Classification of Hepatocellular Carcinoma Based on N6-Methylandenosine–Related lncRNAs Profiling

HCC is one of the most common types of malignancies worldwide and the fourth-leading cause of cancer deaths. Thus, there is an urgent need to search for novel targeted therapies in HCC. 186 m6a-related lncRNAs were screened for subsequent analysis. Two distinct m6A modification clusters were identified to be associated with the overall prognosis in TCGA-LIHC based on the m6A-related lncRNAs profiling, followed by univariate Cox regression analysis. In addition, four m6A-related lncRNAs prognostic signatures were developed and validated that could predict the OS of HCC patients, followed by univariate Cox regression, LASSO regression, and multivariate Cox regression analysis. Moreover, four m6A-related lncRNAs were identified to be related to HCC prognosis. ESTIMATE was used to evaluate the stromal score, immune score, ESTIMATE score, and tumor purity of each HCC sample. ssGSEA was performed to identify the enrichment levels of 29 immune signatures in each sample. Finally, quantitative real-time polymerase chain reaction shown that KDM4A-AS1, BACE1-AS, and NRAV expressions were upregulated in HCC patients. We proved that our m6A-related lncRNAs signature had powerful and robust ability for predicting OS of different HCC subgroups.

Introduction

Among primary liver cancer, HCC represents the major histological subtype. According to the global cancer statistics, HCC is one of the most common types of malignancies worldwide and the fourth leading cause of cancer deaths, causing approximately 66,2000 deaths each year worldwide (Venook et al., 2010; Ferlay et al., 2015). The incidence rates of HCC vary globally, but more than 80% of cases occur in low- and middle-resource countries, particularly in East Asia and Africa (Torre et al., 2015; Fitzmaurice et al., 2017; Yang et al., 2019). As one of the countries with the highest incidences of HCC in the world, China accounts for about 50% of newly diagnosed cases and deaths (Yang and Heimbach, 2020). Over the last few decades, intensive investigations and efforts have focused on the role of protein-coding genes in the pathogenesis of HCC to seek appropriate diagnostic markers and therapeutic targets for HCC (Hoshida et al., 2012; Xie et al., 2013). However, on a global scale, effective molecular markers for diagnosis, evaluation, and treatment of advanced HCC are still insufficient. Therefore, there is an urgent need to search for novel targeted therapies in HCC.

m6A is methylation that occurs in the N6-position of adenosine, which is the most prevalent epigenetic internal modification on eukaryotic messenger RNAs (mRNAs) and non-coding RNAs (ncRNAs) (He et al., 2019). m6A is installed by writers (methyltransferases), removed by erasers (demethylases), and recognized by readers (reader proteins) (Batista, 2017; Dai et al., 2018). In physiological processes, m6A regulates RNA metabolism, including translation, splicing, export, and degradation (Liu et al., 2017; Liu and Gregory, 2019). Accumulating evidence showed that m6A was involved in the pathogenesis and progression of cancer by regulating gene expression and affecting cell self-renewal, differentiation, invasion, and apoptosis (Deng et al., 2018; Panneerdoss et al., 2018; Karthiya and Khandelia, 2020; Lan et al., 2021). For instance, Niu et al. reported that FTO was overexpressed and promoted cell proliferation, colony formation, and metastasis in breast cancer by reducing BNIP3 methylation and promoting BNIP3 degradation (Niu et al., 2019). He et al. indicated that ALKBH5 was downregulated in pancreatic cancer with elevated m6A level and decreased expression of KCNK15-AS1, leading to enhanced migration and invasion of pancreatic tumor cells (He et al., 2018a). Meanwhile, increased research has suggested that m6A plays an important role in the occurrence and development of hepatocellular carcinoma, including METTL14 (Ma et al., 2017) and YTHDF2 (Zhong et al., 2019).

LncRNA is a non-coding RNA with a length of more than 200 nucleotides. Recent studies had shown that lncRNA played an important role in the pathophysiology of cancer, such as epigenetic regulation, cell cycle regulation, and cell regulation (Berretta and Morillon, 2009). It has been identified that the dysregulation of lncRNA is closely associated with tumorigenesis, metastasis, prognosis, and diagnosis in HCC (Abbastabar et al., 2018). For instance, overexpression of MALAT-1 could promote the proliferation, migration, and invasion of liver cancer cells by inhibiting the expression of mir-146-5p, thus leading to the recurrence and metastasis of hepatocellular carcinoma (Lai et al., 2012; Li et al., 2017).In addition, Wang et al. demonstrated that UCA1 might act as an endogenous sponge by repressing the expression of FGFR1 and activating an FGFR1/ERK signaling pathway in HCC (Wang et al., 2015). However, the role of m6A modification patterns in dysregulation of lncRNAs in HCC is still unclear.

In this study, we aimed to explore the m6A-related lncRNAs profiling and tumor microenvironment in HCC. Our study revealed two distinct subtypes based on m6A-related lncRNAs profiling and surprisingly found that the tumor microenvironment under these two clusters was highly consistent with overall survival, immunotherapy, and prognosis, respectively, suggesting that m6A-related lncRNAs played a pivotal and robust role in tumor microenvironment characterizations of individual HCC patients. For that, we established a novel m6A-related lncRNAs signature to predict prognosis and guide individual treatment strategies in HCC patients.

Materials and Methods

Data Source and m6A-Related lncRNAs

The mRNA expression files and clinical information of hepatocellular carcinoma were downloaded from TCGA database, including age, gender, tumor grade, pathologic stage, and AJCC-TNM. The annotation file of GRCH38 for long non-coding RNA was downloaded from GENCODE to identify lncRNA and mRNA. Based on the Ensemble IDs of the gene symbols, 14,086 lncRNAs were acquired in TCGA-LIHC. Furthermore, based on previous literature, the expression matrixes of 23 m 6A-related genes were obtained from TCGA-LIHC (Supplementary Table S1). Ultimately, we performed Pearson correlation analysis to obtain 186 m 6A-related lncRNAs for subsequent analysis with |Cor| > 0.5 and p < 0.001.

Clustering, ESTIMATE, and ssGSEA

For TCGA-LIHC, we first performed univariate Cox regression analysis on the 186 m 6A-related lncRNAs to obtain 78 prognostic m6A-related lncRNAs. In addition, based on the 78 m 6A-related lncRNAs, we performed hierarchical clustering using the “ConsensusClusterPlus” package (R implementation, K = 2). ESTIMATE was used to evaluate the stromal score, immune score, ESTIMATE score, and tumor purity of each HCC sample (Yoshihara et al., 2013). ssGSEA was performed to identify the enrichment levels of 29 immune signatures in each sample (Barbie et al., 2009; Hänzelmann et al., 2013). 10-fold cross validation was implemented to evaluate the classification performance with the accuracy and the weighted F-score in TCGA-LIHC.

M6A-Related lncRNAs Risk Signature in HCC

To identify the potential optimal m6A-related lncRNAs prognostic signature, we randomly divided TCGA-LIHC into two parts (training dataset and testing dataset). Followed by univariate Cox regression, LASSO regression, and multivariate Cox regression, we constructed a m6A-related lncRNA prognostic risk signature for HCC patients in the training dataset, which involved four m6A-related lncRNAs. The calculation formula of the risk score was: Risk score = Coef_AC099850.4*Expression_AC099850.4 + Coef_KDM4A-AS1*Expression_KDM4A-AS1+Coef_BACE1-AS*Expression_BACE1-AS + Coef_NRAV*Expression_NRAV, where Coef represented coefficients with the lowest AIC values. A Coef greater than zero meant increased risk and vice versa to protective factors. In addition, we performed ROC curve and AUC values for the risk signature and other clinicopathological features to evaluate the prognostic ability in the training dataset. Similarly, we validated our risk signature in the testing dataset, followed by K-M survival analysis and ROC curve.

Prediction Analysis of m6A-Related lncRNAs Risk Signature

To identify three of the four m6A-related lncRNAs expressions and their associations with HCC prognosis, the HCC tissues’ expressions and OS and DFS plots were downloaded from the GEPIA database. In addition, we put the training dataset and testing dataset together for subsequent analysis. All HCC patients were divided into high- or low-risk groups based on the median of risk scores as a threshold. We evaluated the risk scores in different clinicopathological features, followed by K-M survival analysis.

Additionally, the “limma” package was implemented to uncover the DEGs between the two subgroups. The “Metascape” website was employed for pathway and functional enrichment analysis involving GO Biological Process, KEGG pathway and Reactome Gene sets, and Canonical Pathways (Zhou et al., 2019). Finally, to identify the tumor hallmarkers, Gene Set Enrichment Analysis (GSEA) was performed.

Four m6A-Related lncRNAs Function Prediction

To further identify how m6A-related lncRNAs regulated and influenced the development and progress of HCC, we performed co-expression analysis to predict downstream mRNAs with |Cor|> 0.5 and p < 0.001. 1985 related downstream mRNAs were identified and intersected with the DEGs to earn 97 m 6A-related lncRNAs targets. Cytoscape software 3.7.2 was implemented to visualize the co-expression networks. In addition, the “Metascape” online tool was used to perform functional and pathway enrichment analysis for the 97 m 6A-related lncRNAs targets.

Cell Culture and Samples

To further evaluate the expressions of three of the four m6A-related lncRNAs in cells and tissues, L02, Huh7, HepG2, and Hep3B cells were cultured in DMEM containing 10% FBS. A total of 10 HCC samples and adjacent normal tissues were collected from HCC patients who underwent surgical resection from January 2018 to December 2020 in the Affiliated Hangzhou First People’s Hospital. All patients involved provided written informed consent. This research was approved by the Ethics Committee of Affiliated Hangzhou First People’s Hospital.

For assessment of three of the four m6A-related lncRNAs in cells and tissue, we extracted the total RNA using an RNA-Quick purification kit. A PrimeScriptTM RT reagent kit with gDNA Eraser (Takara, RR047A) was employed to reverse the total RNA into cDNA. Quantitative real-time polymerase chain reaction was conducted using TB Green Premix EX Taq^TM II (Takara, RB820A). M6A-related lncRNAs expressions were analyzed using the 2^−ΔΔC method (Table 1).

TABLE 1

TABLE 1. Primers’ sequences.

Statistical Analysis

All calculations and analysis were implemented using R software 4.0.4 and Perl language. Kaplan–Meier survival analysis and log-rank test were employed to compare overall survival of two subgroups. The “limma” package was used to generate differential expressed genes. The “ConsensusClusterPlus” package was performed for hierarchical clustering analysis. Followed by univariate Cox regression, LASSO regression, and multivariate Cox regression, a m6A-related risk signature was constructed and validated. Independent prognostic factors were identified with univariate Cox regression and multivariate Cox regression analysis. Student’s t-test was performed for statistical comparisons. Values of p < 0.05 were considered statistically significant.

Results

Identification of m6A-Related lncRNAs in TGCA-LIHC

First, we downloaded the expression matrixes of 374 HCC samples and 50 matched normal controls from TCGA database. Using the annotation file downloaded from the GENCODE website, we confirmed 14,806 lncRNAs and 19,604 mRNAs for subsequent analysis. Then the expression matrixes of 23 m 6A-related genes were extracted and analyzed from the TCGA-LIHC cohort. An m6A-related lncRNA was defined as a lncRNA associated with one or more of 23 m 6A-related genes using Pearson correlation analysis with |Cor| > 0.5 and p < 0.001. The workflow of the entire research is shown in Figure 1A, and the associations of 23 m 6A-related genes and 186 lncRNAs are shown in Figure 1B.

FIGURE 1

FIGURE 1. Overview of study. (A) Study flow chart. (B) Network of the correlations between m6A-related genes and 186 m6A-related lncRNAs.

Hierarchical Clustering Identifies Two HCC Subtypes Based on m6A-Related lncRNAs Profiling

For the previously obtained 186 m 6A-related lncRNAs, we performed univariate Cox regression analysis and identified 78 m 6A-related prognostic lncRNAs. On the basis of the 78 m 6A-related lncRNAs, we hierarchically clustered HCC using the “ConsensusClusterPlus” package (R implementation, K = 2). Interesting, HCC patients were clearly separated into two clusters. We defined the two clusters as Cluster 1 (Immune Response Low) and Cluster 2 (Immune Response High), respectively. K-M survival analysis showed that Cluster two was significantly associated with better prognosis (p = 0.012, Figure 2A). In addition, PD-L1 was more highly expressed in tumor and Cluster one than in normal and Cluster 1: PD-L1 might act as a mediator in tumor progression (Figures 2B,C). We performed ESTIMATE to evaluate the stromal score, immune score, ESTIMATE score, and tumor purity of each HCC sample. When comparing the stromal score, immune score, and ESTIMATE score in two clusters, we found scores increasing from Cluster one to Cluster 2 (Figures 2D–F). Tumor purity had the opposite trend (Figure 2G). Above all, these results demonstrated that Cluster two had a higher number of immune cells and stromal cells, while Cluster one contained more tumor cells.

FIGURE 2

FIGURE 2. Hierarchical clustering identifies two HCC subtypes based on m6A-related lncRNAs profiling. (A) Kaplan–Meier survival analysis illustrated that Cluster 2 (Immune Response High) was significantly associated with better prognosis. (B,C) PD-L1 was more highly expressed in tumor and Cluster one than in normal and Cluster 2. (D) Stromal score. (E) Immune score. (F) ESTIMATE score. (G) Tumor purity. (H) Heatmap of the enrichment levels of 29 immune signatures in each sample. *p < 0.05; **p < 0.01; ***p < 0.001.

Moreover, as shown in Figure 2H, ssGSEA was performed to identify the enrichment levels of 29 immune signatures in each sample. 10-fold cross validation was implemented to evaluate the classification performance with an accuracy value of 91% and F-score of 89%. Collectively, these results proved that hierarchical clustering identifies two subtypes significantly associated with the overall prognosis in TCGA-LIHC based on the m6A-related lncRNAs profiling.

Construction and Validation of Four m6A-Related lncRNAs Prognostic Risk Signature in TCGA-LIHC

To construct a potential optimal m6A-related lncRNA to predict the overall survival of HCC patients, we randomly divided TCGA-LIHC into two parts, namely, the training set and the testing set (Supplementary Tables S2,3). For the training set, we constructed a m6A-related lncRNAs prognostic signature based on the 186 m 6A-related lncRNAs, which involved Coefs of each and four m6A-related lncRNAs (Figures 3A,B). We calculated the risk scores of each HCC patient based on coefficients and gene expression (Figure 3C, Supplementary Figure S2; Table 2). We divided the HCC patients of the training set into high- and low-risk subgroups based on the median of risk scores as the threshold. K-M survival analysis showed that HCC patients with high-risk scores suffered shorter overall survival and worse prognosis (Figure 3D). The survival status distribution proved that HCC patients died more as the risk increased (Figure 3E). In addition, we evaluated the prognostic ability of the risk signature compared with age, gender, tumor grade, pathologic stage, and AJCC-TNM. The ROC curve analysis showed that our four m6A-related lncRNA had the optimal prognostic ability to forecast the OS of HCC patients (AUC = 0.758, Figure 3F).

FIGURE 3

FIGURE 3. Construction and validation of four m6A-related lncRNAs prognostic risk signature in TCGA-LIHC. (A,B) LASSO regression was performed, calculating the optimal criteria. (C) Coefficients. (D) Kaplan–Meier survival analysis showed that the high-risk group had poor prognosis and shorter overall survival in the training set. (E) The scatter plot of risk scores in the training set. (F) The scatter plot of survival status in the training set. (G–I) The survival plot, scatter plot of risk scores, and survival status in the testing set.

TABLE 2

TABLE 2. Multivariate Cox regression analysis results during signature construction.

To validate our m6A-related lncRNAs prognostic signature, we calculated the risk scores of each HCC patient in the testing set using the same scheme. The results were similar with the training set: HCC patients with high risk had worse survival outcome with statistical significance (Figure 3G, Supplementary Figure S2). The distributions of risk scores and survival status were shown in Figure 3H. ROC curve analysis demonstrated that the risk score has a stable and robust predictive ability (AUC = 0.787, Figure 3I). All results indicated that our four m6A-related lncRNAs signatures could predict the OS of HCC patients.

Three of Four m6A-Related lncRNAs Expressions and Their Relationships With Prognosis in HCC Patients

To further explore three of the four m6A-related lncRNAs expressions and their associations with prognosis in HCC, the Gepia database was employed to obtain the gene expressions, overall survival, and disease-free survival of three m6a-related lncRNAs in HCC. As shown in Figure 4, compared with the control group, KDM4A-AS1, BACE1-AS, and NRAV were significantly elevated in HCC (Figures 4A–G). Additionally, HCC patients with increased expressions of KDM4A-AS1, BACE1-AS, and NRAV had shorter OS and DFS, as well as worse prognosis (Figures 4B–I). Thus, these results suggested that three of four m6A-related lncRNAs could be implemented as independent biomarkers for predicting prognosis in HCC.

FIGURE 4

FIGURE 4. Evaluating three of the four m6A-related lncRNAs in the GEPIA database. (A–C) KDM4A-AS1. (D–F) BACE1-AS. (G–I) NRAV expression, overall survival plot, and disease-free survival plot from the GEPIA database, respectively.

Association Between Risk Signature and Other Clinicopathological Features in TCGA-LIHC

We attempted to explore the associations between risk signature and other clinicopathological features in TCGA-LIHC. The results indicated that risk scores were associated with AJCC-T, tumor grade, pathologic stage, immune score, and hierarchical clusters (Figures 5A–E). To better evaluate the independent prognostic ability of m6A-related lncRNAs, we conducted a stratified analysis to confirm whether m6A-related lncRNAs could predict OS of different HCC subgroups. Compared with low-risk patients, high-risk HCC patients in various subgroups (age<65, age≥65, male, G1-2, G3, Stage I-II, Stage III-IV, T1-2, T3-4, N0, N1-3 and M0 subgroups) had shorter OS and worse prognosis (Figure 5F–Q).

FIGURE 5

FIGURE 5. Association between risk signature and other clinicopathological features in TCGA-LIHC. (A–E) Risk scores were associated with AJCC-T, tumor grade, pathologic stage, immune score, and hierarchical clusters. (F–Q) Compared with low-risk patients, high-risk HCC patients in various subgroups had shorter OS and worse prognosis.

Moreover, univariate Cox regression analysis indicated that pathologic stage, AJCC-T, and risk score were associated with the prognosis of HCC patients. Multivariate Cox regression analysis proved that risk score is an independent prognostic indicator for HCC patients (Table 3). These data proved that our m6A-related lncRNAs signature had powerful and robust ability for predicting OS of different HCC subgroups.

TABLE 3

TABLE 3. Univariate and multivariate Cox regression analysis of OS in TCGA-LIHC.

Functional Enrichment and GSEA

To identify the potential functional and pathway enrichment involved in gene heterogeneity between high- and low-risk subgroups, R package “limma” was employed to obtain 2054 DEGs, followed by the standards of p < 0.05 and |log2(Foldchange)| > 1. We put 2054 DEGs into the “Metascape” online tool for analysis, which mainly enriched in these terms: Cell Cycle, nuclear division, DNA replication, Retinoblastoma Gene in Cancer, microtubule cytoskeleton organization, and so on (Figures 6A–C). Gene Set Enrichment Analysis (GSEA) was performed to identify the tumor hallmarkers enriched in the high-risk group, such as PI3K-AKT-mTOR signaling, mTORC1 signaling, the P53 pathway, and so on (Figure 6D). These results could provide us with novel insights into molecular biological function related to m6A-related lncRNAs.

FIGURE 6

FIGURE 6. Functional enrichment and Gene Set Enrichment Analysis (GSEA). (A) Heatmap of enriched terms between high- and low-risk groups, colored according to p-value. Network of enriched terms colored according to (B) cluster ID and (C) p-value. (D) Gene Set Enrichment Analysis (GSEA) showed that tumor hallmarkers were enriched in the high-risk group.

Predicting m6A-Related lncRNAs Function Using Co-Expression Analysis

To investigate how m6A-related lncRNAs regulate and influence the development and progress of HCC, we performed co-expression analysis to predict downstream mRNAs with |Cor|> 0.5 and p < 0.001. 1985 related downstream mRNAs were identified and intersected with the DEGs to earn 97 m 6A-related lncRNAs targets (Figure 7A). Metascape was performed for functional and pathway enrichment analysis of 97 mRNAs. The enriched terms were cell division, microtubule cytoskeleton organization involved in mitosis, PID AURORA B PATHWAY, PID PLK1 PATHWAY, and so on. These data mainly provide us with clues to confirm the functions of the four m6A-related lncRNAs.

FIGURE 7

FIGURE 7. Predicting m6A-related lncRNAs function using co-expression analysis. (A) Network of four m6A-related lncRNAs and their target mRNA using co-expression analysis. (B) Heatmap of enriched terms across the 97 target mRNAs. Network of enriched terms across the 97 target mRNAs colored according to (B) cluster ID and (C) p-value.

Validation of Three m6A-Related lncRNAs Expression in Cells and Tissues

For validation of three of the four m6A-related lncRNAs expressions in cells and tissue, we detected three of the four m6A-related lncRNAs expressions in L02, Huh7, HepG2, and Hep3B cells. The results showed that KDM4A-AS1, BACE1-AS, and NRAV expressions were significantly upregulated in HCC cells in contrast to normal liver cells (Figures 8A–C). In addition, we detected three m6A-related lncRNAs expressions in 10 pairs of HCC tissues and adjacent normal tissues we collected by quantitative real-time polymerase chain reaction (Q-PCR). Our results indicated that KDM4A-AS1, BACE1-AS, and NRAV expressions were upregulated in HCC patients (Figures 8D–F).

FIGURE 8

FIGURE 8. Validation of three m6A-related lncRNAs expression in cells and tissues. (A–C) KDM4A-AS1, BACE1-AS, and NRAV expression in normal liver cells and HCC cells. (D–F) KDM4A-AS1, BACE1-AS, and NRAV expression in HCC tissues and adjacent normal tissues.

Discussion

A total of 374 HCC samples from TCGA database were included in our study to explore the m6A-related lncRNAs profiling and tumor microenvironment. 186 m 6A-related lnRNAs were obtained from TCGA-LIHC. Two HCC subtypes were confirmed by hierarchical clustering on the basis of 186 m 6A-related lncRNAs, followed by univariate Cox regression analysis, which were associated with the OS and tumor microenvironment. For that, the TCGA-LIHC database was randomly separated into the training dataset and the testing dataset. Four m6A-related lncRNAs prognostic risk signatures were constructed and validated in TCGA-LIHC. In addition, we performed co-expression analysis to identify the potential function of the four m6A-related lncRNAs. Finally, Q-PCR was implemented to detect three of the four m6A-related lncRNAs expressions in cells and tissue.

Abundant evidence had suggested that m6A modification was involved in cancer pathogenesis and progression with epigenetic regulation of oncogenes or tumor suppressor genes (He et al., 2018b). However, how it works in a lncRNA-dependent manner during hepatocellular carcinoma is still unknown. m6A modification can regulate malignance and invasion of several tumors by modifying specific lncRNAs. Yuan et al. reported that overexpressing METTL3 significantly promoted osteogenic differentiation of primary ligament fibroblasts by increasing m6A methylation of long non-coding RNA (lncRNA) X-inactive specific transcript (XIST) (Yuan et al., 2021). M6A reader YTHDF3 negatively regulates the long non-coding RNA GAS5, inhibiting the progression of colorectal cancer by interacting with and triggering YAP phosphorylation and degradation (Ni et al., 2019). Furthermore, Yang et al. indicated that low expression of TRAF3IP2-AS1 promotes the progression of NONO-TFE3 translocation renal cell carcinoma by stimulating the N6-methyladenosine of PARP1 mRNA and downregulating PTEN (Yang et al., 2021). Meng et al. demonstrated that m6A-mediated upregulation of LINC00857 promotes the occurrence of pancreatic cancer by functioning as a competing endogenous RNA (ceRNA) for sponging miR-150-5p, leading to overexpressing E2F3 and ultimately promoting tumorigenesis in PC (Meng et al., 2021) Studies had shown that the m6A modification of lncRNA may affect cancer occurrence and development, and lncRNA could act as a competitive endogenous RNA to target m6A modulators, thereby affecting tumor invasion and progression. Taken together, we firmly believe that m6A modification could target lncRNA to regulate the occurrence and development of tumors, and we should pay more attention to the function and interaction between m6A modification and lncRNA to identify potential prognostic markers or therapeutic targets for cancers.

In this study, four m6A-related lncRNAs were identified to be related to HCC prognosis, and a few studies had revealed biological functions of these lncRNAs. Elevation of lncRNA BACE1-AS was a potential mechanism to inhibit the proliferation and invasion of human ovarian cancer stem cells and may be a novel target of anisomycin in the treatment of ovarian cancer (Chen et al., 2016). LncRNA NRAV negatively regulated the initial transcription of multiple key interferon-stimulated genes and was important in regulating the antiviral interferon response (Ouyang et al., 2014). In the meantime, lncRNA NRAV has also been proven to be an independent prognostic factor for patients with low-grade glioma (Maimaiti et al., 2021). A recent meta-analysis reported that KDM4A-AS1 was screened as an important prognostic factor to predict overall survival in HCC patients. However, there were few reports on the modification and interaction between lncRNAs and m6A-related genes in HCC. Thus, we hope that our results would help identify m6a-related prognostic lncRNAs, thereby providing novel insights into their potential role and serving as novel therapeutic targets in HCC.

However, there were several limitations and challenges in our study. The modification and interaction between lncRNA and m6A-related genes should be confirmed by in vivo and in vitro experiments. Indeed, it is essential to further investigate the molecular mechanisms and signaling pathways of m6A-related lncRNAs involved in HCC.

In conclusion, our study revealed two distinct subtypes based on m6A-related lncRNAs profiling and identified a novel m6A-related lncRNAs signature to predict prognosis and guide individual treatment strategies in HCC patients.

Data Availability Statement

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

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committee of Affiliated Hangzhou First People’s Hospital. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

LY and LZ conceived and designed the study. LY, LZ, SG, and YF conducted the study HZ contributed to the acquisition of data. LY and LZ analyzed the data. LY, LZ, and SG interpreted the data. LY, LZ, JX, and RX reviewed and edited the manuscript. All authors read and approved the manuscript.

Funding

This work was supported by the Natural Science Foundation of Zhejiang Province (LY21H160016), the Medical Science and Technology Plan of Zhejiang Province (2018KY569), and the Medical Science and Technology Project of Hangzhou (A20210073).

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/fmolb.2022.807418/full#supplementary-material

Abbreviations

AUC, area under the curve; HCC, hepatocellular carcinoma; HR, hazard ratio; m6A-related lnRNAs, N6-Methylandenosine-Related lncRNAs; OS, overall survival; ROC curve, receiver operating characteristic curve; ssGSEA, single sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; TNM: tumor-node-metastasis.

References

Abbastabar, M., Sarfi, M., Golestani, A., and Khalili, E. (2018). lncRNA Involvement in Hepatocellular Carcinoma Metastasis and Prognosis. Excli j 17, 900–913. doi:10.17179/excli2018-1541

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA Interference Reveals that Oncogenic KRAS-Driven Cancers Require TBK1. Nature 462 (7269), 108–112. doi:10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

Batista, P. J. (2017). The RNA Modification N 6 -methyladenosine and its Implications in Human Disease. Genomics, Proteomics & Bioinformatics 15 (3), 154–163. doi:10.1016/j.gpb.2017.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Berretta, J., and Morillon, A. (2009). Pervasive Transcription Constitutes a New Level of Eukaryotic Genome Regulation. EMBO Rep. 10 (9), 973–982. doi:10.1038/embor.2009.181

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Q., Liu, X., Xu, L., Wang, Y., Wang, S., Li, Q., et al. (2016). Long Non-coding RNA BACE1-AS Is a Novel Target for Anisomycin-Mediated Suppression of Ovarian Cancer Stem Cell Proliferation and Invasion. Oncol. Rep. 35 (4), 1916–1924. doi:10.3892/or.2016.4571

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, D., Wang, H., Zhu, L., Jin, H., and Wang, X. (2018). N6-methyladenosine Links RNA Metabolism to Cancer Progression. Cell Death Dis 9 (2), 124. doi:10.1038/s41419-017-0129-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, X., Su, R., Weng, H., Huang, H., Li, Z., and Chen, J. (2018). RNA N6-Methyladenosine Modification in Cancers: Current Status and Perspectives. Cell Res 28 (5), 507–517. doi:10.1038/s41422-018-0034-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferlay, J., Soerjomataram, I., Dikshit, R., Eser, S., Mathers, C., Rebelo, M., et al. (2015). Cancer Incidence and Mortality Worldwide: Sources, Methods and Major Patterns in GLOBOCAN 2012. Int. J. Cancer 136 (5), E359–E386. doi:10.1002/ijc.29210

CrossRef Full Text | Google Scholar

Fitzmaurice, C., Fitzmaurice, C., Allen, C., Barber, R. M., Barregard, L., Bhutta, Z. A., et al. (2017). Global, Regional, and National Cancer Incidence, Mortality, Years of Life Lost, Years Lived with Disability, and Disability-Adjusted Life-Years for 32 Cancer Groups, 1990 to 2015: A Systematic Analysis for the Global Burden of Disease Study. JAMA Oncol. 3 (4), 524–548. doi:10.1001/jamaoncol.2016.5688

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

He, L., Li, H., Wu, A., Peng, Y., Shu, G., and Yin, G. (2019). Functions of N6-Methyladenosine and its Role in Cancer. Mol. Cancer 18 (1), 176. doi:10.1186/s12943-019-1109-9

PubMed Abstract | CrossRef Full Text | Google Scholar

He, L., Li, J., Wang, X., Ying, Y., Xie, H., Yan, H., et al. (2018). The Dual Role of N6‐methyladenosine Modification of RNAs Is Involved in Human Cancers. J. Cel Mol Med 22 (10), 4630–4639. doi:10.1111/jcmm.13804

CrossRef Full Text | Google Scholar

He, Y., Hu, H., Wang, Y., Yuan, H., Lu, Z., Wu, P., et al. (2018). ALKBH5 Inhibits Pancreatic Cancer Motility by Decreasing Long Non-coding RNA KCNK15-AS1 Methylation. Cell Physiol Biochem 48 (2), 838–846. doi:10.1159/000491915

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoshida, Y., Moeini, A., Alsinet, C., Kojima, K., and Villanueva, A. (2012). Gene Signatures in the Management of Hepatocellular Carcinoma. Semin. Oncol. 39 (4), 473–485. doi:10.1053/j.seminoncol.2012.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Karthiya, R., and Khandelia, P. (2020). m6A RNA Methylation: Ramifications for Gene Expression and Human Health. Mol. Biotechnol. 62 (10), 467–484. doi:10.1007/s12033-020-00269-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Lai, M.-c., Yang, Z., Zhou, L., Zhu, Q.-q., Xie, H.-y., Zhang, F., et al. (2012). Long Non-coding RNA MALAT-1 Overexpression Predicts Tumor Recurrence of Hepatocellular Carcinoma after Liver Transplantation. Med. Oncol. 29 (3), 1810–1816. doi:10.1007/s12032-011-0004-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Lan, Y., Liu, B., and Guo, H. (2021). The Role of M6A Modification in the Regulation of Tumor-Related lncRNAs. Mol. Ther. - Nucleic Acids 24, 768–779. doi:10.1016/j.omtn.2021.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C., Miao, R., Liu, S., Wan, Y., Zhang, S., Deng, Y., et al. (2017). Down-regulation of miR-146b-5p by Long Noncoding RNA MALAT1 in Hepatocellular Carcinoma Promotes Cancer Growth and Metastasis. Oncotarget 8 (17), 28683–28695. doi:10.18632/oncotarget.15640

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, N., Zhou, K. I., Parisien, M., Dai, Q., Diatchenko, L., and Pan, T. (2017). N 6-methyladenosine Alters RNA Structure to Regulate Binding of a Low-Complexity Protein. Nucleic Acids Res. 45 (10), 6051–6063. doi:10.1093/nar/gkx141

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Q., and Gregory, R. I. (2019). RNAmod: an Integrated System for the Annotation of mRNA Modifications. Nucleic Acids Res. 47 (W1), W548–w555. doi:10.1093/nar/gkz479

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J. z., Yang, F., Zhou, C. c., Liu, F., Yuan, J. h., Wang, F., et al. (2017). METTL14 Suppresses the Metastatic Potential of Hepatocellular Carcinoma by Modulating N 6 ‐methyladenosine‐dependent Primary MicroRNA Processing. Hepatology 65 (2), 529–543. doi:10.1002/hep.28885

PubMed Abstract | CrossRef Full Text | Google Scholar

Maimaiti, A., Jiang, L., Wang, X., Shi, X., Pei, Y., Hao, Y., et al. (2021). Identification and Validation of an Individualized Prognostic Signature of Lower-Grade Glioma Based on Nine Immune Related Long Non-coding RNA. Clin. Neurol. Neurosurg. 201, 106464. doi:10.1016/j.clineuro.2020.106464

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, X., Deng, Y., He, S., Niu, L., and Zhu, H. (2021). m6A-Mediated Upregulation of LINC00857 Promotes Pancreatic Cancer Tumorigenesis by Regulating the miR-150-5p/E2F3 AxisA-Mediated Upregulation of LINC00857 Promotes Pancreatic Cancer Tumorigenesis by Regulating the miR-150-5p/E2F3 Axis. Front. Oncol. 11, 629947. doi:10.3389/fonc.2021.629947

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 (1), 143. doi:10.1186/s12943-019-1079-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Niu, Y., Lin, Z., Wan, A., Chen, H., Liang, H., Sun, L., et al. (2019). RNA N6-Methyladenosine Demethylase FTO Promotes Breast Tumor Progression through Inhibiting BNIP3. Mol. Cancer 18 (1), 46. doi:10.1186/s12943-019-1004-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouyang, J., Zhu, X., Chen, Y., Wei, H., Chen, Q., Chi, X., et al. (2014). NRAV, a Long Noncoding RNA, Modulates Antiviral Responses through Suppression of Interferon-Stimulated Gene Transcription. Cell Host & Microbe 16 (5), 616–626. doi:10.1016/j.chom.2014.10.001

CrossRef Full Text | Google Scholar

Panneerdoss, S., Eedunuri, V. K., Yadav, P., Timilsina, S., Rajamanickam, S., Viswanadhapalli, S., et al. (2018). Cross-talk Among Writers, Readers, and Erasers of m6A Regulates Cancer Growth and Progression. Sci. Adv. 4 (10), eaar8263. doi:10.1126/sciadv.aar8263

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 Cancer J. Clin. 65 (2), 87–108. doi:10.3322/caac.21262

CrossRef Full Text | Google Scholar

Venook, A. P., Papandreou, C., Furuse, J., and de Guevara, L. L. (2010). The Incidence and Epidemiology of Hepatocellular Carcinoma: a Global and Regional Perspective. Oncologist 15 Suppl 4 (Suppl. 4), 5–13. doi:10.1634/theoncologist.2010-S4-05

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F., Ying, H.-Q., He, B.-S., Pan, Y.-Q., Deng, Q.-W., Sun, H.-L., et al. (2015). Upregulated lncRNA-UCA1 Contributes to Progression of Hepatocellular Carcinoma through Inhibition of miR-216b and Activation of FGFR1/ERK Signaling Pathway. Oncotarget 6 (10), 7899–7917. doi:10.18632/oncotarget.3219

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, H., Ma, H., and Zhou, D. (2013). Plasma HULC as a Promising Novel Biomarker for the Detection of Hepatocellular Carcinoma. Biomed. Res. Int. 2013, 136106. doi:10.1155/2013/136106

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J. D., Hainaut, P., Gores, G. J., Amadou, A., Plymoth, A., and Roberts, L. R. (2019). A Global View of Hepatocellular Carcinoma: Trends, Risk, Prevention and Management. Nat. Rev. Gastroenterol. Hepatol. 16 (10), 589–604. doi:10.1038/s41575-019-0186-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J. D., and Heimbach, J. K. (2020). New Advances in the Diagnosis and Management of Hepatocellular Carcinoma. Bmj 371, m3544. doi:10.1136/bmj.m3544

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Chen, Y., Liu, N., Shi, Q., Han, X., Gan, W., et al. (2021). Low Expression of TRAF3IP2-AS1 Promotes Progression of NONO-TFE3 Translocation Renal Cell Carcinoma by Stimulating N6-Methyladenosine of PARP1 mRNA and Downregulating PTEN. J. Hematol. Oncol. 14 (1), 46. doi:10.1186/s13045-021-01059-5

CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, X., Shi, L., Guo, Y., Sun, J., Miao, J., Shi, J., et al. (2021). METTL3 Regulates Ossification of the Posterior Longitudinal Ligament via the lncRNA XIST/miR-302a-3p/USP8 Axis. Front. Cel Dev. Biol. 9, 629895. doi:10.3389/fcell.2021.629895

CrossRef Full Text | Google Scholar

Zhong, L., Liao, D., Zhang, M., Zeng, C., Li, X., Zhang, R., et al. (2019). YTHDF2 Suppresses Cell Proliferation and Growth via Destabilizing the EGFR mRNA in Hepatocellular Carcinoma. Cancer Lett. 442, 252–261. doi:10.1016/j.canlet.2018.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat. Commun. 10 (1), 1523. doi:10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar