Prognostic Value of Genomic Instability of m6A-Related lncRNAs in Lung Adenocarcinoma

Background: Genomic instability of N6-methyladenosine (m⁶A)–related long noncoding RNAs (lncRNAs) plays a pivotal role in the tumorigenesis of lung adenocarcinoma (LUAD). Our study identified a signature of genomic instability of m⁶A-associated lncRNA signature and revealed its prognostic role in LUAD.

Methods: We downloaded RNA-sequencing data and somatic mutation data for LUAD from The Cancer Genome Atlas (TCGA) and the GSE102287 dataset from the Gene Expression Omnibus (GEO) database. The “Limma” R package was used to identify a network of regulatory m⁶A-related lncRNAs. We used the Wilcoxon test method to identify genomic-instability–derived m⁶A-related lncRNAs. A competing endogenous RNA (ceRNA) network was constructed to identify the mechanism of the genomic instability of m⁶A-related lncRNAs. Univariate and multivariate Cox regression analyses were performed to construct a prognostic model for internal testing and validation of the prognostic m⁶A-related lncRNAs using the GEO dataset. Performance analysis was conducted to compare our prognostic model with the previously published lncRNA models. The CIBERSORT algorithm was used to explore the relationship of m⁶A-related lncRNAs and the immune microenvironment. Prognostic m⁶A-related lncRNAs in prognosis, the tumor microenvironment, stemness scores, and anticancer drug sensitivity were analyzed to explore the role of prognostic m⁶A-related lncRNAs in LUAD.

Results: A total of 42 genomic instability–derived m⁶A-related lncRNAs were differentially expressed between the GS (genomic stable) and GU (genomic unstable) groups of LUAD patients. Four differentially expressed lncRNAs, 17 differentially expressed microRNAs, and 75 differentially expressed mRNAs were involved in the genomic-instability–derived m⁶A-related lncRNA-mediated ceRNA network. A prediction model based on 17 prognostic m⁶A-associated lncRNAs was constructed based on three TCGA datasets (all, training, and testing) and validated in the GSE102287 dataset. Performance comparison analysis showed that our prediction model (area under the curve [AUC] = 0.746) could better predict the survival of LUAD patients than the previously published lncRNA models (AUC = 0.577, AUC = 0.681). Prognostic m⁶A-related-lncRNAs have pivotal roles in the tumor microenvironment, stemness scores, and anticancer drug sensitivity of LUAD.

Conclusion: A signature of genomic instability of m⁶A-associated lncRNAs to predict the survival of LUAD patients was validated. The prognostic, immune microenvironment and anticancer drug sensitivity analysis shed new light on the potential novel therapeutic targets in LUAD.

Introduction

Lung cancer is the most frequently diagnosed cancer and the leading cause of cancer-related deaths globally (Torre et al., 2015; Bray et al., 2018). Lung adenocarcinoma (LUAD), a type of non–small cell lung cancer (NSCLC), accounts for more than 40% of lung cancer (Jordan et al., 2017). Despite great progress in drugs, including tyrosine kinase inhibitors and immune checkpoint inhibitors, the 5-year survival rate of patients with lung adenocarcinoma was less than 15% (Zhao et al., 2018; Zhang et al., 2021). Therefore, there is an urgent need to identify new prognostic molecular biomarkers to predict survival and serve as new therapeutic targets in LUAD patients.

Genomic instability, a hallmark of cancer, results from mutations in DNA repair genes and promotes cancer development (Negrini et al., 2010; Andor et al., 2017; Seton-Rogers, 2018; Duijf et al., 2019; Hypoxic Tumors Share Geno, 2019). Recent studies have demonstrated that the characteristics of genomic instability in cancers are associated with clinical implications and prognosis (Sheffer et al., 2009; Sahin et al., 2016). Research demonstrated that human long noncoding RNA (lncRNA) LINC00657 that is induced after DNA damage maintains genomic instability by sequestering PUMILIO proteins (Lee et al., 2016). Genomic instability–related lncRNAs have a critical role in the tumorigenesis of cancers (Chua et al., 2017; Munschauer et al., 2018; Tracy et al., 2018; Petti et al., 2019; Bao et al., 2020). Recently, a study demonstrated that EZH2 mediates ribosomal DNA stability via silencing of lncRNA PHACTR2-AS1, which promotes breast cancer cell proliferation and metastasis (Chu et al., 2020). N⁶-methylandenosine (m⁶A)–related lncRNAs play a significant role in multiple cancers (Zhang et al., 2019a; Ni et al., 2019; Wang et al., 2019), and an m⁶A-related lncRNA was developed to predict the diagnosis and prognosis of cancers (Tu et al., 2020). Recently, a study developed mutator-derived lncRNA signatures of genome instability for predicting clinical outcomes in breast cancer (Bao et al., 2020). Although some lncRNAs have been shown to be involved in genomic instability, the specific role of genomic instability–associated m⁶A-related lncRNAs and their clinical implications in LUAD remains largely unexplored.

In this study, we identified differentially expressed m⁶A-related lncRNAs associated with genomic instability and constructed a competing endogenous RNA (ceRNA) network and then combined m⁶A-related lncRNA expression and somatic mutation profiles based on the tumor genome to reveal a prognostic m⁶A-associated lncRNA signature for predicting the survival of LUAD patients. We also further elucidated the pivotal role of prognostic m⁶A-associated lncRNA in LUAD and validated the expression of prognostic m⁶A-related lncRNAs in LUAD cells and normal bronchial epithelioid cells; meanwhile, prognostic m⁶A-related lncRNAs of the tumor microenvironment, stemness scores, anticancer drug sensitivity, and immune subtype in LUAD were also explored, providing novel therapeutic targets based on RNA modification in LUAD.

Materials and Methods

Data Acquisition

RNA sequencing (RNA-seq) transcription expression data, clinical characteristic data, and simple nucleotide variation (Masked Somatic Mutation from VarScan2 Variant Aggregation and Masking) data of lung adenocarcinoma (LUAD) patients were obtained from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). Correlation analysis was conducted to identify the m⁶A-associated lncRNAs in LUAD according to the following criteria: absolute value of Pearson’s coefficient >0.4, p-value < 0.001. The Gene Expression Omnibus (GEO) dataset GSE102287 was downloaded and used to validate the prognostic mutator-derived m⁶A-associated lncRNAs related to genomic instability.

Identification of Genomic Instability–Associated m⁶A-Related lncRNAs

We used the Perl language to calculate mutation counts in LUAD patients. Patients were classified into two groups comprising the top 25% (genomic unstable [GU] group) and the bottom 25% (genomic stable [GS] group) with respect to mutation frequency. The Wilcoxon test method was used to perform the differential expression analysis of genomic instability–associated m⁶A-related lncRNAs between the GS and GU groups. Furthermore, we explored the differential expression levels of m⁶A-associated genes and the differential somatic mutation counts between the GS and GU groups using the “Limma” and “Ggpubr” R packages. Euclidean distances and Ward’s linkage method were used for hierarchical cluster analyses.

Construction of ceRNA Network Mediated by Mutator-Derived m⁶A-Related lncRNAs

To further explore the mechanisms of mutator-derived m⁶A-related lncRNAs, we constructed a mutator-derived, m⁶A-related, lncRNA-mediated, ceRNA network. First, we used the Wilcoxon test method to perform differential expression analysis of mutator-derived m⁶A-related lncRNAs between the GS and GU groups with thresholds of |log fold change (FC)| > 1 and p-value < 0.05. The MiRcode database was used to determine the relationships between microRNAs (miRNAs) and lncRNAs. The miRTarBase, MiRDB, and TargetScan databases were used to determine the relationships between miRNAs and mRNAs. Finally, we obtained a mutator-derived, m⁶A-related, lncRNA-mediated, ceRNA network. We used Cytoscape 3.6.1 to visualize the genomic instability of the lncRNA–miRNA–mRNA–ceRNA network.

7Functional Enrichment and ConsensusPathDB Pathway Analysis

To further explore the functions of genome instability–related m⁶A-associated lncRNAs, we conducted Gene Ontology (GO) and ConsensusPathDB analysis. We used the Database for Annotation, Visualization, and Integrated Discovery to determine the functional enrichment of mutator-derived m⁶A-associated lncRNAs. Then, we used the “prepareplot.pl” package to visualize the biological processes and molecular functions of genome instability–related m⁶A-associated lncRNAs. The ConsensusPathDB human database (http://cpdb.molgen.mpg.de/) integrates interaction networks in Homo sapiens and currently incorporates genetic, signaling, and gene-regulatory interactions from 32 public resources (Kamburov et al., 2013). We used ConsensusPathDB to analyze the pathway enrichment of mutator-derived m⁶A-related lncRNAs.

Construction and Validation of a Risk Prediction Model of Mutator-Derived m⁶A-Related lncRNAs

We used the “survival,” “caret,” “glmnet,” “survminer,” and “timeROC” packages to construct the prediction model of m⁶A-related lncRNAs. First, TCGA data grouping is cycled 1 time, and the groups were determined with significant training and testing at a ratio of 7:3. Then, we conducted univariate Cox regression analysis to identify the significant m⁶A-associated lncRNAs using the criterion of Cox p-value less than 0.05. Finally, multivariate Cox regression analysis was performed to establish the Cox prediction model. Receiver operating characteristic (ROC) curves were constructed to validate the effects of the m⁶A-related lncRNA signature on clinical outcomes. The area under the ROC curve (AUC) values were calculated to assess the ability of the prognostic m⁶A-associated lncRNAs to predict the survival of LUAD patients. We used time-dependent ROC curves to assess the performance of the m⁶A-related lncRNA signature. All statistical analyses were performed using R version 4.0.3.

Expression, Clinical Characteristics, Immune Microenvironment, and Anticancer Drug Sensitivity Analyses of Prognostic m⁶A-Associated lncRNAs

To investigate the role of prognostic m⁶A-associated lncRNAs in LUAD, we performed correlation analyses involving expression levels of m⁶A-related lncRNAs, LUAD immune environment, stemness scores based on RNA methylation and DNA expression (RNAss and, DNAss, respectively), and drug sensitivity. We downloaded LUAD transcription expression data, RNAss and DNAss data, and TCGA LUAD phenotype-immune subtype data from UCSC-Xena (https://xenabrowser.net/). Then, we combined the expression profiles of 17 prognostic m⁶A-related lncRNAs in LUAD with the LUAD phenotype-immune subtype data, transcription expression data, and stemness scores to perform Spearman’s correlation analysis for immune subtype and stemness score. Using RNA-seq data and compound activity (DTP NCI-60) data from the CellMiner database, combined with FDA (Food and Drug Administration)-approved anticancer drugs and expression of 16 prognostic m⁶A-related lncRNAs, we performed Pearson’s correlation analysis of anticancer drug sensitivity in LUAD. All statistical analyses used R software version 4.0.3.

Validation of Expression of 17 Prognostic m⁶A-Associated lncRNAs

To further explore the expression of 17 prognostic m⁶A-associated lncRNAs between LUAD and non-LUAD tissues, we used the paired sample t-test to perform differential expression analysis in 57 paired LUAD and adjacent non-LUAD tissues from TCGA database. Meanwhile, we used the Wilcoxon rank sum test to perform the differential expression analysis in 535 LUAD tissues and 59 adjacent non-LUAD tissues from TCGA database.

Immune Microenvironment Analysis of Prognostic m⁶A-Associated lncRNAs

To verify the relationship of the m⁶A-associated lncRNAs and immune microenvironment, we used the CIBERSORT (Cell Type Identification by Estimating Relative Subsets of RNA Transcripts) algorithm to obtain the fraction of 22 immune cell types from six microarray public datasets based on the expression file. The Perm was set as 1,000. According to the median expression of 17 prognostic m⁶A-associated lncRNAs, we divided TCGA-LUAD patients into the high-expression group and low-expression group. A p-value < 0.05 was selected for further analysis in CIBERSORT results. The Mann–Whitney U test was used to compare differences in immune subtypes in high-expression and low-expression groups.

Cell Culture

Beas-2B (human bronchial epithelioid cells) was purchased from Saibai Kang Biotechnology Co., Ltd. Human LUAD cells A549, H1299, and H1975 were purchased from Procell Life Sciences& Technology Co., Ltd. We used the short tandem repeat (STR) method to identify these cells. Beas-2B was cultivated in DMEM (PM150210) with 10% fetal bovine serum (Lonsera, Shanghai Shuangru Biotechnology Co., Ltd.) and penicillin/streptomycin (100 mg/ml), and A549, H1299, and H1975 were cultivated in RPMI-1640 medium (Gibco, Invitrogen, Carlsbad CA, United States ) with 10% fetal bovine serum (Lonsera, Shanghai Shuangru Biotechnology Co., Ltd.) and penicillin/streptomycin (100 mg/ml). All cells were cultured at 37°C in a sterile humid incubator containing 5% carbon dioxide.

RNA Isolation, Reverse Transcription, and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

Total RNA was isolated from A549, H1299, H1975, and Beas-2B using TRIzol (Invitrogen, Carlsbad, CA, United States ). Then, RNA was reverse-transcribed onto cDNA using the Evo M-MLV RT Kit (Accurate Biotechnology (Hunan) Co., Ltd.) and subsequently subjected to qRT-PCR using the SYBR Green Premix Pro Taq HS qPCR Kit (Accurate Biology) according to the manufacturer’s instructions. We used GAPDH as the reference gene for 17 prognostic m⁶A-related lncRNAs. The primer sequences were synthesized by BioSune Co., Ltd. (Shanghai, China) as detailed in Supplementary Table S1. The value of the relative expression was calculated by the 2^−ΔΔCT method.

Statistical Analysis

Each experiment was repeated at least thrice. All data were analyzed using GraphPad Prism 7.0 and R software. Significant differences between A549, H1299, H1975, and Beas-2B cell were evaluated by one-way ANOVA. P-value less than 0.05 was considered as statistically significant.

Results

Identification of m⁶A-Associated lncRNAs in LUAD Patients

A flow diagram of the whole study is shown in Figure 1. We downloaded RNA-seq transcription expression data for TCGA-LUAD patients. Then, we extracted the expression of 19 m⁶A-related regulators (METTL3, VIRMA, METTL14, ZC3H13, RBM15, RBM15B, WTAP, HNRNPC, YTHDC1, HNRNPA2B1, YTHDC2, IGF2BP1, YTHDF1, IGF2BP2, YTHDF2, IGF2BP3, YTHDF3, ALKBH5, and FTO) from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). Next, according to the standard filters of p < 0.001 and correlation coefficient >0.4, 2,710 correlated pairs of m⁶A-related regulators and m⁶A-related lncRNA correlations were identified (Supplementary Table S2) (Figure 2A).

FIGURE 1

FIGURE 1. Flow chart of the whole study.

FIGURE 2

FIGURE 2. m6A regulator–lncRNA coexpression network construction, differential analysis, and hierarchical cluster analyses in LUAD. (A) Coexpression network between m⁶A regulators and m⁶A-related lncRNAs. Red circles indicate m6A regulators; blue diamonds indicate m⁶A-related lncRNAs. (B) Heatmap showing the top 20 upregulated and top 20 downregulated differentially expressed m⁶A-related lncRNAs between the GS-like group and GU-like group. (C) Hierarchical cluster analyses based on the expression of differentially expressed m⁶A-related lncRNAs.

Identification of Genome Instability–Derived m⁶A-Related lncRNAs in LUAD Patients

To identify genome instability–derived m⁶A-related lncRNAs, we first downloaded TCGA somatic mutation data from the GDC Data Portal. Then, we calculated the cumulative number of somatic mutations per LUAD patient and sorted patients in descending order based on this number; the lowest 25% (n = 134) and the top 25% (n = 139) patients were defined as the GS-like group and GU-like group, respectively (Figure 2B). Next, 42 genome instability–derived m⁶A-related lncRNAs were differentially expressed between the GS and GU groups based on |logFC|≥1 and false discovery rate adjusted p-value < 0.05. In total, 12 genome instability–derived m⁶A-related lncRNAs were upregulated and 30 were downregulated in the GU group (Table 1). Unsupervised hierarchical clustering analysis was performed for 535 LUAD patients from TCGA dataset using the set of 42 differentially expressed genome instability–derived m⁶A-related lncRNAs (Figure 2C). All 535 LUAD patients were clustered into two cohorts according to the expression levels of differentially expressed genome instability–derived m⁶A-related lncRNAs. The somatic mutation count was significantly different between the two cohorts. We defined the cohort with the lower cumulative somatic mutation count as the GS group and the one with the higher cumulative somatic mutation count as the GU group. Then, we analyzed the differential expression of 19 m⁶A-related regulators between the GU and GS groups. As shown in Figure 3, the expression levels of FTO, HNRNPC, HNRNPA2B1, IGF2BP3, IGF2BP1, IGF2BP2, YTHDC2, YTHDF1, YTHDF3, METTL14, RBM15, VIRMA, RBM15B, and WTAP differed significantly between GU and GS patients (p < 0.05). To explore the role of immune checkpoint molecules, we performed differential expression analysis for CTLA4 (T-lymphocyte–associated protein 4), HAVCR2 (hepatitis A virus cellular receptor 2), PDCD1, and 4-1 BB (TNFRSF9) between the two groups. CTLA4 showed higher expression in the GU group than in the GS group (p = 0.015), whereas the reverse was true for HAVCR2 (p = 8.2e-05). However, PDCD1 and 4-1 BB (TNFRSF9) had no significant difference between the two groups (p = 0.08, p = 0.058). Furthermore, the somatic mutation counts of LUAD patients showed significant differences between the two groups (p < 2.22e-16).

TABLE 1

TABLE 1. Differentially expressed genomic instability–derived m⁶A-related lncRNAs in the GS-like cohort and GU-like cohort.

FIGURE 3

FIGURE 3. Correlation analysis of LUAD clusters, expression of m⁶A regulators, CTLA4, HAVCR2, PDCD1, 4-1BB (TNFRSF9), and somatic mutation counts of LUAD. (A–N) Correlation analysis of LUAD clusters and expression of m⁶A regulators. (O–S) Correlation analysis of LUAD clusters and expression of CTLA4, HAVCR2, PDCD1, and 4-1BB (TNFRSF9). (Q) Correlation analysis of LUAD clusters and somatic mutation counts.

CeRNA Network Construction

To further identify the mechanism of genome instability–associated m⁶A-related lncRNAs in LUAD patients, we established the genomic instability of the lncRNA-mediated ceRNA network. First, we searched for differentially expressed genome instability–derived m⁶A-related lncRNAs in the miRcode database and presented 306 pairs of interacting miRNAs and lncRNAs using Perl language. Then, we identified 17 of 298 identified differentially expressed miRNAs (DEmiRNAs) interacting with four differentially expressed genome instability–derived m⁶A-related lncRNAs. The interactions between differentially expressed miRNAs and lncRNAs are shown in Table 2. Next, we identified 17 DEmiRNAs targeting 436mRNAs using TargetScan ((http://www.targetscan.org/), miRTarbase ((http://mirtarbase.mbc.nctu.edu.tw/), and miRDB (http://www.mirdb.org/). The correlations between mRNAs and miRNAs are shown in Table 3. Finally, 75 mRNAs were used to construct the network of genome instability–derived m⁶A-related lncRNAs. Furthermore, we constructed a network of genome instability–derived m⁶A-related lncRNA-mediated ceRNA network incorporating four differentially expressed lncRNAs (DElncRNAs), 17DEmiRNAs, and 75 differentially expressed mRNAs (Figure 4).

TABLE 2

TABLE 2. Interaction between differentially expressed m⁶A-related lncRNAs and miRNAs associated with genomic instability.

TABLE 3

TABLE 3. Interaction between differentially expressed miRNAs and mRNAs associated with genomic instability.

FIGURE 4

FIGURE 4. Construction of the ceRNA network. Red and blue circles indicate upregulated and downregulated mRNAs, respectively. Red and blue rectangles indicate upregulated and downregulated miRNAs, respectively. Red and blue diamonds indicate upregulated and downregulated genomic instability–derived differentially expressed m⁶A-related lncRNAs, respectively.

Functional and ConsensusPathDB Analysis of Genome Instability–Derived m⁶A-Related lncRNAs

To further explore the potential functions and pathways involving the genome instability–derived m⁶A-related lncRNAs, we conducted functional enrichment and ConsensusPathDB pathway analyses. According to the genome instability–derived m⁶A-related lncRNA-mediated ceRNA network, we obtained 75 lncRNA-correlated mRNAs. As shown in Figure 5, the 75 mRNAs were most enriched in eight GO terms: the cellular response to fatty acids, lipid metabolic process, protein kinase binding, protein binding, signal transduction, cellular response to amino acid stimulus, negative regulation of gene expression, and receptor signaling protein serine/threonine kinase activity. The ConsensusPathDB pathway analyses demonstrated that 75 mRNAs were most enriched in 14 pathways: sudden infant death syndrome (SIDS) susceptibility pathways, regulation of the microtubule cytoskeleton, ncRNAs involved in STAT3 signaling in hepatocellular carcinoma, integrins in angiogenesis, SMAD2/3 phosphorylation motif mutants in cancer, SMAD2/3 MH2 domain mutants in cancer, loss of function of SMAD2/3 in cancer, signaling by TGF-beta receptor complex in cancer, epithelial–mesenchymal transition in colorectal cancer, syndecan-1–mediated signaling events, AGE-RAGE signaling pathway in diabetic complications – Homo sapiens (human), collagen chain trimerization, beta3 integrin cell surface interactions, and HIF-1 signaling pathway – Homo sapiens (human) (Figure 6) (Table 4).

FIGURE 5

FIGURE 5. Functional enrichment analysis. (A–C) GO analysis indicating top eight GO terms for which differentially expressed mRNAs were most enriched.

FIGURE 6

FIGURE 6. ConsensusPathDB analysis indicating that differentially expressed mRNAs were most enriched in 14 pathways.

TABLE 4

TABLE 4. ConsensusPathDB pathway analysis.

Construction of the Cox Prediction Model of m⁶A-Associated lncRNAs in LUAD

First, we tried constructing a Cox prediction model of genomic instability of m⁶A-related lncRNAs in LUAD patients. However, the results showed that there was only one group in LUAD patients. Then, we further explored whether m⁶A-associated lncRNAs could serve as prognostic biomarkers for predicting the clinical outcomes of LUAD patients and established a Cox prediction model. First, we combined the expression level of m⁶A-related lncRNAs and the complete survival status and survival time data of 535 LUAD patients, Then, the patients were divided into two cohorts in the ratio of 7:3: training cohort (n = 356) and testing cohort (n = 148). Next, univariate Cox regression analysis of the training TCGA cohort showed that 39 m⁶A-related lncRNAs were significantly related to overall survival (p < 0.05) (Figure 7A) (Table 5). We selected these 39 prognosis-related lncRNAs for multivariate Cox regression analysis. This yielded 17 m⁶A-related lncRNAs which were used to construct the Cox prediction model (Figure 7B) (Table 6). The weighted relative coefficients in the Cox prediction model were as follows: riskscore = (-0.1938)×AL122010.1 expression + (-0.2721) × STIM2-AS1 expression + (-0.3821) × LINC00654 expression + (-1.6785) × AL133445.2 expression +0.300 × LINC01137 expression + (-0.1602) × GAS6-AS1 expression + (-0.2401) × AC090617.5 expression +0.8548 × AC093495.1 expression + (-0.4769) × AC026202.2 expression +0.0733 × AL590666.2 expression + (-0.6618) × AC123595.1 expression + (-0.5149) × AL590226.1 expression +0.0734 × AC245041.1 expression +0.0648 × LINC02555 expression + (-0.0900) × AL049555.1 expression + (-0.2295) × AC024075.1 expression +0.2374 × AC079949.2 expression). According to the median risk scores (1.62, 1.69, and 1.41), we divided our entire TCGA set, training TCGA set, and testing TCGA set into high-risk and low-risk groups. Survival analysis on these three sets, (entire, training, and testing) showed that patients in the high-risk groups had worse prognosis than those in the low-risk groups (p < 0.05) (Figures 8A,E,I). Time-dependent ROC curves for the entire, training, and testing sets at 1 year, 2 years, and 3 years showed that the Cox prediction model could moderately well-predict the overall survival of LUAD patients (Figures 8B,F,J). Independent prognosis analyses showed that pathological N stage and riskscore were independent prognostic factors for overall survival in the entire TCGA-LUAD cohort (p < 0.05) (Figure 8D). Riskscore was an independent prognostic factor for survival in TCGA training LUAD cohort (p < 0.05) (Figure 8H), and pathological T and N stages were independent prognostic factors for clinical outcomes in the testing TCGA-LUAD cohort (p < 0.05) (Figure 8L).

FIGURE 7

FIGURE 7. Univariate and multivariate Cox regression analyses. (A) Forest map of univariate Cox regression analysis. (B) Forest map of multivariate Cox regression analysis. Red and green circles indicate high-risk and low-risk factors, respectively.

TABLE 5

TABLE 5. Univariate Cox regression analysis of training TCGA sets of prognostic m⁶A-related lncRNAs in LUAD.

TABLE 6

TABLE 6. Multivariate Cox regression analysis of m⁶A-related lncRNAs in LUAD.

FIGURE 8

FIGURE 8. Prognostic value of m⁶A-related lncRNAs of entire, training, and testing cohorts of LUAD patients. (A,E,I) Survival analysis of the prognostic model of the entire, training, and testing cohorts of LUAD patients. (B,F,J) ROC AUC at 1, 2, and 3 years in the entire, training, and testing cohorts. (C,G,K) Forest map of univariate Cox independent regression analyses of the prognostic model for the entire, training, and testing cohorts. (D,H,L) Forest map of multivariate Cox independent regression analyses of the prognostic model for the entire, training, and testing cohorts. (M) ROC curve of AUC analysis for our prediction model and other published prognostic models.

Performance Comparison of the LncSig With the Signature of 17 Prognostic m⁶A-Related lncRNAs

We compared the performance of our m⁶A-related lncRNA signature with the previously published lncRNA signatures: a five-lncRNA signature from Sun’s study (SunlncSig) (Sun et al., 2020) and a four-lncRNA signature from Shukla’s study (ShuklalncSig) (Shukla et al., 2017). We extracted the published lncRNA expression and combined it with our risk profile and complete survival information for the performance comparison analysis. As shown in Figure 8M, the AUC for overall survival for the signature of 17 m⁶A-related lncRNAs was 0.746, significantly higher than that for SunlncSig (AUC = 0.577) and ShuklalncSig (AUC = 0.681). Performance comparison analysis showed that our signature could better predict the clinical outcomes of LUAD patients than the two recently published lncRNA signatures.

Correlation Analysis of the Expression of m⁶A Regulators and the Risk Model

To investigate the relationships between m⁶A regulators and the risk model, we performed correlation analysis. As shown in Figure 8, eight m⁶A regulators (HNRNPC, YTHDC1, METTL3, IGF2BP1, METTL14, YTHDC2, IGF2BP2, and IGF2BP3) were significantly associated with the risk model (p < 0.05). In the entire LUAD patient cohort, HNRNPC, IGF2BP1, IGF2BP2, and IGF2BP3 had higher expression in the high-risk group than the low-risk group (p < 0.05; Figures 9A–D), whereas YTHDC1, YTHDC2, METTL3, and METTL14 had higher expression in the low-risk group (p < 0.05; Figures 9E–H). In the training cohort (n = 356), the expression of HNRNPC, IGF2BP1, IGF2BP2, and IGF2BP3 was significantly higher in the high-risk group than in the low-risk group (p < 0.05; Figure 8I-L), whereas the expression of YTHDC1, YTHDC2, METTL3, and METTL14 was significantly higher in the low-risk group (p < 0.05; Figures 9M–P). In the testing cohort (n = 148), the expression levels of HNRNPC, IGF2BP2, and IGF2BP3 were higher in the high-risk group (p < 0.05; Figure 9Q, S, T), whereas IGF2BP1 and METTL14 showed no significant difference in expression between the two risk groups (p > 0.05; Figures 9R, X), and YTHDC1, YTHDC2, and METTL3 expression was higher in the low-risk group (p < 0.05; Figures 9U–W).

FIGURE 9

FIGURE 9. Correlation analysis of m⁶A regulators and the prognostic risk model. (A–X) Differential analysis of m⁶A regulators in the prognostic model for the high- and low-risk cohorts.

Correlation Analysis Between Somatic Mutation Count and the Risk Model in the LUAD Cohort

We performed correlation analysis between somatic mutation count and the risk model. As shown in Figures 10A–C, the somatic mutation counts differed significantly between the high- and low-risk groups in the entire training and testing LUAD cohorts (p < 0.05). As shown in Figure 10, the somatic mutation counts were significantly higher in the high-risk groups than in the low-risk groups in all three cohorts (p < 0.05).

FIGURE 10

FIGURE 10. Correlation analysis of the prognostic risk model and somatic mutation counts of LUAD. (A–C) Differential analysis of somatic mutation counts in the prognostic model for the high- and low-risk cohorts.

Riskplot Analyses

To reveal the relationships between the riskplot and expression and genomic instability of HNRNPC, IGF2BP1, METTL14, IGF2BP2, YTHDC1, IGF2BP2, YTHDC2, and METTL3, we performed riskplot analyses. In the training group, LUAD patient risk increased and the expression levels of LINC01137, AC245041.1, AL049555.1, AL590666.2, and AC079949.2 increased, whereas those of AL590226.1, LINC02555, STIM2-AS1, AL133445.2, AC123595.1, AC090617.5, AC026202.2, LINC00654, GAS6-AS1, AL122010.1, AC093495.1, and AC024075.1 decreased (Figure 11A). In the testing group of prognostic m⁶A-related lncRNA signature the heatmap was consistent with that in the training group; LINC01137, AC245041.1, AL049555.1, AL590666.2, AC079949.2 were high-risk factors, whereas AL590226.1, LINC02555, STIM2-AS1, AL133445.2, AC123595.1, AC090617.5, AC026202.2, LINC00654, GAS6-AS1, AL122010.1, AC093495.1, and AC024075.1 were low-risk factors (Figure 11B). In the training andtesting groups, LUAD patient risk increased and somatic mutation counts decreased (Figures 11C,D). There were differences in the expression of genomic instability of HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, METTL3, METTL14, YTHDC1, and YTHDC2 between high- and low-risk patients in the training set (Figures 11E,G,I,K, M, O, Q, S). There were also differences in the expression of genomic instability of HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, METTL3, METTL14, YTHDC1, and YTHDC2 between high- and low-risk patients in the testing set (Figures 11F,H,J, L, N, P, R, T).

FIGURE 11

FIGURE 11. Riskplot analysis of m⁶A-related regulators (A,C) Riskplot heatmap and somatic mutation counts for the prognostic model in the training cohort of LUAD patients. (B,D) Riskplot heatmap and somatic mutation counts for the prognostic model in the testing cohort of LUAD patients. (E,G,I,K,M,O,Q,S) Differences in expression of HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, METTL3, METTL14, YTHDC1, and YTHDC2 in the high- and low-risk groups of the training cohort. (F,H,J,L,N,P,R,T) Differences in the expression of HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, METTL3, METTL14, YTHDC1, and YTHDC2 in the high- and low-risk groups of the testing cohort.

Prognostic Model Validation for Clinical Grouping

We combined the clinical data of 522 LUAD patients from TCGA with complete survival information and the risk file of the prognosis model for survival verification of clinical samples. Patients with LUAD in the high-risk cohort had worse prognosis than those in the low-risk cohort, regardless of whether they were older than 65 years or younger than 65 years (p < 0.001) (Figures 12A,B), or whether they were male or female (p < 0.001) (Figures 12C,D). Patients with pathological M0 stage in the high-risk cohort had poorer prognosis than those in the low-risk cohort (p < 0.001) (Figure 12E). However, there was no such significant difference with respect to pathological M1 stage between the high- and low-risk cohorts (p = 0.169) (Figure 12F). Among LUAD patients with pathological N0/N1-3 stage, the low-risk cohort had better prognosis than the high-risk cohort (p Figures 12G,H< 0.001) (Figures 12G,H). LUAD patients with pathological stages I–II/III–IV in the high-risk cohort had poorer prognosis than those in the low-risk cohort (p < 0.001) (Figures 12I,J). LUAD patients with pathological T1–2/3–4 stages in the high-risk cohort had poorer survival than those in the low-risk cohort (p < 0.05) (Figures 12K,L)

FIGURE 12

FIGURE 12. Prognostic model validation of clinical grouping. (A–L) Survival analysis of LUAD patients by age, gender, pathological M0/M1 stage, pathological T1–2/3–4 stage, pathological N0/N1–3 stage, and pathological stage I–II/III–IV in the high- and low-risk cohorts.

Prognosis Validation of m⁶A-Related lncRNAs

We performed survival analysis of the 17 prognosis m⁶A-related lncRNAs (AL122010.1, STIM2-AS1, LINC00654, AL133445.2, GAS6-AS1, AC090617.5, AC093495.1, AC026202.2, AL590666.2, AC123585, AL590226.1, AC245041.1, LINC02555, AL049555.1, AC024075.1, AC079949.2, and LINC01137) from TCGA (https://portal.gdc.cancer.gov/) LUAD (lung adenocarcinoma) project in level 3 HTSeq-FPKM (fragments per kilobase per million) format RNA-seq data and survival information of 594 LUAD patients from TCGA database to verify the prognosis of the 17 lncRNAs. As shown in Supplementary Figure S1, patients with high expression levels of STIM2-AS1, LINC00654, AL133445.2, GAS6-AS1, AC090617.5, AC123595.1, and AC024075.1 had better prognosis than those with low expression levels, showing that these m⁶A-related lncRNAs are protective factors regarding the prognosis of LUAD patients, while the patients with high expression of AL590666.2, AL049555.1, and LINC01137 had poor prognosis than those with low expression, indicating that AL590666.2, AL049555.1, and LINC01137 are risk factors for predicting the survival of LUAD.

Expression, Clinical Correlation, Immune Microenvironment, and Anticancer Drug Sensitivity Analyses of 17 Prognostic m⁶A-Related lncRNAs in LUAD

To investigate the role of the 17 m⁶A-associated lncRNAs in LUAD, we downloaded the LUAD data from UCSC-Xena TCGA GDC to conduct differential expression, LUAD immune subtype, LUAD immune microenvironment, LUAD stemness score (DNAss and RNAss), and LUAD drug sensitivity analyses. As shown in Figure 13A, LINC01137 had higher expression levels than the other prognostic lncRNAs in LUAD. The expression levels of AC093495.1 and AC024075.1, and those of AL133445.2 and AC026202.2, were positively correlated with LUAD (Figure 13B). Expression levels of AL122010.1, STIM2-AS1, AL133445.2, GAS6-AS1, AC093495.1, AC026202.2, AC123595.1, AL590226.1, AC245041.1, LINC02555, AL049555.1, AC024075.1, and AC079949.2 were significantly different in immune subtypes C1 (wound healing), C4 (lymphocyte depleted), C3 (inflammatory), C2 (IFN-gamma dominant), and C6 (TGF-beta dominant) (Figure 13C). There was no statistically significant difference in the expression of the 17 prognostic m⁶A-related lncRNAs between pathological types T1/T2/T3/T4 (Figure 13D). However, the expression of AL590666.2 was significantly different in pathological stage M0/M1 of LUAD (p < 0.01) (Figure 13E). The expression levels of m⁶A-related lncRNAs STIM2-AS1, LINC00654, LINC01137, GAS6-AS1, AC090617.5, AC093495.1, AL590666.2, AC123595.1, AC024075.1, and AC079949.2 were significantly different in pathological stages N0/N1/N2/N3 (p < 0.05) (Figure 13F). As shown in Figure 13G, most of the m⁶A-related lncRNAs had negative correlations of their expression with RNAss and DNAss in LUAD, that is, higher expression of m6A-related lncRNAs, lower RNAss and DNAss stemness scores, weaker activity of LUAD stem cells, and greater degree of LUAD differentiation. The expression of LINC01137 was negatively correlated with RNAss in LUAD. The expression levels of most m⁶A-related lncRNAs were positively correlated with immune score, estimate score, and stromal score, which indicated that the content of immune and stromal cells was higher in LUAD, whereas tumor cell contents were lower in LUAD, indicating that most of the m⁶A-related lncRNAs were protective factors with respect to the survival of LUAD patients. As shown in Figure 14, LINC00654 expression was positively correlated with sensitivity to anticancer drugs procarbazine (p < 0.001), simvastatin (p < 0.001), testolactone (p = 0.003), and calusterone (p < 0.001) but negatively associated with sensitivity to anticancer drug byproduct of CUDC-305 (p < 0.001). The expression of GAS6-AS1 had a positive relationship with sensitivity to anticancer drugs chelerythrine (p < 0.001), hydroxyurea (p < 0.001), melphalan (p < 0.001), nelarabine (p < 0.001), cyclophosphamide (p = 0.001), pipobroman (p = 0.001), idarubicin (p = 0.001), cytarabine (p = 0.001), BML-277 (p = 0.002), imexon (p = 0.002), ABT-199 (p = 0.002), cladribine (p = 0.002), and uracil mustard (p = 0.003) but was negatively correlated with kahalidef sensitivity (p < 0.001). The expression of LINC01137 was negatively correlated with sensitivity to anticancer 3-bromopyruvate (acid) in LUAD (p = 0.002).

FIGURE 13

FIGURE 13. Expression, clinical correlation, and immune microenvironment analyses of prognostic m⁶A-related lncRNAs in LUAD. (A) Expression of 17 prognostic m⁶A-related lncRNAs in LUAD. (B) Correlation analysis of 17 prognostic m⁶A-related lncRNAs. (C) Boxplot of 17 prognostic m⁶A-related lncRNAs in six LUAD immune subtypes. (D–F) Boxplots of 17 prognostic m⁶A-related lncRNAs and clinical features in LUAD. (G) Correlations between the expression of 17 prognostic m⁶A-related lncRNAs and tumor stemness scores (based on RNA methylation and DNA methylation), immune score, estimate score, and stromal score in LUAD.

FIGURE 14

FIGURE 14. Correlation analysis between m⁶A-related lncRNAs and sensitivity to anticancer drugs in CellMiner. Top 20 statistically significant correlations between prognostic m⁶A-related lncRNAs and sensitivity to anticancer drugs.

Verification of Expression of 17 Prognostic m⁶A-Associated lncRNAs in LUAD

The expression levels of LINC01137, AL590666.2, AL049555.1, AC245041.1, LINC02555, AL590226.1, AL122010.1, AC093495.1, and AC090675.1 were significantly different between 57 paired LUAD and adjacent non-LUAD tissues (p < 0.001). The expressions of AC123595.1 and AC026202.2 in paired LUAD tissues were higher than those in paired non-LUAD tissues (p < 0.01). The expressions of AL133445.2 and AL123595.1 in paired LUAD tissues were lower than those in paired non-LUAD tissues (p < 0.01). The expressions of STIM2-AS1, LINC01137, AL590666.2, AL049555.1, AC079949.2, and AC026202.2 in unpaired LUAD tissues were higher than those in unpaired non-LUAD tissues (p < 0.001). The expressions of LINC02555, AL590226.1, AL133445.2, AL122010.1, AC245041.1, AC123595.1, AC093495.1, AC090617.5, and AC024075.1 in unpaired LUAD tissues were lower than those in unpaired non-LUAD tissues (p < 0.001). The expressions of 17 prognostic m⁶A-associated lncRNAs in paired LUAD, adjacent non-LUAD tissues, and unpaired LUAD and non-LUAD tissues are shown in Supplementary Figures S2, 3. The differential expression of 17 prognostic m⁶A-associated lncRNAs in LUAD and adjacent non-LUAD tissues is shown in Supplementary Table S3. Supplementary Table 3 showed that the log(fold change (FC)) of AL590226.1, AC245041.1, LINC02555, AC024075.1, AC090617.5, AC123595.1, and AC093495.1 was less than 0, which indicated that the expression of AL590226.1, AC245041.1, LINC02555, AC024075.1, AC090617.5, AC123595.1, and AC093495.1 was highly expressed in non-LUAD tissues compared with LUAD tissues (p < 0.05). The log FC of AL133445.2, AC079949.2, LINC00654, STIM2-AS1, GAS6-AS1, AC026202.2, LINC01137, AL049555.1, and AL590666.2 was greater than 0, which demonstrated that the expression of AL133445.2, AC079949.2, LINC00654, STIM2-AS1, AC026202.2, LINC01137, AL049555.1, and AL590666.2 was highly expressed in LUAD tissues than non-LUAD tissues (p < 0.05).

Validation of 17 Prognostic m⁶A-Related lncRNAs in Human LUAD Cells and Normal Bronchial Epithelial Cells

To further evaluate whether the expression of 17 lncRNAs is significantly different between normal bronchial epithelial cells (Beas-2B) and LUAD cell lines (A549, H1299, and H1975), we performed the RT-qPCR. The results showed that the expression of LINC00654, AC090617.5, AC093495.1, AL590226.1, AC245041.1, AC123595.1, AL049555.1, AC079949.2, and LINC01137 was significantly higher in A549 than in Beas-2B (p < 0.05). The expression of AC079949.2, LINC01137, and AC024075.1 was significantly higher in H1299 than in Beas-2B (p < 0.05). The expression of AL122010.1, LINC00654, AC090617.5, GAS6-AS1, AC093495.1, AC026202.2, AL590226.1, AC245041.1, AL049555.1, AC079949.2, AL590666.2, LINC01137, and AC024075.1 was significantly different between Beas-2B and H1975 (p < 0.05). However, the expression of STIM2-AS1, AL133445.2, and LINC02555 had no significant difference between Beas-2B and A549, H1299, and H1975 (p > 0.05). The validation results are shown in Supplementary Figure S4.

Validation of 22 Immune Cell Subtypes in High and Low Expression of 17 Prognostic m⁶A-Associated lncRNAs

To further validate the relationship between 17 prognostic m⁶A-associated lncRNAs and the immune microenvironment, we performed the CIBERSORT algorithm. The results showed that the fraction of naïve B cells was significantly different in the high- and low-expression groups of AL122010.1 and LINC00654 (p = 0.008, p = 0.031); the fraction of memory B cells was significantly different in the low- and high-expression groups of GAS6-AS1 (p = 0.017); the fraction of plasma cells was significantly different in the high- and low-expression groups of AC026202.2, AC090617.5, AC245041.1, AL049555.1, AL122010.1, LINC02555, and STIM2-AS1; the fraction of CD8 T cells was significantly different in the low- and high-expression groups of AC245041.1, AL049555.1, and LINC01137; the fraction of resting memory CD4 T cells was significantly different in the low- and high-expression groups of AC024075.1,AC093495.1, AL590226.1, GAS6-AS1, LINC00654, and LINC02555; the fraction of memory activated T cells was significantly different in the high- and low-expression groups of AC093495.1, AL049555.1, AL122010.1, AL590226.1, GAS6-AS1, and LINC01137; the fraction of follicular helper T cells was significantly different in the high- and low-expression groups of LINC01137 and AC024075.1; the fraction of regulatory T cells (Tregs) was significantly different in the high- and low-expression groups of AC090617.5, AC093495.1, AC245041.1, AL049555.1, GAS6-AS1, LINC01137, and LINC02555; the fraction of gamma delta T cells was significantly different in the high- and low-expression groups of AC026202.2, AL122010.1, and AL590666.2; the fraction of resting NK cells was significantly different in the high- and low-expression groups of LINC01137 and AL590226.1; the fraction of activated NK cells was significantly different in the high- and low-expression groups of AC123595.1, AL590666.2, and LINC01137; the fraction of monocytes was significantly different in the high- and low-expression groups of AC026202.2, GAS6-AS1, and LINC02555; the fraction of M0 macrophages was significantly different in the high- and low-expression groups of AC026202.2, AL590226.1, LINC01137, and LINC02555; the fraction of M1 macrophages was significantly different in the high- and low-expression groups of AC123595.1, AC245041.1, AL049555.1, LINC01137, and STIM2-AS1; the fraction of macrophages was significantly different in the high- and low-expression groups of AC026202.2, AC245041.1, AL122010.1, AL590226.1, LINC02555, and STIM2-AS1; the fraction of resting dendritic cells was significantly different in the low- and high-expression groups of AC093495.1, AL133445.2, AL590226.1, GAS6-AS1, LINC01137, and LINC02555; the fraction of activated dendritic cells was significantly different in the high- and low-expression groups of AC024075.1, AC093495.1, AL590666.2, LINC01137, and LINC02555; the fraction of resting mast cells was significantly different in the low- and high-expression groups of AC024075.1, AC093495.1, AL590226.1, GAS6-AS1, LINC01137, and LINC02555; the fraction of activated mast cells was significantly different in the high- and low-expression groups of AC093495.1, AL590226.1, and LINC02555; the fraction of eosinophils was significantly different in the low- and high-expression groups of AC245041.1; and the fraction of neutrophils was significantly different in the high- and low-expression groups of AC026202.2, AC079949.2, AC245041.1, AL590666.2, LINC01137, and LINC02555. The CIBERSORT immune microenvironment analysis result is shown in Supplementary Figure S5.

Discussion

In this study, we first identified m⁶A-related lncRNAs by correlation analysis. Then, we identified differentially expressed genomic instability–associated m⁶A-related lncRNAs and constructed a ceRNA network and combined m⁶A-related lncRNA expression with somatic mutation profiles based on the tumor genome to explore the ability of the prognostic m⁶A-associated lncRNA signature to predict outcomes in LUAD. We also identified differences in the expression of immune checkpoint blockade CTLA4, HAVCR2, PDCD1, and 4-1 BB (TNFRSF9) in the GS-like and GU-like groups, which provided novel insight into the correlation relationship between immune checkpoint blockade and genomic instability of LUAD. Furthermore, we constructed an m⁶A-related lncRNA risk model and compared it with other reported lncRNA models, which revealed novel insights for predicting the outcome of LUAD. Eventually, pan-cancer analysis of prognostic m⁶A-related lncRNAs provided new prospects for identifying pan-cancer therapeutic targets. Previous studies have demonstrated the pivotal role of m⁶A modification of lncRNAs in cancers (Fazi and Fatica, 2019). Recently, a study demonstrated the role of oncogenic lncRNA THOR in m⁶A modification (Liu et al., 2020a). Additionally, the m⁶A-induced lncRNA RP11 was identified as a new predictive biomarker that triggers the dissemination of colorectal cancer metastasis through upregulation of ZEB1 (Wu et al., 2019). Some studies showed that the m⁶A RNA methyltransferase METTL3/14 could enhance response to anti–PD-1 treatment in pMMR-MSI-L melanoma and CRC (Wang et al., 2020a). However, the pivotal role of genomic instability of m⁶A-related lncRNAs in LUAD remained unclear.

In this study, we first identified genomic instability of the m⁶A-related lncRNA-associated ceRNA network, revealing novel insights related to new therapeutic targets involving RNA epigenetic mechanisms. We constructed a genomic instability–associated DElncRNA-mediated ceRNA network and found that four DElncRNAs (SRGAP3-AS2, AC110619.1, ATP13A4-AS1, and PSORS1C3) were significantly different between the GS-like group and GU-like groups. Recently, studies have revealed that SRGAP3-AS2 could play an important role in predicting the prognosis of LUAD (Yang et al., 2018; Jin et al., 2020). However, whether the expression of AC110619.1, ATP13A4-AS1, PSORS1C3, and SRGAP3-AS2 was associated with genomic instability was still unclear. We first constructed a ceRNA network of genomic instability–associated m⁶A-related lncRNAs, revealing novel insights for exploring new RNA epigenetic regulatory mechanisms in LUAD. Furthermore, we combined the expression of m⁶A regulators and immune checkpoint blockades CTLA4 and HAVCR2 with somatic mutations of LUAD, providing a new approach to exploring novel immune therapeutic targets in LUAD. Recently, a study showed that CTLA4 can serve as a prognostic biomarker for predicting the survival of LUAD (Wang et al., 2020b). Furthermore, an immunogenic gene signature associated with immune checkpoint blockade provided novel therapeutic targets for immunology in LUAD (Ahluwalia et al., 2021).

In our study, we identified a signature of 17 prognostic m⁶A-related lncRNAs based on all LUAD patients, a training group, and a testing group for predicting the outcomes of LUAD. Furthermore, we conducted a performance comparison analysis of our prognostic model and two recently published lncRNA signatures; the results indicated that our prognostic model of m⁶A-associated lncRNAs performed better in predicting the clinical outcomes of LUAD patients than the other two lncRNA signatures. In addition, validation using the GEO dataset GSE102287 confirmed that high expression of LINC01137 was associated with better prognosis of LUAD patients than low expression of LINC01137. Some studies have identified AL122010.1 as a predictor of survival of breast cancer patients (Li et al., 2020; Ma et al., 2020). The short-lived lncRNA LINC01137 serves as a useful indicator of chemical stress responses (Tani et al., 2019). Low expression of lncRNA GAS6-AS1 is a biomarker for predicting survival in NSCLC (Han et al., 2013). The expression of AC079949.2 is associated with clinical outcomes of patients with esophageal cancer (Liu et al., 2020b).

In our study, we identified the correlations between the expression of 17 prognostic m⁶A-related lncRNAs and the LUAD immune microenvironment, microenvironment, stemness scores, and immune subtype using the UCSC-Xena website. The immune microenvironment analysis showed that the 17 prognostic m⁶A-related lncRNAs were protective factors regarding clinical outcome on a pan-cancer basis. Recently, the m⁶A-related lncRNA LINC00958 was identified as a therapeutic target in hepatocellular carcinoma (Zuo et al., 2020). Other studies showed that ALKBH5 could promote the metastasis and invasion of gastric cancer by demethylating the expression of lncRNA NEAT1 (Zhang et al., 2019b) and that METTL14 acts as a prognostic biomarker in colorectal cancer by downregulating oncogenic lncRNA XIST (Yang et al., 2020). We constructed an m⁶A-related lncRNA prognostic risk model based on all LUAD patients, a training group, and a testing group from UCSC-Xena GDC TCGA-LUAD. Performance comparison analysis confirmed that the m⁶A-related lncRNA risk model could better predict clinical outcomes of LUAD patients than the previously published lncRNA signatures. Our study first verified the relationship between prognostic m⁶A-associated lncRNAs and the immune microenvironment based on the CIBERSORT algorithm, which provided a novel insight for revealing potential immune therapeutic targets in RNA epigenetics.

We first identified prognostic genomic instability of m⁶A-related lncRNAs in LUAD and constructed and validated a prognostic model based on the internal datasets (all TCGA-LUAD, training TCGA-LUAD, and testing TCGA-LUAD); furthermore, we validated the prognostic value of m⁶A-related lncRNA LINC01137 in the external dataset GSE102287 and analyzed the expression, clinical correlation, immune microenvironment, and anticancer drug sensitivity of the prognostic value of the 17 m⁶A-related lncRNAs in LUAD; the results provided novel insights for exploring new therapeutic targets in LUAD. Our study subjects were extracted from TCGA, the UCSC-Xena database, and the GEO database; it is not known whether our findings apply to other groups. We further validated the expression of prognostic m⁶A-related lncRNAs in normal bronchial epithelial (Beas-2B) and LUAD cells (A549, H1299, and H1975), and the results showed that the expression of AL122010.1, AC026202.2, AC123595.1, AL049555.1, AC079949.2, and LINC01137 was significantly higher in LUAD cells than normal bronchial epithelial cells (p < 0.05), which is consistent with TCGA database. However, the reason why the expression levels of the remaining lncRNAs are contrary to the results in the database may be related to ethnicity. In future, many molecular biology experiments will be performed, and many clinical samples will be collected to verify the 17 prognostic m⁶A-related lncRNAs.

In this study, we further explored the relationships between the expression levels of 17 prognostic m⁶A-associated lncRNAs and pathological N1-3/N0 and pathological M0/M1 stages. AL590666.2 was significantly different in pathological stages M0/M1, and 10 prognostic m⁶A-associated lncRNAs were differentially expressed in pathological N0/N1-3 stages. Our study first explored the association of the expression of m⁶A-associated lncRNAs with anticancer drug sensitivity using the CellMiner database, which provided novel insights regarding new therapeutic targets based on RNA epigenetics in LUAD.

Conclusion

In conclusion, we identified the genomic instability of m⁶A-related lncRNAs and constructed a ceRNA network and then constructed and validated the Cox prediction model based on m⁶A-related lncRNAs using all TCGA-LUAD patients, a training group, and a testing group. Finally, clinical, prognostic, immune microenvironment, stemness scores, and anticancer drug sensitivity analyses of m⁶A-related lncRNAs in LUAD shed new light on potential therapeutic targets based on RNA epigenetics.

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.

Author Contributions

Y-QQ and RL designed the study; J-PL and T-TL collected the data; CH, JY, and X-LJ analyzed the data; RL wrote the manuscript; Y-QQ revised the manuscript. All authors approved the submitted version.

Funding

This work was funded by grants from the Major Scientific and Technological Innovation Project of Shandong Province (Grant No. 2018CXGC1212), the CSCO-Qilu Cancer Research Fund (Grant No. Y-Q201802-014), the Medical and Health Technology Innovation Plan of Jinan City (Grant No. 201805002), and the Special Fund for Clinical Research of Jinan City (201912011).

Conflict of Interest

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

Publisher’s Note

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

Supplementary Material

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

References

Ahluwalia, P., Ahluwalia, M., Mondal, A. K., Sahajpal, N., Kota, V., Rojiani, M. V., et al. (2021). Immunogenomic Gene Signature of Cell-Death Associated Genes with Prognostic Implications in Lung Cancer. Cancers 13, 155. doi:10.3390/cancers13010155

CrossRef Full Text | Google Scholar

Andor, N., Maley, C. C., and Ji, H. P. (2017). Genomic Instability in Cancer: Teetering on the Limit of Tolerance. Cancer Res. 77, 2179–2185. doi:10.1158/0008-5472.CAN-16-1553

PubMed Abstract | CrossRef Full Text | Google Scholar

Bao, S., Zhao, H., Yuan, J., Fan, D., Zhang, Z., Su, J., et al. (2020). Computational Identification of Mutator-Derived lncRNA Signatures of Genome Instability for Improving the Clinical Outcome of Cancers: a Case Study in Breast Cancer. Brief Bioinform 21, 1742–1755. doi:10.1093/bib/bbz118

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: A Cancer J. Clinicians 68, 394–424. doi:10.3322/caac.21492

CrossRef Full Text | Google Scholar

Chu, W., Zhang, X., Qi, L., Fu, Y., Wang, P., Zhao, W., et al. (2020). The EZH2-PHACTR2-AS1-Ribosome Axis Induces Genomic Instability and Promotes Growth and Metastasis in Breast Cancer. Cancer Res. 80, 2737–2750. doi:10.1158/0008-5472.can-19-3326

PubMed Abstract | CrossRef Full Text | Google Scholar

Chua, M. L. K., Lo, W., Pintilie, M., Murgic, J., Lalonde, E., Bhandari, V., et al. (2017). A Prostate Cancer " Nimbosus ": Genomic Instability and SChLAP1 Dysregulation Underpin Aggression of Intraductal and Cribriform Subpathologies. Eur. Urol. 72, 665–674. doi:10.1016/j.eururo.2017.04.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Duijf, P. H. G., Nanayakkara, D., Nones, K., Srihari, S., Kalimutho, M., and Khanna, K. K. (2019). Mechanisms of Genomic Instability in Breast Cancer. Trends Mol. Med. 25, 595–611. doi:10.1016/j.molmed.2019.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Fazi, F., and Fatica, A. (2019). Interplay between N6-Methyladenosine (m6A) and Non-coding RNAs in Cell Development and Cancer. Front. Cel Dev. Biol. 7, 7. doi:10.3389/fcell.2019.00116

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, L., Kong, R., Yin, D.-D., Zhang, E.-B., Xu, T.-P., De, W., et al. (2013). Low Expression of Long Noncoding RNA GAS6-AS1 Predicts a Poor Prognosis in Patients with NSCLC. Med. Oncol. 30, 694. doi:10.1007/s12032-013-0694-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hypoxic Tumors Share Genomic Instability, Cancer Discov..(2019); 9: 314. doi:10.1158/2159-8290.CD-NB2019-012

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, D., Song, Y., Chen, Y., and Zhang, P. (2020). Identification of Three lncRNAs as Potential Predictive Biomarkers of Lung Adenocarcinoma. Biomed. Res. Int. 2020, 1–13. doi:10.1155/2020/7573689

PubMed Abstract | CrossRef Full Text | Google Scholar

Jordan, E. J., Kim, H. R., Arcila, M. E., Barron, D., Chakravarty, D., Gao, J., et al. (2017). Prospective Comprehensive Molecular Characterization of Lung Adenocarcinomas for Efficient Patient Matching to Approved and Emerging Therapies. Cancer Discov. 7, 596–609. doi:10.1158/2159-8290.CD-16-1337

PubMed Abstract | CrossRef Full Text | Google Scholar

Kamburov, A., Stelzl, U., Lehrach, H., and Herwig, R. (2013). The ConsensusPathDB Interaction Database: 2013 Update. Nucleic Acids Res. 41, D793–D800. doi:10.1093/nar/gks1055

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S., Kopp, F., Chang, T.-C., Sataluri, A., Chen, B., Sivakumar, S., et al. (2016). Noncoding Rna Norad Regulates Genomic Stability by Sequestering Pumilio Proteins. Cell 164, 69–80. doi:10.1016/j.cell.2015.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Li, Y., Yu, X., and Jin, F. (2020). Identification and Validation of Stemness-Related lncRNA Prognostic Signature for Breast Cancer. J. Transl Med. 18, 331. doi:10.1186/s12967-020-02497-4

CrossRef Full Text | Google Scholar

Liu, H., Xu, Y., Yao, B., Sui, T., Lai, L., and Li, Z. (2020). A Novel N6-Methyladenosine (m6A)-dependent Fate Decision for the lncRNA THOR. Cell Death Dis 11, 11. doi:10.1038/s41419-020-02833-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Zhang, Q., Lou, Q., Zhang, X., Cui, Y., Wang, P., et al. (2020). Differential Analysis of lncRNA, miRNA and mRNA Expression Profiles and the Prognostic Value of lncRNA in Esophageal Cancer. Pathol. Oncol. Res. 26, 1029–1039. doi:10.1007/s12253-019-00655-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, W., Zhao, F., Yu, X., Guan, S., Suo, H., Tao, Z., et al. (2020). Immune-related lncRNAs as Predictors of Survival in Breast Cancer: a Prognostic Signature. J. Transl Med. 18, 442. doi:10.1186/s12967-020-02522-6

CrossRef Full Text | Google Scholar

Munschauer, M., Nguyen, C. T., Sirokman, K., Hartigan, C. R., Hogstrom, L., Engreitz, J. M., et al. (2018). The NORAD lncRNA Assembles a Topoisomerase Complex Critical for Genome Stability. Nature 561, 132–136. doi:10.1038/s41586-018-0453-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Negrini, S., Gorgoulis, V. G., and Halazonetis, T. D. (2010). Genomic Instability - an Evolving Hallmark of Cancer. Nat. Rev. Mol. Cel Biol 11, 220–228. doi:10.1038/nrm2858

PubMed Abstract | CrossRef Full Text | Google Scholar

Ni, W., Yao, S., Zhou, Y., Liu, Y., Huang, P., Zhou, A., et al. (2019). Long Noncoding RNA GAS5 Inhibits Progression of Colorectal Cancer by Interacting with and Triggering YAP Phosphorylation and Degradation and Is Negatively Regulated by the m6A Reader YTHDF3. Mol. Cancer 18, 143. doi:10.1186/s12943-019-1079-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Petti, E., Buemi, V., Zappone, A., Schillaci, O., Broccia, P. V., Dinami, R., et al. (2019). SFPQ and NONO Suppress RNA:DNA-hybrid-related Telomere Instability. Nat. Commun. 10, 1001. doi:10.1038/s41467-019-08863-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahin, I. H., Lowery, M. A., Stadler, Z. K., Salo-Mullen, E., Iacobuzio-Donahue, C. A., Kelsen, D. P., et al. (2016). Genomic Instability in Pancreatic Adenocarcinoma: a New Step towards Precision Medicine and Novel Therapeutic Approaches. Expert Rev. Gastroenterol. Hepatol. 10, 1–13. doi:10.1586/17474124.2016.1153424

PubMed Abstract | CrossRef Full Text | Google Scholar

Seton-Rogers, S. (2018). The Sting of Metastasis. Nat. Rev. Cancer 18, 137. doi:10.1038/nrc.2018.16

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheffer, M., Bacolod, M. D., Zuk, O., Giardina, S. F., Pincas, H., Barany, F., et al. (2009). Association of Survival and Disease Progression with Chromosomal Instability: a Genomic Exploration of Colorectal Cancer. Pnas 106, 7131–7136. doi:10.1073/pnas.0902232106

PubMed Abstract | CrossRef Full Text | Google Scholar

Shukla, S., Evans, J. R., Malik, R., Feng, F. Y., Dhanasekaran, S. M., Cao, X., et al. (2017). Development of a RNA-Seq Based Prognostic Signature in Lung AdenocarcinomaJNCI J. Natl. Cancer Inst. 109, djw200. doi:10.1093/jnci/djw200

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, J., Zhang, Z., Bao, S., Yan, C., Hou, P., Wu, N., et al. (2020). Identification of Tumor Immune Infiltration-Associated lncRNAs for Improving Prognosis and Immunotherapy Response of Patients with Non-small Cell Lung Cancer. J. Immunother. Cancer 8, 1–12. doi:10.1136/10.1136/jitc-2019-000110

CrossRef Full Text | Google Scholar

Tani, H., Numajiri, A., Aoki, M., Umemura, T., and Nakazato, T. (2019). Short-lived Long Noncoding RNAs as Surrogate Indicators for Chemical Stress in HepG2 Cells and Their Degradation by Nuclear RNases. Sci. Rep. 9, 20299. doi:10.1038/s41598-019-56869-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Torre, L. A., Bray, F., Siegel, R. L., Ferlay, J., Lortet-Tieulent, J., and Jemal, A. (2015). Global Cancer Statistics, 2012. CA: A Cancer J. Clinicians 65, 87–108. doi:10.3322/caac.21262

CrossRef Full Text | Google Scholar

Tracy, K. M., Tye, C. E., Ghule, P. N., Malaby, H. L. H., Stumpff, J., Stein, J. L., et al. (2018). Mitotically-Associated lncRNA (MANCR) Affects Genomic Stability and Cell Division in Aggressive Breast Cancer. Mol. Cancer Res. 16, 587–598. doi:10.1158/1541-7786.MCR-17-0548

PubMed Abstract | CrossRef Full Text | Google Scholar

Tu, Z., Wu, L., Wang, P., Hu, Q., Tao, C., Li, K., et al. (2020). N6-Methylandenosine-Related lncRNAs Are Potential Biomarkers for Predicting the Overall Survival of Lower-Grade Glioma Patients. Front. Cel Dev. Biol. 8, 642. doi:10.3389/fcell.2020.00642

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Hui, H., Agrawal, K., Kang, Y., Li, N., Tang, R., et al. (2020). M 6 A Rna Methyltransferases Mettl3/14 Regulate Immune Responses To Anti‐Pd‐1 Therapya Rna Methyltransferases Mettl3/14 Regulate Immune Responses To Anti-Pd-1 Therapy. EMBO J. 39, e104514. doi:10.15252/embj.2020104514

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Luo, X., Cheng, C., Amos, C. I., Cai, G., and Xiao, F. (2020). A Gene Expression-Based Immune Signature for Lung Adenocarcinoma Prognosis. Cancer Immunol. Immunother. 69, 1881–1890. doi:10.1007/s00262-020-02595-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Lu, J.-H., Wu, Q.-N., Jin, Y., Wang, D.-S., Chen, Y.-X., et al. (2019). Lncrna Linris Stabilizes Igf2bp2 And Promotes The Aerobic Glycolysis In Colorectal Cancer. Mol. Cancer 18, 174. doi:10.1186/s12943-019-1105-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y., Yang, X., Chen, Z., Tian, L., Jiang, G., Chen, F., et al. (2019). m6A-induced lncRNA RP11 Triggers the Dissemination of Colorectal Cancer Cells via Upregulation of Zeb1. Mol. Cancer 18, 87. doi:10.1186/s12943-019-1014-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Zhang, S., He, C., Xue, P., Zhang, L., He, Z., et al. (2020). METTL14 Suppresses Proliferation and Metastasis of Colorectal Cancer by Down-Regulating Oncogenic Long Non-coding RNA XIST. Mol. Cancer 19, 46. doi:10.1186/s12943-020-1146-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z., Li, H., Wang, Z., Yang, Y., Niu, J., Liu, Y., et al. (2018). Microarray Expression Profile of Long Non-coding RNAs in Human Lung Adenocarcinoma. Thorac. Cancer 9, 1312–1322. doi:10.1111/1759-7714.12845

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Zhang, M., Ge, S., Huang, W., Lin, X., Gao, J., et al. (2019). Reduced m6A Modification Predicts Malignant Phenotypes and Augmented Wnt/PI3K‐Akt Signaling in Gastric Cancer. Cancer Med. 8, 4766–4781. doi:10.1002/cam4.2360

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Guo, S., Piao, H.-y., Wang, Y., Wu, Y., Meng, X.-y., et al. (2019). ALKBH5 Promotes Invasion and Metastasis of Gastric Cancer by Decreasing Methylation of the lncRNA NEAT1. J. Physiol. Biochem. 75, 379–389. doi:10.1007/s13105-019-00690-8

CrossRef Full Text | Google Scholar

Zhang, L., Luo, Y., Cheng, T., Chen, J., Yang, H., Wen, X., et al. (2021). Development and Validation of a Prognostic N6-Methyladenosine-Related Immune Gene Signature for Lung Adenocarcinoma. Pgpm Vol 14, 1549–1563. doi:10.2147/PGPM.S332683

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Varn, F. S., Cai, G., Xiao, F., Amos, C. I., and Cheng, C. (2018). A P53-Deficiency Gene Signature Predicts Recurrence Risk of Patients with Early-Stage Lung Adenocarcinoma. Cancer Epidemiol. Biomarkers Prev. 27, 86–95. doi:10.1158/1055-9965.EPI-17-0478

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuo, X., Chen, Z., Gao, W., Zhang, Y., Wang, J., Wang, J., et al. (2020). M6A-mediated Upregulation of LINC00958 Increases Lipogenesis and Acts as a Nanotherapeutic Target in Hepatocellular Carcinoma. J. Hematol. Oncol. 13, 5. doi:10.1186/s13045-019-0839-x

CrossRef Full Text | Google Scholar