Endogenous DNA methylation reprogramming in mammalian germline

We first segregated the ACG and TCG trinucleotides of the human genome, and calculated methylation level of the cytosine in the center position to estimate the endogenous DNA methylation levels of each sample. We found that in the early postimplantation embryos (heart tissue of 5-week embryo), the average DNA methylation level is as high as 72.9% (n = 2), whereas in the 7-week male FGCs, the methylation level decreases drastically to 22.2% (n = 2); in 12-week male FGCs, the DNA methylation level reaches its lowest point at only 6.7% on average (n = 2), confirming that germline DNA demethylation occurs within the first 12 weeks of human development, consistent with our previous results⁵. After 12 weeks, there is a small but significant increase in DNA methylation level in male FGCs, reaching to 10.7% (n = 2) at 17 weeks, and 15.1% (n = 2) at 26 weeks, indicating a global remethylation in the male germline takes place after 12 weeks of human gestation (Figure 1A and ​and1B1B).

Open in a separate window

Figure 1

Endogenous DNA methylation reprogramming in mammalian germline. (A) Boxplots showing the average endogenous DNA methylation dynamics in human germline. Endogenous DNA methylation level was calculated using the WCG sites (ACG and TCG trinucleotides). (B) Endogenous DNA methylation distribution along the gene body regions and their flanking regions in human fetal germ cells and somatic cells, with the DNA methylation level decreasing at TSS and increasing along the gene body region, and then decreasing at transcription end site (TES). (C) Boxplots showing the average endogenous DNA methylation dynamics in mouse germline. (D) Endogenous DNA methylation distribution along the gene body regions and their flanking regions in mouse PGCs and somatic cells. (E, F) The heatmap views of a representative section of chromosome 1 showing the dynamics of DNA methylome in human (E) and mouse (F) germline. Color key from white to dark blue indicates the endogenous DNA methylation level from low to high.

We also analyzed the DNA methylation dynamics in mouse PGCs in parallel. The average DNA methylation level falls from 77.7% in E6.5 epiblast to 5.0% in E11.5 PGCs, and further decreases to 4.2% and 3.8% in E13.5 male and female PGCs, respectively. This extremely hypomethylated state is maintained in E16.5 female PGCs (5.6%), whereas in E16.5 male PGCs, the DNA methylation level recovers to 61.4% (Figure 1C and ​and1D1D).

Despite drastic changes in DNA methylation in the human and mouse germline, the average DNA methylation level of neighboring somatic cells maintains at a high level (64.7%), which is consistent with previous reports^5,6,7,21,22 and verifies the accuracy of our analysis (Figure 1).

We further investigated DNA methylation dynamics and chromatin accessibility landscapes on chromosome X (chrX) to obtain a better understanding of epigenetic reprogramming of the sex chromosome, and we included one of the autosomes, chromosome 1 (chr1) as the control. Intriguingly, both chrX and chr1 clearly undergo global DNA methylation resetting during human and mouse germ cell development, with chrX having a slightly higher DNA methylation level than autosomes in all human male germ cells and most mouse PGCs, except in E16.5 male PGCs. ChrX is slightly less accessible than autosomes in both human and mouse germ cells. Higher DNA methylation level and relatively closed chromatin state suggest that epigenetic reprogramming is slightly less extensive on chrX than on autosomes during the period of germ cell development (Supplementary information, Figure S4).

Chromatin accessibility of the promoter regions in mammalian germline

To investigate comprehensively nucleosome occupancy and chromatin accessibility of the mammalian germline we applied a customized bisulfite sequencing analysis pipeline on our NOMe-seq data (Materials and Methods) and identified, in total, 116 887 and 137 557 nucleosome-depleted regions (NDRs) in human and mouse FGCs, respectively. 32 053 and 28 293 NDRs are located in the promoter regions of RefSeq genes (proximal NDRs) in human and mouse FGCs, respectively (Supplementary information, Table S2). As expected, the nucleosome-depleted signals are strongly enriched in the promoter regions, especially at the transcription start site (TSS), indicating high accessibility of the promoters of actively transcribing genes to facilitate the recruitment of transcriptional machinery and some regulatory transcription factors (TFs) to the TSS region (Figure 2A and ​and2B).2B). More importantly, downstream of the TSS, we can clearly identify at least three nucleosome positioning signals with periodical spacing (+1, +2 and +3 nucleosome position), whereas upstream of the TSS, the nucleosomes are weakly phased. This asymmetrical distribution of the nucleosome spacing surrounding the TSSs indicates that nucleosome patterning may be involved in regulating transcriptional orientation of RNA polymerase stalling and elongation. It is also possible that these regularly positioned nucleosomes immediately downstream of the TSSs are important for repressing undesired transcription initiation from alternative TSSs of the same gene^23,24,25.

Open in a separate window

Figure 2

Chromatin accessibility of the promoter regions in mammalian fetal germ cells. (A, B) Chromatin accessibility patterns of the promoter regions in human (A) and mouse (B) samples. Chromatin accessibility is calculated using the DNA methylation level of the GCH sites (GCA, GCT and GCC trinucleotides). (C) Relationship between chromatin accessibility of the promoter regions and corresponding protein-coding RefSeq gene expression in human and mouse samples. Genes are classified into four groups according to their expression level. Human fetal germ cells and mouse PGCs are presented in the middle panel as the average signals of all the human germ cells from 7 to 26 weeks of gestation and all of the mouse PGCs from E11.5 to E16.5, respectively. (D) The Spearman correlation (r1) between the gene expression level (red) and chromatin accessibility of the promoter regions (light blue), and the Spearman correlation (r2) between the gene expression level (red) and endogenous DNA methylation level of the promoter regions (green) in human fetal germ cells and somatic cells; the horizontal axis from left to right represents the genes with their expression level from low to high.

To further investigate the relationship between promoter accessibility and transcriptional activity, we grouped genes into highly expressed (RPKM > 10), intermediately expressed (1 < RPKM ≤ 10), lowly expressed (0.1 < RPKM ≤ 1) and silenced genes (RPKM ≤ 0.1) according to their expression levels. We found that promoter accessibility is highly correlated with the expression level of corresponding genes, with the promoter regions of highly expressed genes exhibiting the highest chromatin accessibility, whereas those of lowly expressed or silenced genes having lower chromatin accessibility and relatively closed chromatin states (Figure 2C and Supplementary information, Figures S1, S2). The highly expressed genes also have strongly phased nucleosomes at +1 and +2 positions downstream of the TSS, whereas genes with low expression levels are likely to have weakly phased or even evicted nucleosomes at +1 position (Figure 2C). All cells, including FGCs and somatic cells in human and mouse share similar global organization of nucleosomes in the promoter regions described above, which demonstrates evolutionary conservation of the nucleosome patterning at TSSs and its potential impact to the gene expression divergence between germ cells and somatic cells, in agreement with the previous publications^25,26,27. Besides the protein-coding RefSeq genes, these patterns are also observed in long-non-coding RNAs (lncRNAs). We quantified the expression levels of known lncRNAs in the human and mouse FGCs using Gencode datasets, and found highly expressed lncRNAs also have the more open promoter and more strongly phased nucleosomes downstream of the TSS regions, whereas lncRNAs with lower transcriptional activities tend to have closed promoters and weakly phased nucleosomes surrounding the TSS regions (Supplementary information, Figure S5A). These results collectively indicate that the openness of the promoter regions can indicate the transcriptional activity of genes (regardless of whether they have protein coding potential or not) in the mammalian germline during the period when the genome is nearly free of DNA methylation.

In parallel, we compared the endogenous CpG methylation levels in the promoter regions and the expression levels of the corresponding genes, and found, expectedly, a negative correlation in both FGCs and somatic cells at different developmental stages, including stage when the endogenous CpG methylation is reset to a very low level (Figure 2D and Supplementary information, Figure S5B and S5C).

Chromatin accessibility of the distal regulatory elements in mammalian germline

More than half of the NDRs we identified are distal NDRs (located at least 2 kb away of the annotated TSS). Different from the proximal NDRs, accessibility of the distal NDRs is weakly correlated with the gene expression levels of their neighboring genes (data not shown). Distal NDRs are symmetrically surrounded by 3-5 regularly spaced nucleosomes with intervals ranging from 160 to 280 bp (Figure 3A and ​and3B).3B). Most of the distal regulatory elements, such as enhancers, are likely to recruit tissue-specific TFs to regulate expression of their target genes^28,29,30. We took full advantage of the distal NDRs we identified to analyze the preferential enrichment patterns of TF motifs in FGCs and their surrounding somatic cells in human and mouse. We found that distal NDRs in human FGCs are more likely to be enriched for binding motifs of TFs such as NANOG, SOX17, AP2γ (also called TFAP2C) and OCT4 (also called POU5F1) (Figure 3C). This indicates that these pluripotency-specific or early germ cell-specific TFs are functionally important for human FGC development potentially by setting an open chromatin state in a large set of downstream target genes. On the other hand, distal NDRs in human gonadal somatic cells have a strong enrichment for the binding motifs of FOXA2, STAT4, NKX6-1, OLIG2 and TEAD (Figure 3C). This suggests that these TFs play crucial roles in initiating and maintaining open chromatin states for somatic cell-specific genes and consequently impact on gonadal somatic cell development. Furthermore, Ap2γ, Sox2, Oct4, Smad4, Esrrb and Klf4 binding motifs are enriched in the distal NDRs of mouse PGCs and E6.5 epiblast, whereas Gata2, Gata4 and Hoxa2 binding motifs are enriched in the distal NDRs of mouse gonadal somatic cells (Figure 3D-3F). This is compatible with the fact that Sox2 is highly expressed in early mouse PGCs and is crucial for their development, whereas SOX17 is highly expressed in early human PGCs and is important for their development^7,31,32. This also suggests that beyond the general evolutionarily conserved characteristics of the chromatin states, there are also distinct features between human and mouse germ cell development. More importantly, human FGCs from 7- to 26-week male embryos share similar TF-enrichment patterns of distal NDRs and are clustered together in an unsupervised hierarchical clustering analysis, whereas gonadal somatic cells and heart tissue at early postimplantation stage are clustered together (Figure 3C). Similarly, most of the mouse PGCs are clustered together with E6.5 epiblast, except the E16.5 female PGCs, which possess a relatively unique TF-enrichment pattern, suggesting that after entering meiotic arrest, mouse E16.5 female PGCs acquire very striking meiosis-related features in chromatin states (Figure 3D).

Open in a separate window

Figure 3

Chromatin accessibility of the distal regulatory elements in mammalian germline. (A, B) Chromatin accessibility patterns of NDRs and the flanking regions showing 4-5 symmetrically positioned nucleosomes in human (A) and mouse (B) samples. (C, D) Heatmaps showing the enrichment patterns, the corresponding relative gene expression levels and the endogenous DNA methylation levels of the promoter regions of the known transcription factors of the human (C) and mouse (D) distal NDRs. Color key from green to red indicates the enrichment from weak to strong; from dark blue to yellow indicates the gene expression level from low to high; from blue to red indicates the DNA methylation level from low to high. (E, F) Representative genome browser snapshots of the NOMe-seq signal in Pou5f1 (E) and SYCP3 (F) loci. E11.5 H3K27ac ChIP-seq data are from¹⁰. Blue bars in (E) and (F) indicate DNA methylation levels of GCH sites calculated based on NOMe-seq data sets.

Chromatin accessibility of imprinted genes and germline-specific genes

We and others previously reported that DNA methylation in the known imprinted differentially methylated regions (DMRs) is erased during the global DNA demethylation event in human and mouse FGCs^5,6,7, but little is known about chromatin accessibility of these DMRs. Our NOMe-seq results now show that unlike DNA methylation, chromatin accessibility of the imprinted DMRs is comparable before and after DNA demethylation (Figure 4A-4C).

Open in a separate window

Figure 4

Chromatin accessibility landscapes of imprinted genes and germline-specific genes. (A, B) Heatmaps showing endogenous DNA methylation (WCG DNA methylation level), chromatin accessibility of the promoter regions (GCH DNA methylation level) and relative gene expression (z-score of the FPKM values) of human (A) and mouse (B) imprinted genes and germline-specific genes. Notably, in the heatmaps, the endogenous DNA methylation levels of human and mouse imprinted genes were calculated based on the WCG DNA methylation levels of the imprinting DMRs, whereas the endogenous DNA methylation levels for the germline-specific genes were calculated based on the WCG DNA methylation levels of their promoter regions. (C) DNA methylation graphs showing endogenous DNA methylation dynamics (upper panel) and chromatin accessibility (bottom panel) of human H19 imprinted DMR across stages. In upper panel, white and open cycles indicate the unmethylated CpG sites in WCG context (endogenously unmethylated state), whereas the black and filled cycles indicate the methylated CpG sites in WCG context (endogenously methylated state). In the bottom panel, green and open cycles (the unmethylated GCH sites) indicate closed chromatin, whereas the red and filled cycles (the methylated GCH sites) indicate opened chromatin. Only the pair-ended reads with no less than four consecutive WCG or GCH trinucleotides covered are plotted.

Moreover, we found that DNA methylation levels of the promoter regions of germline-specific genes, which are highly and specifically expressed in FGCs, are significantly lower in FGCs than in the neighboring somatic cells in the gonad, whereas chromatin accessibility of most of these promoter regions is significantly higher in FGCs than in the somatic cells. Interestingly, the promoter of some germ cell-specific genes is also in an open state in 5-week heart tissue, which indicates that chromatin of some lineage-restricted genes is already made accessible at an earlier developmental stage preparing for the flexible regulation of their expression (Figures 3F and 4A-4C).

Distinct features of chromatin accessibility in mammalian germline

We next explored unique features of chromatin accessibility existing in the mammalian FGCs by systematically comparing the accessibility of distal and proximal NDRs in human and mouse FGCs and somatic cells. From the hierarchical clustering analysis and principle component analysis (PCA), we find that both the distal and proximal NDRs can accurately separate FGCs and somatic cells, with FGCs at different time points clustered together whereas somatic cells grouped together and separately from FGCs, implying that chromatin accessibility patterns of FGCs are distinct from somatic cells. Furthermore, during the DNA demethylation process, chromatin status of FGCs is kept relatively stable, compatible with the fact that during early FGC development from 4 to 11 weeks, the transcriptome of the germ cells is relatively stable despite the global demethylation of the genome^5,6,7 (Figure 5A-5D and Supplementary information, Figure S6A-S6D).

Open in a separate window

Figure 5

Distinct features of chromatin accessibility in mammalian germline. (A-D) Unsupervised hierarchical clustering analysis of chromatin accessibility of distal and proximal NDRs in human and mouse samples across replicates. (E) Multidimensional scaling (MDS) analysis of chromatin accessibility of the proximal NDRs of human-to-mouse homologous regions.

We also calculated signal intensity of the distal and proximal NDRs and only took the human-mouse homologous regions for comparison. Multidimensional scaling analysis of human and mouse proximal NDRs shows that dimension 1 (Dim.1) mainly separates these two species, while dimension 2 (Dim.2) reflects the major differences between FGCs and somatic cells, together demonstrating both species specificity and functional significance of the chromatin states of proximal NDRs for mammalian germ cell development (Figure 5E and Supplementary information, Figure S6E).

We subsequently identified several hundreds of genes with the accessibility of the promoter regions most positively correlated with their expression levels, and further categorized them into heart (or epiblast in mouse)-specific, FGC-specific and gonadal somatic cell-specific proximal NDR-associated genes. As shown in Figure 6A and ​and6B,6B, cell type-specific proximal NDR-associated genes tend to have more open promoters and higher expression levels than those in the remaining cell types. Human heart-specific proximal NDR-associated genes, such as MYL9, NKX2-5, FOXP4, MYH6, ACTC1, TNNT2 and HAND1, are significantly enriched for terms of cardiovascular system or heart development in our gene ontology (GO) analysis (Materials and Methods), whereas human FGC-specific proximal NDR-associated genes, such as SYCP3, TDRD9, DDX4, PIWIL2, BRDT and DND1, show a strong enrichment for GO terms of piRNA metabolic process, meiotic cell cycle or spermatogenesis (Figure 6C). Similarly, mouse epiblast-specific proximal NDR-associated genes are preferentially enriched for GO terms of cellular metabolic processes, whereas mouse PGC-specific proximal NDR-associated genes are enriched for reproduction, meiotic cell cycle and gamete generation-related terms (Figure 6D). These findings suggest that the accessibility of the promoter regions of cell fate determination genes contribute to the precise spatial-temporal regulation of their expression.

Open in a separate window

Figure 6

Tissue-specific proximal NDR identification. (A, B) Heatmaps showing tissue-specific proximal NDRs and the associated genes in human (A) and mouse (B) samples. (C, D) Gene ontology (GO) analysis of tissue-specific proximal NDR-associated genes in human (C) and mouse (D) samples.

Chromatin accessibility at annotated genomic elements and repetitive elements

Next, we explored the genome-wide distribution of accessible chromatin with respect to various genomic elements. We found that CpG islands (CGIs) and promoters are the most highly accessible regions. Intragenic regions are more accessible than intergenic regions, with exons more accessible than introns, and enhancers are also open as expected. These patterns are consistent in FGCs and somatic cells in both human and mouse (Figure 7A and Supplementary information, Figure S7A). Further classification of promoters into high-density CpG promoters (HCPs), intermediate-density CpG promoters (ICPs) and low-density CpG promoters (LCPs) shows that HCPs are more open than ICPs, and LCPs have relatively lower accessibility; and genes with HCP tend to have the highest expression level, whereas genes with LCP are expressed at relatively low levels (Figure 7B and Supplementary information, Figure S7B).

Open in a separate window

Figure 7

Chromatin accessibility at annotated elements and repetitive elements in human. (A) The relative enrichment analysis of chromatin accessibility at the genomic regions. (B) Chromatin accessibility of high-density CpG promoter (HCP), intermediate-density CpG promoter (ICP), low-density CpG promoter (LCP) and promoter of germline-specific genes in human germline. (C, D) Chromatin accessibility at the repetitive elements (C) and their subfamilies (D) with different evolutionary ages in human germline.

When we examined the repeat elements, we found that SINE/variable number of tandem repeats/Alu elements (SVAs) are specifically open in human germline, and also have significantly higher accessibility than SINEs, LINEs and long-terminal repeats (LTRs) (Figure 7C). SVAs also have the most abundant transcripts in human FGCs compared with neighboring somatic cells (Supplementary information, Figure S8A). In mouse, SINEs and LTRs are more accessible than LINEs, and intracisternal A particles (IAPs) are more likely to be open in mouse PGCs than in somatic cells and epiblast (Supplementary information, Figure S7C).

Furthermore, when we categorized the repetitive elements into subfamilies according to their evolutionary ages, we found a more interesting pattern. Both LINE-1 (L1) and LINE-2 (L2) belong to LINE family, with L1 evolutionarily younger than L2. And both Alu and MIR belong to SINE family, with Alu evolutionarily younger than MIR. We found that evolutionarily younger L1 is less open than L2, and evolutionarily younger Alu is less accessible than MIR in human FGCs (Figure 7D and Supplementary information, Figure S7D). In addition, the evolutionarily younger subfamilies, L1 and Alu, tend to have more active transcription activities than the evolutionarily older ones, and also retain higher levels of residual DNA methylation than the older ones after the global DNA demethylation in human and mouse FGCs, consistent with our previous reports^5,33,34 (Supplementary information, Figure S8). Taken together, these findings indicate that during the DNA methylation reprogramming process, mammalian FGCs tend to retain more residual DNA methylation and adopt a less accessible chromatin state in the evolutionarily younger and more hazardous transposable elements to, potentially, repress their transcription and transposition.

The relationship between histone modification and chromatin accessibility

Accompanying global reprogramming of the endogenous DNA methylation, the histone modifications in mammalian FGCs were also reported to be reprogrammed genome-wide^8,9,10. We analyzed the relationship between histone modifications and chromatin accessibility in mammalian germ cells at fine resolution. We incorporated the previously published histone ChIP-seq data set in mouse PGCs¹⁰, and called H3K4me3, H3K27me3 and H3K27ac peaks independently using our previously reported ChIP-seq analysis pipeline³⁵. We find that the endogenous DNA methylation levels of the genomic regions within H3K27ac, H3K4me3 and H3K27me3 peaks and 1 kb upstream/downstream flanking regions are very low in E11.5 and E13.5 mouse PGCs, which is consistent with the global hypomethylation patterns in E11.5 and E13.5 PGCs (Supplementary information, Figure S9A), while from the GCH DNA methylation levels in E11.5 and E13.5 PGCs, we can see that the H3K4me3-marked regions are more accessible than H3K27me3-marked regions at these developmental stages. When we systematically compared signal intensities of the ChIP peaks with chromatin accessibility or endogenous DNA methylation levels of the genome using non-overlapped 1 kb bins, we found that, in general, H3K4me3, H3K27ac and H3K27me3 histone modifications are clearly correlated with chromatin accessibility, with H3K4me3 having the highest correlation with chromatin accessibility than H3K27ac, and H3K27me3 being marginally correlated with chromatin accessibility. These data collectively suggest that chromatin regions with H3K4me3 modification are most likely in open states in mouse PGCs. Notably, we also found anti-correlation relation between the endogenous DNA methylation levels and H3K4me3/H3K27ac/H3K27me3 ChIP-seq signals as expected (Supplementary information, Figure S9B and S9C).

Human and mouse sample preparation

The de-identified fetus at 7-26 weeks of gestation was dissected under the microscopes, and the gonad was extensively washed with DPBS to remove any blood and other contaminants, then digested with 500 μl of Accutase Cell Detachment Solution (Millipore #SCR005) for 10-25 min at 37 °C. After digestion, the single-cell suspension was obtained by filtering through 30 μm Pre-Separation Filters (MiltenyiBiotec #130-041-407) before centrifuging at 300× g for 10 min at 4 °C. The cell pellet was washed with DPBS and finally resuspended with 300 μl L15 medium (Gibco # 21083027) with 10% FBS. Accordingly to the previous publication⁴⁵, we chose CD117 (also known as KIT) surface marker to isolate KIT-positive FGCs. Briefly, 100 μl of FcR Blocking Reagent and 100 μl of CD117 MicroBeads (MiltenyiBiotec #130-091-332) were added to the 300 μl gonad cell suspension and mixed well by gently pipetting. Then, 10 μl of PE Mouse Anti-Human CD117 antibody (BD Pharmingen #555714, clone YB5.B8) was added and mixed well before incubation at 4 °C for 30 min. After centrifuge at 300× g for 10 min at 4 °C, the cell pellet was resuspended in 500 μl L15 medium with 10% FBS, and then subjected to the MACS (MiltenyiBiotec). After MACS enrichment, the KIT-positive fraction was further sorted by FACS (BD FACSAriaII). All the human FGCs used in this study were sorted by MACS coupled with FACS, which are restricted to the KIT-positive population of the germ cells, and all the somatic cells were collected as the KIT-negative fraction after FACS. Human heart tissue was obtained from 5-week fetus with careful dissection, and all the blood was squeezed out and heart was further washed several times with DPBS before digestion with Accutase Cell Detachment Solution for 15 min at 37 °C.

Mouse E11.5, E13.5 and E16.5 PGCs were isolated from timed mated females carrying the Oct4-Gfp transgene expressed in the developing gonad with CD-1 background. For each stage, 10-30 embryos were dissected and gonads were pooled, and E13.5 and E16.5 male and female embryos were distinguished morphologically and collected separately. The pooled gonads were digested with 500 μl of Accutase Cell Detachment Solution (Millipore #SCR005) for 10-20 min at 37 °C, and the digestion reaction was quenched after adding 500 μl L15 medium (Gibco # 21083027) with 10% FBS, the cells were filtered through 30 μm Pre-Separation Filters (MiltenyiBiotec #130-041-407) and then pelleted by centrifuging at 300× g for 10 min at 4 °C. The mouse PGCs at each time point were collected after BD FACSAriaII cell sorter as GFP-positive fraction, the gonadal somatic cells were collected as the GFP-negative fraction. Mouse epiblasts were obtained at E6.5 with manually dissected under microscopy, and further separated away with the extraembryonic tissues. And single-cell suspension of mouse epiblasts was obtained after digestion with Accutase Cell Detachment Solution for 10 min at 37 °C.

All the cells isolated from the FACS cell sorter or just digested with Accutase solution were further washed with DPBS, pelleted and freezed at −80 °C for temporary storage.

Nucleosome occupancy and methylome sequencing

Approximately 1 000-10 000 cells were used for the NOMe-seq in this study. Briefly, cell pellet was first lysed in 1× lysis buffer (part of the NOMe-Seq Kit, Active motif #54000) with Protease inhibitor and PMSF (part of the NOMe-Seq Kit, Active motif #54000) on ice for 1 h with intermittent vortex. The nuclei was released, washed with DPBS, and spiked in with 0.1% unmodified lambda DNA (ThermoFisher #SD0011), and then subjected to the 15 U M.CviPI GpC Methyltransferase (NEB #M0227L) treatment for 1 h at 37 °C, supplemented with 160 nM fresh SAM (NEB #B9003S), then followed by a boost with an additional 15 U M.CviPI and 160 nM fresh SAM for another 2 h at 37 °C. The reaction was stopped by adding EDTA followed by Proteinase K digestion at 55 °C overnight, then genomic DNA was column-purified using Genomic DNA Clean and Concentrator Kit (VisTech #DC2008). For NOMe-seq, the libraries were prepared using PBAT protocol with some minor modifications^46,47. Briefly, the isolated genomic DNA was first bisulfite converted using MethylCode Bisulfite Conversion Kit (ThermoFisher #MECOV-50), and first strand was synthesized using random nonamer primers with biotin-tagged truncated Illumina P5 adapter (5′-biotin-CTACACGACGCTCTTCCGATCTNNNNNNNNN-3′), and second strands were synthesized using random nonamer primers containing a truncated P7 Illumina adapter (5′-AGACGTGTGCTCTTCCGATCTNNNNNNNNN-3′), the final library was amplified with 4-6 cycles of PCR using 1 U Kapa HiFi HS DNA Polymerase (KAPA Biosystems), and then AMpure XP Beads (Beckman Coulter) purified, quantified and pooled on Illumina HiSeq 2500 platform with 100 or 150 bp paired-end mode (sequenced by Novogene).

RNA-seq library preparation and sequencing

Total RNAs of human 5-week heart, mouse E6.5 epiblasts, mouse PGCs and corresponding somatic cells at different stages are isolated and RNA from these cells are extracted using RNeasy mini kit (Qiagen), and bulk RNA-seq libraries are constructed under the instruction of NEBNext Ultra RNA Library Prep Kit for Illumina (NEB).

For other human samples, including KIT-positive FGCs and KIT-negative gonadal somatic cells at different development stages, we first picked the single cells using mouth pipette after FACS sorting, and obtained the single-cell cDNA library using the single-cell RNA-seq protocol we developed in 2009^48,49. The amplified cDNAs were fragmented into 200-300 bp using Covaris S2 system, and the final single-cell RNA-seq library was prepared using NEBNext Ultra DNA Library Prep Kit for Illumina (NEB).

The quality-ensured RNA-seq libraries were also pooled and sequenced on HiSeq 2500 platform with 100 or 150 bp paired-end mode (sequenced by Novogene).

Data processing and bioinformatic analysis

NOMe-seq data processing

1. Data quality control

Adapter or other artificial sequence polluted reads, such as the illumina sequencing adapters or reads with more than 10% N (undetermined bases), were filtered out prior to the analysis using customized Perl scripts.

2. Sequencing alignment

Cleaned reads were mapped to the mouse or human reference genomes (mm9 or hg19 version) following our modified bisulfite sequencing pipeline⁵. PCR duplicates are removed using SAMtools⁵⁰, and the DNA methylation level of every cytosine was calculated according our previous publications^33,51. According to previous reports, GCG and CCG context should be excluded in the subsequent analysis¹⁹. Because cytosine methylation can be designated as endogenous CpG methylation or artificial M.CviPI GpC methyltransferase treatment within GCG trinucleotides, while CCG methylation can be attributed to the slight 'off-target' M.CviPI GpC methyltransferase activity. That is to say, we only focus on ACG/TCG trinucleotides (WCG, W indicates A and T nucleotides) to estimate the endogenous CpG methylation level, and GCA/GCT/GCC (GCH, H represents A, T and C nucleotides) to deduce the nucleosome position and chromatin accessibility.

3. Data reproducibility

In each human and mouse stage, we performed 2-3 replicates if materials available. The WCG and GCH methylation levels along the transcripts and their flanking regions are highly reproducible across replicates. Also, we spiked in 5% unmethylated lambda DNA to estimate the in vitro methylation efficiency of the M.CviPI, also the non-conversion rate during the bisulfite treatment. Only the samples with no less than 90% in vitro methylation efficiency, and also higher than 98% bisulfite conversion rate were retained for the downstream analysis.

4. Nucleosome-depleted regions and nucleosome position identification

Only the DNA methylation information in GCH sites was used to call NDR and nucleosome. First, the genome was sized to 100 bp sliding windows with 20 bp steps, and the C and T read counts of every GCH site with no less than 3 depths in each window were summed up and p-values (χ²-test) for the enrichment of unmethylated GCH sites of each window were calculated as the differences to the genome background. Only the significant windows with p-values passed the cutoffs (−log10(p-value) > 5) also with a minimum size of 140 bp were retained for the downstream analysis. And the NDRs were further classified as the proximal NDRs and distal NDRs depending on their distances between NDR centers and the transcriptional start sites.

Similarly, the nucleosome position was identified by using the similar strategy with 40 bp sliding windows, 20 bp steps. Only windows with p-values (χ²-test) for the enrichment of methylated GCH sites passed the cutoffs (−log10(p-value) > 3) and with a minimum size of 60 bp were kept for downstream analysis.

5. NDR-enriched motifs analysis

HOMER (version 4.7.2) was applied to identify the binding motifs and further calculate the enrichment scores for different TFs for TSS NDRs and distal NDRs with the following command:

findMotifsGenome.pl input.bed hg19 (or mm9) output_dir -size 2000 -len 8 -S 100^35,52.

6. Evolutionary comparison between human and mouse

In order to compare the endogenous DNA methylation information (WCG methylation level) and the chromatin accessibility (GCH methylation level) between human and mouse species with minimal biases, the mouse genome coordinate (mm9) was lifted over to the human assembly (hg19) using UCSC LiftOver tool (https://genome.ucsc.edu/cgi-bin/hgLiftOver), together with the genomic annotations. The endogenous DNA methylome and chromatin accessibility of mouse samples were re-analyzed using the mouse-to-human homology DNA segments, and spearman coefficients were computed and hierarchical clustering was performed using the 'hclust' function, and 'ward.D2' method with 'pairwise.complete.obs' in R software. PCA was performed using pca and prcomp package in Bioconductor^53,54.

To ensure proper transcriptome comparison between these two species, only the 16 391 well-annotated human-mouse orthologous genes downloaded from the Vertebrate Homology Database (http://www.informatics.jax.org/homology.shtml) were used for further analysis.

7. Stage-specific proximal NDR identification

Both the endogenous DNA methylation level (WCG DNA methylation level) of the promoter regions and the chromatin accessibility (GCH DNA methylation level) of these regions contribute to the regulation of corresponding gene expression level. To identify the stage-specific (or cell type-specific) proximal NDRs and also their related genes, we systematically computed the spearman correlation coefficients of endogenous DNA methylome, chromatin accessibilities and gene expression levels, and only the proximal NDRs with correlation higher than 0.6 and p-values < 0.1 were retained and further classified into heart (or epiblast), FGC and somatic cell-specific proximal NDRs, respectively. GO analysis of these proximal NDR-related genes was performed using R GOstats package⁵⁵.

RNA-seq data processing

For the RNA-seq data, the cleaned reads were aligned to the hg19 or mm9 reference genome using TopHat (version 2.0.9) with the default parameters⁵⁶. And the transcript annotations in GTF format downloaded from the genecode project (http://www.gencodegenes.org) with the version v19 for human and vM1 for mouse. The regions for protein-coding genes and lncRNAs were also defined in the GTFs. The gene expression levels (FPKM) of each sample were calculated using the cuffinks⁵⁷.

LY, JQ and FT were supported by grants from the National Natural Science Foundation of China (31230047, 31322037, 31522034, 31271543, 31571544 and 81521002) and the National Basic Research Program of China (2014CB943200 and 2012CB966704). JQ and FT were also supported by a grant from the Beijing Municipal Science and Technology Commission (D151100002415000). LY was supported by a grant from the National High Technology Research and Development Program Grant (2015AA020407). The work was supported by Beijing Advanced Innovation Center for Genomics at Peking University. We thank Mr Zhonglin Fu, Ms Xuefang Zhang and Ms Shujing Wang from the National Center for Protein Sciences Beijing (Peking University) for their excellent assistance with the FACS.