In vivo temporal resolution of acute promyelocytic leukemia progression reveals a role of Klf4 in suppressing early leukemic transformation.

Genome organization plays a pivotal role in transcription, but how transcription factors (TFs) rewire the structure of the genome to initiate and maintain the programs that lead to oncogenic transformation remains poorly understood. Acute promyelocytic leukemia (APL) is a fatal subtype of leukemia driven by a chromosomal translocation between the promyelocytic leukemia (PML) and retinoic acid receptor α (RARα) genes. We used primary hematopoietic stem and progenitor cells (HSPCs) and leukemic blasts that express the fusion protein PML-RARα as a paradigm to temporally dissect the dynamic changes in the epigenome, transcriptome, and genome architecture induced during oncogenic transformation. We found that PML-RARα initiates a continuum of topologic alterations, including switches from A to B compartments, transcriptional repression, loss of active histone marks, and gain of repressive histone marks. Our multiomics-integrated analysis identifies Klf4 as an early down-regulated gene in PML-RARα-driven leukemogenesis. Furthermore, we characterized the dynamic alterations in the Klf4 cis-regulatory network during APL progression and demonstrated that ectopic Klf4 overexpression can suppress self-renewal and reverse the differentiation block induced by PML-RARα. Our study provides a comprehensive in vivo temporal dissection of the epigenomic and topological reprogramming induced by an oncogenic TF and illustrates how topological architecture can be used to identify new drivers of malignant transformation.

The 3D organization of the genome, ranging from nucleosomes to heterochromatin/euchromatin compartments and chromosome territories, provides a fundamental mechanism for genome regulation (Schoenfelder and Fraser 2019; Zheng and Xie 2019). Transcriptional regulatory elements, including enhancers and promoters, are in physical contact to fine-tune the timing and magnitude of gene expression, and perturbation of this contact can profoundly affect cell identity, differentiation, and tumorigenesis (Gröschel et al. 2014; Northcott et al. 2014; Lupiáñez et al. 2015; Flavahan et al. 2016; Hnisz et al. 2016; Akdemir et al. 2020). Epigenetic changes often drive the initiation, maintenance, and progression of cancer, and their reversibility and plasticity make them attractive targets in the clinical field (Dawson 2017). However, little is known about how the genome structure is rewired during the acquisition of oncogenic features, or whether structural changes are functionally linked to epigenome and transcriptome alterations during oncogenic transformation.

Acute promyelocytic leukemias (APLs) represent 10%–15% of acute myeloid leukemias (AMLs) and are characterized by the presence of the t(15;17) chromosomal translocation between PML and RARα (de Thé et al. 1990; Goddard et al. 1991). Expression of the oncofusion protein PML-RARα in hematopoietic stem/progenitor cells (HSPCs) results in a differentiation block at the promyelocytic stage and in malignant transformation, as these cells are able to recapitulate most clinical and morphological features of human APL in transplantation mouse models (Brown et al. 1997; Grisolano et al. 1997; He et al. 1997; Grignani et al. 2000; Westervelt et al. 2003; Guibal et al. 2009; Wojiski et al. 2009). Mechanistic studies using APL cell lines and primary blasts have shown that the PML-RARα oncofusion protein competes with normal RARα functions and recruits histone deacetylases (HDACs), the NuRD chromatin remodeling complex, and Polycomb-repressive complexes (PRCs) to constitutively repress target genes of the TFs RARα and PU.1 (Pandolfi 2001; Villa et al. 2007; Morey et al. 2008; Martens et al. 2010; Wang et al. 2010; Mas and Di Croce 2016). These epigenetic complexes mediate long-range interactions to instruct gene expression programs during development and in tumor cells (Denholtz et al. 2013; Schoenfelder et al. 2015; Mas et al. 2018; Oksuz et al. 2018; Basu et al. 2020). In addition to its repressive functions, PML-RARα binds superenhancer regions to directly transactivate genes that encode key myeloid-determining TFs or enzymes, including GFI1, MPO, WT1, and MYC (Tan et al. 2021). PML-RARα also disrupts PML nuclear bodies, which are structures involved in the control of cell cycle, apoptosis, senescence, DNA damage, and antiviral immunity (Bernardi and Pandolfi 2007; Scherer and Stamminger 2016; Chang et al. 2018). Induction of DNA damage by PML-RARα results in increased mutability, favoring the occurrence of cooperating secondary mutations and development of full-blown leukemia (di Masi et al. 2016; Voisset et al. 2018). Recent reports using APL cell lines and primary blasts showed that PML-RARα mediates the formation of long-range interactions to repress the expression of genes controlling myeloid differentiation and maturation (Li et al. 2018; Wang et al. 2020) and activate the GFI1 superenhancer (Tan et al. 2021). However, these studies did not provide a dynamic perspective of how PML-RARα remodels the genome to impair the function and differentiation of normal primary HPSCs to generate fully transformed leukemic blasts. Here, we used the PML-RARα model system as a paradigm to temporally dissect the dynamics of epigenomic and transcriptomic reprogramming occurring at the onset, during progression, and in full-blown APL leukemias in animal models that faithfully recapitulate human APL clinical features. Our global profiling identified the Klf4 locus as one of the most extensively reorganized genes during PML-RARα-driven APL progression. Klf4 encodes a TF with important roles in myeloid differentiation (Feinberg et al. 2007; Park et al. 2016, 2019a). Klf4 expression has been shown to be lower in samples from AML patients than in those from healthy individuals (Faber et al. 2013b; Morris et al. 2016). However, the function of Klf4 in APL has remained controversial, with a few studies reporting that Klf4 overexpression induces differentiation using the APL cell line HL-60 (Feinberg et al. 2007; Alder et al. 2008; Morris et al. 2016), and others showing that Klf4 expression supports cell growth and survival of the APL cell line NB4 (Lewis et al. 2021). Using our integrative multiomics analysis, we temporally resolved the genomic alterations induced by PML-RARα and showed that the Klf4 locus undergoes extensive reprogramming of enhancer–promoter interactions, transcriptional down-regulation, and gain of repressive histone modifications. We further showed that ectopic overexpression of Klf4 partially restored the phenotypic defects induced by the expression of PML-RARα. This work provides a dynamic model of the genomic reprogramming triggered by an oncogenic TF in vivo and highlights the use of topological information to identify new drivers of malignant transformation.

Results

PML-RARα induces a progressive reorganization of genome architecture

To dissect the dynamic changes induced by PML-RARα in genome architecture and transcription during leukemia progression, we infected primary mouse bone marrow hematopoietic stem/progenitor cells (lineage negative [Lin−]) with lentiviruses carrying an empty vector control or a Flag-tagged human PML-RARα and harvested cells at different stages of APL transformation (Fig. 1A). Stage 0 and stage I corresponded to sorted GFP+ cells transformed with empty vector or PML-RARα-3xFlag vector, respectively (Supplemental Fig. S1A,B). We followed the progressive transformation of cells carrying PML-RARα by culturing them in semisolid media and harvesting after 2 or 4 wk of serial replating (equivalent to stage II or III, respectively). The final stage of APL transformation (stage IV) corresponded to blasts isolated from mice that were transplanted with cells carrying PML-RARα; mice developed leukemia after ∼6 mo (Fig. 1A). We verified that cells expressing PML-RARα-3xFlag showed increased serial replating capacity, impaired differentiation, and promyelocytic morphology, as compared with cells expressing empty vector control (Supplemental Fig. S1C,D). We then used multiple biological replicates of cells from stages 0–IV to generate in situ Hi-C libraries (Rao et al. 2014), RNA-seq libraries, and ChIP-seq libraries in order to comprehensively characterize the genome architecture, the transcriptome, and the epigenome, respectively, during the process of leukemic transformation.

Figure 1.

PML-RARα induces a continuum of A-to-B compartment switches. (A) Experimental strategy and definition of stages. Lineage-negative cells (Lin−) were isolated from bone marrow of adult mice (6 to 8 wk old) and infected with lentiviruses carrying an empty vector (Empty) or a vector containing a Flag-tagged PML-RAR < fusion gene; PML-RARα-3xFlag). Successfully infected Lin− cells (GFP+) were sorted and corresponded to stage 0 (carrying empty vector) or stage I (carrying PML-RARα-3xFlag vector). Stage I cells were then cultured in methylcellulose media and harvested at the second replating (for stage II) or fourth replating (for stage III). Stage I cells were also transplanted into lethally irradiated recipient mice. Blast cells were harvested from bone marrow of mice developing leukemias, at ∼6 mo after transplant (stage IV). (B) Principal component analysis (PCA) based on Hi-C eigenvectors for all autosomes. (C) Scatter plots of first eigenvectors along the time course for chromosome 7. The Y-axis represents the first eigenvector associated with replicate 1 of stage 0, and the X-axis represents the first eigenvector of the different stages (replicate 2 for stage 0, and replicate 1 for the remaining stages). Pearson correlations are highlighted in red. (D) Proportion of the genome that changed compartment during the time course. We assumed that a region was A (or B) if all replicates at the same stage were flagged as A (or B). Regions not consistent between replicates were considered ambiguous and represented 2% of the genome. About 5% of the genome was excluded due to low mapability. (E) Alluvial plot showing the dynamic A-to-B compartment switching of bins during the time course. Stages are represented along the X-axis, and the genomic size is represented on the Y-axis as well as by the width of the ribbons. Bins that did not switch compartments or that were flagged as “ambiguous” at any point were excluded. (F) Stacked bar plots of the number of genes in bins that were (1) in compartment A at stage 0 and switched to compartment B at another stage (left), or (2) in compartment B at stage 0 and switched to compartment A at another stage (right). Only genes that switched compartments from one stage to the next and were stably maintained in the new compartment were considered (i.e., genes in bins that switched compartments more than once during the time course were excluded). (G) KEGG analysis and WikiPathways analysis of 724 genes that switched from A to B compartments. (H) Example of A-to-B compartment switching of chromosome 4. The left panel shows the first eigenvector (compartments; Y-axis) along the genomic position in megabases (X-axis). Each row corresponds to one independent biological replicate of the indicated stage. A compartments are depicted in yellow, and B compartments are shown in blue. The right panel corresponds to an 11.4-Mb zoomed region of chromosome 4 that contains the Klf4 locus. (I) Aggregate genome-wide contact profiles centered on TAD borders defined in stage 0. Data are the log2 ratio of observed and expected contacts in 10-kb bins, pooling biological replicates.

We obtained high-quality maps of the 3D genome organization across all stages (Supplemental Fig. S1E), which allowed us to examine the segregation of active (A, gene-rich) and inactive (B, gene-poor) compartments (Lieberman-Aiden et al. 2009; Imakaev et al. 2012). Principal component analysis (PCA) of the eigenvectors of all autosomes revealed that PML-RARα induced genome-wide, cumulative changes in A/B compartments (Fig. 1B,C). Overall, 7.1% of the genome changed compartment at some point during APL transformation, with 3% of the genome stably switching from the A to B compartment, and 1.1% switching from B to A (Fig. 1D,E). A greater proportion of switching events from stage 0 to III occurred from the A to B compartment (Fig. 1D,E; Supplemental Fig. S1F); this is in agreement with the known role of PML-RARα as a transcription repressor (Di Croce et al. 2002; Segalla et al. 2003; Villa et al. 2004, 2006, 2007; Carbone et al. 2006; Morey et al. 2008; Saumet et al. 2009; Martens et al. 2010; Subramanyam et al. 2010; Saeed et al. 2011, 2012; Cole et al. 2016). We found a cumulative total of 724 genes in bins that stably switched from compartment A to B, and 223 genes in bins that switched from B to A (Fig. 1F). The gene set that switched from A to B was enriched for genes associated with MAPK signaling (including those encoding Ras, Rap1, and cAMP), immune signaling via TGF-β, cellular differentiation, apoptosis, and transcriptional misregulation in cancer, as shown by KEGG analysis (Fig. 1G). In addition, this A-to-B gene set was significantly enriched for SMAD4 targets, as shown by ChEA analysis (adjusted P = 0.00045) and Polycomb targets (enriched in H3K27me3; adjusted P = 0.00034); specific genes included Mef2c, Flt3, Hmga2, Maf, Pax7, Met, Igf1, Wnt16, Aff1, Ptk2, Runx2, Rel, and Prom1. Of note, the A-to-B gene set also included several genes with known roles in leukemia development at compartment boundaries, such as Klf4, Setbp1, Efl1, and Hhip (Fig. 1H; Supplemental Fig. S1H; Alder et al. 2008; Kobune et al. 2012; Faber et al. 2013a; Schoenhals et al. 2013; Huang et al. 2014; Filarsky et al. 2016; Morris et al. 2016; Seipel et al. 2016; Makishima 2017; Park et al. 2019b; Tan et al. 2019). In contrast, the set of genes that switched from the B to A compartment was enriched for immune system processes, as shown by KEGG pathway analysis (Supplemental Fig. S1G); this included Il33, Il9r, Klf5, and Mcm10. When examining Hi-C interactions with intra-TAD regions, we observed minimal changes in TAD border strength (Fig. 1I). Overall, our data showed that PML-RARα expression induced a dynamic reorganization of the genome, affecting a large set of actively transcribed regions of the genome and causing their interaction patterns to switch toward those in the inactive chromatin compartment.

PML-RARα promotes dynamic changes in gene expression that are linked to changes in genome topology

Our in situ Hi-C data indicated that PML-RARα reorganized long-range interactions in a cumulative manner across the genome and was potentially accompanied by dynamic transcription alterations. To confirm this hypothesis, we performed RNA-seq on independent biological replicates harvested in duplicate at all stages. Based on PCA of the RNA-seq data sets, we observed a trajectory of transcriptome alterations concurrent with PML-RARα expression; of note, full-blown leukemias (stage IV) showed extensive transcriptome reprogramming as compared with the other stages (Supplemental Fig. S2A). We used two differential gene expression analyses to identify (1) genes significantly deregulated during APL transformation with respect to stage 0 (control) cells, and (2) genes uniquely deregulated (i.e., excluding genes that were also deregulated at other stages) at each stage of APL progression as compared with stage 0, which identified genes that are “transiently” altered during the kinetic analysis (Supplemental Table S1). The first analysis revealed an increasing number of significantly deregulated genes (Q-value < 0.05) from stage 0 during leukemic transformation (Fig. 2A). In addition, these differentially regulated genes progressively increased or decreased in expression along the time course (Fig. 2B), suggesting that PML-RARα expression induced early alterations in expression (e.g., at stage I) that were maintained during APL transformation. The second analysis revealed that a relatively small subset of genes was uniquely up-regulated or down-regulated at early stages of APL transformation, and that a larger number of genes was transcriptionally deregulated specifically at stages III and IV (Supplemental Fig. S2B; Supplemental Table S1). These results put forward the hypothesis that early expression of PML-RARα impaired the expression of a relatively few genes encoding for key hematopoietic TFs, which subsequently altered the transcriptional landscape genome-wide. Indeed, several genes encoding TFs or enzymes were either significantly up-regulated (e.g., Bcl2, Bmp2, Hes1, Mycn, Twist1, and Id2) or down-regulated (e.g., Cdh1, Lef1, Rxrb, and Rarg) at stages I and II. Overall, more genes were found to be down-regulated than up-regulated (Fig. 2A; Supplemental Fig. S2B), in line with previous reports of PML-RARα driving transcriptional repression (Morey et al. 2008; Gaillard et al. 2015; Li et al. 2018; Wang et al. 2020).

Figure 2.

PML-RARα promotes dynamic changes in expression of gene pathways that control cell cycle progression, immune signaling, and DNA repair. (A) Number of differentially expressed genes (DEGs) obtained by DESeq2 (Q-value > 0.05) at each of the indicated stages, using stage 0 as baseline. (B) Heat maps showing unsupervised clustering of expression levels (Z-score values) of genes at each stage, using stage 0 as reference. The left heat map depicts all genes that were up-regulated with respect to stage 0, and the right heat map depicts all down-regulated genes. (C) GSEA signatures of DEGs at each indicated stage with respect to stage 0. Gene expression signatures related to cell cycle control and p53/DNA damage repair were positively enriched during the time course, while signatures related to immune signaling were negatively enriched. (D) Heat map depicting the dynamic gene expression alterations (log RPKM values) of key hematopoietic transcription factors and leukemia-associated genes at each stage. (E, left panel) Expression of genes in regions that were in compartment A at stage 0 and switched to compartment B at another stage (724 genes). P = 2 × 10−4 between stages 0 and II; P = 1.1 × 10−8 between stages 0 and III. (Right panel) Same as the left panel, but for genes in the B compartment at stage 0 that switched to the A compartment at another stage (223 genes). P = 0.79 between stages 0 and II; P = 0.11 between stages 0 and III. Genes in bins that switched compartments more than once during the kinetic assay were excluded. P-values were computed using the Wilcoxon test (two-sided).

We next performed gene ontology and KEGG pathway analyses to dissect the pathways perturbed by PML-RARα. Genes with an increased expression in early stages showed enrichment in pathways related to MAPK signaling, regulation of cell proliferation and/or adhesion, or pathways in cancer. In turn, genes with a reduced expression in early stages were mostly related to hematopoietic cell lineage or immune response (Supplemental Table S1; Supplemental Fig. S2C). GSEAs of genes during APL progression revealed increased expression of genes involved in pathways related to cell cycle (E2F targets, G2M checkpoint, and mitotic spindle) or DNA repair, and decreased expression of genes involved in apoptosis and immune signaling pathways (Fig. 2C). In the leukemic stage (IV), a large number of genes was deregulated with respect to their status in stage 0 (Supplemental Fig. S2D); however, a large proportion of these genes already showed altered expression at stage III (Supplemental Fig. S2E). Included in the top transcriptionally deregulated genes were genes that encode TFs or enzymes that play fundamental roles in myeloid differentiation (Rosenbauer and Tenen 2007) and HSPC function, including Gata2, Cebpα, Bcl2, Hoxa10, Irf8, Myc, Spi1, and Klf4 (Fig. 2D). These results were validated in independent biological samples using qRT-PCR (Supplemental Fig. S2F).

Overall, our Hi-C and transcriptomic data indicated that cells expressing PML-RARα undergo progressive and profound alterations in genome architecture that may be correlated to changes in gene transcription. To confirm this hypothesis, we examined the transcriptional status of genes located in bins that switch compartments during APL transformation. Indeed, expression of genes in regions that switched from the A compartment at stage 0 to the B compartment at any later stage was significantly down-regulated (Fig. 2E). Although not statistically significant, the opposite trend was observed for genes that switched from the B to A compartment. Genes that stably switched from one compartment to the other were frequently found at pre-existing boundaries between A and B compartments at stage 0 (22.28% of the total genes are located at ±100 kb from A/B boundaries, whereas this proportion increases to 66.58% and 66.97% for stable A-to-B and stable B-to-A genes, respectively; P-value < 2.2 × 10−16). Together, our data indicated that PML-RARα induced early chromatin topological alterations, and in particular switched the interaction patterns of active regions of the genome to inactive chromatin compartments, which correlated with transcription repression.

PML-RARα induces epigenomic alterations at enhancers, which correlate with changes in expression of nearby genes

PML-RARα induces important expression changes of key genes involved in hematopoietic stem cell function and differentiation (Fig. 2D; Supplemental Fig. S2F; Tan et al. 2021). Given that transcriptional output is controlled by the activity of distal regulatory enhancers, we hypothesized that PML-RARα influences transcription of these genes by modulating enhancer activation. To comprehensively examine alterations of enhancer activity during APL progression, we collected samples at all experimental stages and performed ChIP-seq to map the genome-wide distribution of active enhancers (H3K4me1 and H3K27ac), active promoters (H3K4me3 and H3K27ac), and Polycomb-mediated repression (H3K27me3). These experiments revealed interesting patterns in both the number of peaks and their genomic distribution (Supplemental Fig. S3A,B). First, while the global number of H3K4me1-enriched regions was very similar between stages, the number of regions enriched in H3K4me3 and H3K27ac—hallmarks of active promoters—substantially decreased during APL progression (stages III and IV) (Supplemental Fig. S3A). Second, regions decorated by H3K27me3 increased along the four stages (Supplemental Fig. S3A), suggesting that PML-RARα led to a cumulative repression of the epigenome. Third, the reduced H3K27ac and increased H3K27me3 levels mostly occurred outside promoters of coding genes (i.e., in intergenic and intragenic regions), suggesting that PML-RARα had a primary role in epigenetic repression of putative enhancers (Supplemental Fig. S3B). Following these results, we next mapped the dynamic loss of enhancer activity during leukemic transformation (Fig. 3A). We identified 27,341 active enhancers at stage 0, of which 5%, 17%, 20%, and 21% lost H3K27ac at stages I, II, III, and IV, respectively (Fig. 3A; Supplemental Fig. S3C). We observed a striking progressive reduction of H3K27ac levels at enhancers with reduced levels in one stage during the subsequent stages. For example, enhancers with reduced H3K27ac levels at stage II continued to lose H3K27ac levels at stages III and IV (Fig. 3A). Importantly, loss of H3K27ac was accompanied by a gain in the repressive histone mark H3K27me3 (Fig. 3B). Interestingly, motif analysis of sequences of these enhancers revealed significant hits for the PU.1/SPI1 and the myeloid-determining transcription factor GFI1B (Supplemental Fig. S3D). These data align with previous literature (Wang et al. 2010; Tan et al. 2021) and confirm the role of the PML-RARα–PU.1 and PML-RARα–GFI1B axes during APL progression. Moreover, the KLF4 motif was enriched at enhancers that are inactivated during leukemia progression, suggesting that KLF4 down-regulation might be one of the key events that could induce decommissioning of enhancers at a later time point, although some of the observed changes might be indirect. Together, our results indicate that PML-RARα induced a vast reprogramming of the epigenome that involved the repression of active enhancers concomitant with a gain of Polycomb-mediated repression.

Figure 3.

PML-RARα induces epigenomic alterations at distal regulatory elements correlated with changes in expression of nearby genes. (A) Intensity of H3K27ac signal at 27,341 stage 0 active enhancers (non-TSS regions with overlapping H3K4me1 and H3K27ac peaks) that lost H3K27ac during the time course. Box plots below correspond to the normalized H3K27ac ChIP-seq signal intensity of enhancers lost at each stage. For H3K27ac, P = 2.2 × 10−16 between stages 0 and I, between stages 0 and II, between stages 0 and III, between stages 0 and IV. For H3K27me3, P = 0.33 between stages 0 and I, P = 2.2 × 10−16 between stages 0 and II, P = 2.2 × 10−16 between stages 0 and III, and P = 2.2 × 10−16 between stages 0 and IV. (B) Intensity of H3K27me3 signal at the same enhancers shown in A. Box plots show the normalized H3K27me3 ChIP-seq signal intensity of enhancers lost at each stage. (C) UCSC genome browser snapshots of differential H3K27ac (purple) and H3K27me3 (green) ChIP-seq profiles at promoters and putative distal regulatory elements of the indicated genes. Each row corresponds to the ChIP-seq signal intensity at each indicated stage subtracted from the signal at stage 0 as baseline. (D) Expression of genes within 5 kb around active enhancers at stage 0 that are lost in stage I (1464 enhancers; left graph), stage II (4668 enhancers; middle graph), or stage III (5406 enhancers; right graph). P-values were computed using Wilcoxon test (two-sided). (E) Dynamic changes of overall contact enrichment (intra-TAD) of promoters depending on activation status from stage 0 to stage III are as follows: Active (blue dots) maintained H3K27ac in both stages, inactive (red dots) were not marked by H3K27ac in either stage, gain (purple dots) gained H3K27ac at stage III, and loss (green dots) lost H3K27ac at stage III. Contact enrichments were measured as log2 of observed contacts over expected (log2 Obs/Exp) and were corrected against background. (F) Box plots showing the changes in H3K27ac, H3K27me3, and RNA levels per TAD for the top and bottom 10% of TADs with higher changes in domain score.

To closely examine the dynamic alterations of the epigenome during leukemic transformation, we next subtracted the normalized signal intensity of H3K27ac and H3K27me3 at each stage from the baseline signal at stage 0; we then inspected the regulatory landscape near key hematopoietic transcription factors. We observed that PML-RARα expression induced a progressive loss of H3K27ac at enhancers near Klf4 and Spi1, which became transcriptionally repressed during APL transformation (Figs. 2D, 3C). In contrast, Gata2 and its putative enhancers showed progressively increased H3K27ac and reduced H3K27me3 levels (Fig. 3C), in line with its increased expression (Fig. 2D). These examples suggested that the epigenetic alterations induced by PML-RARα at enhancers were associated with changes in nearby gene expression. To address this question genome-wide, we examined the levels of expression of genes located within 5 kb from enhancers that presented a significant decrease in the H3K27ac levels at each stage (Supplemental Fig. S3C). Our data confirmed that loss of enhancer activity correlated with a significant decrease in expression of nearby genes (Fig. 3D).

Next, we used Hi-C to examine whether the overall physical contacts within the same TAD (topologically associating domain; i.e., contacts with other promoters or enhancers) were affected in promoters that lost or gained H3K27ac during APL transformation. Notably, promoters that had decreased levels of H3K27ac—and thus had become repressed—showed decreased contacts during the early phases of leukemic transformation (Fig. 3E). In contrast, activated promoters with increased H3K27ac levels showed the opposite trend, whereas the contacts of stably active or inactive promoters were maintained (Fig. 3E). Examples of intra-TAD reorganizations for a repressed gene (KLF4) and an activated gene (GATA2) are shown in Figure 4 and Supplemental Figure S4, respectively. These results suggest that changes in intra-TAD interactions may be required for transcriptional activation but not for transcriptional deactivation. To generalize those observations, we ranked TADs according to their changes in domain score, which reflect internal reorganization and compartmentalization of TADs (Krijger et al. 2016; Stadhouders et al. 2018) between stage 0 and stage III. The 10% of TADs with a higher increase in domain score at stage III showed increased levels of H3K27ac (Fig. 3F, left panel), reflecting intra-TAD reorganization and establishment of regulatory contacts in TADs that become active, confirming previous observations (Krijger et al. 2016; Stadhouders et al. 2018). Changes in domain score did not correlate with changes in H3K27me3 levels, and both TADs with a higher increase or decrease in gene expression showed an increased domain score at stage III, indicating a link between TAD reorganization and gene expression changes and suggesting complex reorganization of TADs upon PML-RARα expression. Collectively, our data indicated that PML-RARα promoted extensive epigenomic reprogramming of enhancer regions and induced transcriptional changes by rewiring promoter–promoter and promoter–enhancer contacts, thus affecting genes that encode for critical regulators of hematopoietic differentiation.

Figure 4.

The Klf4 locus undergoes progressive changes in long-range interactions driven by PML-RARα expression. (A, top panels) Hi-C interaction matrices showing normalized interaction counts at the region between 54 and 57 Mb of chromosome 4. Each stage is indicated by the top gray bar. (Middle panels) Hi-C interaction matrices following subtraction of the signal at stage 0. Blue depicts interactions that decreased after stage 0, and brown depicts interactions that increased after stage 0. Red arrows indicate the location of the Klf4 locus. (Bottom panels) Same as the middle panels, but zooming in at the region between 55 and 56 Mb of chromosome 4, depicting a progressive loss of interactions of Klf4 with downstream genomic elements. (B) Differential matrix at 5-kb resolution, showing normalized interaction counts at the region between 55 and 56 Mb of chromosome 4 at stage III after subtraction of stage 0 signal. Blue depicts interactions that decreased at stage III, and brown depicts interactions that increased at stage III. The insulation score track is shown below for stages 0 and III. RNA-seq tracks (black) show a progressive decrease of Klf4 gene expression. Differential H3K27ac (purple) and K3K27me3 (green) ChIP-seq tracks show a sequential loss of H3K27ac signal and gain of H3K27me3. A/B compartment bins showed progressive compartment switching of the Klf4 locus from A (orange) to B (blue). RefSeq genes of this genomic region are shown at the bottom.

The Klf4 genomic locus undergoes progressive rearrangement of long-range interactions during PML-RARα-induced transformation

KLF4 is a master hematopoietic transcription factor that acts as a tumor suppressor in leukemia by activating the expression of genes that promote myeloid differentiation, apoptosis, and cell cycle arrest (Feinberg et al. 2007; Alder et al. 2008; Huang et al. 2014; Filarsky et al. 2016; Morris et al. 2016). Our data showed that the Klf4 locus underwent an extensive regulatory reprogramming during APL transformation with an A-to-B compartment switch (Fig. 1H) and transcriptomic and epigenomic repression (Figs. 2D, 3C; Supplemental Fig. S2F). To investigate whether this reprogramming was accompanied by alterations in long-range interactions, we inspected the temporal changes in interactions centered around the Klf4 gene. The Klf4 locus is located at the boundary between two well-defined TADs (Fig. 4A). During APL transformation, we observed a progressive loss of long-range interactions of the Klf4 locus with the downstream TAD, with multiple interaction loops profoundly decreased at stage III as compared with stage 0 (Fig. 4A, middle Hi-C map). Simultaneously, the insulation between the two TADs (Crane et al. 2015) decreased during the kinetic (Fig. 4B), leading to increased interactions with the upstream TAD linked to the change of compartment of the Klf4 locus. Interestingly, such changes in long-range interactions were accompanied by epigenetic and transcriptional reprogramming, as shown by decreases in H3K27ac levels and gene expression toward the downstream TAD (Fig. 4B).

Similar topological alterations were observed in the Etv1 locus, a gene that is recurrently transposed in acute leukemia (Sacchi et al. 1986) and that was also repressed during APL progression (Supplemental Fig. S4A, left panel). The Etv1 gene progressively lost contacts, H3K27ac signal strength, and transcriptional output, concomitant with a gain of H3K27me3 levels. We also inspected the pattern of interactions around the Gata2 locus as an example of a gene encoding a master regulator of myeloid differentiation that is up-regulated in APL (Fig. 2D; Zhang et al. 2008a; Li et al. 2018). The Gata2 gene showed conspicuously increased contacts with neighboring genes in a region of ∼0.5 Mb (Supplemental Fig. S4A, right panel). In addition, we observed an enrichment in H3K27ac levels and transcriptional output at the Gata2 locus during APL transformation (Supplemental Fig. S4A, right panel).

Given the remarkable topological rearrangements observed at the Klf4 locus, we next sought to identify potential cis-regulatory elements that interacted with the Klf4 locus and to examine their contact profiles during APL transformation. To this end, we generated virtual 4C-seq maps at stage 0 and stage III that were centered at the Klf4 locus (Fig. 5A). These maps revealed that, in normal hematopoietic stem/progenitor cells, the Klf4 promoter had strong interactions with potential enhancers located at ∼119, 198, and 274 kb upstream of the promoter (Fig. 5A). These interactions were markedly reduced at stage III of the time course, while interactions downstream from the Klf4 promoter tended to increase. Notably, the Klf4 putative enhancers identified at stage 0 were enriched in H3K4me1 and H3K27ac in normal cells, and these marks were reduced in stage III. In addition, the region spanning the +119-kb enhancer showed a conspicuous increase in the levels of the Polycomb-repressive mark H3K27me3 (Fig. 5A; Di Carlo et al. 2019). The virtual 4C-seq map around the Etv1 locus also confirmed that transcriptional repression was accompanied by loss of contacts between the Etv1 promoter and its downstream enhancers, which decreased their activity from stage 0 to stage III, as shown by the loss of active histone modifications (Supplemental Fig. S4B). In contrast, the Gata2 locus (which is up-regulated during APL progression) showed a marked gain in interactions both upstream of and downstream from the gene, including at Gata2 putative enhancers (Supplemental Fig. S4C). Furthermore, we found a marked decrease in the Polycomb-mediated repressive mark H3K27me3 around Gata2. Altogether, our data showed that PML-RARα expression induced extensive rearrangements in long-range interactivity at loci encoding for master hematopoietic transcription factors.

Figure 5.

Klf4 overexpression restores the normal function of PML-RARα-expressing HSCs. (A) Virtual 4C-seq signal around the Klf4 genomic region. The top panel illustrates overall contact profiles at the gene promoter at stage 0 (red line) and stage III (black line) for the region between 54.5 and 56.5 Mb of chromosome 4. The zoomed-in panel corresponds to the region between 55 and 56 Mb of chromosome 4. ChIP-seq tracks of the indicated histone modifications at stage 0 and stage III are shown. Shadowed regions highlight the regions spanning the Klf4 gene and the putative enhancers upstream of the Klf4 promoter (+119 kb, +198 kb, and +274 kb from the promoter). (B) Immunophenotyping analyses of lineage-negative bone marrow cells overexpressing KLF4-GFP, PML-RARα-hCD4, or both (KLF4 OE + PML-RARα). Cells were sorted and analyzed by flow cytometry using the indicated cell surface markers. P-values were calculated using a Student's t-test between PML-RARα and KLF4 OE + PML-RARα conditions. (**) P = 0.016, (*) P = 0.003. (C) Quantification of colony-forming units (CFUs) during four consecutive replatings of sorted cells overexpressing KLF4-GFP, PML-RARα-hCD4, or both (KLF4 OE + PML-RARα).

Klf4 overexpression inhibits self-renewal and promotes differentiation of PML-RARα-expressing cells

Our data indicated that PML-RARα progressively down-regulated Klf4 expression by remodeling long-range interactions at the Klf4 locus. Both tumor suppressor and oncogenic roles have been reported for Klf4 in the context of APL (Feinberg et al. 2007; Alder et al. 2008; Morris et al. 2016; Lewis et al. 2021). Klf4 expression appears to be lower in APL patient samples carrying the t(15;17) translocation as compared with other AML subtypes or healthy bone marrow (Supplemental Fig. S5A). To determine whether down-regulation of this TF contributes to leukemic phenotypes, we examined whether ectopic Klf4 expression was able to reverse the phenotypic alterations driven by PML-RARα. To this end, we generated lineage-negative (Lin−) cells expressing PML-RARα only, Klf4 only, or both using a retroviral strategy (Supplemental Fig. S5B). Expression of PML-RARα arrested cellular differentiation, as shown by the decreased frequency of CD11b+ cells and the increased frequency of cKit+ cells, as compared with an empty vector control (Supplemental Fig. S5C). Klf4 overexpression alone did not significantly change the frequency of cKit+ or CD11b+ cells as compared with empty vector control (Supplemental Fig. S5C). However, ectopic Klf4 expression in cells simultaneously expressing PML-RARα substantially increased differentiation (Fig. 5B; Supplemental Fig. S5C). In addition, the self-renewal capacity of PML-RARα-expressing cells was completely abrogated when Klf4 was overexpressed (Fig. 5C; Supplemental Fig. S5D,E). These results were further supported by RNA-seq analysis of cells overexpressing KLF4 in the absence and presence of PML-RARα (Supplemental Fig. S6A,B). By comparing the transcriptome of cells co-overexpressing PML-RARα and KLF4 with the one from cells overexpressing PML-RARα only, we observed alterations in different processes that are essential for cell growth and leukocyte function. Interestingly, KLF4 expression in PML-RARα cells results in up-regulation of cellular senescence programs (Supplemental Fig. S6C). Together, these results suggested that Klf4 down-regulation, which is induced by PML-RARα, is indeed a leukemia-promoting event that can be reversed by ectopic Klf4 expression.

Discussion

We have shown that the expression of the chimeric protein PML-RARα in primary HPSCs induces a rapid and extensive remodeling of contacts genome-wide as well as reprogramming of both the epigenome and the transcriptome. We showed that this process is, at least partially, dynamic and continuous, impacting transcription and the enhancer landscape around genes that encode for key transcription factors, which control the differentiation and function of HPSCs. Among these alterations, we identified major changes and transcriptional repression at the Klf4 gene and neighboring enhancers and showed that ectopic overexpression of Klf4 restored the differentiation capacity of HPSCs that expressed PML-RARα. Our findings add to recent studies addressing the role of PML-RARα (Li et al. 2018; Wang et al. 2020; Tan et al. 2021) and of other oncofusion proteins, such as RUNX1-ETO (Ptasinska et al. 2019), in genome architecture by (1) using a primary cellular and animal model system that closely recapitulates clinical and morphological features of human APL, (2) providing the first temporal multiomics analysis of the alterations driven by the chimeric protein at both promoter and enhancer regions, and (3) identifying specific changes that occur at Klf4 putative enhancers and demonstrating a tumor suppressor role of this transcription factor in promyelocyte leukemic transformation. While our experimental system has been validated extensively and is known to maintain two key functional aspects of PML-RARα function (inhibiting differentiation and enhancing self-renewal), it is important to acknowledge that using in vitro cultured cells to characterize early stages of PML-RARα function may not provide the complete view of the molecular events or genetic mutations required to develop leukemogenesis in vivo.

A few studies have mapped PML-RARα occupancy genome-wide using human APL cell lines or APL blasts (Hoemme et al. 2008; Martens et al. 2010; Mikesch et al. 2010; Wang et al. 2010, 2020; Singh et al. 2018; Tan et al. 2021). Despite multiple attempts, we were unable to map the fusion protein in primary mouse HSPCs by ChIP-seq; however, our temporal analysis allowed us to identify early (and potentially direct) alterations in the topology, epigenome, and transcriptome driven by PML-RARα expression. Such earlier changes are more likely driven by the direct effect of the fusion protein and its primary targets, while stage IV alterations might reflect the occurrence of additional genetic alterations and other potentially cooperative effects that promote the fully transformed leukemic phenotype. Among the early alterations, we focused on those occurring at the Klf4 locus, which we characterized in depth. Furthermore, we cross-validated our RNA-seq data by overlapping our differentially expressed gene lists with known PML-RARα target genes in the NB4 APL cell line (Tan et al. 2021), and with known target genes of the TF PU.1 in HSCs (Pundhir et al. 2018), which are reported to be coregulated by PML-RARα (Wang et al. 2010; Yang et al. 2014). These analyses showed that ∼30% of differentially expressed genes between stages 0 and II, and between stages 0 and III, are bona fide PML-RARα targets in the patient-derived NB4 cells. We also found that >52% or 60% of up-regulated or down-regulated genes, respectively, were PU.1 targets in HSCs. In addition, genes down-regulated during the kinetic analysis showed significant enrichment in PRC2 (Suz12) as well as SMAD4 targets (adjusted P-values 1 × 2.5–8 and 1 × 1.05−6, respectively), confirming previous reports (Lin et al. 2004; Villa et al. 2007; Morey et al. 2008) and further validating our analyses. Interestingly, up-regulated genes were enriched not only for Suz12 and SMAD4 targets (adjusted P-values 10 × 4.88−8 and 10 × 9.4−3, respectively), but also for Gata2 targets (adjusted P-value 10 × 4.88−8). Expression of Gata2 increased very early during the kinetic analysis (Fig. 2D; Supplemental Fig. S4A,C) and in human APL patient samples (Sukhai et al. 2008; Katerndahl et al. 2021). Notably, Gata2 up-regulation was recently reported to suppress PML-RARα-induced leukemic transformation, indicating that, in addition to Klf4, Gata2 modulation might contribute to suppressing malignant transformation (Katerndahl et al. 2021).

The temporal resolution of our data revealed that the previously reported repressive functions of PML-RARα occurred in a progressive manner as cells underwent transformation (Figs. 2B,D, 3A). We found that the most pronounced alterations in long-range interactions, the epigenome, and the transcriptome occurred from stage III to stage IV (full-blown leukemia). These observations integrated other studies (Gaillard et al. 2015) that have shown that the initial changes in gene expression driven by PML-RARα are relatively subtle and related to metabolism, cell cycle, and DNA damage response signatures (Fig. 2C), which may be insufficient to terminally arrest differentiation. This model is further reinforced by DNA methylation analyses that report only modest epigenome alterations at early stages of APL development (Schoofs et al. 2013; Gaillard et al. 2015). It is important to note that stage IV leukemic blasts were isolated from live animals and thus were likely to be influenced by microenvironmental cues in the bone marrow. Furthermore, the development of full-blown APL blasts requires secondary mutations, which could further contribute to the divergence in our stage IV samples. Future studies are warranted to identify the contribution of secondary lesions and the bone marrow microenvironment in the cellular phenotypes observed during later stages of APL development.

Here, we focused on uncovering early alterations occurring at regulatory enhancers encoding for key hematopoietic TFs, including Klf4, that can subsequently have major impacts in the transcriptome, genome architecture, and methylome of APL blasts. The role of Klf4 in hematopoietic malignancies has remained controversial. Our study sheds new light on this issue by demonstrating a tumor suppressor role of Klf4 in the context of APL, showing that the gene is progressively down-regulated during APL progression and that its ectopic overexpression counteracts oncogenic transformation. Given that Klf4 is required for mesoderm lineage commitment (Aksoy et al. 2014), we speculate that Klf4 down-regulation rewires gene regulatory networks that promote HSPC differentiation, thus contributing to leukemogenesis. We found that the cis-regulatory landscape within the Klf4 locus substantially changed its pattern of long-range interactions and histone modifications concomitant with a reduction in Klf4 expression. The increase in long-range interactions observed around the Klf4 locus could be triggered by the switch from the A to B compartment: As Klf4 is progressively embedded into a larger B compartment, B-to-B interactions might be facilitated. Although treatment of APL with all-trans retinoic acid in combination with chemotherapy results in remission in >90% of patients, our data suggest that ATRA combined with enhanced KLF4 expression may open a novel avenue of therapeutic intervention. Indeed, we have observed synergistic effects of ATRA in inducing apoptosis, enhanced G1-phase arrest, and differentiation of cells coexpressing PML-RARα and KLF4 (data not shown), confirming the tumor suppressor role of KLF4 overexpression and its molecular effects in the context of APL.

Our work delineates the dynamic mechanisms whereby the oncogenic TF PML-RARα builds a network of chromosome interactions that repress transcription of master hematopoietic regulators. We propose that the dynamic changes in the genome architecture mediated by PML-RARα may serve as a general paradigm for other oncogenic proteins that act as transcriptional repressors, bringing new light to the molecular mechanisms by which these transcriptional repressors drive malignant transformation, and possibly leading to the identification of novel transformative therapeutic strategies.

Materials and methods

Murine APL model, bone marrow harvest, and cell culture

Bone marrow lineage-negative hematopoietic stem/progenitor cells from 8- to 10-wk-old female 129SvEv mice were harvested and infected with high-titer retroviruses expressing either an empty PINCO-3xFlag vector or a PINCO–PML-RARα-3 × Flag vector carrying human PML-RARα. PINCO plasmids expressing human PML-RARα from the 5′ viral long terminal repeat (LTR) and GFP from an internal promoter (cytomegalovirus [CMV]) were described previously (Grignani et al. 1998; Minucci et al. 2002) and were modified by cloning three copies of a Flag tag (69 bp) at the C-terminal of the human PML-RARα sequence. GFP+ cells transformed with empty vector or PML-RARα-3xFlag vector were sorted by FACS and correspond to stage 0 and stage I, respectively. Stage I cells were then plated in methylcellulose supplemented with cytokines and stem cell factor and serially replated for 2 wk (stage II) and 4 wk (stage III). The GFP+ cells that were passaged on methylcellulose were not resorted at each passage. In parallel, ∼1 million GFP+ PML-RARα-3xFlag transduced lineage-negative cells (stage I) were transplanted via tail vein injection into lethally irradiated (9 Gy) syngeneic mice (129SvEv) as previously described (Minucci et al. 2002). The animals were monitored periodically for signs of disease and the presence of blasts as evaluated by complete blood counts (CBC) and peripheral blood smears. Leukemic mice were humanely euthanized, and leukemic blasts were isolated from the spleen (with >95% of leukemic cell infiltration) for subsequent experiments (stage IV).

Bone marrow lineage-negative cells were obtained and transduced as described previously (Minucci et al. 2002). Serial replating assays of GFP+ cells were performed by seeding 10,000 cells/well in methylcellulose medium (Methocult, Stem Cell Technologies M3434) and replating every 7 d. Flow cytometry analyses were performed by staining cells with antimouse CD11b (eBioscience 25-0112-82) and using BD FACSCalibur 2.

Animal handling was performed following Italian laws (D.L.vo 116/92 and subsequent additions), which enforce EU Council Directive 86/609/EEC of November 24, 1986, on the approximation of laws, regulations, and administrative provisions of the Member States regarding the protection of animals used for scientific purposes. Mice were housed according to guidelines from the Commission Recommendation 2007/526/EC, June 18, 2007. The protocol was approved by the Italian Ministry of Health (authorization October 2013).

Western blotting

Whole-cell lysates of 293T cells infected with empty PINCO-3xFlag vector, PINCO-PML-RARα (Minucci et al. 2002), or PINCO-PML-RARα-3xFlag were obtained using RIPA buffer containing protease inhibitors (Roche). Sixty micrograms of total protein was loaded per lane on an 8% SDS-PAGE. After blocking in 5% milk–TBST-1X, the following antibodies were incubated overnight at 4°C: anti-Flag (mouse monoclonal, 1:500; Sigma F1804) and anti-Tubulin (mouse monoclonal, 1:5000; Abcam ab7291). Immunodetection was performed using ECL.

In situ Hi-C experiments

For in situ Hi-C experiments, 5 million to 10 million cells of each stage were harvested at two independent biological replicates per stage. Cells were cross-linked for 10 min at room temperature with 1% formaldehyde and quenched during a 5-min incubation at room temperature with 125 mM glycine, followed by a 15-min incubation on ice and two washes with cold PBS; samples were then pelleted and frozen at −80°C.

In situ Hi-C libraries were generated as previously described (Rao et al. 2014) with minor modifications (Mas et al. 2018). Two biological replicates were sequenced for all stages, with one additional technical replicate for stage II, giving between 70 million and 400 million valid reads per replicate. Supplemental Table S2 summarizes the statistics and reads obtained for all in situ Hi-C samples.

RNA-seq and quantitative real-time PCR (qRT-PCR)

For RNA-seq experiments, 1 million to 3 million cells at each stage were resuspended in 350 µL of Qiazol (Qiagen) and frozen at −80°C. RNA was obtained by thawing the samples, adding an additional 350 µL of Qiazol, and using the miRNEAsy mini kit as recommended (Qiagen). After RNA extraction, contaminating genomic DNA was eliminated with DNase I digestion. Two independent biological replicates per stage were used to generate RNA-seq libraries.

RNA samples were quantified using Nanodrop, and RNA quality was evaluated with an Agilent Bioanalyser (RIN > 9.9). Total RNA (1 µL) was used to generate RNA-seq libraries with rRNA depletion using TruSeq stranded total RNA library preparation kit (Illumina RS-122-2201). Libraries were sequenced in a HiSeq 2000 (75-bp, paired-end reads) to obtain ∼300 million raw reads per sample.

For qRT-PCR, 0.5–1 µg of total RNA obtained from independent biological samples at each stage was converted to cDNA, and qRT-PCR was conducted using SYBR Green (LightCycler Roche) and the following primer sequences (5′ to 3′): Klf4 (Fwd-CGGGAAGGGAGAAGACA, Rev- GAGTTCCTCACGCCAAC), Spi1 (Fwd- GCGTGCAAAATGGAAGGGTT, Rev- GTGTGCGGAGAAATCCCAGT), Irf8 (Fwd-CAATCAGGAGGTGGATGCTT, Rev-AGCACAGCGTAACCTCGTCT), Myc (Fwd-CCTAGTGCTGCATGAGGAGA, Rev-TCCACAGACACCACATCAATTT), Flt3 (Fwd-ATCTCCGAGGGTGTTCCAGA, Rev-TGAACAGCTTGGTGCATTCG), Gata2 (Fwd-GCTTCACCCCTAAGCAGAGA, Rev-TGGCACCACAGTTGACACA), Gata1 (Fwd-ACGACCACTACAACACTCTGGC, Rev- TTGCGGTTCCTCGTCTGGATTC), c-Kit (Fwd-GATCTGCTCTGCGTCCTGTT, Rev-CTTGCAGATGGCTGAGACG), and Bcl2 (Fwd- GAACTGGGGGAGGATTGTGG, Rev- GGCCATATAGTTCCACAAAGGC). Rplp0 (Fwd- TTCATTGTGGGAGCAGAC, Rev-CAGCAGTTTCTCCAGAGC) was used as housekeeping control for Klf4, Gata1, c-Kit, and Bcl2. For those genes, final values were multiplied by 1000. For the rest of genes, β-actin (Fwd- GGCCCAGAGCAAGAGAGGTATCC, Rev-ACGCACGATTTCCCTCTCAGC) was used as housekeeping control.

ChIP-seq experiments

For ChIP-seq experiments, 5 million to 10 million cells were cross-linked as described above. Experiments were performed as previously published (Mas et al. 2018). Chromatin complexes were immunoprecipitated using anti-H3K27me3 (Millipore 07-449), anti-H3K4me1 (Abcam ab8895), anti-H3K4me3 (Diagenode C15410003), and anti-H3K27ac (Millipore 07-360). A small aliquot of ChIP DNA was used for ChIP-qPCR validations using primers of transcriptionally active and repressed genes (Nucleolin, Sox2, and Gapdh) to verify enrichment of the histone modifications. About 2–10 ng of ChIP or input DNA material was used to prepare ChIP-seq libraries following the NEBNext Ultra DNA library preparation kit for Illumina (NEB E7370L) as per the manufacturer's instructions. Final ChIP-seq libraries were size-selected to remove fragments <100 bp and then amplified for 10 PCR cycles. Libraries were sequenced on a HiSeq 2000 platform (Illumina) to obtain ∼30 million reads per library (50 bp, single end).

RNA-seq and ChIP-seq data analyses

RNA-seq replicate samples were mapped against the mm10 mouse genome assembly using TopHat (Trapnell et al. 2009) with the option –g 1 to discard reads that could not be uniquely mapped to just one region. DESeq2 (Love et al. 2014) was run to quantify the expression of every annotated transcript using the RefSeq catalog of exons and to identify each set of differentially expressed genes. Expression values shown in the box plots correspond to the averaged FPKMs across the two replicates in each stage. The rows of the heat maps of gene expression were scaled to have mean 0 and a standard deviation of 1 (Z-score). To define the set of unique differentially expressed genes (up or down), only genes reported to significantly change expression in a single stage as compared with stage 0 were included in the heat maps. Gene set enrichment analysis of the preranked lists of genes by DESeq2 stat value was performed with the GSEA software (Subramanian et al. 2005).

ChIP-seq raw reads were mapped against the mm10 mouse genome assembly using Bowtie (Langmead et al. 2009) with the option -m 1 to discard reads that did not map uniquely to one region. MACS (Zhang et al. 2008b) was run with the default parameters but with the shift size adjusted to 100 bp to perform the peak calling of ChIP-seq experiments. The genome distribution of each set of peaks was calculated by counting the number of peaks fitted on each class of region according to RefSeq annotations (O'Leary et al. 2016). “Promoter” wass the region within ±2.5 kb of the transcription start site (TSS), intragenic regions corresponded to the rest of the gene not classified as promoter, and the rest of the genome was considered to be intergenic. Peaks that overlapped with more than one genomic feature were counted in multiple categories. Active enhancers were defined by the presence of overlapping peaks of H3K4me1 and H3K27ac at stage 0 within intronic and intergenic regions. To define the set of active enhancers at stage 0 that were lost along the rest of stages (I, II, III, and IV), the ChIP-seq signal of H3K27ac was subtracted in stage 0 from the rest of H3K27ac profiles, and enhancers were identified in which the final value was less than one normalized read. To identify examples of enhancers gaining H3K27ac signal, the subtraction was inversely performed. The same procedure was used to determine the gain or loss of H3K27me3 in the same enhancer collection. Heat maps displaying the density of H3K27ac and H3K27me3 ChIP-seq reads around the center of each enhancer were generated by counting the number of reads for each individual enhancer and normalizing this value with the total number of mapped reads of the sample. The rows of the heat maps were scaled to have mean 0 and standard deviation 1 (Z-score), and plots were generated using SeqCode (Blanco et al. 2021).

In all analyses, we used release 68 of the RefSeq annotations (O'Leary et al. 2016) as provided by the UCSC genome browser on the refGene.txt file (Tyner et al. 2017). This RefSeq version contains 34,904 transcripts corresponding to 24,338 mouse genes. No preprocessing filtering steps were performed on this file. The UCSC genome browser was used to generate screenshots of the genomic landscape of selected genes (Tyner et al. 2017). Enrichr (Kuleshov et al. 2016) was used to perform gene ontology (GO), KEGG, and other functional analysis (such as ChEA) of the gene sets obtained from RNA-seq and genes in bins that switched A/B compartments. Supplemental Table S1 lists all differentially expressed genes and Enrichr results in each comparison. Graphical treatment and quantification of the ChIP-seq and the RNA-seq experiments was performed using SeqCode (Blanco et al. 2021).

Hi-C data analysis

Hi-C data were processed using an in-house pipeline based on TADbit (Serra et al. 2017). Reads were mapped according to a fragment-based strategy: Each side of the sequenced read was mapped in full length to the reference genome mouse December 2011 (GRCm38/mm10). TADbit filtering module was used to remove noninformative contacts and to create contact matrices as previously described (Serra et al. 2017). PCR duplicates were removed, and the Hi-C filters applied corresponded to potential nondigested fragments (extradangling ends), nonligated fragments (dangling ends), self-circles, and random breaks. Contact matrices were normalized for sequencing depth and genomic biases using OneD (Vidal et al. 2018). A and B chromatin compartment analysis was performed at 100-kb resolution as previously described (Lieberman-Aiden et al. 2009; Serra et al. 2017). Differential Hi-C matrices were computed from normalized Hi-C matrices at 5-kb resolution. Matrices of the specific regions were corrected for read coverage, a Gaussian filter was applied for noise reduction, and the difference between maps at stage III and stage 0 was plotted. Virtual 4C-seq profiles were generated from local coverage-normalized Hi-C matrices at 5-kb resolution, and Gaussian filter was applied for smoothing. The 5-kb bin containing the TSS of the gene of interest was used as the viewpoint. The domain score of consensus TADs was computed as previously described (Krijger et al. 2016; Stadhouders et al. 2018). TADs were ranked according to the ratio stage III/stage 0 of this score. The normalized level of H3K27ac, H3K27me3, and RNA per TAD was obtained using respective ChIP-seq and RNA-seq data sets. The ratio of the levels of these marks between stage III and stage 0 was compared between the 10% of TADs with higher changes (higher at stage 0 or higher at stage III).

Klf4 overexpression experiments

The KLF4 (mouse) and PML-RARα (human) cDNAs were cloned into MSCV-GFP and MSCV-hCD4 vectors (Addgene vector 35712) under the control of the EV promoter. The ecotropic phoenix packaging cell line was transiently transduced with the retroviral vectors cited above, and the retroviral supernatant was collected and filtered. Bone marrow lineage-negative cells were obtained from C57/Bl6 wild-type mice and transduced with retroviruses carrying either MSCV-GFP-KLF4, MSCV-hCD4-PE-PML-RARα, or both by two rounds of spinfection in nontissue culture-treated plates (Corning 351147) coated with retronectin (Takara T100A). Transduced Lin− cells were sorted and serially replated in methylcellulose medium (Methocult, Stem Cell Technologies M3434) by seeding 10,000 cells/well and replating every 7 d. Flow cytometry analyses were performed by staining cells with PE anti-hCD4 (BD Pharmigen 555347), APC antimouse Cd11b (BD 53312), and APC-fluo780 antimouse c-Kit (Invitrogen 47-1171-82) using FACSAria (BD).

Data availability

All sequencing data sets are available at GEO under accession number GSE151837.