Multiscale 3D genome rewiring during PTF1A-mediated somatic cell reprogramming into neural stem cells

Morphological and gene expression changes during PTF1A-driven reprogramming of human somatic cells into neural stem cells

In our previous work, we have demonstrated that by direct somatic cell reprogramming, the non-neural progenitor TF PTF1A alone is able to directly convert both mouse and human fibroblasts, which do not normally express PTF1A, into induced neural stem cells (iNSCs)³³. To further determine whether PTF1A exerts such a strong “master TF” effect on cell transdifferentiation by altering the pattern of 3D genome dynamics and the 3D chromatin architecture during the reprogramming process, we introduced doxycycline (Dox)-inducible PTF1A-GFP or GFP lentiviruses³³ into the HFFs. Subsequently, we carried out the reprogramming procedure and conducted integrated multiomics experiments and analyses including 3D genome, epigenome and transcriptome analyses at different stages of the reprogramming process: HFFs, pre-NSCs (pre-induced NSCs, cells at day 7 after reprogramming initiation), and hiNSCs (cells at day 14 after reprogramming initiation) (Fig. 1a).

a Schematic of the experimental design. Human foreskin fibroblasts (HFFs) were infected with doxycycline-inducible PTF1A-GFP lentiviruses and cultured in the N3 medium for 7 days to generate pre-neural stem cells (pre-NSCs), which were further transdifferentiated into human induced neural stem cells (hiNSCs) by culture for additional 7 days. As a control, HFFs were infected with GFP lentiviruses. Experiments were always performed at these three timepoints (HFF, pre-NSC and hiNSC). The colored circle, hexagon, triangle, square, and rhombus below the cells represent the experiments or analyses of Hi-C, CTCF CUT&Tag, RNA-seq, PTF1A CUT&Tag, and H3K27ac CUT&Tag in corresponding cells. The 3D genome structures display dynamic changes over time, referred to as the 4DN (4D nucleome) pattern. The 4DN pattern and the impact of PTF1A on this pattern during somatic cell reprogramming were obtained by integrating the multiomics experiments and analyses. b Representative images of HFFs infected with either GFP (+GFP) or PTF1A-GFP (+PTF1A) lentiviruses at days 0, 7 (pre-NSC) and 14 (hiNSC) post-infection, showing morphological changes over induction time. HFFs infected with PTF1A-GFP lentiviruses exhibited a rounder and smaller morphology at day 7, which progressed to form spheroid neurosphere-like structures by day 14. In contrast, control HFFs displayed no such changes in morphology over time. c Representative images of pre-NSCs and hiNSCs under bright-field microscopy. d Quantification of PTF1A-induced neurospheres. HFFs (5 × 10⁴) were seeded into each well of 12-well plates, infected with GFP or PTF1A-GFP viruses, and neurosphere-like structures in each well were then counted on days 7 and 14 following virus infection. There were significantly more neurospheres in PTF1A-induced samples on day 14. Data are presented as mean ± SD (n = 4). The asterisk indicates significance in unpaired two-tailed Student’s t-test: *p < 0.01. e hiNSCs (+PTF1A) and corresponding HFFs infected with GFP lentiviruses (+GFP) were immunolabeled for NES (nestin), OLIG2, PAX6, or SOX2, and counterstained with nuclear DAPI. f hiNSCs were immunostained for OCT4 or MAP2 and counterstained with DAPI. g Principal component analysis (PCA) of the HFF, pre-NSC and hiNSC samples based on gene expression determined by RNA-seq. h Expression heatmaps of a set of marker genes for HFFs and NSCs based on RNA-seq experiments of HFFs, pre-NSCs and hiNSCs. i, j NSC-related GO-BP (biological process) and KEGG terms enriched by upregulated DEGs (differentially expressed genes) with log₂FC > 2 and FDR < 0.01 between hiNSCs and HFFs (cutoff of p value = 0.01 and cutoff of q value = 0.05). k Percentage of PTF1A CUT&Tag peaks at promoters of DEGs (purple) and all other genomic regions (gray). Scale bars, 80 μm in (b), (c), (e) and (f).

Full size image

When HFFs were infected with Dox-inducible PTF1A-GFP lentiviruses and cultured in N3 medium supplemented with EGF, bFGF and Dox, they expressed PTF1A protein and underwent overt morphological changes (Figs. 1b, c; S1a, b). By day 7, the cells underwent a reduction in size and acquired a round morphology, characteristic of pre-NSCs. By day 14, these cells aggregated to form spheroid neurospheres, providing clear evidence for successful reprogramming into hiNSCs (Fig. 1b–d). In contrast, control HFFs infected with GFP lentiviruses retained their typical spindle-shaped morphology with clearly oval nuclei, hardly formed any neurospheres, and did not express PTF1A (Figs. 1b, d; S1a, b). Immunofluorescence staining confirmed that the PTF1A-reprogrammed hiNSCs expressed the NSC marker proteins, including NES (nestin), OLIG2, PAX6 and SOX2 (Fig. 1e), but not the pluripotent stem cell marker OCT4 or mature neuron marker MAP2 (Fig. 1f). These results suggest that the reprogrammed cells have successfully acquired the characteristic features of NSCs.

To elucidate the reprogramming route, we validated PTF1A-mediated reprogramming of HFFs at the transcriptome level by RNA-seq using three biological replicates at each time point. High correlation was detected among the replicates of each cell type (Fig. S1c). Principal component analysis (PCA) of the gene expression data among HFFs, pre-NSCs and hiNSCs revealed that hiNSCs were significantly divergent from HFFs and pre-NSCs (Fig. 1g). By identifying and analyzing the differentially expressed genes (DEGs), we observed a progressive transition from fibroblasts toward hiNSCs (Fig. 1h). Genes with high expression levels in hiNSCs were associated with the regulation of neurogenesis and DNA-replication (proliferation), confirming the acquisition of NSC identity. Conversely, the downregulated genes were involved in fibroblast functions, including extracellular matrix/structure organization and cell–substrate adhesion. Notably, genes exhibiting transient activation in reprogramming intermediates (pre-NSCs) were associated with regulation of inflammatory response, wound healing and cell migration (Fig. S1d). This result is in line with previous studies demonstrating that the activation of innate immunity and the MET (mesenchymal-to-epithelial transition) process is crucial for effective reprogramming^35,36.

Expression heatmap analysis revealed a large set of DEGs during the reprogramming process. Notably, NSC marker genes, including PAX6, OLIG2, SOX2, FOXG1, POU3F2 and GBX2, exhibited upregulation in the late stage of reprogramming (pre-NSC → hiNSC), while other NSC markers including NCAM1, HES1, HEY1 and SOX4 were upregulated during the early stage (HFF → pre-NSC)^37,38 (Fig. 1h). Conversely, the expression of HFF marker genes such as TWIST1, LUM, FBLN1, COL8A1 and TNC was gradually downregulated during reprogramming^39,40 (Fig. 1h). Gene ontology (GO) analysis of the upregulated DEGs (hiNSC vs HFF) showed that they were associated with neuron fate commitment and axon development. Additionally, KEGG analysis of these DEGs further revealed that PTF1A played a critical role in committing HFFs toward the neural lineage by activating various signaling pathways including the NOTCH, WNT and TGFβ signaling pathways (Fig. 1i, j). Taken together, these results suggest that HFFs are reprogrammed by PTF1A into hiNSCs with molecular properties similar to those of native NSCs.

Only a small fraction of PTF1A binding sites are associated with local transcriptional regulation

To understand how PTF1A may regulate the expression of the large number of DEGs revealed by RNA-seq, we first performed PTF1A CUT&Tag analysis in pre-NSCs and hiNSCs during the cell reprogramming process (Fig. 1a). By combining the CUT&Tag dataset with the RNA-seq data, we aimed to investigate potential functional links between the PTF1A occupancy and transcriptional activation/repression.

Our analyses revealed pervasive genome-wide binding of PTF1A, with 27,278 and 17,014 binding peaks identified in pre-NSCs and hiNSCs, respectively (Fig. S1e). The majority of the peaks were located in the promoters, introns and intergenic regions (Fig. S1e). Notably, we observed a strong PTF1A binding preference for the prototypical E-box motif SSGGTCACGTGA (p value: 1e − 49 in pre-NSCs and 1e − 70 in hiNSCs), consistent with previous reports^41,42. Intriguingly, we found that only ~5–10% of PTF1A binding sites were located within the promoter regions of DEGs during cell reprogramming (Fig. 1k), indicating that the occupancy of PTF1A on the genome does not directly affect gene transcription in the great majority of cases. Therefore, we speculated that beyond its canonical role as a transcriptional activator or repressor during cell reprogramming, PTF1A may have a broader function, potentially in higher-order genome organization.

PTF1A-bound TAD boundaries are associated with strength-changed TAD reorganization during cell reprogramming

To directly investigate the potential involvement of PTF1A in 3D chromatin and genome organization during cell reprogramming, we conducted in situ DNase Hi-C analysis^43,44 of HFFs, pre-NSCs and hiNSCs (Supplementary Data 1). The resulting Hi-C contact maps exhibited the characteristic checkerboard patterns, indicative of the presence of distinct local genomic domains (Fig. 2a). In addition, they manifested a gradual transition from fibroblasts toward NSCs such that the contact pattern of hiNSCs most closely resembled that of neural progenitor cells (NPCs) derived from iPSCs but highly diverged from that of HFFs (Fig. 2a–f).

a Contact matrices at 100-kb resolution for chr1:0–25 Mb in four cell types: HFFs, pre-NSCs, hiNSCs and neural progenitor cells (NPCs). The Hi-C data used for NPCs are publicly available (GSE158382). Note that the chromatin organization pattern in hiNSCs most closely resembles that in NPCs. b Cluster heatmap of the correlation of contact matrices among the three indicated sample types and NPC at the resolution of 40 kb. c–e Contact matrices at 10-kb resolution, along with PTF1A CUT&Tag tracks over ~2 Mb in four cell types: HFFs, pre-NSCs, hiNSCs and NPCs. f Boxplot showing the stratum correlation coefficients between Hi-C contact matrices of the four cell types at 10-kb resolution.

Full size image

We then conducted an analysis of the TAD and chromatin loop dynamics during the cell reprogramming process⁴⁵, and observed an elevation in the abundance of TAD boundaries during reprogramming (Fig. 3a), with a corresponding decrease in their size (Fig. S2a). Using cLoops⁴⁶, we identified chromatin loops at the three stages of cell reprogramming and found that the loop number and size also varied across these stages (Figs. 3b; S2b), indicating that TADs and loops may undergo reorganization during cell reprogramming.

a The number of TAD boundaries was quantified at each stage of reprogramming using Hi-C data. b Chow–Ruskey diagrams of chromatin loops in HFFs, pre-NSCs and hiNSCs. c Metaplots are used to visualize the enrichment of CTCF CUT&Tag signals at the TAD boundaries in both pre-NSCs and hiNSCs. The results indicate that CTCF CUT&Tag signals are significantly enriched at TAD boundaries at both stages. d Metaplots showing PTF1A CUT&Tag and IgG signal enrichment around PTF1A peak regions in pre-NSCs and hiNSCs. e Venn plots showing overlap between PTF1A and CTCF CUT&Tag peak summits using a 100-bp window in pre-NSCs and hiNSCs, which are presented separately. f CTCF and PTF1A CUT&Tag tracks near HES4 at the indicated cell stages. g hiNSCs were immunostained for PTF1A and CTCF, and counterstained with nuclear DAPI, showing co-localization of the two proteins. Scale bar, 20 μm. h Validation of CTCF interaction with PTF1A. Western blot analysis showed the presence of CTCF protein in the products immunoprecipitated (IP) by an anti-PTF1A antibody, indicating the interaction between CTCF and PTF1A. Lane 1: molecular weight marker (Marker), Lane 2: input, Lane 3: IP with anti-PTF1A, Lane 4: IP with IgG. i A high-resolution average profile of PTF1A and CTCF CUT&Tag signals around CTCF motifs at shared sites in hiNSCs. j Metaplots are used to visualize the enrichment of PTF1A CUT&Tag signals in the CTCF peaks at TAD boundaries in both pre-NSCs and hiNSCs. k A stacked bar plot showing the number of TAD boundaries bound (PTF1A+) or unbound (PTF1A−) by PTF1A, which are either co-bound (CTCF+) or not (CTCF−) by CTCF in pre-NSCs and hiNSCs in a 20-kb window. The results indicate that most PTF1A-bound boundaries are co-occupied by CTCF. l Heatmaps depicting the gradual increase from pre-NSCs to hiNSCs in intensity of PTF1A CUT&Tag signals in the CTCF peaks at TAD boundaries of hiNSCs. m A pie chart showing the proportion of CTCF + PTF1A + TAD boundaries in hiNSCs that existed exclusively in hiNSCs or were already occupied by both CTCF and PTF1A in pre-NSCs, indicating that most CTCF + PTF1A + TAD boundaries in hiNSCs had been bound by CTCF and PTF1A at the pre-NSC stage. n The percentage of dynamic TADs (orange) and unchanged TADs (green) during reprogramming is determined, and the details of dynamic TADs are presented in a pie diagram. The dynamic TADs are classified into three categories based on the stages of reprogramming at which they undergo changes: “Early change only,” refers to TADs that are reorganized in pre-NSCs; “Late change only,” referring to TADs that are reorganized at late stage; and “Continuous change,” refers to TADs that are reorganized in both pre-NSCs and hiNSCs. o The number and percentage of reorganized TADs are assessed at each stage during reprogramming. TAD reorganization results in complex, strength-changed, shifted, merged and split TADs. p Comparison of the proportions of distinct TAD reorganization categories is performed at PTF1A-bound TAD boundaries (PTF1A+) and unbound boundaries (PTF1A−) at early and late stages. q The expression dynamics of DEGs associated with [in close proximity to (<20 kb)] PTF1A-bound and strength-changed TAD boundaries in pre-NSCs. DEGs are classified into upregulated and downregulated ones between hiNSCs and HFFs. The p value is calculated using the Wilcoxon rank-sum test, and significance is denoted as follows: *p < 0.05, **p < 0.01, ***p < 0.001, n.s. no significance. The results suggest that TAD reorganization precedes changes in gene expression.

Full size image

Given that TADs are often demarcated by CTCF, we performed CTCF CUT&Tag analysis for pre-NSCs and hiNSCs (Fig. 1a), and additionally obtained the CTCF CUT&Tag dataset of HFFs from the 4D Nucleome database⁴⁷. These analyses revealed that CTCF binding sites were enriched at TAD boundaries (Fig. 3c), consistent with previous studies^12,48,49. Moreover, our results showed that the proportion of the CTCF-bound TAD boundaries increased during reprogramming (Fig. S2c, d). Interestingly, we also found that a large fraction of PTF1A binding sites (65% in pre-NSCs and 77% in hiNSCs) were colocalized with CTCF binding peaks (Figs. 3d, e; S2e). Integrative analysis revealed a positive correlation between PTF1A binding and CTCF levels (Fig. S2f). For instance, near the HES4 gene, PTF1A binding sites overlap with those of CTCF (Fig. 3f). In addition, through immunofluorescence co-localization and co-immunoprecipitation in hiNSCs, we were able to show a direct interaction between PTF1A and CTCF (Figs. 3g, h; S2g), suggesting that PTF1A may directly recruit CTCF to TAD boundaries. Furthermore, a high-resolution average profile of PTF1A and CTCF CUT&Tag signals around CTCF motifs at shared sites indicates that PTF1A predominantly binds downstream of the CTCF motif (Fig. 3i). We also found an enrichment for PTF1A at the CTCF binding sites at TAD boundaries (Fig. 3j), with essentially all PTF1A-bound TAD boundaries also bound by CTCF (Figs. 3k; S2h). Notably, the majority of the boundaries bound by both PTF1A and CTCF in hiNSCs were already occupied by both proteins at the early reprogramming stage (Figs. 3l, m; S2i). Thus, during somatic cell reprogramming, PTF1A exhibits early binding to TAD boundaries, with increased binding signals at the late stage, highlighting its potential role in shaping TAD organization and chromatin architecture during the reprogramming process.

To determine whether PTF1A is involved in TAD reorganization, we examined the TAD organization patterns during the reprogramming process using TADCompare⁵⁰. 64% of TAD boundaries remained unchanged during reprogramming (Fig. 3n), which is consistent with previous report that 60–70% of TADs exhibit conservation between human ESCs and differentiated cells¹². Notably, we identified five types of TAD reorganization, including complex, strength-changed, shifted, merged and split, based on significant alterations in TAD boundaries or size (Fig. 3o).

Interestingly, compared to boundaries without PTF1A binding, TAD boundaries with PTF1A binding showed a significantly higher occurrence of strength-changed TAD reorganization, representing 39.3% and 29.9% of all TAD reorganizations at the early and late stages, respectively (Fig. 3p). Additionally, strength-changed boundaries exhibited the highest density of PTF1A signal at both the early and late stages of cell reprogramming (Fig. S2j). For instance, near the SOX2, HES4 and GATA2 loci, PTF1A binding can be observed near the TAD boundaries. In HFFs at the original cellular state, these TAD boundaries exhibited relatively weak insulation strength (Fig. S2k–m). The binding of PTF1A to the TAD boundaries was accompanied by an increase in insulation strength and the disappearance of nesting TADs during reprogramming (Fig. S2k–m), suggesting that PTF1A may be involved in the change of TAD boundary insulation.

The boundaries of strength-changed TADs did not correlate with increased or decreased local gene expression (Fig. S2n, o), suggesting that changes in the transcription level are not the primary driving force of TAD reorganization. However, at PTF1A-bound TAD boundaries where transcriptional changes occurred (≥2 fold change between HFFs and hiNSCs), gene expression changes were not significant until after TAD boundaries were strength-changed (Fig. 3q). We further investigated the fold change of gene expression at PTF1A-bound and strength-changed TAD boundaries. Minimal changes were seen in pre-NSCs, while significant differences were observed in hiNSCs (Fig. S2p). These results suggest that PTF1A may influence strength-changed TAD reorganization, which frequently precedes transcriptional changes.

PTF1A-bound boundaries exhibit increased early-stage CTCF binding and enhanced insulation during cell reprogramming

The local chromatin insulation by TAD boundaries can be quantified using the insulation score⁵¹. High correlation of insulation scores was detected between the experimental replicates of each cell type (Fig. S3a). We found a positive correlation between PTF1A CUT&Tag signals and insulation strength, with stronger PTF1A binding linked to lower insulation scores, indicating stronger insulation strength resulting from stronger Ptf1a binding during cell reprogramming (Fig. 4a). Our analysis revealed that during the early stage of reprogramming, the insulation strength of unchanged TAD boundaries was significantly enhanced, followed by a moderate weakening (Fig. S3b). To investigate the impact of PTF1A on the insulation strength of TAD boundaries, we classified all boundaries as either PTF1A-bound (PTF1A+) or PTF1A-unbound (PTF1A−) in both pre-NSCs and hiNSCs (Fig. 4b). Comparative analysis of the insulation scores of the grouped boundaries revealed that PTF1A exhibited a preference for binding to TAD boundaries with insulation strength that was weaker at the original stage (Figs. 4c, d; S3c). Moreover, PTF1A+ boundaries in both pre-NSCs and hiNSCs exhibited a marked enhancement in insulation strength compared to the same position in HFFs (Figs. 4c–f; S3d). Comparison of TAD boundaries with and without PTF1A binding revealed that 73% of PTF1A-bound boundaries had decreased insulation scores, vs 34% without PTF1A binding (Fig. 4g).

a Density plot of PTF1A CUT&Tag signals at TAD boundaries in pre-NSCs with change in insulation score (pre-NSCs vs HFFs). The Pearson correlation coefficient and p value are shown in the plot. b Percentage of PTF1A-bound (PTF1A+) and -unbound (PTF1A−) TAD boundaries in a 20 kb window in pre-NSCs (inner circle) and hiNSCs (outer circle). c, d Metaplots illustrating the dynamics of insulation scores within a 0.5 Mb genomic region around PTF1A+ and PTF1A− boundaries in pre-NSCs (c) and hiNSCs (d). The dashed line denotes the position where the insulation score reaches zero. e, f Boxplots represent the insulation scores within a 20 kb window around PTF1A+ or PTF1A− boundaries in pre-NSCs (e) and hiNSCs (f). P value is calculated using Wilcoxon signed-rank test, **p < 0.01, ****p < 0.0001, n.s. no significance. g Percentage of insulation score decrease (pre-NSCs vs HFFs) at TAD boundaries with or without PTF1A binding. h Boxplots depicting the density of CTCF CUT&Tag signals within a 20 kb window around PTF1A+ or PTF1A− boundaries in pre-NSCs, indicating that PTF1A+ boundaries gain CTCF binding at the early stage of reprogramming. P value is calculated using the Wilcoxon rank-sum test, ****p < 0.0001. i Schematic showing the criteria by which TAD boundaries in HFFs are divided into high- and low-insulation boundaries based on the insulation score. j The bar plot illustrates the percentage of high- and low-insulation boundaries that were bound by PTF1A in pre-NSCs and hiNSCs. These results suggest that PTF1A may prefer binding to low-insulation boundaries in HFFs. k Representative Hi-C maps, insulation scores, and PTF1A CUT&Tag tracks are shown for the indicated cell stages. The arrows indicate the positions of a low-insulation boundary as well as a high-insulation one. l GO term enrichment analysis of upregulated and downregulated DEGs (hiNSC vs HFF) in TADs with PTF1A+ low-insulation boundaries (PTF1A CUT&Tag peaks in hiNSCs) (cutoff of p value = 0.01 and cutoff of q value = 0.05). m Hi-C maps in ~1 Mb region around HOXA1–9 genes. Shown are also insulation score, PTF1A and CTCF CUT&Tag tracks and RNA-seq tracks at the indicated cell stages.

Full size image

Previous studies have demonstrated that reorganized TAD boundaries exhibit significant changes in CTCF enrichment during reprogramming^12,52, and that the number of architectural proteins at TAD boundaries is directly correlated with their strength in Drosophila^53,54. We therefore investigated whether PTF1A binding is correlated with changes in CTCF binding at boundaries. The results showed that at early stage of reprogramming, CTCF density at PTF1A+ boundaries significantly increased, whereas CTCF density at PTF1A− boundaries decreased (Fig. 4h). Therefore, PTF1A binding correlates with increased CTCF binding and enforced insulation at TAD boundaries during reprogramming.

Architectural proteins like CTCF play a critical role in modulating TAD boundaries, enabling or restricting inter-TAD interactions to establish novel gene expression profiles during cell differentiation^55,56. Therefore, we further studied the characteristics of TADs with low-insulation boundaries, primarily subTAD boundaries, that were bound by PTF1A. We defined boundaries with insulation scores in HFFs greater than 0 as low-insulation boundaries and those with scores less than −0.8 as high-insulation boundaries, as illustrated in Fig. 4i. Our investigation revealed that PTF1A bound to approximately twice as many low-insulation boundaries as the high-insulation ones (Fig. 4j, k). We then identified DEGs in those TADs with PTF1A-bound low-insulation boundaries, and found that the upregulated DEGs were associated with the regulation of cell cycle and neuron differentiation while the downregulated DEGs were related to extracellular matrix organization (Fig. 4l). Overall, these results suggest that PTF1A-bound boundaries are associated with increased CTCF binding and enhanced insulation at the low-insulation boundaries of HFFs, which may contribute to the upregulation of NSC-related genes and downregulation of HFF-related genes during reprogramming.

We further examined in more detail several pivotal genes located within TADs with PTF1A-bound low-insulation boundaries, including neurogenesis genes HOXA1–9 and MAPK15 and a key fibroblast gene FBLN5, which are located within TADs with PTF1A-bound low-insulation boundaries (Figs. 4m; S3e, f). The TAD boundaries surrounding these genes exhibited a reinforcement of insulation strength during reprogramming, accompanied by an increase in PTF1A and CTCF binding. Importantly, these changes occurred prior to their transcriptional alterations (Figs. 4m; S3e, f). Thus, PTF1A appears to prefer binding to subTAD boundaries, where CTCF binding is elevated and insulation is enhanced, particularly during the early stage of cell reprogramming.

Reorganization of chromatin loops near low-insulation TAD boundaries during cell reprogramming

Chromatin loops denote robust interactions between two distal genomic regions, which critically contribute to the regulation of gene expression⁴³. Given their importance, we further investigated chromatin loop changes during PTF1A-mediated cell reprogramming, and found that numerous chromatin loops (5794) decreased at the early stage, which we refer to as “early-decreased loops”, while many loops (3384) increased at the late stage, designated as “late-increased loops” (Fig. 5a, b). CTCF plays a key role in facilitating chromatin looping in conjunction with cohesin, possibly via a loop-extrusion mechanism^52,57. Given our finding that PTF1A seems to cooperate with CTCF in modulating the insulation of TADs, we asked whether PTF1A may be involved in the establishment or loss of chromatin loops during cell reprogramming. To explore this, we conducted a motif enrichment analysis focusing on loop anchors and detected significant enrichment for both PTF1A and CTCF binding motifs in pre-NSCs and hiNSCs^42,58 (Fig. 5c). This enrichment was further supported by our observation that PTF1A-binding peaks were enriched at loop anchors (Fig. S4a), and that PTF1A binding signals exhibited a weak positive correlation with the interaction strength of CTCF-CTCF loops (Fig. S4b, c). Notably, an increase in Ptf1a binding is associated with a stronger rise in CTCF-CTCF loop signal (75%) compared to the smaller change observed when Ptf1a binding decreases (29%) (Fig. S4d). Given the enrichment of PTF1A at loop anchors, we investigated whether PTF1A is associated with the reorganization of loops. Our analysis revealed that a substantial portion of early-decreased (78%) and late-increased (56%) chromatin loops were bound by PTF1A at anchors (Fig. 5d), suggesting that PTF1A may play a role in chromatin loop reorganization.

a Number of chromatin loops in HFFs, pre-NSCs and hiNSCs. b Aggregate Hi-C maps around pairs of the decreased chromatin loop regions at the early stage (referred to as early-decreased chromatin loops), as well as around pairs of the increased chromatin loop regions at the late stage (referred to as late-increased chromatin loops). c Enriched CTCF and PTF1A binding motifs at the anchors of early-decreased and late-increased chromatin loops. P values of motif enrichment are determined using cumulative binomial distributions, as applied by Homer. d Percentage of early-decreased chromatin loops, late-increased chromatin loops and static loops that are occupied (PTF1A+) or unoccupied (PTF1A−) by PTF1A at anchors. e Violin plot representing the gene expression dynamics associated with PTF1A-bound, late-increased loops. The p value is calculated using the Wilcoxon rank-sum test, ****p < 0.0001, n.s. no significance. f Violin plot showing the distance between PTF1A-bound, late-increased loops and the closest TAD boundaries. P value is calculated using the Wilcoxon rank-sum test, ****p < 0.0001. g The SOX2 locus as an example of a PTF1A-bound, late-increased loop during cell reprogramming. The gray boxes highlight the two anchors of the PTF1A-bound late, increased loops. The circles indicate the positions of hiNSC-specific contacts. h A schematic indicating that PTF1A preferentially binds at low-insulation boundaries, which are associated with increased CTCF binding and reinforced insulation at the early stage of reprogramming. This may affect gene expression within those TADs. Near low-insulation boundaries, PTF1A binding correlates with new loop formation at the late stage.

Full size image

To gain deeper insights, we focused on analyzing the PTF1A-bound loops and the expression changes of their associated genes (Fig. 5e). In particular, we examined PTF1A-bound, late-increased loops, because we speculated that these loops may be directly regulated by PTF1A binding. Indeed, the genes associated with PTF1A-bound, late-increased chromatin loops were upregulated at the late stage of cell reprogramming, and enriched for GO terms related to neuron fate commitment and differentiation (Figs. 5e; S4f). Notably, genes like SOX2 and POU3F2, known to regulate neurogenesis⁵⁹, were found to be enclosed within PTF1A-bound loops and exhibited upregulation in hiNSCs (Figs. 5g; S4h). In addition, the PTF1A-bound, late-increased chromatin loops are located closer to low-insulation TAD boundaries (mainly subTADs) than to high-insulation ones (Fig. 5f). However, the genes associated with PTF1A-bound, early-decreased loops showed no significant changes during cell reprogramming (Fig. S4e). This lack of expression change might be attributed to the presence of many inhibitory regulatory elements, such as super-silencers⁴⁹, looped to these genes. Nevertheless, the downregulated genes associated with PTF1A-bound, early-decreased loops were enriched for known fibroblast-related functions (Fig. S4g), suggesting that this type of loop is necessary for downregulating gene expression related to fibroblast identity. For instance, TWIST1 and FBLN2 are known to play a key role in fibroblasts^39,40, and their loop interactions were disrupted in pre-NSCs, accompanied by their expression downregulation (Fig. S4i, j). Together, these results suggest that PTF1A may regulate gene expression perhaps by re-organizing the chromatin loops near low-insulation TAD boundaries during cell reprogramming (Fig. 5h).

PTF1A binding is correlated with H3K27ac deposition in TADs with low-insulation boundaries

The maintenance of cell identity has been shown to rely heavily on particular chromatin environments established by lineage TFs to regulate the transcriptional landscape^20,60. To further dissect the molecular mechanism underlying PTF1A-mediated regulation of neurogenic genes within TADs, we examined tracks of active histone modification H3K27ac in the entire genome. A total of 12,111, 51,857 and 12,440 CUT&Tag peaks were identified at the three timepoints of cell reprogramming (Fig. 6a), with most of these peaks distributed in the intron, intergenic and promoter regions, like the PTF1A peaks (Figs. S1e; S5a). Remarkably, we observed an enrichment of H3K27ac signals at the PTF1A binding sites (Fig. 6b). Moreover, the proportion of H3K27ac peaks overlapping with PTF1A peaks increased gradually during cell reprogramming, ultimately reaching 86% in hiNSCs (Figs. 6c; S5b), suggesting a correlation between PTF1A binding and H3K27ac deposition. Furthermore, only about one fourth of the H3K27ac CUT&Tag peaks in pre-NSCs were retained in hiNSCs (Fig. 6a), and PTF1A binding appears to be associated with the retention of these H3K27ac peaks during reprogramming, with transient peaks in pre-NSCs likely contributing to dynamic gene expression changes important for reprogramming (Fig. S6).

a Quantification of the number of H3K27ac CUT&Tag peaks at the three stages of cell reprogramming. b Metaplots are used to visualize the enrichment of H3K27ac CUT&Tag signals around the PTF1A CUT&Tag peak regions in both pre-NSCs and hiNSCs. The results indicate that there is significant enrichment of H3K27ac CUT&Tag signals around PTF1A peak regions at both stages. c Stacked bar plots depicting the fraction of H3K27ac CUT&Tag peaks overlapping (PTF1A+) or not overlapping (PTF1A−) with PTF1A CUT&Tag peaks at the three stages of cell reprogramming. d Heatmaps of the PTF1A and H3K27ac CUT&Tag signals and RNA-seq signals around increased PTF1A CUT&Tag peak regions (upper panel) or decreased PTF1A CUT&Tag peak regions (lower panel) in pre-NSCs and hiNSCs. Each row represents a 1-kb region centered on differentially remodeled PTF1A regions, sorted by PTF1A signal enrichment. e Line plots showing H3K27ac CUT&Tag signals around increased PTF1A (left panel) and decreased PTF1A (right panel) regions. f Pairwise correlations of fold change (FC) between altered PTF1A and H3K27ac signals at late stage of reprogramming. The Pearson correlation coefficient (R) was used for the calculation. g Venn plot showing the overlapped peaks (61%) between increased PTF1A peaks and increased H3K27ac peaks, which are designated as BPHI peaks. h–j Genome browser tracks of the indicated PTF1A and H3K27ac CUT&Tag signals and RNA-seq signals across the NFYA, PAX6 and FOXG1 loci at different stages of reprogramming. k Comparison of the distance between BPHI peak regions and the indicated TAD boundaries. P value is calculated using Wilcoxon rank-sum test, ****p < 0.0001. l Comparison of the density of BPHI peaks in TADs with the indicated boundaries reveals that these peaks are more enriched within TADs that have low-insulation boundaries. P value is calculated using Wilcoxon rank-sum test, ****p < 0.0001. m The genome browser tracks around OLIG2. Shown are Hi-C maps, insulation score, PTF1A and H3K27ac CUT&Tag tracks, and RNA-seq tracks at the indicated cell stages. n A schematic indicating that BPHI peaks near low-insulation boundaries are associated with the upregulation of nearby gene expression.

Full size image

To investigate whether PTF1A shapes a cell stage-specific epigenetic landscape, we compared H3K27ac binding profiles in pre-NSCs and hiNSCs. 38,858 differential H3K27ac peaks were identified at the late stage of reprogramming (Fig. S5c, d), 24% of which showed overlap with the PTF1A differential peaks; and at the early stage, ~30% of increased H3K27ac peaks overlapped with the PTF1A peaks (Fig. S5c). The genomic distribution of H3K27ac and PTF1A differential peaks was highly similar (Fig. S5d, e). We recovered 9705 and 12,981 peak regions with increased and decreased PTF1A binding, respectively (Fig. 6d), and found that increased PTF1A peak regions were accompanied by increased H3K27ac and RNA-seq signals and the opposite was true for the decreased PTF1A peak regions (Fig. 6d, e). The alterations in PTF1A binding and the modifications of H3K27ac show a pronounced correlation (Fig. 6f). Notably, 61% of increased H3K27ac peaks overlapped with increased PTF1A peaks (BPHI peaks) when hiNSCs were compared to pre-NSCs (Fig. 6g), suggesting that increased PTF1A binding correlates with most of the H3K27ac deposition during the late stage of reprogramming. Integration with RNA-seq data further revealed that 74% of promoter-associated, upregulated DEGs with increased H3K27ac peaks also showed increased PTF1A signals nearby (Fig. S5f). This implies that PTF1A binding is correlated with H3K27ac deposition and gene transcription during cell reprogramming. For example, a dramatic increase in H3K27ac levels occurred in the increased PTF1A peak regions across the NFYA, PAX6, FOXG1 and HOX 1-9 loci (Figs. 6h–j; S5g). In agreement, genes whose enhancers were newly bound by PTF1A were upregulated in both pre-NSCs and hiNSCs compared to HFFs, whereas the expression of genes with newly active H3K27ac peaks unbound by PTF1A were unaltered in pre-NSCs and downregulated or upregulated in hiNSCs (Fig. S7).

Given the preferential binding of PTF1A to low-insulation TAD boundaries, we examined the relationship between BPHI peaks and TAD strength, and found that BPHI peaks were located closer to low-insulation boundaries than to high-insulation ones, and that their distribution abundance in TADs with low-insulation boundaries was greater than that in TADs with high-insulation ones (Figs. 6k, l; S5h). For example, BPHI peaks across the OLIG2 locus were near a low-insulation TAD boundary (Fig. 6m). These results suggest that PTF1A binding is correlated with H3K27ac deposition in TADs with low-insulation boundaries, thereby forming BPHI peaks associated with upregulation of nearby gene expression (Fig. 6n).

PTF1A preferentially activates super-enhancers in TADs with low-insulation boundaries to promote neural stem cell fate

Given the importance of super-enhancers (SEs) in the specification of cell identity¹⁷, we aimed to identify SEs at each stage of cell reprogramming. To achieve this, we utilized the ROSE algorithm¹⁹ to rank the clustered H3K27ac peaks and defined H3K27ac-rich regions as SEs (Fig. S8a). The results showed that the number of SEs varied during reprogramming (Fig. 7a). Additionally, SEs were found to span large chromatin domains, with a median size ranging from 12 to 47 kb (Fig. S8b). Of particular interest, we found a significant enrichment of PTF1A signals within SEs, particularly in hiNSCs (Fig. 7b). Furthermore, we noted that PTF1A had already occupied binding sites within hiNSC-specific SEs at the pre-NSC stage (Fig. 7c), indicating its early involvement in SE regulation during reprogramming. Consistent with this, there were 22% of PTF1A CUT&Tag peaks overlapping with SEs at the pre-NSC stage (Fig. S8c). Moreover, PTF1A displayed stronger binding strength at SE-overlapping CTCF-CTCF loop anchors compared to other anchors in both pre-NSCs and hiNSCs (Fig. S8d).

a Chow–Ruskey diagrams of SE (super-enhancer) regions in HFFs (blue border), pre-NSCs (orange border) and hiNSCs (green border). The color of the border surrounding each intersection corresponds to the cell type whose SE regions overlap. The central circle represents the SE regions that overlap in all three cell types. The inner areas with two- or three-color borders represent the overlap of fewer cell types. The peripheral areas with a unique color border represent the cell type-specific SE regions. The area of each intersection is proportional to the number of SE regions within the intersection. b Metaplots depicting the enrichment of PTF1A CUT&Tag signals around pre-NSC-specific SEs and hiNSC-specific SEs. c Percentage of hiNSC-specific SEs bound by PTF1A in pre-NSCs or hiNSCs. d The distance between SEs and the closest TAD boundaries suggests that SEs are more likely to be located closer to low-insulation boundaries than to high-insulation ones. P value is calculated using the Wilcoxon rank-sum test, ****p < 0.0001. e The distribution percentage of TADs with SEs indicates that SEs are more likely to be located in TADs with low-insulation boundaries (pink) than in TADs with high-insulation ones (red). f, g Comparison of the enrichment of SEs in hiNSCs within TADs with high- or low-insulation boundaries reveals that SEs in hiNSCs are more enriched within TADs with low-insulation boundaries than those with high-insulation ones. P value is calculated using the Wilcoxon rank-sum test, ****p < 0.0001. h Gene expression dynamics of SE-associated genes within TADs with low-insulation boundaries in hiNSCs, at the three stages of cell reprogramming. P value is calculated using the Wilcoxon rank-sum test, ****p < 0.0001, n.s. no significance. i Volcano plot of the SE-associated genes within TADs with low-insulation boundaries. They display significant expression change between hiNSCs and HFFs and most are upregulated. j GO-biological process terms enriched by SE-associated genes within TADs with low-insulation boundaries (cutoff of p value = 0.01 and cutoff of q value = 0.05). k Genome browser tracks are displayed to visualize the PTF1A and H3K27ac CUT&Tag signals, identified SEs, and RNA-seq data across the DLX1, DLX2, HOXD11, HOXD13 and MYBL2 loci at the indicated stages of reprogramming.

Full size image

Notably, similar to the binding location preference of PTF1A, SEs in hiNSCs tended to be located near the low-insulation boundaries instead of the high-insulation ones (Fig. 7d). Upon closer examination, we found that TADs with low-insulation boundaries contained a higher number of SEs and had a greater enrichment for SEs compared to TADs with high-insulation boundaries (Fig. 7e–g). This observation aligns with previous findings⁶¹ and further supports the idea that PTF1A may play a role in activating SEs within subTADs during cell reprogramming. Next, we investigated the SEs located in TADs with low-insulation boundaries and their associated genes. A substantial proportion of these genes, including several critical neural TF genes, were upregulated at the late stage of reprogramming (Fig. 7h, i). These upregulated genes likely played a vital role in cell fate transition since they were enriched for GO terms related to neuron fate commitment, regulation of neurogenesis, and positive regulation of cell cycle (Fig. 7j). Specifically, SEs proximal to genes involved in neuron fate commitment, such as DLX1, DLX2, HOXD11, HOXD13, DLL1, POU2F1, GATA4, HES4 and WNT3A loci, as well as the MYBL2 locus, which is involved in cell cycle progression, were all activated by PTF1A binding in hiNSCs (Figs. 7k; S8e–i). These results suggest that the activation of SEs in TADs with low-insulation boundaries is critical for cell fate transition and that PTF1A plays a crucial role in this process.

We subsequently investigated how the chromatin loops endow the regulatory function of SEs. We identified 1311, 858 and 433 SE-anchored loops at the HFF, pre-NSC and hiNSC stages of cell reprogramming, respectively (Figs. 8a; S9a). More than half of the SE-anchored loops involved SE-typical enhancer (TE) and SE–SE interactions in pre-NSCs and hiNSCs (Figs. 8a; S9a). Pre-NSCs and hiNSCs were enriched for SEs associated with shorter chromatin loops compared to HFFs (Fig. S9b), and the SEs in pre-NSCs and hiNSCs displayed lower interaction density than those in HFFs (Fig. S9c). Most SE-associated chromatin loops exhibited binding of PTF1A at one or both anchors (Figs. 8b; S9d), indicating potential roles of PTF1A in SE-anchored loops. Notably, SE-looped genes displayed cell-stage specificity (Fig. 8c). For instance, the promoters of ID3 and POU2F1 were directly looped to SEs bound by PTF1A in hiNSCs (Fig. 8d, e). These results suggest that PTF1A may be involved in SE clustering to regulate the expression of key reprogramming genes (Fig. 8f).

a Pie chart showing the fractions of the indicated SE (super-enhancer)-anchored loop types in hiNSCs. TE typical enhancer. b Fractions of the indicated loop types that are occupied (PTFA+) or unoccupied (PTFA−) by PTF1A at anchors in hiNSCs. c Chow–Ruskey diagram of SE-looped genes during reprogramming. d, e The ID3 and POU2F1 loci as examples of PTF1A-bound loops that interact with one or more SEs in hiNSCs. The gray boxes highlight the regions proximal to the anchors of the PTF1A-bound, SE-anchored loops, which encompass or are near either the promoter or SE regions. Hi-C heatmaps are generated using pyGenomeTracks (top). The circles indicate the positions of hiNSC-specific contacts. The blue boxes represent loop anchors. f Working model depicting the relationship between PTF1A binding and reorganization of 3D genome architecture during somatic cell reprogramming. There are changes at multiple hierarchical organization levels of the 3D genome during reprogramming. PTF1A shows a preference for binding at low-insulation boundaries (mainly boundaries of subTADs) in HFFs, which correlates with increased CTCF binding at these sites. During reprogramming, subTAD demarcation becomes more pronounced, while the nesting TADs gradually disappear. These changes are associated with the downregulation of HFF (human foreskin fibroblast)-related genes and the upregulation of NSC (neural stem cell)-related genes. Within subTADs in HFFs, PTF1A also binds to and activates enhancers to regulate the expression of genes within the same TAD, while directly and/or indirectly facilitating the deposition of H3K27ac. Furthermore, PTF1A is involved in the dynamic reorganization of chromatin loops, such as the formation of SE (super-enhancer)-promoter and TE (enhancer)-promoter loops, which are coupled with the upregulation of NSC transcription factor (TF) genes, including SOX2, POU2F1, ID3 and POU3F2.

Full size image

• SEO Powered Content & PR Distribution. Get Amplified Today.
• PlatoData.Network Vertical Generative Ai. Empower Yourself. Access Here.
• PlatoAiStream. Web3 Intelligence. Knowledge Amplified. Access Here.
• PlatoESG. Carbon, CleanTech, Energy, Environment, Solar, Waste Management. Access Here.
• PlatoHealth. Biotech and Clinical Trials Intelligence. Access Here.
• Source: https://www.nature.com/articles/s42003-024-07230-1