MicroRNA-363-3p promote the development of acute myeloid leukemia with RUNX1 mutation by targeting SPRYD4 and FNDC3B.

BACKGROUND: Runt-related transcription factor 1 (RUNX1) is one of the most frequently mutated genes in most of hematological malignancies, especially in acute myeloid leukemia. In the present study, we aimed to identify the key genes and microRNAs based on acute myeloid leukemia with RUNX1 mutation. The newly finding targeted genes and microRNA associated with RUNX1 may benefit to the clinical treatment in acute myeloid leukemia. MATERIAL/
METHODS: The gene and miRNA expression data sets relating to RUNX1 mutation and wild-type adult acute myeloid leukemia (AML) patients were downloaded from The Cancer Genome Atlas database. Differentially expressed miRNAs and differentially expressed genes (DEGs) were identified by edgeR of R platform. Gene ontology and the Kyoto Encyclopedia of Genes and Genomes enrichment analyses were performed by Metascape and Gene set enrichment analysis. The protein-protein interaction network and miRNA-mRNA regulatory network were performed by Search Tool for the Retrieval of Interacting Genes database and Cytoscape software.
RESULTS: A total of 27 differentially expressed miRNAs (25 upregulated and 2 downregulated) and 561 DEGs (429 upregulated and 132 downregulated) were identified. Five miRNAs (miR-151b, miR-151a-5p, let-7a-2-3p, miR-363-3p, miR-20b-5p) had prognostic significance in AML. The gene ontology analysis showed that upregulated DEGs suggested significant enrichment in MHC class II protein complex, extracellular structure organization, blood vessel development, cell morphogenesis involved in differentiation, embryonic morphogenesis, regulation of cell adhesion, and so on. Similarly, the downregulated DEGs were mainly enriched in secretory granule lumen, extracellular structure organization. In the gene set enrichment analysis of Kyoto Encyclopedia of Genes and Genomes pathways, the RUNX1 mutation was associated with adherent junction, WNT signaling pathway, JAK-STAT signaling pathway, pathways in cancer, cell adhesion molecules CAMs, MAPK signaling pathway. Eleven genes (PPBP, APP, CCR5, HLA-DRB1, GNAI1, APLNR, P2RY14, C3AR1, HTR1F, CXCL12, GNG11) were simultaneously identified by hub gene analysis and module analysis. MicroRNA-363-3p may promote the development of RUNX1 mutation AML, targeting SPRYD4 and FNDC3B. In addition, the RUNX1 mutation rates in patient were obviously correlated with age, white blood cell, FAB type, risk(cyto), and risk(molecular) (P < .05).
CONCLUSION: Our findings have indicated that multiple genes and microRNAs may play a crucial role in RUNX1 mutation AML. MicroRNA-363-3p may promote the development of RUNX1 mutation AML by targeting SPRYD4 and FNDC3B.

Introduction

Acute myeloid leukemia (AML), the most common acute leukemia in adults,[ is a malignant clonal disorder that inhibits differentiation of cells and induces proliferation or accumulation of blasts, instead of producing healthy hematopoietic cells.[ While acute myeloid leukemia mostly invades elderly patients with poor survival prognosis.[ As the leading cause of acute leukemia in the adult population, the number of deaths is more than other type of leukemia, such as acute lymphocytic leukemia, chronic lymphocytic leukemia, and chronic myeloid leukemia.[ Abnormal gene expression and epigenetic changes in the genome promote the development of AML.[ A large number of gene mutations related to AML patients have been discovered, such as DNMT3A, ASXL1, TET2, IDH1, IDH2, and FLT3, studied for improving the efficiency of early diagnosis as well as the effect of initial treatment.

Regulating the expression of genes[ and various cellular activities, microRNAs (miRNAs) is regarded as the short noncoding RNA molecules, such as cell growth, development, differentiation, as well as apoptosis.[ Most of genes can be influenced by just a single miRNA that involved in the functional pathway by targeting relational mRNAs.[

Including Runt-related transcription factor 1 (RUNX1), RUNX2, RUNX3, the family of RUNX transcription factors play a major role in the regulation of cell identities and functions. RUNX1, spanned ∼261 kb on the long arm of chromosome 21,[ encodes a sequence-specific transcription factor, which is one of the most common mutated genes in most of hematological malignancies and is essential for the formation of the differentiation of lymphoid, myeloid, and megakaryocytic cells.[ Specially, located at the signalized point of t(8,21), RUNX1 is regarded as the gene that is associated with acute myeloid leukemia.[RUNX1 is involved in hematopoietic differentiation, ribosome biogenesis, cell cycle regulation, and is one of the p53 and transforming growth factor β signaling pathways.[ Fangxiao Zhu et al[ have investigated the key genes associated with RUNX1 mutations in AML. But the differentially expressed miRNAs and relevant pathways have not been analyzed. Therefore, in order to better understand the primary biological processes associated with RUNX1 mutation in adult AML, we aim to identify the key miRNAs and pathway through bioinformatics analysis.

Materials and methods

Data collection

The gene and miRNA expression data sets and clinical information of patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://gdc-portal.nci.nih.gov/).[

Identification of differentially expressed miRNAs (DE-miRNAs) and differentially expressed genes (DEGs)

Limma, an R package for examining gene expression microarray data, was utilized to screen differential expression of miRNAs and differential expression of genes between RUNX1 mutation and wild-type AML patients according to the user's guide.[ The DE-miRNAs were identified by P value <.05 and |log2fold change (FC)|≥1, while P value <.05 and |log2FC|≥1.5 were set as the threshold values for DEG identification.

Functional enrichment analysis of DEGs

Metascape (http://metascape.org) was utilized to carry out gene ontology (GO) term analysis, providing an extensive gene list annotation and analysis resource for experimental biologists.[ While the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was analyzed by the gene set enrichment analysis (GSEA) (https://www.gsea-msigdb.org/), a method that was used to evaluate whether a set of genes showed statistically significant and consistent differences between 2 biological states.[ Significant gene sets with an FDR < 0.25 and a P value < .05 were identified.

Prediction of target genes of DE-miRNAs

The target genes of DE-miRNAs were identified by Target Scan,[ miRDB,[ miRPathDB,[ and miRWalk,[ all of which were the target prediction databases. Simultaneously, the total target genes, indicated by these 4 databases, were corresponded to DGEs (P < .05, not limited fold change) for further research, while the target genes of upregulated miRNA matched with downregulated genes and target gens of downregulated miRNA matched with upregulated genes.

Protein–protein interaction (PPI) network, miRNA-mRNA regulatory network

To further explore the relationships between DEGs at the protein level, DEGs were uploaded to Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/).[ Subsequently, we used Cytoscape software to construct the PPI network of DEGs and the miRNA-mRNA regulatory network, in which the hub genes and screening modules were also identified by Molecular Complex Detection plugin and cytoHubba plugin.[ Moreover, all the parameters of plug-in were kept as the default values.

Statistical analysis

Evaluating the prognostic value of DE-miRNAs, the Kaplan–Meier method was utilized to construct survival curves. Among them, according to the Cox proportional hazards regression model, the hazard ratio (HR) with 95% confidence intervals (CIs) and log rank P value were calculated. χ analysis and t test were used to evaluate the clinical characteristic and expression between the RUNX1 mutations and RUNX1 wild type in adult AML. All the statistical analyses were conducted with SPSS version 20.0 and GraphPad Prism version 8.0. A value of P<.05 was considered statistically significant.

Ethical approval

All the data in this study were obtained from open, public databases; therefore, ethical approval was not necessary.

Results

Procedure of bioinformatics analysis

First, the DE-miRNAs and DEGs of AML with RUNX1 mutations were identified, especially, the valuable DE-miRNAs were found. Second, GO and GSEA enrichment analyses of DEGs were accomplished. Subsequently, we built PPI network of DEGs and proposed hub genes and module analysis. In addition, the valuable DE-miRNAs were evaluated by the survival analysis, and the miRNA-mRNA regulatory network was constituted by those target genes. Ultimately, miR-363-3p and its target genes were selected for further analyses (Fig. 1).

Figure 1

The flow chart of bioinformatics analysis. AML = acute myeloid leukemia, DEGs = differentially expressed genes, DE-miRNAs = differentially expressed miRNAs, GO = gene ontology, GSEA = gene set enrichment analysis, PPI = protein–protein interaction, TCGA = The Cancer Genome Atlas.

Identification of differentially expressed genes (DEGs) and differentially expressed miRNAs (DE-miRNAs)

The study has shown that RUNX1 mutations were remarkably related to poor clinical outcomes as well as increased risk of death in AML patients (log-rank P = .012, HR = 2.145, 95% CI = 1.156–3.982).[ Furthermore, the expression of RUNX1 between RUNX1 mutation and wild-type patients was dissimilar (Fig. 2F). In a word, the crucial mRNA and miRNA played significant roles in RUNX1 mutation.

Figure 2

The expression level of RUNX1 (mutation and wild-type) and prognostic significance of 5 DE-miRNAs. Kaplan–Meier survival curves of overall survival (OS) in AML were performed using TCGA dataset. The top 50% miRNA expression was defined as high expression group and the rest was defined as low expression group. The hazard ratio (HR) with 95% confidence intervals (CIs) was analyzed by the Cox proportional hazards regression model. AML = acute myeloid leukemia.

In this study, according to the sequencing data of the TCGA database, the DE-miRNAs (P value <.05 and |log2fold change (FC)|≥1) and DEGs (P value <.05 and |log2FC|≥1.5) were identified between RUNX1 mutation and wild-type AML. As a result, 27 DE-miRNAs (25 upregulated and 2 downregulated) and 561 DEGs (429 upregulated and 132 downregulated) were totally obtained from TCGA (Table 1). Five miRNAs (miR-151b, miR-151a-5p, let-7a-2-3p, miR-363-3p, miR-20b-5p) with prognostic significance in AML were chosen for further researches after overall survival analyses (Fig. 2, A–E and Supplementary Figure 1).

Table 1

Identification of DE-miRNAs and DEGs.

Type	Upregulated	Downregulated
miRNA	hsa-miR-455-3p, hsa-miR-30a-5p, hsa-miR-214-5p, hsa-miR-363-3p, hsa-miR-151b, hsa-miR-20b-5p, hsa-miR-203b-3p, hsa-miR-199a-5p, hsa-miR-130a-3p, hsa-miR-30a-3p, hsa-miR-99b-3p, hsa-miR-196a-5p, hsa-miR-126-3p, hsa-miR-483-5p, hsa-miR-10b-5p, hsa-miR-99b-5p, hsa-miR-5579-3p, hsa-miR-125a-5p, hsa-miR-4754, hsa-miR-126-5p, hsa-miR-139-5p, hsa-miR-153-5p, hsa-miR-151a-3p, hsa-miR-151a-5p, hsa-miR-551b-3p (n = 25)	hsa-miR-6718-5p, hsa-let-7a-2-3p (n = 2)
mRNA	DNTT, OPALIN, COBL, C3orf50, SDK2, TLL2, L3MBTL4, CHRDL1, HMGA2, DLGAP2, SCN3A, THSD7A, BAALC, KIAA1462, SETBP1, GPR173, DOCK1, PROM1, SLCO5A1, AOX1, LRP6, TSPAN7, BLNK, ZNF521, CD34, GNG11, FAM171B, NPR3, C5orf23, ABO, PTK2, PEAR1, CLEC4C, PACSIN1, ABCA9, MST1R, FAM196A, KIRREL, CD109, KCNK10, KCNA5, ABCA6, AR, FAM171A1, COBLL1, MAGED4B, HTR1F, CREG2, DSG2, DDR1, EPHB1, JAKMIP2, LIMCH1, LAMA3, F2RL1, APP, LILRA4, RASD1, C7orf58, BCL2L14, CRHBP, DUSP27, ATP10A, PPM1J, DLG3, tAKR, S100A16, SPRED1, CLEC9A, PGM5, CALN1, ROBO4, GNAI1, MYCT1, GLI3, GJA1, CDH9, CYP46A1, MN1, LMO1, IGF2, NRP1, MYEF2, LHX6, SPON1, C1orf21, HLA-DOA, DACH1, CALCRL, APCDD1, AKT3, DSCR6, AASS, ABCG2, LDOC1, CDC42BPA, HOPX, TNFRSF21, H1F0, MLLT3, LILRB4, BEND7, C10orf58, CD200, CD1E, FAT3, SHD, CNKSR2, IL12RB2, WWTR1, LOC100240735, MGAT5B, NOS2, C5orf20, RXFP1, CXCL12, XIRP2, NR5A2, LOC283761, ST8SIA5, ENAH, SMOC1, PSD3, ITGA9, SLITRK5, DNASE1L3, CCND1, MYLK, NOVA1, MPP6, IFI44L, EVC, SLC41A2, WFS1, KCNH8, ESR1, HS3ST1, IRF8, GBP4, GPR110, PPP1R9A, KIAA0125, EBF1, BMP3, GSTM5, CMKLR1, DOCK6, CUX2, OTOA, PFKFB2, IL15, CCR2, DCHS1, MFAP4, TRIM9, C13orf18, FLNB, MMP28, SALL2, NRP2, AHNAK2, ABCA8, GCOM1, TIFAB, ADAMTS10, MACROD2, DLK1, FLJ22536, CLEC14A, GGTA1, TGFBI, SPARC, ESM1, MYCL1, PTPRM, ABCA12, CPNE7, CXCR2P1, LRRC16A, FHL1, CCDC48, TCTEX1D1, TCF4, PLXNA4, EGF, LRP2BP, STL, C9orf106, PBX1, COL6A1, IL28RA, KCNA3, SNX25, CD1C, TLR7, ZNF462, PCSK6, KLHL13, EXPH5, SLC4A3, PTPRS, FAM160A1, RPH3A, STK32B, COL1A2, MPPED2, FAM66D, MS4A4A, LGMN, PHEX, SPIB, KIF17, MSR1, LAG3, NFIA, RAG2, LOC646982, PRTFDC1, BCAR3, FNBP1L, NECAB1, CLEC1B, P2RY14, ENAM, HLA-DQA2, PPP1R3G, PAWR, MAN1A1, RPS6KA2, EPSTI1, MYO5C, SVEP1, PCDH18, HLA-DPA1, PAR5, COL8A2, HLA-DQB2, SEMA3C, MRC1, KIAA1257, APOD, PCDHGC3, SHANK3, GLYATL1, SPRY1, HCRTR1, SLC41A1, AFF3, FMO2, PLEKHH2, PAR-SN, VPREB1, AJAP1, FRZB, DLL1, FAM81A, PARD3, GARNL3, SNORD116-4, PCDHGC4, PMEPA1, EMID2, ZNF608, HLA-DQA1, AMZ1, LRRC2, SH3BP4, NBEA, EVC2, PCDHGB7, VCAM1, SLAMF7, PKIA, HLA-DRA, HEMGN, ZNF135, PTGIR, EYA1, MTMR11, SYT6, C6orf192, TM4SF1, CYP4F2, IRF4, NOG, SCRN1, FAM129C, FOXO1, LPHN2, SNORD116-28, SERPINF1, FCRLA, SULF1, ZNF704, TMPRSS13, SLC4A4, CETP, ASAP3, SGCE, ADIPOQ, SGSM1, LGR6, DPPA3, LPAR6, BTLA, SMAD1, GABBR1, SLC8A1, TRPS1, LRP1, TET1, ITGA8, CDKL5, LY9, GZMB, HLA-DRB5, CHL1, OPTN, CECR7, IGFBP3, ARL5C, SHANK1, FGFR2, PLXNC1, APLNR, CARD11, AMPD1, MYO1E, DFNA5, SLC37A3, LOC284551, FLJ37543, PTPRO, BMP5, SLC38A1, SLC47A1, GFRA1, IGJ, ENHO, FAM64A, IPW, VANGL2, COL6A2, SMO, CIITA, ANGPTL2, SLC12A3, KIAA1984, SLC6A1, HLA-DPB1, EPB41L2, MPEG1, VGLL3, HLA-DRB1, ARNT2, CYP1B1, HLA-DMB, MX1, FAM134B, ABHD12B, LRRN1, MPDZ, BCAT1, MMP11, DCN, PPBP, TMEM117, CCL3L1, EPHA2, DFNB31, JAM2, PAPPA, COL5A1, CCR5, FLJ40330, MAGED4, ZNF99, DTX1, PAR1, RTN1, TTTY10, FAT4, PLEKHG4B, HFE, MDFIC, CCDC80, SYT1, CPLX2, SNORD116-20, SCT, HLA-DQB1, TMEM132B, CLIP4, HOMER2, CDH11, PIP5K1B, CACNA1H, EDNRB, C12orf42, PRO1768, PLOD2, DIRAS1, IGF2BP3, C6orf145, ME3, TMEM130, TINAGL1, FMO3, PNMAL2, ITGB3, C7, NRIP1, GRAMD1B, GREM1, LDHAL6B, PCDHGA11, NCRNA00200, MCF2L, CRMP1, DMD, PRKACG, DAGLA, APOB, HLF, FAM84A, IFITM3, MPP3, LRRK2, HOXC4 (n = 429)	LOC646851, IL5RA, PTGER3, WT1, PRICKLE1, MMP15, ITPKA, LOC441208, PTGER1, ANXA8L2, C1orf150, TFCP2L1, GPT2, TIE1, RNASE2, ADAMTS14, LGALS3BP, COL23A1, DYX1C1, MPV17L, SLC27A6, TOM1L1, C12orf59, LOC254559, C8orf79, ZNF114, RGS9BP, ADAMTS3, MAOA, KCTD15, WDR35, COL9A2, MAP1LC3A, CEBPE, GSDMC, WIT1, PRSSL1, SAGE1, HYAL3, C1orf106, NR1I2, DLC1, CCDC24, PAQR5, DHCR24, IL17RE, ZNF503, ANXA8, IGFBP2, CFD, CCL1, MS4A3, JAG1, MACC1, WWC2, C10orf114, CILP2, KDELC1, IRX5, POMC, HOXB8, TTLL10, S100P, KCTD1, SLC22A20, ZNF334, KIAA1958, KRT8, TBX1, KCNQ5, GPR27, WNT7B, FAM83H, LTC4S, CT45A1, DSC2, ANO7, MAMDC2, ABP1, C2, NTRK1, BIK, RETN, ARTN, REXO1L1, PROK2, ROPN1L, MOSC2, LRG1, S100B, FANK1, HOXB9, C3AR1, GTSF1, DEFB1, LPO, CYP7B1, RNASE3, C10orf11, KCNE1L, DPP10, PVR, SV2B, CLTCL1, GPR85, C21orf56, SPRY2, CLEC11A, SLC28A3, KRT18, LGALS12, SEC16B, C19orf51, KRT17, IGLL1, DPY19L2, NDST3, NTNG2, MPO, IRX3, PTX4, NXF3, CEACAM19, LTK, MOSC1, LPPR3, AZU1, ZNF804A, CTSG, CCNA1, APOC2, LOC728606 (n = 132)

DEGs = differentially expressed genes, DE-miRNAs = differentially expressed miRNAs.

GO and GSEA enrichment analysis of DEGs

The functional enrichment of candidate DEGs was estimated at the Metascape website. As demonstrated in a cluster heat map of the GO analysis (Fig. 3A), the upregulated DEGs suggested the principal enrichment in MHC class II protein complex, extracellular structure organization, blood vessel development, embryonic morphogenesis, cell morphogenesis involved in differentiation, regulation of cell adhesion, and etc. Similarly, the downregulated DEGs were mainly enriched in secretory granule lumen, extracellular structure organization. In addition, an enrichment network that demonstrated the relationships between the terms was built. In this network, each node symbolized a rich term that was colored by cluster ID (Fig. 3B), P value (Fig. 3C), and gene lists identification (Fig. 3D).

Figure 3

GO enrichment analysis of DEGs. A, Heatmap of enriched terms for downregulated and upregulated DEGs colored by P values was visualized by Metascape. B–D, Network of enriched terms: (B) colored by enriched terms, (C) colored by P value, (D) color-coded based on the identities of the gene lists. DEGs = differentially expressed genes, GO = gene ontology.

Moreover, top 20 clusters were performed in Table 2, showing their representative enriched terms. The GSEA analyzed the influence of RUNX1 mutations on diverse biological functional gene sets. In the GSEA analysis of KEGG pathways, the RUNX1 mutation related to adherent junction, WNT signaling pathway, MAPK signaling pathway, JAK-STAT signaling pathway, cell adhesion molecules CAMs (Table 3 and Supplementary Figure 2).

Table 2

Top 20 clusters with their representative enriched terms.

Type	GO	Category	Description	Count	%	Log10(P)	Log10(q)
↑	GO:0042613	GO cellular components	MHC class II protein complex	11	2.62	−15.81	−11.46
↓↑	GO:0043062	GO biological processes	Extracellular structure organization	40	7.27	−13.29	−9.24
↑	GO:0001655	GO biological processes	Urogenital system development	34	6.18	−12.61	−8.86
↓↑	GO:0001568	GO biological processes	Blood vessel development	54	9.82	−12.43	−8.86
↑	GO:0000904	GO biological processes	Cell morphogenesis involved in differentiation	45	10.71	−12.15	−8.27
↓	GO:0048598	GO biological processes	Embryonic morphogenesis	42	7.64	−9.85	−6.83
↑	GO:0001501	GO biological processes	Skeletal system development	39	7.09	−9.84	−6.83
↑	GO:0030155	GO biological processes	Regulation of cell adhesion	39	9.29	−9.44	−6.49
↑	GO:0048589	GO biological processes	Developmental growth	37	8.81	−8.97	−6.14
↑	GO:0003094	GO biological processes	Glomerular filtration	8	1.9	−8.52	−5.76
↑	GO:0008285	GO biological processes	Negative regulation of cell proliferation	39	9.29	−8.5	−5.75
↑	GO:0061005	GO biological processes	Cell differentiation involved in kidney development	12	2.18	−8.32	−5.56
↑	GO:0098794	GO cellular components	Postsynapse	34	8.1	−8.07	−5.39
↑	GO:0009880	GO biological processes	Embryonic pattern specification	11	2.62	−8.05	−5.38
↑	GO:0045177	GO cellular components	Apical part of cell	26	6.19	−7.87	−5.22
↑	GO:0050808	GO biological processes	Synapse organization	31	5.64	−7.68	−5.07
↓↑	GO:0009792	GO biological processes	Embryo development ending in birth or egg hatching	40	7.27	−7.56	−4.97
↑	GO:0030425	GO cellular components	Dendrite	32	7.62	−6.97	−4.44
↓	GO:0034774	GO cellular components	Secretory granule lumen	12	9.23	−6.76	−4.24
↑	GO:0009611	GO biological processes	Response to wounding	33	7.86	−6.72	−4.21

↓ = downregulated genes, ↑ = upregulated genes, GO = gene ontology.

Table 3

GSEA results of RUNX1 mutations in AML patients.

NAME	ES	NES	NOM p-val	FDR q-val
KEGG_ADHERENS_JUNCTION	0.49826324	1.7916677	0.001834862	0.18580072
KEGG_DORSO_VENTRAL_AXIS_FORMATION	0.5595871	1.7560349	0.005714286	0.15923245
KEGG_WNT_SIGNALING_PATHWAY	0.38487324	1.5906051	0.005747126	0.21179874
KEGG_TASTE_TRANSDUCTION	0.5330845	1.7734187	0.006465518	0.16741094
KEGG_ONE_CARBON_POOL_BY_FOLATE	-0.64856833	-1.7047765	0.010799136	0.24881025
KEGG_GAP_JUNCTION	0.45841148	1.7137706	0.010869565	0.17216179
KEGG_COLORECTAL_CANCER	0.39059794	1.536395	0.012544803	0.24306875
KEGG_LONG_TERM_DEPRESSION	0.41905212	1.6034809	0.01532567	0.24336743
KEGG_JAK_STAT_SIGNALING_PATHWAY	0.39050087	1.576051	0.016627079	0.22054975
KEGG_CALCIUM_SIGNALING_PATHWAY	0.3339285	1.4759699	0.018404908	0.22962376
KEGG_PATHWAYS_IN_CANCER	0.32847777	1.4595866	0.018832391	0.24317616
KEGG_GLYOXYLATE_AND_DICARBOXYLATE_METABOLISM	-0.6353147	-1.6713699	0.022177419	0.24820589
KEGG_ASTHMA	0.69502634	1.7460698	0.0251938	0.14858118
KEGG_ALDOSTERONE_REGULATED_SODIUM_REABSORPTION	0.4761012	1.5969595	0.028355388	0.21735093
KEGG_AUTOIMMUNE_THYROID_DISEASE	0.56841636	1.6766285	0.02955665	0.2084749
KEGG_TYPE_I_DIABETES_MELLITUS	0.6203307	1.6638426	0.030075189	0.18521671
KEGG_REGULATION_OF_ACTIN_CYTOSKELETON	0.36554304	1.4903752	0.030828517	0.24551603
KEGG_CELL_ADHESION_MOLECULES_CAMS	0.47629744	1.6032443	0.031189084	0.22495992
KEGG_INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION	0.6331871	1.6653942	0.034026466	0.20303084
KEGG_MAPK_SIGNALING_PATHWAY	0.33732373	1.4583732	0.037924152	0.23690787
KEGG_ALLOGRAFT_REJECTION	0.66069037	1.6163293	0.038910504	0.24292912
KEGG_DILATED_CARDIOMYOPATHY	0.37348428	1.4308964	0.049586777	0.24564306

AML = acute myeloid leukemia, ES = enrichment score, FDR = false discovery rate, GSEA = gene set enrichment analysis, KEGG = Kyoto Encyclopedia of Genes and genomes.

Construction of PPI network

To further explore the hub genes and interaction of those DEGs, protein–protein interaction was created through Search Tool for the Retrieval of Interacting Genes, extracted for visualization by Cytoscae (Fig. 4A). In this protein–protein interaction network, nodes were color-coded according to the expression of DEGs, for example, red represented upregulated DEGs and green represented downregulated DEGs. Including 12 algorithms, cytoHubba plug-in of Cytoscae detected those top 20 hub genes (Supplementary Table 1). Recognized by more than 7 algorithms, 11 hub genes (PPBP, APP, CCR5, HLA-DRB1, GNAI1, APLNR, P2RY14, C3AR1, HTR1F, CXCL12, GNG11) with higher degree of connectivity were chosen to construct the hub gene protein–protein interaction network (Fig. 4B). The enrichment analysis demonstrated that those hub genes enriched in positive regulation of leukocyte migration, side of membrane, positive regulation of leukocyte migration, second-messenger-mediated signaling, as well as chemical synaptic transmission. After that, imported into Cytoscape software, the PPI networks were analyzed via plug-ins Molecular Complex Detection. Ultimately, 9 modules in PPI network were detected, then we selected 3 significant modules (Fig. 4C). The enrichment analysis illustrated that those DEGs in modules primarily enriched in chemokine-mediated signaling pathway, MHC class IIprotein complex, positive regulation of cytokine production, regulation of cytosolic calcium ion concentration, regulation of T-cell activation, blood circulation, and extracellular matrix structural constituent conferring tensile strength, positive regulation of T-cell migration, and so on. Eleven genes (PPBP, APP, CCR5, HLA-DRB1, GNAI1, APLNR, P2RY14, C3AR1, HTR1F, CXCL12, GNG11) were simultaneously distinguished through module analysis and hub gene analysis.

Figure 4

PPI network, hub genes network, and modules analyses of DEGs. A, PPI network of DEGs. B, The network of 11 hub genes with a higher degree of connectivity and enrichment analysis of these genes. C, Genes of top 3 modules were performed GO enrichment analysis by Metascape. Nodes were color-coded based on the expression of DEGs (red, upregulated; green, downregulated). DEGs = differentially expressed genes, GO = gene ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes, PPI = protein–protein interaction.

MiRNA-mRNA regulatory network analysis

The target genes of 5 candidate DE-miRNAs were identified by Target Scan, miRWalk, miRDB, and miRPathDB (Fig. 5A), then matched with DGEs (P < .05, not limited fold change) (the target genes of upregulated miRNA were adapted to downregulated DEGs, and the target gens of downregulated miRNA were adapted to upregulated DEGs). As a result, we identified 190 target genes that expressed differently between RUNX1 mutation and wild type in AML. Here, these DE-miRNA (n = 5) and target genes (n = 190) were utilized to construct the miRNA-mRNA regulatory network, showing that let-7a-2-3p and miR-363-3p had rich external connections (Fig. 5B).

Figure 5

miRNA-mRNA regulatory network. A, Diagrams illustrating in prediction of target genes of 5 candidate miRNAs. B, miRNA-mRNA regulatory network. Nodes were color-coded based on components (red, DE-miRNAs; green, target DEGs). The red lines indicated the regulation relationship between DE-miRNAs and their targets DGEs. The silver lines indicated the relationship between DEGs. DEGs = differentially expressed genes, DE-miRNAs = differentially expressed miRNAs.

Target genes analysis of miR-363-3p

In a previous study, the high expression of let-7a-2-3p could be potentially used as favorably prognostic biomarkers in cytogenetically normal AML patients,[ while the differential expression of miR-151b was regarded as a prognostic signal in upper tract urinary carcinoma,[ and a blood-based biomarker for diagnosing ischemic stroke patients.[ MiR-151a-5p and miR-151b were significantly downregulated in the amyotrophic lateral sclerosis[ and also were regarded as prognostic biomarkers in the blood of primary CNS lymphoma patients.[ MiR-151a-5p was significantly over-expressed in colorectal cancer[ and lung cancer.[ MiR-20b-5p played an significant role in prostate cancer(Pca),[ lung cancer,[ laryngeal squamous cell carcinoma,[ renal cell carcinoma,[ cancer stem cells,[ and so on. Nevertheless, there were a few studies involving miR-363-3p in AML. Therefore, miR-363-3p was selected for a further study.

According to the miRNA-mRNA regulatory network, these were 19 target genes associated with miR-363-3p: SPRY-domain containing protein 4 (SPRYD4), PLEKHB2, ZNF385B, PIAS4, ASPH, ZDHHC5, ITGA5, GDF11, TECPR2, PIP4K2C, E2F3, Fibronectin type III domain containing 3B (FNDC3B), DPP10, FAM16A2, WWC1, PLXDC2, DOCK5, DSC2, CDK16, all of them were downregulated genes in RUNX1 mutation AML of TCGA. Simultaneously, performed between 19 target genes and miR-363-3p, the correlation analysis showed that miR-363-3p had negative correlations with SPRYD4, FNDC3B, E2F3 (P < .05, Fig. 6, A–C). Moreover, SPRYD4 and FNDC3B had prognostic significance in overall survival of AML patients (P < .05, Fig. 6, D–E), except E2F3 (Fig. 6F). The expression of SPRYD4 and FNDC3B were both decreased in AML with RUNX1 mutation (Fig. 6, G–I).

Figure 6

SPRYD4 and FNDC3B may be the target genes of miR-363-3p involved in AML with RUNX1 mutation. The correlation analyses showed miR-363-3p had significant negative correlations with SPRYD4 (Pearsons’ rho = −0.17), FNDC3B (Pearsons’ rho = −0.33), and E2F3 (Pearsons’ rho = −0.21) in AML with RUNX1 mutation from TCGA. Survival analyses showed that SPRYD4 and FNDC3B had better overall survival (OS) of AML patients from TCGA and had different mRNA expression levels between RUNX1 mutation and wild-type AML (P < .05), except E2F3. OS = overall survival, TCGA = The Cancer Genome Atlas.

Association between RUNX1 mutation and clinicopathological parameters

A total of 197 AML patients were obtained from TCGA, including 18 RUNX1 mutation and 179 RUNX1 wild type. χ2analysis was utilized to estimate the relationship between RUNX1 mutation and clinicopathological parameters. As a result, the differences of RUNX1 mutation rates in sex, disease-free status, bone marrow blast percentage, and peripheral blood blast percentage were not significant (P > .05) between the 2 groups. However, RUNX1 mutation was obviously correlated with age (P = .027), WBC (P = .021), FAB type (P = .008), risk(cyto) (P = .009), and risk(molecular) (P = .02) (Table 4).

Table 4

The analysis of clinical and pathological variables between RUNX1 mutation and wild-type patients.

	RUNX1
Variable	Mutation (18)	Wild type (179)	Frequency (197)	P value
Age (yr)
≤60	6	108	114	.027
>60	12	71	83
Sex
Male	11	94	105	.486
Female	7	85	92
FAB
M0	7	11	18	.008
M1	4	42	46
M2	2	42	44
M3	0	20	20
M4	5	34	39
M5	0	22	22
M6	0	3	3
M7	0	3	3
Unknown	0	2	2
Disease free status
Disease free	10	91	101	.275
Recurred/progressed	8	85	93
Unknown	0	3	3
risk (cyto)
Good	0	36	36	.009
Intermediate	16	98	114
Poor	1	41	42
Unknown	1	4	5
Risk (molecular)
Good	0	38	38	.02
Intermediate	16	90	106
Poor	1	48	49
Unknown	1	3	4
WBC
≥10	6	110	116	.021
<10	12	69	91
BM blast percentage
≥50	13	145	158	.561
<50	5	34	39
PB blast percentage
≥50	6	68	74	.633
<50	12	106	118

BM = bone marrow, PB = peripheral blood, WBC = white blood cell.

Discussion

Regarded as the majority common acute leukemia in adults,[ AML is a malignant clonal disease that promotes the growth of malignant cells, instead of producing healthy hematopoietic cells.[ As a transcription factor, RUNX1 regulates critical processes in various aspects of hematopoiesis, while t(8,21) is identified as a chromosomal translocation in acute myeloid leukemia.[ Our study showed that RUNX1 expression was higher in DNMT3A mutation AML compared with wild-type AML. The previous study has shown that RUNX1 mutations were remarkably associated with poor clinical outcomes as well as increased risk of death in AML patients (log-rank P = .012, HR = 2.145, 95% CI = 1.156–3.982).[ Stengel et al[ indicated that RUNX1-mutated AML showed significant genetic abnormalities and poor prognosis. Wang et al[ showed that RUNX1 mutations in AML were associated with distinct clinical features. Thus, it was significant to further understand the biological roles of RUNX1 mutations in AML.

According to our study, a comprehensive bioinformatics analysis was conducted to explore the crucial genes and miRNA of acute myeloid leukemia with RUNX1 mutation. In this study, 27 DE-miRNAs (25 upregulated and 2 downregulated) and 561 DEGs (429 upregulated and 132 downregulated) were identified between RUNX1 mutation and wild-type AML patients who were obtained from TCGA database. The GO analysis showed that upregulated DEGs suggested significant enrichment in MHC class II protein complex, extracellular structure organization, blood vessel development, cell morphogenesis involved in differentiation, embryonic morphogenesis, regulation of cell adhesion, and so on. These were associated with the cancer. Similarly, the downregulated DEGs were mostly enriched in secretory granule lumen, extracellular structure organization. It was suggested that RUNX1 mutations regulated the development and prognosis of AML by blood vessel development, cell morphogenesis involved in differentiation, embryonic morphogenesis, regulation of cell adhesion. In the GSEA analysis of KEGG pathways, the RUNX1 mutation related to adherent junction, WNT signaling pathway, MAPK signaling pathway, JAK-STAT signaling pathway, cell adhesion molecules CAMs. Consistent with previous studies, those pathways had a crucial role in tumor progression.[

In addition, we constructed the PPI network, then PPBP, APP, CCR5, HLA-DRB1, GNAI1, APLNR, P2RY14, C3AR1, HTR1F, CXCL12, GNG11 were selected as the hub genes according to the degree of connectivity. The enrichment analysis showed that those hub genes enriched in cellular calcium ion homeostasis, side of membrane, positive regulation of leukocyte migration, second-messenger-mediated signaling, as well as chemical synaptic transmission. Eventually, 9 modules in PPI network were detected, among which we selected 3 significant modules, mainly enriched in chemokine-mediated signaling pathway, MHC class II protein complex, positive regulation of cytokine production, regulation of cytosolic calcium ion concentration, regulation of T-cell activation, blood circulation and extracellular matrix structural constituent conferring tensile strength, positive regulation of T-cell migration, and so on.

In addition, the results of log rank test by overall survival analysis showed that 5 miRNAs (miR-151b, miR-151a-5p, let-7a-2-3p, miR-363-3p, miR-20b-5p) had a prognostic significance in AML. The miRNA-mRNA regulatory network was built by those miRNAs and their target genes. In a previous study, the high expression of let-7a-2-3p could be potentially used as favorably prognostic biomarkers in cytogenetically normal AML patients,[ while the differential expression of miR-151b was regarded as a prognostic signal in upper tract urinary carcinoma,[ a blood-based biomarker for diagnosing ischemic stroke patients.[ MiR-151a-5p and miR-151b were significantly downregulated in the amyotrophic lateral sclerosis[ and were regarded as prognostic biomarkers in the blood of primary CNS lymphoma patients.[ MiR-151a-5p was significantly over-expressed in colorectal cancer[ and lung cancer.[ The miR-20b-5p played an significant role in prostate cancer,[ lung cancer,[ laryngeal squamous cell carcinoma,[ renal cell carcinoma,[ cancer stem cells,[ and so on. Nevertheless, there were few studies involving miR-363-3p in AML. So, we chose miR-363-3p for further study, which had 19 target genes: SPRYD4, PLEKHB2, ZNF385B, PIAS4, ASPH, ZDHHC5, ITGA5, GDF11, TECPR2, PIP4K2C, E2F3, FNDC3B, DPP10, FAM16A2, WWC1, PLXDC2, DOCK5, DSC2, CDK16, all of them were downregulated genes in RUNX1 mutation AML of TCGA. After correlation analysis and survival analysis, SPRYD4 and FNDC3B had prognostic significance in AML patients.

Moreover, inducing the death of apoptotic cell, SPRYD4 restricted the progression in hepatocellular carcinoma.[ Nevertheless, there was less study on SPRYD4, which required a further research.

Fibronectin type III domain containing 3B (FNDC3B/ FAD104) was the member of FNDC3 family, promoting epithelial–mesenchymal transition of tongue squamous cell carcinoma,[ colorectal cancer progression,[ hepatocellular carcinoma,[ and glioblastoma.[ Moreover, a study showed that upregulated expression of miR-143 repressed FNDC3B expression, which promoted the metastasis of prostate cancer.[ When the FNDC3B was knockdown in NB4 cells, the NB4 cell proliferated faster. It meant that FNDC3B inhibited the proliferation of acute myeloid leukemia cell. Here, FNDC3B was regarded as a fusion gene in AML.[ The results were consistent with our expectation that downregulated FNDC3B promoted the development of AML patients.

According to the AML patients from TCGA, RUNX1 mutation rates in patients over 60 years of age were significantly higher than those under 60 years of age (P < .05). While the white blood cell was at high level, the RUNX1 mutation rates in patients were lower than those with low level (P < .05). Moreover, the RUNX1 mutation rates in M0 were higher than other types in FAB, the differences were statistically significant (P < .05). The RUNX1 mutation rates in patients at a good rick rank (cyto and molecular) were lower than other rick ranks (P < .05). Nevertheless, the differences of RUNX1 mutation rates in sex, disease-free status, bone marrow blast percentage, and peripheral blood blast percentage were not significant (P > .05). In summary, this indicated that RUNX1 mutation was obviously correlated with age, WBC, FAB type, and risk(cyto and molecular). In conclusion, the identified DE-miRNAs and DEGs might play a key role in acute myeloid leukemia with RUNX1 mutation. MicroRNA-363-3p may promote the development of RUNX1 mutation AML, targeting SPRYD4 and FNDC3B.

Furthermore, further studies are required to support our results, such as cell culture and patient samples.

Conclusion

RUNX1 is regarded as one of the most frequently mutated genes in acute myeloid leukemia. Our study has pointed that multiple genes and microRNAs may play a crucial role in RUNX1 mutation AML. RUNX1 mutations regulated the development and prognosis of AML by blood vessel development, cell morphogenesis involved in differentiation, embryonic morphogenesis, regulation of cell adhesion. In the GSEA analysis of KEGG pathways, the RUNX1 mutation related to some pathways, which had a crucial role in tumor progression. Moreover, we find that microRNA-363-3p may promote the proliferation of RUNX1 mutation AML by targeting SPRYD4 and FNDC3B. However, further experiments, like the cell culture and patient samples, are still needed to confirm our results.

Acknowledgments

The authors thank the TCGA database for their sharing of the AML sequencing dataset.

Author contributions

Conceptualization: Huo Tan, Lihua Xu.

Data curation: Yimin Chen, Shuyi Chen.

Formal analysis: Yimin Chen, Shuyi Chen, Danyun Yuan.

Funding acquisition: Huo Tan, Lihua Xu.

Investigation: Jielun Lu, Pengfei Qin.

Methodology: Lang He.

Project administration: Lihua Xu.

Resources: Lang He, Pengfei Qin.

Software: Yimin Chen.

Supervision: Lihua Xu.

Validation: Huo Tan, Lihua Xu.

Visualization: Huo Tan, Lihua Xu.

Writing – original draft: Yimin Chen, Shuyi Chen.

Writing – review & editing: Jielun Lu, Danyun Yuan.