Skip to main content

ORIGINAL RESEARCH article

Front. Oncol., 09 April 2021
Sec. Molecular and Cellular Oncology
This article is part of the Research Topic Therapeutic Opportunities and Innovative Biomarkers in Tumor Microenvironment View all 23 articles

Tumor Microenvironmental Competitive Endogenous RNA Network and Immune Cells Act as Robust Prognostic Predictor of Acute Myeloid Leukemia

Yaqi Cheng&#x;Yaqi Cheng1†Xiaoran Wang*&#x;Xiaoran Wang1*†Peiyan QiPeiyan Qi2Chengxiu LiuChengxiu Liu3Shoubi WangShoubi Wang1Qi WanQi Wan1Yurun LiuYurun Liu1Yaru SuYaru Su1Lin JinLin Jin1Ying LiuYing Liu1Chaoyang LiChaoyang Li1Xuan SangXuan Sang1Liu YangLiu Yang1Chang LiuChang Liu1Hucheng DuanHucheng Duan4Zhichong Wang*Zhichong Wang1*
  • 1State Key Laboratory of Ophthalmology, Zhongshan Ophthalmic Center, Sun Yat-sen University, Guangzhou, China
  • 2Guangzhou International Travel Health Care Center, Guangzhou, China
  • 3Department of Ophthalmology, Affiliated Hospital of Qingdao University Medical College, Qingdao, China
  • 4Department of Ophthalmology, The Second People’s Hospital of Foshan, Foshan, China

Acute myeloid leukemia (AML) is malignant hematologic tumors with frequent recurrence and cause high mortality. Its fate is determined by abnormal intracellular competitive endogenous RNA (ceRNA) network and extracellular tumor microenvironment (TME). This study aims to build a ceRNA network related to AML TME to explore new prognostic and therapeutic targets. The RNA expression data of AML were obtained from The Cancer Genome Atlas (TCGA) database. First, we used the ESTIMATE algorithm to calculate the immune cells and stromal cells infiltration scores in the TME and found that all scores were highly correlated with AML’s prognostic characteristics. Subsequently, differentially expressed mRNAs and lncRNAs between high and low score groups were identified to construct a TME-related ceRNA network. Further, the Cox-lasso survival model was employed to screen out the hub prognostic ceRNA network composed of two mRNAs (EPB41L3, COL2A1), three miRNAs (hsa-mir-26a-5p, hsa-mir-148b-3p, hsa-mir-148a-3p), and two lncRNAs (CYP1B1-AS1, C9orf106), and construct nomograms. Finally, we used CIBERSORT algorithm and Kaplan-Meier survival analysis to identify the prognostic TME immune cells and found that naive B cells, M2-type macrophages, and helper follicular T cells were related to prognosis, and the hub ceRNAs were highly correlated with immune cell infiltration. This study provided a new perspective to elucidate how TME regulates AML process and put forward the new therapy strategies combining targeting tumor cells with disintegrating TME.

Introduction

Acute myeloid leukemia (AML) is one of the most malignant hematologic tumors, characterized by the massive expansion of abnormally differentiated hematopoietic precursor cells (1, 2). AML shows a high degree of heterogeneity due to its complex genetic mutations and variable molecular phenotypes, which results in the failure of traditional treatment with a single mode of action in achieving an ideal effect (3). Although 50% of the AML patients could obtain a complete remission (CR) after receiving therapies, such as intensive induction chemotherapy and post-remission, still more than 20% of AML cases remain unresponsive and refractory. Even among patients who achieve CR, there is still a recurrence rate of up to 50% (4). Hematopoietic stem cell transplantation, the only method to cure AML currently, still faces severe challenges, such as tumor immune escape, recurrence, and graft-versus-host response (5). Targeted chemotherapy and chimeric antigen receptor cell therapy have brought a new breakthrough to the treatment of AML. But off-target, drug resistance, and even treatment-related death, still pose insurmountable obstacles to its safe application (6). Old age, high white blood cell counts, abnormal genotypes of leukemia cells, and complicating other myeloid diseases have significant effect on the prognosis of AML. To promote the transformation of prognostic markers, exploration into clinical application is the driving force for the leap-forward development of all tumor, including AML treatment strategies.

The generation and progression of tumor are determined by abnormal molecular aberrations inside tumor cells and the tumor microenvironment (TME) outside the cells. AML originates from abnormal bone marrow (BM) microenvironment. Tumor cells hijack and reshape the BM microenvironment ecology to transform it into tumor protective phenotype, mediating the immune escape and therapeutic tolerance of AML (7). Under the protection of abnormal BM microenvironment, AML stem cells and initiation cells can maintain their potential of regeneration, and small residual lesions can obtain effective incubation and induce the recurrence of AML (8). As a hotbed of AML, the abnormal BM microenvironment employs extracellular matrix (ECM) as a scaffold and contains cellular components, such as immune cells, stromal cells, endothelial cells, and various soluble molecules, such as exosomes, cytokines, and hormones. Among them, stromal cells and immune cells are key components that affect the progression of AML and therapeutic response.

The emergence of competitive endogenous RNA (ceRNA) network theory provides new ideas for exploring the internal mechanisms of tumorigenesis and development. Long non-coding RNA (lncRNA), circular RNA (circRNA), and other non-coding RNA (ncRNA) form an endogenous RNA competitive regulatory network when competing with mRNA for the binding to microRNA (miRNA). The competitive effect of ncRNAs is crucial to the structural stability and translation function of mRNAs. Abnormal ceRNA networks are involved in various tumor processes, including AML (9, 10). A number of studies have clarified the role of lncRNA and other ncRNAs or TME in AML and revealed their potential value in predicting the prognosis of AML (11, 12). Conducting in-depth exploration of the internal and external mechanisms of AML progress, looking for highly efficient prognostic targets, and combining precise targeting tumor cells with completely disintegrating tumor protective microenvironment may bring new breakthroughs for the treatment of AML.

In this study, we obtain the clinical information and transcriptome expression data of AML patients from The Cancer Genome Atlas (TCGA) database. The Estimate of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) was employed to quantify the scores of immune and stromal cells infiltrations in TME. DEmRNAs and DElncRNAs between high and low scores groups were identified, and the TME-related ceRNA network was constructed. We screened hub prognostic genes and established nomograms to quantify their predictive power. Meanwhile, The Cell Type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) algorithm was used to calculate the composition and proportion of immune cells in TME, and the correlation between ceRNAs and immune cell infiltration was verified, so as to provide more reliable biological targets for AML therapy. Finally, we identified two lncRNAs (CYP1B1-AS1, C9orf106), three miRNAs (hsa-mir-26a-5p, hsa-mir-148b-3p, hsa-mir-148a-3p), and two mRNAs (EPB41L3 and COL2A1) to construct a prognostic lncRNA-miRNA-mRNA ceRNA network. Furthermore, our study identified that several immune cells, such as M2 macrophages, naive B cells, and helper follicular T cell, have dramatic prognosis value.

Materials And Methods

AML Transcriptome Data and Clinical Information

The RNA sequencing data and corresponding clinical information of 173 AML samples were obtained from the TCGA database (https://portal.gdc.cancer.gov/). The miRNA expression array GSE142699 (GPL26945 NanoString nCounter Human miRNA) downloaded from the Gene Expression Omnibus (GEO) database contains 24 AML and 24 normal control samples, and mRNA expression array GSE71014 (GPL10558 Illumina humanht-12 V4.0 expression beadchip) contains mRNA data and clinical information of 104 AML samples.

Analysis of Differentially Expressed Genes Related to Microenvironment

The ESTIMATE algorithm can evaluate the non-tumor cell components in TME based on the gene expression characteristics of the tumor, and quantify the TME immune cell and stromal cell infiltration scores (13). ESTIMATE algorithm was used to analyze the AML mRNA sequencing data and calculate the BM microenvironment scores. The AML samples were divided into high, low immune infiltration groups and stromal infiltration groups with the median score as the boundary. The “limma” R package was used to identify differentially expressed mRNAs (DEmRNAs) and lncRNAs (DElncRNAs) between high- and low-stromal or immune score groups with |log2FC| > 1 and adjusted P value < 0.05 as the cutoff criteria, and then intersection was taken (14, 15). Meanwhile, differentially expressed miRNAs (DEmiRNAs) between AML and normal samples were identified after normalization, |log2FC|>1 and adjusted P value <0.05 was defined as the cutoff criteria (16, 17). The “pheatmap” and “ggplot2” R packages were involved in drawing heatmap and volcanoes of the differential genes.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Analysis

The “clusterProfiler” (18) R package was used to perform gene ontology (GO) enrichment analysis of DEmRNAs to reveal the biological processes (BP), cellular components (CC) and molecular functions (MF) that they are involved in. Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis was used to annotate the signaling pathways DEmRNAs involved. The adjusted P value <0.05 was statistically significant.

CeRNA Network Construction

MiRcode (19) database was employed to predict the miRNAs (lnc-pre-miRNAs) targeted by DElncRNAs. The intersection of lnc-pre-miRNAs and DEmiRNAs was obtained, and then the mRNAs (mi-pre-mRNAs) targeted by the intersecting miRNAs were predicted through Targetscan (20)and miRDB (21) databases. The intersection was taken within mi-pre-mRNAs and DEmRNAs. miRNAs and lncRNAs related to common mRNAs were identified according to the targeting relationship, and Cytoscape v3.7.2 was used to construct the initial TME-related ceRNA network.

Cox, Lasso Regression Analysis, and Construction of Nomogram

The Cox-lasso survival model was constructed to screen the hub prognostic mRNAs. Univariate Cox proportional hazards regression analysis was employed to identify the relationship between the mRNAs expression and overall survival (OS) of patients using the “survival” R package, and the forest map was drawn using the “forestplot” R package. The mRNAs with P<0.05 entered lasso regression and multivariate Cox proportional hazards regression, which were also performed with the “survival” R package. Subsequently, with the median of the risk scores obtained by multivariate Cox regression as the boundary, the AML patients were divided into high and low risk groups. The “survival” R package was employed to perform survival analysis and draw survival curves of the two groups. At the same time, the hub prognostic mRNAs nomogram was established based on the results of multivariate Cox regression to quantitatively predict the prognosis of AML, using the “rms” R package.

TME Immune Cell Infiltration Analysis

CIBERSORT algorithm analyzes gene expression data to identify cell abundance and proportion in mixed cell populations. CIBERSORT can easily recognize each cell type and its count in each sample (22). We used CIBERSORT algorithm to identify the proportion of the 22 types of immune cells infiltrated in AML samples. Analyses were performed with 100 permutations with the default statistical parameters to improve the credibility of the results. Samples with a CIBERSORT output of P < 0.05 were considered to be eligible for further analysis (16, 23).

Statistical Analysis

The Kaplan‐Meier survival analysis was performed using the R package “Survival” to analyze the correlations between patients’ OS and other variates, such as TME scores, hub prognostic mRNAs’ expression, and immune cell infiltration proportion, and also employed in identifying the OS’ difference between high and low risk groups after Cox regression, two-sided P <0.05 was the statistically significant cutoff. The statistical significance of the correlation was tested by the log-rank test. Pearson correlation analysis was done for each prognostic biomarker in the ceRNA network and the proportion of each microenvironment related immune cell using the “limma” package (23). All statistical analyses were performed by R software (v.3.6.3) and corresponding program packages.

Results

TME Scores Were Highly Related to the AML Patients Prognosis

The results of ESTIMATE analysis showed that the immune scores of 173 AML samples ranged from 1329.53 to 3971.97, and the stromal scores ranged from −1888.81 to 435.75. The ESTIMATE scores were the sum of immune scores and stromal scores. The AML patients were divided into high and low groups according to stromal scores, immune scores, and ESTIMATE scores to perform survival analysis. Results showed that the survival time of the patients with high immune scores (P = 0.021; Figure 1A) and high ESTIMATE scores (P = 0.011; Figure 1C) is much shorter than those with the low scores. Meanwhile, the overall survival time of patients with high stromal scores was shortened, but there was no significant statistical difference (P = 0.69; Figure 1B).

FIGURE 1
www.frontiersin.org

Figure 1 Identification of the relationship between TME scores and prognosis of AML. (A) Survival curve of high and low immune groups: the survival rate of patients with high immune scores was significantly reduced (P = 0.021). (B) Survival curves of high and low stromal groups: the survival rate of patients with high stromal score decreased, with no significant statistical difference (P = 0.69). (C) Survival analysis showed that the survival rate of patients with high ESTIMATE score was significantly reduced (P = 0.011). (D–F) The distribution of the three scores in different FAB classifications are significantly different (all P values < 0.0001); (G–H) The relationship between the survival status of AML patients and three scores: the immune scores of the dead cases (P = 0.0093) and the ESTIMATE scores (P = 0.029) were significantly higher than those of the surviving patients, and the patients with high stromal scores were more likely to be in the state of death, but there was no significant statistical difference (P = 0.2).

At the same time, we analyzed the relationship between three types of scores and prognosis-related clinical traits. The results showed that immune, stromal, and ESTIMATE scores were all significantly different among different FAB classification of leukemia (P < 0.001; Figures 1D–F); The three scores were also correlated with the survival status of the patients. Each score in the death cases was significantly increased, among which stromal score were not statistically significant (Figures 1G–I). In addition, immune scores were significantly different in different cytogenetic risk classifications (P = 0.011; Supplementary Figures 1A–C), while the three scores were not related to gender (P >0.05; Supplementary Figures 1D–G). All the above results indicate that TME immune cells and stromal cells are of great significance in the prognosis of AML, especially in the diagnosis and classification of AML.

Identification of TME-Related DEmRNAs and DElncRNAs

In order to evaluate the possible impact of stromal and immune scores on breast cancer, we investigated the expression patterns in different stromal and immune groups. DEmRNAs and DElncRNAs were compared between high- and low-stromal or immune score groups with |log2FC| > 1 and adjusted P value < 0.05 as the cutoff criteria through R package “limma.” For comparison based on immune scores, there were 414 DEmRNAs (Figures 2A, C) and 294 DElncRNAs (Figures 3A, B) in the high immune-score group. Similarly, for the high and low groups based on stromal scores, 367 DEmRNAs (Figures 2B, D) and 222 DElncRNAs (Figures 3C, D) were obtained. After taking the intersection of mRNAs and lncRNAs separately, we got 285 mRNAs (Figure 2E) and 172 lncRNAs (Figure 3E) that were differentially expressed in in both stromal and immune groups. We believed that these genes were highly correlated with TME and used them for further analysis.

FIGURE 2
www.frontiersin.org

Figure 2 Identification of TME-related mRNAs. (A) Heatmap of DEmRNAs between high and low immune groups. (B) Heatmap of DEmRNAs between high and low stromal groups. (C) Volcano map of mRNAs expression between high and low immune groups. (D) Volcano map of mRNAs expression between high and low stromal groups. (E) Venn diagram was used for the intersection of DEmRNAs within immune and stromal groups, we considered the common DEmRNAs to be TME-related mRNAs. (F) GO enrichment analysis was performed to annotate TME-related DEmRNAs: the larger the bubble and longer columns represent the more genes enriched in this function, the deeper the color of the bubble and bars, the smaller the P value. (G) TME-related DEmRNAs KEGG enrichment analysis results.

FIGURE 3
www.frontiersin.org

Figure 3 Identification of TME-related lncRNAs. (A–D) Heatmaps and volcano maps of DElncRNAs between high and low immune and stromal groups. (E) Venn diagram was used for the intersection of DElncRNAs within immune and stromal groups, and we considered the common DElncRNAs to be TME-related mRNAs. (F) miRNAs predicted by TME-related lncRNAs.

285 common DEmRNAs were subjected to perform GO and KEGG enrichment analysis to further verify the biological processes involved. The results of GO analysis showed that they were mainly enriched in cytokine production and secretion (GO: 0050663, GO: 0001819, GO: 0050707), neutrophil activation (GO: 0042119, GO: 0043312), lymphocytes (GO: 0050900), and T Cell (GO: 0042110) function and immune response (GO: 0050727, GO: 0045088) and other immune processes (Figure 2F, Table 1). KEGG analysis showed they were enriched in signaling pathways, such as osteoclast differentiation (hsa04380), cytokines and cytokine receptor interaction (hsa04060), phagosome (hsa04145), and B cell receptor (hsa04662) (Figure 2G, Table 2). Both are related to tumor immune response and TME remodeling.

TABLE 1
www.frontiersin.org

Table 1 GO enrichment analysis results of the TME-related DEmRNAs.

TABLE 2
www.frontiersin.org

Table 2 KEGG enrichment analysis results of the TME-related DEmRNAs.

Identification of DEmiRNAs and Preliminary Construction of TME-Related ceRNA Network

We obtained the DEmiRNAs between the 24 AML and 24 normal control samples in GSE142699 miRNA microarray. A total of 61 DEmiRNAs were identified with |log2FC|>1 and adjusted P value <0.05 as the cutoff criteria (Figure 4A).

FIGURE 4
www.frontiersin.org

Figure 4 Construction of TME-related ceRNA network. (A) Heatmap of DEmiRNAs between AML and normal control samples in GSE142699 microarray. (B) Venn diagram of the intersecting miRNAs between DEmiRNAs and lnc-pre-miRNAs. (C) Venn diagram of the intersecting of mRNAs between mi-pe-mRNAs and TME-related DEmRNAs, a total of 43 common mRNAs. (D) Preliminary establishment of TME-related ceRNA network, including 43 mRNAs, 23 miRNAs, and 7 lncRNAs.

We predicted the targeting relationships among lncRNAs, miRNAs and mRNAs through a variety of databases and preliminarily established TME-related ceRNA network.

First, 55 lnc-pre-miRNAs were obtained by predicting the 172 common DelncRNAs targeting genes (Figure 3F). The intersection of lnc-pre-miRNAs and DEmiRNAs yielded 23 miRNAs (Figure 4B). Subsequently, 4187 mi-pre-mRNAs targeted by the 23 intersecting miRNAs were predicted by both Targetscan and miRDB databases (Figure 4C). In order to improve the correlation between ceRNA network and TME, the intersection of mi-pre-mRNAs and DEmRNAs was taken, and 43 hub mRNAs were obtained. Finally, we preliminarily constructed a TME-related ceRNA network consisting of 43 mRNAs, 23 miRNAs and 7 lncRNAs (Figure 4D). GO and KEGG enrichment analysis were performed on 43 mRNAs to verify their involvement in TME and immune-related processes (Supplementary Figure 2).

Survival Model to Screen Hub Prognostic mRNAs and Construct Nomograms

In order to ensure the accuracy and sensitivity of AML prognosis prediction, we constructed a survival model consisting of univariate Cox proportional hazard regression-lasso regression-multivariate Cox proportional hazard regression to further screen the hub prognostic ceRNAs. Firstly, 16 of 43 common mRNAs which were highly related to overall survival (OS) of AML patients were screened by univariate Cox regression (P <0.05; Figure 5A). To prevent overfitting of multivariate Cox regression, Lasso regression was used to screen eight mRNAs from 16 mRNAs (Figures 5B, C). Three hub prognostic mRNAs: KCNK10, EPB41L3, and COL2A1 were finally screened out by multivariate Cox regression (Figure 5D, Table 3). The receiver operating characteristic (ROC) curve was drawn to check the accuracy of the model. The area under the curve (AUC) of 5-year survival was 0.777, indicating the high accuracy of this model (Figure 5E). Kaplan-Meier survival analysis of high and low risk groups showed that the survival rate of the high-risk group was significantly lower than that of the low-risk group (P = 0.001; Figures 5F, G). Simultaneously, heatmap of the differential expression of three hub mRNAs between the high and low risk groups were drawn (Figure 5H), and a nomogram was constructed according to the expression of the three mRNAs (Figure 5I).

FIGURE 5
www.frontiersin.org

Figure 5 Construction of the survival model to screen hub prognostic mRNAs. (A) Forest map of univariate Cox proportional hazards regression: 16 mRNAs with P < 0.05 were screened out of 43 common mRNAs (marked by *). (B, C) Lasso regression screening of mRNAs into multivariate Cox regression: λ and cv graphs show that eight hub mRNAs were selected. (D) Forest map of multivariate Cox proportional hazards regression: three hub prognostic mRNAs (KCNK10, EPB41L3, COL2A1) were screened out. (E) ROC curve of multivariate Cox regression: AUC of 5-year survival was 0.777, indicating the high accuracy of multivariate Cox regression model. (F) The AML samples were divided into high and low risk groups according to the risk scores obtained from Cox regression. (G) Survival curve of high and low risk groups: the survival rate of the high-risk group was significantly reduced (P = 0.0011). (H) Heatmap of KCNK10, EPB41L3, and COL2A1 expressions among the high and low risk group samples. (I) Nomogram of KCNK10, EPB41L3, and COL2A1 based on the results of the multivariate Cox regression.

TABLE 3
www.frontiersin.org

Table 3 Identification of hub prognostic mRNAs through univariate Cox regression-lasso regression-multivariate Cox regression.

At the same time, we performed Kaplan-Meier survival analysis of the three hub mRNAs separately. The results showed that the survival rate of patients with high EPB41L3 expression was significantly reduced (P = 0.011; Figure 6A), and low expression of COL2A1 predicted poor prognosis (P = 0.013; Figure 6B), while there was no significant relationship between the KCNK10 expression and the prognosis of AML (P = 0.681; Figure 6C). The survival analysis results of EPB41L3, COL2A1, and KCNK10 in GSE71014 array were consistent with the above, but there was no significant statistical difference in the survival analysis of the three mRNAs (Supplementary Figure 3). Therefore, we identified EPB41L3 and COL2A1 as the hub prognostic mRNAs and plotted the prognostic nomogram of the two genes (Figure 6D). Finally, the hub ceRNA network consisting of two mRNAs (EPB41L3 and COL2A1), three miRNAs (hsa-mir-26a-5p, hsa-mir-148b-3p, hsa-mir-148a-3p), and two lncRNAs (CYP1B1-AS1, C9orf106) were constructed (Figure 6E).

FIGURE 6
www.frontiersin.org

Figure 6 Verifying the prognostic value of the three core mRNAs. (A–C) Survival curve of the three hub mRNAs: the survival rate of AML patients with high expression of EPB41L3 significantly decreased (P = 0.011), and the survival rate of AML patients with low expression of COL2A1 significantly decreased (P = 0.013). There is no significant difference in survival rate between different KCNK10 expression levels (P = 0.681), so we believe that EPB41L3 and COL2A1 are core prognostic mRNAs with predictive value. (D) Nomogram of EPB41L3 and COL2A1 based on the results of the multivariate Cox regression. (E) Hub prognostic ceRNA network based on EPB41L3 and COL2A1; (F) Forest map of multivariate Cox regression analysis: the risk score of prognostic model maintained independence in predicting the OS of AML patients (HR = 1.8469, 95% CI = 1.2652–2.6961, P = 0.00147). At the same time, age and cytogenetic risk groups of AML patients also have independent prognosis value. (G) The nomogram of the risk score of prognostic model and clinic characters.

Furthermore, in order to evaluate the independent predictive value of the two hub mRNAs prognostic model in patients with complete clinical information from the TCGA cohort, we constructed univariate and multivariate Cox regression analyses according to the risk scores from prognostic model and clinical covariates including age, gender, cytogenetic risk group, FAB classification, history of hydroxyurea treatment. The result of univariate Cox regression analysis revealed that risk scores have obvious predictive value (P < 0.0001), and some clinic characters, such as age (P < 0.0001) and cytogenetic risk group (P < 0.001) also had prognosis significance (Table 4). The multivariate Cox regression analysis indicated that the risk score of prognostic model maintained independence in predicting the OS of AML patients (HR = 1.8469, 95% CI = 1.2652–2.6961, P = 0.00147) (Figure 6F, Table 4). We plotted the prognostic nomogram of the risk score of prognostic model and clinic characters (Figure 6G).

TABLE 4
www.frontiersin.org

Table 4 Identification the independent prognostic value of survival model and clinic characters through univariate Cox regression and multivariate Cox regression.

AML Microenvironment Immune Cell Infiltration and Its Prognostic Value

In this study, the CIBERSORT algorithm was used to analyze the 22 immune cell subtypes infiltration in the AML microenvironment (Figure 7A) and the infiltrating differences between the high and low immune and stromal score groups obtained by the ESTIMATE algorithm were analyzed (Supplementary Figure 4). The results showed that the resting dendritic cells, resting mast cells and neutrophils were significantly different between different immune and stromal groups (P <0.05), and the activated mast cells (P = 0.013) were significantly different between the high and low stromal groups (P = 0.013). Then Pearson test was performed to identify the correlation between TME immune cells and the correlation heatmap was drawn (Figure 7B). The results showed that there were strong correlations among various immune cells such as M2 macrophages, activated NK cells, memory B cells, and helper follicular T cells (Figure 7C, Supplementary Figure 5).

FIGURE 7
www.frontiersin.org

Figure 7 Correlation analysis of TME immune cells and prognosis. (A) Histogram of the infiltration proportion of 22 types of immune cells between different samples. (B) Correlation heat map between immune cells. (C) Plots of correlation between immune cells: memory B cells and monocytes (P = 0.013) was significantly positively correlated; there was a significant negative correlation between resting mast cells and M2 macrophages (P = 0.007), activated NK cells and M0 macrophages (P = 0.044), naive CD4 memory T cells and helper follicular T cells (P = 0.026). (D) Survival curves of TME immune cell: the survival rate of AML patients with high proportion of M2 macrophages infiltration was significantly reduced (P = 0.002), patients with low proportion of naive B cells (P <0.001) and helper follicular T cells (P = 0.016) had decreased survival rate. In addition, the increase in the ratio of γδ T cells, activated CD4 memory T cells and memory B cells also predicted the adverse outcome of AML, with no statistically significant difference.

The results of Kaplan-Meier survival analysis showed that the M2 macrophages (P = 0.002), naive B cells (P <0.001), helper follicular T cell (P = 0.016) were significantly correlated with the survival rate of patients. In addition, the increased proportion of γδT cells, activated CD4 memory T cells and memory B cells also predicted a poor outcome, but the results were not statistically significant (Figure 7D, Supplementary Figure 6).

Correlation Between TME-Related Hub Prognosis ceRNA Network and Immune Cell Infiltration

In order to fully elucidate the mechanism of the hub prognostic ceRNA network in AML, and at the same time look for the potential way for TME to determine the fate of AML, we performed Pearson test on the correlation between EPB41L3, COL2A1, CYP1B1-AS1, C9orf106 and TME immune cells and plotted heatmap (Figure 8A). The results showed that EPB41L3, COL2A1, CYP1B1-AS1, and C9orf106 were all related to TME immune cells. EPB41L3 was significantly positively correlated with resting mast cells, neutrophils and M1 macrophages, while being negatively correlated with plasma cells; COL2A1 was positively correlated with naive CD4 T cells; CYP1B1-AS1 was positively correlated with resting CD4 memory T cells, and C9orf106 was positively correlated with activated dendritic cells (all P values < 0.05; Figure 8B). In addition, we found that the expression of CYP1B1-AS1 was significantly positively correlated with that of the EPB41L3 and C9orf106 (Supplementary Figure 7).

FIGURE 8
www.frontiersin.org

Figure 8 Correlation analysis of hub prognostic ceRNA network and TME immune cells. (A) Correlation heat map between EPB41L3, COL2A1, CYP1B1-AS1, C9orf106, and TME infiltrating immune cells; (B) Verifying the correlation between EPB41L3, COL2A1, CYP1B1-AS1, C9orf106, and immune cells: EPB41L3 was significantly positively correlated with resting mast cells, neutrophils, and M1 macrophages, and had strong negative correlation with plasma cells, COL2A1 was positively correlated with naive CD4 T cells. CYP1B1-AS1 was positively correlated with resting CD4 memory T cells. C9orf106 was positively correlated with activated dendritic cells (all P values <0.05).

Discussion

AML progresses rapidly and causes a very high mortality rate. Its initiation, development, drug resistance and recurrence all depend on the abnormal molecular and genetic changes inside the cell and the protection of the extracellular TME (24). More and more scholars have realized that current immune therapies and targeted therapies cannot fully deal with the extremely complex heterogeneity of AML, and are committed to excavating new prognostic biomarkers to promote the innovation of AML precision diagnosis and treatment (25, 26). This study identified the TME related hub prognostic lncRNA-miRNA-mRNA ceRNA network, while focusing on the correlations between specific TME immune cell populations and ceRNA network of AML. A multiple prognostic verification model based on TME score survival analysis-ceRNA survival model-TME immune cell infiltration survival analysis was constructed, which highly improved the prognosis predicting accuracy of the biomarkers in this study while ensuring the high correlation between ceRNA and TME. We tried to elucidate the mechanism of TME regulating the fate of AML, and provide a unique perspective for finding new targets for AML.

AML cells manipulate TME through a complex interaction network, domesticate and reshape it into a pro-leukemia phenotype. The modified AML protective TME in turn promotes AML progression, providing tumor cells with an increasingly strong bastion, and forming a positive feedback loop of tumor-TME mutually promoting. Stromal cells and immune cells in TME have important regulatory and protective effects on AML. Abnormally expressed adhesion molecules, cell cycle regulators and angiogenic factors in tumor-associated stromal cells promotes the angiogenesis in TME, and meanwhile enable AML cells to accelerate proliferation, resist apoptosis and malignantly invade (27). TME immune cell-mediated inflammatory response is considered to be an important driving force for the remodeling of the AML microenvironment (28). In this study, ESTIMATE algorithm was employed to identify the infiltrating immune and stromal cells in the AML TME, and it was confirmed that the infiltration degree of these two types of cells was significantly correlated with the prognosis of AML, and it had a prominent effect on assisting the diagnosis of AML FAB typing. At the same time, the TME related mRNAs and lncRNAs were obtained according to the TME score grouping and the ceRNA network was constructed. These results and opinions are consistent with the findings of Yan (11), Huang (29) and others who utilized the ESTIMATE algorithm to identify hub prognostic mRNAs related to the AML microenvironment.

As an important post-transcriptional regulatory mechanism, ceRNA network is of great significance in the progression of AML. Several studies have proved that HOXA-AS2, RPPH1, and other lncRNAs regulate mRNAs expression as ceRNA and participate in AML proliferation, differentiation, and invasion (10, 30, 31). We constructed a survival prediction model of univariate Cox regression-lasso regression-multivariate Cox regression, supplemented by Kaplan-Meier survival analysis and nomograms, to identify hub mRNAs (EPB41L3, COL2A1) with both high TME correlation and high prognostic efficacy and construct the hub prognostic ceRNA network by layer-upon-layer screening. As early as 2003, Celal et al. reported the high expression of EPB41L3 in AML cell line HL-60 (32), which was consistent with the poor prognosis of AML patients with up-regulated EPB41L3 expression in this study. As the gene involved in cytoskeletal construction, EPB41L3 has also been shown to promote tumor metastasis by promoting epithelium-mesenchymal transformation in advanced lung cancer (33). The dynamic balance of osteoblasts and osteoclasts is of great significance for maintaining the normal BM microenvironment. Increased osteoclast activity in the BM of AML patients results in bone demineralization and the destruction of the normal BM structure (34). COL2A1 encodes type II collagen α1 chain, which participates in the composition of ECM, is an important osteogenic protein. Decreased expression of COL2A1 reduces the osteogenic protein, destroys normal ECM structure, promotes the remodeling of the BM microenvironment toward the direction of tumor promotion, and participates in the mutual promotion between TME and AML (35, 36). Ganapathi et al. demonstrated that the low expression of COL2A1 was significantly associated with the rapid recurrence of high-grade serous breast cancer, and proposed that the tumor suppressive effect of COL2A1 might be achieved by depleting oncogene miR-301 as ceRNA (37), revealing the important role of COL2A1 related ceRNA network in tumor development. Meanwhile, miRNAs in this study have been proved to be closely related to AML. The decreased expression of hsa-mir-26a-5p in AML can cause high expression of peroxiredoxin III, thereby promptly clearing the reactive oxygen species within cells and protecting tumor cells from oxidative stress injury (38, 39); In another study, researchers found that enforced expression of hsa-mir-26a-5p in AML cells was able to inhibit cell cycle progression by downregulating cyclin E2 expression, potentiated the antiproliferative effects of 1,25-dihydroxyvitamin D (3) (VitD) and stimulated myeloid differentiation by targeting E2F7 (40). Furthermore, hsa-mir-26a-5p was also proved as a target of c-Myc, which revealed the vital role this miRNA played in ceRNA network (41). Huang et al. demonstrated that MLL-fusion/MYC⊣miR-26a⊣TET1 signaling circuit played an important role in AML, in which hsa-mir-26a-5p functioned as an essential tumor-suppressor mediator and its transcriptional repression was required for the overexpression and oncogenic function of TET1 in MLL-rearranged AML (42). miRNAs sequencing of exosomes derived from bone marrow mesenchymal stromal cells and bone marrow specimens also found hsa-mir-26a-5p to be significantly associated with overall survival of AML patients and was closely related to HOX-related genes (43). Wang et al. found that the expression of hsa-mir-148 family decreased in AML patients, and it was highly correlated with FAB typing of AML (44).They also verified that DNMT1 was identified to be a downstream target of hsa-mir-148, and was negatively regulated by miR-148a in AML cell lines. It was a strong evidence that hsa-mir-148 was involved in the regulation of AML as a core member of ceRNA network (45). However, the role of two hub lncRNAs (CYP1B1-AS1, C9orf106) in tumors has not been reported. In our previous study, we have established an AML related circRNA-lncRNA-miRNA-mRNA ceRNA network based on the differentially expressed RNAs and the target relationship among circRNA-miRNA, lncRNA-miRNA, and miRNA-mRNA. Through comparing the previous results and the new discovery in this study, we noticed that lncCYP1B1-AS1 were involved in both studies, indicating that lncCYP1B1-AS1 may perform an important function in regulating the intracellular biological process and extracellular microenvironment of AML cells, and it had great potential and research value as a powerful prognostic marker (46). Zhang at al. identified 10 RNAs (LINC00471, hsa-mir-100, hsa-mir-150, ANP32E, ERMP1, MYO1B, PAPD7, PTGIS, TERF1, and VEGFA) to be ceRNAs closely related to childhood AML. Even prognostic RNAs in Zhang’s study are different from those in our study. We speculate that it may be due to the childhood AML remarkably differing from adult AML in karyotype, therapeutic strategy, and therapeutic effects. Still and all, it reminds us to consider more general markers for both types of AML and excavate more accurate markers for various subtypes (47).The abovementioned confirms that the ceRNA network established in this study plays an important role in the remodeling of the AML microenvironment, and has high prognostic value.

Chimeric antigen receptor immune cell therapy brings new hope for the treatment of leukemia and myeloma, and confirms that the tumor treatment using normal immune cells to penetrate and change the abnormal TME immune cell populations is effective and feasible. To systematically reveal the immune microenvironment of AML, we also used CIBERSORT algorithm to identify TME infiltrating immune cells, verify their prognostic value, and explore the relationship between TME immune cells and hub prognosis ceRNA. The results showed that the higher the proportion of M2 macrophages infiltration was, the poorer the prognosis (P = 0.002), and the lower the proportion of naive B cells (P <0.001) and helper follicular T cells (P = 0.016) would be, resulting in the lower the survival rate of the patients. In addition, the increased proportion of γδT cells, activated CD4 memory T cells and memory B cells also predicted the poor outcome of AML. Macrophages can be divided into M1 type that suppresses tumors and M2 type that promotes tumors. TME can induce the differentiation of monocytes and mesenchymal stem cells into M2 macrophages and promote the malignant process of tumor. M2 macrophages can not only suppress the anti-tumor immune response, but also promote tumor proliferation, invasion, migration and angiogenesis (48). Through analyzing a variety of blood tumors by CIBERSORT algorithm, Xu et al. found that the proportion of M2 macrophages in the AML microenvironment was much higher than that of the normal control and even other tumors, and it also predicted the reduced survival rate of patients and the rapid recurrence (49). In another study which identifies infiltrating lymphocytes in TME by CIBERSORT algorithm, γδT cells were found significantly increased in various hematological tumors, such as M3-AML, chronic myeloid leukemia, B-cell acute lymphoblastic leukemia, and could promote tumor progression by promoting tumor-related inflammatory response and mediating the formation of immunosuppressive TME (50). It is consistent with the results that proportions of TME infiltrating of M2 macrophages and γδT cells determines the outcome of AML in our study. Besides, we found that EPB41L3, COL2A1, CYP1B1-AS1, and C9orf106 are all highly correlated with immune infiltration, which is consistent with various studies. For example, Wendell et al. found that COL2A1 was bound up with cytotoxic lymphocyte immune signature and T-cell trafficking (51). It coincides with the finding of this study that there is a positive correlation between COL2A1 and naive CD4 T cells, confirming that COL2A1 is involved in the regulation of AML immune microenvironment.

Nevertheless, there are still some inescapable deficiencies in this study. On the one hand, due to the limitation of the existing data amount, the comparative study of different FAB classification and cytogenetic risk groups could not be completed, suggesting that more expression data is urgently needed for more accurate analysis. On the other hand, although we correlated ceRNA network with AML microenvironment through layers of verification, the role of multiple hub ceRNAs in AML is still not completely clear, and the exploration and verification of the functions of these RNAs is just beginning.

Conclusion

In summary, our research combined ESTIMATE and CIBERSORT algorithms to identify AML TME related genes and immune cells, screened out ceRNA network that have high prognostic predictive power and highly TME association through multiple rigorous survival models and survival analysis. This study proposed a new way to reveal the role of TME in AML from the perspective of post transcriptional regulation, and contributed to exploring more diverse and effective AML biological markers.

Data Availability Statement

The datasets analyzed for this study can be found in The Cancer Genome Atlas database (https://portal.gdc.cancer.gov/) and GEO (GSE142699–GPL26945, GSE71014–GPL10558) (https://www.ncbi.nlm.nih.gov/geo).

Author Contributions

Conceptualization, YC and XW; methodology, PQ, CXL, SW, and LJ; software, XW, HD, QW and XS; validation, XW, YRL, and LY; formal analysis, YC, PQ and YS; investigation, YL, CYL, and CL; resources, XW and ZW; writing—original draft preparation, YC, CXL, and XW; writing—review and editing, XW and ZW. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by The National Key R&D program of China (2018YFC1106000).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.584884/full#supplementary-materials

References

1. Torre LA, Bray F, Siegel RL, Ferlay J, Lortet-Tieulent J, Jemal A. Global cancer statistics, 2012. CA: Cancer J Clin (2015) 65(2):87–108. doi: 10.3322/caac.21262

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Dohner H, Weisdorf DJ, Bloomfield CD. Acute Myeloid Leukemia. N Engl J Med (2015) 373(12):1136–52. doi: 10.1056/NEJMra1406184

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Papaemmanuil E, Gerstung M, Bullinger L, Gaidzik VI, Paschka P, Roberts ND, et al. Genomic Classification and Prognosis in Acute Myeloid Leukemia. N Engl J Med (2016) 374(23):2209–21. doi: 10.1056/NEJMoa1516192

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Burnett A, Wetzler M, Lowenberg B. Therapeutic advances in acute myeloid leukemia. J Clin Oncol (2011) 29(5):487–94. doi: 10.1200/JCO.2010.30.1820

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Christopher MJ, Petti AA, Rettig MP, Miller CA, Chendamarai E, Duncavage EJ, et al. Immune Escape of Relapsed AML Cells after Allogeneic Transplantation. N Engl J Med (2018) 379(24):2330–41. doi: 10.1056/NEJMoa1808777

PubMed Abstract | CrossRef Full Text | Google Scholar

6. van Galen P, Hovestadt V, Wadsworth Ii MH, Hughes TK, Griffin GK, Battaglia S, et al. Single-Cell RNA-Seq Reveals AML Hierarchies Relevant to Disease Progression and Immunity. Cell (2019) 176(6):1265–81.e24. doi: 10.1016/j.cell.2019.01.031

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Duarte D, Hawkins ED, Lo Celso C. The interplay of leukemia cells and the bone marrow microenvironment. Blood (2018) 131(14):1507–11. doi: 10.1182/blood-2017-12-784132

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Shafat MS, Gnaneswaran B, Bowles KM, Rushworth SA. The bone marrow microenvironment - Home of the leukemic blasts. Blood Rev (2017) 31(5):277–86. doi: 10.1016/j.blre.2017.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Karreth FA, Pandolfi PP. ceRNA cross-talk in cancer: when ce-bling rivalries go awry. Cancer Discov (2013) 3(10):1113–21. doi: 10.1158/2159-8290.CD-13-0202

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Liu Y, Cheng Z, Pang Y, Cui L, Qian T, Quan L, et al. Role of microRNAs, circRNAs and long noncoding RNAs in acute myeloid leukemia. J Hematol Oncol (2019) 12(1):51. doi: 10.1186/s13045-019-0734-5

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Yan H, Qu J, Cao W, Liu Y, Zheng G, Zhang E, et al. Identification of prognostic genes in the acute myeloid leukemia immune microenvironment based on TCGA data analysis. Cancer Immunol Immunother CII (2019) 68(12):1971–8. doi: 10.1007/s00262-019-02408-7

CrossRef Full Text | Google Scholar

12. Wang JD, Zhou HS, Tu XX, He Y, Liu QF, Liu Q, et al. Prediction of competing endogenous RNA coexpression network as prognostic markers in AML. Aging (2019) 11(10):3333–47. doi: 10.18632/aging.101985

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Xu M, Li Y, Li W, Zhao Q, Zhang Q, Le K, et al. Immune and Stroma Related Genes in Breast Cancer: A Comprehensive Analysis of Tumor Microenvironment Based on the Cancer Genome Atlas (TCGA) Database. Front Med (Lausanne) (2020) 7:64. doi: 10.3389/fmed.2020.00064

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Jia D, Li S, Li D, Xue H, Yang D, Liu Y. Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging (2018) 10(4):592–605. doi: 10.18632/aging.101415

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Wang Z, Liu T, Li G, Cao Z. The exploration of new therapeutic targets for HPV-negative head and neck squamous cell cancer through the construction of a ceRNA network and immune microenvironment analysis. J Cell Biochem (2020) 121(5-6):3426–37. doi: 10.1002/jcb.29615

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Huang R, Zeng Z, Li G, Song D, Yan P, Yin H, et al. The Construction and Comprehensive Analysis of ceRNA Networks and Tumor-Infiltrating Immune Cells in Bone Metastatic Melanoma. Front Genet (2019) 10:828. doi: 10.3389/fgene.2019.00828

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Metsalu T, Vilo J. ClustVis: a web tool for visualizing clustering of multivariate data using Principal Component Analysis and heatmap. Nucleic Acids Res (2015) 43(W1):W566–70. doi: 10.1093/nar/gkv468

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Jeggari A, Marks DS, Larsson E. miRcode: a map of putative microRNA target sites in the long non-coding transcriptome. Bioinformatics (2012) 28(15):2062–3. doi: 10.1093/bioinformatics/bts344

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. eLife (2015) 4:e05005. doi: 10.7554/eLife.05005

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wong N, Wang X. miRDB: an online resource for microRNA target prediction and functional annotations. Nucleic Acids Res (2015) 43(Database issue):D146–52. doi: 10.1093/nar/gku1104

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Huang R, Meng T, Chen R, Yan P, Zhang J, Hu P, et al. The construction and analysis of tumor-infiltrating immune cell and ceRNA networks in recurrent soft tissue sarcoma. Aging (2019) 11(22):10116–43. doi: 10.18632/aging.102424

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Shallis RM, Wang R, Davidoff A, Ma X, Zeidan AM. Epidemiology of acute myeloid leukemia: Recent progress and enduring challenges. Blood Rev (2019) 36:70–87. doi: 10.1016/j.blre.2019.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Xu Q, Li Y, Lv N, Jing Y, Xu Y, Li Y, et al. Correlation Between Isocitrate Dehydrogenase Gene Aberrations and Prognosis of Patients with Acute Myeloid Leukemia: A Systematic Review and Meta-Analysis. Clin Cancer Res an Off J Am Assoc Cancer Res (2017) 23(15):4511–22. doi: 10.1158/1078-0432.CCR-16-2628

CrossRef Full Text | Google Scholar

26. Thol F, Schlenk RF, Heuser M, Ganser A. How I treat refractory and early relapsed acute myeloid leukemia. Blood (2015) 126(3):319–27. doi: 10.1182/blood-2014-10-551911

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Corre J, Mahtouk K, Attal M, Gadelorge M, Huynh A, Fleury-Cappellesso S, et al. Bone marrow mesenchymal stem cells are abnormal in multiple myeloma. Leukemia (2007) 21(5):1079–88. doi: 10.1038/sj.leu.2404621

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell (2011) 144(5):646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Huang S, Zhang B, Fan W, Zhao Q, Yang L, Xin W, et al. Identification of prognostic genes in the acute myeloid leukemia microenvironment. Aging (2019) 11(22):10557–80. doi: 10.18632/aging.102477

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Dong X, Fang Z, Yu M, Zhang L, Xiao R, Li X, et al. Knockdown of Long Noncoding RNA HOXA-AS2 Suppresses Chemoresistance of Acute Myeloid Leukemia via the miR-520c-3p/S100A4 Axis. Cell Physiol Biochem Int J Exp Cell Physiol Biochem Pharmacol (2018) 51(2):886–96. doi: 10.1159/000495387

CrossRef Full Text | Google Scholar

31. Lei B, He A, Chen Y, Cao X, Zhang P, Liu J, et al. Long non-coding RNA RPPH1 promotes the proliferation, invasion and migration of human acute myeloid leukemia cells through down-regulating miR-330-5p expression. EXCLI J (2019) 18:824–37. doi: 10.17179/excli2019-1686

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Ulger C, Toruner GA, Alkan M, Mohammed M, Damani S, Kang J, et al. Comprehensive genome-wide comparison of DNA and RNA level scan using microarray technology for identification of candidate cancer-related genes in the HL-60 cell line. Cancer Genet Cytogenet (2003) 147(1):28–35. doi: 10.1016/S0165-4608(03)00155-9

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yu F, Yang H, Zhang Z, Wang Z, Xiong J. DAL-1/4.1B contributes to epithelial-mesenchymal transition via regulation of transforming growth factor-beta in lung cancer cell lines. Mol Med Rep (2015) 12(4):6072–8. doi: 10.3892/mmr.2015.4217

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Bataille R, Chappard D, Marcelli C, Dessauw P, Sany J, Baldet P, et al. Mechanisms of bone destruction in multiple myeloma: the importance of an unbalanced process in determining the severity of lytic bone disease. J Clin Oncol (1989) 7(12):1909–14. doi: 10.1200/JCO.1989.7.12.1909

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Narakornsak S, Poovachiranon N, Peerapapong L, Pothacharoen P, Aungsuchawan S. Mesenchymal stem cells differentiated into chondrocyte-Like cells. Acta Histochem (2016) 118(4):418–29. doi: 10.1016/j.acthis.2016.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Alvarez J, Balbin M, Santos F, Fernandez M, Ferrando S, Lopez JM. Different bone growth rates are associated with changes in the expression pattern of types II and X collagens and collagenase 3 in proximal growth plates of the rat tibia. J Bone Miner Res (2000) 15(1):82–94. doi: 10.1359/jbmr.2000.15.1.82

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Ganapathi MK, Jones WD, Sehouli J, Michener CM, Braicu IE, Norris EJ, et al. Expression profile of COL2A1 and the pseudogene SLC6A10P predicts tumor recurrence in high-grade serous ovarian cancer. Int J Cancer (2016) 138(3):679–88. doi: 10.1002/ijc.29815

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Antonov SM, Dudel J, Franke C, Hatt H. Argiopine blocks glutamate-activated single-channel currents on crayfish muscle by two mechanisms. J Physiol (1989) 419:569–87. doi: 10.1113/jphysiol.1989.sp017887

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Jiang W, Min J, Sui X, Qian Y, Liu Y, Liu Z, et al. MicroRNA-26a-5p and microRNA-23b-3p up-regulate peroxiredoxin III in acute myeloid leukemia. Leuk Lymph (2015) 56(2):460–71. doi: 10.3109/10428194.2014.924115

CrossRef Full Text | Google Scholar

40. Salvatori B, Iosue I, Mangiavacchi A, Loddo G, Padula F, Chiaretti S, et al. The microRNA-26a target E2F7 sustains cell proliferation and inhibits monocytic differentiation of acute myeloid leukemia cells. Cell Death Dis (2012) 3:e413. doi: 10.1038/cddis.2012.151

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Salvatori B, Iosue I, Djodji Damas N, Mangiavacchi A, Chiaretti S, Messina M, et al. Critical Role of c-Myc in Acute Myeloid Leukemia Involving Direct Regulation of miR-26a and Histone Methyltransferase EZH2. Genes Cancer (2011) 2(5):585–92. doi: 10.1177/1947601911416357

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Huang H, Jiang X, Wang J, Li Y, Song CX, Chen P, et al. Identification of MLL-fusion/MYC dash, verticalmiR-26 dash, verticalTET1 signaling circuit in MLL-rearranged leukemia. Cancer Lett (2016) 372(2):157–65. doi: 10.1016/j.canlet.2015.12.032

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Wang Y, Li Z, He C, Wang D, Yuan X, Chen J, et al. MicroRNAs expression signatures are associated with lineage and survival in acute leukemias. Blood Cells Mol Dis (2010) 44(3):191–7. doi: 10.1016/j.bcmd.2009.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Wang XX, Zhang R, Li Y. Expression of the miR-148/152 Family in Acute Myeloid Leukemia and its Clinical Significance. Med Sci Monit Int Med J Exp Clin Res (2017) 23:4768–78. doi: 10.12659/MSM.902689

CrossRef Full Text | Google Scholar

45. Wang XX, Zhang H, Li Y. Preliminary study on the role of miR148a and DNMT1 in the pathogenesis of acute myeloid leukemia. Mol Med Rep (2019) 19(4):2943–52. doi: 10.3892/mmr.2019.9913

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Cheng Y, Su Y, Wang S, Liu Y, Jin L, Wan Q, et al. Identification of circRNA-lncRNA-miRNA-mRNA Competitive Endogenous RNA Network as Novel Prognostic Markers for Acute Myeloid Leukemia. Genes (2020) 11(8):868. doi: 10.3390/genes11080868

CrossRef Full Text | Google Scholar

47. Zhang N, Chen Y, Shen Y, Lou S, Deng J. Comprehensive analysis the potential biomarkers for the high-risk of childhood acute myeloid leukemia based on a competing endogenous RNA network. Blood Cells Mol Dis (2019) 79:102352. doi: 10.1016/j.bcmd.2019.102352

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Condeelis J, Pollard JW. Macrophages: obligate partners for tumor cell migration, invasion, and metastasis. Cell (2006) 124(2):263–6. doi: 10.1016/j.cell.2006.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Xu ZJ, Gu Y, Wang CZ, Jin Y, Wen XM, Ma JC, et al. The M2 macrophage marker CD206: a novel prognostic indicator for acute myeloid leukemia. Oncoimmunology (2020) 9(1):1683347. doi: 10.1080/2162402X.2019.1683347

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Wu D, Wu P, Qiu F, Wei Q, Huang J. Human gammadeltaT-cell subsets and their involvement in tumor immunity. Cell Mol Immunol (2017) 14(3):245–53. doi: 10.1038/cmi.2016.55

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Jones WD, Michener CM, Biscotti C, Braicu I, Sehouli J, Ganapathi MK, et al. RNA Immune Signatures from Pan-Cancer Analysis Are Prognostic for High-Grade Serous Ovarian Cancer and Other Female Cancers. Cancers (2020) 12(3):620. doi: 10.3390/cancers12030620

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: acute myeloid leukemia, tumor microenvironment, ceRNA network, prognosis, immune microenvironment

Citation: Cheng Y, Wang X, Qi P, Liu C, Wang S, Wan Q, Liu Y, Su Y, Jin L, Liu Y, Li C, Sang X, Yang L, Liu C, Duan H and Wang Z (2021) Tumor Microenvironmental Competitive Endogenous RNA Network and Immune Cells Act as Robust Prognostic Predictor of Acute Myeloid Leukemia. Front. Oncol. 11:584884. doi: 10.3389/fonc.2021.584884

Received: 18 July 2020; Accepted: 08 March 2021;
Published: 09 April 2021.

Edited by:

Zhenyu Jia, University of California, Riverside, United States

Reviewed by:

Kuang-Ming Liao, Chi Mei Medical Center, Taiwan
Xu Ye, The Second Affiliated Hospital of Guangzhou Medical University, China

Copyright © 2021 Cheng, Wang, Qi, Liu, Wang, Wan, Liu, Su, Jin, Liu, Li, Sang, Yang, Liu, Duan and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Xiaoran Wang, xiaoranw1@126.com; Zhichong Wang, wangzhichong@gzzoc.com

These authors have contributed equally to this work

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.