Identification of m6A-Related lncRNAs Associated With Prognoses and Immune Responses in Acute Myeloid Leukemia

Background: Acute myeloid leukemia (AML) remains the most common type of hematopoietic malignancy in adults and has an unfavorable outcome. Herein, we aimed to construct an N6-methylandenosine (m6A)-related long noncoding RNAs (lncRNAs) signature to accurately predict the prognosis of patients with AML using the data downloaded from The Cancer Genome Atlas (TCGA) database. Methods: The RNA-seq and clinical data were obtained from the TCGA AML cohort. First, Pearson correlation analysis was performed to identify the m6A-related lncRNAs. Next, univariate Cox regression analysis was used to determine the candidate lncRNAs with prognostic value. Then, feature selection was carried out by Least absolute shrinkage and selection operator (LASSO) analysis, and seven eligible m6A-related lncRNAs were included to construct the prognostic risk signature. Kaplan–Meier and receiver operating characteristic (ROC) curve analyses were performed to evaluate the predictive capacity of the risk signature both in the training and testing datasets. A nomogram was used to predict 1-year, 2-year, and 3-year overall survival (OS) of AML patients. Next, the expression levels of lncRNAs in the signature were validated in AML samples by qRT-PCR. Functional enrichment analyses were carried out to identify probable biological processes and cellular pathways. The ceRNA network was developed to explore the downstream targets and mechanisms of m6A-related lncRNAs in AML. Results: Seven m6A-related lncRNAs were identified as a prognostic signature. The low-risk group hold significantly prolonged OS. The nomogram showed excellent accuracy of the signature for predicting 1-year, 2-year and 3-year OS (AUC = 0.769, 0.820, and 0.800, respectively). Moreover, the risk scores were significantly correlated with enrichment in cancer hallmark- and malignancy-related pathways and immunotherapy response in AML patients. Conclusion: We developed and validated a novel risk signature with m6A-related lncRNAs which could predict prognosis accurately and reflect the immunotherapy response in AML patients.


INTRODUCTION
Acute myeloid leukemia (AML) is one of the most aggressive hematological malignancies, with the highest incidence in adults (Short et al., 2018). Despite advancements in the biological understanding of AML, the standard treatment has not significantly improved over 40 years (Choi et al., 2020). The prognosis remains dismal, and most patients encounter disease recurrence or die of the disease several months after their initial remission. The long-term overall survival (OS) rate of young patients (<60 years) is less than 40%, and that of elderly patients (≥60 years) is only 15% (Koenig and Mims, 2020). Therefore, it is essential to identify novel molecules to accurately predict the prognosis and guide treatment in AML.
RNA methylation remains the frontier and focus of current epigenetic research, and N6-methyladenosine (m6A) is the main internal epigenetic modification in eukaryotic messenger RNAs (mRNAs) and noncoding RNAs (ncRNAs) (Deng et al., 2018). m6A affects the stability, alternative splicing, nuclear exit and translation efficiency of mRNA, and therefore plays important roles in various disease processes, including cancer (Cao et al., 2016;Wei et al., 2017;Dai et al., 2018). Recent research has revealed that the inactivation of m6A RNA methylases and demethylases in AML suppresses malignant cells through multiple m6A-dependent mechanisms. m6A RNA methylases METTL3/14 (Vu et al., 2017;Weng et al., 2018) and demethylases FTO/ALKBH5 (Huang et al., 2019;Shen et al., 2020), which are aberrantly expressed in specific leukemia subtypes, play a critical role in leukemogenesis by regulating specific gene targets and signaling pathways. Likewise, the m6A reader YTHDF2, which is highly expressed across AML subgroups, has been found to promote leukemogenesis by inhibiting gene targets through a YTHDF2-mediated mRNA decay mechanism (Paris et al., 2019). Therefore, m6A mRNA modification has promising therapeutic and prognostic potential in AML. Long noncoding RNAs (lncRNAs), which possess no protein-coding capacity, are defined as transcripts longer than 200 nt and regarded as "noise" in genome transcription (Gibb et al., 2011). Previous studies have revealed that lncRNAs play essential roles in many important regulatory processes, such as transcription regulation, genome imprinting, chromatin modification, and nuclear transport, which have attracted widespread attention (Groff et al., 2016;Long et al., 2017). Some lncRNAs, including CASC15 (Fernando et al., 2017), UCA1 , H19 (Zhang et al., 2018), HOTAIRM1 (Chen et al., 2020), LINC00152 (Cui et al., 2021), have been identified to be associated with the recurring mutations, clinical features and prognosis of AML. Accumulating studies have indicated that aberrant expression of certain specific lncRNAs in tumor cells can be used as a diagnostic marker or a potential drug target. Additionally, lncRNAs can easily be detected in serum and saliva, urine, blood, or tissue biopsy, which makes them attractive for clinical diagnosis and prognostic prediction (Chandra Gupta and Nandan Tripathi, 2017).
Extensive studies have showed that some m6A modifications could be directly or indirectly regulated by lncRNAs. The interaction of m6A modification and lncRNAs plays crucial roles in tumor progression, metastasis, response to immune, drug resistance, and offers new insights for early diagnosis and new treatment strategies of cancer (Chen et al., 2020). However, specific roles of the m6A-related lncRNAs in AML remains to be elucidated. Therefore, understanding how m6A-related lncRNAs work may enable to identify potential molecules as therapeutic targets. In this study, we first constructed a prognostic risk signature with seven m6A-related lncRNAs and further validated the reliability and sensitivity of the signature. Moreover, we also explored the correlation of the risk score and tumor microenvironment. Finally, we built a ceRNA and PPI networks in order to further study the potential mechanisms of m6A-related lncRNAs in AML.

Data Collection
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.

Generation of a Prognostic Signature Using LASSO Regularization
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 m6Arelated 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.

Survival Analysis
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.

Construction and Evaluation of the Nomogram
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.

Functional Enrichment Analysis
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.

CeRNA Network
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.

Sample Collection
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.

qRT-PCR Analysis
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.

Statistical Analysis
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.

Construction of the m6A-Related lncRNAs Prognostic Signature for AML Patients
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 Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 770451 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).

Validation of the Signature Performance in the TCGA Dataset
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 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.

Independent Prognostic and Clinicopathological Correlation Analyses
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-, 2and 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.

Validation of the Expression Levels of Seven m6A-Related lncRNAs in AML Samples
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.

Cellular Biological Effects Related to the Signature
t-distributed stochastic neighborhood embedding (t-SNE) was used to investigate the differences between the low-risk and highrisk 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.

Correlation of the Signature and Immunotherapy Response
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 Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 770451 6 group (p＜0.05) ( Figure 9A), and the risk scores in the noresponse 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.

Construction of the ceRNA Network and Functional Enrichment Analysis
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 Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 770451 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.

DISCUSSION
To the best of our knowledge, AML is a deadly disease associated with poor outcomes despite advancements in targeted molecular and immunotherapy (Short et al., 2018). The m6A is the most abundant epigenetic modification of mRNA and lncRNAs and plays important roles in many biological processes (Wei et al., 2017). The lncRNAs do not possess protein coding ability but participate in a plethora of cellular functions (Ponting et al., 2009;Morlando and Fatica, 2018). The m6A modifications of lncRNAs are associated with the development, occurrence, and prognosis of a variety of tumors and other diseases, including AML (Coker et al., 2019), lncRNAs can also affect tumor invasive progression by targeting m6A regulators as competitive endogenous RNAs (ceRNAs) (Tu et al., 2020). Both m6A and lncRNAs are important regulators of AML (Zimta et al., 2019;Zheng and Gong, 2021). However, the potential roles of m6A-related Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 770451 8 lncRNAs in AML remains unclear. Thus, based on the data from the TCGA dataset, we developed a novel m6A-related lncRNAs signature to accurately predict the prognoses in AML.
A total of 144 AML samples were included to explore the prognostic value of m6A-related lncRNAs. Seven of 15 m6Arelated lncRNAs with prognostic value were used to establish a signature to predict the OS of AML patients. The high-risk group was significantly associated with a poor prognosis. Furthermore, multivariate Cox regression analysis illustrated that the signature acts as an independent prognostic factor in AML, which was consistent with the results of ROC curve analysis.
Among the lncRNAs in the prognostic signature, USP30-AS1 has been reported to be related to autophagy and immunity in bladder cancer (Wan et al., 2021), cervical cancer (Chen et al., 2020) and melanoma (Ding et al., 2021). MIR133A1HG was also reported to be one of the lncRNAs in the autophagy-related  signature for AML (Zhao et al., 2021). However, the underlying molecular biological mechanism has not been elucidated. Notably, AC114271.2, AF064858.8, RP11-22L13.1, RP11-544A12.4 and MIR181A1HG were revealed for the first time.
The correlation analysis revealed that two of seven lncRNAs, MIR133A1HG and RP11-544A12.4, have significant correlation with the expression of almost all m6A regulators, and that MIR181A1HG, AC114271.2 and USP30-AS1 were related to  Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 770451 10 some m6A regulators, including m6A writers, erasers, and readers. Furthermore, AF064858.8 is related to some m6A RNA writers and readers, while RP11-22L13.1 is related to m6A RNA readers IGF2BP2 and IGF2BP3. However, the specific mechanism on how lncRNA regulates m6A remains to be further studied (Supplementary Figure S3). The results of the pan-cancer analysis showed that USP30-AS1, AF064858.8, MIR133A1HG, RP11-544A12.4 and MIR181A1HG were specifically expressed in AML, which indicated that these lncRNAs play a crucial role in its carcinogenesis (Supplementary Figures S4, S5).
GO and KEGG revealed that cancer hallmark-related and malignancy-related pathways were more enriched in the highrisk group. In addition, we found that T cell activation was also enriched, indicating that m6A-related lncRNAs may affect the prognosis of patients through immune mechanisms, which is consistent with the results of the previous studies Zhang et al., 2021). The infiltrating immune cells in the tumor microenvironment (TME), particularly T cells, are key mediators of tumor destruction and play important roles in immunotherapy (Lei et al., 2020). m6A has been reported to regulate the maturation and neoantigen presentation involved in the immunotherapeutic response . Nevertheless, the potential function and prognostic value of m6A-related lncRNAs in mediating immunotherapeutic response and prognosis in AML remain to be characterized. The TIDE prediction score is a successfully validated computational framework for immunotherapy prediction (Jiang et al., 2018). Correspondingly, the TIDE algorithm was also performed to predict the correlation between the risk score and the immunotherapeutic response, the results indicated that risk scores of the no-response group were higher than those of the response group. Obviously, this signature was significantly correlated with immunity and a potential biomarker for predicting the response to immunotherapy in AML. Seven of the 15 m6A-related lncRNAs with prognostic value were used to construct the prognostic signature. This does not necessarily mean that these 7 lncRNAs are more important than other m6A-related lncRNAs; rather, it indicates that the combination of these 7 m6A-related lncRNAs can adequately predict AML prognosis. The results of the ceRNA network and functional enrichment analysis showed that other m6A-related lncRNAs (SENCR, PSMB8-AS1, ITGB2-AS1, RP11-89N17.1, RP11-483I13.5 and LINC01547) also play critical roles during tumor progression by regulating the expression of important genes such as PTEN, VEGFA, MAPK1, IGF1, etc. These genes may provide us with new therapeutic targets for AML patients.
In conclusion, we first constructed a 7 m6A-related lncRNA risk signature to predict the prognosis of AML and verified the predictive reliability and sensitivity of the signature. We also explored the distinct molecular landscape classified by the signature, including biological processes, pathways, correlation with immune therapies, potential targets, and provided new insights into the potential roles and mechanisms of m6Arelated lncRNAs in AML.
However, our study has several limitations. First, we used the data from the public dataset TCGA to construct and validate the signature. There was no suitable external database on AML to assess the reliability of the signature. Additionally, functional analyses and mechanistic studies of the signature were not carried out.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Medical Ethics Committee of The Affiliated Cancer Hospital of Zhengzhou University. The patients/ participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
WZ, XW, and DL conceived and designed the study. DL, JL, CC, WG, SL, WS, ZS, YB, and YZ participated in the coordination of data acquisition and data analysis and reviewed the manuscript.

ACKNOWLEDGMENTS
We are grateful to the contributors to the public databases used in this study.
Supplementary Figure S2 | The most pronounced differential biological processes between the high-and low-risk groups explored by GSEA.
Supplementary Figure S3 | The correlation analysis between seven lncRNAs and 21 m6A genes.
Supplementary Figure S4 | The expressions of 7 signature lncRNAs in pan-cancer from TCGA.
Supplementary Figure S5 | The expressions of 7 signature lncRNAs in AML and normal tissues analyzed by GEPIA.