Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

  • Article
  • Published:

Epigenetic modulation of a hardwired 3D chromatin landscape in two naive states of pluripotency

A Publisher Correction to this article was published on 16 May 2019

This article has been updated

Abstract

The mechanisms underlying enhancer activation and the extent to which enhancer–promoter rewiring contributes to spatiotemporal gene expression are not well understood. Using integrative and time-resolved analyses we show that the extensive transcriptome and epigenome resetting during the conversion between ‘serum’ and ‘2i’ states of mouse embryonic stem cells (ESCs) takes place with minimal enhancer–promoter rewiring that becomes more evident in primed-state pluripotency. Instead, differential gene expression is strongly linked to enhancer activation via H3K27ac. Conditional depletion of transcription factors and allele-specific enhancer analysis reveal an essential role for Esrrb in H3K27 acetylation and activation of 2i-specific enhancers. Restoration of a polymorphic ESRRB motif using CRISPR–Cas9 in a hybrid ESC line restores ESRRB binding and enhancer H3K27ac in an allele-specific manner but has no effect on chromatin interactions. Our study shows that enhancer activation in serum- and 2i-ESCs is largely driven by transcription factor binding and epigenetic marking in a hardwired network of chromatin interactions.

This is a preview of subscription content, access via your institution

Access options

Buy this article

Prices may be subject to local taxes which are calculated during checkout

Fig. 1: The transcriptional dynamics and the chromatin interactions in serum- and 2i-ESCs.
Fig. 2: Differential gene expression and enhancer–promoter rewiring in different states of ESCs.
Fig. 3: Epigenetic marking of enhancers and the association with differential gene expression.
Fig. 4: Time-resolved enhancer activation during serum-to-2i ESCs transition.
Fig. 5: Motif enrichment and TF occupancy at enhancers.
Fig. 6: ESRRB is associated with enhancer H3K27ac in the 2i state.
Fig. 7: Allele-specific enhancer H3K27 acetylation.

Similar content being viewed by others

Data availability

RNA-, ChIP-, ATAC- and 4C-sequencing data that support the findings of this study have been deposited in the Gene Expression Omnibus (GEO) under the accession code GSE92412. Previously published sequencing data that were reanalysed here are available under the accession codes GSE23943 (ref. 16), GSE72164 (ref. 18), GSE56138 (ref. 3), GSE56312 (ref. 22) and GSE110164 (ref. 15). Source data for Figs. 17 and Supplementary Figs. 17 have been provided as Supplementary Table 6. All other data supporting the findings of this study are available from the corresponding author on reasonable request.

Code availability

All computational codes used in this study are available on the Bitbucket repository: https://bitbucket.org/wmegchel/ncb_2019.

Change history

  • 16 May 2019

    In the version of the article originally published, extra lines were displayed in Fig. 7. Fig. 7a contained a solid black line that extended into panel b, and Fig. 7c contained two extra scale bars on the left. These have been removed from the figure. The errors have been corrected in the HTML and PDF versions of the article.

References

  1. Long, H. K., Prescott, S. L. & Wysocka, J. Ever-changing landscapes: transcriptional enhancers in development and evolution. Cell 167, 1170–1187 (2016).

    Article  CAS  Google Scholar 

  2. Li, W., Notani, D. & Rosenfeld, M. G. Enhancers as non-coding RNA transcription units: recent insights and future perspectives. Nat. Rev. Genet. 17, 207–223 (2016).

    Article  CAS  Google Scholar 

  3. Buecker, C. et al. Reorganization of enhancer patterns in transition from naive to primed pluripotency. Cell Stem Cell 14, 838–853 (2014).

    Article  CAS  Google Scholar 

  4. Factor, D. C. et al. Epigenomic comparison reveals activation of “seed” enhancers during transition from naive to primed pluripotency. Cell Stem Cell 14, 854–863 (2014).

    Article  CAS  Google Scholar 

  5. Ghavi-Helm, Y. et al. Enhancer loops appear stable during development and are associated with paused polymerase. Nature 512, 96–100 (2014).

    Article  CAS  Google Scholar 

  6. Comoglio, F. et al. Thrombopoietin signaling to chromatin elicits rapid and pervasive epigenome remodeling within poised chromatin architectures. Genome Res. 28, 295–309 (2018).

    Article  CAS  Google Scholar 

  7. Ji, X. et al. 3D chromosome regulatory landscape of human pluripotent cells. Cell Stem Cell 18, 262–275 (2016).

    Article  CAS  Google Scholar 

  8. Jin, F. et al. A high-resolution map of the three-dimensional chromatin interactome in human cells. Nature 503, 290–294 (2013).

    Article  CAS  Google Scholar 

  9. Palstra, R. J. et al. The beta-globin nuclear compartment in development and erythroid differentiation. Nat. Genet. 35, 190–194 (2003).

    Article  CAS  Google Scholar 

  10. Phanstiel, D. H. et al. Static and dynamic DNA loops form AP-1-bound activation hubs during macrophage development. Mol. Cell 67, 1037–1048.E6 (2017).

    Article  CAS  Google Scholar 

  11. Weinberger, L., Ayyash, M., Novershtern, N. & Hanna, J. H. Dynamic stem cell states: naive to primed pluripotency in rodents and humans. Nat. Rev. Mol. Cell Biol. 17, 155–169 (2016).

    Article  CAS  Google Scholar 

  12. Atlasi, Y. & Stunnenberg, H. G. The interplay of epigenetic marks during stem cell differentiation and development. Nat. Rev. Genet. 18, 643–658 (2017).

    Article  CAS  Google Scholar 

  13. Ying, Q. L. et al. The ground state of embryonic stem cell self-renewal. Nature 453, 519–523 (2008).

    Article  CAS  Google Scholar 

  14. Boroviak, T., Loos, R., Bertone, P., Smith, A. & Nichols, J. The ability of inner-cell-mass cells to self-renew as embryonic stem cells is acquired following epiblast specification. Nat. Cell Biol. 16, 516–528 (2014).

    Article  CAS  Google Scholar 

  15. Senft, A. D. et al. Combinatorial Smad2/3 activities downstream of nodal signaling maintain embryonic/extra-embryonic cell identities during lineage priming. Cell Rep. 24, 1977–1985 (2018).

    Article  CAS  Google Scholar 

  16. Marks, H. et al. The transcriptional and epigenomic foundations of ground state pluripotency. Cell 149, 590–604 (2012).

    Article  CAS  Google Scholar 

  17. Habibi, E. et al. Whole-genome bisulfite sequencing of two distinct interconvertible DNA methylomes of mouse embryonic stem cells. Cell Stem Cell 13, 360–369 (2013).

    Article  CAS  Google Scholar 

  18. Joshi, O. et al. Dynamic reorganization of extremely long-range promoter-promoter interactions between two states of pluripotency. Cell Stem Cell 17, 748–757 (2015).

    Article  CAS  Google Scholar 

  19. Ter Huurne, M., Chappell, J., Dalton, S. & Stunnenberg, H. G. Distinct cell-cycle control in two different states of mouse pluripotency. Cell Stem Cell 21, 449–455 e444 (2017).

    Article  Google Scholar 

  20. Choi, J. et al. Prolonged Mek1/2 suppression impairs the developmental potential of embryonic stem cells. Nature 548, 219–223 (2017).

    Article  CAS  Google Scholar 

  21. Tosolini, M. et al. Contrasting epigenetic states of heterochromatin in the different types of mouse pluripotent stem cells. Sci. Rep. 8, 5776 (2018).

    Article  Google Scholar 

  22. Galonska, C., Ziller, M. J., Karnik, R. & Meissner, A. Ground state conditions induce rapid reorganization of core pluripotency factor binding before global epigenetic reprogramming. Cell Stem Cell 17, 462–470 (2015).

    Article  CAS  Google Scholar 

  23. Schoenfelder, S. et al. Polycomb repressive complex PRC1 spatially constrains the mouse embryonic stem cell genome. Nat. Genet. 47, 1179–1186 (2015).

    Article  CAS  Google Scholar 

  24. Sahlen, P. et al. Genome-wide mapping of promoter-anchored interactions with close to single-enhancer resolution. Genome Biol. 16, 156 (2015).

    Article  Google Scholar 

  25. Domcke, S. et al. Competition between DNA methylation and transcription factors determines binding of NRF1. Nature 528, 575–579 (2015).

    Article  CAS  Google Scholar 

  26. Hackett, J. A. & Surani, M. A. Regulatory principles of pluripotency: from the ground state up. Cell Stem Cell 15, 416–430 (2014).

    Article  CAS  Google Scholar 

  27. Dixon, J. R. et al. Chromatin architecture reorganization during stem cell differentiation. Nature 518, 331–336 (2015).

    Article  CAS  Google Scholar 

  28. Rao, S. S. et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell 159, 1665–1680 (2014).

    Article  CAS  Google Scholar 

  29. Dixon, J. R., Gorkin, D. U. & Ren, B. Chromatin domains: the unit of chromosome organization. Mol. Cell 62, 668–680 (2016).

    Article  CAS  Google Scholar 

  30. Apostolou, E. et al. Genome-wide chromatin interactions of the Nanog locus in pluripotency, differentiation, and reprogramming. Cell Stem Cell 12, 699–712 (2013).

    Article  CAS  Google Scholar 

  31. Wei, Z. et al. Klf4 organizes long-range chromosomal interactions with the oct4 locus in reprogramming and pluripotency. Cell Stem Cell 13, 36–47 (2013).

    Article  CAS  Google Scholar 

  32. de Wit, E. et al. The pluripotent genome in three dimensions is shaped around pluripotency factors. Nature 501, 227–231 (2013).

    Article  Google Scholar 

  33. Denholtz, M. et al. Long-range chromatin contacts in embryonic stem cells reveal a role for pluripotency factors and polycomb proteins in genome organization. Cell Stem Cell 13, 602–616 (2013).

    Article  CAS  Google Scholar 

  34. Krijger, P. H. et al. Cell-of-origin-specific 3D genome structure acquired during somatic cell reprogramming. Cell Stem Cell 18, 597–610 (2016).

    Article  CAS  Google Scholar 

  35. Beagan, J. A. et al. Local genome topology can exhibit an incompletely rewired 3D-folding state during somatic cell reprogramming. Cell Stem Cell 18, 611–624 (2016).

    Article  CAS  Google Scholar 

  36. Dowen, J. M. et al. Control of cell identity genes occurs in insulated neighborhoods in mammalian chromosomes. Cell 159, 374–387 (2014).

    Article  CAS  Google Scholar 

  37. Cairns, J. et al. CHiCAGO: robust detection of DNA looping interactions in Capture Hi-C data. Genome Biol. 17, 127 (2016).

    Article  Google Scholar 

  38. Novo, C. L. et al. Long-range enhancer interactions are prevalent in mouse embryonic stem cells and are reorganized upon pluripotent state transition. Cell Rep. 22, 2615–2627 (2018).

    Article  CAS  Google Scholar 

  39. Loh, Y. H. et al. The Oct4 and Nanog transcription network regulates pluripotency in mouse embryonic stem cells. Nat. Genet. 38, 431–440 (2006).

    Article  CAS  Google Scholar 

  40. Dunn, S. J., Martello, G., Yordanov, B., Emmott, S. & Smith, A. G. Defining an essential transcription factor program for naive pluripotency. Science 344, 1156–1160 (2014).

    Article  CAS  Google Scholar 

  41. Adachi, K. et al. Esrrb unlocks silenced enhancers for reprogramming to naive pluripotency. Cell Stem Cell 23, 266–275 e266 (2018).

    Article  CAS  Google Scholar 

  42. Luikenhuis, S., Wutz, A. & Jaenisch, R. Antisense transcription through the Xist locus mediates Tsix function in embryonic stem cells. Mol. Cell Biol. 21, 8512–8520 (2001).

    Article  CAS  Google Scholar 

  43. Panning, B., Dausman, J. & Jaenisch, R. X chromosome inactivation is mediated by Xist RNA stabilization. Cell 90, 907–916 (1997).

    Article  CAS  Google Scholar 

  44. Furlong, E. E. M. & Levine, M. Developmental enhancers and chromosome topology. Science 361, 1341–1345 (2018).

    Article  CAS  Google Scholar 

  45. Hnisz, D., Day, D. S. & Young, R. A. Insulated neighborhoods: structural and functional units of mammalian gene control. Cell 167, 1188–1200 (2016).

    Article  CAS  Google Scholar 

  46. Cho, W. K. et al. Mediator and RNA polymerase II clusters associate in transcription-dependent condensates. Science 361, 412–415 (2018).

    Article  CAS  Google Scholar 

  47. Gu, B. et al. Transcription-coupled changes in nuclear mobility of mammalian cis-regulatory elements. Science 359, 1050–1055 (2018).

    Article  CAS  Google Scholar 

  48. Stadhouders, R. et al. Transcription factors orchestrate dynamic interplay between genome topology and gene regulation during cell reprogramming. Nat. Genet. 50, 238–249 (2018).

    Article  CAS  Google Scholar 

  49. Ostuni, R. et al. Latent enhancers activated by stimulation in differentiated cells. Cell 152, 157–171 (2013).

    Article  CAS  Google Scholar 

  50. Parkhomchuk, D. et al. Transcriptome analysis by strand-specific sequencing of complementary DNA. Nucleic Acids Res. 37, e123 (2009).

    Article  Google Scholar 

  51. Langmead, B., Trapnell, C., Pop, M. & Salzberg, S. L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25 (2009).

    Article  Google Scholar 

  52. Wu, T. D. & Watanabe, C. K. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics 21, 1859–1875 (2005).

    Article  CAS  Google Scholar 

  53. Anders, S., Pyl, P. T. & Huber, W. HTSeq—A Python framework to work with high-throughput sequencing data. Bioinformatics 31, 166–169 (2015).

    Article  CAS  Google Scholar 

  54. Anders, S. & Huber, W. Differential expression analysis for sequence count data. Genome Biol. 11, R106 (2010).

    Article  CAS  Google Scholar 

  55. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).

    Article  Google Scholar 

  56. Ernst, J. & Kellis, M. ChromHMM: automating chromatin-state discovery and characterization. Nat. Methods 9, 215–216 (2012).

    Article  CAS  Google Scholar 

  57. Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y. & Greenleaf, W. J. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat. Methods 10, 1213–1218 (2013).

    Article  CAS  Google Scholar 

  58. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, R137 (2008).

    Article  Google Scholar 

  59. Georgiou, G. & van Heeringen, S. J. fluff: exploratory analysis and visualization of high-throughput sequencing data. PeerJ 4, e2209 (2016).

    Article  Google Scholar 

  60. Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 38, 576–589 (2010).

    Article  CAS  Google Scholar 

  61. Gupta, S., Stamatoyannopoulos, J. A., Bailey, T. L. & Noble, W. S. Quantifying similarity between motifs. Genome Biol. 8, R24 (2007).

    Article  Google Scholar 

  62. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate - a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289–300 (1995).

    Google Scholar 

  63. Poser, I. et al. BAC TransgeneOmics: a high-throughput method for exploration of protein function in mammals. Nat. Methods 5, 409–415 (2008).

    Article  CAS  Google Scholar 

  64. Marks, H. et al. Dynamics of gene silencing during X inactivation using allele-specific RNA-seq. Genome Biol. 16, 149 (2015).

    Article  Google Scholar 

  65. Splinter, E., de Wit, E., van de Werken, H. J., Klous, P. & de Laat, W. Determining long-range chromatin interactions for selected genomic sites using 4C-seq technology: from fixation to computation. Methods 58, 221–230 (2012).

    Article  CAS  Google Scholar 

  66. Nagano, T. et al. Comparison of Hi-C results using in-solution versus in-nucleus ligation. Genome Biol. 16, 175 (2015).

    Article  Google Scholar 

  67. Klein, F. A. et al. FourCSeq: analysis of 4 C sequencing data. Bioinformatics 31, 3085–3091 (2015).

    Article  CAS  Google Scholar 

  68. Raviram, R. et al. 4C-ker: a method to reproducibly identify genome-wide interactions captured by 4C-Seq experiments. PLoS Comput. Biol. 12, e1004780 (2016).

    Article  Google Scholar 

Download references

Acknowledgements

We thank S. van Heeringen, G. Georgiou and A. Brinkman for assistance with bioinformatics analysis; and E. Janssen-Megens, K. Berentsen and N. Sharifi for sequencing. The computational work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. This research was supported and funded by the European Union grant ERC-2013-ADG-339431 ‘SysStemCell’.

Author information

Authors and Affiliations

Authors

Contributions

Y.A., W.M. and H.G.S. conceived the paper. Y.A. performed the experiments. Y.A., W.M. and H.G.S. analysed the data. T.P. generated the ATAC-seq data. C.W. and C.L. assisted with the 4C-seq. I.P. generated the BAC DNA, O.J. and S.-Y.W. generated the cHiC data, E.H. assisted with bioinformatics analysis and H.M. assisted with SNP analysis. Y.A., W.M. and H.G.S. wrote the paper.

Corresponding author

Correspondence to Hendrik G. Stunnenberg.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Integrated supplementary information

Supplementary Fig. 1 cHiC in serum- and 2i-ESCs.

a) Gene ontology analysis of genes assigned to PC1 (n=945 genes) and PC2 (n=599 genes). b) Robust K-means clustering (PAM) of differentially expressed genes in the serum-to-2i and the reverse transition (left part). Heatmap showing enrichment of genes assigned to each K-means cluster in different stages of pre-implantation embryos (right part). N number indicated in the panel represents genes. c) Schematic representation of the cHiC analysis for the 340,000 pairwise interactions. Interactions were constrained to have a regulatory element (ATAC-seq peak or CTCF binding site) within 2 kb of both interaction anchors. The interaction network, constrained to interacting regulatory elements, contained ~135,000 distinct pairwise interactions. d) The promoter interaction intensity (sum of reads of all interactions with that promoter) in our cHiC data is in good agreement with published promoter capture in serum ESCs, n=11,842 promoters, R-values represent Pearson’s correlation coefficient. e) ChromHMM-based segmentation of the mouse genome using six chromatin states revealed CTCF binding sites, weak- and strong regulatory elements, H3K27me3-marked regions and actively transcribed gene bodies. The other state has none of the studied marks. f) Assignment of each interacting anchor site to corresponding ChromHMM states. Sites within 2 kb of a TSS were assigned as promoters. Promoter regions with H3K27me3 are colored in yellow. Strong and weak distal elements were combined and assigned as enhancers. g) Pairwise combinations of anchors assigned to one of the main chromatin states: CTCF, weak-element, strong-element or H3K27me3 region located in promoter or distal sites. Circle sizes correspond to numbers of interactions and colors represent the interaction intensities. In a-b FDR was computed using Hypergeometric tests with Benjamini-Hochberg correction.

Supplementary Fig. 2 Characterizing the chromatin interactions in serum- and 2i-ESCs.

a) The percentage of differential chromatin interactions (y-axis) detected per category using different fold change cutoffs (x-axis). Decreasing the fold-change cutoff does not significantly increase the number of differential interactions for non-H3K27me3-associated loci. b) Insulated neighborhoods are highly similar between serum- and 2i-ESCs. The cumulative sum of all cHiC interaction-reads covered within an insulated neighborhood is shown. N number indicated in the panel represents insulated neighborhoods. R-values represent Pearson’s correlation coefficient. Differential insulated neighborhood was computed using two sided Wald tests (DEseq2) and adjusted for multiple comparison using Benjamini-Hochberg correction. c) The insulation scores of strong and week CTCF loops are similar between serum- and 2i-ESCs. p=2.1e-107 (ser) and p =8.5e−113 (2i), Wilcoxon rank-sum test (two-sided). n=11,441 strong-strong and n=5,222 weak-weak boundaries. Box=5-75th percentile; bar=median; whiskers=5-95th percentile. d) The sum of interaction reads of all connected enhancers per promoter is highly similar between serum- and 2i-ESCs with the exception of H3K27me3-associated promoters. Circle size: number of interactions per promoter. n=250 up-regulated and n=250 down-regulated genes, R-values represent Pearson’s correlation coefficient. e) Scatter plots showing the intensities of promoter interactions for the 250 most induced genes in serum- and 2i-ESCs. Even after increasing the CHICAGO-score and read count cutoff, the vast majority of non-H3K27me3-associated genes show similar interactions in serum- and 2i-ESCs. Genes with differential interactions (FDR < 0.1) mainly encompass genes decorated by H3K27me3 (shown in blue). f) Grid analysis of the number of significant interactions found by increasing the interaction read count (x-axis) and CHICAGO significance threshold (y-axis). While the number of interactions rapidly decrease (panel 1, blue), the number of differential interactions remains low, but may increase due to a lower multiple correction penalty (panel 2, green). Increasing the read and score cutoff resulted in a minor increase in the number of differential non-H3K27me3 interactions (panel 3, yellow). Increasing the CHICAGO significance threshold biased towards long-range interactions, dominated by CTCF-CTCF loops (panel 4, red). Increasing the read count cutoff biased the analysis towards short-range interactions, which have a higher proximity ligation background. g) Scatter plots showing the intensities of promoter interactions for the 250 most induced genes in serum- and 2i-ESCs across different biological replicates. Similar interaction intensity in both culture conditions is observed even after increasing the sequencing depth five-fold. R-values represent Pearson’s correlation coefficient. n=250 genes.

Supplementary Fig. 3 Association of chromatin interactions and the occupancy of CTCF, SMC1, YY1 and H3k4me1.

a) Representative examples of chromatin interactions for a gene upregulated in 2i-ESCs (Ablim2) and a gene upregulated in serum-ESCs (c-Myc). Promoter-mediated cHiC interactions are shown in blue. 4C-seq profiles display the running median of 7 DpnII fragments. Dots represent the mean and error bars show the SEM based on n=4 biological replicates in serum-ESCs, 2i-ESCs and EpiLSCs. Fold-changes and significance were determined in bins of 5 kb, 1 Mb up- and downstream of the viewpoint. Differential interactions consistent between replicates are marked with asterisks (p < 0.1, * p < 0.05, ** p < 0.001, *** p < 1e-5; Two sided Wald test (DEseq2)). The viewpoint position is marked with an eye symbol. b) Normalized promoter cHiC interaction intensities of 10,211 genes in serum-ESC and EpiSCs (see text). The majority of the genes (85%) have no differential interactions. The 250 genes most up- or downregulated in the serum-to-2i transition also display few differential interactions between serum-ESCs and EpiSCs. Moreover, most differential interactions are found at promoters associated with H3K27me3 in serum-ESCs. c) Pearson correlation between biological duplicates employed in ChIP-seq analysis. Numbers in the graph indicate Pearson correlation coefficients. n=41630 CTCF_peaks, 35873 SMC1_peaks,19924 YY1_peaks. d) Change in H3K27ac level across interaction-anchors does not correlate with change in cHiC interactions. n=8,016 differential E2 peaks (involving 56,451 significant pairwise interactions), R-value represents Pearson’s correlation coefficient. e) Violin plot showing the occupancy of CTCF, SMC1, YY1 and H3k4me1 across differential enhancers that strongly change H3K27ac between serum- and 2i-ESCs (n=4,303 2i-specific E2 enhancers). Violin=kernel density estimation; center=median; bar=interquartile range; thin line=95% confidence interval. f) Interaction frequency between strongly differential (fold-change >= 3) H3K4me1 sites in 2i- vs serum ESCs (left) and serum-ESCs vs EpiSCs (right). Change in H3K4me1 is significantly associated with change of interaction intensity (One-sided Wilcoxon signed-rank test). Few promoters- and enhancers (n=236 enhancers) show differential H3K4me1 between serum- and 2i-ESCs, and this number is significantly higher (n=399 enhancers) when ESCs are compared with EpiSCs. P-values from left to right: P=0.03; P=0.008; P=7.2e-6; P=2.6e-9, Wilcoxon signed-rank test (one-sided); Box=5-75th percentile; bar=median; whiskers=5-95th percentile. n numbers represent enhancers and are stated in the graph for each group.

Supplementary Fig. 4 Features of serum- and 2i-enhancers.

a) Pearson correlation between biological duplicates employed in ChIP-seq analysis. Numbers in the graph indicate Pearson correlation coefficients. n=32,448 H3K27ac peaks. b) Pearson correlation between biological duplicates employed in ATAC-seq analysis. Numbers in the graph indicate Pearson correlation coefficients. n=54,488 ATAC-seq peaks. c) Scheme of the number of total and differential peaks identified in ChIP-seq or ATAC-seq experiments. d) Scatterplot showing the correlation between ATAC-seq and H3K27ac signal in total ATAC-seq peaks (n=54,488 peaks) identified in serum or 2i-ESCs. ATAC-seq signal was measured under the ATAC-seq peaks and H3K27ac signal was measured in 4 kb windows flanking the ATAC-seq peak-centers. Black, green and blue colors represent E1, E2 and E3 regions as identified in materials and methods and Fig. 2. Each dot represents one ATAC-seq region. e) Number of differential enhancers detected at different time points of the serum-to-2i transition. f) Average profile of different histone marks in different classes of distal regulatory elements. The peak center ± 5 kb is shown.

Supplementary Fig. 5 Motif enrichment in dynamic 2i-enhancers.

a) Promoter connectivity of different enhancer types based on cHiC data. Only dynamic enhancers that differ between the 2i- and serum-ESCs were used. n=1,645 2i-specific E1 enhancers, n=2,941 2i-specific E2 enhancers, n=814 2i-specific E3 enhancers Random baited-DpnII-fragments (n = 5,000) or random genomic regions (n=5,000) were employed as positive and negative controls, respectively. P-values were computed using two-tailed ‘N-1’ Chi-squared test with 95% confidence interval. Asterisks represent P-values < 0.0001. b) RNA-expression of genes connected to different classes of 2i-specific enhancers. cHiC was used to assign enhancer–promoter connections. Box=5-75th percentile; bar=median; whiskers=5-95th percentile. n=1,645 dynamic E1 sites; n=2,941 dynamic E2 sites, n=814 E3 sites. Only dynamic enhancers that differ between the 2i- and serum-ESCs were used. c) Bar plot showing the number of differential distal regulatory elements at different stages of serum-to-2i transition. d) Boxplots represent the temporal change in chromatin accessibility for 2i-specific E1 and E2 regions and during 2i-to-serum transition. Box=5-75th percentile; bar=median; whiskers=5-95th percentile. n=3,815 dynamic E1 sites; n=4,303 dynamic E2 sites. e) Motif enrichment in 2i-specific E1 and E2 distal regulatory elements. Top enriched motifs are shown in the graph. f) RNA expression of different members of the nuclear receptor family at different time points of the serum-to-2i transition. g) Linear regression analysis connects all motifs to change in ATAC-seq signals of all detected peaks. Motifs highly associated to the epigenetic mark (both in regression coefficient and P-value) are annotated in the figure. n= 7,340 differential distal ATAC-seq peaks. P-values derived from t-statistic and adjusted for multiple comparison using Benjamini–Hochberg correction. h) Motif co-occurrence for all discriminative motifs (from regression analysis, motif enrichment and random forest classification) in total distal regulatory elements. Graph clustering: 4 groups (green, yellow, blue and red), edge width relates to significance (padj. ranging for 1e-10 – 1e-300), ball size relates to the number of peaks that are differential between serum and 2i (differential E1 and E2). n=15,960 differential E1/E2 peaks. P-values were derived from a Chi-squared statistic and adjusted for multiple comparisons using Benjamini–Hochberg correction.

Supplementary Fig. 6 TF occupancy at dynamic 2i-enhancers.

a) Motif enrichment in serum-specific and EpiLSC-specific E1 and E2 distal regulatory elements. Top enriched motifs are shown in the graph. b) Correlation between biological replicates employed in ChIP-seq analysis for ESRRB, NR5A2 and KLF4 transcription factors. Two independent clones were used for each transcription factor and ChIP-seq experiments were performed in serum or 2iDay3 ESCs. Numbers represent the Pearson correlation coefficient. n=18,135 NR5A2-GFP peaks, n=3,743 ESRRB-GFP peaks, n=13,921 KLF4-GFP peaks. c) Motif analysis on ChIP-seq binding sites of ESRRB, NR5A2 and KLF4 transcription factors. High enrichments for known motifs and de novo motifs were found in the corresponding binding sites confirming the specificity of ChIP-seq signals. n=18,135 NR5A2-GFP peaks, n=3,743 ESRRB-GFP peaks, n=13,921 KLF4-GFP peaks. P-values were derived from hypergeometric tests (two-sided). d) Co-occurrence of different transcription factor motifs in E1 and E2 enhancers. e) Average profiles of NANOG, OCT4 and SOX2 for different types of 2i-specific enhancers. As a control, 2,000 random ATAC-seq peaks were used. f) Genomic distribution of differential binding sites (defined in Galonska et al.) of NANOG, OCT4 and SOX2 in serum- and 2i-ESCs. The genomic distribution of these transcription factors was linked to E1-E3 distal regulatory elements identified in this study.

Supplementary Fig. 7 ESRRB is linked to enhancer H3K27ac in 2i-ESCs.

a) ESRRB occupancy at ESRRB-binding sites that are activated during EpiSCs to 2i-ESCs reprogramming41 (Adachi et al.). ESRRB binding sites were classified into E2- and non-E2 type enhancers. P-values (Wilcoxon rank sum test; two-sided); left to right: p=0.007, p=0.35, p=2.7e-07, p=9.4e-29, p=8.9e-64. n=1,648 ESRRB binding sites from Adachi et al.41 that overlap E2 enhancers, n=1,648 sampled non-overlapping enhancers. b) Western blot analysis of Tamoxifen-inducible Esrrb-knockout ESCs cultured in serum or 2i-medium. Experiment was repeated twice with similar results. c) Representative microscopy pictures of Tamoxifen-inducible Nr5a2-knock out ES-cells cultured in serum or 72 h 2i-medium and treated with Tamoxifen (4-OH) for Cre-activation and LoxP excision. The experiment was repeated twice with similar results. Scale bar in the image=100 µm. d) Luciferase activity of representative E2 enhancers. Tamoxifen-inducible Esrrb-knockout ESCs were transfected with the reporter constructs and cultured in 2i-medium. After 12 h of transfection, medium was changed to culture-medium supplemented or not with Tamoxifen (4-OHT) and cells were harvested after 48 h and assayed for luciferase activity. Values represent the means of n = 3 biological replicates. Renilla-luciferase was employed for normalization and signal intensities represent Firefly-luciferase / Renilla-luciferase values. e) qRT-PCR analysis of Tamoxifen-inducible Esrrb-knockout ESCs cultured in serum or 2i-medium. Two biological replicates were used per condition. f) Average profile of H3K27ac at 2i-specific E2-enhancers in long term serum- or 2i-cultured ESCs. Cells were treated with Tamoxifen (4-OHT) for 48 h followed by 72 h culture in serum- or 2i-medium (similar to Fig. 8b). Signal in a ±5 kb window flanking the peak center is shown. Reads from two biological samples per condition were pooled. g), Boxplot showing the ChIP-seq levels of H3K27ac and ESRRB-occupancy at E2 mono-allelic enhancers in F1-ESCs. Biallelic enhancers were used for comparison. Reads in input DNA were measured to exclude the bias in allele mapping. n=882 mono-allelic enhancers with higher signal at 129 alleles, n=366 mono-allelic enhancers with higher signal in Castaneus alleles, n=3,527 bi-allelic enhancers. Wilcoxon rank-sum test (two-sided). From left to right: P-value < 2.2e-16, P-value < 2.2e-16, P-value < 2.2e-16, P-value = 5.9e-12. Box=5-75th percentile; bar=median; whiskers=5-95th percentile.

Supplementary Figure 8

Raw western blot image.

Supplementary information

Supplementary Information

Supplementary Figs. 1–8, Supplementary table titles/legends

Reporting Summary

Supplementary Table 1

List of new generated samples and public data used in this study.

Supplementary Table 2

Normalized RNA expression at different time points of serum–2i transitions.

Supplementary Table 3

Total and differential cHiC interactions in serum- and 2i-ESCs and the statistics of 4C-seq libraries.

Supplementary Table 4

Total ATAC-seq and H3K27a sites, dynamic enhancers (E1–E3), CTCF, YY1, SMC1, ESRRB, NR5A2, KLF4 peaks.

Supplementary Table 5

Table summarizing the primers used in this study.

Supplementary Table 6

Statistics source data.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Atlasi, Y., Megchelenbrink, W., Peng, T. et al. Epigenetic modulation of a hardwired 3D chromatin landscape in two naive states of pluripotency. Nat Cell Biol 21, 568–578 (2019). https://doi.org/10.1038/s41556-019-0310-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1038/s41556-019-0310-9

This article is cited by

Search

Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing