Emergence of replication timing during early mammalian development

Embryo collection and culture

All experiments were performed under the authorization of the authorities from Upper Bavaria (Tierversuchsantrag von Regierung von Oberbayern). The temperature, humidity and light cycle of mouse cages were maintained at 20–24 °C, 45–65% and 12/12 h dark/light, respectively. F₁ female mice (C57BL/6J × CBA) under 10 weeks of age were superovulated by intraperitoneal injection of 10 U of pregnant mare serum gonadotropin, followed by 10 U of hCG 48 h later, and were then mated with DBA/2J male mice. Zygotes were collected from the oviduct and cumulus cells removed following brief incubation in M2 medium containing hyaluronidase (Sigma-Aldrich). Zygotes were placed in drops of KSOM (potassium simplex optimized medium) and cultured at 37 °C with 5% CO₂ as previously described. For induction of parthenogenetic embryos, MII-stage oocytes were collected, as described above, from superovulated females without mating. Following removal of cumulus cells, oocytes were treated with 10 mM Sr²⁺ for 2 h in Ca²⁺-free CZB medium and then incubated in KSOM. For generation of IVF-derived zygotes, MII oocytes from F₁ female mice (C57BL/6J × CBA) were inseminated with activated spermatozoa obtained from the caudal epididymides of adult DBA/2 J male mice.

Detection of 5-ethynyl-2′-deoxyuridine incorporation

Cells were incubated with 50 μM 5-ethynyl-2′-deoxyuridine (EdU) for 1 h for each time window, as indicated, and processed for quantification of signal intensity. Incorporated EdU was visualized by Click-iT chemistry (Thermo Fisher Scientific) followed by permeabilization as described in the manufacturer’s instructions. Images were acquired on a SP8 confocal laser-scanning microscope (Leica). EdU was coupled to Alexa 594 and images acquired with a Plan-Apochromat ×63/1.4 numerical aperture 1.4 oil-immersion objective (Leica) at 561 nm excitation.

Analysis of EdU incorporation

To quantify EdU incorporation we manually cropped confocal stacks containing several embryos so that each image contained only one single embryo. Only embryos that looked fertilized and with normal pronuclei following visual inspection were included in this analysis. From embryo images we then automatically obtained the maximum intensity value in the EdU channel of the whole stack by ImageJ (v.1.53k) with a custom-made ImageJ macro. We plotted and analysed the resulting EdU intensity values for each time bin with R.

Inhibition of ZGA

For inhibition of both minor and major ZGA, embryos were treated with either 0.1 mg ml⁻¹ α-amanitin or 100 μM DRB from the zygote stage at 17 h after hCG injection until their collection for single-cell Repli-seq at the 2-cell stage. Validation of the α-amanitin effect on transcriptional silencing was done using a Click-iT RNA Alexa Fluor 594 Imaging Kit (Thermo Fisher Scientific) at the 2-cell stage (at 40 h after hCG injection).

Gene expression analyses following treatment with α-amanitin and DRB

Twelve embryos were treated with either 0.1 mg ml⁻¹ α-amanitin or 100 μM DRB from 17 to 40 h after hCG to inhibit both minor and major ZGA, then flash-frozen in liquid nitrogen in 5 μl of 2× reaction buffer (CellsDirect One-Step qRT–PCR kit, no. 11753100, Thermo Fisher). Next, 0.5 μl of a 1:200 dilution of ERCC spike-in mix (Thermo Fisher) was added to each group and TaqMan Gene Expression assays were performed according to previous work³⁸. Complementary DNA was diluted tenfold before analysis with Universal PCR Master Mix and TaqMan Gene Expression assays (Applied Biosystems). All raw C_t values were normalized by those acquired from the ERCC spike-in specific primer set, and relative expression levels of each gene were determined by the ddCt method. We assigned C_t values below the detection range as expression level 0. Primers and probes for ribosomal DNA (Hsa1) were produced by TIB MolBiol (custom design)⁴⁵. Primers and probes for Zscan4 cluster and ERCC spike-in were purchased from Applied Biosystems.

Immunostaining following either treatment by α-amanitin and DRB or expression of KDM5B

Embryos were treated with either 0.1 mg ml⁻¹ α-amanitin^55,56 or 100 μM DRB from 17 to 40 h after hCG and fixed with 4% paraformaldehyde (PFA) for 20 min at room temperature. For KDM5B expression, 2 μg μl⁻¹ KDM5B of in vitro synthesized messenger RNA was microinjected into zygotes at 18 h after hCG and fixed with 4% PFA for 20 min at room temperature at 48 h after hCG, similar to previous experiments^13,33. Embryos were then permeabilized with 0.5% Triton X-100 containing PBS for 20 min. For immunostaining following Triton pre-extraction, embryos were first permeabilized with pre-extraction buffer (50 mM NaCl, 3 mM MgCl₂, 300 mM sucrose, 25 mM HEPES, pH adjusted to 7.4) with 0.5% Triton X-100 for 10 min on ice and washed three times in pre-extraction buffer before fixing in 4% PFA at room temperature for 20 min. Following blocking for 1 h at room temperature in blocking solution (5% normal goat serum in PBS), embryos were incubated with either anti-RNA polymerase II (no. sc-899, 1:100), anti-RNA polymerase II CTD repeat YSPTSPS (phospho S2, no. ab5095, 1:1,000) or anti-H3K4me3 (Diagenode, no. C15410003, 1:250) antibody in blocking solution overnight at 4 °C. Embryos were incubated for 1.5 h at room temperature in blocking solution containing goat anti-rabbit IgG highly cross-adsorbed secondary antibody, Alexa Fluor 488 (Thermo Fisher Scientific, no. A11034, 1:1,000). After washing, embryos were mounted in Vectashield (Vector Laboratories). Confocal microscopy was performed using a ×40 oil objective on an SP8 confocal microscope (Leica) and images acquired with LAS X software.

Repli-seq

Single-cell Repli-seq was performed as previously described¹⁹ based on ref. ⁵. In brief, early-stage zygotes were collected and cultured until they reached the S phase at each developmental stage, based on their time following hCG injection. Embryos were collected at different time points at each developmental stage to achieve sampling over the entire S phase. Collection times are indicated in Supplementary Table 1. For parthenogenetic embryos and IVF-derived zygotes, the timing of S phase was calculated based on the time elapsed since activation and insemination, respectively. For KDM5B experiments, 2 μg μl⁻¹ KDM5B of in vitro synthesized mRNA was microinjected into zygotes at 18 h after hCG as previously described¹³. For each developmental stage, embryos were obtained from several litters and embryos from different litters were collected across different dates to ensure robust data collection. The number of mice used for collection of samples for each developmental stage is indicated in parentheses, as follows: zygote (20), 2-cell (30), 4-cell (27), 8-cell (20), 16-cell (15), morula (16), ICM (19), parthenotes (14), IVF zygotes (14), 2-cell + α-amanitin (14), 2-cell + DRB (24) and 2-cell + KDM5B (24). Zona pellucida was removed by exposure to acid Tyrode, and each blastomere was dissociated by gentle pipetting following trypsin treatment. For Repli-seq with physically isolated pronuclei we distinguished maternal and paternal pronuclei based on their size and relative position to the second polar body, and isolated them using micromanipulation. The remaining zygote containing a single pronucleus was also collected following removal of the polar body so that both pronuclei from the same zygote were further processed for Repli-seq. ICM cells were collected following trypsin digestion as previously described⁵⁷, with repeated oral pipetting in 0.5% trypsin and 1 mM EDTA; collection times are indicated in Supplementary Table 1. To distinguish ICM from trophectoderm cells, blastocysts were labelled with Fluoresbrite YG Microspheres (0.2 μm, Polysciences) before incubation with trypsin, and individual cells were sorted according to either positive (trophectoderm) or negative (ICM) fluorescence under a fluorescence microscope following disaggregation. Individual blastomeres or pronuclei were placed in eight-strip PCR tubes containing lysis buffer, and extracted DNA was fragmented by heat incubation. Fragmented DNA was tagged by the universal primer 5′-TGTGTTGGGTGTGTTTGGKKKKKKKKKKNN-3′ and amplified with whole-gene amplification primer sets, which have individual barcodes. This whole-genome amplification procedure was successfully used for single-cell Repli-seq in cell culture^4,5. Amplified DNA was purified using the QIAquick 96 PCR Purification Kit (QIAGEN), and concentration determined by NanoDrop (Thermo Scientific). Equal amounts of DNA from each sample (up to 96 samples) were pooled and 1 μg of each was ligated with Illumina adaptors using the NEBNext Ultra II DNA Library Prep Kit (NEB). Illumina sequences (NEBNext Multiplex Oligos for Illumina, NEB) were added to adaptor-ligated samples by PCR. Clean-up and size selection of the PCR product was done using SPRIselect (Beckman Coulter), and the quality of the library was confirmed using a 2100 Bioanalyzer with the High Sensitivity DNA Kit (Agilent).

Single-cell Repli-seq read alignment and quality control filtering

An overview of sample collection, mapping statistics and quality control is included in Supplementary Table 1. The quality control parameters we used were (1) the number of reads, which we set as 750,000 aligned reads as minimum; and (2) a coefficient of variation, which we established as a measure of equal/balanced coverage between chromosomes, thus filtering out potential cells with aneuploidy. At early stages, the reason for failure was equally the low number of reads or a high coefficient of variation (typically due to either lack of reads on a complete chromosome or in fragments of the genome; for example, zygotes 13 and 8 were excluded due to low number of reads and zygote 56 to a high coefficient of variation). At later stages, chromosome imbalances were the most common reason for failure (59 cells with high coefficient of variation versus three with low reads in the blastocyst stage), which reflects the known aneuploidy of cells at this embryonic stage. Sequencing reads were aligned to the mm10 genome using bowtie2 (v.2.3.5)⁵⁸ with the ‘–local’ option. Duplicates were marked using SAMtools (v.1.9) ‘markdup’ as described by SAMtools⁵⁹ documentation (the commands ‘fixmate’ and ‘sort samtools’ were used for this purpose accordingly). Using SAMtools view, reads were filtered by retaining only properly paired reads, removing duplicates and selecting those whose mapping quality was higher than or equal to 20. BED files of the read coordinates were generated with the BEDtools⁶⁰ (v.2.29.0) command ‘bamtobed’. Using BEDtools intersect, read counts were obtained for contiguous 50 kb genomic bins. For each cell the average of the bin counts was calculated for chromosomes 1–19; these 19 values were then next used to calculate the coefficient of variation as standard deviation divided by the mean. Cells with a coefficient of variation greater than 0.1 were removed from analyses due to chromosome imbalance. To maximize the number of samples used, the coefficient of variation was recalculated, excluding chromosomes one at a time. Cells were considered for further analysis if they passed the threshold when only one specific chromosome was removed. This chromosome was subsequently masked in downstream analyses; this filter removes abnormal genotypes and cells with aneuploidy.

Assignment of replication status

Using the read counts obtained for contiguous 50 kb genomic bins, we used the single-cell Repli-seq bioinformatic pipeline previously described⁵, which we followed with some modifications for each embryonic stage as summarized below. Window counts were first normalized to reads per million, and then each bin by its respective average of all samples within the same stage, aiming to correct for mappability biases intrinsic to genomic regions. Outlier regions were then masked, specifically the windows of the lower fifth percentile and upper first percentile values. To correct for low mappability, windows were segmented with the R package copy number (v.1.28.0, R v.4.0.0)⁶¹ to retain segments with the highest 95% of values. We did not perform the G1/G2 normalization described previously⁵, but we verified that this did not impact the results of these analyses. In brief, we used the validated mouse ES cell scRepli-seq datasets in ref. ⁵ and ran the analysis pipeline as described in their methods section with and without G1 control cells. Subsequently we compared the generated matrix of ones and zeros (that is, bins replicated and not replicated, respectively) by determining the percentage of windows that remained the same (for example, their 1 or 0 replication state did not change) after running the pipeline versus without G1 control. These analyses showed a high concordance between the two pipelines, with over 91% identity of genomic bins with zeros and ones on average across cells (Extended Data Fig. 1b). Importantly, those cells classified as outliers based on our analysis correspond to those that were removed in the original publication⁵ based on their ‘Removing outlier cells’, and were not considered for further analyses. Data were centred by the mean, scaled by the IQR for each cell and smoothed using a median filter with a running width of 15 windows, followed by segmentation with the R package copynumber. Finally, using the function normalmixEM in the R package mixtools (v.1.2.0)⁶², segmented values were used to fit a mixture model with two components to identify replicated and non-replicated window populations. To do this, two normal distribution functions were used to select a cutting threshold that better separated distributions; this value is located where the two individual normal distribution functions intersect. If no intersection was found between the means of the two normal distribution functions, the mid-point of the means was used as a threshold.

Computing replication scores, RT values and variability scores

Genome-wide replication score was defined as the percentage of replicated genomic bins for each cell. Throughout the manuscript we have used a 50 kb bin size, but we obtained similar results when using 25 and 100 kb bin size. Cells with a replication score greater than 90% and less than 10% were excluded from downstream analyses. We used the replication score to rank cells by S-phase progression for visualization of their replication status on heatmaps (Fig. 1c). Next we calculated raw RT values as the fraction of cells that replicated the given genomic bin for each stage, respectively. A RT value indicates earlier RT, because a higher proportion of cells replicated the bin. To correct for potential sampling bias of cells, we calculated the fraction of replicated cells in overlapping intervals of the genome-wide replication score with interval size of 35% and increment of 4.33% (for example, 0–35%, 4.33–39.33% and so on) for each genomic bin. The average of these 16 intervals served as the interval RT value that was used for both visualization of RT profiles (Fig. 1e) and downstream analyses. Raw and interval-averaged RT values looked similar overall (Extended Data Fig. 1c; RT raw versus interval), except for some stages in which the number of cells within replication score intervals showed a different distribution. Variability score was calculated using the following formula: score = 1 − (abs(p − 0.5)/0.5), where p is the fraction of replicated cells (ones) for the given bin; note that p is corrected for sampling (as described above). The variability score is therefore a measure of variation in the RT programme across cells, because it represents the number of cells that either replicated or did not replicate a given bin. A value of 1 means that one-half of the cells replicated a given bin and corresponds to the highest variance; likewise, a value of 0 means that either all cells replicated or did not replicate a given bin, which corresponds to the lowest variance and/or no variance.

Identification of initiation zones (referred to as RT peaks), TTRs and termination zones (referred to as RT troughs)

To distinguish the features of RT, initiation zones, TTRs and termination zones were defined based on RT values. Genomic bins were grouped into 15 clusters by their RT values using the Mclust function from the R package mclust (v.5.4.10, R v.4.1.2). Clusters were ranked by their average RT values following analysis similar to that described previously¹⁰, except that we used RT values for clustering as opposed to the 16 Repli-seq fractions. Initiation zones and termination zones were defined as consecutive bins with local maxima or minima of their cluster ranks, respectively, in sliding windows of 21 genomic bins using the rollappy function from the R package zoo (v.1.8-10). Regions between initiation zones and termination zones were defined as TTRs (Extended Data Fig. 3b). The number of initiation zones, which we refer to as RT peaks, recorded previously¹⁰ (approximtely 2,200 in neuronal progenitor cells) is similar to that reported here. To determine the significance of the changes in the number or region size of initiation zones, TTRs and termination zones throughout development, a linear model was fitted using the lm function in R (v.4.1.2). The rank of the developmental stages (that is, 1–7) served as the independent variable. The dependent variable was either the number of regions or the upper quartile of region sizes (75th percentile) for each region type. The P value of the coefficient corresponding to the slope indicates the significance of the linear trend. For composite plots, RT values were centred at the middle point of RT peak coordinates in 2 Mb windows and the median of RT values was calculated per position (Fig. 1h). To visualize relative RT compared with the neighbouring region, the minimum value of the 2 Mb window was subtracted for each stage.

Analysis of RT heterogeneity

Heterogeneity analysis was performed using the sigmoidal model formula as described previously^5,63. A sigmoidal curve was fitted for each genomic bin by the nls function from the R package stats (v.4.1.2), such that nls(y ~ 100/(1 + exp(−g × (x − M))), start = list(g = 0.1, M = m0)) (Extended Data Fig. 6a). The average genome-wide replication score of each of the 16 overlapping intervals (see above) served as the independent variable (x), with the percentage of cells that replicated the bin within the same replication score interval as dependent variable (y). Model parameters were M = mid-point, g = slope (gain) and m0 = initial value for M (100 minus the mean of y values). By this method, the replication status of the given genomic bin was related to the overall S-phase progression of cells (measured in intervals of replication score). To anchor the start and end points of the curve, 16 data points of 0 and 100 values were added to the x and y variable, respectively. Two parameters were calculated from the curve fitting, M-value and T_width. The M-value (RT mid-point, sometimes also referred to as T_rep in the literature¹⁰) is the replication score (roughly S-phase time) at which 50% of the cells replicated the given bin. A higher M-value indicates later RT. T_width is a measure of RT heterogeneity and is defined as the replication score difference (approximate S-phase time difference) of between 25 and 75% of the cells that replicated the given genomic bin. A higher T_width value indicates higher heterogeneity, because the transition from non-replicated to replicated status is greater.

Allele-specific analyses

To address any bias that could have been caused by SNPs during alignment, reads were realigned to a SNP-masked genome sequence containing an ‘N’ anywhere in which a SNP between any of the paternal (DBA) or maternal genomes (C57BL/6 × CBA) is located. The bam files were subsequently divided into paternal and maternal reads. Importantly, not all potential SNPs between strains were used. Splitting considered only SNPs that were different for the three genomes or those whose nucleotide was the same for both maternal genomes but different compared with the paternal one. Both reference preparation and splitting were performed with SNPsplit⁶⁴ (v.0.5.0). Reads were filtered using the same tools and thresholds as described above for non-allelic analyses—that is, taking into account read duplication, properly paired criteria and a mapping quality filter. Finally, as previously described, BEDtools intersect was used to count the number of reads for each contiguous 50 kb window. All subsequent analyses were performed on genomic bins, with at least five reads assigned either to the maternal or paternal genome of the same sample.

To determine allelic bias, the log₂ ratio of maternal:paternal read counts was calculated for each bin. The majority of physically separated maternal or paternal pronuclei showed a high positive (over +2) or negative (below −2) log₂ ratio, respectively. Pronuclei with a log₂ ratio of the opposite sign were exchanged for downstream analyses. We identified several parthenogenic examples among IVF zygotes (log₂ ratio above 1), which were excluded from further analyses. Finally we calculated Spearman’s correlation coefficients on log₂ maternal:paternal ratios pairwise across single zygotes and visualized these as a correlation heatmap (Extended Data Fig. 4f). A high correlation value between two zygotes indicates that, if a genomic bin has a high allelic bias in one of the zygotes it also has a high bias in the other.

Analysis of imprinted genes

Lists of maternally and paternally imprinted genes were downloaded from the Geneimprint database (https://www.geneimprint.com/site/genes-by-species.Mus+musculus). RT values were extracted for genomic bins overlapping imprinted genes. If multiple bins overlapped the same gene, RT values were averaged. For expression level and allelic bias analysis, supplementary data were downloaded from Gene Expression Omnibus (GEO) (GSE38495 and GSE45719)⁶⁵. A gene was considered expressed when its average fragments per kilobase exon per million mapped reads value in the given stage was greater than zero. Allelic bias was calculated as the log₂-transfomed ratio between read counts assigned to Cast or C57BL/6 genomes. A gene was considered maternally biased if the average log² allelic ratio was greater than zero, and paternally biased if less than zero. RT values at imprinted genes were visualized on heatmaps and ordered by their expression and allelic bias status. In total we analysed 49 maternally and 37 paternally imprinted genes, corresponding to 98 and 100 genomic bins, respectively.

Analysis of transposable elements

Transposable element annotation for the mm10 genome was obtained from Hammell’s laboratory repository (https://labshare.cshl.edu/shares/mhammelllab/www-data/TEtranscripts/TE_GTF/mm10_rmsk_TE.gtf.gz).

Enrichment of transposable elements in RT peaks, TTRs or RT troughs was estimated by calculating the log₂ ratio of the number of transposable elements of the given type overlapping with RT peaks, TTRs or RT troughs relative to the overlap of randomly shifted transposable elements with RT peaks, TTRs or RT troughs, respectively. The final enrichment value was the average of 1,000 iterations.

Statistical and genome-wide enrichment analysis

For statistical analyses of single-cell RT data we established a bootstrapping approach and calculated 95% confidence intervals to judge statistical significance⁶⁶. We chose this method to avoid the inflation of P values when n is large due to a large number of genomic bins (n = approximately 49,000) and thus we applied bootstrapping to samples, in this case single cells (n = approximately 30–70), rather than to genomic bins. Namely, we iteratively resampled individual cells with replacement 1,000 times for each stage or condition. For each iteration we recalculated RT values and any subsequent statistic—for example, Spearman’s correlation coefficient or ΔRT between conditions, as described above. We constructed confidence intervals from the bootstrap distribution using the percentile method. The 95% confidence interval is the interval between the 2.5th and 97.5th percentiles of the distribution; when 95% confidence intervals do not include zero or two intervals do not overlap, they are significantly different from zero or different from each other, respectively. For enrichment analysis of overlapping regions or gene classes, genomic bins were grouped by significantly differential RT values to increasing (earlier), decreasing (later) or non-significant (no change) bins. Enrichments were visualized on heatmaps by calculating the ratio of the observed number of overlapping bins relative to the expected value, which is the product of the row and column sums divided by the total number of bins in the corresponding contingency table.

Analysis of public chromatin datasets

Published datasets were downloaded from GEO with accession numbers GSE66581, GSE101571 (ATAC-seq³⁶), GSE71434 (H3K4me3 chromatin immunoprecipitation sequencing (ChIP)³⁴), GSE112834 (H3K36me3 ChIP⁶⁷), GSE98149 (H3K9me3 ChIP⁶⁸), GSE73952 (H3K27me3ChIP³⁹) GSE76687 (H3K27me3 ChIP⁶⁹) and GSE135457 (Pol2 Stacc-seq⁵²) andGSE76642 (DNase I hypersensitive sites sequencing⁷⁰). Paired-end reads were trimmed by cutadapt (v.3.4) with parameters -a CTGTCTCTTATA -A CTGTCTCTTATA -a AGATCGGAAGAGC -A AGATCGGAAGAGC –minimum-length=20. Following trimming, reads were aligned to the mouse reference (GRCm38) using bowtie2 (v.2.3.5) with parameters –end-to-end –very-sensitive –no-unal –no-mixed –no-discordant -I 10 -X 500. Reads were filtered by mapping quality score using SAMtools (v.1.3) with the parameter -q 12. Read pairs were read into R (v.3.6.3) using the readGAlignmentPairs function from the GenomicAlignment package (v.1.22.0) and were filtered for unique fragments. Fragments aligned to the mitochondrial genome or small scaffolds were not considered in analyses. Fragments were counted in 50 kb consecutive genomic bins (same bins as for RT profiles), normalized by the sum of fragment counts and multiplied by 1 million. Finally, normalized counts were log₂ transformed following the addition of a pseudocount of 1. Note that, for the analysis of H3K27me3 in Extended Data Fig. 10b,c the dataset used was that of Liu et al. (GSE73952)³⁹ whereas in Fig. 5f the dataset used was that of Zheng et al.⁶⁹ (GSE76687). For the correlation analysis shown in Fig. 5f we used the following stages when the actual stage was not available: early 2-cell ATAC-seq for zygote, morula DNase I hypersensitive sites sequencing for ICM and ES cell LmnB1 DamID for ICM. Differential genomic bins between conditions (for example, ATAC-seq following α-amanitin treatment) were called by DESeq2 (v.1.34.0) with an adjusted P value cutoff of 0.05. For ATAC-seq analysis in α-amanitin-treated embryos, 2-cell-stage embryos administered α-amanitin treatment by Wu et al.³⁷ (GSE101571) were compared with untreated 2-cell-stage embryos derived from Wu et al.³⁶ (GSE66581).

Analysis of public HiC and LAD datasets

HiC compartment coordinates and scores (GSE82185)¹⁶, as well as LAD coordinates (GSE112551)¹³, were analysed as previously described¹³.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.