Genome-wide dynamics of replication timing revealed by in vitro models of mouse embryogenesis

Abstract

Differentiation of mouse embryonic stem cells (mESCs) is accompanied by changes in replication timing. To explore the relationship between replication timing and cell fate transitions, we constructed genome-wide replication-timing profiles of 22 independent mouse cell lines representing 10 stages of early mouse development, and transcription profiles for seven of these stages. Replication profiles were cell-type specific, with 45% of the genome exhibiting significant changes at some point during development that were generally coordinated with changes in transcription. Comparison of early and late epiblast cell culture models revealed a set of early-to-late replication switches completed at a stage equivalent to the post-implantation epiblast, prior to germ layer specification and down-regulation of key pluripotency transcription factors [POU5F1 (also known as OCT4)/NANOG/SOX2] and coinciding with the emergence of compact chromatin near the nuclear periphery. These changes were maintained in all subsequent lineages (lineage-independent) and involved a group of irreversibly down-regulated genes, at least some of which were repositioned closer to the nuclear periphery. Importantly, many genomic regions of partially reprogrammed induced pluripotent stem cells (piPSCs) failed to re-establish ESC-specific replication-timing and transcription programs. These regions were enriched for lineage-independent early-to-late changes, which in female cells included the inactive X chromosome. Together, these results constitute a comprehensive "fate map" of replication-timing changes during early mouse development. Moreover, they support a model in which a distinct set of replication domains undergoes a form of "autosomal Lyonization" in the epiblast that is difficult to reprogram and coincides with an epigenetic commitment to differentiation prior to germ layer specification.

Figure 1.

Synchronous differentiation of pluripotent ESCs to neurectoderm. (A) A 9-d neural differentiation scheme. See text for details. (B) Analyses of steady-state mRNA levels by microarrays (line charts; y-axis, arbitrary units) and reverse transcriptase (RT)-PCR (gel images) at 3-d intervals. See text for details. Actb (beta actin), loading control. (C) A heatmap showing 17,311 RefSeq genes present on the gene expression microarray, where each horizontal line represents a single gene. Log₂ transformed ratios of steady state mRNA levels in EBM3, EBM6, and EBM9 to those of ESC are displayed using the color code shown. Genes are ordered based on a K-means clustering into 11 expression clusters (Tx1–11) with distinct patterns, where gray (Tx1–2), red (Tx3–7), and green (Tx8–11) color spectra show unchanged, up-regulated and down-regulated clusters, respectively. (D) Synchronous elimination of POU5F1 from nuclei. (Left) ESC and EBM3–9 were scored for POU5F1-positive nuclei. Similar results were obtained from three experiments and their average is shown. We scored over 330 nuclei for EBM5 and EBM6, while 54–302 were scored for others. Error bars represent standard error of the mean. (Right) Representative images from EBM4, EBM5, EBM6, and EBM7 are shown. Top, POU5F1 (green); bottom, DAPI (blue); Percentages, POU5F1-positive nuclei. (E,F) RNA-FISH reveals synchronous changes in transcription of tissue-specific genes. (E) Representative images of cells with and without RNA-FISH signals (arrowheads). (F) RNA-FISH analysis at 3-d intervals. One representative experiment for each gene is shown. At least 147 nuclei were scored for each state. Similar results were obtained in at least two replicate differentiation series. Note that even highly expressed genes are not 100% positive due to the probabilistic nature of interactions with the transcription machinery (Mitchell and Fraser 2008).

Figure 2.

Kinetics of replication-timing changes during neural differentiation of ESCs. (A) Replication-timing profiling. Exemplary profile of EBM3 is shown. (Gray dots) Probe log ratios [=log₂(Early/Late)] along chromosome 16. A local polynomial smoothing (loess) curve is overlaid (blue). Replication domains (red lines) and their boundaries (dotted lines) identified by segmentation are also overlaid (Hiratani et al. 2008). (B) Pearson's R²-values for pairwise comparisons of replication timing. Smoothed data for all probes were used for the calculation. (C,D) Overlaid replication-timing profiles of ESC, EBM3, EBM6, and EBM9, at selected EtoL domains (C) and LtoE domains (D), using the color code at the bottom. Representative genes within the domains are shown with their chromosomal positions in red squares. (E) Pearson's R²-values are shown for the relationship between average replication-timing ratios of replication domains vs. their %GC and %LINE-1 values in ESC, EBM3, EBM6, and EBM9. (F) A heatmap showing replication-timing ratios [=log₂(Early/Late)] of 18,679 RefSeq genes based on the microarrays using the color code shown. Genes are ordered based on a K-means clustering into 20 replication-timing clusters with distinct patterns. Orange, blue, red, and gray boxes represent EtoE (RT1–8), LtoE (RT9–12), EtoL (RT13–16), and LtoL (RT17–20) clusters, respectively (kinetics in G,H and Supplemental Fig. 4). (G,H) Magnified views of LtoE (G, RT9–12) and EtoL (H, RT13–16) clusters, respectively. Line charts show mean replication-timing ratios of each cluster.

Figure 3.

Subnuclear repositioning associated with replication-timing changes occurs during the EBM3–EBM6 transition, in parallel with chromatin fiber reorganization. (A–C) Analysis of subnuclear positions of eight genomic regions by 3D DNA-FISH in EBM3 and EBM6. ESC and EBM9 data (Hiratani et al. 2008) are also shown for comparison. Box plots show the distribution of relative radial distance to the nuclear periphery, where 0 and 1 represents the periphery and the center of the nucleus, respectively. Horizontal bars represent the 10th, 25th, 50th (median), 75th, and 90th percentiles. P-values were obtained from a two-sample Kolmogorov-Smirnov test. Below the box plots are the overlays of the four replication-timing profiles (using the color code at the top) with the probe gene positions (red squares). All EtoL domains (A; Zfp42, Rex2, and Dppa2) and LtoE domains (B; Ptn, Ephb1, and Akt3) analyzed make the most significant movement toward and away from the nuclear periphery, respectively, during the EBM3–EBM6 transition. During this period, replication-timing changes traverse the mid-late S-phase (pink shades), which is when a dramatic interior to peripheral transition in the spatial patterns of DNA replication foci occurs (Hiratani et al. 2009). In contrast, two EtoE domains (C; Pou5f1 and Nanog) do not change subnuclear positioning or replication time. Comparable results were obtained from two to four biological replicates, and the sum of all experiments is shown. Seventy-one to 223 FISH signals were measured per state. The identical data sets are also displayed using cumulative frequency plots (Supplemental Fig. 5). (D) Representative FISH images of Dppa2. Dotted lines represent the rim of nuclear DAPI signals. Arrowheads point to green FISH signals. (E) Electron spectroscopic imaging (ESI) analysis of nuclei from ESC, EBM3, EBM6, and EBM9. Relative levels of phosphorus and nitrogen levels were used to delineate chromatin (yellow) vs. protein and ribonucleoprotein (blue) (Bazett-Jones et al. 2008).

Figure 4.

Relationship between replication-timing changes, transcription changes, and their reversibility (A) Cross-tabulation of K-means clusters of replication timing and transcription. All 17,311 RefSeq genes present on the gene expression array were analyzed. Each of the 20 replication-timing clusters, RT1–20 (Fig. 2F), were presented as a stacked bar consisting of the 11 expression clusters, Tx1–11 (Fig. 1C). N, number of genes in RT1–20. (B) Line chart showing the kinetics of transcription changes of Tx1–11 during differentiation. Log₂ transformed ratios to ESC levels are plotted. (C,D) The same cross-tabulation graph as in A, for 10,586 HCP genes (C) and 2650 LCP genes (D). HCP and LCP classification is based on Mikkelsen et al. (2007). (E) A scheme for the reversal of differentiation. EBM6R is a condition in which EBM6 cells were placed back into ESC medium containing LIF and cultured for 6 d. (F,G) Box plots showing the reversibility of different classes of genes. Log₂ transformed expression level in EBM6R relative to EBM6 is plotted. Twofold down-regulated genes during the ESC–EBM6 transition were selected from the down-regulated expression clusters, Tx8–11. Among them, genes in the EtoL clusters (RT13–16) are presented in F, while genes in the EtoE clusters (RT1–8) are presented in G. P-values were obtained from a two-tailed t-test for comparison of two paired groups (i.e., EBM6R vs. EBM6). (H) LCP is overrepresented in EtoL genes that are not reversed in EBM6R. Twofold down-regulated EtoL and EtoE genes that did not show reversal (i.e., EBM6R/EBM6 < 1) were calculated for their LCP-to-HCP ratio. Note that there are four times as many HCPs as LCPs in the genome (=All genes). Numbers of genes in each category are 49 (EtoL) and 313 (EtoE). A P-value was obtained from a χ² test.

Figure 5.

Lineage-independent replication-timing changes are completed at the EpiSC stage, coincident with spatial repositioning of early-to-late domains. (A) Pearson's R²-values for pairwise comparisons of gene expression microarray profiles of 17,311 RefSeq genes. EpiSC shows the highest correlation to EBM3. (B) Immunostaining confirms POU5F1 expression in EpiSC5 (89% positive; 88% for EpiSC7 [data not shown]). The identical EpiSC batches were used for replication analysis. ESCs and NPCs are positive and negative controls, respectively. (C) Hierarchical clustering of 12 cell lines. The whole genome was divided into 10,974 ∼200-kb segments and their average replication-timing ratios were compared between cell lines. (Left) A whole-genome heatmap and the percentages of segments that stay early (EtoE), stay late (LtoL), shift earlier (LtoE), and shift later (EtoL), based on a K-means clustering (K = 20). LtoE and EtoL clusters were defined as those that show a replication-timing differential of above 1.0 between any cell types. With this stringent cutoff, 15% of the genome exhibited differences (as opposed to 20% with a less stringent cutoff; Hiratani et al. 2008). Clusters 9–14 correspond to LtoE and EtoL segments used for hierarchical clustering on the right. The dendrogram draws a clear distinction between EPL/EBM3 and EpiSCs, which corresponds to early and late epiblast stages, respectively. Cell lines not described in this study are from Hiratani et al. (2008). (D–I) Two-color 2D DNA-FISH in ESCs, EPL cells, and EpiSCs (EpiSC7 line). (D) Pou5f1, which is early replicating in all three cell lines, maintains its internal positioning. Box plots show the distribution of relative radial distance to the periphery, where 0 and 1 represents the periphery and the center of the nucleus, respectively. Relative radial distance medians were 0.40 (ESC), 0.38 (EPL), and 0.42 (EpiSC). N = 222–430. (E–G) Subnuclear positioning of EtoL loci (Zfp42, Rex2, and Dppa2) relative to Pou5f1. Radial distance of Pou5f1 loci to the nuclear periphery was subtracted from that of Zfp42 (E), Rex2 (F), or Dppa2 (G) loci for all four combinations of allele pairs within a given nucleus, divided by the radius and their distribution was plotted as box plots. Most peripheral, identical, and most internal positioning of these loci relative to Pou5f1 are represented by –1, 0 and 1, respectively. P-values were obtained from a two-sample Kolmogorov-Smirnov test. All three EtoL loci exhibited significant repositioning toward the periphery in EpiSCs, but not others. N = 120–384. The identical data sets are also shown using cumulative frequency plots (Supplemental Fig. 10). (H) Representative images of 2D DNA-FISH. (I) Summary of 2D DNA-FISH. Positions shown for each locus (red, EtoL loci; green, Pou5f1) represent the median relative radial distance to the nuclear periphery, showing repositioning of EtoL loci toward the periphery in EpiSCs, but not EPL cells. In EPL cells, median values were: Zfp42 (0.40), Rex2 (0.30), and Dppa2 (0.42). In EpiSCs, median values were: Zfp42 (0.23), Rex2 (0.11), and Dppa2 (0.27).

Figure 6.

Replication-timing changes across different lineages collectively affect nearly half the genome and include a set of lineage-independent early-to-late changes. (A) Hierarchical clustering with the addition of Gsc+Sox17– mesoderm, Gsc+Sox17+ endoderm, embryonic fibroblasts (f, female; m, male), and fetal myoblast cells. Their characterization is provided in Supplemental Figures 11 and 12. K-means clustering was applied first (K = 20; these K-means clusters are different from those defined in Fig. 5C) and for hierarchical clustering, K-means clusters 4–17 (framed in blue) were used, which showed a differential of >0.80 between any cell types. (B) Centroids (a set of average replication-timing ratios for a given cluster) of K-means clusters 1–20 in A presented in a heatmap format, which shows the average ratios for six cell states. Asterisks indicate lineage-independent EtoL clusters. (C) Sequence properties of K-means clusters 1–20 in A and B. Based on the average replication-timing ratios, we categorized clusters into those that stay early (EtoE), late (LtoL), middle-early (ME), middle-late (ML), shift earlier (LtoE), and shift later (EtoL). LtoE and EtoL clusters were defined using a stringent cutoff (those with a differential of >1.0). (D) Correlation between replication timing and GC/LINE-1 content in different cell states. Clusters 4–17 were used for calculation of Pearson's R²-values.

Figure 7.

Lineage-independent early-to-late replication-timing changes are difficult to reprogram. (A) Hierarchical clustering of 22 cell lines, with the addition of piPSCs (1A2, 1B3, V3) and iPSCs (1D4, 2D4) (Maherali et al. 2007), using K-means clusters 4–17 defined in Figure 6A (framed in blue). (B) Properties of seven K-means clusters (identical clusters as in Fig. 6A–C) that showed large differential (>0.95) between MEFs and ICM. The fifth column (% Reversed Segments) shows the difficulty of regaining ESC-specific replication timing in piPSCs derived from MEFs. The numbers represent the percentage of 200-kb segments that showed more than 75% recovery in replication timing in piPSCs. (C) Average replication-timing differential of K-means clusters 1–20 relative to ICM in a heatmap format (identical clusters as in Fig. 6A–C). Here, red and green represent earlier and later shifts, respectively. (Asterisks) Clusters analyzed in B. (D) Pearson's R²-values for pairwise comparisons of gene expression in iPSCs, piPSCs, and MEFs to ESCs within each of the 20 clusters. At the bottom of the table is the whole genome comparison. Note that clusters 15 and 16 deviate considerably from the rest of the genome when piPSCs and ESCs are compared, indicating that gene expression program in these clusters is particularly resistant to reprogramming, in contrast to other clusters and the whole genome comparison. Gene expression levels are based on the work of Sridharan et al. (2009) (Supplemental Table 2).

Figure 8.

Analysis of late replicating inactive X chromosome. (A) Heatmap representation of replication timing of all 200-kb segments along the X chromosome in 22 cell lines. Below are the chromosome average replication-timing ratios. Note that piPSCs and female MEFs show markedly later replication timing (**). Female EpiSC7 is also later replicating (*). Female cell lines include iPSC1D4, iPSC2D4, EpiSC7, piPSC1A2, piPSC1B3, piPSCV3, and MEF (f); all others are male lines. (B–G) Pairwise comparisons of smoothed replication-timing profiles of the X chromosome. Note the difference between cell types, particularly at the early replicating peaks along the chromosome. (Note: Coordinates 23–32 Mb exhibit large gaps in probes due to highly repetitive sequences.) (H) H3K27me3 staining of EpiSC7. Forty-four percent of the cell population showed Barr body-like nuclear staining (arrowheads) of H3K27me3.

Figure 9.

A model: Lineage-independent and lineage-dependent replication-timing changes before and after the post-implantation epiblast stage, respectively. Schematic diagrams of mouse embryos at different stages of embryogenesis are shown along with cell types that model different tissues and stages analyzed in this study (red). We propose that lineage-independent EtoL changes are completed by the post-implantation epiblast stage (E5.0–E6.0), prior to germ layer specification and down-regulation of key pluripotency transcription factors (POU5F1/NANOG/SOX2). These changes are accompanied by several events listed, which collectively represent a form of genome reorganization (“autosomal Lyonization”) that is coincident with EpiSCs having lost the ability to easily revert back to the ESC state. Lineage-dependent EtoL and LtoE changes occur continuously after the late epiblast stage to create cell-type-specific profiles. piPSCs are trapped at an epigenetic state that has not yet reprogrammed this epiblast stage genome reorganization.