Transcriptional network dynamics during the progression of pluripotency revealed by integrative statistical learning

Abstract

The developmental potential of cells, termed pluripotency, is highly dynamic and progresses through a continuum of naive, formative and primed states. Pluripotency progression of mouse embryonic stem cells (ESCs) from naive to formative and primed state is governed by transcription factors (TFs) and their target genes. Genomic techniques have uncovered a multitude of TF binding sites in ESCs, yet a major challenge lies in identifying target genes from functional binding sites and reconstructing dynamic transcriptional networks underlying pluripotency progression. Here, we integrated time-resolved 'trans-omic' datasets together with TF binding profiles and chromatin conformation data to identify target genes of a panel of TFs. Our analyses revealed that naive TF target genes are more likely to be TFs themselves than those of formative TFs, suggesting denser hierarchies among naive TFs. We also discovered that formative TF target genes are marked by permissive epigenomic signatures in the naive state, indicating that they are poised for expression prior to the initiation of pluripotency transition to the formative state. Finally, our reconstructed transcriptional networks pinpointed the precise timing from naive to formative pluripotency progression and enabled the spatiotemporal mapping of differentiating ESCs to their in vivo counterparts in developing embryos.

Figure 1.

Identification of dynamically regulated TFs in pluripotency progression. (A) Schematics of TF expression and their binding at promoter-proximal (P1, P2 and P3) and distal (D1, D2 and D3) sites during pluripotency progression from naive to formative states. Dash lines (red) represent interactions between distal TF binding sites and their target genes. (B) Schematic summary of the time-course trans-omic dataset utilised in this study for reconstructing and characterising the transcriptional networks in pluripotency progression from naive ESCs to EpiLCs that represent formative state. (C) Volcano plot of genes profiled on both transcriptome and proteome levels in the trans-omic dataset (27). TFs that are dynamically regulated (DR) during the ESC to EpiLC transition are highlighted in blue and within these DR TFs, those that have been profiled previously using ChIP-seq in ESCs are highlighted in red. (D) Pie charts showing the distribution of promoter-proximal target genes (TSS ± 1 kb) and putative distal target genes (TSS > 1 kb) with chromatin loops (Pol2-ChIA-PET) for each TF according to its ChIP-seq profile in ESCs.

Figure 2.

Prediction and validation of TF target genes in transition from naive to formative pluripotency using AdaEnsemble and trans-omic data. (A) Time-course showing the log₂ fold change (compared to time = 0) in expression for mRNA (green) and protein (blue) for Sox2 and c-Myc. Bars represent standard deviation among biological replicates (n = 2 for mRNA and n = 4 for protein). Pearson's correlation coefficient for concordance between protein and mRNA are shown. (B) Time-course expression profiles of putative target genes (i.e. putative candidates) supported by chromatin loops from Pol2-ChIA-PET (top) and AdaEnsemble predicted target genes (bottom) for Sox2 and c-Myc, respectively. Profiles are divided into those from promoter-proximal and distal target genes, and then further divided into those from mRNA and protein levels. (C, D) Cumulative distribution showing degree of correlation (Pearson's) between time-course expression profiles of each TF with its putative target genes (green and black) and those predicted by AdaEnsemble (red and yellow), on mRNA (C, red) and protein (D, yellow) levels. P-values were computed using Wilcoxon Mann-Whitney U test (two-sided). (E) Log₂ fold change in mRNA after Sox2, Nr5a2, Nanog, Klf4, Esrrb knockdown, or Otx2 knockout in ESCs compared to WT ESCs. Wilcoxon Mann–Whitney U test (one-sided) are performed on AdaEnsemble-identified target genes (iv) versus all quantified genes (i), proximal and distal target genes by nearest TSS assignment of all TF binding sites identified in ChIP-seq data (ii), and putative target genes by chromatin loop support (Pol2-ChIA-PET) (iii), respectively, and the largest p-value is displayed as an upper bound for all pairwise comparisons for each TF. (F) Similar to (E) but in c-Myc overexpressed ESCs compared to WT ESCs.

Figure 3.

Evaluation of predictability of histone marks and Pol2 ChIP-seq on TF target genes. (A) Overall predictive power (in terms of area under the ROC curve [AUC]) of histone marks and Pol2 ChIP-seq data, averaged across all seven TFs, when assessed against the intermediate predictions from the initial training phase of AdaEnsemble. (B) Bi-clustered heatmap showing relationship among histone modifications and Pol2 ChIP-seq in terms of predictiveness (scaled AUC) in TF target genes. The dashed box highlights H3K4me1 and H3K4me3 marks for their contrast between naive pluripotency TFs (Sox2, Nanog, Esrrb, Nr5a2, Klf4) and formative pluripotency TFs (Otx2 and c-Myc). (C) Levels of H3K4me1 (6 h) and H3K4me3 (48 h) signal at the promoters of the target genes of naive and primed pluripotency TFs.

Figure 4.

Characterisation of TF target genes and transcriptional networks in naive and formative pluripotency. (A) Heatmap showing proportion (as quantified by Jaccard index) of genes regulated by pairs of TFs. (B) Venn diagram showing a three-way overlap of target genes of naive TFs (union), Otx2 and c-Myc. (C) Over-representation of gene ontology (GO) of target genes unique to naive TFs (naive-specific), Otx2 (Otx2-specific) and c-Myc (c-Myc-specific). (D) Contribution of target genes that are themselves TFs toward the total transcriptional regulation. (E) Correlation between (D) and the percentage of target genes that are themselves TFs (y-axis) for each of the seven TFs.

Figure 5.

Formative target genes in naive pluripotency. (A) Boxplot for comparing H3K27me3 and H3K27ac signal at the promoter of AdaEnsemble-identified TF target genes of naive and formative TFs during pluripotency progression. (B) Time-course of average log₂ fold change (relative to 0 min) in mRNA expression for the target genes of naive and formative TFs. Shaded areas represent standard deviation among target genes. (C) Chromatin accessibility (ATAC-seq) of c-Myc, Otx2, and naive TF binding sites in ESCs and EpiLCs. (D) Schematic illustration of the poised chromatin associated with formative target genes in naive pluripotency.

Figure 6.

Dynamic rewiring of transcriptional networks from naive to formative pluripotency. (A) TF networks prior to differentiation (naive pluripotent state). AdaEnsemble-identified TF target genes (prediction probability > 0.95) that are themselves TFs are included in the network reconstruction. Edges measure mean expression of each TF–gene pair. (B) Dynamic change of TF networks during the transition from naive to formative pluripotent states. Edge colour across times reflects the change of expression of each TF–gene pair in the transcriptional networks. (C) Relative change in all pairwise expressions from naive to formative states (see ‘Materials and Methods’ section for details). (D) Mapping the differentiating ESCs to the equivalent epiblast cell populations by the activity of reconstructed transcriptional networks at different time points. Colour of the corn plots denotes enrichment from low (green) to high (red). (E) Correlation between the expression of transcription networks from single cells of the E5.5 epiblast and the ESC to EpiLC transcriptional networks at different time points.