Construction of circRNA–miRNA–mRNA Network Reveal Functional circRNAs and Key Genes in Acute Myeloid Leukemia

Introduction

Acute myeloid leukemia (AML) is a group of highly aggressive hematological malignancies with abnormal proliferation and differentiation of hematopoietic stem cells.¹ Abnormal blood cells proliferate in a considerable number and inhibit normal hematopoiesis, with the clinical manifestations of anemia, bleeding, hepatosplenomegaly and other organ infiltration. It is one of the most common hematologic malignancies, and the mortality rate is the highest compared with other types of leukemia.² The treatment is based on conventional chemotherapy and hematopoietic stem cell transplantation, and molecular targeted therapy and cellular immunotherapy are promising for treating AML patients. However, they are still at risk for drug resistance and disease recurrence, thus resulting in shorter survival and poorer quality of life.³ Studying the molecular mechanism of AML development takes on a great significance in finding novel therapeutic targets and facilitating AML prognosis.

Circular RNAs (circRNAs) are a special type of closed-loop non-coding RNA molecule without 5’-caps or 3’ -polyadenylate tails. Existing research has suggested that circRNAs exhibit complex biological functions in the body.⁴ An increasing number of studies have suggested that circRNAs play a vital role in the development of tumors.⁵ Existing research has also suggested that circRNA is correlated with the emergence of drug-resistance in osteosarcoma,⁶ breast cancer,⁷ as well as chronic myeloid leukemia.⁸ CircRNA is stably expressed and is present at high level in body fluids, making it an ideal candidate for noninvasive liquid biopsy biomarkers in tumors.⁹ Existing research has suggested that circRNAs are capable of playing a role in the post-transcriptional regulation of AML by binding miRNAs, activating downstream signaling cascades, and regulating the expression of related genes.¹⁰ CircRNA is closely correlated with a wide variety of processes of AML and may serve as a novel target for AML diagnosis and treatment. With the development of bioinformatics technology, there have been increasing circRNA-based bioinformatics studies, whereas the molecular mechanism and the specific function of the novel circRNAs in AML should be explored in depth.

In this study, the microarray chip data in AML were analyzed and mined from Gene Expression Omnibus (GEO) database, differentially expressed circRNAs (DEcircRNAs), differentially expressed miRNAs (DEmiRNAs), and differentially expressed mRNAs (DEmRNAs) were selected using R software. circRNA/miRNA/mRNA interaction networks were built using Cytoscape. Moreover, the main functional pathways were evaluated through gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Lastly, validation and survival analysis of hub genes were conducted using The Cancer Genome Atlas (TCGA) database and genotype-tissue expression (GTEx) database. Subsequently, we performed real-time quantitative polymerase chain reaction (RT-qPCR) to confirm that these circRNAs and miRNAs are expressed differently in AML and healthy controls. The above results provide new directions for the pathogenesis and treatment of AML.

Materials and Methods

Data Acquisition and Identification of DEcircRNAs

Gene Expression Omnibus (GEO) database contains high-throughput gene expression data submitted by research institutions around the world. Microarray chip data from AML patients (GSE94591) in the GEO database were selected. R language was used to screen differential genes between AML patients and healthy controls for significant differential circRNAs (P<0.05, logFC absolute value > 2).

Prediction of miRNA Response Elements (MREs)

Cancer-Specific CircRNA Database (CSCD) explores the function of circRNAs by predicting the miRNA response elements (MREs) loci. We predicted the target genes of circRNAs by CSCD. The microarray data from AML patients in the GEO database (GSE142700-GPL26945) were selected, and the R language was employed to screen for meaningful miRNAs (P<0.05, logFC absolute value>1). By intersecting the circRNA-predicted target miRNAs with them, we selected 8 overlapping miRNAs for further mRNA prediction.

Forecast of miRNA Targeted Genes, Collection of AML Differentially Expressed Genes (DEG) and Obtaining Overlapping Genes

MiRanda was developed by researchers at the prestigious Memorial Sloan-Kettering Cancer Research Center to predict the binding sites of miRNAs interacting with their target genes. MiRTarBase is a manually collected database of experimentally verified miRNA target genes. For each miRNA target gene record, a unique miRNA-target interaction (referred to as MTs) number will be assigned. The target gene mRNA of miRNA was predicted by miRanda and miRTarBase software. Similarly, the microarray data from AML patients in the GEO database (GSE142700-GPL19956) were selected, and the R language was employed to screen for meaningful mRNAs (P<0.05, logFC absolute value>1). The genes between the DEGs and predicted miRNA target genes were obtained through the Venn diagram.

Functional Enrichment Analysis of Overlapped Genes

DAVID V6.8 (https://david.abcc.ncifcrf.gov/), an online bioinformatics resource, provides tools for functional interpretation of a large number of genes, which was used for Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of overlapping genes. The setting was adopted to perform the R package, with p<0.05 and the count is>2.

Overall Survival (OS) Analysis of the Hub Genes

GEPIA 2 (Gene Expression Profiling Interactive Analysis) mainly focuses on The Cancer Genome Atlas (TCGA) database and Genotype-Tissue Expression (GTEx) project. Analysis includes tumor/normal differential expression profiling, expression analysis, survival analysis, similar genes, gene expression correlation and dimensionality reduction analysis. Overall survival (OS) analysis was used to explore the effect of differentially expressed hub genes on OS between healthy controls and AML patients.

Construction of circRNA–miRNA–mRNA Network

A circRNA-miRNA-mRNA regulatory network was constructed by mapping the regulatory networks using Cytoscape software (Figure 1).

Figure 1 Flowchart of this study about constructing a circRNA/miRNA/mRNA network.

AML Patients and Healthy Controls

We enrolled 43 patients with confirmed AML and 35 healthy controls who were admitted to the Second Affiliated Hospital of Zhejiang Chinese Medical University from Jan 2021 to Sep 2022. There was no significant difference in gender and age between patients with AML and healthy controls. Bone marrow aspirate was performed on two groups to collect BM cells. Bone marrow mononuclear cells (BMMNCs) were isolated using lymphocyte isolation medium (TBD Science).

RNA Extraction and Quantitative Real-Time PCR Validation

Total RNA were extracted from BMMNCs of AML patients and controls with TRIzol reagent (Life Technologies, CA, USA). The quantitative analysis of RNAs was performed by NanoDrop 2000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA). For circRNA, cDNA was synthesized using the PrimeScript™ RT reagent kit (Takara, Tokyo, Japan) under the guidance of the manufacturer’s instructions. Realtime-qPCR was performed using the TB Green® Premix Ex Taq™ II kit (Takara, Tokyo, Japan). PCR primers are listed in Table 1. For microRNA, cDNA was synthesized using the Mir-X miRNA First-Strand Synthesis Kit (Takara, Tokyo, Japan). The special sequences of hsa-miRNAs were amplified with the forward primer in Table 1. The universal reverse primer and U6 primers in Mir-X miRNA qRT-PCR TB Green Kit (Takara, Tokyo, Japan). Indicated genes were normalized to β-actin or U6. Thereafter, the folding change of genes was calculated by the 2-ΔΔCt analytic method.

Table 1 Sequences of Primers for RT-qPCR and miRNA Related Sequence

Ethical Statement

This study were approved by the Ethics Committee of The Second Affiliated Hospital of Zhejiang Chinese Medical University (the Institutional Review Board approval number: 2021-KL-125-01) and were carried out in accordance with the ethical principles of the Declaration of Helsinki. All patients signed the informed consent form.

Statistical Analysis

Bioinformatic statistical analyses were performed using R software (v.3.6.1) and the above software packages. Other statistical analyses were done using GraphPad Prism software. Correlations were analyzed by Pearson’s correlation coefficient. P < 0.05 was statistically significant.

Result

Acquiring 10 DEcircRNAs in Acute Myeloid Leukemia

To explore the function of circRNAs and construct a circRNA-miRNA-mRNA interaction network, we obtained the microarray dataset GSE94591 from the GEO database to identify the differentially expressed circRNAs in bone marrow specimens from healthy controls and AML patients. The gene chip was obtained from the Agilent-069978 Arraystar Human CircRNA microarray V1 platform. Using the Limma package, we obtained 4520 differentially expressed CircRNAs (Figure 2a). After screening at P<0.05 and differential fold > 4, we finally found that 72 CircRNAs matched and used the top 50 CircRNAs with the most significant differences for the heat map (as shown in Figure 2b). The top 5 circRNAs with the most significant rise and fall were used for the follow-up study. 10 circ RNAs were differentially expressed in bone marrow specimens from healthy controls and AML patients (Figure 2). The basic characteristics of the 10 circRNAs are shown in Table 2.

Table 2 Basic Characteristics of These 10 circRNAs

Figure 2 Acquired DEcircRNAs of AML. (a) Volcano plots of DEcircRNAs, the blue dots and red dots represent down- and up-expressed circRNAs. (b) Heatmap of circRNAs. Type: (C) healthy controls. (T) AML patients.

Identification of 10 circRNA-miRNA Interactions

To explore circRNA-miRNA interactions, we previously analyzed the structures of 10 circRNAs (Figure 3a) which have the structures of MRE, RBP and ORF. We also predicted the target miRNAs to which CSCD binds. Meanwhile, we identified differentially expressed miRNAs between healthy controls and AML patients. The microarray dataset GSE142700 was obtained from the GEO database. Gene chips were obtained from NanoString nCounter Human miRNA platform. Using the Limma package, 777 differential circRNAs were screened (Figure 3b). After screening with P<0.05 and the differential fold > 2-fold, 54 miRNAs were finally found to be concordant, and the top 50 miRNAs with the most significant differences were mapped to a heat map (Figure 3c). circRNA-predicted target genes intersected with miRNAs screened by the AML patient microarray (Figure 3d), and 6 miRNAs intersections were identified, namely: hsa-miR-382-5p, hsa-miR-3190-3p, hsa-miR-150-5p, hsa-miR-1193, hsa-miR-3613-5p, hsa-miR-454-3p. Six miRNAs were used as follow-up studies.

Figure 3 Differentially expressed miRNAs in AML patients. (a) Basic structures of the 10 circRNAs. (b) Volcano plots for differentially expressed miRNAs, the blue and red dots represent down- and up-expressed miRNAs. (c) Heatmap of differentially expressed miRNAs. Type: (C) healthy controls. (T) AML patients. (d) Venn graph showing the 6 overlapped genes by intersecting circRNA targeted genes and miRNA response elements from AML patients.

Identification of 6 miRNA-mRNA Interactions

We used the miRanda and miRTarBase tools to predict the target mRNAs of the 6 miRNAs in order to investigate miRNA-mRNA interactions. At the same time, differentially expressed mRNAs between healthy controls and AML patients were identified. The microarray dataset GSE142700 was obtained from the GEO database. The gene chip was obtained from the NanoString nCounter human PanCancer Pathways Panel platform. Using the Limma package, 779 differential mRNAs were screened (as shown in Figure 4a). After screening with P<0.05 and the differential fold >2, 137 matched mRNAs were finally found, with the top 50 mRNAs with the most significant differences used to plot a heat map (Figure 4b). The target genes predicted by miRNA were intersected with the mRNAs screened by the AML patient microarray 33 mRNAs were found to intersect (Figure 4c) and used as follow-up studies.

Figure 4 Differentially expressed mRNAs in AML patients. (a) Volcano plots of differentially expressed mRNAs, the blue and red dots represent down- and up-expressed miRNAs. (b) Heatmap for differentially expressed miRNAs. Type: (C) healthy controls. (T) AML patients. (c) Venn graph showing the 33 overlapped genes by intersecting the miRNA targeted genes and mRNA targets from AML patients.

Function Enrichment Analyses

For the potential functions and biological roles of the 33 overlapping genes, we performed GO analysis, including biological process (BP), molecular function (MF) and cellular component (CC) and KEGG signaling pathway enrichment analysis. GO analysis showed that, in terms of BP, these overlapping genes were mainly enriched in response to oxygen levels, regulation of cell-cell adhesion, response to alcohol, regulation of leukocyte cell-cell adhesion, and cellular response to the abiotic stimulus. For MF, overlapping genes are mainly enriched in cytokine receptor binding, protein serine/threonine kinase activity, phosphatase binding, kinase regulator activity and transforming growth factor beta receptor binding. KEGG signaling pathway enrichment analysis revealed that it was mainly enriched in PI3K-Akt signaling pathway, Pancreatic cancer, colorectal cancer, AGE-RAGE signaling pathway in diabetic complications, hepatocellular carcinoma and other pathways (Figure 5).

Figure 5 Dot plot of function enrichment analysis. (a) Biological process. (b) Cellular component. (c) Molecular function. (d) KEGG pathway.

CircRNA/miRNA/mRNA Competing for Endogenous RNA (ceRNA) Network Construction

Based on the circRNAs, miRNAs, and mRNAs screened in the previous stage, an association network between circRNA, miRNA and hub genes were established by Cytoscape (Figure 6). circRNA/miRNA/mRNA regulatory axis, including hsa_circ_0009581/hsa-miR-150-5p/RAD51 axis, hsa_circ_0009581/hsa-miR-150-5p/CDK4 axis, hsa_circ_0009581/hsa-miR-150-5p/PKMYT1 axis, hsa_circ_0009581/hsa-miR-150-5p/MYB axis, hsa_circ_0009581/hsa-miR-150-5p/HSP90B1axis, hsa_circ_0009581/hsa-miR-150-5p/ZNF346 axis, hsa_circ_0009581/hsa-miR-150-5p/PIK3R1 axis, hsa_circ_0001947/hsa-miR-454-3p/CREB5 axis, hsa_circ_0001947/hsa-miR-454-3p/BMP6 axis, hsa_circ_0001947/hsa-miR-454-3p/TGFBR2 axis, hsa_circ_0007765/hsa-miR-382-5p/CASP3 axis, hsa_circ_0004277/hsa-miR-3190-3p/GSK3B axis, hsa_circ_0004277/hsa-miR-3190-3p/SETBP1axis, hsa_circ_0004277/hsa-miR-3190-3p/SMAD2 axis, were extracted from this network.

Figure 6 ceRNA network. circRNAs, and miRNAs, hub-genes are shown in purple, yellow, and red respectively. circRNAs: dark purple: rising, light purple: falling; miRNAs: dark yellow: rising, light yellow: falling; hub-genes: dark red: rising, light red: falling.

Effects of Hub-Genes on Overall Survival

Hub genes have been discovered. The overall survival (OS) time of patients with AML will now be shown, which is the biggest concern for all doctors. As shown in Figure 7a, the hub genes, MYB, HSP90B1, PIK3R1, SMAD2, PKMYT1, RAD51 and TGFBR2 were significantly differentially expressed between normal controls and AML patients from TCGA database and GTEx project. The overall survival of RAD51 is shown in Figure 7b, and higher expression of RAD51 is associated with poor overall survival of AML patients (HR = 2.2, P = 0.044).

Figure 7 Survival analysis and validation of the 7 genes. (a) mRNA expression levels of 7 hub genes.*p<0.05. (b) The overall survival analysis of RAD51 related to AML patients. P < 0.05.

Validation of Differentially Expressed circRNAs and miRNAs by qPCR

To investigate the expression of circRNA and miRNA from bioinformatic analysis and determine whether they may be important to leukemia, we used RT-qPCR to confirm that some circRNA and miRNA are expressed differently in AML and healthy controls. The different expression levels of 10 circRNAs and 6 miRNAs from AML patients (n=43) and healthy controls (n=35) were analyzed using RT-qPCRs. The results showed that hsa_circ_0009581 and hsa_circ_0005273 were significantly upregulated in AML patients; hsa_circ_0000497 and hsa_circ_0001947 were significantly downregulated in AML patients compared to that in the controls (Figure 8a, p <0.05). Hsa-miR-150-5p was significantly downregulated; hsa-miR-454-3p was upregulated in AML patients (Figure 8b, p<0.05). To explore the relationship between circ RNAs and miRNAs, hsa_circ_0009581 and hsa-miR-150-5p; hsa_ circ_0001947 and hsa-miR-454-3p were inversely correlated using Pearson’s correlation coefficient (Figure 8c). Hence, they might be regarded as prognostic indicators, but also could be designed as therapeutic targets of AML.

Figure 8 Quantitative real-time PCR results of circRNAs and miRNAs expression. (a): RT-qPCR of 10 circRNAs (b): RT-qPCR of 6 miRNAs. (c): Relationships of circRNAs and miRNAs. **p<0.05.

Discussion

AML has been confirmed as a highly heterogeneous malignant hematologic tumor in cytogenetics and molecular biology, which has been precisely diagnosed by MIMC (ie, morphology, immunology, molecular biology, and cytogenetics). Molecular mechanism research of AML has facilitated the development of several novel molecular targeted therapies over the past few years,¹¹ whereas most AML patients will relapse, and the overall prognosis remains disappointing. CircRNAs are novel non-coding RNA family members with characteristic covalent closed loops. As high-throughput sequencing technology has been leaping forward, circRNAs, a novel class of RNA molecules, has been extensively investigated over the past few years. In AML, Li et al¹² analyzed circRNA microarray from 10 samples (four healthy controls and six newly diagnosed AML patients included) and suggested that hsa_circ_0004277 may be a promising diagnostic marker and treatment target in AML. Chen et al¹³ have suggested that circ-ANAPC7 up-regulated significantly in AML, will be a potential biomarker for AML. Circ-ANAPC7 may play a certain role in AML pathogenesis by serving as a sponge for the miR-181 family. The above two studies involved few samples and used a single chip, so the mentioned results may not be ideally validated in large sample data. With the rapid advance of gene sequencing technology, there have been increasing biological data, and some research has employed multiple microarrays for analysis. For instance, Zhang et al¹⁴ adopted GSE116617 and GSE94591 to obtain the microarray data of circRNA, employed GSE142699 to obtain the microarray data of miRNA, and used GSE114868 to obtain the microarray data of mRNA. They built a circRNA-miRNA-hub gene subnetwork and revealed four key circRNAs (including hsa_circ_0074371, hsa_circ_0029407, hsa_circ_0001394, and hsa_circ_0000994) involved in the pathogenesis of AML. With one of the GEO microarrays consistent with to ours, this study is a better addition to Zhang et al study. The reason for the better addition is that the rest of the microarray data were different, such that the results obtained by overlapping were different.

In the current study, we selected three microarray data (GSE94591, GSE142700-GPL26945, GSE142700-GPL19956) from GEO database, identified 10 DEcircRNA: hsa_circ_0017446, hsa_circ_0012152, hsa_circ_0005273, hsa_circ_0001947, hsa_circ_0071375, hsa_circ_0009581, hsa_circ_0077765, hsa_circ_0000497, hsa_circ_0058058, hsa_circ_0004277. Among the 10 identified DEcircRNAs, hsacirccirc0001947, Circcirc0058058 and hsacirc0004277 have been found to be correlated with AML.^15–18 Existing research has suggested that hsacirc0001947 inhibits AML cell proliferation via the hsa-miR-329-5p/CREBRF axis.¹⁵ Circ0058058 knockdown accelerated leukocytes apoptosis and inhibited AML progression by targeting miR-4319 to down-regulated EIF5A2 expression, and circ0058058 was considered a possible therapeutic target for AML.¹⁶ Li et al¹² has suggested that the level of hsacirc0004277 is down-regulated in newly diagnosed AML patients, while the expression of hsacirc0004277 is up-regulated in complete remission cases, and in relapsed/refractory cases, the level was again down-regulated. The above results suggested that hsacirc0004277 can be a diagnostic marker and therapeutic target in AML. Liu et al¹⁷ further suggested that circ0004277 hinders the development of AML by targeting miR-134-5p and up-regulating SSBP2. Of the remaining seven circRNAs. hsacirc0000497 can promote peritoneal metastasis of ovarian cancer by regulating epithelial-mesenchymal transition (EMT) in vitro.¹⁸ The remaining six circRNAs have not been found to be correlated with tumors.

To further identify and validate key circRNAs, a network of the correlations between circRNAs, miRNAs, and hub genes was established through Cytoscape. The network contains four circRNAs that have been mentioned above, four miRNAs and fourteen hub genes. Wang et al¹⁹ found that circSPI1 contributes to proliferation and apoptosis by interacting with hsa-miR-382-5p in AML. hsa-miR-150-5p²⁰ has been found to be associated with the proliferation and growth of a variety of solid tumors, such as non-small cell lung cancer,²¹ gastric cancer,²² and cervical cancer.²³ hsa-miR-454-3p²⁴ was associated with metastatic invasion of solid tumors such as breast and esophageal cancers.

We further performed RT-qPCR to confirm that some circRNA and miRNA are expressed differently in AML and healthy controls. The results showed that hsa_circ_0009581 and hsa_circ_0005273 were significantly upregulated in AML patients; hsa_circ_0000497 and hsa_circ_0001947 were significantly downregulated in AML patients. Hsa-miR-150-5p was significantly downregulated and hsa-miR-454-3p was upregulated in AML patients. hsa_circ_0009581 and hsa-miR-150-5p, hsa_circ_0001947 and hsa-miR-454-3p expression were significantly negatively correlated. According to the results of bioinformatics analysis and experiments, we speculated that hsa_circ_0009581 through hsa-miR-150-5p and hsa_circ_0001947 through hsa-miR-454-3p promoted leukemia development. hsa_circ_0001947 and hsa_circ_0009581 may provide new directions in the pathogenesis of AML.

Data from the TCGA database and the GTEx project were compared using GEPIA 2, and validation of the above hub genes (MYB, HSP90B1, PIK3R1, SMAD2, PKMYT1, RAD51 and TGFBR2) revealed differential expression in both AML patients and controls. Survival analysis revealed that high RAD51 expression was associated with poor prognosis. RAD51 whose core function is DNA repair is an essential gene for cell survival it can protect cells from genetic instability. In normal cells, it is regulated tightly and its expression is low while in most tumor cells, its expression level is high. RAD51 expression has been found to be significantly increased in breast and pancreatic cancers and is strongly associated with poor prognosis.^25,26 Seedhouse et al²⁷ investigated the role of RAD51 polymorphisms in modulating susceptibility to the development of AML. Park et al²⁸ suggested that hypermethylation of RAD51, MLH1, and ATM shows better overall survival than the non-hypermethylation group in a univariate analysis in AML patients. Inhibition of RAD51 expression enhances cytarabine and daunorubicin sensitivity in AML cells.²⁹ In our study, survival analysis of RAD51 revealed that high expression was associated with a poor prognosis, consistent with existing literature findings.

Four of the above hub genes were reported to be correlated with tumorigenesis in leukemia. Proto-oncogene MYB can be an important transcription factor to promote the development of AML. When its inhibitor monensin was used to effectively reduce MYB expression, it resulted in induction of cell differentiation and promotion of apoptosis.³⁰ Hyperactivation of the PI3K/AKT pathway is correlated with drug resistance and relapse in AML cells.³¹ Disease-free survival and overall survival were significantly reduced in AML patients with up-regulated PI3K/AKT signaling pathway.³² Glioma-related oncogene 1 (GLI1) can directly bind PI3KR1 and activate the PI3K/AKT/GSK3β/CDK signaling pathway, resulting in AML cell growth and reduced chemosensitivity.³³ Hypoxia can activate the TGF-β/Smad2 signaling pathway and induce multidrug resistance in leukemia cells.³⁴ Chen et al³⁵ have suggested that miR-425 is a prognostic factor in AML, and SMAD2 is one of the predicted targets of miR-425. Dou et al³⁶ also suggested that DARS-AS1 up-regulates TGFB1 expression by binding to miR-425 to regulate AML progression through the Smad2/3 pathway. Yu et al³⁷ have identified significant differences in miRNA data for ILCregs (regulatory innate lymphoid cells)-related genes (TGFBR2 included) in AML patients. These previous studies confirmed that our constructed ceRNAs are involved in the development of AML and may be therapeutic targets.

This study suggests that DEcircRNAs take on a critical significance in AML development and may be the effective therapeutic targets. The data from this study provide a novel perspective for pathogenesis of AML. Furthermore, our further research will examine the exact role of the above circRNAs in AML development and verify the hub genes identified in this study.

Conclusions

In this study, differentially expressed circRNAs, miRNAs and mRNAs were selected from the published microarray data of AML by bioinformatics analysis, and CircRNA/miRNA/hub gene networks were built, which confirmed by experiment and provided novel ideas for genetic research of AML. Despite rigorous bioinformatics analysis of this study, due to the limitations of sample number, Further clinical experiments are required to confirm the key circRNAs and seven hub genes including RAD51.