Expression of m6A-regulatory genes in AML, and identification of related lncRNAs. The m6A regulatory genes are important regulators in tumorigenesis and progression of AML. In our present study, we performed a differential expression analysis using AML and normal samples. Compared with that of normal samples, AML samples showed significantly higher expression of METTL3, METL14, METTL16, ZC3H13, RBM15, RBM15B, YTHDC1, YTHDC2 YTHDF1, YTHDF2, YTHDF3, HNRNPC, LRPPRC, HNRNPA2B1, RBMX, and FTO. the expression of WTAP, VIRMA, FMR1, IGFBP1, IGFBP2, and IGFBP3 in normal samples was significantly higher than that in AML samples (P<0.05) (Fig.1A). We then extracted the expression profiles of 14,086 lncRNAs and 23 m6A regulatory genes from AML samples obtained from the TCGA database. Using correlation coefficient to evaluate the relationship between m6A regulatory genes and lncRNAs, we show that there were 525 m6A-regulated lncRNAs and 680 interactions (Fig.1B) (absolute correlation coefficient>0.5, P<0.001).
Prognostic analyses of m6A-related lncRNAs. Using univariate Cox proportional hazard regression analysis, 15 lncRNAs, shown significant in AML prognosis, were screened out from the abovementioned 525 m6A-regulated lncRNAs (P < 0.001, Fig. 2A). Among these 15 lncRNAs, AP003498.2 was a survival risk factor in AML patients, while the remaining lncRNAs were protection factors. We then analyzed the expression levels of these m6A-related prognostic lncRNAs. Our results show that the expression levels of AC025430.1, AFF2-IT1, LINC02593, AC000120.2, AP003498.2, AL158163.1, AC048382.1, AL391834, AC008770.3, and AL133492.1 in AML samples were significantly higher than those in normal samples. The expression levels of AC020916.2 and AJ239328.1 in AML samples were significantly lower than those in normal samples (Fig. 2B). Heatmap shows single-sample expression of m6A-related prognostic lncRNAs in AML and normal-tissue samples (Fig. 2C).
m6A-related lncRNAs are associated with clinicopathologic characteristics of AML. To explore the relationship between m6A-related lncRNAs and clinicopathological characteristics of AML, we used unsupervised cluster analysis to analyze 151 AML samples based on TCGA RNA-seq data of 15 m6A-related prognostic lncRNAs. The 151 AML samples were divided into two clusters having the highest stabilities (Fig. 3A–3D). First, we evaluated the two clusters using survival analysis. Our results indicate that patients in cluster1 showed significantly worse OS that that of patients in cluster2 (P < 0.001, Fig. 3E). Then, we investigated whether there were differences in clinicopathologic factors (i.e., gender, age, and FAB classification), as well as in the expression of the immune checkpoint molecule PD-L1, between the two clusters, and analyzed correlations between the expression levels of PD-L1 and m6A-related lncRNAs. Cluster2 had more elderly patients, and the expression levels of m6A-related prognostic lncRNAs were generally higher in Cluster2 than in Cluster1 (Fig. 3F). PD-L1 showed differential expression in the two clusters, as well as in AML and normal samples (Fig. 3G, 3H). Gene correlation analysis shows that PD-L1 expression was negatively correlated with LINC02593 and AC020916.2 co-expression, and positively correlated with AP003498.2 co-expression (Fig. 3I). AP003498.2 was shown a high-risk lncRNA in AML prognosis. High expression of AP003498.2 was associated with increased risk of death in patients with AML. High expression of PD-L1 may promote the immune escape of AML cells. These results indicate that m6A-related lncRNAs affected the clinicopathologic characteristics of AML and were related to the age and survival of AML patients. m6A-related lncRNAs also affected the expression of PD-L1.
m6A-related lncRNAs are associated with biological characteristics of AML. Next, we analyzed differences in the biological responses of the subgroups generated using consensus clustering in order to further explore the relationship between m6A-related lncRNAs and AML biological characteristics, and impact of m6A-related lncRNAs on AML survival, tumorigenesis, and progression. GSEA was used to explore the main KEGG signaling pathways in the two subgroups. Our results show that cluster1 was mainly involved in toll-like receptor signaling, NOD-like receptor signaling, and B-cell receptor signaling pathways (Fig. 4A–4C), which are related to immunomodulation. Cluster2 was mainly involved in the metabolism of various substances, such as histidine metabolism, and heparan sulfate and glycosylphosphatidylinositol (GPI) anchored protein biosynthesis (Fig. 4D–4F), which are important for adhesion, proliferation, invasion, and metastasis of cancer cells. These results indicate that m6A-related lncRNAs were involved in the biological regulation of AML, and affected tumorigenesis and progression of AML, and survival of AML patients.
m6A-related lncRNAs affect the AML tumor immune microenvironment. To evaluate the characteristics of tumor immune microenvironment in the two clusters, we analyzed the levels of immune cell infiltration and tumor purity in the clusters. The CIBERSORT algorithm was used to calculate infiltration ratio scores for 22 different immune cells in 151 AML samples. The median scores for different types of immune cells in the two clusters were analyzed and compared (Fig. 5A). The ESTIMATE algorithm was used to score the ratio of immune and stromal cells in the tumor microenvironment of patients in cluster1 and cluster2. Our results indicate that infiltration levels of monocytes and M2 macrophages in cluster1 were significantly higher than those in cluster2 (Fig. 5B-5C). The infiltration levels of naïve B cells, plasma cells, resting natural killer (NK) cells, and activated mast cells were significantly higher in cluster2 than in cluster1 (Fig. 5D-5G). In addition, the immune score, stromal score, and ESTIMATE score of cluster1 were higher than those of cluster2 (Fig. 5H–5J). These results indicate that the levels of immune cell infiltration differed considerably between the two groups. Additionally, cluster1, with a high ESTIMATE score, had a lower tumor purity than that of cluster2, indicating that m6A-related lncRNAs affected the tumor immune microenvironment in AML.
Risk signaling of m6A-related lncRNAs and risk prediction model show prognostic value in AML. Using LASSO regression analysis, we screened out the six most representative combinations of m6A-related lncRNAs from the 15 m6A-related prognostic lncRNAs, and established a risk prediction model to evaluate the prognostic value of m6A-related lncRNAs in AML (Fig. 6A- 6B). To ensure the accuracy of the model, 140 AML patients were randomly assigned to a training set (n = 72) and a test set (n = 68) for the construction and verification of the model. After calculating the risk scores of individual patients in the training set, patient samples were divided into a high-risk group and a low-risk group based on median risk score (Fig. 6C). Survival curves show that the number of patient deaths increased as the risk score increased, indicating that the risk score was related to the survival status of AML patients (Fig. 6D). Heatmap shows changes in the risk score and expression levels of lncRNAs. The expression levels of lncRNAs used in the construction of the risk prediction model in the high-risk group were generally lower than those in the low-risk group (Fig. 6E). Analysis of patient survival indicates that the OS of the low-risk group was higher than that of the high-risk group (P < 0.001, Fig. 6F). Area under the curve (AUC) of the ROC curve was 0.852, indicating that the model showed high accuracy in predicting the prognosis of AML patients (Fig. 6G). Next, we used the model to calculate risk scores for individual patients in the test set in order to verify the risk prediction model established using the training set. Our results show that the characteristics of the test set were consistent with those of the training set (Fig. 7A–7E), indicating that m6A-related lncRNAs had prognostic value, and that the risk prediction model showed satisfactory performance in predicting the prognosis of patients with AML.
Independent predictive effect of AML prognosis in the risk prediction model. To further verify whether the risk prediction model was applicable in patients with varying clinicopathological factors, we used the univariate and multivariate independent prognosis analyses in the training and test sets. Both analyses show that risk score was a prognostic factor independent of clinicopathologic factors (P < 0.001, Fig. 8A–8D). Subsequently, patients were stratified according to their clinical pathologic factors including gender, age, and FAB classification, and their risk scores were calculated. Median score was used to distinguish the high-risk and low-risk groups. Survival analysis shows that in each stratified subgroup, patients in the low-risk group showed higher OS than that of patients in the high-risk group (Fig. 9A-9F). These results indicate that the m6A-related lncRNA risk prediction model could predict the prognosis of AML patients without being confounded by gender, age, and FAB classification.
Risk score is correlated with clinicopathological factors and immune infiltration levels. To evaluate the characteristics of risk signaling in m6A-related lncRNA, we analyzed the relationship between risk scores, clinicopathological factors, and immune infiltration levels. Patients were grouped according to their age, gender, clustering subgroups of m6A-related prognostic lncRNAs, and scores obtained for the immune cell ratio, in order to study risk correlation with respect to the grouping of each clinical factor. Our results show that risk scores in the elderly group (≥ 60 years old) were significantly increased compared with those of the young group (< 60 years old, P < 0.05, Fig. 10A). No significant difference in risk scores was found between the genders and between the FAB classifications (Fig. 10B-10C). In the clustering subgroups of m6A-related prognostic lncRNAs, the risk score of cluster1 was significantly higher than that of cluster2 (P < 0.001, Fig. 10D). Comparison of groups with high and low immune scores, categorized based on median immune score, showed that the risk score of the group with high immune scores was significantly higher than that of the group with low immune scores (P < 0.001, Fig. 10E). Heatmap confirmed the above results, and shows that the expression levels of the six lncRNAs used in model construction were low in the high-risk group (Fig. 10F). Further analysis of the relationship between various types of immune cells and risk scores shows that risk scores were significantly related to naïve B cells (R = − 0.3, P = 0.00065), activated dendritic cells (R = − 0.23, P = 0.01), M2 macrophages (R = 0.19, P = 0.036), resting mast cells (R = − 0.24, P = 0.0063), monocytes (R = 0.47, P = 4.8e-08), resting NK cells (R = − 0.18, P = 0.042), plasma cells (R = − 0.27, P = 0.0026), resting memory CD4 T cells (R = − 0.37, P = 2.6e-05), and regulatory T cells (Tregs) (R = 0.21, P = 0.018), and were positively correlated with M2 macrophages, monocytes, and Tregs (Fig. 11A-11I).
Expression of TRAF3IP2-AS1 is regulated by SRSF10. TRAF3IP2-AS1 was selected as a candidate lncRNA to explore its regulatory genes. We first made predictions on the website ENCORI (http://starbase.sysu.edu.cn/)39. These showed that the SRSF10 protein has binding regions for nine transcripts of TRAF3IP2-AS1 in chr6: 111821214–111821244[+] (GSE71096). Compared with normal samples from GTEx databases, the levels of expression of SRSF10 and TRAF3IP2-AS1 were significantly upregulated in AML samples from TCGA databases (Fig. 12A). The expression level of SRSF10 also positively correlated with TRAF3IP2-AS1 expression in TCGA AML samples (Fig. 12B). Next, we verified that the levels of expression of SRSF10 and TRAF3IP2-AS1 were also upregulated in AML samples by RT-PCR (Fig. 12C), and the GEO chip data (GSE65263, GSE114868) showed a similar trend(Fig. 12D-12E). Next, we tested the ability of SRSF10 to regulate the expression of TRAF3IP2-AS1 by siRNA targeting SRSF10 in THP-1 cells (Fig. 12F). RT-PCR showed that the expression of TRAF3IP2-AS1 was significantly downregulated after SRSF10 knockdown (Fig. 12G). These results indicated that the expression of TRAF3IP2-AS1 is regulated by SRSF10 in AML.