MAPK Signaling Pathway was Dysregulated in Acute Myeloid Leukemia with RUNX1 Mutations

The purpose of this study was to investigate the relationship between RUNX1 mutations and MAPK signaling pathway in acute myeloid leukemia (AML). In this study, we analyzed miRNA expression with 5 mutant RUNX1 and 9 wild-type RUNX1 cases from TCGA database of AML. Six miRNAs were differently expressed with signicance, and three of them were related to overall survival. Predicted target genes of these 3 miRNAs were highly enriched in MAPK signaling pathway by functional enrichment with miRWalk3.0. Besides, genes among RUNX1 associated genes directly regulated by RUNX1 were involved in MAPK signaling pathway, too. Taken together, we demonstrate three DEmiRNAs and three genes correlated to RUNX1 were correlated with prognosis in AML, and RUNX1 modulated MAPK signaling pathway in AML. found in AML(37). study, with transcriptome data of miRNAs and mRNAs from GDC-TCGA-LAML, to distinguish targets and uncover the underlying mechanism of RUNX1, differentially expressed targets and enriched pathways of these targets were revealed. hsa-miR-873, hsa-miR-4473, and hsa-miR-551b were suggested to be up-regulated by RUNX1, and their expression levels were positively correlated to OS. MEF2C, TGFBR2, CACNB3 were correlated to RUNX1, and their expression levels were positively correlated to OS, too. And target genes of these miRNAs modulated MAPK signaling pathway.


Introduction
Acute myeloid leukemia (AML) is a clonal proliferation disease of immature myeloid precursors, caused by series of genetic alterations in hematopoietic stem cells/progenitors(HSPCs) (1). Despite remarkable progress has been made in AML treatment, 5-year survival rates of ∼10% (2). To improve treatment e ciency, it is urgently required to understand the reveal mechanism of AML. Transcription factors (TF) controlling gene expression are vital for hematopoietic stem cells (HSCs) maintenance, differentiation, apoptosis, and aging, such as RUNX1, GATA2, ETV6, and CEBPA (3)(4)(5). Bunches of evidence showed changes in these TFs were critical for leukemogenesis(6). Illustrating the dependency of these altered TFs in AML would provide therapeutic targets and improve treatment e ciency of patients.
The Runt-related transcription factor 1 (RUNX1/AML1) gene located on chromosome 21q22 is one of the most frequently altered genes in AML, including chromosomal translocation and loss of function mutations. In FAB classi cation, M0 is the subtype with the most frequently mutant RUNX1(65.2%), followed by M1(32.4%) and M2(30.2%) (7). According to WHO classi cation of myeloid neoplasms and acute leukemia, AML with RUNX1 mutations is recognized as a new subtype (or category) of AML.
Usually, RUNX1 forming heterodimer with CBF regulates target gene expression by activation or suppression(8, 9). RUNX1 is essential for hematopoietic stem cell generation, and largely dispensable in adult hematopoiesis. But heterozygous germline mutations in RUNX1 cause familial platelet disorder with a predisposition to AML (10). Different from chromosome translocation, loss of function (LOF) mutations in RUNX1 conferred stress-resistance to pre-leukemic stem cells (11). AML patients with RUNX1 mutations usually have poor treatment outcomes, including lower complete remission (CR) rate, inferior disease-free survival (DFS), and overall survival (OS) (12,13). With the advance in high-through sequencing technology, cytogenetic and molecular genetic abnormalities were identi ed and co-occurred with RUNX1 mutations in AML, such as TP53, PHF6, DNMT3A, IDH1, TET2, ASXL1, FLT3, NPM1, WT1, RAS, and MLL (14,15). But the molecular mechanism of chemo-resistance caused by RUNX1 mutations is still elusive. As a transcription factor, RUNX1 can activate or repress its target gene expression, including microRNAs and mRNAs. MicroRNAs (miRNAs), are small short non-coding RNAs, which can regulate post-transcriptional gene expression through binding to 3'-untranslated regions (3'-UTRs) of target mRNAs(16, 17). As miRNAs are critical regulators in cell proliferation, differentiation, migration, and apoptosis, they have been con rmed that they were involved in tissue development, carcinogenesis, and tumor metastasis. In normal hematopoiesis, the formation of hematopoietic stem cells (HSCs), maintenance of HSCs, and differentiation of HSCs to lineages are regulated by miRNAs. In leukemia, dysregulation of miRNAs also has been demonstrated in pathogenesis. miRNA expression pro les were correlated to prognosis, which can be used for diagnosis, treatment, and prognosis prediction (18-20).
As a master hematopoietic transcription factor, RUNX1 can control miRNA expression, as well as being targeted by miRNAs (21). In megakaryopoiesis, RUNX1 determined the expression of the miRNA-144/451 (22), and its expression was controlled by miR-142 in HSCs(23). Further, many target genes of RUNX1 were screened out, which were essential for hematopoiesis. And lots of evidence of these genes participating in leukemogenesis were accumulated. Hmga2, one of RUNX1 target genes, regulated expansion of myeloid progenitors (24). CSF1R and MPO controlled by RUNX1 were involved in myelogenesis (25, 26). But the relationship of RUNX1 and their targets in AML was not fully illustrated. In this study, we found that RUNX1 speci cally modulated MAPK signaling pathway by controlling miRNAs and mRNAs expression.

Data preparation
Use R package 'TCGAbiolinks-GDCquery', LAML-somatic mutation, miRNA expression, and gene expression data from GDC TCGA Acute Myeloid Leukemia (LAML) database were downloaded. There were 200 clinically annotated adult cases of de novo AML, including 19 patients with RUNX1 mutations.
Patients classi ed in M0/M1 according to the FAB classi cation, 'DECEASED 'statue, the mutations occurred in 'Runt homology domain' (amino acid change position<300) (27), no fusion-gene mutations (AML1-ETO or TEL-AML1), and available mRNA and miRNA expression data were used as criteria. There were only 5 cases with RUNX1 mutations according to these ltering criteria. 'LIVING' statue, normal karyotype, no fusion-gene mutations, and availability of mRNA and miRNA data were found in 9 cases with wild-type RUNX1. The clinical characteristics of AML samples were summarized in Supplemental 1, sample characteristics diagnosis after lter in AML. The GSE106291 dataset from GEO database was used as the external validation cohort.
Statistical methods R package 'DESeq2' (28) was used to perform the differential analysis. The differently expressed miRNAs with absolute fold changes over 0.5 and p-value less than 0.05 were identi ed with signi cance.
Spearman rank correlation was used to analyze correlations between RUNX1 with miRNAs or mRNAs, and a p-value less than 0.05 was signi cant. Cox regression model and log-rank test used for survival analysis by the R packages 'survival' and 'ggplot2'.
Prediction of miRNA target genes and identi cation of RUNX1 target genes MiRBase database was used to check miRNA two different matures, ('-5p' and '-3p'), and con rm previous IDs of differently expressed miRNAs. miRNA target gene was predicted by miRWalk3.0(29-31) which integrates resources of Targetscan, Mirdb, and experimental veri cation of Mirtarbase. Target region in 3'-UTR position and score more than 0.95 of miRNAs were used as selection criteria in at least one of these three databases.
Functional annotation and pathway enrichment analysis GSEA genesets reactome pathway, KEGG pathway, and gene ontology can be analyzed by miRWalk3.0. R packages 'clusterPro ler' was also used to draw the graph of GO functional annotation and KEGG pathway enrichment. Signi cantly enriched pathways were sorted by p values less than 0.05. With ENCODE RUNX1 CHIP-Seq on human K562(ENCSR414TYY)(32, 33) optimal and conservative peaks annotation was analyzed by R packages 'TxDb.Hsapiens. UCSC.hg38. knownGene' and 'ChIPseeker', and 'annotatePeak' function were used to de ne gene promoter, which included 10kb upstream to 0.5kb downstream of the transcription start site (TSS).

RUNX1 regulates miRNA expression in AML
To identify potential roles of RUNX1 in AML, the differentially expressed miRNAs were identi ed from 14 TCGA AML samples by using DESeq2 algorithms. Compared the expression in AML samples with wildtype RUNX1, in those with mutant RUNX1 28 miRNAs were differently expressed (DEmiRNAs)with signi cance, in which 16 DEmiRNAs were up-regulated and 12 DEmiRNAs were down-regulated ( Fig. 1 and Supplemental 2, signi cant differentially expression miRNAs and raw miRNAs expression array).

Signi cant Correlation of RUNX1 miRNAs with the Characteristics and Survival of Patients with AML
To explore the correlation of DEmiRNAs with RUNX1 in AML, the association of miRNAs and RUNX1 in expression was determined by Spearman algorithms. (Supplemental 3, signi cant correlation miRNAs to RUNX1, p < 0.005). We found 6 miRNAs were signi cantly correlated to RUNX1 in DEmiRNAs. hsa-miR-551b, has-miR-4473, and has-miR-873 were positive to RUNX1, hsa-miR-6877, has-miR-30c-2, and has-miR-30c-1 were negative to RUNX1. As shown in table I, by Cox regression analyses, the expression of 3 miRNAs (hsa-miR-873, hsa-miR-551b, and hsa-miR-4473) of these six miRNAs was correlated to patients' overall survival (OS, p < 0.005), and patients with low expression had better treatment outcome (Fig. 2). These ndings revealed that these three miRNAs could be involved in chemo-resistance by interfering expression of target mRNAs, and confer poor prognosis in AML.

RUNX1 regulates the MAPK signaling pathway through miRNAs
Check the previous IDs of hsa-miR-4473, hsa-miR-873 and hsa-miR-551b in miRBase database, we used hsa-miR-873 as hsa-miR-873-5p, and hsa-miR-551b as hsa-miR-551b-3p. miRNAs regulated gene expression at the post-transcriptional level by binding to 3' UTR of the target mRNAs. To illustrate the mechanism on chemo-resistance of those 3 miRNAs associated with OS, miRWalk3.0 was applied to predict target genes of hsa-miR-4473, hsa-miR-551b-3p, and hsa-miR-873-5p. By target regions within 3'UTR and more than 0.95 of score, 4633 target genes of these three miRNAs were identi ed. The KEGG pathway enrichment of these target genes was analyzed with miRWalk3.0 (Supplemental 4, predictions result of DEmiRNAs from miRWalk3.0.). As presented in Fig. 3A-B, with the highest scores, MAPK signaling pathway(hsa04010), Ras signaling pathway(hsa04014), and Rap1 signaling pathway(hsa04015) were found by these enrichment assays. The target genes of three miRNAs were most signi cantly enriched in the MAPK signaling pathway. The target genes in MAPK pathway signaling were shown in Supplemental 5, graph of miRNA-target genes mapping in the MAPK signaling pathway.
Total 95 miRNA-target genes by the red font mapping in the MAPK signaling pathway from KEGG Mapper, and these genes were involved in ERK, JNK, and p38 MAP kinase pathways (red color).
Bioinformatics target prediction and pathway analysis results indicated RUNX1 could regulate MAPK signaling pathway through miRNAs, which implied that MAPK signaling could modulate chemo-response of AML patients with RUNX1 mutations.

RUNX1 regulates the MAPK signaling pathway directly
To investigate the control of RUNX1 on mRNA expression, genes associated with RUNX1 expression in AML carrying wild-type RUNX1 were gathered from GEO dataset GSE106291, and the most signi cant associated genes(p<0.05) were used for KEGG enrichment analysis (Supplemental 6, KEGG enrichment of GSE106291 and row wild-type RUNX1 expression array and Supplemental 7, dot plot of KEGG enrichment of RUNX1 correlate genes in GSE106291). It was noted that the MAPK signaling pathway was the highest enriched pathway, too. One hundred and eighty-two RUNX1-related genes were regulating MAPK signaling pathway (hsa04010, pvalue=0.0071). With RUNX1 ChIP-seq data of human K562 cell line in ENCODE (Fig. 4A), it was found that 35 genes of the optimal and conserved peaks with annotation were enriched to MAPK signaling pathway (Fig. 4B). Meanwhile, in human CD34+ cells there were 9727 genes within whose promoters there were RUNX1 binding motifs (TGTGGT/TGTGTT) (38). Among these genes, 172 genes were participating in MAPK signaling pathway ( gure 4B). These results supported RUNX1 could transcriptionally control genes' expression directly, which were components of MAPK signaling pathway.
As it was shown above, by CHIP-seq data of K562 and human CD34 cells, many genes of MAPK signaling pathway were recognized and potentially controlled by RUNX1. In the RUNX1-related genes, 61.90% genes mediated MAPK signaling pathway. More than 60% genes, which had RUNX1 peaks within 10kb of their promoters (CD34 cell and K562 cell) were involved in MAPK signaling pathway. In comparison, there was less than 30% genes genome-wide with RUNX1 peaks in their promoters. In DemiRNAs target genes, there were 95 genes involved in MAPK signaling pathway, more than 60% of genes (59/95) of which were overlapped with RUNX1-related genes. Further, 63 of 95 target genes of DemiRNAs had RUNX1 binding peaks in their promoters. It meant that RUNX1 could regulate their 63 genes' expression indirectly and/or directly. Among these four subsets, there were 44 shared genes ( Fig. 5 and Supplemental 8, RUNX1 related genes target to MAPK signaling pathway). Log-rank test of these 44 genes found three of them were signi cantly related to survival (Fig. 5D, p<0.01), including MEF2C, TGFBR2, and CACNB3. It indicated that RUNX1 regulated MAPK signaling pathway indirectly (miRNAs) and directly (mRNA). Discussion RUNX1, as a critical transcription factor, played important roles in hematopoiesis, including maintaining the pool of hematopoietic stem cells, megakaryocyte differentiation, and lymphocyte development. With high-through technology, advances in the research of RUNX1 had been made. However, the mechanism of loss of function mutation in Runx1 required further investigation urgently. In this study, with transcriptome data of miRNAs and mRNAs from GDC-TCGA-LAML, to distinguish targets and uncover the underlying mechanism of RUNX1, differentially expressed targets and enriched pathways of these targets were revealed. hsa-miR-873, hsa-miR-4473, and hsa-miR-551b were suggested to be up-regulated by RUNX1, and their expression levels were positively correlated to OS. MEF2C, TGFBR2, CACNB3 were correlated to RUNX1, and their expression levels were positively correlated to OS, too. And target genes of these miRNAs modulated MAPK signaling pathway.
MAPKs is highly conserved protein kinases, including nineteen species MAPKKK (also named MAP3K or MEKK) were found, seven kinds of MAPKK (MAP2K) and fourteen kinds of MAPK. They transmit extracellular signals to intracellular to regulate cell growth, differentiation, apoptosis, migration, and so on. There are four distinct MAPK signaling cascades, including extracellular signal-related kinases (ERK1/2 and ERK5) Jun amino-terminal kinases(JNK1/2/3) and p38 kinase (39). JNK and p38 MAPK share roles in in ammation, apoptosis, and growth. The kinases of this pathway were dysregulated in various diseases, including cancers so that they could be used as biomarkers or therapeutic targets. MEK inhibitor (PD0325901) was tested in a phase II clinical trial for advanced non-small cell line lung cancer (NSCLC). Moreover, LGX818 is a BRAF inhibitor, which also has shown potential bene t for NSCLC patients. Inhibitor for p35 MAPK improved behavior and synaptic function of mouse with Alzheimer's disease. Inhibition of JNK pathway provided a potent regulator for brain in ammation in neurological disorders (40,41).
In AML, MAPK pathway was activated, and ERK phosphorylation was associated with cell growth (42), Flt3 ligand stimulation caused rapid activation of MAPK (43). It has been reported that combinatorial blockade of the MAPK pathways promoted cell death and suppresses proliferation in the majority of primary AML cells (44). Further, activated MAPK signaling contributed to primary therapeutic resistance to an IDH2 inhibitor in AML (45). pan-RAF inhibition, had been proved that it was a promising treatment strategy for AML(46). And in preclinical models, pexmetinib as a p38 MAPK inhibitor showed its potential for MDS/AML treatment (47). In AML, MAPK signaling pathway could be regulated by miR-181a, miR-497, and miR-143. And it has been identi ed that RUNX1 could affect the activity of MAPK signaling pathway miRNAs and mRNAs(48-50).
In this study, we found Six miRNAs were differently expressed with signi cance, and three of them were related to overall survival. Predicted target genes of these 3 miRNAs were highly enriched in MAPK signaling pathway by functional enrichment with miRWalk3.0. Besides, genes among RUNX1 associated genes directly regulated by RUNX1 were involved in MAPK signaling pathway, too. Taken together, RUNX1 modulated MAPK signaling pathway in AML. These results could explain the poor treatment outcome with mutant RUNX1, and enhance our understanding of the TF-miRNA regulatory associated with cell proliferation, apoptosis, and survival with MAPK signaling pathway.

Conclusions
Here, in GDC-TCGA-AML, we found new targets, including miRNAs and mRNAs, enriched in MAPK signaling pathway, and presumably proved RUNX1 regulated MAPK signaling pathway. It was concluded that this regulatory network was of great signi cance to the phenotype of AML with RUNX1 LOF mutations, and this also provided a potential target for clinical treatment in AML. At the same time, there were some shortcomings, the initial samples were relatively small and these research conclusions need to be veri ed by clinical specimens and animal models.      Evidences of RUNX1 related genes target to MAPK signaling pathway (A) Distribution of peak's closest genetic distance, in other words, distribution transcription factor binding loci relative to TSS. gene promoter was 10kb upstream to 0.5kb downstream of the TSS. (B) K562 peaks annotation genes and CD34 RUNX1 binding coding genes match to MAPK signaling pathway genes.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.