We obtained the RNA-seq data and relevant clinical profiles of AML patients from The Cancer Genome Atlas (TCGA) database. After screening, cases with missing clinical data and/or OS ≤ 30 days were excluded from the study, and a total of 144 AML cases were included in the analysis. We subsequently transformed the probe IDs of each AML cohort into gene symbols based on the annotation files. The clinical characteristics for AML cases were summarized in Supplementary Table S1.

We retrieved 21 m6A gene expression matrixes from the TCGA AML cohort, including expression data on methyltransferases, binding proteins, and demethylases. The 21 m6A genes included erasers (ALKBH5 and FTO), writers (ZC3H13, KIA1429, RBM15, RBM15B, METTL3, METTL14, METTL16 and WTAP), and readers (RBMX, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, YTHDC1, YTHDC2, YTHDF1, YTHDF2 and YTHDF3). 11,904 lncRNA expression data were extracted from the TCGA AML cohort. Subsequently, Pearson correlation analysis was applied to investigate the relationship between m6A genes and lncRNA expression to determine the m6A-related lncRNAs with a correlation coefficient >0.3 and p < 0.05.

To establish a reliable signature, the 144 cases were randomly divided into training (100 cases) and testing datasets (44 cases). The log-rank test was performed to extract the m6A-related lncRNAs that were closely associated with survival time for development of the prognostic signature. Subsequently, least absolute shrinkage and selection operator (LASSO) was performed to exclude the candidate lncRNAs significantly associated with each other to restrict overfitting with the R package “glmnet.” Then, we identified seven m6A-related lncRNAs to build the signature through multivariate Cox proportional hazards regression analysis. The risk scores were calculated using the following formula based on the included m6A-related lncRNAs. Based on the median risk score, the samples were classified to high- and low-risk groups. R i s k   S c o r e = 0.1906265 ∗ USP 30 − AS 1 + 0.10995083 ∗ AC 114271.2 + 0.0704641 ∗ AF 064858.8 + 0.02511635 ∗ RP 11 − 22 L 13.1 + ( − 0.129286 ) ∗ MIR 181 A 1 HG + ( − 0.0512868 ) ∗ RP 11 − 544 A 12.4 + ( − 0.0383528 ) ∗ MIR 133 A 1 HG

Kaplan–Meier curves analysis and log-rank tests were used to compare the discrepancy of OS between the predicted high-risk and low-risk groups. p ≤ 0.05 was defined as statistical difference. All survival analyses and log-rank tests were carried out using the R package survival, while the R package “surviminer” was used to plot the Kaplan–Meier curve.

R package “rms” was performed to construct a nomogram based on clinical stage, T stage and the risk score using the TCGA AML cohort. To evaluate the utility of the nomogram, the R package “ROC survival” was used to construct ROCs for the prediction of the 1-, 2- and 3-year OS by the nomogram. The R package “ggDCA” was used to construct a decision analysis curve to evaluate the clinical utility. Finally, package “rms” was used to establish a calibration curve to assess the precision for prediction of 1-, 2- and 3-year OS.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the differentially expressed genes (DEGs) in the two risk groups were performed by “cluster Pofolier” package. (Yu et al., 2012). Further, Gene Set Enrichment Analysis (GSEA) was used to analyze the biological functions and pathways associated with the high and low-risk groups. FDR <0.05 and p < 0.05 were regarded statistically significant.

MiRcode was used to predict the miRNAs that interact with m6A-related lncRNAs. A total of 25 pairs of interactions between 9 lncRNAs and 10 miRNAs were identified by the miRcode database. A total of 271 mRNAs interacting with the miRNAs were identified to construct a PPI network through TargetScan, miRDB, and miR TarBase. The ceRNA network was visualized using Cytoscape software. The Cyto Hubba plugin was performed to obtain hub genes from the PPI network. GO and KEGG enrichment analyses were conducted to identify biological processes and potential signaling pathways.

We totally collected 21 bone marrow samples, including 14 AML samples (derived from 7 AML patients at the time of first diagnosis and at first relapse) and 7 healthy controls in the Hematological Department of The Affiliated Cancer Hospital of Zhengzhou University. This research was approved by the Medical Ethics Committee of The Affiliated Cancer Hospital of Zhengzhou University (approval no. 2020239). Informed consent and approval were provided by all participants.

Total RNA was isolated from 21 patients’ samples. cDNA synthesis was conducted with a reverse transcription kit (TransGen Biotech, #AU311-02). Then, real-time PCR was performed on the ABI 7500 Fast System (Applied Biosystems, United States) with TB Green Premix Ex Taq (Takara Bio, #RR420A). Relative expression of lncRNAs were normalized to GAPDH and calculated by 2-ΔΔCt method. Primers sequences are listed in Supplementary Table S2.

All statistical data were analyzed using R version 4.0.2. Pearson’s correlation analyses were applied to evaluate the correlation between the risk score and specific gene expression. Kaplan–Meier curves and log-rank tests were conducted for survival analysis. Univariate and multivariate Cox regression were applied to assess the prognostic independence. The reliability and sensitivity of the signature were evaluated using ROC curve analysis. Student’s t-test was used to determine significance between two groups. p < 0.05 was defined as statistical significance.

The flowchart summarized the construction and subsequent analysis of the risk signature (Figure 1). The expression profiles of 21 well known m6A genes and 11,904 lncRNAs were extracted from TCGA AML cohort. The m6A-related lncRNAs were defined as ones which were associated with one or more of these m6A genes (|Pearson R| > 0.3 and p < 0.05), and a total of 3,812 lncRNAs were included (Figure 2A). Next, 15 m6A-related lncRNAs significantly related to the OS of AML patients were screened out through univariate Cox regression analysis (p < 0.001). Then, feature selection was performed using LASSO analysis, and 7 m6A-related lncRNAs were ultimately identified to develop the prognostic risk signature (Figures 2B–D). Among them, USP30-AS1, AC114271.2, AF064858.8 and RP11-22L13.1 are detrimental factors with a hazard ratio (HR) > 1, whereas MIR181A1HG, RP11-544A12.4 and MIR133A1HG are protective factors with a HR < 1 (Figure 2E). Sankey diagram shows the relationship between the lncRNAs and m6A genes and the risk types (Figure 2F).

Based on the median risk score, the samples were classified to high- and low-risk groups. The scatter plot showed that mortality increased with a higher risk score (Figure 3A). Furthermore, the heatmap showed that the expression level of AF064858.8, RP11-22L13.1, USP30-AS1 and AC114271.2 were higher in the high-risk group, whereas MIR181A1HG, RP11-544A12.4 and MIR133A1HG were at lower expressed levels (Figure 3B). Moreover, Kaplan–Meier curve analyses indicated that the low-risk group had prolonged OS (Figure 3C). In addition, the AUC values for the 1-year, 2-year, and 3-year OS were 0.74, 0.765 and 0.702, respectively (Figure 3D), suggesting the high predictive capacity of the signature in the training dataset.

To further validate the predictive capacity of the signature, we used the same algorithm to calculate the risk scores both in the testing dataset and the overall dataset. The risk score distribution plot, scatter plot, and heatmaps were consistent with those in the training dataset. Furthermore, Kaplan–Meier curve analyses showed the consistent results in the testing dataset and overall dataset (p < 0.001). In addition, the time-ROC curves and their AUC values also displayed good performance for predicting prognosis. The AUCs for 1-year, 2-year and 3-year OS in the testing dataset were 0.798, 0.813 and 0.784, and in the overall dataset, the AUCs were 0.748, 0.778 and 0.719, respectively (Figure 4). In summary, the high-risk score based on the signature could indicate a poor prognosis, accurately.

A stratified analysis was carried out according to the clinicopathological characteristics, including age (<60/≥60 years, Figures 5A,B), sex (female/male, Figures 5C,D), abnormalities (normal/abnormal, Figures 5E,F), chromosomal FLT3 mutations (no/yes, Figures 5G,H), IDH1 mutations (no/yes, Figures 5I,J), NPM1 mutations (no/yes, Figure 5K,L), and RAS mutations (no/yes, Figures 5M,N). Kaplan–Meier curve analyses revealed that high-risk group had worse survival outcome than low-risk group when stratified by the different clinical features, except for IDH1, NPM1 or RAS mutations.

Univariate and stepwise multivariate Cox regression analyses were applied to explore whether the m6A-related lncRNAs signature and clinical characteristics, such as age, sex, cytogenic abnormalities, FLT3 mutation, IDH1 mutation, NPM1 mutation and RAS mutation, may serve as independent prognostic factors. Finally, age and risk score were selected as independent prognostic factors for the survival prediction (Figures 6A,B). Next, to quantitatively predict the survival probability of each case, a prognostic nomogram incorporating the clinical features and the signature was plotted (Figure 6C). Furthermore, calibration curves of 1-year, 2-year, and 3-year OS were plotted to confirm the predictive probability of the nomogram. The results showed that the predicted survival probability by the nomogram was consistent to the actual one (Figures 6D–F). Moreover, time-dependent ROC curves revealed that the nomogram showed remarkable accuracy to predict 1-, 2- and 3-year OS (AUC = 0.769, 0.82, and 0.8, respectively) (Figure 6G). In addition, the decision curve analysis of the LNM nomogram showed that even if the threshold probability of the patient is very small, the use of the LNM nomogram in predicting LNM brings more benefit than treating either all or no patients (Figure 6H). Above all, this signature displayed good performance in predicting the survival of AML.

Besides, Kaplan-Meier curve analyses were used to evaluate the prognostic role of each gene from the signature, and the results showed that higher expression levels of USP30-AS1, AC114271.2, AF064858.8 and RP11-22L13.1, and lower expression levels of MIR181A1HG, RP11-544A12.4 and MIR133A1HG were associated with poorer survival outcomes (Supplementary Figure S1).

To further demonstrate the feasibility of the prognostic signature and measure the potential in clinical practice, we performed qRT-PCR assays in our collected bone marrow samples. Our results showed that seven of the m6A-related lncRNAs could be easily detected both in AML patients and healthy controls. And, the expression levels of these lncRNAs were relatively higher in AML patients than healthy controls. It is worthy of note, compared with samples at first diagnosis, the detrimental factors of USP30-AS1, AC114271.2, AF064858.8 and RP11–22L13.1 were upregulated, and the protective factors of MIR181A1HG, RP11–544A12.4 and MIR133A1HG exhibited a decreased tendency in first relapsed AML patients (Figure 7). Generally speaking, the relapsed AML patients always encounter therapy resistance, resulting in poor survival outcomes, which consistent with our experimental validation in AML samples.

t-distributed stochastic neighborhood embedding (t-SNE) was used to investigate the differences between the low-risk and high-risk groups. The results obtained based on entire genes, 21 m6A genes and 7 m6A-related lncRNAs showed that the distributions were relatively scattered (Figures 8A–C). However, the result obtained based on the signature showed that the two risk groups have different distributions, which suggested that the prognostic signature can easily distinguish the low-risk and high-risk groups (Figure 8D).

In view of the excellent predictive capacity of this signature, we further investigated the biological effect related to the molecular heterogeneity. We first identified 1,186 DEGs between the two risk groups. And the DEGs were primarily enriched in the following terms, including cellular defense response, mononuclear cell differentiation, chemokine-mediated signaling pathway and positive regulation of T cell activation (GO Biological Processes) (Figure 8E). Antigen processing and presentation, cytokine-cytokine receptor interaction, intestinal immune network for IgA production as well as Th1 and Th2 cell differentiation (KEGG Pathway) (Figure 8F). GSEA showed several tumor hallmarks were enriched in the high-risk group, such as interferon γ response, HEME metabolism, allograft rejection, interferon α response, complement, myogenesis, inflammatory response, KRAS signaling up, TNFα signaling via NF-κB and so on (Supplementary Figure S2). Above all, these evidences indicated that the signature is correlated with immune cell-related biological pathways and may be associated with the immune microenvironment in AML.

We subsequently evaluated the immunotherapy response in patients with different risk scores. The analysis indicated that the low-risk group had a higher response rate than the high-risk group (p<0.05) (Figure 9A), and the risk scores in the no-response group were higher than those in the response group (Figure 9B). These results suggested that the signature might serve as a good tool to evaluate the immunotherapy response in AML.

LncRNAs can regulate the expression of downstream mRNAs by combining shared miRNAs as ceRNAs (Tay et al., 2014). Therefore, a ceRNA network was constructed to view the potential roles of the m6A-related lncRNAs in AML. A total of 9 lncRNAs, 10 microRNAs and 271 mRNAs were used to construct the network (Figure 10A). Two significant modules, dominated by FOS and NOTCH1 nodes, were identified from the PPI network, and 10 hub genes were extracted (Figure 10B). Moreover, the functional enrichment analysis with 271 target mRNAs showed that the target genes were enriched in regulation of protein serine/threonine kinase activity, protein autophosphorylation, positive regulation of catabolic process, response to transforming growth factor β and epithelial cell proliferation and migration (GO Biological Processes, Figure 10C). MAPK, PI3K-Akt, Rap1, Ras, TGF-β, Estrogen, FoxO signaling pathway, and microRNA in cancer (KEGG Pathway, Figure 10D). The results showed other m6A-related lncRNAs also play critical roles during tumor progression, and provide new insights to study the potential roles and mechanism of m6A-related lncRNAs in AML.