Analysis of m6A-related signatures associated with the tumor immune microenvironment and predict survival in acute myeloid leukemia
Original Article

Analysis of m6A-related signatures associated with the tumor immune microenvironment and predict survival in acute myeloid leukemia

Shushu Yuan1,2#, Zhirong Cong2#, Jiali Ji2#, Li Zhu2, Qi Jiang2, Ying Zhou2, Qian Shen2, Daniela Damiani3, Xiaohong Xu2, Bingzong Li1

1Department of Hematology, The Second Affiliated Hospital of Soochow University, Suzhou, China; 2Department of Oncology, Affiliated Tumor Hospital of Nantong University, Nantong, China; 3Division of Hematology and Stem Cell Transplantation, Azienda Sanitaria Universitaria Friuli Centrale-Ospedale S. M. Misericordia, Udine, Italy

Contributions: (I) Conception and design: S Yuan, X Xu, B Li; (II) Administrative support: X Xu, B Li; (III) Provision of study materials or patients: S Yuan, Z Cong, J Ji; (IV) Collection and assembly of data: L Zhu, Q Jiang; (V) Data analysis and interpretation: Y Zhou, Q Shen; (VI) Manuscript writing: All authors; (VII) Final approval of manuscript: All authors.

#These authors contributed equally to this work and should be considered as co-first authors.

Correspondence to: Bingzong Li, MD, PhD. Department of Hematology, The Second Affiliated Hospital of Soochow University, San Xiang Road 1055, Suzhou 215006, China. Email: lbzwz0907@hotmail.com; Xiaohong Xu. Department of Oncology, Affiliated Tumor Hospital of Nantong University, Nantong 226001, China. Email: xhx107@163.com.

Background: Most previous studies have focused on the intrinsic carcinogenic pathways of tumors; however, little is known about the potential role of N6-methyladenosine (m6A) methylation in the tumor immune microenvironment (TIME). To better diagnose and treat acute myeloid leukemia (AML), we sought to examine the correlation between m6A regulatory factors and immune infiltration in cases of AML. At the same time, a prognostic model was constructed to predict the survival of AML.

Methods: We extracted data from The Cancer Genome Atlas (TCGA) database, including ribonucleic acid sequencing (RNA-seq) transcriptome data and data on the corresponding clinical characteristics of AML patients. We identified two m6A modification patterns with distinct clinical outcomes and found a significant relationship between them. Simultaneous discovery of distinct m6A clusters associated with the tumor immune microenvironment [immune cell types and Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm] are closely related. Next, we implemented Lasso (Least Absolute Shrinkage and Selection Operator) Cox regression to build a predictive model in the 2-m6A regulator TCGA dataset to further explore m6A prognostic features in AML, and perform correlation validation.

Results: We identified 2 molecular subtypes (Clusters 1 and 2) by the consistent clustering of significant m6A regulators in AML. Cluster 2 was associated with a higher immune score and obvious immune cell infiltration, and thus patients in Cluster 2 had a poorer prognosis than those in Cluster 1 (P<0.05). Additionally, the 2 m6A-related signatures representing the independent prognostic factors in AML were screened to construct a prognostic risk-score model. We found that patients with low-risk scores had higher immune scores than those with high-risk scores (P<0.05).

Conclusions: Our research confirmed that m6A methylation plays an important role in AML. Further provide new directions for the prognosis and treatment of AML.

Keywords: Tumor immune microenvironment (TIME); m6A; survival; acute myeloid leukemia (AML)


Submitted Jun 09, 2022. Accepted for publication Aug 18, 2022.

doi: 10.21037/atm-22-3858


Introduction

Acute myeloid leukemia (AML) is a group of diseases derived from the malignant transformation of hematopoietic progenitor cells at different stages during the differentiation and development of myeloid cells (1). The number of cases of AML increased from 63.84×103 in 1990 to 119.57×103 cases in 2017; thus, there has been an 87.3% increase in the incidence of AML over the past 27 years worldwide (2). AML can infiltrate various organs, such as the liver, spleen, and lymph nodes. Its clinical manifestations include anemia, bleeding, infection, and fever, which are critical conditions with poor prognosis. The survival of AML patients is short; AML patients have a 2-year survival rate of 32% and a 5-year survival rate of 24% (3). The treatment of AML is divided into general treatment, chemotherapy, targeted therapy, autologous and allogeneic hematopoietic stem cell transplantation, etc. (4,5). Recently, encouraging progress has been made due to the development of existing technologies; however, the treatment results of patients with AML are still not satisfactory, and more than half of the patients die from the disease (6,7). The development of AML is not only related to gene mutation and chromosomal variation, but also abnormal epigenetic regulation such as DNA methylation and histone modification play an important role in the occurrence of AML (1).

N6-methyladenosine (m6A) is the most common post-transcriptional modification of eukaryotic messenger ribonucleic acid (mRNA), and accounts for 80% of RNA methylation modifications (8). The execution of the function is mainly completed by 3 types of regulators (i.e., writers, erasers, and readers) (9). As early as the 1970s, m6A modifications were discovered in the mRNA and long non-coding RNA of eukaryotes. In most eukaryotes, the methylation modification that occurs in the 3'untranslated region (UTR) of mRNA plays an important role in mRNA splicing, editing, stability, degradation, and polyadenylation (10). The modification of the 3'UTR methylation contributes to the nuclear transport of mRNA, the initiation of translation, and the maintenance of the structural stability of mRNA combined with the polyA binding protein (11). Further, research has shown that m6A modifiers are correlated with AML; for example, Huang et al. found that fat mass and obesity-associated protein (FTO) is a mRNA m6A demethylase that can be targeted by small molecule inhibitors and has the potential to treat AML (12). Methyltransferase-like 14 (METTL14) modulates its mRNA target through m6A modification to exert its carcinogenic effect, and METTL14 and m6A modification have key roles in normal and malignant hematopoiesis (13,14).

Immune checkpoints refer to interactions that inhibit cytotoxic T lymphocyte (CTL) activation. These checkpoints act as a “brake” for CTL, allowing CTL to activate for a while but not allowing the reaction to proceed indefinitely. This mechanism is one of many regulatory mechanisms of the immune response. Programmed death-ligand 1 (PD-L1) is one of the most important immune checkpoints (15-17). Yang et al. showed that programmed death-1 (PD-1)/PD-L1 signaling might be involved in the pathogenesis of myelodysplastic syndromes and resistance to hypomethylation drugs, which provide potential treatments for myelodysplastic syndromes and AML (18). PD-1(−/−) mice challenged with C1498 cells produced enhanced anti-tumor T cell responses, had a reduced burden of AML in the blood and other organs, and survived significantly longer than wild-type mice, which shows the importance of the PD-1/PD-L1 pathway in the immune evasion of hematological malignancies (19,20). PD-L1 is closely related to the tumor immune microenvironment (TIME). However, very few studies have examined the relationship between the m6A-related signatures and TIME. There is also a study that enhances the understanding of the heterogeneity and complexity of TIME through a comprehensive analysis of specific m6A modulators to guide related immunotherapy (14). Unfortunately, these studies did not extend to AML.

Although current gene-targeted mutation testing and minimal residual disease (MRD) testing can provide preliminary prognostic stratification for AML. However, according to the clinical application results, it is far from enough to accurately predict the prognosis. Therefore, it is critical to uncover the genomic properties of AML, construct prognostic models for developing effective treatments and predicting individual survival and risk of relapse. This study explored the genes related to m6A in AML, examined the correlation between cluster typing and immune cell infiltration, and performed clinical modeling to evaluate the survival prognosis of AML patients. Our findings provide a basis for the molecular diagnosis and targeted therapy of AML in the future. We present the following article in accordance with the TRIPOD reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-22-3858/rc).


Methods

Data selection

We extracted all data, including the ribonucleic acid sequencing (RNA-seq) transcriptome data and data on the corresponding clinical characteristics of AML patients, from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). A total of 151 specimens and relevant clinicopathological information embodying sex and age were downloaded for further analysis. From all samples, 129 AML patient samples with complete prognosis and relevant clinical data were selected. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

M6A signatures

Previous studies have identified 23 m6A-related regulators (21-23). We evaluated the differential expression of m6A signatures between the normal samples and AML samples.

Cancer clustering

The “ConsensusClusterPlus” package in R software was used for the consistent clustering analysis in the 1,000 repeat iterations to determine the optimal number of clusters, K, and sort the AML patients into subgroups. This method was used to verify the rationality of clustering by resampling. A gene set enrichment analysis (GSEA) was conducted to evaluate the distribution trends of genes in the gene table ranked by phenotype correlation and to examine their contribution to the phenotype. We also identified the enrichment pathways of 2 subtypes with a false discovery rate <0.05 and calculated the normalized enrichment score (NES).

Immune score analysis

The “estimate” R package was used to calculate the stromal scores, estimate scores, and immune scores of the AML patients using the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm. CIBERSORT is a tool used to deconvolute the expression matrix of immune cell subtypes based on the principle of linear support vector regression. The RNA-seq data was used to estimate immune cell infiltration.

Construction of a prognosis model

A univariate Cox regression analysis was conducted to select 5 m6A signatures. All the AML specimens were categorized into the modeling group and the verification group at a ratio of 1:1 by applying the caret package. Next, a least absolute shrinkage and selection operator (LASSO) regression was conducted to filter out 2 candidate m6A-related variables, while fitting a generalized linear model in TCGA modeling set. The following formula was used to calculate the risk score: Risk score = (−0.0397) × YTHDF3 + 0.0222 × ALKBH5. Using the above formula, patients were divided into low-risk and high-risk groups based on the mean value of the risk scores. To improve the predictive accuracy and interpretability of statistical models, Lasso Cox regression analysis was performed to examine the relationship between m6A prognostic features and AML risk. Cox proportional hazards regression models included age, sex, and risk score. The hazard ratio (HR) from Cox regression analyses was used to distinguish between positive and negative prognostic factors. Genes with HR >1 were considered risk genes, and genes with HR <1 were considered protection genes. Subsequently, the availability of prognostic models was assessed using the Kaplan-Meier survival method, and the prognostic accuracy of signature construction was assessed using the sensitivity and specificity of receiver operating characteristic (ROC) curves.

Statistics

R version 4.1.0 software was mainly used in the statistical analysis and to generate images. We used the Kaplan-Meier (K-M) method to account for differences in the survival rates of the 2 groups. A subgroup analysis was conducted to evaluate the stability of risk characteristics, and patients were divided into 2 subtypes according to their age (≤65 or >65 years) and sex (female or male). The ROC curve image and the corresponding area under the ROC curve was used to assess the predicted value of the risk model. The value of area under the curve (AUC) ranges between 0.5 and 1. The closer the AUC is to 1.0, the higher the authenticity of the detection method; when it is equal to 0.5, the authenticity is the lowest and has no application value. The differences between the two groups were analyzed using the Student’s t-test and the one-way analysis of variance test was used for the comparison of multiple conditions. The “Spearman” method was used to calculate the gene expression correlations. In this study, a bilateral P value <0.05 was considered statistically significant.


Results

Expression of m6A-related signatures in AML patients

A total of 8 writer regulators (i.e., METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15, and RBM15B), 13 reader regulators [i.e., YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTH M6A RNA binding protein 3 (YTHDF3), HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3, and RBMX], and 2 eraser regulators [i.e., FTO and alkB homolog 5 (ALKBH5)] were identified between the tumor samples and the normal samples. Additionally, a total of 5 m6A-related signatures were examined based on the univariate Cox regression analysis (applying the cutoff criteria of a P value <0.1), and their expression levels in AML were determined (see Figure 1). Higher expression levels of METTL14, YTHDC2, YTHDF, and HNRNPA2B1 were significantly correlated with better survival. Conversely, higher expression levels of ALKBH5 were significantly correlated with a worse prognosis (see Table 1).

Figure 1 The 5 M6A-related markers of differentially expressed genes between AML patients and normal individuals. (A) Single-factor Cox regression analysis. (B) Expression of related differentially expressed genes in AML. (C) Heatmap of related genes in AML. AML, acute myeloid leukemia.

Table 1

Univariate Cox regression analyses of the 5 m6A-related signatures associated with the OS of the AML patients

Gene HR HR.95L HR.95H P value
METTL14 0.851978 0.705265 1.029209 0.096638
YTHDC2 0.88141 0.774514 1.00306 0.055665
YTHDF3 0.94964 0.913381 0.987338 0.009281
HNRNPA2B1 0.99129 0.982991 0.99966 0.041424
ALKBH5 1.050996 1.011429 1.09211 0.011073

OS, overall survival; AML, acute myeloid leukemia; HR, hazard ratio.

Association between m6A-related genes and consensus clustering analysis

A consensus clustering analysis was performed on all AML patient samples using the combination of the identified probes. Following 1,000 clustering iterations at cluster counts (k) ranging from 2 to 8, k=2 was found to be the best k with the lowest proportion of ambiguous clustering (PAC) measures. This framework generated 2 distinct clusters, both with consensuses >0.99 (see Figure 2). A total of 129 specimens with AML were divided into Cluster 1 (n=59) and Cluster 2 (n=70). To further explore the relationship between the 5 m6A-related genes and clinical features, including age, sex, and cluster, we generated a heatmap showing the expression of these signatures (see Figure 3A). The K-M curve analysis showed that there was no significant difference between the 2 clusters (see Figure 3B). Additionally, as Figure 3C shows, the 5 m6A genes were correlated to each other. PD-L1 was positively correlated with YTHDC2, while PD-L1 was negatively correlated with HNRNPA2B1 and ALKBH5. The expression of PD-L1 was differentially increased in Cluster 2 compared to Cluster 1 (see Figure 3D).

Figure 2 Consensus clustering analysis of the AML patients. (A) Consensus CDF plot. This figure illustrates the CDF plot of the consensus matrix for each k (stratified by colors). In any curve of a consensus matrix, the lower left portion represents the samples that are rarely clustered together, the upper right portion represents those that are always clustered together, and the middle segment represents those with ambiguous assignments in different iterations. The PAC measure is defined as the CDF value of the sample pairs with a consensus index =1 minus the CDF value for those close to 0 (top right vs. lower left). A low PAC value indicates a flat middle segment and a low rate of discordance across permuted clustering runs. The PAC measure for clustering at k=2 was the smallest (almost zero), which suggested that k=2 was the best cluster count in this analysis. (B) Delta area plot. This figure shows the relative increase in the clustering consensus following a change from k−1 to k. The k providing the largest change in the consensus represents the best cluster count. In this plot, after k=2, none of the other k values resulted in a considerable increase in the consensus. (C-I) Consensus cluster heatmaps at k=2–8. Both rows and columns are samples. The intensity of the blue represents how often every 2 samples from the rows and columns clustered together in the 1,000 repeat iterations. The pane on top demonstrates the consensus cluster for each sample. AML, acute myeloid leukemia; CDF, cumulative distribution function; PAC, proportion of ambiguous clustering.
Figure 3 Association with m6A-related genes. (A) Evaluation of heatmaps of the 5 m6A-related genes associated with AML patients based on different clinical characteristics. (B) K-M curve of AML OS. (C) Correlation between PD-L1 and the 5 m6A-related genes. (D) PD-L1 expression level in AML. *, P<0.05; ns, P>0.05. AML, acute myeloid leukemia; K-M, Kaplan-Meier; OS, overall survival; PD-L1, programmed death-ligand 1.

Figure 4A shows the results of the differential analysis in immune cell infiltration between the two clusters. These boxplots also show that the expression levels of memory-activated cluster of differentiation (CD)4 T cells, naive CD4 T cells, and eosinophils in Cluster 1 were higher than those in Cluster 2. Conversely, the expression levels of CD8 T cells and regulatory T cells in Cluster 1 were lower than those in Cluster 2 (see Figure 4B-4F). Additionally, the ESTIMATE score of Cluster 2 was higher than that of Cluster 1 (P<0.05; see Figure 4G) when grouped by stromal score and immune score, respectively, and differences were observed between Cluster 1 and Cluster 2 (see Figure 4H,4I).

Figure 4 Differential analysis of immune cell infiltration in the 2 clusters. (A) The proportion of the 22 immune cell subgroups. (B) The proportion of the activated memory CD4 T cells displayed by box plot (P=0.056). (C) The proportion of CD4 naive T cells displayed by box plot (P=0.0082). (D) The proportion of CD8 T cells displayed by box plot (P=0.04). (E) The proportion of regulatory T cells displayed by box plot (P=0.0029). (F) The proportion of eosinophils displayed by box plot (P=9.1e–05). (G) ESTIMATE score level displayed by box plot (P=0.045). (H) Stromal score level displayed by box plot (P=0.23). (I) Immune score level displayed by box plot (P=0.022). NK, natural killer.

To clarify the underlying regulatory mechanism causing the difference between the 2 subtypes, we described the top pathways in the different clusters based on the GSEA. The outcome illustrated that spliceosome (NES =−1.85, normalized P=0.01), protein export (NES =−1.71, normalized P=0.02), and basal transcription (NES =−1.69, normalized P=0.025) factors were significantly correlated with Cluster 1 (see Figure 5A-5C). Additionally, Cluster 2 was enriched in cell adhesion molecules (CAMs) (NES =2.03, normalized P<0.001), arrhythmogenic right ventricular cardiomyopathy (ARVC) (NES =1.94, normalized P<0.001) and the notch signaling pathway (NES =1.91, normalized P<0.001) (see Figure 5D-5F).

Figure 5 Enrichment plots from the GSEA. Gene set enrichment plots of (A) spliceosome; (B) protein export; (C) basal transcription factors; (D) CAMs; (E) ARVC; and (F) the notch signaling pathway. GSEA, gene set enrichment analysis; CAM, cell adhesion molecule; ARVC, arrhythmogenic right ventricular cardiomyopathy; NES, normalized enrichment score; NOM, nominal; FDR, false discovery rate.

Construction of the prognostic model of the m6A signatures

A total of 129 patients with AML were categorized into a training set (n=65) and a validation set (n=64) at a rough ratio of 1:1. Next, 2 m6A signatures (i.e., YTHDF3 and ALKBH5) were selected to build a prognostic model according to the LASSO regression analysis of the training set (see Figure 6A,6B). The patients were then divided into a high-risk score group and a low-risk score group based on the median risk scores. To further examine the stability and effectiveness of the established model, we used K-M plots to analyze both the training and validation sets and found that the low-risk score cohort had a better survival rate than the high-rick score cohort (both P<0.05; see Figure 6C,6D). The AUC of the training set reached 0.774, and that of the validation set was 0.665 (see Figure 6E,6F). These results suggest that the risk model can be used as an important indicator for evaluating the prognosis of AML, but it is not the only criterion, and we can increase the sample size and repeat the validation in future clinical data analysis. We also evaluated the relationship between the risk score subgroups, overall survival (OS), OS status, and the level of expression, and the 2 m6A signatures. We found that the expression of YTHDF3 in the training cohort was lower in the high-risk group and higher in the low-risk group, while the expression of ALKBH5 showed the opposite trend. The same conclusion was drawn in relation to YTHDF3 in the validation cohort, but the expression of ALKBH5 did not change significantly (see Figure 7).

Figure 6 Prognostic value of the m6A-based prognostic model in patients stratified by AML subgroup. (A) LASSO regression analysis of the training set. (B) Cross-validation plot for the penalty term. (C) K-M curves for the training set of the AML patients. (D) K-M curves for the validation set of the AML patients. (E) Time-dependent ROC curves showed the predictive efficiency of the m6A-based prognostic model in the training set of the AML patients. (F) Time-dependent ROC curves showing the predictive efficacy of the m6A-based prognostic model in the validation set of the AML patients. AML, acute myeloid leukemia; LASSO, least absolute shrinkage and selection operator; K-M, Kaplan-Meier; ROC, receiver operating characteristic; AUC, area under the curve.
Figure 7 The signature-based risk score was a promising marker in the training and validation cohorts. (A) Survival overview of the training set. (B) Risk score distribution of the training set. (C) Heatmap of the training set showing the expression profiles of the signature in the low- and high-risk groups. (D) Survival overview of the validation set. (E) Risk score distribution of the validation set. (F) Heatmap of the validation set showing the expression profiles of the signature in the low-risk and high-risk groups.

Clinical characteristics, such as age and sex, were examined in the subgroup analysis. From the univariate and multivariate Cox regression analyses, we found that age [hazard ratio (HR) 2.187 (1.155−4.142); P=0.016] and risk scores [HR 6.253 (1.097−35.634); P=0.039] were independent risk factors in the training cohort (see Figure 8A-8D). Notably, in relation to OS, the low-risk score group had a longer prognosis than the high-risk score group in AML patients aged ≤65 years and males or females (all P<0.05) (see Figure 8E-8H). The heatmap shows the expression of YTHDF3 and ALKBH5 in the different clinical subgroups of AML. YTHDF3 was more lowly expressed in the AML high-risk group, while ALKBH5 was more highly expressed in the AML high-risk group (see Figure 9A). Additionally, patients aged >65 years had a higher risk score than those aged ≤65 years (P=0.042; see Figure 9B). The results of the differential analysis revealed that Cluster 1 had a lower risk score than Cluster 2 (P=0.00029) (see Figure 9C). Notably, a higher immune score was associated with an increased risk score (P=0.026) (see Figure 9D). There was no statistically significant difference in PD-L1 between the low-risk and high-risk score groups (P=0.29; see Figure 9E).

Figure 8 Analysis of the clinical parameters of the AML patients and the prognosis of OS. (A) Univariate regression analysis of the training set of the AML patients. (B) Multivariate regression analysis of the training set of the AML patients. (C) Univariate regression analysis of the validation set of the AML patients. (D) Multivariate regression analysis of the validation set of the AML patients. (E) K-M curve showing the OS of the AML patients aged ≤65 years at high or low risk. (F) K-M curve showing the OS of the AML patients aged >65 years at high or low risk. (G) K-M showing the OS of the male AML patients at high or low risk. (H) K-M curve showing the OS of the female AML patients at high or low risk. AML, acute myeloid leukemia; OS, overall survival; K-M, Kaplan-Meier.
Figure 9 Analysis of the difference between the risk score of the AML patients and PD-L1. (A) Heatmap showing the expression of YTHDF3 and ALKBH5 in different clinical subgroups. (B) Box plot assessing the risk of the AML patients according to age. (C) Box plot assessing the risk of the AML patients in different groups. (D) Box plot assessing the risk of the AML patients based on the immune score. (E) Box plot assessing PD-L1 expression of the AML patients based on the risk score. **, P<0.01. AML, acute myeloid leukemia; PD-L1, programmed death-ligand 1.

The relationship between the immune microenvironment and disease prognosis

To further explore the relationship between the risk scores and immune cell infiltration, we conducted a correlation analysis, and found that the expression levels of the resting CD4 memory T cells (R=−0.18, P=0.046) and eosinophils (R=−0.21, P=0.02) were positively correlated to the risk scores, while the expression levels of the M2 macrophages (R=0.42, P=2.3e−06) and memory B cells (R=0.26, P=0.0041) were positively correlated to the risk scores (see Figure 10).

Figure 10 Associations of the risk score with immune cell infiltration in AML. (A) Correlation between the risk score and resting CD4 memory cells. (B) Correlation between the risk score and CD8 T cells. (C) Correlation between the risk score and the M2 macrophages. (D) Correlation between the risk score and the B cell memory. (E) Correlation between the risk score and the eosinophils. AML, acute myeloid leukemia.

Discussion

M6A methylase is heavily involved in mammalian development, immunity, tumor generation and metastasis, stem cell renewal, fat differentiation, and other life processes (12,24-26). Li et al. showed that FTO, which is a m6A demethylase, has a key carcinogenic effect in AML (27). AML is an aggressive clonal disease of hematopoietic stem cells and primitive progenitor cells that can prevent bone marrow differentiation and produce self-renewing leukemia stem cells (LSCs) (28). Paris et al. revealed that the mRNA m6A reader YTHDF2 was overexpressed in a broad spectrum of human AML and was a unique therapeutic target that inhibited the selective targeting of LSCs and promoted hematopoietic stem cell amplification (29,30). To date, previous studies have mostly focused on the internal mechanisms of tumors, and little has been done to investigate the potential regulatory relationship between m6A methylation and TIME, including the unclear role of m6A regulators in AML prognosis. Based on the TCGA dataset, we identified two clusters according to optimal k-means clustering, and observed two clusters with significant differences, suggesting that m6A-related expression regulators are closely related to the prognosis of AML.

The present study showed the expression pattern, prognostic value, and effect of m6A regulatory factors in AML. The expression level of YTHDF3 was lower in the high-risk AML samples than the low-risk AML samples, while the expression level of ALKBH5 did not change significantly between the AML samples. An analysis of the difference between the m6A regulators and PD-L1 in the AML patient samples revealed their close relationship. We identified 2 subtypes of AML through the consistent clustering of m6A regulatory factors (i.e., Cluster 1 and Cluster 2). The Cluster 1 and Cluster 2 subtypes were found to affect the prognosis and different clinicopathological characteristics of AML and to be closely related to PD-L1, immune score, and immune cell infiltration level. We also identified 2 prognostic risk characteristics related to m6A, and divided the AML patients into high-risk and low-risk groups. The high-risk and low-risk groups were also associated with different cluster subtypes and immune scores. The AML prognostic risk assessment model constructed in this study includes treatment-naive, secondary, or treatment-related AML. The AML patients were divided into high-risk and low-risk groups. The survival time of the high-risk group was significantly shorter than that of the low-risk group; the survival time of patients aged ≤65 years was significantly longer. Compared with the prognostic predictive ability of other clinical features, this model showed better predictive performance. In conclusion, the risk signature we constructed may be regarded as a new biomarker with special potential and even more effective in clinical treatment. However, due to the limited number of databases, the model needs to be repeatedly confirmed in more clinical data.

The expression of the m6A demethylase ALKBH5 is regulated by changes in the chromatin state during the development of human AML, and ALKBH5 is necessary to maintain the function of LSCs. The dynamics of the chromatin state are correlated with the expression regulation of m6A-related modifiers, revealing the selective and key role of ALKBH5 in AML, which can be used as a specific therapeutic target for LSC (31). Shen et al. also found that ALKBH5 plays a key role in the occurrence of leukemia and LSC/(leukemia initiating cells) LIC self-renewal/maintenance and highlighted the therapeutic potential of targeting the ALKBH5/m6A axis (32). However, YTHDF3 had not previously been examined. These 2 m6A-related genes were used to construct a prognostic model that may provide a foundation for the future diagnosis and treatment of AML.

The TIME is mainly composed of tumor cells, surrounding immune and inflammatory cells, tumor-related fibroblasts, and nearby interstitial tissues, capillaries, and various cytokines and chemokines. Its synthesis is complex. The system can be divided into an immune microenvironment dominated by immune cells and a non-immune microenvironment dominated by fibroblasts (33). We found that resting memory CD4 T cells, CD8 T cells, M2 macrophages, memory B cells, and eosinophils were positively correlated to the risk scores of patients. Additionally, compared to the patients with high-risk scores, those with low-risk scores had higher immune scores. These results suggest that a comprehensive analysis of m6A is useful for understanding and facilitating the exploration of the relationship between cellular infiltration and disease risk associated with the tumor immune microenvironment, ultimately leading to individual immunotherapy. In this study, it is unclear whether the risk score is influenced by the biological characteristics of AMLs. In other words, what effect does genetic mutation or (World Health Organization/European LeukemiaNet) WHO/ELN disease risk score have on 6A-related modifications? This still needs further research.

It is undeniable that our research has some limitations. First, the extrapolation of our results was based on an internal cohort from TCGA. Due to the lack of sufficient data available in our cohort, our risk-score model was not externally validated. Thus, further external verification is required with a multicenter cohort to test the results. Additionally, it is necessary to conduct further experimental research.


Conclusions

Our study clarified the important role of m6A methylation in the TIME of AML. We also constructed a prognostic model using the m6A-related signatures associated with TIME to predict the survival of the AML patients. Changes in m6A-related signatures associated with TIME not only reflect disease progression, but may also be promising therapeutic targets for improving the efficacy of immunotherapy in AML.


Acknowledgments

The authors appreciate the academic support from the AME Oncology Collaborative Group.

Funding: This work was supported by grants from the Nantong Health Commission Project (No. QA2021031) and the Nantong Science and Technology Project (No. HS2019003).


Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://atm.amegroups.com/article/view/10.21037/atm-22-3858/rc

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://atm.amegroups.com/article/view/10.21037/atm-22-3858/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.


References

  1. Koenig KL, Sahasrabudhe KD, Sigmund AM, et al. AML with Myelodysplasia-Related Changes: Development, Challenges, and Treatment Advances. Genes (Basel) 2020;11:845. [Crossref] [PubMed]
  2. Yi M, Li A, Zhou L, et al. The global burden and attributable risk factor analysis of acute myeloid leukemia in 195 countries and territories from 1990 to 2017: estimates based on the global burden of disease study 2017. J Hematol Oncol 2020;13:72. [Crossref] [PubMed]
  3. Shallis RM, Wang R, Davidoff A, et al. Epidemiology of acute myeloid leukemia: Recent progress and enduring challenges. Blood Rev 2019;36:70-87. [Crossref] [PubMed]
  4. McNerney ME, Godley LA, Le Beau MM. Therapy-related myeloid neoplasms: when genetics and environment collide. Nat Rev Cancer 2017;17:513-27. [Crossref] [PubMed]
  5. Molica M, Breccia M, Foa R, et al. Maintenance therapy in AML: The past, the present and the future. Am J Hematol 2019;94:1254-65. [Crossref] [PubMed]
  6. Short NJ, Rytting ME, Cortes JE. Acute myeloid leukaemia. Lancet 2018;392:593-606. [Crossref] [PubMed]
  7. Böhme M, Kayser S. Immune-Based Therapeutic Strategies for Acute Myeloid Leukemia. Cancers (Basel) 2021;14:105. [Crossref] [PubMed]
  8. Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. Biomed Pharmacother 2019;112:108613. [Crossref] [PubMed]
  9. Chen M, Wong CM. The emerging roles of N6-methyladenosine (m6A) deregulation in liver carcinogenesis. Mol Cancer 2020;19:44. [Crossref] [PubMed]
  10. Deng LJ, Deng WQ, Fan SR, et al. m6A modification: recent advances, anticancer targeted drug discovery and beyond. Mol Cancer 2022;21:52. [Crossref] [PubMed]
  11. Reichel M, Köster T, Staiger D. Marking RNA: m6A writers, readers, and functions in Arabidopsis. J Mol Cell Biol 2019;11:899-910. [Crossref] [PubMed]
  12. Huang Y, Su R, Sheng Y, et al. Small-Molecule Targeting of Oncogenic FTO Demethylase in Acute Myeloid Leukemia. Cancer Cell 2019;35:677-691.e10. [Crossref] [PubMed]
  13. Weng H, Huang H, Wu H, et al. METTL14 Inhibits Hematopoietic Stem/Progenitor Differentiation and Promotes Leukemogenesis via mRNA m6A Modification. Cell Stem Cell. 2018;22:191-205.e9. [Crossref] [PubMed]
  14. Yin T, Zhao L, Yao S. Comprehensive characterization of m6A methylation and its impact on prognosis, genome instability, and tumor microenvironment in hepatocellular carcinoma. BMC Med Genomics 2022;15:53. [Crossref] [PubMed]
  15. Xu-Monette ZY, Zhang M, Li J, et al. PD-1/PD-L1 Blockade: Have We Found the Key to Unleash the Antitumor Immune Response? Front Immunol 2017;8:1597. [Crossref] [PubMed]
  16. Wu Y, Chen W, Xu ZP, et al. PD-L1 Distribution and Perspective for Cancer Immunotherapy-Blockade, Knockdown, or Inhibition. Front Immunol 2019;10:2022. [Crossref] [PubMed]
  17. Shen X, Zhao B. Efficacy of PD-1 or PD-L1 inhibitors and PD-L1 expression status in cancer: meta-analysis. BMJ 2018;362:k3529. [Crossref] [PubMed]
  18. Yang H, Bueso-Ramos C, DiNardo C, et al. Expression of PD-L1, PD-L2, PD-1 and CTLA4 in myelodysplastic syndromes is enhanced by treatment with hypomethylating agents. Leukemia 2014;28:1280-8. [Crossref] [PubMed]
  19. Zhang L, Gajewski TF, Kline J. PD-1/PD-L1 interactions inhibit antitumor immune responses in a murine acute myeloid leukemia model. Blood 2009;114:1545-52. [Crossref] [PubMed]
  20. Yang X, Ma L, Zhang X, et al. Targeting PD-1/PD-L1 pathway in myelodysplastic syndromes and acute myeloid leukemia. Exp Hematol Oncol 2022;11:11. [Crossref] [PubMed]
  21. Wang S, Sun C, Li J, et al. Roles of RNA methylation by means of N6-methyladenosine (m6A) in human cancers. Cancer Lett 2017;408:112-20. [Crossref] [PubMed]
  22. Chen Y, Lin Y, Shu Y, et al. Interaction between N6-methyladenosine (m6A) modification and noncoding RNAs in cancer. Mol Cancer 2020;19:94. [Crossref] [PubMed]
  23. Kwok CT, Marshall AD, Rasko JE, et al. Genetic alterations of m6A regulators predict poorer survival in acute myeloid leukemia. J Hematol Oncol 2017;10:39. [Crossref] [PubMed]
  24. Niu Y, Zhao X, Wu YS, et al. N6-methyl-adenosine (m6A) in RNA: an old modification with a novel epigenetic function. Genomics Proteomics Bioinformatics 2013;11:8-17. [Crossref] [PubMed]
  25. Tang Y, Chen K, Song B, et al. m6A-Atlas: a comprehensive knowledgebase for unraveling the N6-methyladenosine (m6A) epitranscriptome. Nucleic Acids Res 2021;49:D134-43. [Crossref] [PubMed]
  26. Wang T, Kong S, Tao M, et al. The potential role of RNA N6-methyladenosine in Cancer progression. Mol Cancer 2020;19:88. [Crossref] [PubMed]
  27. Li Z, Weng H, Su R, et al. FTO Plays an Oncogenic Role in Acute Myeloid Leukemia as a N6-Methyladenosine RNA Demethylase. Cancer Cell 2017;31:127-41. [Crossref] [PubMed]
  28. Strickland SA, Vey N. Diagnosis and treatment of therapy-related acute myeloid leukemia. Crit Rev Oncol Hematol 2022;171:103607. [Crossref] [PubMed]
  29. Paris J, Morgan M, Campos J, et al. Targeting the RNA m6A Reader YTHDF2 Selectively Compromises Cancer Stem Cells in Acute Myeloid Leukemia. Cell Stem Cell 2019;25:137-148.e6. [Crossref] [PubMed]
  30. Ayyadurai VAS, Deonikar P, McLure KG, et al. Molecular Systems Architecture of Interactome in the Acute Myeloid Leukemia Microenvironment. Cancers (Basel) 2022;14:756. [Crossref] [PubMed]
  31. Wang J, Li Y, Wang P, et al. Leukemogenic Chromatin Alterations Promote AML Leukemia Stem Cells via a KDM4C-ALKBH5-AXL Signaling Axis. Cell Stem Cell 2020;27:81-97.e8. [Crossref] [PubMed]
  32. Shen C, Sheng Y, Zhu AC, et al. RNA Demethylase ALKBH5 Selectively Promotes Tumorigenesis and Cancer Stem Cell Self-Renewal in Acute Myeloid Leukemia. Cell Stem Cell 2020;27:64-80.e9. [Crossref] [PubMed]
  33. Quail DF, Joyce JA. Microenvironmental regulation of tumor progression and metastasis. Nat Med 2013;19:1423-37. [Crossref] [PubMed]

(English Language Editor: L. Huleatt)

Cite this article as: Yuan S, Cong Z, Ji J, Zhu L, Jiang Q, Zhou Y, Shen Q, Damiani D, Xu X, Li B. Analysis of m6A-related signatures associated with the tumor immune microenvironment and predict survival in acute myeloid leukemia. Ann Transl Med 2022;10(16):902. doi: 10.21037/atm-22-3858

Download Citation