Introduction Despite our rapidly expanding knowledge of the structure and function of eukaryotic chromatin at the individual nucleosome level, little is known about the higher-order organization of chromosomes [1]. DNA replication provides an excellent forum with which to investigate these levels of chromosome organization. Large segments of chromosomes replicate coordinately, mediated by the nearly synchronous firing of clusters of replication origins (“replicon clusters”) at specific times during S-phase ([2] and references therein). Replicon clusters can be visualized in living cells as discrete foci by pulse labeling with fluorescent nucleotide analogs. When followed through multiple cell divisions, labeled foci do not mix, separate, or change in shape, indicating that the DNA that replicates coordinately derives from a single chromosome segment [3–6]. Moreover, replicon clusters that fire at different times during S-phase occupy different subnuclear compartments, with early-replicating foci showing enrichment in the nuclear interior, whereas foci replicating later during S-phase are enriched in perinucleolar regions and the nuclear periphery [4,7]. The order in which these segments replicate is established during early G1-phase, coincident with the establishment of their subnuclear positions after nuclear reassembly [4,8]. Together, these results support the hypothesis that coordinately replicated segments of chromosomes form stable units of chromosome structure and nuclear architecture that are maintained from cell cycle to cell cycle. However, the data supporting this model are mainly cytogenetic; molecular evidence for stable replication-timing boundary sequences has been difficult to obtain. Several studies have found that regions of coordinate replication timing correspond to regions of alternating GC content, or isochores, with GC-rich regions replicating early and AT-rich regions replicating late [9–11], leading some to conclude that replication timing is a relatively static chromosomal feature conserved in all cell types [12,13]. However, replication timing cannot be dictated by sequence alone, as both genomic imprinting and X chromosome inactivation are accompanied by the asynchronous replication of homologs [14,15]. Moreover, the replication-timing program for at least some regions is different in different tissues [16–21], and changes in replication timing can be detected during the course of differentiation [22,23]. These differences appear to be related to differential gene expression since there is a strong positive correlation between early replication and transcriptional activity in cultured Drosophila [24,25] and human cell lines [9,10,26], and genes generally replicate earlier when transcriptionally active. This relationship is not direct; instead, early replication appears to be associated with a chromatin state that is permissive for transcription [16]. Nonetheless, it is not clear how much of the genome ever changes replication timing. A comparison of human Chromosome 22 between fibroblast and lymphoblastoid cells revealed that only 1% of this chromosome differed in replication time [9], whereas analyses of individual genes revealed that changes during differentiation were restricted to a subset of genes residing within AT-rich isochores [22]. Hence, determining the extent to which replication timing changes occur during differentiation is a fundamental unresolved question. Major cell fate changes occur early in development when pluripotent cells commit to specific germ layers. Loss of pluripotency can be recapitulated through the differentiation of embryonic stem cells (ESCs), which has been associated with changes in the dynamics of chromatin [27], posttranslational modifications of histones, and nuclear architecture [28]. By focusing on this specific developmental period, we wished to analyze the extent of replication-timing changes genome-wide to address whether replication timing is a static or dynamic property of chromosomes during the course of differentiation. Here, we performed genome-wide analyses of three mouse ESC (mESC) lines before and after differentiation to neural precursor cells (NPCs). Replication domain organization was highly conserved between all three ESC lines and with fibroblasts that were induced to the pluripotent state (induced pluripotent stem cells [iPS cells]). However, 20% of the genome showed substantial changes in replication timing upon neural differentiation. Intriguingly, differentially replicating domains in ESCs consolidated to generate larger coordinately regulated units in NPCs that showed a considerably higher degree of alignment of replication timing to isochore sequence composition. We conclude that replication domain organization is highly dynamic, including its relationship to isochores, and we provide evidence that smaller replication domains that disrupt the alignment of replication timing to isochores define a novel characteristic of the pluripotent state. Furthermore, our findings suggest that DNA replication may provide a molecular handle on the study of previously impenetrable levels of higher-order chromosomal organization. Results Replication Domain Structure in Embryonic Stem Cells The genome-wide analysis of replication timing in mammalian cells has been reported for only one cell type at a density of one probe per megabase [10], which was not sufficient to evaluate the extent to which the genome is organized into coordinately replicating domains. Hence, we mapped replication timing in mESCs using high-density oligonucleotide arrays, adapting a previously developed retroactive synchronization method [24,29]. ESCs were chosen because they provide the opportunity to directly evaluate dynamic changes in the replication program in response to changes in growth conditions [22,23], in contrast to comparisons of separately isolated cell lines that may harbor genetic differences or long-term epigenetic adaptations. Cells were pulse labeled with BrdU and separated into early and late S-phase fractions by flow cytometry (Figure 1A). BrdU-substituted DNA from each fraction was immunoprecipitated with an anti-BrdU antibody, differentially labeled, and cohybridized to a mouse whole-genome oligonucleotide microarray (Nimblegen Systems) (Figure 1A). The ratio of the abundance of each probe in the early and late fraction [“replication timing ratio” = log2(Early/Late)] was then used to generate a replication-timing profile for the entire genome at a density of one probe every 5.8 kb. Replicate experiments in which early- and late-replicating DNA were reciprocally labeled (“dye-switch”) showed a high degree of correlation and were averaged (R 2 values ranged between 0.86 and 0.95 after loess smoothing). Datasets were confirmed by PCR analysis of 18 genes (100% consistent) and by comparison to two previously published replication-timing analyses of 90 individual genes in mESCs (91% consistent) [22,23] (Figure S1A–S1C). Figure 1 Genome-Wide Analysis of Replication Timing in mESCs (A) Protocol for genome-wide replication timing analysis using oligonucleotide microarrays with one probe every 5.8 kb. (B and C) Generating replication-timing profiles. An exemplary mESC replication-timing profile of a Chromosome 1 segment is shown. Raw values for probe log ratios [i.e., log2(Early/Late)] along the chromosome revealed a clear demarcation between regions of coordinate replication (B), which is highlighted upon overlaying a local polynomial smoothing (loess) curve (C). (D) Analyses at a density of one probe per 5.8 kb or 100 bp show essentially identical smoothed replication-timing profiles. Figure 1B shows the mean replication-timing ratio for each probe plotted as a function of chromosomal coordinate for an exemplary 50-Mb segment of Chromosome 1, and Figure 1C shows a loess-smoothed curve fit for the same region. This profile revealed a surprisingly clear demarcation between regions of coordinate replication that we heretofore refer to as “replication domains.” To address whether 5.8-kb probe density was sufficient to provide a complete profile of replication domains, we hybridized the same duplicate preparations of replication intermediates to tiling microarrays (one probe every 100 bp) of Chromosomes 6 and 7. Despite the nearly 60-fold–higher probe density, results showed an almost indistinguishable smoothed profile (Figure 1D). This is consistent with known properties of DNA replication; a 2-h BrdU pulse is expected to label 200–400-kb stretches of DNA (fork rate 1–2 kb/min; [6,30,31]), and since multiple replicons across hundreds of kilobases fire synchronously (reviewed in [2]), probes spaced 5.8 kb apart would be expected to replicate at very similar times. Indeed, high autocorrelation of replication timing between neighboring probes extends over 1 Mb (Figure S2). Hence, replication timing across the entire genome can be reliably profiled on a single oligonucleotide chip. Replication profiles for all chromosomes are displayed on our Web site and are available for downloading: http://www.replicationdomain.org. To quantify the numbers and positions of replication domains and their boundaries genome-wide, we adapted a segmentation algorithm—originally developed to identify copy number differences for comparative genomic hybridization [32]—to identify regions of uniform y-axis values (see Materials and Methods), which are illustrated in Figure 2A. This algorithm generates a dataset consisting of the nucleotide map positions for the boundaries of each replication domain. Domain sizes ranged from 200 kb to 2 Mb, with some considerably larger domains (Figures 2B and S3A). These domain sizes provide a logical explanation as to why existing ENCODE replication-timing data for HeLa cells [33] does not reveal replication domains; the ENCODE regions cover 1% of the genome and consist primarily of scattered 500-kb genomic segments, which would be too small to discern replication domain–level chromosome organization. It is also possible that the genetic and epigenetic instability of HeLa cells contributed to the blurring of domain boundaries. Domains were found to replicate at all times during S-phase, however, domains larger than 2.5 Mb were either very early or very late replicating, suggesting that coordinately replicating regions larger than a certain threshold size tend to replicate at one extreme or another of S-phase (Figure S3D). Our results were not an artifact of probe density, segmentation algorithm, or synchronization method; similar distributions were obtained with a density of one probe per 100 bp, using different segmentation parameters, and using an alternative protocol [10] that determines replication timing by probe copy number in S-phase versus G1-phase, without fractionation of S-phase (unpublished data). Similar results were also obtained with human ESCs (hESCs; J. Lu, I. Hiratani, T. Ryba, and D. M. Gilbert, unpublished data). Figure 2 Replication Domain Structure and Its Conservation between Three Independent mESC Lines (A) Identification of replication domains (red lines) and their boundaries (dotted lines) by a segmentation algorithm [32]. (B) Box plots of early (E; log ratio > 0) and late (L; log ratio 0.65 kb) to generate the final list of 18,679 RefSeq genes with replication-timing ratios matched. Complete replication-timing datasets for all (384,849) probes are downloadable from our Web site (http://www.replicationdomain.org ) and are graphically displayed on the Web site. Transcription analysis by microarrays. Total cellular RNA was isolated from D3 ESCs or NPCs (three biological replicates per cell state), and steady-state transcript levels were determined by Affymetrix GeneChip microarrays (Mouse Genome 430 2.0), which were highly reproducible (R 2 > 0.98 between all replicates). After quality control tests [79], datasets were subject to normalization by the Probe Logarithmic Intensity Error algorithm (PLIER) developed by Affymetrix for calculating probe signals. For each Affymetrix “probe set,” signal intensity of the three biological replicates were averaged (i.e., average intensity). Genes are often represented by multiple probe sets. In such cases, the one with the highest total intensity (i.e., sum of ESC and NPC average intensity) was defined as the representative probe set, and the other probe sets were not used. We did so because such highest-intensity probe sets were empirically most consistent with reverse transcriptase (RT)-PCR analysis and can be defined in an objective way. Present (transcriptionally active) and absent (inactive) calls are generated by MAS5.0 (Affymetrix) per replicate per probe set, which results in multiple present–absent calls for a given gene [=3 × (total number of probe sets for a gene)]. We defined “present” genes as those with more than 50% of all their probe set calls being present. A total of 15,143 (81%) of the 18,679 RefSeq genes, for which replication-timing ratios were obtained, were represented on the Affymetrix GeneChip microarrays and were assigned transcription levels and present–absent calls. Validation of transcription array results was evident from previously published transcription analysis under the same condition [80]. Identification of replication domains and domains that change replication timing. DNAcopy (R/Bioconductor) is a segmentation algorithm for the analysis of microarray-based DNA copy number data [32]. For identification of replication domains, we applied this method directly to datasets containing mean replication-timing ratios for all probes before loess smoothing. The parameters, nperm (number of permutation) and alpha (the significance level for the test to accept change-points), were set at 10,000 and 1 × 10−15, respectively, which were empirically determined based on how well the resultant segmentation profile traced the loess-smoothed profile. Once determined, these parameters were fixed and used for objective segmentation of all datasets. Although DNAcopy segmentation infrequently fails to identify segments discernible by visual inspection of loess-smoothed profiles, it does so equally for ESC and NPC profiles. Thus, its performance was sufficient to provide objective evidence for replication domain consolidation. Others have reported the superiority of its overall performance to current alternatives [81,82]. The same strategy was used to identify chromosomal domains that change replication timing, except in this case, datasets consisting of replication timing ratio differential (i.e., NPC ratio − ESC ratio) for all probes were used for segmentation. Among the resultant 2,042 segments, we selected 102 EtoL, 102 LtoE, 232 EtoE, and 96 LtoL domains based on the criteria described in Figure 4G. Analysis of transitions between replication domains. We selected the three chromosomes 2, 11, and 16 because we reasoned that they were representative of the genome. Chromosomes 2, 11, and 16 are large-, intermediate-, and small-size chromosomes with gene density that is intermediate, high, and low, respectively. For random selection of replication domain boundaries, we focused on the middle portion of each chromosome, counting all transitions after nucleotide position 40,000,000 on each chromosome until we counted 25 boundaries. As a result, the following regions were analyzed: chr2:40,000,000–75,000,000; chr11:40,000,000–68,000,000; and chr16:40,000,000–65,000,000. Transition regions were defined as regions with large and unidirectional changes in replication timing along the chromosomes on the loess-smoothed curve. The positions at which this unidirectionality stopped were defined as the two “ledges” of a transition region. GC and LINE-1 content calculation. GC and LINE-1 content was calculated based on the UCSC Genome Browser database (gc5base.txt and chrN_rmsk.txt, mm8 assembly; http://genome.ucsc.edu) using the Table Browser function of the UCSC Genome Browser as well as an R/Bioconductor script. DNA-FISH. DNA-FISH was performed essentially as described [62], with some modifications. Briefly, preparation and fixation of cells were done as described [83] to preserve 3-D structure. BAC probes were used for all genes tested, with some genes additionally tested by PCR probes of 8.9–10.2 kb. Digoxigenin (DIG)-labeled probes were generated using the DIG-nick translation mix (Roche, Cat#11745816910). Primary and secondary antibodies used to detect the DIG-labeled probes were sheep anti-DIG-fluorescein (Roche Applied Science, Cat# 11207741910) and rabbit fluorescein anti-sheep IgG (Vector, Cat#FI-6000), respectively. Images were captured with a DeltaVision Image Restoration Microscope System (Applied Precision) attached to an Olympus IX-71 fluorescence microscope equipped with an Olympus PlanApo 100× 1.42 NA oil objective lens. Optical sections were taken with 0.2-μm spacing and were subsequently enhanced using constrained iterative deconvolution process by softWoRx software (Applied Precision). We defined the radius of each nucleus as one half of the largest diameter of DAPI staining, which decreased slightly upon differentiation but was comparable in ESC and NPC (average radius: 5.3 μm in ESC, n = 1,250 vs. 4.9 μm in NPC, n = 1,339). We then measured the distance from FISH signals to the nearest nuclear periphery, and divided it by the radius to obtain relative radial distance to the nuclear periphery. RNA-FISH. LINE-1 RNA-FISH was performed essentially as described [84]. LINE-1 primer sequences were 5′-TAATACGACTCACTATAGGGGGCTCAGAACTGAACAAAGA-3′ (forward; underline, T7 promoter) and 5′-GCTCATAATGTTGTTCCACCT-3′ (reverse), which amplifies a 1,041-bp fragment of LINE-1 corresponding to portions of ORF2 and the 3′-UTR (L1MdA2; accession number M13002; 7,713 bp). Importantly, this sequence is conserved in other subfamilies of LINE-1. We used genomic DNA for PCR, and the amplified DNA fragment was purified and used for in vitro transcription, followed by reverse transcription to generate a DIG-labeled, single-stranded DNA probe. Supporting Information Figure S1 Validation of Microarray-Based Replication-Timing Analysis by PCR (A) Validation of microarray experiments by PCR. Pairs of immunoprecipitated BrdU DNA samples from early and late S fractions were subject to PCR and mean percent early S-phase values (i.e., [intensity of early fraction]/[intensity of early and late fractions combined]) from six to seven pairs of DNA samples were calculated, as previously described [22]. Genes above and below 50% were classified as early (E) and late replicating (L), respectively. From microarray data, replication-timing ratios of genes were obtained from the loess-smoothed curve at the transcription start sites. Replication-timing ratios above and below zero were classified as early (E) and late replicating (L), respectively. The resultant binary datasets for 18 genes showed a 100% match (18/18) in ESCs and a 94% match (17/18) in NPCs. Note that this binary classification of PCR results forces some genes that actually change replication timing to be not classified as such: for instance, Crisp1 (later shift), Cdh2, Postn, and Mash1 (earlier shift). However, even such subtle changes were detected on the microarray, as shown by the changes in replication timing ratios from ESCs to NPCs. (B and C) Comparison to two previously published replication-timing analyses using 46C ESCs [22] (B) and OS25 ESCs [23] (C). (B) PCR results from Hiratani et al. [22] were classified as early (E) and late (L) based on the same criteria as in Figure S1A. (C) Genes called E, ME, and M by Perry et al. [23] were classified as early (E), whereas genes called ML and L were classified as late (L). Both studies combined, 91% (82/90) of the PCR results matched those from the microarray. (51 KB PDF) Click here for additional data file. Figure S2 Autocorrelation Analysis of Replication-Timing Data Autocorrelation analysis of replication timing in ESCs. The autocorrelation function (ACF) indicates the degree of similarity between neighboring data points. The result illustrates that relatively uniform replication timing extends over large regions of approximately 1 Mb. The x-axis shows the distance on the chromosome in megabases as calculated from the interprobe distance. (36 KB PDF) Click here for additional data file. Figure S3 Replication Domain Size Distribution (A–C) Size distribution of replication domains in ESCs (A) and NPCs (B), as well as domains that change replication timing ([C] EtoL and LtoE domains). Domains with replication-timing ratios above and below zero were defined as early- and late-replicating domains, respectively. Domains were categorized by their sizes into bins of equal intervals (0.2 Mb) starting from 0–0.2 Mb as the first bin. Insets in (A–C) show domains below 0.4 Mb in bins of equal intervals (40 kb) starting from 0–40 kb as the first bin. (D) Scatter plots of replication timing versus domain size in ESCs and NPCs. (67 KB PDF) Click here for additional data file. Table S1 Replication-Timing Ratio, mRNA Expression, Histone Modifications [51,53,54], and Promoter CpG Classification [51] of 18,702 RefSeq Genes Analyzed in This Study (6.25 MB XLS) Click here for additional data file.
