Comprehensive analyses of molecular features, prognostic values, and regulatory functionalities of m6A-modified long non-coding RNAs in lung adenocarcinoma
Clinical Epigenetics volume 15, Article number: 60 (2023)
Lung adenocarcinoma (LUAD) has a high incidence and recurrence rate. N6-methyladenosine (m6A) modification of RNA has become a promising epigenetic marker in tumors. The dysregulation of both RNA m6A levels and m6A regulator expression levels reportedly affects essential biological processes in various tumors. Long non-coding RNAs (lncRNAs), a subgroup of RNAs over 200 nucleotides in length that do not code for protein, can be modified and regulated by m6A, but the relevant profile in LUAD remains unclear.
The m6A levels of total RNA were decreased in LUAD tumor tissues and cells. Multiple m6A regulators were abnormally expressed at both the RNA and protein levels, and were related in expression patterns and functionally synergistic. Our microarray revealed 2846 m6A-modified lncRNA transcripts as well as its molecular features, 143 of which were differentially m6A-modified and manifested a negative correlation between expression levels and m6A modification levels. More than half of the differentially m6A-modified lncRNAs associated with dysregulated expression. The 6-MRlncRNA risk signature was a reliable indicator for assessing survival time of LUAD patients. The competitive endogenous regulatory network suggested a potential m6A-induced pathogenicity in LUAD.
These data have demonstrated that differential RNA m6A modification and m6A regulator expression levels were identified in LUAD patients. In addition, this study provides evidence increasing the understanding of molecular features, prognostic values, and regulatory functionalities of m6A-modified lncRNAs in LUAD.
According to GLOBOCAN 2020 statistics, lung cancer was the second most common tumor type and had the highest mortality worldwide . Despite the great improvements in antitumor treatments over the past three decades, the 5-year survival rates for patients with regional and widespread disseminated lung cancer were only 27% and 4%, respectively . Notably, as the most prevalent lung cancer histological subtype, lung adenocarcinoma (LUAD) accounted for about 40% of cases  and has had increasing morbidity in recent years , making it the main research focus in this field. To improve the therapeutic effect and prolong patient survival, novel effective biomarkers are required for the early detection, prognosis, and monitoring of LUAD, as well as to explore the molecular mechanisms underlying this disease.
N6-methyladenosine (m6A), the most prevalent post-transcriptional modification in human mRNAs, was discovered within the consensus motifs DRm6ACH (D = G/A/U, R = G/A, H = A/U/C) and is highly enriched within long internal exons, around stop codons, and in the 3′ untranslated region (3′ UTR) of mRNAs [5, 6]. Such cellular modifications are part of a highly dynamic state that is regulated by three types of RNA m6A regulators, including methyltransferases (METTL3, METTL14, METTL16, VIRMA, RBM15, RBM15B, WTAP, ZC3H13, and CBLL1), demethylases (FTO and ALKBH5), and binding functional proteins (YTHDF1/2/3, YTHDC1/2, HNRNPC/A2B1, RBMX, IGF2BP1/2/3, FMR1, EIF3A, ABCF1, ELAVL1, G3BP1/2, ZNF217, and LRPPRC). These regulator types mediate the occurrence, removal, and functionality of m6A modifications, respectively [7,8,9,10]. Several studies have found that an abnormal abundance of RNA m6A modifications, likely mediated by the mutation or aberrant expression of m6A regulators, is a new epigenetic gene regulatory mechanism in cell differentiation and proliferation in certain tumors [11, 12].
Long non-coding RNAs (lncRNAs), defined as a type of RNA molecule that is more than 200 nucleotides in length and has limited or non-existent protein coding ability, are poorly conserved among species and show dynamic and specific expression patterns in various tissues and cells [13, 14]. Because lncRNAs were found to be generally far lower as abundant as mRNAs in a wide range of human organs, they were initially believed to be a kind of transcriptional byproduct [15, 16]. However, researchers later demonstrated that lncRNAs can significantly impact cell proliferation and differentiation in tumor cells because of their multiple roles in transcriptional and post-transcriptional regulation of certain oncogenes . Currently, the m6A modification was found in many lncRNAs and could support their functionality in regulating the expression levels of tumor-related genes via affecting self-stability , nuclear accumulation , RNA–protein interactions , and affinity of binding microRNAs (miRNAs) . Nevertheless, the m6A-modified lncRNA profile in LUAD remains unclear.
In this study, we first evaluated and confirmed the relationship between the m6A modification and LUAD by detecting the m6A modification levels in total RNA and exploring the m6A regulator expression levels. Importantly, the Arraystar Human M6A-modified LncRNA Epitranscriptomic Microarray was utilized to analyze the characteristics of m6A-modified lncRNAs and screen out lncRNAs with differential methylation level in LUAD. Combined analysis with The Cancer Genome Atlas (TCGA) database, we developed and validated the six m6A-regulated lncRNAs (6-MRlncRNA) risk signature to predict the overall survival (OS) of LUAD patients. At last, we predicted the competitive endogenous regulatory role of differentially m6A-modified lncRNAs in LUAD.
LUAD is associated with low levels of total m6A-modified RNA
To explore the potential role of m6A modification in LUAD, we first examined the m6A levels in total RNA samples. We found that the majority of early-stage LUAD tumor tissues exhibited reduced RNA m6A methylation levels compared with adjacent non-tumor lung tissues (N = 20, P < 0.01) using the EpiQuik™ m6A RNA methylation quantification kit (Fig. 1A). Consistently, the m6A modification abundance of RNA from the six LUAD cell lines was lower (P < 0.05) compared with that from the normal lung epithelial cell BEAS-2B. This was particularly clear in the H838, PC9, and A549 cell lines (Fig. 1B). These results indicated that the RNA m6A modification level was decreased in LUAD tissues and cells.
Disordered expression of m6A regulators in LUAD tissues and cell lines
We hypothesized that RNA m6A modification could be more broadly associated with the altered expression levels of m6A methyltransferases, demethylases, and functional proteins in LUAD. Thus, we evaluated the gene expression levels of 30 reported m6A regulators. Through the conjoint analysis of TCGA and GTEx databases, the transcript expression levels of these regulators were all disordered in LUAD tumor tissues (all P < 0.05, except for METTL14), but all differences were not significant (all |log2FC|< 1) (Fig. 2A). Verification by RT-qPCR indicated that METTL3 and FTO were decreased (P < 0.01), and IGF2BP3 was increased (P < 0.01) in LUAD. METTL14, ALKBH5, and IGF2BP1 didn’t show a significant statistical difference (Fig. 2B). Overall, the disordered expression patterns of six regulators in eight paired clinical LUAD samples showed similar trends as the database analysis results. Next, western blot analysis revealed that the six regulators’ protein levels were all significantly increased in LUAD tumor tissues (Fig. 2C), which were consistent with the results in LUAD cell lines with BEAS-2B or HBE cells as the normal control (Fig. 2D). Only IGF2BP3 showed coincident differences in transcript and protein levels, both of which were elevated in LUAD. These results suggested that the disordered expression patterns of m6A regulators were likely involved in the development of LUAD.
To further exploring the potential relationship of the 30 m6A regulators, we investigated their RNA expression relationship and protein interaction by Pearson coefficient and PPI network analysis, respectively. Beyond METTL16, most m6A regulators had correlations with at least one other regulator at the transcript level in lung tissues (|r|> 0.3, P < 0.001), and the strongest negative and positive co-expressions were METTL3 and RBM15B (r = − 0.76, P < 0.001), and YTHDF1 and RBM15B (r = 0.82, P < 0.001) (Fig. 2E). Additionally, most of m6A regulators had frequent protein interactions with each other, especially with WTAP, METTL14, ALKBH5, and HNRNPC (Fig. 2F), suggesting a cooperative interaction pattern in LUAD.
Overview of m6A-modified LncRNA microarray in LUAD tissues
The above experiments demonstrated that m6A regulator expression levels were abnormal in LUAD, making us speculate that lncRNAs modulated by these m6A regulators may also have characteristics changes. To examine the m6A-modified lncRNA profile, we performed the Arraystar Human m6A-lncRNA Epitranscriptomic Microarray on six paired LUAD tumor tissues and adjacent normal tissues. In 12,496 lncRNA-specific probes, the microarray hybridized to and identified a total of 2846 m6A-modified lncRNA transcripts, whose source distribution were mainly in exon sense-overlapping (48%), intergenic (24%), natural antisense (17%), and bidirectional (13%) (Fig. 3A). The length distributions of 2142 m6A-modified lncRNA transcripts (≤ 3000 bp) were mostly 500 to 1000 bp (29%), while there was no significant proportion difference in the other length ranges, which were all 10% to 17% (Fig. 3B). The m6A-modified lncRNA transcripts were localized to all chromosomes, with chromosomes 1 and Y having the highest (13.5%) and lowest (0.5%) abundance, respectively (Fig. 3C). In addition, the amounts of transcripts originating from sense and antisense strands had no apparent difference (Fig. 3C).
To improve biological repetition and eliminate intra-group error, we more specifically examined the data of samples with great heterogeneity (Pearson.r < 0.9) according to the sample correlation analysis of lncRNA expression levels before statistical analysis (Additional file 1: Fig. S1). Finally, the data of Tumor_2/6 were omitted, and the data of other four tumor tissues (Tumor_1/3/4/5) and six adjacent normal tissues (Normal_1/2/3/4/5/6) were used for follow-up difference analysis. Collectively, 123 hypermethylated and 20 hypomethylated lncRNA transcripts had different m6A methylation levels in LUAD tumor tissues (|log2FC|> 0.585, P < 0.05) (Fig. 3D and Additional file 2: Table S1), which were displayed in the clustering heatmap (Fig. 3E). We observed that most changes in m6A methylation levels were within a twofold difference. Subsequently, the expression analysis of 143 differentially modified lncRNA transcripts indicated that 32 were upregulated and 48 were downregulated (|log2FC|> 1, P < 0.05) (Fig. 3F, G and Additional file 2: Table S2). These expression differences ranged from twofold to 37-fold, much higher than the observed fold changes for the m6A modification level. The additional 63 differentially modified lncRNA transcripts did not show significantly differential expression levels (|log2FC|< 1 or P > 0.05) (Fig. 3F, G and Additional file 2: Table S2).
Combined analysis of the m6A modification profile and expression profile in LUAD tissues
Because epigenetic modifications often affect expression levels of the gene, we speculated that the difference in lncRNA expression levels may be caused by its altered m6A modification level. To investigate these, the UpSet diagram was exploited to divide the 143 differentially modified lncRNA transcripts into six categories on the basis of the alteration of both modification and expression (Fig. 4A and Additional file 2: Table S3). We found that “hyper–not different” was the most numerous category (N = 54), followed by “hyper–down” (N = 42), “hyper–up” (N = 27), “hypo–not different” (N = 9), “hypo–down” (N = 6) and “hypo–up” (N = 5). Next, we performed Spearman correlation analysis between the expression levels and m6A levels of the 143 lncRNA transcripts in normal tissues, tumor tissues, and both in normal and tumor tissues, respectively. The results showed a negative relationship, suggesting that in LUAD, a differentially m6A-modified lncRNA with a higher m6A methylation level may have a lower expression level (Fig. 4B). Subsequently, four differentially m6A-modified lncRNAs, SOCAR, PCAT19, SNHG8 and COLCA1, were selected to evaluate the accuracy of the microarray result. As shown in Fig. 4C, RT-qPCR on clinical tissue samples revealed enhanced expression of SNHG8, decreased expression of SOCAR and PCAT19, and unaltered expression of COLCA1 in LUAD. Additionally, MeRIP-qPCR assays validated that SOCAR and SNHG8 transcripts were hypermethylated, and PCAT19 and COLCA1 transcripts were hypomethylated in LUAD tumor tissues (Fig. 4D). The experimental validation results of SNHG8 and SOCAR were consistent with the results in the microarray, which indicated that the microarray had accuracy and its analyses had reference significance. Overall, these results suggested that disordered expression patterns of lncRNA associated with abnormally altered m6A modification on lncRNA.
Construction and evaluation of the m6A-regulated LncRNA prognostic signature
Through the co-expression analysis between the 2846 m6A-modified lncRNAs and 30 m6A regulators from the criteria of |Pearson.r| > 0.3 and P < 0.001 in TCGA-LUAD set, 215 lncRNAs were defined as m6A-regulated lncRNAs (Fig. 5A and Additional file 2: Table S4). Then, we recognized that 13 m6A-regulated lncRNAs were correlated with the OS of patients in the TCGA-LUAD set by univariate Cox regression analysis (P < 0.01) (Fig. 5B). The co-expressions of 13 m6A-regulated lncRNAs and m6A regulators were showed in Fig. 5C. The differences of clinical characteristics between the training and validation sets were ruled out (Table 1), implying that the random division of TCGA-LUAD samples had no effect on subsequent analysis. Next, LASSO Cox analysis and lambda value were used to eliminate 13 m6A-regulated lncRNAs that were highly correlated with each other to avoid overfitting in the training set (Fig. 6A, B). Finally, this produced a m6A-regulated lncRNA prognostic signature containing six lncRNAs and the coefficient of each (Fig. 6C). This used the following formula: risk score = 0.109 * AL590666.2 + (− 0.024) * CH17-340M24.3 + 0.291 * MIR31HG + (− 0.493) * MIR99AHG + (− 0.186) * LINC01936 + 0.111 * LINC02802. Then, RT-qPCR assays on LUAD clinical tissue samples were used to preliminarily explore the expression levels of the six lncRNAs in the signature. As shown in Fig. 6D, the expression of AL590666.2, a risk factor of OS (HR = 1.232, 95% CI: 1.104–1.375), was increased in LUAD tumor tissues (P < 0.05), while the expression of LINC01936, a protective factor (HR = 0.721, 95% CI: 0.567–0.917), was decreased in LUAD tumor tissues compared with corresponding normal tissues (P < 0.05). The expression levels of other prognostic m6A-regulated lncRNAs had no change (all P > 0.05).
Subsequently, for evaluating the performance of the prognostic risk-related signature, the TCGA-LUAD patients were assigned to low- and high-risk subgroups according to the median value of risk scores of the training set. KM survival curves depicted that the survival probability of the high-risk subgroup was lower (P < 0.001 in training set; P = 0.005 in validation set) (Fig. 6E, H). ROC curves showed that the AUC values of the 1-, 2-, and 3-year OS prediction for the signature were all greater than 0.66 (Fig. 6F, I). In addition, the risk score distribution plots and scatter plots showed that the high-risk subgroup had shorter survival times than the low-risk subgroup (upper and middle parts of Fig. 6G, J). The cluster heatmaps showed expression differences of the six prognosis-related m6A-regulated lncRNAs between the two subgroups (lower part of Fig. 6G, J).
Independence assessment of the signature and stratification analysis based on clinical features
As revealed by univariate Cox regression analyses, risk score predicted undesirable OS (HR = 2.694, 95% CI: 1.999–3.630 in training set; HR = 1.412, 95% CI: 1.091–1.828 in validation set; and all P < 0.01) as well as AJCC stage (HR = 1.712, 95% CI: 1.386–2.115 in training set; HR = 1.566, 95% CI: 1.296–1.891 in validation set; and all P < 0.001) (left parts of Fig. 7A, B). Moreover, the multivariate analyses simultaneously confirmed the independence of our constructed risk score signature for predicting the prognosis of LUAD patients (HR = 2.353, 95% CI: 1.710–3.239 in training set; HR = 1.362, 95% CI: 1.064–1.743 in validation set; and all P < 0.05) as well as AJCC stage (HR = 1.530, 95% CI: 1.228–1.907 in training set; HR = 1.567, 95% CI: 1.290–1.902 in validation set; and all P < 0.001) (right parts of Fig. 7A, B).
We employed stratification analysis to confirm the prognostic value of the 6-MRlncRNA prognostic signature. The 482 TCGA-LUAD set patients with complete required clinical features were stratified into different groups based on age (≥ 65 and < 65 years old), gender (female and male) and tumor stage (I-II and III-IV), and were classified into high- and low-risk subgroups in each group according to the median risk score. As expected, the KM survival curves showed that the high-risk subgroup had worse OS compared with the low-risk subgroup in all strata of clinical characteristics (P < 0.05; Fig. C–H), suggesting that the risk score could accurately distinguish the prognosis of the populations with different characteristics. Overall, the stratification analyses indicated that the lncRNA prognostic signature was effective.
Construction of the m6A-induced CeRNA network
In the 43 recorded cellular localizations of 143 differentially m6A-modified lncRNAs based on RNALocate, 40 lncRNAs were detected in extracellular exosomes. Four lncRNAs were present only in the cytoplasm, and 12 were found only in the nucleus, and 12 were present in both the cytoplasm and nucleus, the names of which were separately showed in Fig. 8A. All in all, a total of 16 lncRNAs were in the cytoplasm. Next, we used miRcode to identify 201 highly conserved miRNAs interacting with 16 cytoplasmic lncRNAs. Subsequently, three databases combined to predict 1269 miRNAs-targeted mRNAs. Simultaneously, a total of 4702 differentially expressed mRNAs (|log2FC|> 1, P < 0.05) in LUAD were screened out by the combined analysis of TCGA and GTEx. Ultimately, after taking the intersection, 11 correlative m6A regulators, seven differentially m6A-modified lncRNAs, 30 sponged miRNAs, and 110 targeted mRNAs formed the ceRNA network (Fig. 8B, C). Furthermore, we listed 11 mRNAs with the number of connecting nodes greater than or equal to 3 in the regulatory network (Fig. 8D), which had a high probability to be regulated by the ceRNA network induced by differential m6A modification on lncRNA in LUAD. Simultaneously, we labeled the dysregulated expression levels of 11 mRNAs in LUAD tumor tissues, which had significant fold changes (|log2FC|> 1, P < 0.001) (Fig. 8D).
After studies on histone acetylation and DNA 5-methylcytosine modification, RNA m6A methylation has become another epigenetics focus that has attracted much attention in cancer research in recent years. The overall RNA m6A modification abundances in tumors were reportedly frequently abnormal. Ma et al. observed that the m6A% level of total RNA displayed a gradient descending trend in normal liver tissues, adjacent tissues, and tumor tissues of patients with hepatocellular carcinoma . Ge et al. validated that the m6A% in total RNA of gastric tumor cells was significantly higher than that in normal gastric epithelial cells . Here, we observed a decreased tendency of m6A modification abundance in LUAD tumor tissues at stage I and in most LUAD cell lines. These results indicated that m6A was vital for tumorigenesis in early-stage LUAD. Using statistical analysis and validation experiments, we found that multiple m6A regulators were abnormally expressed at both the RNA and protein levels in LUAD tissues, which may lead to the dysregulation of the m6A%. Furthermore, a majority of m6A regulators were related in expression patterns and functionally synergistic in LUAD, which were consistent with other reported tumor types [24, 25]. Overall, the m6A modification level was dysregulated and potentially involved in the occurrence and development of LUAD.
In recent years, studies have demonstrated that the dysregulation of m6A modifications on lncRNAs mediated by m6A regulators could change the lncRNA expression levels and further facilitate tumor malignant progression. For example, in lung cancer, the stability of MALAT1 was increased by its hyper m6A modification level induced by METTL3. Then, the upregulated MALAT1 could bind miR-1914-3p to further regulate the expression of YAP, a co-transcription factor upregulating chemo-resistant genes . In gastric cancer, METTL14-mediated m6A modification led to LINC01320 upregulation, which promoted an aggressive phenotype of increased cancer cell proliferation, migration, and invasion via regulating the miR-495-5p/RAB19 axis . However, relevant studies on the m6A-modified lncRNA profile in LUAD are still in a preliminary stage and are not fully comprehensive. Among the 2846 m6A-modified lncRNA transcripts detected in our epitranscriptomic microarray, we observed that the m6A-modified lncRNA with the highest proportion in LUAD tumor tissues had the following characteristics: exon sense-overlapping and intergenic transcript type, derived from chromosome 1 or 2, and was within the 500 to 1000 bp range. Among the 143 differentially m6A-modified lncRNA transcripts, about 55% concurrently associated with an altered expression level, which may be worthy of more in-depth studies of functional mechanisms. This phenomenon was confirmed by the experimental validation of PCAT19, SOCAR and SNHG8. According to the fold change analyses, we inferred that in LUAD, small variations of m6A modification levels on lncRNAs may cause a significant alteration of its expression levels. Additionally, we noted that in both normal and tumor tissues, differentially modified lncRNAs with higher m6A modification levels had lower expression levels. Overall, this study screened out specific m6A-modified lncRNAs and complemented existing research about the relationship between m6A modification and expression levels.
Accumulating evidence indicated that m6A modification was an important activator or stabilizer for lncRNAs [26, 28,29,30]. Therefore, we considered that the functionality of m6A-modified lncRNA may be stronger than that of its unmodified form, and was more likely to be an accurate and effective tumor biomarker. We identified the m6A-modified lncRNAs co-expressed with m6A regulators as the “m6A-regulated lncRNAs”, from which we constructed and verified a new prognostic signature that could successfully predict the OS of LUAD patients from their risk scores. The signature was subsequently demonstrated to be an independent prognostic factor, have favorable stability in different clinicopathological subgroups, and be a reliable metric to assess survival time of LUAD patients. Of the six m6A-regulated lncRNAs included in the signature, we verified that AL590666.2 was increased and LINC01936 was decreased in LUAD tumor tissues, while the expression levels of the other lncRNAs were not statistically different. AL590666.2 and CH17-340M24.3 had not been previously reported in tumor. MIR31HG participated in an immune-related lncRNA signature to predict survival and immunotherapy of LUAD patients . In addition, MIR31HG targeted HIF1A and P21 to facilitate head and neck cancer cell proliferation and tumorigenesis by promoting cell cycle progression . In pancreatic tumors, aberrantly upregulated MIR99AHG could modulate notch receptor 2 (NOTCH2) expression and stimulate the Notch signaling pathway to accelerate malignant progression through sponging miR-3129-5p and recruiting ELAV like RNA binding protein 1 (ELAVL1) . In LUAD, LINC01936 displayed adequate performance in distinguishing LUAD patients from healthy people (AUC of ROC = 95.3%). Increased expression levels of LINC01936 were strongly related to a decreased risk of death . LINC02802 was identified as a member of a prognostic signature for cervical tumors .
M6A modifications not only affected the lncRNA expression levels, but also directly modulated the ceRNA model of lncRNAs. A recent study demonstrated that decreased m6A modification level of LINC1281 attenuated its interaction with let7, suggesting that the m6A modification was necessary for the LINC1281-mediated ceRNA model . Similarly, the ZFAS1-miR-647 interaction was regulated by METLL3-mediated m6A modification on ZFAS1 . To explore the underlying mechanism, we constructed a ceRNA model originating from differentially m6A-modified lncRNAs localized in the cytoplasm, the axis of which included 11 m6A regulators, seven lncRNAs, 30 miRNAs, and 110 abnormally expressed mRNAs. 11 mRNAs, SOX4, MAP3K9, MATR3, TFAP2C, RRM2, SALL3, E2F2, C17orf49, MTMR3, SEMA6A and SIK1, were most likely modulated by hyper or hypo m6A modification levels on lncRNAs in LUAD tumor tissues. Interestingly, the subcellular localization results showed that most differentially m6A-modified lncRNAs in LUAD were observed in exosomes, a type of extracellular vesicles. In line with the hypothesis put forward by other researchers , m6A modifications may favor specific packaging of lncRNAs into exosomes. Given that m6A can regulate many RNA metabolic processes, including stability, translation, splicing, and structure switch , it would be meaningful to investigate if exosomes serve as intercellular transport carriers of tumor-associated m6A-modified RNAs in the transformation of target cells from normal to malignant.
The present study had some limitations. Although the microarray analysis revealed specific m6A-modified lncRNAs benefiting from specific hybridization probes, we could not determine the exact sites and motif sequences of m6A modification. Additionally, the overall m6A modification level of a full-length lncRNA may cover up the methylation difference at one specific site because multiple m6A sites were often on one lncRNA transcript. The reliability of the proposed prognostic signature needs to be further evaluated using external data. Furthermore, the m6A-induced ceRNA model requires rigorous molecular experiments to be validated. Even so, our study is meaningful and pioneering, as we have revealed the molecular features, prognostic values, and regulatory functionalities of m6A-modified lncRNAs in LUAD, providing references for the future research in this field.
Clinical tissue sample collection and cell culture
Twenty pairs of tumor and adjacent normal tissues were collected from LUAD patients undergoing surgical treatment in Shanghai Tongji Hospital from 2018 to 2021. The patients were screened to have no history of other malignancies and not receive any preoperative radiotherapy or chemotherapy. The clinicopathological characteristics of the study subjects were summarized in Additional file 2: Table S5. The study was approved by the hospital’s medical ethics committee with written informed consent received from all patients. H1792, H1975, H838, HCC827, and H1944 cell lines were cultured in RPMI-1640 medium (Hyclone, South Logan, UT, USA). A549, PC9, and BEAS-2B cell lines were cultured in DMEM medium (Hyclone, South Logan, UT, USA). HBE cells were cultured in Keratinocyte medium (ScienCell, San Diego, CA, USA). All medium contained 10% fetal bovine serum (FBS; Gibco, Grand Island, NY, USA) and all cells were incubated at 37 °C and 5% CO2.
RNA extraction and RT-qPCR
Total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA), and genomic DNA was removed using Recombinant DNase I (Takara, Dalian, China). 5 × Evo M-MLV RT Master Mix (AG, Changsha, China) was used for cDNA synthesis. SYBR® Green Premix Pro Taq HS qPCR Kit II (AG, Changsha, China) was used for qPCRs, which were performed on a QuantStudio® 5 real-time quantitative PCR instrument. The relative expression levels of targeted genes were calculated using the 2(−∆∆CT) method using β-actin or GAPDH as the internal reference gene. The primers were listed in Additional file 2: Table S6.
Quantification of total m6A RNA
The m6A content of approximate 200 ng RNA extracted from the indicated cells or tissues were analyzed using the EpiQuik m6A RNA Methylation Quantification Kit (Colorimetric) (Epigentek, Farmingdale, NY, USA) following the manufacturer’s instructions. The absorbance values were measured at a wavelength of 450 nm using a microplate reader within 15 min. The amount of m6A-modified RNA (ng) was calculated based on the standard curve. The value of m6A% was obtained by dividing the m6A amount by the input RNA amount.
Western blot analysis and antibodies
Total protein was extracted using RIPA lysis buffer and was quantified using a BCA protein assay kit (Beyotime Biotechnology, Shanghai, China). Equal amounts of protein samples were separated by 10% SDS-PAGE and then transferred onto 0.45 μm nitrocellulose filter membranes (Cytiva, Amesham, UK). The membranes were blocked with 1 × protein free rapid blocking buffer (Epizyme Biotech, Shanghai, China) for 20 min at room temperature and then incubated with primary antibodies at 4 °C overnight. After three washes in Tris-buffered saline and Tween 20 (TBST), the membranes were incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies (AP132P, 1:2000, Millipore, Billerica, MA, USA) for 1 h at room temperature. An ECL detection system (Tanon 4600SF, Shanghai, China) was used for visualization. β-Actin served as an internal control. The primary antibodies were as follows: β-actin (4970S, 1:1000, Cell Signaling Technology (CST), Beverly, MA, USA), METTL3 (86132S, 1:2000, CST, Beverly, MA, USA), METTL14 (26158-1-AP, 1:2000, Proteintech, Wuhan, China), ALKBH5 (80283S, 1:1000, CST, Beverly, MA, USA), FTO (27226-1-AP, 1:2000, Proteintech, Wuhan, China), IGF2BP1 (EPR26408-18, 1:2000, Abcam, Cambridge, UK), and IGF2BP3 (14642-1-AP, 1:2000, Proteintech, Wuhan, China).
Acquisition of LUAD expression profiles and expression analysis of targeted genes
We downloaded the gene expression RNA-seq files from Genotype-Tissue Expression (GTEx) and TCGA, as well as their corresponding clinical phenotype and survival files from UCSC Xena (http://xena.ucsc.edu/). The sample GTEx-SUCS-0626-SM-5CHQE with abnormal value of gene expression was removed, and the data of 346 normal lung tissue samples and 526 LUAD tumor tissue samples were obtained. The gene set of two databases was intersected by the Perl program. The expression value of each gene was uniformly defined as log2FPKM+1. The normalize Between Arrays function of limma R package was used to correct the data. Subsequently, genes were categorized into lncRNA and mRNA groups using the human genome annotation file of GRCH38.p13 version. The Wilcox Test was used to analyze the differential expression of targeted genes.
Correlation analysis and string protein interaction network of m6A regulators
Pearson correlation analysis was performed on the expression levels of 30 m6A regulators in 872 lung tissues. The Cor.mtest function was used to calculate the significance of correlation coefficients and P < 0.001 was considered to be statistically significant. The correlation heatmap was plotted by the “corrplot” R package. In addition, we analyzed the String protein interaction network (https://www.string-db.org/) of 30 m6A regulators with a confidence score cutoff greater than 0.6 as the screening condition. The active interaction sources included textmining, experiments, databases, co-expression, neighborhood, and co-occurrence. Cytoscape was utilized to visualize the protein–protein interaction (PPI) network.
M6A-modified LncRNA epitranscriptomic microarray
We performed the Arraystar Human m6A-lncRNA Epitranscriptomic Microarray detection on six paired LUAD tumor tissues and adjacent normal tissues confirmed by pathological examination. The entire experimental procedure performed by Shanghai Aksomics Biotechnology Co., Ltd was shown in Additional file 1: Fig. S2. The m6A methylation levels and expression levels of lncRNA transcripts were calculated according to the following formulas:
Three µg total RNA was added to 60 µL 5 × IP buffer (50 mM Tris HCl pH 7.4, 750 mM NaCl, 0.5% IGEPAL CA-630) containing 2 µg affinity purified anti-m6A rabbit polyclonal antibody (Synaptic Systems, Goettingen, Germany), and then incubated in a rotary shaker for 2 h at 4 °C. A rabbit anti-human IgG antibody (Abcam, Cambridge, UK) was used as negative control. The samples were then mixed with 20 µL sheep anti-rabbit IgG magnetic beads (Invitrogen, Carlsbad, CA, USA) that were blocked in advance with 0.5% BSA for 2 h. The mixture was then incubated at 4 °C for 14 h. The beads were washed three times with 300 µL 1 × IP buffer and twice with 300 µL 1 × wash buffer (100 mM Tris HCl pH 7.4, 50 mM NaCl, 0.1% IGEPAL CA-630). The m6A-modified RNA was eluted from the magnetic beads using 300 µL elution buffer (100 mM Tris HCl pH 7.4, 1 mM EDTA, 0.05% SDS, 4 µL proteinase K, 2 µL RNase inhibitor) in the rotary shaker for 1 h at 50 °C. The RNA was extracted by phenol/chloroform/isoamylol (25:24:1) reagent (Tinadz, Beijing, China). The input RNA and m6A-modified RNA were both analyzed via RT-qPCR in a 1:1 ratio. The primers are listed in Additional file 2: Table S6. The m6A methylation level of lncRNA was calculated using the following formula:
Construction and validation of m6A-regulated LncRNA prognostic signature
According to previous publications, 30 proven m6A regulators were included, including methyltransferases, demethylases, and binding functional proteins. Through the Pearson correlation analysis (|r|> 0.3 and P < 0.001) of expression level, a total of 2017 lncRNAs related to m6A regulators were verified, including 215 m6A-modified lncRNAs screened by m6A microarray. Samples without survival values were removed out, leaving 500 TCGA-LUAD tumor samples to be included in the subsequent survival analysis. Next, 500 samples were randomly divided into a training set (N = 252) and validation set (N = 248) for the construction and validation of the prognostic signature, respectively. Later, the least absolute shrinkage and selection operator (LASSO) regression analysis was applied to eliminate those prognostic-associated lncRNAs highly correlated with each other to avoid overfitting. The risk score of LUAD patients was calculated using the following formula:
The LUAD patients were classified into high-risk and low-risk subgroups using the median risk score as the cutoff value. Kaplan–Meier (KM) survival curve analysis was performed to compare the survival outcomes of the two subgroups. The receiver operating characteristic curve (ROC) and its area under the curve (AUC) values were utilized to evaluate the specificity and sensitivity of the signature using the “timeROC” R package. The information and risk score for each sample in the training and validation set were supplemented in the Additional file 2: Table S7.
Independence assessment of the signature and stratification analysis
We excluded the TCGA-LUAD tumor samples with unknown clinical variables of age, gender and AJCC stage. A total of 482 tumor samples (Additional file 2: Table S7) were included to evaluate the independence of the signature using univariate and multivariate Cox regression analyses, and to assess the signature’s ability to predict OS in different groups of age (≥ 65 and < 65 years old), gender (female and male) and AJCC tumor stage (I-II and III-IV).
Construction of m6A-induced competitive endogenous RNA (CeRNA) network
To clarify the competitive endogenous regulation of differentially m6A-modified lncRNAs, we firstly analyzed their cellular localization using the RNALocate v2.0 database (http://www.rna-society.org/rnalocate/). Next, miRcode (http://www.mircode.org/) was used to identify highly conserved miRNAs interacting with lncRNAs. Subsequently, miRNAs-targeted mRNAs were predicted by three databases, which were TargetScan (http://www.targetscan.org/), miRDB (http://www.mirdb.org/), and miRTarBase (http://mirtarbase.mbc.nctu.edu.tw). Only mRNAs predicted by all three databases were eligible. Simultaneously, differentially expressed mRNAs in LUAD were screened out with the combined analysis of TCGA and GTEx. The combination of lncRNA–m6A regulator protein detected by CLIP-seq was recorded in the POSTAR3 database (http://184.108.40.206/RBP.html). Cytoscape was utilized to visualize the ceRNA regulatory relationship.
Availability of data and materials
All data in the current study are available from the corresponding author on reasonable request.
Long non-coding RNAs
Receiver operating characteristic curve
Area under the curve
The Cancer Genome Atlas
American Joint Committee on Cancer
False discovery rate
Competitive endogenous RNA
Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71:209–49. https://doi.org/10.3322/caac.21660.
Dong N, Shi L, Wang DC, Chen C, Wang X. Role of epigenetics in lung cancer heterogeneity and clinical implication. Semin Cell Dev Biol. 2017;64:18–25. https://doi.org/10.1016/j.semcdb.2016.08.029.
Denisenko TV, Budkevich IN, Zhivotovsky B. Cell death-based treatment of lung adenocarcinoma. Cell Death Dis. 2018;9:117. https://doi.org/10.1038/s41419-017-0063-y.
Shi JF, Wang L, Wu N, Li JL, Hui ZG, Liu SM, et al. Clinical characteristics and medical service utilization of lung cancer in China, 2005–2014: overall design and results from a multicenter retrospective epidemiologic survey. Lung Cancer. 2019;128:91–100. https://doi.org/10.1016/j.lungcan.2018.11.031.
Linder B, Grozhik AV, Olarerin-George AO, Meydan C, Mason CE, Jaffrey SR. Single-nucleotide-resolution mapping of m6A and m6Am throughout the transcriptome. Nat Methods. 2015;12:767–72. https://doi.org/10.1038/nmeth.3453.
Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature. 2012;485:201–6. https://doi.org/10.1038/nature11112.
Khan R, Malla WA. m(6)A modification of RNA and its role in cancer, with a special focus on lung cancer. Genomics. 2021;113:2860–9. https://doi.org/10.1016/j.ygeno.2021.06.013.
Chokkalla AK, Mehta SL, Vemuganti R. Epitranscriptomic regulation by m(6)A RNA methylation in brain development and diseases. J Cereb Blood Flow Metab. 2020;40:2331–49. https://doi.org/10.1177/0271678X20960033.
Wang S, Chai P, Jia R, Jia R. Novel insights on m(6)A RNA methylation in tumorigenesis: a double-edged sword. Mol Cancer. 2018;17:101. https://doi.org/10.1186/s12943-018-0847-4.
Fazi F, Fatica A. Interplay between N (6)-methyladenosine (m(6)A) and non-coding RNAs in cell development and cancer. Front Cell Dev Biol. 2019;7:116. https://doi.org/10.3389/fcell.2019.00116.
Liu J, Eckert MA, Harada BT, Liu SM, Lu Z, Yu K, et al. m(6)A mRNA methylation regulates AKT activity to promote the proliferation and tumorigenicity of endometrial cancer. Nat Cell Biol. 2018;20:1074–83. https://doi.org/10.1038/s41556-018-0174-4.
Lan T, Li H, Zhang D, Xu L, Liu H, Hao X, et al. KIAA1429 contributes to liver cancer progression through N6-methyladenosine-dependent post-transcriptional modification of GATA3. Mol Cancer. 2019;18:186. https://doi.org/10.1186/s12943-019-1106-z.
Peng WX, Koirala P, Mo YY. LncRNA-mediated regulation of cell signaling in cancer. Oncogene. 2017;36:5661–7. https://doi.org/10.1038/onc.2017.184.
Jiang C, Li Y, Zhao Z, Lu J, Chen H, Ding N, et al. Identifying and functionally characterizing tissue-specific and ubiquitously expressed human lncRNAs. Oncotarget. 2016;7:7120–33. https://doi.org/10.18632/oncotarget.6859.
Derrien T, Johnson R, Bussotti G, Tanzer A, Djebali S, Tilgner H, et al. The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genome Res. 2012;22:1775–89. https://doi.org/10.1101/gr.132159.111.
Deveson IW, Hardwick SA, Mercer TR, Mattick JS. The dimensions, dynamics, and relevance of the mammalian noncoding transcriptome. Trends Genet. 2017;33:464–78. https://doi.org/10.1016/j.tig.2017.04.004.
Fang Y, Fullwood MJ. Roles, functions, and mechanisms of long non-coding RNAs in cancer. Genomics Proteomics Bioinform. 2016;14:42–54. https://doi.org/10.1016/j.gpb.2015.09.006.
Zuo X, Chen Z, Gao W, Zhang Y, Wang J, Wang J, et al. M6A-mediated upregulation of LINC00958 increases lipogenesis and acts as a nanotherapeutic target in hepatocellular carcinoma. J Hematol Oncol. 2020;13:5. https://doi.org/10.1186/s13045-019-0839-x.
Wu Y, Yang X, Chen Z, Tian L, Jiang G, Chen F, et al. m6A-induced lncRNA RP11 triggers the dissemination of colorectal cancer cells via upregulation of Zeb1. Mol Cancer. 2019;18:87. https://doi.org/10.1186/s12943-019-1014-2.
Liu N, Dai Q, Zheng G, He C, Parisien M, Pan T. N(6)-methyladenosine-dependent RNA structural switches regulate RNA-protein interactions. Nature. 2015;518:560–4. https://doi.org/10.1038/nature14234.
Yang Z, Ma J, Han S, Li X, Guo H, Liu D. ZFAS1 exerts an oncogenic role via suppressing miR-647 in an m(6)A-dependent manner in cervical cancer. Onco Targets Ther. 2020;13:11795–806. https://doi.org/10.2147/OTT.S274492.
Ma JZ, Yang F, Zhou CC, Liu F, Yuan JH, Wang F, et al. METTL14 suppresses the metastatic potential of hepatocellular carcinoma by modulating N(6) -methyladenosine-dependent primary MicroRNA processing. Hepatology. 2017;65:529–43. https://doi.org/10.1002/hep.28885.
Ge L, Zhang N, Chen Z, Song J, Wu Y, Li Z, et al. Level of N6-methyladenosine in peripheral blood RNA: a novel predictive biomarker for gastric cancer. Clin Chem. 2020;66:342–51. https://doi.org/10.1093/clinchem/hvz004.
Chai RC, Wu F, Wang QX, Zhang S, Zhang KN, Liu YQ, et al. m(6)A RNA methylation regulators contribute to malignant progression and have clinical prognostic impact in gliomas. Aging. 2019;11:1204–25. https://doi.org/10.18632/aging.101829.
Zhao X, Cui L. Development and validation of a m(6)A RNA methylation regulators-based signature for predicting the prognosis of head and neck squamous cell carcinoma. Am J Cancer Res. 2019;9:2156–69.
Jin D, Guo J, Wu Y, Du J, Yang L, Wang X, et al. m(6)A mRNA methylation initiated by METTL3 directly promotes YAP translation and increases YAP activity by regulating the MALAT1-miR-1914-3p-YAP axis to induce NSCLC drug resistance and metastasis. J Hematol Oncol. 2019;12:135. https://doi.org/10.1186/s13045-019-0830-6.
Hu N, Ji H. N6-methyladenosine (m6A)-mediated up-regulation of long noncoding RNA LINC01320 promotes the proliferation, migration, and invasion of gastric cancer via miR495-5p/RAB19 axis. Bioengineered. 2021;12:4081–91. https://doi.org/10.1080/21655979.2021.1953210.
Yang D, Qiao J, Wang G, Lan Y, Li G, Guo X, et al. N6-Methyladenosine modification of lincRNA 1281 is critically required for mESC differentiation potential. Nucl Acids Res. 2018;46:3906–20. https://doi.org/10.1093/nar/gky130.
Yang Z, Ma J, Han S, Li X, Guo H, Liu D. ZFAS1 Exerts an oncogenic role via suppressing miR-647 in an m6A-dependent manner in cervical cancer. Onco Targets Ther. 2020;13:11795–806. https://doi.org/10.2147/OTT.S274492.
Zheng ZQ, Li ZX, Zhou GQ, Lin L, Zhang LL, Lv JW, et al. Long noncoding RNA FAM225A promotes nasopharyngeal carcinoma tumorigenesis and metastasis by acting as ceRNA to sponge miR-590-3p/miR-1275 and upregulate ITGB3. Cancer Res. 2019;79:4612–26. https://doi.org/10.1158/0008-5472.CAN-19-0799.
Zhuang J, Chen Z, Chen Z, Chen J, Liu M, Xu X, et al. Construction of an immune-related lncRNA signature pair for predicting oncologic outcomes and the sensitivity of immunosuppressor in treatment of lung adenocarcinoma. Respir Res. 2022;23:123. https://doi.org/10.1186/s12931-022-02043-4.
Wang R, Ma Z, Feng L, Yang Y, Tan C, Shi Q, et al. LncRNA MIR31HG targets HIF1A and P21 to facilitate head and neck cancer cell proliferation and tumorigenesis by promoting cell-cycle progression. Mol Cancer. 2018;17:162. https://doi.org/10.1186/s12943-018-0916-8.
Xu J, Xu W, Yang X, Liu Z, Zhao Y, Sun Q. LncRNA MIR99AHG mediated by FOXA1 modulates NOTCH2/Notch signaling pathway to accelerate pancreatic cancer through sponging miR-3129-5p and recruiting ELAVL1. Cancer Cell Int. 2021;21:674. https://doi.org/10.1186/s12935-021-02189-z.
Hou J, Yao C. Potential prognostic biomarkers of lung adenocarcinoma based on bioinformatic analysis. Biomed Res Int. 2021;2021:8859996. https://doi.org/10.1155/2021/8859996.
Ye J, Chen X, Lu W. Identification and experimental validation of immune-associate lncRNAs for predicting prognosis in cervical cancer. Onco Targets Ther. 2021;14:4721–34. https://doi.org/10.2147/OTT.S322998.
Mateescu B, Kowal EJ, van Balkom BW, Bartel S, Bhattacharyya SN, Buzás EI, et al. Obstacles and opportunities in the functional analysis of extracellular vesicle RNA—an ISEV position paper. J Extracell Vesicles. 2017;6:1286095. https://doi.org/10.1080/20013078.2017.1286095.
Huang H, Weng H, Chen J. m(6)A modification in coding and non-coding RNAs: roles and therapeutic implications in cancer. Cancer Cell. 2020;37:270–88. https://doi.org/10.1016/j.ccell.2020.02.004.
We thank Susan Zunino, PhD, and J. Iacona, PhD, from Liwen Bianji (Edanz) (www.liwenbianji.cn), for editing the English text of a draft of this manuscript.
This work was supported by the National Natural Science Foundation of China (82072362, 81974314, 82272630).
Ethical approval and consent to participate
All procedures related to patients were carried out under the approval of the Internal Review and the Ethics committee of Shanghai Tongji Hospital. Written informed consent was obtained from all patients.
Consent for publication
All authors agreed to publish this manuscript.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Ping, Y., Huang, J., Zhu, J. et al. Comprehensive analyses of molecular features, prognostic values, and regulatory functionalities of m6A-modified long non-coding RNAs in lung adenocarcinoma. Clin Epigenet 15, 60 (2023). https://doi.org/10.1186/s13148-023-01475-z