The Temporal Order of DNA Replication Shaped by Mammalian DNA Methyltransferases

Multiple epigenetic pathways underlie the temporal order of DNA replication (replication timing) in the contexts of development and disease. DNA methylation by DNA methyltransferases (Dnmts) and downstream chromatin reorganization and transcriptional changes are thought to impact DNA replication, yet this remains to be comprehensively tested. Using cell-based and genome-wide approaches to measure replication timing, we identified a number of genomic regions undergoing subtle but reproducible replication timing changes in various Dnmt-mutant mouse embryonic stem (ES) cell lines that included a cell line with a drug-inducible Dnmt3a2 expression system. Replication timing within pericentromeric heterochromatin (PH) was shown to be correlated with redistribution of H3K27me3 induced by DNA hypomethylation: Later replicating PH coincided with H3K27me3-enriched regions. In contrast, this relationship with H3K27me3 was not evident within chromosomal arm regions undergoing either early-to-late (EtoL) or late-to-early (LtoE) switching of replication timing upon loss of the Dnmts. Interestingly, Dnmt-sensitive transcriptional up- and downregulation frequently coincided with earlier and later shifts in replication timing of the chromosomal arm regions, respectively. Our study revealed the previously unrecognized complex and diverse effects of the Dnmts loss on the mammalian DNA replication landscape.


Introduction
Mammalian chromosomes consist of multiple replication origins. Coordinate activation of multiple origins across the scale of several hundred kilobases to several megabases forms replication domains, the regulatory units of DNA replication [1][2][3]. Genome-wide replication timing analysis has revealed the periodic distribution of early and late replication domains along the chromosome arms [4]. Replication domains are reorganized in the contexts of development and disease, giving rise to cell type-specific organization of replication domains [4]. For example, the chromosomal region that contains pluripotencyassociated genes such as Dppa2/Dppa4 displays early replication in mouse ES cells, while in closely related epiblast stem cells, the same region switches to late replication [5]. Replication timing of some genomic regions is commonly dysregulated in pediatric leukemia patient cells [6]. Given these observations, epigenetic mechanisms are thought to play important roles in the temporal regulation of DNA replication. To date, however, only a handful of epigenetic modifiers have been tested to determine whether and to what extent they contribute towards the DNA replication program [7].
DNA methylation is an epigenetic modification found in gene regulatory regions, heterochromatin, DNA repeats, and transposable elements [8,9]. DNA methylation patterns in the genome are stably inherited over many cell divisions by the cooperative action of the DNA methyltransferases (Dnmts) [10,11]. Regulation of gene transcription by DNA methylation has been well studied. Mechanistically, methylated CpGs are recognized by a set of methyl CpG-binding proteins (MBDs and the BTB/POZ domain family of proteins) that recruit chromatin-modifying complexes and bring about repressive chromatin environments near gene promoter regions [12]. Mouse pericentromeric heterochromatin (PH) forms a microscopically observable nuclear domain that contains highly methylated DNA [13]. Not only Dnmts, but also MDBs are localized at PH, thus providing a model to examine the relationship between DNA methylation and histone modifications [14,15]. Genetic inactivation of genes responsible for DNA methylation leads to aberrant gene expression patterns [16][17][18]. Several biochemical studies have shown the interaction of DNA methyltransferases with various histone modifiers, suggesting that non-catalytic functions of Dnmts also contribute to the regulation of gene expression [19,20]. However, the roles of the Dnmts other than gene regulation are poorly understood. DNA hypomethylation is associated with the promotion of tumorigenesis and other disease states [21,22], and a causal relationship between DNA hypomethylation and chromosome instability has been reported; however, the underlying mechanisms remain unclear [23][24][25].
Although genome-wide replication timing analysis has not been performed in Dnmtdeficient cells, a global correlation between DNA methylation and early replication has been observed [26]. This is somewhat unexpected since the Dnmts and DNA methylation generally induce repressive chromatin environments that favor late replication. To directly examine the extent to which the loss of Dnmts and DNA methylation affect the DNA replication program, we performed detailed replication timing analysis using genetically engineered Dnmt-mutant mouse ES cell lines.
Dnmt3a −/− Dnmt3b −/− double-knockout (DKO) ES cell clones stably expressing WT or catalytic-defective (C487S) Dnmt3a2 were generated by electroporation-mediated introduction of cDNA expression plasmids driven by the CAG promoter, followed by selection with blasticidin. An empty expression plasmid was used to generate control mock clones.
Dnmt3a −/− Dnmt3b −/− DKO ES cell clones stably expressing RU486 (mifepristone, a progesterone receptor antagonist)-inducible Dnmt3a2 were generated by electroporationmediated introduction of an expression plasmid for Dnmt3a2 fused with a modified ligandbinding domain of the human progesterone receptor (termed "GPRh") driven by the CAG promoter, followed by selection with blasticidin. A mammalian expression vector with the RU486-inducible GPRh was provided by Hitoshi Niwa, Kumamoto University. The GPRh is a chimeric protein comprising the human glucocorticoid receptor (513-554; UniProt ID, P04150) and the human progesterone receptor (710-914; UniProt ID, P06401), in which the N-terminal helices 1 and 2 in the ligand-binding domain of the progesterone receptor are replaced with those of the glucocorticoid receptor to increase the capacity for cytosolic retention in the absence of ligand [28], and the C-terminal 19 amino acids of the progesterone receptor are deleted allowing response to RU486 but not to progesterone [29,30]. Five amino acids (MHAYG) derived from the cloning site are added at the C-terminal end. To induce expression of the Dnmt3a2 protein, ES cells (clone 1723) were cultured in ES medium containing 1 µM RU486 for 14 days on feeder mouse embryonic fibroblast (MEF) cells. Control cells were treated with ethanol only.

Genome-Wide Replication Timing Analysis
Replication timing analysis with BrdU labeling was performed as described previously [31][32][33] with modifications. In brief, cells were labeled with 50 µM BrdU for 2 h, washed twice with ice-cold PBS, trypsinized, and then fixed in 75% (v/v) ethanol. These cells were resuspended in PBS containing 1% FBS, stained with propidium iodide (20 µg/mL) for 30 min in the presence of RNase A (250 µg/mL), and then sorted into early and late S phase fractions by flow cytometry (SH800 Cell Sorter; Sony Biotechnology, Tokyo, Japan). After phenol-chloroform extraction of the cellular DNA, immunoprecipitation with anti-BrdU mouse monoclonal antibody (Cat#555627; BD Biosciences, San Jose, CA, USA) was performed for each fraction to enrich BrdU-substituted replicating DNA. Isolated early and late replicating DNA were amplified by whole-genome amplification (WGA) (SeqPlex Enhanced DNA Amplification Kit; Sigma-Aldrich, St. Louis, MO, USA) and used to construct libraries for sequencing with an Ion Proton next generation sequencer (Life Technologies, Waltham, MA, USA). Library preparation and data acquisition were performed according to standard Ion Proton procedures. Data analyses were performed with SeqMonk (http://www.bioinformatics.babraham.ac.uk/projects/seqmonk/) and R/Bioconductor (http://www.r-project.org; http://www.bioconductor.org). The number of total reads in each sample were normalized per million reads in SeqMonk. The ratio of sequence reads between early and late S phase samples was determined in windows of 5 kb and LOESS-smoothed over a 300-kb window size in R/Bioconductor. These smoothed datasets were used to generate the replication timing profiles.
For some analyses, datasets were averaged into 200-kb windows (fixed position) and the replication timing differential (i.e., TKO ratio − WT ratio) was determined for each 200-kb segment. To determine the significant replication timing switching domains that were independent of changes between replicates, we calculated the Euclidian distance at 13,222 200-kb segments between groups (i.e., WT vs. TKO) and within groups (i.e., WT replicate-1 vs. WT replicate-2), which was used to calculate the p-values at each 200-kb genomic segment [34]. Statistical significance was then calculated using the qvalue package in R/Bioconductor, which yields a q-value for each segment that reflects the proportion of false positives (i.e., the false discovery rate; FDR) among segments deemed to have significant replication timing changes. High-confidence replication timing switching domains were selected with a q-value cutoff of 0.01, corresponding to an overall FDR of 1%. A q-value cutoff of 0.05 was also used to identify a set of lower-confidence domains.

DNA Combing
DNA combing and immunofluorescence detection of replication-labeled DNA were performed as described previously [35]. Cells were first labeled with 100 µM 5-iodo-2deoxyuridine (IdU; Sigma-Aldrich, St. Louis, MO, USA) for 20 min. After washing with PBS, the cells were labeled with 100 µM 5-chloro-2 -deoxyuridine (CldU; Sigma-Aldrich, St. Louis, MO, USA) for 20 min. The labeling reaction was stopped by washing the cells with ice-cold PBS. The cells were collected by centrifugation at 400× g for 3 min at 4 • C, and the supernatant was removed before adding buffer A (250 mM sucrose, 20 mM HEPES, pH 7.0, 10 mM KCl, 1.5 mM MgCl 2 , 1 mM EDTA, 1 mM EGTA) to the cell pellets. Then, agarose LM-MP (Roche, Mannheim, Germany) was added to a final concentration of 0.7% (w/v), and the resultant cell suspension was transferred into a plug mold. After 2 h of incubation on ice, the solidified gel containing the cells was incubated in ESP solution (10 mM Nlauroylsarcosine sodium salt, 10 mM Tris-HCl, pH 8.0, 0.5 M EDTA, 2 mg/mL Proteinase K) for 16 h at 50 • C to digest the proteins. Next, the gel was washed with 0.5 M EDTA three times and then washed twice with 1× TE (10 mM Tris-HCl, pH 8.0, 1 mM EDTA). After 1 h of incubation at room temperature, the gel was washed with 1× TE and incubated with 50× β-agarase buffer (Lonza Japan, Tokyo, Japan) for 1 h. Then, the gel was melted at 70 • C, and β-agarase I (Lonza Japan, Tokyo, Japan) was added to complete the gel dissolution. Finally, the DNA was diluted with 0.5 M MES, pH 5.5 and stored at 4 • C until use. DNA molecules isolated in the MES buffer were combed at a rate of 300 µm/s onto silanated coverslips (Matsunami Glass, Osaka, Japan) as described previously [36]. To estimate the degree of DNA molecule stretching on the coverslips, λ-DNA molecules were similarly combed and stained with YOYO-1 for length measurement. Because the λ phage virus genome is 48.5 kb in length, the rate of DNA stretching was calculated as 2.32 ± 0.11 kb/mm. Combed DNA molecules were denatured in 50% formamide and 2× SSC (saline-sodium citrate) buffer for 12 min at 72 • C. Denatured DNA molecules were incubated in blocking solution (3% Blocking Reagent from Sigma-Aldrich, and 0.05% Tween 20 in PBS) for 30 min at 37 • C to reduce nonspecific binding, followed by incubation in a detection solution (1% Blocking Reagent, 0.05% Tween 20 in PBS) containing the primary antibodies for 1 h at 37 • C. A mouse anti-BrdU monoclonal antibody (1:5 dilution, clone B44, mouse IgG1, Cat#347580; BD Biosciences, San Jose, CA, USA) and rat anti-BrdU monoclonal antibody (1:50 dilution, clone BU1/75, rat IgG2a, Cat#ab6326; Abcam, Cambridge, UK) were used for immunodetection of IdU-and CldU-labeled DNA, respectively. After washing with 0.05% Tween 20/PBS, the DNA molecules were incubated in detection solution containing Alexa Fluor 555-conjugated goat anti-mouse IgG (1:100 dilution; Life Technologies, Carlsbad, CA, USA) and Alexa Fluor 488-conjugated rabbit anti-rat IgG (1:100 dilution; Life Technologies, Carlsbad, CA, USA) for 30 min at 37 • C. After washing with 0.05% Tween 20/PBS, the coverslips were mounted in Vectashield (Vector Laboratories, Burlingame, CA, USA).

Visualization of Replication Sites in the Cell Nucleus
Cells were labeled with 100 µM EdU (Life Technologies, Carlsbad, CA, USA) for 10 min, aliquoted onto coverslips and fixed with 4% paraformaldehyde in PBS, and then permeabilized with 0.5% Triton X-100 in PBS. Blocking was performed by incubation in 3% BSA/0.1% Tween 20/Tris-buffered saline, pH for 30 min. EdU-labeled DNA was detected with the Click-iT EdU Alexa Fluor 488 Imaging Kit (Life Technologies, Carlsbad, CA, USA) according to the manufacturer's protocol. Digoxigenin (DIG)-labeled DNA was detected with an anti-DIG antibody conjugated with FITC (1:200 dilution; Roche, Mannheim, Germany). DAPI was used to counterstain the DNA. In the experiment in Figure 1B, digoxigenin-11-dUTP (DIG-dUTP; Roche, Mannheim, Germany) was loaded into cells with the hypotonic shift method for replication labeling [37]. Incorporated DIG-dUTP was detected with an anti-DIG antibody conjugated with rhodamine (Cat#11207750910; Roche, Mannheim, Germany).

Microarray Analysis
Total cellular RNA from two independent experiments was reverse transcribed and hybridized to the Affymetrix Mouse Genome 430 2.0 microarray, as described previously [38]. Data processing, quality control, and output file generation was performed with Affymetrix Gene Expression Console. Raw data (CEL files) were used to perform normalization through the Robust Multiarray Analysis (RMA) algorithm.

Methylation Analysis by Southern Blotting
Genomic DNA was digested with CpG methylation-sensitive restriction enzymes (HpaII or MaeII), blotted, and hybridized with major satellite (pSAT) probes as described previously [39].

Imaging System and Measurements
Images were collected using a Leica DM RA2 fluorescence microscope equipped with a cooled charge-coupled device camera (C4742-95-12ER; Hamamatsu Photonics, Shizuoka, Japan), controlled by an Apple Macintosh G4 computer running the software program IPLab (Signals Analytics, Vienna, VA, USA). The images were captured at different stage positions and processed with the Huygens Essential deconvolution software (Scientific Volume Imaging, Hilversum, Netherlands). For some experiments, a fluorescence microscope (Axioplan 2 MOT; Carl Zeiss, Jena, Germany) equipped with an ORCA R2 camera (Hamamatsu Photonics, Shizuoka, Japan) was used for imaging.

H3K27me3 Foci Formed in Mouse Embryonic Stem Cells with Severely Hypomethylated DNA Coincide with Later Replication of Pericentromeric Heterochromatin
Intranuclear distribution patterns of DNA replication sites are known to change dynamically during S phase: DNA replication initially occurs throughout the interior part of the nucleus, then in the region of the pericentromeric heterochromatin (PH) that is densely stained with DAPI, and last in the nuclear periphery [40,41]. To examine the effect of Dnmt loss on the DNA replication program, we first investigated the spatial organization of DNA replication sites in Dnmt1 −/− Dnmt3a −/− Dnmt3b −/− triple-knockout (TKO) mouse ES cells. We observed that the S phase stage-specific distribution patterns of DNA replication sites (patterns I-V) and their frequency in TKO ES cells were comparable to those in wild-type ES cells ( Figure 1A), suggesting that spatial regulation of the DNA replication program is largely maintained in the absence of the Dnmts. We also did not detect a significant difference between wild-type and TKO ES cells in the rate of replication fork progression and ori-to-ori distance as measured by a DNA combing assay ( Figure S1).
A previous report demonstrated that H3K27me3 foci are formed in the PH of extensively DNA hypomethylated cells [42]. Consistent with this, we observed H3K27me3 foci formation at major satellite repeats, which are a main structural component of PH ( Figure S2A,B), and regional exclusion of normally enriched constitutive heterochromatin markers such as HP1β, H3K9me3, and H4K20me3 ( Figure S2C,D). Although H3K27me3 foci are formed in a DNA methylation-dependent manner ( Figure S3), a portion of PH in TKO ES cells still maintains the constitutive heterochromatin marks without H3K27me3 enrichment ( Figure S2C,D), suggesting that loss of DNA methylation destabilizes constitutive chromatin structure and, as a consequence, indirectly induces regional H3K27me3 foci formation within PH. Interestingly, these H3K27me3-enriched PH regions replicate later than neighboring PH regions where H3K27me3 signal intensity is low ( Figure 1B, colocalization of H3K27me3 foci with PCNA-positive later replicating PH but not with DIG-dUTP-positive earlier replicating PH). In H3K27me3-enriched PH regions, we also observed loss of Aurora B recruitment and H3S10 phosphorylation, which are known to coincide with replication timing delay ( Figure 1C) [43]. Because the temporal order of replication within PH is closely correlated with the level of H3K27me3, it is suggested that sporadic chromatin reorganization induced by DNA methylation loss, rather than DNA methylation loss itself, affects the temporal order of DNA replication within PH.

Genome-Wide Replication Timing Analysis Identified a Subset of Chromosomal Regions Sensitive to Dnmt Loss
To systematically explore the effect of Dnmt loss on DNA replication, we next profiled the genome-wide replication timing of TKO ES cells in comparison with the parental control, wild-type ESC cells. To this end, BrdU-immunoprecipitated DNA from FACS-sorted early and late replicating fractions was analyzed using next generation sequencing (Figure 2A), and the ratio of early and late read enrichment (Log 2 early/late) was determined against the chromosomal position ( Figure 2B). Although the obtained replication timing profile of the TKO ES cells was globally similar to that of control ES cells, our detailed analysis identified a set of chromosomal regions that underwent switches from early to late (EtoL) and from late to early (LtoE) after Dnmt loss ( Figure 2B). To evaluate the statistical significance of these replication timing changes, we averaged datasets into 200-kb windows (total 13,222 200-kb segments/dataset) and calculated the proportion of differences greater than those observed between biological replicates. These proportions were converted to false discovery rate (FDR; see details in Materials and Methods) estimates to account for multiple testing. Using this approach, we identified 628 and 300 200-kb segments undergoing either EtoL or LtoE switching, respectively, with an FDR of 1% FDR ( Figure 2B).
To refine the genomic regions affected by the Dnmts, we also examined genomewide replication timing in a Dnmt3a −/− Dnmt3b −/− double-knockout (DKO) cell line stably expressing RU486-inducible Dnmt3a2. Although we observed only partial restoration of DNA methylation upon Dnmt3a2 induction ( Figure S4A), the restored DNA methylation level was sufficient to suppress H3K27me3 foci formation within the PH ( Figure S4B,C). Using this inducible expression system, we found that induction of Dnmt3a2 also induced changes in replication timing in a subset of chromosomal regions ( Figure 2C). As expected, regions undergoing EtoL switching upon Dnmt loss (WT vs. TKO) significantly overlapped with regions undergoing LtoE switching upon Dnmt3a2 induction (DKO − RU486 vs. DKO + RU486) ( Figure 2D top). LtoE switching regions upon Dnmt loss also significantly overlapped with EtoL switching regions upon Dnmt3a2 induction, but to a lesser extent (Figure 2D bottom). These changes in replication timing were not as dramatic as changes seen in other replication domain mutants [44][45][46] but were highly reproducible between experimental replicates ( Figure S5). A previous study showed that the pluripotency of ES cells is maintained in the absence of the Dnmts [27], which was further confirmed by the fact that the mutant cells used in our experiments retained pluripotent cell-specific replication timing profiles ( Figure S6). Taken together, these results clearly demonstrate that the Dnmts are required to maintain the temporal order of DNA replication in a number of chromosomal arm regions in addition to PH.

Redistribution of the H3K27me3 Mark Is not Associated with Replication Timing Changes in the Chromosomal Arm Regions
It has been shown that genome-wide redistribution of H3K27me3 occurs in Dnmt TKO ES cells [42]. In these cells, both abnormal enrichment and loss of H3K27me3 were observed in regions where this modification is normally present [42,47]. This, together with the observation that H3K27me3 enrichment coincides with later replicating PH ( Figure 1B), led us to examine whether replication timing changes in the chromosomal arm regions of TKO ES cells were associated with changes in H3K27me3 enrichment. We analyzed H3K27me3 redistribution with replication timing changes using H3K27me3 ChIP-seq data from control wild-type and TKO ES cells [48] and found no appreciable correlation between the two ( Figure 3A, R = 0.022). Genomic regions undergoing EtoL replication timing switching after Dnmt loss (regions identified in Figure 2D) exhibited no enrichment of the H3K27me3 mark but rather lost this modification, and this was also the case with the LtoE and LtoL domains ( Figure 3B). Further investigation of replication timing in genomic regions that scored in the top 1%, 5%, and 10% in terms of H3K27me3 gain and loss showed poor coordination between H3K27me3 redistribution and the replication timing changes ( Figure 3C). These results suggest that the association of the H3K27me3 mark with the temporal regulation of DNA replication is restricted to specific regions of the genome such as the PH and that other factors account for replication timing changes in the chromosomal arm regions.

Replication Timing Changes Frequently Coincide with Transcriptional Changes in the Chromosomal Arm Regions
Since the Dnmts regulate the expression of a number of genes, we next examined the relationship between replication timing changes and transcriptional changes. We performed genome-wide gene expression analysis in DKO cells with and without Dnmt3a2 induction. As expected, several genes within the MageA and Rhox gene clusters known to be regulated by DNA methylation [18,49] were downregulated upon Dnmt3a2 induction (data not shown). In addition, a subset of genes was upregulated, which might have been partly due to the indirect effects of Dnmt3a2 induction. We found that EtoL switching regions upon Dnmt3a2 induction tended to harbor downregulated genes ( Figure 4A,B). In LtoE switching regions, coordinate changes in replication and transcription were not evident when taking all the expressed genes within these regions into account ( Figure 4B). However, a subset of highly upregulated (>two-fold) genes was found to be enriched in the LtoE switching regions ( Figure S7). This led us to focus on the genes that were highly upand downregulated upon Dnmt3a2 induction to examine their relationship with replication timing changes ( Figure 4C). Since most of the highly up-and downregulated genes were not localized in the EtoL and LtoE switching regions ( Figure S7), we extended the analysis to the chromosomal regions harboring the genes scoring in the top 1%, 5%, and 10% in terms of up-and downregulation, regardless of their statistical significance in replication timing switching. Interestingly, we found that highly up-and downregulated transcription was associated with shifts in replication timing towards earlier and later timing, respectively. A similar but less pronounced result was obtained from the comparison between WT and TKO ( Figure S8). Collectively, these results suggest that transcriptional changes, either directly or indirectly induced by Dnmt loss, could be a cause or consequence of replication timing changes in the chromosomal arm regions.

Discussion
Using Dnmt-deficient mouse ES cells as a model, our study uncovered widespread effects of Dnmt loss on the mammalian DNA replication landscape and also highlighted the complexity of the mechanisms behind induced changes in replication timing.
Within PH regions, loss of DNA methylation destabilizes constitutive heterochromatin structure, leading to regional chromatin reorganization that involves abnormal enrichment of H3K27me3 foci. Loss of DNA methylation along CG-rich sequences is known to allow targeting of polycomb protein KDM2B, which promotes the formation of the H3K27me3 mark via PRC2 recruitment [50,51], and this is likely also the case for H3K27me3 foci formation within the PH of TKO cells. Consistent with this, we observed enrichment of polycomb group proteins at H3K27me3 foci (data not shown). Because H3K27me3-enriched regions replicate later than neighboring PH regions where constitutive heterochromatin marks are normally maintained without DNA methylation, we concluded that while the loss of DNA methylation is necessary, it is not alone responsible for the later replication of PH-chromatin reorganization is also required. Although we do not exclude the possibility that H3K27me3 affects the temporal order of PH replication, it is also possible that either gain or loss of other epigenetic modifications is involved in this phenomenon. For example, Suv4-20h-mediated H4K20me3 has been shown to be necessary for origin activation in a certain genome context [52], and we observed loss of this modification at H3K27me3 foci ( Figure S2C,D).
A previous study has shown that partial loss of DNA methylation within PH advances replication timing [53]. For example, Dnmt1 KO ES cells in which 30% of the methylation is still maintained show earlier replicating major satellite sequences than control wildtype ES cells [53]. Further extensive loss of DNA methylation induces H3K27me3 foci formation within PH and may shift this region back towards late replication. In other words, replication timing of PH may fluctuate depending on the DNA methylation level. It should be noted that H3K27me3 foci are not formed in Dnmt1 single KO cells (data not shown). Abnormal enrichment of the H3K27me3 mark was also observed in histone H3K9 methyltransferase Suv39h1/h2 double-knockout ES cells [54], and delayed replication timing of PH was induced in these cells [53]. Although these lines of evidence suggest the close relationship between the H3K27me3 mark and replication timing regulation, the detailed mechanistic link behind this relationship is currently unclear.
Interestingly, the relationship between H3K27me3 and the replication timing changes discussed above seems to be restricted to PH. Our genome-wide analysis revealed that altered enrichment of the H3K27me3 mark in TKO ES cells is largely uncorrelated with changes in replication timing, suggesting that the mechanism by which the loss of DNA methylation affects the temporal order of PH replication is not involved in the replication timing changes seen in the chromosomal arm regions. Thus, multiple mechanisms are thought to be involved in the observed replication timing changes in Dnmt-deficient ES cells.
We showed that earlier and later shifts in replication timing are often associated with transcriptional up-and downregulation, respectively. As this trend is evident when major changes in transcription occur in response to Dnmt3a2 reintroduction into DKO ES cells, transcription itself or the associated chromatin changes could be a driving force behind the replication timing changes. However, we do not exclude the possibility that replication timing changes affect transcriptional activity. It has been shown that plasmid DNA injected into early and late S phase cell nuclei preferentially forms active and inactive chromatin structures, respectively, via an as yet unknown mechanism [55]. This replication timing-dependent assembly of different types of chromatin could have an impact on gene transcription. The relationship between transcription and DNA replication could, thus, be interdependent through a positive feedback loop, making it difficult to determine which is the cause and which the consequence. Coordinate changes between replication timing and transcription are also observed during mouse ES cell differentiation [4,5]. Genes in EtoL and LtoE switching regions are often down-and upregulated, respectively. Analysis at multiple intermediate differentiation stages has revealed no consistent relationship between transcription and replication that definitively illustrated which change is an upstream event: Some genes undergo transcriptional changes prior to replication timing switching, whereas others undergo transcriptional changes that follow replication timing switching [5].
Considering that significant changes in epigenetic modifications and gene transcription are induced in Dnmt-deficient cells, it may be surprising that the observed changes in replication timing are relatively subtle and that ES cell-specific Mb-sized replication domain organization is globally maintained. However, this was supported by recent electron spectroscopic and Hi-C experimental findings that higher-ordered chromatin structures and A/B compartment structures are largely maintained in the absence of the Dnmts [56,57]. The observed robustness of the replication domain structure could be partly explained by the replication timing program being regulated at two different levels (i.e., at the replication domain level and the replication origin level) in mammalian cells [3,7]. Unlike in yeast, replication timing in mammalian cells is first determined at the level of Mb-sized replication domains during the early G1 phase, which determines the domains to be replicated during the early or late S phase. In the subsequent step during the G1 phase, replication origin selection occurs to determine the local order of replication within each individual (early or late) domain. If regional changes in chromatin structure or transcriptional activity that affect local origin activity occur in response to Dnmt loss, the effects would tend to be buffered by the replication forks of neighboring origins within the same domain (the distances between adjacent origins are estimated to range between 40-150 kb based on DNA fiber experiments), resulting in the robustness of Mb-sized replication domain organization.

Conclusions
In summary, our study elucidated the extent to which DNA methylation and the Dnmts contribute to replication timing regulation. This will provide a basis for understanding diseases associated with chromatin states and replication timing dysregulation.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4 409/10/2/266/s1, Figure S1: DNA replication fork progression rate and ori-to-ori distances as measured in Dnmt-depleted cells, Figure S2: Extensive loss of DNA methylation results in the formation of H3K27me3 foci within pericentromeric major satellite repeats, Figure S3: H3K27me3 foci are formed in a DNA methylation-dependent manner, Figure S4: Establishment of a DKO ES cell clone stably expressing RU486-inducible wild-type Dnmt3a2, Figure S5: Reproducibility of the replication timing changes upon Dnmt depletion, Figure S6: Pluripotent cell-specific replication timing is retained in Dnmt-depleted ES cells, Figure S7: Summary of the number of genes showing more than two-fold change within the LtoE, EtoL, EtoE, and LtoL regions defined in Figure 4A, Figure S8: Relationship between replication timing changes and transcription changes in TKO cells (the same analysis done in Figure 4C).