RUNX1 C-terminal mutations impair blood cell differentiation by perturbing specific enhancer-promoter networks

Hematopoiesis is a vastly complex process, involving many signaling pathways, intricate transcriptional programs, in addition to further epigenetic and RNA splicing regulation. At the top of the hematopoietic hierarchy lie several master regulators which play critical roles throughout the proliferation and differentiation process. RUNX1 is among these master regulators and is required for definitive hematopoiesis.^1-4 

Pathogenic mutations occur in RUNX1. Many mutations have been detected within the DNA-binding runt homology domain (RHD) which disrupt protein binding to DNA and act as loss-of-function mutations. Outside RHD, mutations affect the regulatory regions located at the C-terminus of RUNX1.^5-7 These mutations are associated with sporadic myelodysplastic syndrome/acute myeloid leukemia (AML)^8,9 and an increased risk of AML transformation in patients with chronic myelomonocytic leukemia.¹⁰ Germ line mutations in this region have also been identified as pathogenic drivers of familial platelet disorder with associated myeloid malignancy.^11-13 

Although C-terminal mutations in RUNX1 have been proven to be pathogenic, the underlying mechanisms of this class of mutation remain poorly understood in hematopoietic disorders. We revealed that most C-terminal mutations are nonsense and frameshift mutations that are exempt from nonsense-mediated messenger RNA decay (NMD). Furthermore, to understand the mechanisms and impacts of these C-terminal mutations, we generated an isogenic knock-in human cell line model of RUNX1^R320∗ (RUNX1c notation was used in this study).^13,14 Our study established that RUNX1^R320∗ does not elicit NMD and produces a truncated protein. We examined the effects of RUNX1^R320∗ on transcription, DNA binding, and promoter-enhancer interactions using a combination of RNA-sequencing, chromatin immunoprecipitation (ChIP)-sequence, and Global RNA Interactions with DNA followed by deep sequencing (GRID-seq). Our analysis revealed a RUNX1^R320∗ transcriptional signature that is distinct from that induced by RUNX1 depletion. Interestingly, although we detected similar genome-wide binding between RUNX1 and RUNX1^R320∗, we identified extensive remodeling of enhancer-promoter networks in RUNX1^R320∗ cells. Analysis of RUNX1^R320∗ regulated enhancer-promoter pairs detected significant enrichment of FOXK2 motifs, suggesting a novel role for FOXK2 at enhancers in conjunction with RUNX1^R320∗. At the MYC locus, we found that RUNX1^R320∗ and FOXK2 both exhibited increased binding to hematopoietic MYC enhancers, whereas RNA-seq detected significant MYC upregulation in RUNX1^R320∗ cells. Collectively, we demonstrated that non-RHD RUNX1 mutants can produce proteins that do not act as simple loss of function and dysregulate hematopoiesis through distinct mechanisms and cofactor interactions.

A complete description of all methods is provided in the supplemental Methods section.

RNA extraction was performed in triplicate using TRIzol (Invitrogen #15596026) in accordance with the manufacturer’s protocol, and library preparation and sequencing were performed by Novogene. ChIP-seq samples were prepared as described previously¹⁵ with minor modifications, using an anti-RUNX1 antibody from Abcam (#23980). GRID-seq libraries of RUNX1 wild-type (WT) and RUNX1^R320∗ cells were prepared as previously described.^16,17 

As a transcription factor, RUNX1 binds DNA through the RHD which lies within the N-terminal region of the protein. Pathogenic mutations detected within RHD frequently disrupt DNA binding and act to prevent RUNX1 function. However, the effects of mutations other than those in RHD remain poorly understood. We sought to investigate mutations beyond RHD in more detail. To achieve this, we first assessed the distribution of hematopoietic RUNX1 mutations in the Catalogue of Somatic Mutations in Cancer (COSMIC) database.¹⁴ We found that a significant portion of RUNX1 mutations (27%, n = 387) lie outside the RHD and are distributed throughout the C-terminal region (Figure 1A [top]). A high proportion of C-terminal mutations were revealed to be nonsense or frameshift (78.6%, 304/387) (Figure 1B [left]). This was in stark contrast to our analysis of N-terminal/RHD regions with frameshift and nonsense mutations, which accounted for only 37.3% (384/1030) of mutations. RUNX1 germ line mutations in the RUNX1db (RUNX1 Database) also showed a similar trend to C-terminal nonsense and frameshift mutations.¹⁸ This led us to hypothesize that this disruption of the C-terminus through frameshift or truncation was linked to the pathogenicity of these mutations and warranted further investigation.

Mutations causing frameshifts and early termination codons typically elicit nonsense-mediated decay (NMD), whereas premature termination codons (PTCs) lead to transcript degradation. We reasoned that the high rate of C-terminal nonsense and frameshift mutations might either elicit NMD, causing RUNX1 haploinsufficiency or be exempt from NMD, thereby producing pathogenic protein variants. The mechanisms of NMD are well-defined at the transcript level and enable its prediction.^19-22 Briefly, only PTCs in the last exon and within 50 nucleotides of exon-exon junctions were not subjected to NMD. In the context of the RUNX1 transcript, PTCs beyond residue 305 of 480 (RUNX1c NM_001754.5) were predicted to be exempt from NMD (Figure 1A [bottom]). We found that most (232/304; 76.3%) C-terminal frameshift and nonsense mutations were predicted to be exempt from NMD and produce proteins (Figure 1B [right]). Together, our analysis demonstrated that the most common C-terminal RUNX1 mutations were frameshift or nonsense (78.6%), and most of these mutations were projected to produce proteins with truncated or novel C-termini (76.3%), representing a class of pathogenic RUNX1 mutations that are distinctly different from those found within RHD.

To study C-terminal RUNX1 mutations, we elected to generate a homozygous knock-in of the pathogenic RUNX1^R320∗ mutation (ClinVar: VCV000618862.12) using CRISPR-Cas9 technology in the K562 human leukemia cell line derived from a chronic myeloid leukemia patient at blast crisis (Figure 1C-D). RUNX1^R320∗ generates a premature stop codon predicted to be exempt from NMD mechanisms, representing the most C-terminal RUNX1 mutations in our analysis. The RUNX1^R320∗ transcript was confirmed to produce a protein product via immunoblot (Figure 1E). We observed a 2.37-fold increase in RUNX1^R320∗ protein relative to the wild-type as well as increased transcript expression (Figure 1F), which we hypothesize may be due to reported autoregulation of RUNX1.²³ The remaining RUNX family members, RUNX2 and RUNX3, have been reported to compensate for RUNX1 loss.^24-26 Although we detected a significant change in the transcript expression of RUNX3, the protein levels were barely detectable (supplemental Figure 1A), which enabled us to study the effects of RUNX1^R320∗ as the predominate RUNX protein in our model. These data demonstrate that the endogenous knock-in of RUNX1^R320∗ is not subject to transcript degradation and results in the production of a truncated RUNX1 protein, which is expressed at a level higher than that of the wild-type RUNX1.

RUNX1 plays roles throughout hematopoiesis and has been well described as an essential factor in MK development as well as platelet production and function.^27-29 To investigate whether RUNX1^R320∗ dysregulates K562 cells differentiation into MKs upon induction with 12-O-tetradecanoylphorbol-13-acetate (TPA),³⁰ both RUNX1 and RUNX1^R320∗ cells were treated with TPA and changes in cell morphology and surface markers were assessed after 48 hours. RUNX1 wild-type cells showed characteristic megakaryocytic differentiation after TPA treatment, including the appearance of large cells with lobated nuclei. However, RUNX1^R320∗ cells produced dysplastic MK-like cells that were smaller with dyslobated nuclei (Figure 2A). Undifferentiated K562 cells express CD235a, a maker presented on megakaryocytic-erythroid progenitors (MEPs), and lack the megakaryocytic lineage marker CD61. Upon TPA treatment, RUNX1 wild-type cells lost CD235a and gained CD61, demonstrating megakaryocytic differentiation (Figure 2B-C). RUNX1^R320∗ cells showed significantly less differentiation, confirming our results in Figure 2A. Together, our data show that endogenous expression of RUNX1^R320∗ results in a partial megakaryocytic differentiation block.

DNA damage, as an oncogenic driver, plays an important role in hematologic malignancies. RUNX1 aberrations have been shown to increase DNA damage.^25,31-34 We next sought to examine whether endogenously expressed RUNX1^R320∗ affects these DNA damage pathways. We treated RUNX1 and RUNX1^R320∗ cells with the DNA-damaging agent etoposide (ETOP) and assessed DNA damage sensing using γ-H2AX imaging (supplemental Figure 1B). RUNX1^R320∗ cells were significantly more sensitive to ETOP treatment compared with wild-type cells (Figure 2D); similar results were observed upon camptothecin (CPT)-induced damage (data not shown). We hypothesized that RUNX1^R320∗ induced DNA damage sensitivity may lead to increased apoptosis. Interestingly, extended ETOP treatment for >48 hours did not result in significantly increased apoptosis in RUNX1^R320∗ cells relative to that in wild-type cells (supplemental Figure 1C). Together, these results demonstrated that RUNX1^R320∗ cells become sensitized to DNA damage while evading cell death via apoptosis, suggesting that DNA damage sensitivity is a pathogenic attribute of RUNX1 C-terminal mutants.³² 

Next, we sought to investigate the impact of RUNX1^R320∗ on gene expression and disease pathways. RUNX1 influences gene expression through both direct DNA binding and protein-protein interactions with an abundance of cofactors (reviewed in¹). As a transcriptional master regulator, mutations in RUNX1 lead to aberrant gene expression. We hypothesized that RUNX1^R320∗ may uniquely disrupt transcription, as RUNX1^R320∗ retains the DNA-binding RHD. RNA-seq followed by principal component analysis (PCA) clearly indicated RUNX1^R320∗ samples generated different transcriptome signatures than that of wild-type RUNX1 samples (Figure 3A). Subsequent differential expression analysis revealed 1013 upregulated and 1663 downregulated genes (FDR ≤ 0.05; |log2FC| ≥ 1.5), demonstrating significant transcriptional reprogramming by RUNX1^R320∗ (Figure 3B).

To elucidate how RUNX1^R320∗ may alter transcription differently than RUNX1 RHD loss-of-function mutants, we compared RUNX1^R320∗ dysregulated genes with our previously generated RUNX1 knockdown RNA-seq data set in K562 cells³⁵ (Figure 3C; supplemental Figure 2A). Remarkably, most (74.63%) genes dysregulated in RUNX1^R320∗ were unique and not perturbed in RUNX1-depleted cells. Furthermore, among this small subset of commonly dysregulated genes, only 62.0% of these overlapping genes were dysregulated in the same manner (both up or downregulated) between the RUNX1^R320∗ and RUNX1 depleted data sets. These data demonstrate that RUNX1^R320∗ results in significant changes in transcription, dysregulating 2676 genes, which represents a unique transcriptional signature that differs from that induced by RUNX1 depletion.

Exploring the distinct RUNX1^R320∗ gene expression signature further, we performed an overrepresentation pathway analysis on both our RUNX1^R320∗ and RUNX1 KD data sets. In line with our observed phenotypic changes, we detected a significant enrichment of pathways related to MK and platelet function (Figure 3D). These pathways relate to the known role of RUNX1 in hematopoietic disease as well as MK differentiation and function.^27,29,36-38 We also performed gene set enrichment analysis (GSEA)³⁹ and revealed negative enrichment of RUNX1-regulated megakaryocytic and hematopoietic stem cell (HSC) differentiation gene sets in RUNX1^R320∗ cells (Figure 3E). We identified specific hematopoietic genes that were dysregulated in RUNX1^R320∗ cells using Reactome and gene ontology (GO) databases (supplemental Figure 2B). Furthermore, several oncogenic pathways, such as PI3K/AKT and MAPK signaling, were uniquely enriched in our RUNX1^R320∗ data set (Figure 3D). GSEA also uncovered the enrichment of MYC oncogenic signaling (Figure 3F). c-MYC, a well-established leukemogenic driver,^40-43 was significantly upregulated in RUNX1^R320∗ cells (Figure 3F-G). Together, these data suggest that RUNX1^R320∗ disrupts key MK and HSC differentiation pathways while upregulating oncogenic signaling through the dysregulation of both unique genes and known RUNX1 targets.

We demonstrated that RUNX1^R320∗ dysregulates hematopoietic gene expression, including genes directly related to disease phenotypes (Figure 3). As a master regulator, RUNX1 has been demonstrated to regulate gene expression through promoter and enhancer regulation as well as through chromatin remodeling.^30,42,44-46 We hypothesized that the transcriptional changes we observed in RUNX1^R320∗ cells might result from a combination of altered DNA binding and changes in cofactor interactions.

To explore changes in RUNX1 and RUNX1^R320∗ DNA binding, we performed ChIP-seq. Detailed peak annotation revealed that 39.8% of RUNX1 peaks were within promoter regions, 29.8% intronic, 20.6% intergenic, 5.5% exonic, and 4.3% in the 3’ untranslated region (UTR) and 5’UTR (Figure 4A). Our findings were consistent with previously reported RUNX1 ChIP-seq data sets^30,47 (supplemental Figure 3A-D). As RHD is retained in RUNX1^R320∗, we hypothesized that the loss of the C-terminus would dysregulate binding at a subset of RUNX1 sites through alterations in cofactor interactions and that DNA binding may also be changed as regions of the C-terminus have been reported to have autoinhibitory functions.^6,48 We first compared RUNX1 and RUNX1^R320∗ ChIP-seq data sets and uncovered similar binding annotation patterns: 42.7% of RUNX1^R320∗ peaks at promoter regions, 27.5% intronic, 19.3% intergenic, and 5.6% exonic, with the remaining 5’/3’ UTR and downstream accounting for 5.0% of the peaks (Figure 4B). Despite the increased expression of RUNX1^R320∗ (Figure 1E), we found that most peaks (50 596) were detected in both the RUNX1 and RUNX1^R320∗ data sets, demonstrating that both proteins exhibit similar genomic binding. We detected 1061 sites with significantly downregulated RUNX1^R320∗ binding, while 129 sites showed an increase in RUNX1^R320∗ presence (Figure 4C). These data suggested that RUNX1 and RUNX1^R320∗ bind similarly throughout the genome, displaying differential binding at a small subset of sites. These data also point toward further RUNX1^R320∗ mediated gene regulation through altered interactions with coactivators/corepressors, enhancers, and chromatin modifiers.

To investigate RUNX1^R320∗ transcriptional regulation at promoters we integrated our ChIP-seq and RNA-seq data sets to examine RUNX1 and RUNX1^R320∗ bound genes (Figure 4D). RUNX1^R320∗ promoter binding was correlated with gene expression. However, as we detected RUNX1^R320∗ binding beyond the promoter regions, we hypothesized that additional regulatory elements played a role in the RUNX1^R320∗ transcriptional changes that we observed. To annotate RUNX1^R320∗ binding, we divided the genome into 5 major categories: enhancers, promoters, transcribed regions, repressed regions, and heterochromatin, using publicly available K562 histone modification ChIP-seq data sets (Figure 4E).^49,50 Both RUNX1 and RUNX1^R320∗ differentially bound sites (Figure 4E “up/down”) and shared sites (Figure 4E “nc”) were enriched in promoter and enhancer regions. We conducted further analysis of RUNX1 motif density in differentially bound promoters and enhancers. Enhancers with altered binding were more strongly associated with the RUNX1 DNA-binding motif than with promoter regions (Figure 4F). These analyses demonstrate that the RUNX1 C-terminal region is required for binding to a subset of RUNX1 target sites and these dysregulated sites are strongly enriched in enhancer regions. Collectively, our data suggest a role for RUNX1^R320∗ in enhancer regions of transcription regulation, in addition to canonical promoter binding.

Enhancers have been shown to play critical roles in both normal and abnormal hematopoiesis^41,51-54 and we hypothesized that RUNX1^R320∗ may dysregulate critical enhancer-promoter connections based on our RNA- and ChIP-seq analyses. To uncover these connections, we used GRID-seq to map genome-wide RNA-DNA interactions and generated enhancer-promoter (E-P) network maps in RUNX1 and RUNX1^R320∗ cells.^16,17 GRID-seq detects RNA-DNA interactions using a bivalent linker to capture RNA and DNA molecules in close proximity. Nascent RNAs proximal to their endogenous promoter region, as well as any associated enhancers, were detected as enhancer-promoter pairs (Figure 5A).

We separated GRID-seq interactions into “local,” “cis,” and “trans” interactions. RNA is most likely proximal to the DNA it is transcribed from, typically the gene body, these interactions we define as “local.” Beyond the gene body, “cis” interactions are between RNA and DNA regions within the same chromosome and are most likely to represent enhancer-promoter pairs whereas “trans” interactions are interchromosomal. As shown in Figure 5B, local interactions were the most readily detected, followed by cis interactions. Trans-interactions are significantly more rare and typically weaker by orders of magnitude.^55,56 Generally, chromosomal interactions follow power law scaling, enabling mathematical modeling of the probability of DNA contacts, as described in detail by Lieberman-Aiden et al and others.^52,57 Our GRID-seq data sets successfully recapitulated these findings (Figure 5C) and allowed us to apply this model to GRID-seq detected local, cis, and trans interactions, ranking them and generating a Z-score scale in RUNX1 and RUNX1^R320∗ cells. Examining cis interactions between RUNX1 and RUNX1^R320∗ data sets, we detected 30 365 interactions unique to RUNX1 and 32 903 RUNX1^R320∗ specific interactions, with 52 089 occurring in both (Figure 5D). E-P pairs that were up and downregulated in RUNX1^R320∗ cells correlated with respective increases and decreases in gene expression (supplemental Figure 3E). Furthermore, our data revealed extensive interaction remodeling at differentially expressed hematopoietic and platelet gene loci such as KIT, DIAPH1, NFE2, and STIM1 (supplemental Figure 4A-D). Together, our GRID-seq data set, in combination with ChIP-seq and gene expression analysis, established that RUNX1^R320∗ broadly alters enhancer-promoter networks, leading to significant transcriptional dysregulation.

We hypothesized that RUNX1^R320∗ may remodel enhancer-promoter networks through cofactor interactions, either gained or lost upon truncation of the RUNX1 C-terminus. RUNX1^R320∗ specific enhancer-promoter pairs were examined for cofactor motifs. Using this approach, we successfully detected the enrichment of the RUNX motif in addition to ETS1/PU.1, which are factors known to cooperate with RUNX1 at enhancers and promoters. We also identified forkhead box (FOX) family motifs that were significantly enriched in RUNX1^R320∗ regulated E-P pairs, including those shared with RUNX1 (Figure 5E). FOX proteins are a large family of DNA-binding factors that play a variety of roles in different lineages, including enhancer regulation.^58-60 We next asked which FOX proteins were expressed in our leukemia model. We found that the FOXK subfamily had significantly higher expression in both RUNX1 wild-type and RUNX1^R320∗ cells (Figure 5F). FOXK2 but not FOXK1 showed co-occupancy at RUNX1 sites in the ENCODE data sets (supplemental Figure 5A). Furthermore, the FOXK2 and RUNX1 protein interaction network analysis revealed shared overlapping proteins (supplemental Figure 5B). Together, our analyses suggest a role for FOXK2 in RUNX1^R320∗ regulated enhancer-promoter networks.

To further explore the potential role of FOXK2 in RUNX1^R320∗ bound enhancers, we examined the well-described RUNX1-bound MYC super enhancer locus, where we detected significant E-P remodeling in GRID-seq (Figure 6A) and upregulation of MYC and MYC signaling (Figure 3F-G). RUNX1 has been reported to bind element 3 (E3) of the BENC (Blood ENhancer Cluster) super enhancer.^42,43,61,62 We theorized that RUNX1^R320∗ may dysregulate or hijack this super enhancer to affect MYC expression in conjunction with FOXK2 in the context of our model.

To build upon RUNX1 WT ENCODE H3K27ac and FOXK2 data and assess the presence of these factors in RUNX1^R320∗ cells at MYC enhancer loci, we performed Cleavage Under Targets & Release Using Nuclease (CUT&RUN),⁶³ followed by quantitative polymerase chain reaction (qPCR) in RUNX1 and RUNX1^R320∗ cells. We confirmed RUNX1 binding to E3 and detected significantly higher RUNX1^R320∗ binding at this enhancer (Figure 6B). FOXK2 binding to E3 was also significantly increased in RUNX1^R320∗ cells (Figure 6C). The SWI/SNF component BRG1 has been suggested to play an activating role in the MYC enhancer regions.⁶¹ However, we detected no significant change in BRG1 binding at E3 between RUNX1 WT and RUNX1^R320∗ cells (supplemental Figure 6A). We also observed increased RUNX1^R320∗ and FOXK2 binding to the NOTCH-bound MYC enhancer (N-Me) (supplemental Figure 6B-C). However, H3K27ac signal was not present in this region, indicating the potential requirement for additional cofactors for N-Me activation.

To further investigate the effect of FOXK2 on MYC expression, we performed shRNA-mediated FOXK2 knockdown (Figure 6D). The level of c-MYC was significantly reduced in both WT and RUNX1^R320∗ cells upon FOXK2 knockdown relative to the nontargeting shRNA control (Figure 6E; supplemental Figure 6D-E). Together, these data suggest a potential role for FOXK2 in RUNX1^R320∗ mediated enhancer-promoter networks, as well as in the upregulation of MYC and MYC oncogenic signaling via the BENC super-enhancer (Figure 6F).

In this work, we study how RUNX1 mutations outside the RHD promote abnormal hematopoiesis. We revealed that mutations in the C-terminus of RUNX1 are mainly nonsense or frameshift mutations and remain largely exempt from NMD, producing mutated proteins capable of DNA binding. We note that frameshift mutations may result in novel C-termini but focus on the effects of the retained portion of RUNX1 in this study. Modeling this class of mutation through endogenous gene editing and expression of RUNX1^R320∗, we detected a unique gene expression signature that differs from that induced by RUNX1 depletion. This suggests that truncation of the RUNX1 C-terminus does not function simply as a loss-of-function mutation. We demonstrated that this aberrant transcriptional program contributes to disease phenotypes, including megakaryocytic differentiation blocks and the disruption of hematopoietic and oncogenic pathways. Upon further investigation, we uncovered the remodeling of enhancer-promoter networks in RUNX1^R320∗ cells using GRID-seq. Analysis of altered E-P pairs revealed significant enrichment of the FOX transcription factor motif, which led us to examine FOXK2. Our results suggest a novel potential role for FOXK2 and RUNX1^R320∗in the alteration of enhancer-promoter networks, leading to dysregulated hematopoiesis.

Our work investigates the RUNX1 C-terminus, which has been shown to harbor pathogenic mutations across hematologic malignancies, yet these mechanisms remain incompletely understood. These mutations retain the DNA-binding RHD and therefore exhibit binding to the RUNX motifs. Previous in vitro studies have suggested that the C-terminus of RUNX1 contains multiple intramolecular inhibitory regions that impair DNA binding.^6,48,64 Interestingly, our ChIP-seq data showed that RUNX1 and RUNX1^R320∗ bind to DNA similarly. We hypothesized that our study reflects an endogenous context in which cofactor complexes act to closely regulate RUNX1 DNA binding. Furthermore, RUNX1 frequently interacts with other hematopoietic transcription factors to coregulate critical genes (reviewed in⁶⁵). Additionally, RUNX1 interacts with DNA in the context of chromatin looping and interacts with both the cohesin complex subunit STAG2⁴⁴ and multiple chromatin remodelers, such as PRC1^61,66 and SWI/SNF⁴³ complexes. Our data suggest that in an endogenous environment, RUNX1 DNA binding is modulated through interactions with a combination of factors, which culminates in similar RUNX1 and RUNX1^R320∗ binding on a genome-wide scale. Further studies are required to unravel the combinatorial influences behind RUNX1^R320∗ DNA-binding at the target sites.

Additionally, loss of the multifunctional carboxy-terminus of RUNX1 removes two highly conserved RUNX family domains; the nuclear matrix targeting signal and the terminal VWRPY domain. Independent of DNA binding, the nuclear matrix targeting signal has been reported to be critical for subnuclear localization and cooperation with PU.1, a critical hematopoietic transcription factor.^67,68 We hypothesized that mislocalization of RUNX1^R320∗ alters its subnuclear availability and interactions with nuclear matrix factors, resulting in unique transcriptional perturbations. Furthermore, the conserved VWRPY domain, essential for megakaryopoiesis and HSC maturation,⁶⁹ binds to the TLE1 corepressor and represses RUNX1 activity.^7,70 Lacking TLE1 binding may allow RUNX1^R320∗ to act as an activator at a subset of sites typically repressed by full-length RUNX1, a hypothesis supported by our transcriptome analysis. Taken together, we reason that the truncation of RUNX1 alters its subnuclear localization and ability to interact with various cofactors, resulting in unconventional RUNX1^R320∗ complexes and transcriptional dysregulation of hematopoietic pathways.

The truncation of RUNX1 enhancer-promoter networks are significantly distorted. In dysregulated E-P pairs, we found significant enrichment of the forkhead box (FOX) DNA-binding motif shared among FOX family members. A large family of 44 conserved transcription factors, FOX proteins act to regulate transcription through both direct DNA binding and cooperation with lineage-specific factors. Of the 14 subfamilies, the FOXK family, which consists of FOXK1 and FOXK2, was the most highly expressed in our leukemia model. Although FOXK2 is understudied in the hematopoietic system and unlike FOXK1, ENCODE data sets suggest that RUNX1 and FOXK2 DNA-binding sites frequently overlap. Previous studies have depicted a bivalent role for FOXK2, activating and repressing transcription in a context-dependent manner.^71-73 Our data suggest that loss of the C-terminus of RUNX1 may allow further cooperation between FOXK2 and RUNX1^R320∗, which acts to regulate a subset of enhancer-promoter connections, such as the BENC MYC super-enhancer.

In summary, we established that RUNX1 C-terminal variants consist mostly of nonsense and frameshift mutations, which are largely exempt from nonsense-mediated decay and lead to the production of truncated RUNX1 proteins. These proteins dysregulate hematopoietic transcriptional programs in a manner distinct from RUNX1 depletion. Upon further investigation, we showed that the loss of the domains in the C-terminus of RUNX1 results in the remodeling of enhancer-promoter networks, in which we uncovered a potential role for FOXK2 in cooperation with RUNX1^R320∗ in enhancer regulation.

The authors thank the staff at the Moores Cancer Center Flow Cytometry Core Facility for generous help with fluorescence-activated cell sorting and the staff at the University of California, San Diego Institute for Genomic Medicine for assistance with next-generation sequencing.

This work was supported by the National Institutes of Health (NIH), National Heart, Lung, and Blood Institute grant 5F31HL152652 (N.D.J.), and NIH, National Institute of Diabetes and Digestive and Kidney Diseases grant R01DK098808 (D.-E.Z.). This work was also supported by NIH training grant T32GM007240 (N.D.J.). This publication includes data generated at the UC San Diego IGM Genomics Center using an Illumina NovaSeq 6000, purchased with funding from the NIH SIG grant (#S10 OD026929). The University of California, San Diego Moores Cancer Center Core facilities are supported by NIH Cancer Center support grant CCSG P30CA23100.

Visual abstract and model was created with BioRender.com.

Contribution: N.D.J. and D.-E.Z. devised the study and designed the experimental strategies; N.D.J. performed the research, collected the data, and analyzed the results; D.-H.L. and P.B.C. prepared the GRID-seq and ChIP-seq libraries, respectively; Z.L. performed bioinformatic analyses pertaining to GRID-seq data sets; C.D. and L.X. performed microscopy-based data collection and analysis; K.-I.A. assisted in protein expression/interaction experiments; M.L. assisted in flow cytometry experiments and related statistical analysis; N.D.J. wrote the manuscript; X.-D.F. and B.R. provided resources, imparted expertise, and critically reviewed the manuscript; and D.-E.Z. oversaw the study, supervised manuscript preparation, and secured funding to support the study.

Conflict-of-interest disclosure: B.R. is a cofounder of Epigenome Technologies, Inc and has equity in Arima Genomics, Inc. The remaining authors declare no competing financial interests.

The current affiliation for Z.L. is School of Life Sciences, Southern University of Science and Technology, Shenzhen, Guangdong, China.

The current affiliation for D.-H.L. is School of Systems Biomedical Science, Soongsil University, Seoul, Republic of Korea.

The current affiliation for P.B.C. is Genome Institute of Singapore, Agency for Science, Technology and Research (A∗STAR), Singapore.

Correspondence: Dong-Er Zhang, University of California, San Diego, 9500 Gilman Dr, La Jolla, CA 92037; email: d7zhang@ucsd.edu.