Epigenetic modulation of a hardwired 3D chromatin landscape in two naive states of pluripotency (original) (raw)

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. 1–7 and Supplementary Figs. 1–7 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

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

Author notes

  1. Ehsan Habibi
    Present address: Klarman Cell Observatory, Broad Institute of MIT and Harvard, Cambridge, MA, USA
  2. Onkar Joshi
    Present address: Friedrich Miescher Institute for Biomedical Research, Basel, Switzerland
  3. Shuang-Yin Wang
    Present address: Department of Immunology, Weizmann Institute of Science, Rehovot, Israel
  4. These authors contributed equally: Yaser Atlasi, Wout Megchelenbrink.

Authors and Affiliations

  1. Department of Molecular Biology, Faculty of Science, Radboud University, Radboud Institute for Molecular Life Sciences, Nijmegen, the Netherlands
    Yaser Atlasi, Wout Megchelenbrink, Tianran Peng, Ehsan Habibi, Onkar Joshi, Shuang-Yin Wang, Cheng Wang, Colin Logie, Hendrik Marks & Hendrik G. Stunnenberg
  2. Max Planck Institute of Molecular Cell Biology and Genetics, Dresden, Germany
    Ina Poser

Authors

  1. Yaser Atlasi
  2. Wout Megchelenbrink
  3. Tianran Peng
  4. Ehsan Habibi
  5. Onkar Joshi
  6. Shuang-Yin Wang
  7. Cheng Wang
  8. Colin Logie
  9. Ina Poser
  10. Hendrik Marks
  11. Hendrik G. Stunnenberg

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 toHendrik 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

Rights and permissions

About this article

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