Elucidation of Novel Therapeutic Targets for Acute Myeloid Leukemias with RUNX1-RUNX1T1 Fusion

The RUNX1-RUNX1T1 fusion is a frequent chromosomal alteration in acute myeloid leukemias (AMLs). Although RUNX1-RUNX1T1 fusion protein has pivotal roles in the development of AMLs with the fusion, RUNX1-RUNX1T1, fusion protein is difficult to target, as it lacks kinase activities. Here, we used bioinformatic tools to elucidate targetable signaling pathways in AMLs with RUNX1-RUNX1T1 fusion. After analysis of 93 AML cases from The Cancer Genome Atlas (TCGA) database, we found expression of 293 genes that correlated to the expression of the RUNX1-RUNX1T1 fusion gene. Based on these 293 genes, the cyclooxygenase (COX), vascular endothelial growth factor receptor (VEGFR), platelet-derived growth factor receptor (PDGFR), and fibroblast growth factor receptor (FGFR) pathways were predicted to be specifically activated in AMLs with RUNX1-RUNX1T1 fusion. Moreover, the in vitro proliferation of AML cells with RUNX1-RUNX1T1 fusion decreased significantly more than that of AML cells without the fusion, when the pathways were inhibited pharmacologically. The results indicate that novel targetable signaling pathways could be identified by the analysis of the gene expression features of AMLs with non-targetable genetic alterations. The elucidation of specific molecular targets for AMLs that have a specific genetic alteration would promote personalized treatment of AMLs and improve clinical outcomes.


Introduction
Acute myeloid leukemia (AML) with RUNX1-RUNX1T1 fusion consists of up to 5% of all AML or 10% of AML M2 subtype in the French-American-British classification [1]. The RUNX1-RUNX1T1 fusion transcript results in aberrant products that create nuclear transcriptional co-repressor complexes and suppress the expression of RUNX1 target genes [2][3][4]. Since RUNX1T1 silencing with short hairpin RNAs has shown in vitro therapeutic effects on AML cells with the RUNX1-RUNX1T1 fusion [5], RUNX1T1 overexpression has been suggested as a leukemogenesis factor. Alteration in the specificity of RUNX1 to its target sequences by the fusion might also contribute to leukemogenesis [4]. Moreover, the RUNX1-RUNX1T1 fusion is associated with the differential clinical prognosis of AML patients [6].
Specific genomic alterations not only improve our understanding of leukemogenesis mechanisms, but also are potential therapeutic targets for targeted chemotherapies. A good example would be the inhibition of enzymatic kinase activity of BCR-ABL fusion protein in chronic myeloid leukemia (CML) with imatinib (Gleevec ® ) [7]. However, as it is relatively difficult to prevent the RUNX1-RUNX1T1 fusion protein from mediating non-enzymatic transcriptional repression in vivo, inhibition of signaling pathways modulated in AML cells with the RUNX1-RUNX1T1 fusion cab be considered.
Herein, we utilized the RNA-sequencing (RNA-seq) data of The Cancer Genome Atlas (TCGA) [8], and compared the gene expression profile of seven RUNX1-RUNX1T1 fusion-positive AMLs with those of 86 AMLs that have normal karyotype, in order to systemically identify important cancer signalings and alternative druggable targets of the AMLs with the RUNX1-RUNX1T1 fusion. The analysis revealed specific activation of signaling pathways in the AMLs with the RUNX1-RUNX1T1 fusion, including COX, VEGF, PDGF, and FGFR1, which were pharmacologically inhibited in vitro to show the specific inhibition of the proliferation of AML cells with the RUNX1-RUNX1T1 fusion. Through in silico and in vitro validation, our study provides evidence for repurposing molecular targeting agents approved in different types of cancers or other diseases.
Seven RUNX1-RUNX1T1 fusion-positive and 86 fusion-negative (without detectable genomic structural abnormalities) AMLs were selected from 200 patients available in the Broad DAC Firehose (Figure 1). Group-wise comparison of the clinical and pathological characteristics (Table 1) showed a difference in the mean age at diagnosis; the RUNX1-RUNX1T1 fusion-positive patients were significantly younger than the fusion-negative patients (p < 0.01). In peripheral blood analysis (Table 1), the RUNX1-RUNX1T1 fusion-positive patients had significantly lower blast counts (p < 0.01), hemoglobin measurements (p < 0.01), and platelet counts (p < 0.01) at diagnosis. Significantly lower bone marrow cellularity was also observed in the RUNX1-RUNX1T1 fusion-positive AMLs (p < 0.01, Table 1). Although it was not exclusive, the CD19 expression ratio of cancer cells was significantly higher in the AMLs with RUNX1-RUNX1T1 fusion.
These clinical and pathological features of the RUNX1-RUNX1T1 fusion-positive AMLs were well-matched with those recorded in previous reports [9][10][11][12], indicating the validity of selecting the AMLs with RUNX1-RUNX1T1 fusion in the TCGA database.

Genes and Pathways Specifically Altered in AMLs with RUNX1-RUNX1T1 Fusion
Using the transcriptome data of the 7 RUNX1-RUNX1T1 fusion-positive and 86 fusion-negative AMLs in the TCGA database, 293 genes ( Figure 1 and Table S1) were identified to have expressions significantly correlated with the RUNX1-RUNX1T1 fusion transcript level. Among the 293 genes, 42 genes were cancer-related based on the predefined cancer gene list (2027 genes, http://www. bushmanlab.org/links/genelists) ( Figure S1). CAV1, POU4F1, and ROBO1 were most significantly overexpressed in correlation with the RUNX1T1 expression, which is concordant with the previous report [13]. Interestingly, POU4F1 is known to be up-regulated by RUNX1-RUNX1T1 fusion, and the fusion with POU4F1 has a synergy in driving B-lymphoid gene expression in t (8:21) AML [14,15]. FGFR1, VEGFA, NBL1, TET1, IKZF2, CCND1, and MPL also showed positive correlations (Pearson R > 0.4) with RUNX1T1 in RNA levels. Although not previously shown to be directly related with cancers, several surface markers, such as CD19 and CD34, showed significant positive correlation (Pearson R > 0.5) with RUNX1T1 in RNA levels, whereas CD33 showed negative correlation (Pearson R = −0.36) ( Figure S1).
Subsequently, we applied ConsensusPathDB (CPDB) to the 293 genes to elucidate the specific functional signaling pathways that are altered in the AMLs with RUNX1-RUNX1T1 fusion ( Figure 1). The analysis was performed with the 4,011 curated gene sets provided by CPDB. In the over-representation analysis (ORA) of CPDB, 24 pathways showed statistical significance (q-values < 0.1) (Table 2 and Figure 2). Gene Set Enrichment Analysis (GSEA) was utilized to cross-validate the results from ORA, using CPDB. GSEA showed good normalized enrichment scores (NESs) of over 1.5 in the VEGFR1 specific signals, the PDGFRA signaling pathway, the celecoxib pathway, and the phospholipase C mediated cascade: FGFR1 ( Figure S2). The GSEA results indicated that the affected pathways were cross-validated with two independent in silico pathway analysis methods. Subsequently, we applied ConsensusPathDB (CPDB) to the 293 genes to elucidate the specific functional signaling pathways that are altered in the AMLs with RUNX1-RUNX1T1 fusion ( Figure  1). The analysis was performed with the 4,011 curated gene sets provided by CPDB. In the over-representation analysis (ORA) of CPDB, 24 pathways showed statistical significance (q-values < 0.1) (Table 2 and Figure 2). Gene Set Enrichment Analysis (GSEA) was utilized to cross-validate the results from ORA, using CPDB. GSEA showed good normalized enrichment scores (NESs) of over 1.5 in the VEGFR1 specific signals, the PDGFRA signaling pathway, the celecoxib pathway, and the phospholipase C mediated cascade: FGFR1 ( Figure S2). The GSEA results indicated that the affected pathways were cross-validated with two independent in silico pathway analysis methods. Heat Map of genes with altered expression correlating to RUNX1T1 expression. Genes relevant to the COX, VEGF-vascular endothelial growth factor receptor (VEGFR), PDGF, and FGFR1 related pathways were found to be altered in AMLs with RUNX1-RUNX1T1 fusion, among a total of 293 genes with over-or under-expression, compared to normal karyotype AML patients. Pathways with q value < 0.1 were selected and merged based on the pathway ontology.  Several genes of the 293 genes were involved in multiple pathways in the pathway analyses. Especially, VEGFA, CAV1, FGFR1, PLCG1, and PRKCD were involved in at least three pathways (Figure 3), implying their functional roles in the AMLs with RUNX1-RUNX1T1 fusion. Overall, 14 genes participated in multiple pathways ( Figure 3 and Table S2). Several genes of the 293 genes were involved in multiple pathways in the pathway analyses. Especially, VEGFA, CAV1, FGFR1, PLCG1, and PRKCD were involved in at least three pathways (Figure 3), implying their functional roles in the AMLs with RUNX1-RUNX1T1 fusion. Overall, 14 genes participated in multiple pathways ( Figure 3 and Table S2).

Pathway Prioritization by Targeting Possibility
From the CPDB ORA results, we identified COX-2, VEGF, PDGF, and FGFR1 pathways that can be pharmacologically targeted ( Table 2). The concurrent expression of RUNX1T1 and a gene that is associated with its respective pathway are outlined in Figure 2. Overall, 34 genes were culled from the 293 genes based on their relationship with the COX-2, VEGF, PDGF, and FGFR1 pathways.
Among the 34 genes, 19 were associated with the COX-2 pathway ( Figure 2). As previously suggested, RUNX1-RUNX1T1 fusion mediates leukemogenesis through the COX pathway [16]; this result thus indicates the validity of the analysis of this study.
The VEGF/VEGFR pathway has interesting features, since 4 pathways retrieved from different database sources showed q-values < 0.1 with p-values < 0.004 (Table 2). This suggests the potential for drug repositioning of the VEGF/VEGFR targeting agents from anti-angiogenic treatment in order to direct AML cell-targeting therapy, which is in concordance with previous experimental studies [17]. The PDGF and FGFR1 pathways were also interesting, since they are associated with eosinophilia [18,19]. Eosinophilia has been known to be a characteristic of AMLs with RUNX1-RUNX1T1 fusion [2]. In addition, the PDGF and FGFR1 pathways are associated with the overexpression of CD19, which is the immunophenotypic feature of RUNX1-RUNX1T1 positive AMLs [2]. Considering these pathways can be pharmacologically inhibited by multiple specific targeted agents, the results suggest novel molecular targeted therapies for AMLs with RUNX1-RUNX1T1 fusion.

Pathway Prioritization by Targeting Possibility
From the CPDB ORA results, we identified COX-2, VEGF, PDGF, and FGFR1 pathways that can be pharmacologically targeted ( Table 2). The concurrent expression of RUNX1T1 and a gene that is associated with its respective pathway are outlined in Figure 2. Overall, 34 genes were culled from the 293 genes based on their relationship with the COX-2, VEGF, PDGF, and FGFR1 pathways.
Among the 34 genes, 19 were associated with the COX-2 pathway ( Figure 2). As previously suggested, RUNX1-RUNX1T1 fusion mediates leukemogenesis through the COX pathway [16]; this result thus indicates the validity of the analysis of this study.
The VEGF/VEGFR pathway has interesting features, since 4 pathways retrieved from different database sources showed q-values < 0.1 with p-values < 0.004 (Table 2). This suggests the potential for drug repositioning of the VEGF/VEGFR targeting agents from anti-angiogenic treatment in order to direct AML cell-targeting therapy, which is in concordance with previous experimental studies [17]. The PDGF and FGFR1 pathways were also interesting, since they are associated with eosinophilia [18,19]. Eosinophilia has been known to be a characteristic of AMLs with RUNX1-RUNX1T1 fusion [2]. In addition, the PDGF and FGFR1 pathways are associated with the overexpression of CD19, which is the immunophenotypic feature of RUNX1-RUNX1T1 positive AMLs [2]. Considering these pathways can be pharmacologically inhibited by multiple specific targeted agents, the results suggest novel molecular targeted therapies for AMLs with RUNX1-RUNX1T1 fusion.

Functional Validation of Signaling Pathways
Receptor tyrosine kinase (RTK) array was conducted to confirm the phosphorylation of the RTKs in the VEGF, PDGF, or FGFR1 pathways. The phosphorylation of 71 different RTKs was screened against SKNO-1 and Kasumi-1 human AML cells that have RUNX1-RUNX1T1 fusion [20,21]. THP-1 AML cells were also analyzed as a negative control, as they do not harbor the RUNX1-RUNX1T1 fusion [22]. We further analyzed the phosphorylated status of EGFRs, FGFRs, PDGFRs, and VEGFRs. EGFRs were utilized as negative controls that were not predicted in the transcriptome analyses, to be activated in AMLs with RUNX1-RUNX1T1 fusion.
As expected, the phosphorylation levels of FGFRs, PDGFRs, and VEGFRs were 3~5-fold higher than those of the EGFRs in both SKNO-1 ( Figure 4A) and Kasumi-1 ( Figure 4B). In particular, FGFR, PDGR, and VEGFR were highly phosphorylated in SKNO-1 ( Figure 4A), whereas PDGFR and VEGFR were highly phosphorylated in Kasumi-1 ( Figure 4B). In contrast, the phosphorylation level of EGFRs, FGFRs, and PDGFRs were almost the same in THP-1 ( Figure 4C). Moreover, the phosphorylation level of the EGFRs of THP-1 were higher than those of SKNO-1 and Kasumi-1, while the phosphorylation levels of the FGFRs, PDGFRs, and VEGFRs of THP-1 were lower than those of SKNO-1 and Kasumi-1 ( Figure 4D). These results indicate that the FGFR, PDGFR, and VEGFR signaling pathways may be specifically activated in AMLs with RUNX1-RUNX1T1 fusion. Receptor tyrosine kinase (RTK) array was conducted to confirm the phosphorylation of the RTKs in the VEGF, PDGF, or FGFR1 pathways. The phosphorylation of 71 different RTKs was screened against SKNO-1 and Kasumi-1 human AML cells that have RUNX1-RUNX1T1 fusion [20,21]. THP-1 AML cells were also analyzed as a negative control, as they do not harbor the RUNX1-RUNX1T1 fusion [22]. We further analyzed the phosphorylated status of EGFRs, FGFRs, PDGFRs, and VEGFRs. EGFRs were utilized as negative controls that were not predicted in the transcriptome analyses, to be activated in AMLs with RUNX1-RUNX1T1 fusion.
( Figure 5F), while Kasumi-1 cells rely on the functions of VEGFRs ( Figure 5E) and PDGFRs ( Figure  5F). On the contrary, the IC50s of the EGFR inhibitor in the SKNO-1 (Figure 5A), Kasumi-1 ( Figure  5B), and THP-1 ( Figure 5C) cells did not significantly differ ( Figure 5G). Moreover, the VEGR inhibitor had a lower IC50 than the FGFR and PDGFR inhibitors in SKNO-1 ( Figure S3A) and the PDGFR and VEGFR inhibitors had lower IC50s than the FGFR inhibitor in Katumi-1 ( Figure S3B); this is well-matched with the phosphorylation status of RTKs (Figure 4). In contrast, the effects of the FGFR, PDGFR, and VEGFR inhibitors in THP-1 did not differ compared with the EGFR inhibitors ( Figure S3C) [22].

Discussion
In this study, we investigated cancer-related targetable pathways that are specifically regulated by the RUNX1-RUNX1T1 fusion in AMLs. The unique aspect of the pathway analysis of this study is that the gene sets for ORA were selected based on the expression correlation with the reference gene,

Discussion
In this study, we investigated cancer-related targetable pathways that are specifically regulated by the RUNX1-RUNX1T1 fusion in AMLs. The unique aspect of the pathway analysis of this study is that the gene sets for ORA were selected based on the expression correlation with the reference gene, RUNX1T1. The RUNX1T1 expression level in each RNA-seq sample might be skewed by several factors, such as its tumor purity or subclonal heterogeneity. Expression of the genes of the RUNX1T1 downstream pathways could be affected in a similar pattern. Therefore, the gene set selection based on the expressional correlation with RUNX1T1 has strong potential, since it offers gene selection while considering the extent of the altered gene expression in each sample. Some cancer-related pathways identified in this study were well-matched with those of previous studies of RUNX1-RUNX1T1 fusion-positive AMLs. For example, we identified that COX-related pathways might be activated in AMLs with the RUNX1-RUNX1T1 fusion. Accordingly, it was reported that COX-2 inhibition could reduce the proliferation of RUNX1-RUNX1T1 fusion-positive AML cells in vitro, and in vivo [16]. PDGF-related pathways identified in this study can explain the reason why this leukemia subtype is related to eosinophilia [23] and elevated CD19 expression [24]. These concurrences and consistencies with the well-known features of AMLs with the RUNC1-RUNX1T1 rearrangement might support the validity of the analysis of this study.
Besides the COX-and PDGF-related pathways, the VEGF-and FGFR1-related pathways were also predicted to be specifically activated in AMLs with the RUNX1-RUNX1T1 fusion in this study. The signaling pathways provide novel opportunities for tailored treatments of RUNX1-RUNX1T1 fusion positive AMLs [25][26][27].
In parallel with the in silico analysis of this study, PDGFRs, FGFRs, and VEGFRs were highly phosphorylated in two types of AML cell lines with the RUNX1-RUNX1T1 fusion (SKNO-1 and Kastumi-1). In contrast, the phosphorylation levels were much lower in the THP1 AML cells that have no fusion. The activation status indicated that AMLs with the RUNX1-RUNX1T1 fusion could be more sensitive to PDGFRs, FGFRs, and VEGFRs targeting agents, which were experimentally validated in this study. Although these results are in vitro validations, and not conducted studies done on patient samples, it does shed light on new possibilities for drug reposition for AML with RUNX1-RUNX1T1.
Combinational chemotherapy is a future direction of effective treatments for various cancer types [28]. The US Food and Drug Administration (FDA) has approved several clinical trials on combinational chemotherapy for AML patients such as cytarabine/daunorubicin and PF-04449913(Glasdegib) [29]. Our results indicate that RTK pathway-targeting agents have adverse effects on the proliferation potency of RUNX1-RUNX1T1 fusion cell lines. It showed the possibility that use of RTK targeting agents in combination with the conventional chemotherapeutic drugs for AMLs, like Ara-C, might have better anti-tumor effects. However, it needs more preclinical study to validate which combination strategy is the best, as there are various possible combinations.
Currently, some molecular alterations are used for AML diagnosis, such as FLT3, NPM1, and CEBPA, several of which could be targeted since their protein products have enzymatic activities. In 2017, FLT3 inhibitors were approved to treat AML patients with FLT3-ITD mutations [30]. Despite these efforts, improvement of clinical outcome in AMLs is still limited to a small population of patients [31]. This is partially because many AML-specific genetic alterations cannot be targeted pharmacologically. Therefore, the identification of targetable molecules that are closely associated with AML-specific genetic alterations is important for advancement in the precision medicine of AMLs. The analysis strategy presented in this study might be helpful in the identification of treatment targets and in extending the indication of approved target therapy agents.
In this study, we elucidated gene sets, the expressions of which are specifically changed in AMLs with the RUNX1-RUNX1T1 fusion. Fidelity of identification was upgraded by hiring the expression correlation of candidate genes with the reference gene, RUNX1T1 in the analysis processes. Based on the gene sets, activated signaling pathways in AMLs with the RUNX1-RUNX1T1 fusion were predicted and their functional importance in treatment was validated experimentally. This analysis process could also be adopted in the discovery of molecular targets of AMLs that have no targetable genetic alteration. Further, we expect that combination therapy using the conventional drug (Ara-c) with RTK targetable drugs shows better anti-tumor effects. Despite these encouraging results, these suggestions are required to be further investigated with clinical samples to be fully validated. With future studies possibly validating these results, it can also provoke clinical trials for personalized treatments of AMLs.

Data Acquisition and Inclusion Criteria
Gene level 3(RSEM) mRNA expression for The Cancer Genome Atlas (TCGA) Acute Myeloid Leukemia (AML) study was downloaded from the Broad GDAC Firehose. Clinical information, obtained from the same site, includes chromosome analysis, fluorescence in situ hybridization (FISH) results, immunophenotyping results and other clinical features. Using the results of chromosome analysis and FISH studies, we identified seven samples with RUNX1-RUNX1T1 fusion, which were cross-checked with elevated RUNX1T1 expression level. For the controls, 100 samples with normal karyotype were selected. After filtering 14 samples without RNA expression data, 86 samples were used as controls.

Over-Representation Pathway Analysis and Gene Set Enrichment Analysis
RNA expression data from each TCGA AML sample were merged into a two-dimensional matrix. Each column of the matrix represents the patient and each row represents each gene name. The RNA expression values were described as normalized read counts. To obtain the genes that are correlated in RNA expression with the reference gene (RUNX1T1) in 93 samples, a Pearson correlation test and Spearman correlation test were performed for each gene. Genes with R > 0.3 in both Pearson and Spearman correlation tests were selected, and 294 genes were finally selected as genes closely correlated with RUNX1T1 in RNA expression. With the 294 genes, over-representation analysis (ORA) was performed using ConsensusPathDB (CPDB, http://consensuspathdb.org/) according to current protocols [32]. In this analysis, 4011 pathways, curated from multiple sources including INOH [33], NetPath [34], Reactome [35], HumanCyc [36], KEGG [37], Wikipathways [38], SMPDB [39], PharmGKB [40], EHMN [41] and, Signalink [42] were used.
Gene Set Enrichment Analysis (GSEA) [43] was performed to crosscheck the ORA result. The GSEA analysis software (version 2.2.3) was downloaded from the Broad Institute website (http://www.broadinstitute.org/gsea/index.jsp). The curated gene set provided by CPDB was also downloaded and modified for GSEA analysis.

Network Analysis for Potential Actionable Drugs and Target Genes
To analyze and visualize the relationship between therapeutic agents and their target genes, the expressions of which were specifically altered in the AMLs with RUNX1-RUNX1T1 fusion, Cytoscape version 3.5.1 was used [44]. Knowledge-based databases including CIViC [45] and CancerSCAN [46] were used for the analysis.

Statistical Analysis and Visualization
To select genes that correlate well with RUNX1T1 in RNA expression, the Pearson correlation test and Spearman correlation test were used. To visualize the RNA expression heat map with the related pathways, the R package of the ComplexHeatmap was used [47]. All statistical analysis and visualization were performed based on the open software R version 3.4.3 [48]. We applied a p-value of < 0.05 for statistical relevance, and a q-value of < 0.1 for false detection rate (FDR) control.

Receptor Tyrosine Kinase (RTK) Phosphorylation Array
Cells were cultured with their growth media in 75T flasks. Pellets were prepared upon~80% confluence. A lysis buffer containing a protease and phosphatase inhibitor cocktail (AAH-PRTK-1, RayBiotech, GA) was used to separate the protein from the pellets. Protein concentration was determined by the bicinchoninic acid (BCA) assay (Pierce, Rockford, IL, USA). 700µg protein was reacted with a receptor tyrosine kinase (RTK) phosphorylation array (AAH-PRTK-1) according to the manufacturer's protocol. Spots on the array were analyzed by ImageJ software(National Institutes of Health, Maryland, USA).

Drug Sensitivity Test
The cells were seeded in 384-well plates at a density of 500 cells per well in triplicate for each treatment. The drugs panel consisted of 61 anti-cancer agents (Selleckchem, TX, USA) targeting oncogenic signals. Two hours after plating, the cells were treated with the drugs in a seven-point serial dilution for 6 days at 37 • C in a 5% CO 2 humidified incubator. Cell viability was analyzed using an ATP-monitoring system based on firefly luciferase (ATPLit 1step, PerkinElmer, MA, USA). Dimethyl sulfoxide (DMSO) was included as a negative control in each plate. Controls were used for the calculation of relative cell viability for each plate, and normalization was performed on a per-plate basis. Dose-response curve (DRC) fitting was performed using GraphPad Prism 5 (GraphPad) and evaluated by measuring the half maximal inhibitory concentration (IC 50 ) of the DRC. After normalization, the best-fit lines were determined, and the IC 50 value of each curve was calculated using GraphPad Prism, ignoring the regions defined by fewer than two peaks.

Conflicts of Interest:
The authors declare no conflict of interest.