RA is sufficient to drive SE commitment

An SE differentiation system was developed using a pluripotency growth medium (mTeSR1) supplemented with RA. After 7 days of RA treatment, hESCs lost their embryonic stem cell morphology and were induced into a nearly homogeneous epithelial-like morphology, with robust expression of the SE markers KRT7 and KRT18 and the epidermal master regulator TP63, while the embryonic stem cell pluripotency marker NANOG was undetectable (Fig. 1A, B). To further confirm the SE differentiation phenotype, we performed RNA sequencing (RNA-seq) comparing RA-driven differentiated cells with SE generated according to previously published protocols [16]. Principal components analysis (PCA) revealed that the transcriptome landscape of D7-differentiated cells driven by RA closely resembled that of SE (Fig. 1C). Moreover, gene set enrichment analysis (GSEA) showed that the transcriptome of D7 SE cells in our system was similar to that of mouse SE at embryonic day 9 [23] (Supplementary Fig. 1A). These results suggest that RA is capable of driving hESCs to acquire an SE phenotype.

Fig. 1

RA induces SE commitment. A Phase contrast images of the differentiating hESCs during 7 days of RA induction. Scale bar, 100 μm. B Immunostaining for NANOG, KRT7, KRT18, and TP63 in the differentiated cells on D0 and D7. Scale bar, 100 μm. C PCA of RNA-seq data of cells at each time point during RA-driven differentiation and SE cells. D Heatmap of expression changes of the genes that were differentially expressed between at least two time points during 7 days of RA induction. These differential genes were clustered into four groups based on different variation tendency over time. The color bar shows the relative expression value (Z score of TPM [transcripts per kilobase of exon model per million mapped reads]) from the RNA-seq data. E Representative GO terms (biological process) identified from the cluster II, III, and IV genes. F Schematic diagram of RA-induced SE differentiation

To determine whether RA-driven SE possesses further differentiation potential, we sequentially cultured the cells in keratinocyte maturation medium and terminal differentiation medium. We observed that the RA-driven SE could successfully differentiated into keratinocytes (hereafter “iKCs”) expressing KRT5, KRT14, TP63, ITGB4, LAMB3, and COL7A1, with decreased KRT8 and KRT18 expression (Supplementary Fig. 1B-D). After calcium-induced terminal differentiation, the iKCs highly expressed the terminal differentiation genes KRT1, KRT10, FLG, LOR, and TGM1, indicating their ability to undergo terminal differentiation (Supplementary Fig. 1B-D). Furthermore, gene ontology (GO) enrichment analysis showed that the genes preferentially expressed in iKCs (compared with RA-driven SE) were significantly associated with epidermis and skin development, whereas those enriched in terminally differentiated iKCs (compared with iKCs) were preferentially linked to keratinocyte differentiation, cornification, and establishment of the skin barrier (Supplementary Fig. 1E). Thus, these results validate the functional capacity of the RA-driven SE.

To further understand the process of SE fate determination induced by RA, we analyzed transcriptome dynamics during differentiation. 6,167 genes that were differentially expressed (DEGs) across at least two time points were identified and clustered into 4 groups based on the main temporal expression patterns observed during differentiation (Fig. 1D). Specifically, Cluster I DEGs, typified by pluripotency markers NANOG and POU5F1, were highly expressed in hESCs and sharply downregulated as SE differentiation proceeded (Fig. 1D and Supplementary Fig. 1F). Cluster II DEGs, typified by ribosome assembly genes, UTP14A and RRP36, were progressively upregulated during the first 2 days of differentiation, with enrichment in the ribosome biogenesis, RNA transport, and protein folding, and then gradually declined thereafter (Fig. 1D, E and Supplementary Fig. 1F ). Ribosome biogenesis is known to play key roles in the transition from self-renewal to differentiation [24, 25], implying that during the first 2 days of differentiation, cells prepared for SE initiation by increasing protein synthesis. Cluster III DEGs, typified by SE markers KRT18 and KRT7, were slightly upregulated during the first 2 days of differentiation, abundantly expressed on D3, and peaked on D4 (Fig. 1D and Supplementary Fig. 1F, G ). These DEGs were preferentially linked to epithelium development, BMP signaling pathway, and actin filament organization (Fig. 1E). Finally, cluster IV DEGs, typified by SE maturation genes TP63 and COL3A1, were gradually upregulated as SE differentiation progressed and were significantly associated with epithelium development, response to BMP, and extracellular matrix organization (Fig. 1D, E and Supplementary Fig. 1F).

Taken together, these results confirm that RA is sufficient to drive hESCs to adopt an SE phenotype and that the RA-driven differentiation process can be categorized into four distinct stages: D0 (pluripotency phase), D1-D2 (SE prime phase), D3-D4 (SE initiation phase), and D5-D7 (SE maturation phase) (Fig. 1F).

BMP4 signaling serves as the RA downstream to mediate SE commitment

Given that previous studies have demonstrated the essential role of BMP signaling in SE fate determination and that BMP4 has been used to induce SE fate in vitro [15, 16, 26], we hypothesized that RA may induced SE differentiation by activating endogenous BMP4 signaling. By profiling dynamic changes in BMP family member expression during RA-driven SE differentiation, we discovered BMP4 as the predominantly expressed member during differentiation (Fig. 2A). BMP4 expression is activated on the second day of differentiation, peaks on the third day, and then slightly downregulates while maintaining its expression. (Fig. 2B). Additionally, we observed that RA induced SMAD1/5/8 phosphorylation (pSMAD1/5/8), the downstream nuclear effector of BMP signaling, indicating activation of BMP signaling (Fig. 2C). Through functional validation, we found that BMP4-knockdown cells maintained hESC-like morphology and failed to acquire an epithelial-like morphology, with undetectable expression of pSMAD1/5/8, KRT18 and TP63 (Fig. 2C). Meanwhile, BMP4-knockdown cells showed a significant decrease in the expression of SE genes (KRT8, KRT18, KRT7, KRT19, TP63, TFAP2C, and GRHL3) and an upregulation of pluripotency genes (NANOG, POU5F1, SOX2, and DNMT3B) (Fig. 2D, E). Moreover, GO analysis revealed that genes downregulated upon BMP4 depletion were significantly associated with SE development biological processes, such as response to BMP stimulus, embryonic epithelium morphogenesis, and extracellular matrix organization (Fig. 2F). These results demonstrated that BMP4 signaling as downstream of RA is required to mediate SE fate determination.

Fig. 2

RA drives SE fate determination by activating BMP4 signaling. A Heatmap of the expression changes of BMP family members during SE differentiation. The color bar shows the expression value (TPM) from the RNA-seq. B Immunostaining of BMP4 during 7 days of differentiation. Scale bar, 100 μm. C Phase contrast images and immunostaining for pSMAD1/5/8, KRT18, and TP63 in scrambled shRNA- and shBMP4-treated hESCs after seven days of differentiation. Scale bar, 100 μm. D qRT-PCR analysis of representative genes in scrambled shRNA- and shBMP4-treated hESCs after seven days of differentiation. qRT-PCR values were normalized to the values in scrambled shRNA group. Values are presented as means ± SD (n = 3 biological replicates; **P < 0.01; ***P < 0.001 t test). E Heatmap showing expression levels of representative genes in shBMP4- and scrambled shRNA-treated hESCs after seven days of differentiation. F GO (biological process) analysis of the downregulated genes upon BMP4 knockdown. G Scatterplot of differential accessibility in shBMP4- versus scrambled shRNA-treated hESCs after seven days of differentiation. Sites identified with significant differential accessibility are highlighted in color (red, peaks increased; green, peaks decreased). H Transcription factor motif enrichment in the regions with increased (left) or decreased (right) accessibility in shBMP4- versus scrambled shRNA-treated hESCs after seven days of differentiation. I Genome browser tracks comparing ATAC-seq signal across NANOG, POU5F1, SOX2, KRT8/18, KRT7, and TP63 loci of scrambled shRNA- or shBMP4-treated hESCs after seven days of differentiation

Although the role of BMP4 in SE development is well-established, the molecular regulatory mechanisms through which BMP4 controls this process remain unclear. To elucidate the mechanism underlying BMP4-mediated SE commitment, we performed an assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq). Upon BMP4 depletion, we observed substantial changes in global chromatin accessibility and identified 52,644 differentially accessible regions, in which 64.6% and 35.4% of the loci display a marked decrease and increase in chromatin accessibility, respectively (Fig. 2G). Moreover, TF motif enrichment analysis of these loci revealed that the motifs for SE regulators (TFAP2A, TFAP2C, GATA3, and GRHL2) were specifically enriched in loci with decreased accessibility, whereas the motifs for core pluripotent factors (OCT4, SOX2, and NANOG) were specifically enriched in loci with increased accessibility (Fig. 2H). This suggests a significant repression in the binding of SE regulators and an activation in the binding of hESC critical regulators after BMP4 depletion. Specifically, the loci of SE genes (KRT7, KRT8, KRT18, and TP63) favored a closed chromatin state, and the loci of pluripotency genes (NANOG, POU5F1, and SOX2) exhibited a gain in chromatin accessibility in BMP4-knockdown differentiated cell (Fig. 2I). Collectively, these results demonstrate that BMP4 functions downstream of RA and plays a decisive role in SE development by mediating the establishment of the SE chromatin landscape.

RXRA/B prime lineage-specific epigenetic landscape and regulates SE gene expression

Given that RA played a vital role in SE commitment, we aimed to clarify how RA impels SE determination. RA typically exerts its effects through its nuclear receptors that act as RA-modulated TFs to regulate gene expression [11]. In this study, we focused on elucidating the role of RXRs in RA-driven SE fate determination. We first examined the expression of the three RXR family members (RXRA, RXRB, and RXRG) and observed that RXRA and RXRB were ubiquitously expressed during RA-driven differentiation, whereas RXRG expression was not detectable (Fig. 3A). Using CRISPR/Cas9 technology, we generated RXRA-knockout (KO) and RXRB-KO hESCs and found that either RXRA- or RXRB-KO prevented RA-induced SE differentiation of hESCs and inhibited the expression of SE marker genes in the SE initiation stage (Supplementary Fig. 2A-C). However, after 7 days of RA induction, RXRA- and RXRB-KO-differentiated cells displayed a transition to an epithelial-like morphology, with no difference in SE marker gene expression compared with that in the wild-type (WT) group (Supplementary Fig. 2A-C). Based on these results, we speculated that RXRA and RXRB likely function with partial redundancy.

Fig. 3

RXRA/B primes SE epigenetic and transcriptomic landscape. A RNA-seq analysis of RXRA, RXRB, and RXRG expression during seven days of RA induction. B Left two panels, phase contrast images of the differentiating hESCs on D2 and D7. Scale bar, 100 μm. Right two panels, immunofluorescence staining of NANOG, KRT18, and TP63 in WT and RXRA/B DKO cells after seven days of differentiation. Scale bar, 100 μm. C Representative GO terms (biological processes) identified from the genes highly expressed in WT or RXRA/B DKO cells after seven days of differentiation. D Heatmap of expression of pluripotency, RA and SE-related genes in WT and RXRA/B DKO cells after two days or seven days of differentiation. E Histogram showing the distribution patterns of RXRA and RXRB peaks in the differentiated cells on D2. F Heatmaps of binding signals of H3K4me3, H3K4me1, H3K27ac, and ATAC-seq at the center of RXRA or RXRB peaks. G Pie chart showing the percentage of RXRA or RXRB peaks located in active enhancers or active promoters. H Scatterplot of differential H3K27ac, H3K4me1, and H3K4me3 peaks in WT versus RXRA/B DKO cells after two days of differentiation. Sites identified as differentially bound with significance (FDR < 0.05) are shown in red. I Metaplots of average H3K27ac, H3K4me1, and H3K4me3 density around RXRA or RXRB peaks in WT and RXRA/B DKO-differentiated cells on D2. J Venn diagram showing the overlap among downregulated genes and the target genes of decreased H3K27ac, H3K4me1, and H3K4me3 signals in differentiated WT versus RXRA/B DKO cells on D2. K Pie chart showing percentage of the overlap genes obtained from (J) with or without RXRA and RXRB binding. L Representative GO terms (biological process) identified from the genes bound by both RXRA and RXRB obtained from (K)

Subsequently, we generated RXRA/B double-KO (DKO) hESCs (Fig. 3B and Supplementary Fig. 2D). We found that RXRA/B DKO-differentiated cells on D7 maintained hESC clonal morphology with robust expression of NANOG, and undetectable expression of KRT18 and TP63, validating the critical role of RXRA/B in RA-driven SE commitment (Fig. 3B). PCA demonstrated that RXRA/B DKO-differentiated cells exhibited transcriptome landscapes similar to that of hESCs (Supplementary Fig. 2E). Specifically, RXRA/B DKO-differentiated cells exhibited a remarkable reduction in the expression of SE genes (KRT8, KRT18, KRT7, KRT19, TP63, TFAP2C, and GRHL3) and maintained expression of the pluripotency-related genes (NANOG, POU5F1, and SOX2) (Supplementary Fig. 2F-G). GO analysis showed that the genes downregulated in RXRA/B DKO-differentiated cells were enriched for epithelium development, response to RA, and extracellular matrix organization, whereas genes preferentially expressed in RXRA/B DKO-differentiated cells were associated with stem cell population maintenance and developmental growth regulation (Fig. 3C). These results highlight that RA-driven SE differentiation is governed by RXRA/B.

Given that significant morphological alteration and increased expression of RA and SE -related genes were observed at the SE prime stage (Fig. 3B, D), we focused on investigating the regulatory mechanisms underlying the function of RXRA/B during early RA differentiation. To characterize the global genomic occupancy of RXRA/B, we performed chromatin immunoprecipitation sequencing (ChIP-seq). Intriguingly, the distribution analysis of RXRA and RXRB binding maps demonstrated distinct DNA-binding modes: RXRA binding sites were mainly located in intergenic and intronic regions away from the transcription start site, as well as promoter regions, whereas RXRB binding sites were predominantly in promoter regions (Fig. 3E and Supplementary Fig. 2H). Moreover, epigenetic modifications in RXRA and RXRB binding sites were determined and annotated based on ChromHMM segmentation annotation (Supplementary Fig. 2I). RXRA binding sites were mainly located in active enhancers (47%) and promoters (33%) rich in H3K4me3, H3K27ac, and H3K4me1 histone modifications, whereas RXRB binding sites were mostly located in active promoters (64.4%) rich in H3K4me3 and H3K27ac histone modifications (Fig. 3F, G).

Lineage commitment is highly dependent on the specific epigenetic environment established by the core lineage-determining TFs that dominate transcriptional decisions [27, 28]. To investigate the roles of RXRA and RXRB in rewiring the epigenetic landscape during SE differentiation, we compared the binding maps of H3K27ac, H3K4me1, and H3K4me3 between WT and RXRA/B DKO-differentiated cells at the SE prime stage. We identified 23,313 differential H3K27ac, 37,325 differential H3K4me1, and 11,501 differential H3K4me3 peaks, most of which displayed diminished deposition of H3K27ac (88.7%), H3K4me1 (91.7%), and H3K4me3 (92.6%) (Fig. 3H). Notably, the RXRA and RXRB binding sites showed decreased H3K27ac, H3K4me1, and H4K4me3 signals in RXRA/B DKO cells (compared with that in the WT group), suggesting their essential role in epigenetic remodeling underlying RA-driven SE differentiation (Fig. 3I). By combining the epigenome and transcriptome data, we found that 43.4% (359/827) of the genes that were significantly downregulated in RXRA/B DKO cells were accompanied by simultaneous decreases in H3K27ac, H3K4me1, and H3K4me3 signals near their loci (Fig. 3J). Additionally, most of these genes (65.5%) contained both RXRA and RXRB binding sites near their loci and were significantly enriched in biological processes of epithelium development, RA metabolism, histone modification, and BMP signaling (Fig. 3K, L). Specifically, RXRA and RXRB bound to loci near RA target genes (CRABP2, CYP26A1, and HOXA1) and SE regulators (TFAP2C, GRHL3, TFAP2A, and GATA3), regulating their H3K27ac, H3K4me1, and H3K4me3 deposition and gene expression, whereas in the deficiency of RXRA/B, these processes were remarkably disturbed (Fig. 3K). Taken together, these data demonstrate that RXRA/B modulate SE lineage-specific epigenetic and transcriptomic remodeling.

TCF7 is required for RXRA/B to mediate exit from pluripotency and initiates SE fate determination

To further uncover the RXRA/B regulatory process, we screened the top 10 differentially expressed downstream TFs of RXRA/B (Fig. 4A). Notably, several of these TFs, including TFAP2C and GATA3, are well-known SE development regulators [16], indicating their potential roles in SE fate determination. We then analyzed the trend of expression changes of the top 10 candidate TFs during the early stages of SE differentiation. We focused on TCF7, which has not previously been reported in SE development and exhibited a substantial increase on the first day of differentiation, to further investigate its involvement in SE commitment (Fig. 4B).

Fig. 4

TCF7 is required to mediate the transition from pluripotency to SE initiation. A Heatmap of the expression of the top 10 candidate transcription factors during the first three days of differentiation. B Immunostaining of TCF7 during the first three days of differentiation. Scale bar, 100 μm. C Phase contrast images and immunostaining for KRT18 and TP63 in scrambled shRNA- and shTCF7-treated hESCs after three or seven days of differentiation. Scale bar, 100 μm. D Volcano map of differentially expressed genes of shTCF7- versus scrambled shRNA-treated hESCs after seven days of differentiation. The significant DEGs (fold change ≥ 2 and q-value < 0.05) are highlighted in red (upregulated genes) or green (downregulated genes). E GSEA for Cluster I genes in the gene expression matrix of shTCF7- versus scrambled shRNA-treated hESCs after seven days of differentiation. F GSEA for Cluster III and IV genes in the gene expression matrix of shTCF7- versus scrambled shRNA-treated hESCs after seven days of differentiation. G Pie chart showing the distribution patterns of TCF7 peaks in the differentiated cells on D3. H Heatmaps of binding signals of H3K4me1 and H3K4me3 at the center of TCF7 peaks in the differentiated cells on D3. I Upset plot depicting the overlap between genes bound by TCF7 and the genes that were differentially expressed upon TCF7 depletion. J GO (biological process) analysis of the differentially expressed target genes of TCF7. K Genome browser view of TCF7 peaks and RNA-seq signal across representative loci

We knocked down TCF7 to explore its function during RA-driven SE differentiation. Remarkably, the absence of TCF7 repressed RA-induced SE initiation and prevented cells from acquiring epithelial-like morphology, as indicated by undetectable expression of KRT18 at the SE initiation stage (Fig. 4C). After 7 days of induction, TCF7-knockdown differentiated cells retained the clonal morphology and failed to differentiate into the SE, as evidenced by undetectable KRT8 and TP63 expression, a significant decrease in the expression of SE genes (KRT8, KRT18, KRT7, KRT19, TP63, TFAP2C, and GRHL3), and an upregulation of pluripotency-related genes (POU5F1, SOX2, and NANOG) (Fig. 4C, D). GSEA revealed that Cluster I genes, which are specifically expressed in hESC (shown in Fig. 1D), were upregulated in TCF7-depleted cells, whereas Cluster III and IV genes, which are highly expressed during the SE initiation and maturation stages (shown in Fig. 1D), respectively, were significantly suppressed (Fig. 4E, F).

To further elucidate the mechanism underlying TCF7-mediated SE commitment, we generated the binding map of TCF7 by performing ChIP-seq. The distribution analysis demonstrated that TCF7 binding sites were located in intergenic, intronic, and promoter regions (Fig. 4G). Combined with H3K4me1 and H3K4me3 epigenetic maps, we found that TCF7 binding sites were clustered into two groups: cluster 1 represented promoters with H3K4me3 positivity; cluster 2 were enhancers defined by H3K4me1 enrichment. This indicated that TCF7 regulates downstream genes through promoters and enhancers (Fig. 4H). By integrating TCF7 ChIP-Seq and RNA-seq analyses, we observed that TCF7 targets were enriched among downregulated (854/1682, 50.7%) and upregulated (390/811, 48%) genes, indicating that TCF7 regulates both gene activation and silencing (Fig. 4I). GO analysis showed that the upregulated target genes linked to stem cell population maintenance, response to fibroblast growth factor, and cell growth (Fig. 4J). The downregulated target genes were associated with epithelium development, regulation of the BMP and RA receptor signaling pathways, and extracellular matrix organization (Fig. 4J). Specifically, TCF7 bound to loci near pluripotency genes (NANOG, SOX2, and DPPA4), which were silenced during SE differentiation, and to loci close to SE genes (TFAP2C, GRHL3, and GATA3), which were activated as differentiation progressed (Fig. 4K). Taken together, these results indicate that TCF7, acting as a downstream effector of RA, promotes hESC to exit pluripotency and drives SE differentiation.

MSX2 regulates SE commitment via opening chromatin and mediating SE gene expression

To further decipher the downstream mechanisms by which TCF7 mediates SE commitment, we identified its potential downstream mediators. We screened the top 10 TFs that were bound by TCF7 and significantly affected by its depletion (Fig. 5A). Among these, MSX2 attracted our attention and was chosen for further investigation to explore its role in SE initiation (Fig. 5B). We first examined the expression of MSX2 during SE initiation and observed that MSX2 was absent in hESCs and on the first day of differentiation (Fig. 5C). Its expression remained relatively low on day 2 and underwent a substantial increase during the SE initiation stage (Fig. 5C). Focusing on understanding the role of MSX2 in SE commitment, we knocked down MSX2 and found that MSX2-knockdown differentiated cells exhibited a phenotype resembling that of TCF7-knockdown cells (Fig. 5D). These cells maintained hESC clonal morphology, failed to acquire an epithelial-like morphology, and had undetectable SE marker KRT18 and SE maturation gene TP63 (Fig. 5D, E). Additionally, the loss of MSX2 remarkably suppressed the activation of SE-specific keratins (KRT7, KRT8, KRT18, and KRT19) and SE regulators (GRHL2, GRHL3, TFAP2A, TFAP2C, and GATA3) (Fig. 5E and Supplementary Fig. 3A, B). GSEA showed that loss of MSX2 decreased the expression of Cluster III and IV genes enriched in SE (Fig. 5F and Supplementary Fig. 3C). In contrast, Cluster I genes, which are highly expressed in hESCs, were upregulated when MSX2 was knocked down (Supplementary Fig. 3D). These results indicate that MSX2 is required to mediate SE initiation. Interestingly, despite MSX2 being a downstream target of BMP4 signaling, perturbing MSX2 expression inhibited the expression of BMP4 and its downstream target genes (GRHL3, GATA3, ID2, MSX1, DLX5, and DLX6) (Fig. 5E and Supplementary Fig. 3A-B). Furthermore, transcriptome changes occurring upon loss of BMP4 and MSX2 were highly similar, demonstrating a positive feedback regulatory loop between MSX2 and BMP4 (Fig. 5G and Supplementary Fig. 3E).

Fig. 5

MSX2 opens chromatin and activates SE gene expression. A Volcano map of differentially expressed genes of shTCF7- versus scrambled shRNA-treated hESCs after seven days of differentiation. The candidate downstream activators of TCF7 are highlighted in dark blue text. B Genome browser view of indicated ChIP-seq and RNA-seq signal across MSX2 loci. C Immunostaining of MSX2 during the first three days of differentiation. Scale bar, 100 μm. D Phase contrast images and immunostaining for KRT18 and TP63 in scrambled shRNA- and shMSX2-treated hESCs after three or seven days of differentiation. Scale bar, 100 μm. E Heatmap of expression of represented genes in scrambled shRNA- and shMSX2-treated hESCs after three or seven days of differentiation. F GSEA for Cluster III and IV genes in the gene expression matrix of shMSX2- versus scrambled shRNA-treated hESCs after seven days of differentiation. G GSEA for downregulated genes upon BMP4 knockdown in the gene expression matrix of shMSX2- versus scrambled shRNA-treated hESCs after three or seven days of differentiation. H Pie chart showing the distribution patterns of MSX2 peaks in the differentiated cells on D3. I Heatmaps of binding signals of H3K4me1, H3K27ac, H3K4me3 and ATAC-seq at the center of MSX2 peaks in hESCs and the D3-differentiated cells. J Enrichment of TF motifs identified by HOMER at MSX2 peaks. K Heatmaps of binding signals of GRHL3, GRHL2, TFAP2C, and GATA3 at the center of MSX2 peaks. L Pie chart of the percentages of downregulated genes induced by MSX2 knockdown with or without MSX2 binding in D3-differentiated cells. M GO (biological process) analysis of the target genes of MSX2. N Left panel: Scatterplot of differential accessibility in shMSX2- versus scrambled shRNA-treated hESCs after three days of differentiation. Sites identified as significant differential accessibility are shown in color (red, peaks increased; blue, peaks decreased). Right panel: TF motif enrichment in the regions of differential accessibility in shMSX2- versus scrambled shRNA-treated hESCs after three days of differentiation. O Metaplots of average ATAC-seq density around the MSX2 target genes in shMSX2- and scrambled shRNA-treated hESCs after three days of differentiation. P Genome browser view of MSX2, ATAC-seq and RNA-seq signal across MSX2, TFAP2C, GRHL3, and KRT8/18 loci

Next, we sought to dissect the contribution of MSX2 in SE initiation and performed ChIP-seq to identify genome-wide binding sites of MSX2. Distribution analysis demonstrated that MSX2 binding sites were predominantly located in intergenic and intronic regions distant from the transcription start site of the target genes (Fig. 5H and Supplementary Fig. 3F). Combining this with epigenetic maps, we showed that MSX2 binding sites in SE-initiating cells had increased chromatin accessibility and higher enrichment of active histone marks (H3K27ac and H3K4me1) compared with those in hESCs, implying that MSX2 predominantly occupied enhancers to activate gene transcription during SE initiation (Fig. 5I). This was validated by gene expression analysis, which exhibited significantly higher expression levels of MSX2 target genes in SE-initiating cells compared with those in hESCs (Supplementary Fig. 3G). Furthermore, TF motif enrichment analysis revealed significant co-enrichment of motifs for SE regulators (GATA3, TFAP2A, TFAP2C, and GRHL2) at MSX2 binding sites (Fig. 5J). Notably, integrating the ChIP-seq data of MSX2 with previously published SE regulator (GATA3, TFAP2A, TFAP2C, GRHL2, and GRHL3) ChIP-seq data [26], we found frequent co-binding of MSX2 and SE regulators (Fig. 5K). These findings suggest that MSX2 may cooperate with these factors to regulate SE lineage gene expression and SE commitment. Moreover, by integrating MSX2 ChIP-seq and transcriptome data, we observed that most genes downregulated upon MSX2 knockdown (61%) exhibited MSX2 binding, and these MSX2 target genes were enriched in GO terms associated with pan-epithelial development processes, such as morphogenesis, proliferation, migration, cell junction assembly, and response to BMP, indicating that MSX2 mediated SE differentiation mainly through direct binding (Fig. 5L, M).

To further investigate the mechanisms mediating the effects of MSX2 in SE commitment, we performed genome-wide chromatin accessibility profiling by ATAC-seq. Notably, substantial changes in chromatin accessibility was observed upon MSX2 knockdown (Fig. 5N). Among differentially accessible regions, 93.3% of the loci showed a significant decrease in chromatin accessibility in MSX2-knockdown cells (compared with those in the control cells) and were enriched for motifs of well-known SE regulators, including TFAP2C, TFAP2A, GRHL2, and GATA3. However, only 6.7% of the loci showed a significant increase in chromatin accessibility and were enriched for motifs of core pluripotent factors, including NANOG, OCT4, and SOX2. This suggests a critical role of MSX2 in shaping the SE lineage-specific chromatin accessibility landscape during SE commitment (Fig. 5N). Furthermore, MSX2 knockdown substantially inhibited the opening chromatin accessibility around its target gene loci (Fig. 5O). Specifically, MSX2 bound to loci near itself, SE regulators (TFAP2C and GRHL3), and SE markers (KRT8 and KRT18), opening their chromatin accessibility and activating gene transcription, whereas upon knockdown of MSX2, their chromatin accessibilities were notably interfered, and their gene expression was remarkedly suppressed (Fig. 5P). Collectively, these results demonstrate that MSX2 regulates SE initiation by remodeling nuclear architecture to activate SE lineage-specific genes.

The regulatory hierarchy of RXRA/B, TCF7, and MSX2 in SE lineage regulatory networks

To accurately clarify the regulatory hierarchy among RXRA/B, TCF7, and MSX2 during SE fate determination, we silenced each gene individually and assessed its effect on the expression of the other two genes and BMP4. We found that RXRA/B DKO significantly suppressed the activation of TCF7, MSX2, and BMP4; loss of TCF7 significantly inhibited the expression of MSX2 and BMP4 but did not affect RXRA/B expression; and the absence of MSX2 suppressed BMP4 expression, but had no effect on RXRA/B and TCF7 expression (Fig. 6A). These results establish a regulatory hierarchy governing SE commitment, positioning RXRA/B at the top, followed sequentially by TCF7, MSX2, and BMP4. We next conducted co-immunoprecipitation to further define the roles of these regulators in the RXRA/B-TCF7-MSX2 regulatory network and observed no physical interaction between RXRA/B, TCF7, and MSX2 (Fig. 6B). Taken together, these results demonstrated that RXRA/B, TCF7, and MSX2 control SE commitment by forming the RXRA/B-TCF7-MSX2 regulatory axis (Fig. 6C).

Fig. 6

The hierarchical regulation of RXRA/B, TCF7, and MSX2 during SE commitment. A qRT-PCR analysis of representative genes in the differentiated cells on D3. Top panel, WT and DKO hESCs after three days of differentiation; Middle panel, shTCF7- versus scrambled shRNA-treated hESCs after three days of differentiation; Bottom panel, shMSX2- versus scrambled shRNA-treated hESCs after three days of differentiation. qRT-PCR values were normalized to the values in control group. Values are presented as means ± SD (n = 3 biological replicates; **P < 0.01; ***P < 0.001 t test). B Co-immunoprecipitation analysis of the interactions among endogenous RXRA, RXRB, TCF7, and MSX2 in the differentiated cells on D3. C The graphical summary of the function and mechanism of the RXRA/B-TCF7-MSX2 regulatory axis in RA-driven SE commitment