The regulatory architecture of the primed pluripotent cell state

Abstract

Despite extensive research, the gene regulatory architecture governing mammalian cell states remains poorly understood. Here we present an integrative systems biology approach to elucidate the network architecture of primed state pluripotency. Using an unbiased methodology, we identified and experimentally confirmed 132 transcription factors as master regulators (MRs) of mouse epiblast stem cell (EpiSC) pluripotency, many of which were further validated by CRISPR-mediated functional assays. To assemble a comprehensive regulatory network, we silenced each of the 132 MRs to assess their effects on the other MRs and their transcriptional targets, yielding a network of 1273 MR → MR interactions. Network architecture analyses revealed four functionally distinct MR modules (communities), and identified key Speaker and Mediator MRs based on their hierarchical rank and centrality. Our findings elucidate the de-centralized logic of a "communal interaction" model in which the balanced activities of four MR communities maintain primed state pluripotency.

Fig. 1. Flowchart of experimental design.

a Interactome construction. Two independent EpiSC lines were treated with 5 differentiation conditions and 33 small molecule perturbations, generating 276 expression profiles. These profiles served as input for the ARACNe algorithm to create the “first-generation” EpiSC interactome. b Inference of candidate MRs. 144 expression profiles from 5 distinct lineage-specific differentiation time courses were used to generate signatures. The VIPER algorithm then analyzed these signatures to identify candidate MRs. c Validation of candidate MRs. Two independent screens investigated whether silencing of candidate MRs could modulate the pluripotent state, resulting in 132 validated MRs. The first screen examined whether candidate MR silencing altered endogenous Oct4 protein levels, while the second screen used PLATE-seq to determine if candidate MR silencing recapitulated the EpiSC differentiation signature. d Functional validation of MRs. 70 selected MRs underwent CRISPR-mediated knockout to investigate pluripotency defects through cell culture assays and teratoma formation. e Network modularity analysis. PLATE-seq expression profiles from MR silencing were used to assemble a comprehensive, experimentally-derived regulatory network. Four distinct MR communities were identified within this network. f Regulatory architecture inference. Network hierarchy and centrality were analyzed, followed by further examination of community structure.

Fig. 2. Identification and validation of candidate master regulators of primed state pluripotency.

a Examples of transcriptional regulatory proteins being inactivated, unaffected, or activated during the differentiation treatments. Each heatmap cell represents the VIPER-inferred differential protein activity (DA) expressed as normalized Enrichment Score (NES). Enrichment plots show the distribution of activated (orange) and repressed (purple) transcriptional targets on the 72 h gene expression signatures. The right-most column represents the integrated differential activity of each protein across all time points and treatments (integrated DA). RA: Retinoic acid; WAFB: Wnt3a, Activin A, FGF2, and BMP4; WA: Wnt3a and Activin A. b Oct4 staining of EpiSCs after shRNA-mediated silencing of selected candidate MRs and vehicle or mock-treated controls. Each MR was targeted by two individual hairpins in four independent replicates. Scale bars: 100 microns. c Waterfall plot showing differential Oct4 protein level  (z-score) after silencing of each candidate MR; dashed lines indicate FDR < 0.05 (2-tailed). Raw data are provided in Supplementary Data 5.2. d Scatterplot of the effect of silencing each MR on Oct4 protein level and the differentiation score of EpiSCs, which is determined by the recapitulation of the protein activity signature of EpiSC differentiation by the silencing experiment (see Methods). We observed a significant correlation between these scores (Spearman’s Rho 0.32, p < 10⁻⁴); dotted lines indicate FDR = 0.05 (2-tailed). Candidate MRs showing significant validation (differentiation) score are highlighted in green. Selected MRs labeled in black indicate MRs functionally validated in CRISPR/Cas9 assays, and in blue indicate known key pluripotency regulators. Raw data are provided in Supplementary Data 6.2.

Fig. 3. Functional validation of two master regulators.

a Morphology of KO or control EpiSC lines; representative images for independent lines are shown. b Fold change in cell number after five days of colony formation. c Number of alkaline phosphatase positive colonies. d Time line of neuroectoderm differentiation. e Immunofluorescence staining of biological replicates (clonal lines) for lineage-specific markers. Each MR was characterized by two clonal KO cell lines in at least three independent replicates. f Time line of mesendoderm differentiation. g Staining for lineage-specific markers; number of replicates as in (e). h Western blots for Cbl and Zc3h13. Source data are provided as a Source Data file. i Schematic of the teratoma assay. j Representative images of mutant and control teratomas dissected from contralateral legs of the same animal. k Teratoma weight. l Hematoxylin and eosin (H&E) and immunohistochemical staining of serial sections. m, n Fraction of area showing ectoderm (neural rosette) (m) and endoderm (ciliated epithelium) features (n). Line and whiskers show the mean ± s.d.; statistical significance is shown as FDR, 2-tailed U-test. Scale bars: 100 microns.

Fig. 4. Regulatory network model for primed state pluripotency.

a Regulatory network of primed state pluripotency with 120 nodes and 1273 unidirectional edges. Only significant effects (FDR < 10⁻¹⁰, 2-tailed aREA test) are included. The size of each node is proportional to its Regularized Out-Degree (ROD) score; colors indicate community assignment (strong, FDR < 0.01; weak, 0.01 < FDR < 0.1; unassigned FDR > 0.1). Six outlier MRs showing the highest ROD scores (Speakers) are highlighted. b Heatmap showing the effect of silencing each MR on the activity of all other MRs. The silenced MRs are represented by the columns, and the rows report activity changes for the indicated proteins. The 77 MRs (76 strong community members and Utf1) are grouped by their community membership. The color scale indicates statistical significance (-log₁₀(FDR), 2-tailed aREA test) for each MR protein differential activity; only significant changes in protein activity (FDR < 10⁻¹⁰) are shown.

Fig. 5. Community structure of the primed state pluripotency network.

a Schematic representation of parameters for defining node hierarchy based on their regularized out-degree (ROD) and betweenness score. b Distribution density for ROD across all MRs. The vertical dotted lines indicate the point of equal probability for a mixture of Gaussian models fitted to the data, defining the boundaries for Listeners, Communicators and Speakers. c Heatmap showing network parameters for the 77 MRs, showing scaled values for Degree, ROD and betweenness scores. d–f Violin plots showing distribution of regularized out-degree (d), betweenness (e), and degree (f) (FDR < 0.05, 2-tailed U-test). g Heatmap of the integrated meta-protein activity change for each community across the five differentiation treatments. Each cell represents the weighted average of the relative change in activity vs. untreated EpiSC controls across all MRs of each community, with the weight calculated based on the MR community score (Supplementary Fig. 8). h Heatmap showing GOBP pathways preferentially enriched (-log₁₀(FDR), 1-tailed aREA test) on each of the four communities.