Home
JournalsCollections
For Authors For Reviewers For Editorial Board Members
Article Processing Charges Open Access
Ethics Advertising Policy
Editorial Policy Resource Center
Company Information Contact Us
OPEN ACCESS

Prognostic Role and Potential Mechanisms of N6-methyladenosine-related Long Noncoding RNAs in Hepatocellular Carcinoma

  • Tianxing Dai1,#,
  • Jing Li2,#,
  • Linsen Ye1,
  • Haoyuan Yu1,
  • Mingbin Deng1,
  • Wei Liu3,
  • Hua Li1,
  • Yang Yang1,*  and
  • Guoying Wang1,* 
Journal of Clinical and Translational Hepatology   2022;10(2):308-320

doi: 10.14218/JCTH.2021.00096

Received:

Revised:

Accepted:

Published online:

 Author information

Citation: Dai T, Li J, Ye L, Yu H, Deng M, Liu W, et al. Prognostic Role and Potential Mechanisms of N6-methyladenosine-related Long Noncoding RNAs in Hepatocellular Carcinoma. J Clin Transl Hepatol. 2022;10(2):308-320. doi: 10.14218/JCTH.2021.00096.

Abstract

Background and Aims

Numerous studies have explored the important role of N6-methyladenosine (m6A) in cancer. Nonetheless, the interaction between m6A and long noncoding RNAs (lncRNAs) is poorly investigated. Herein, we systematically analyzed the role and prognostic value of m6A-related lncRNAs in hepatocellular carcinoma (HCC).

Methods

The m6A-related lncRNAs were identified based on the correlation coefficients with m6A-related genes in HCC from The Cancer Genome Atlas. Subsequently, a novel risk score model was determined using the least absolute shrinkage and selection operator Cox regression analyses. Univariate and multivariate Cox analyses were used to identify independent prognostic factors for overall survival (OS) of HCC; thereafter, a prognostic nomogram was constructed.

Results

A total of 259 lncRNAs showed significant correlations with m6A in HCC, while 29 lncRNAs had prognostic significance. Further, six critical m6A-related lncRNAs (NRAV, SNHG3, KDM4A-AS1, AC074117.1, AC025176.1, and AL031985.3) were screened out to construct a novel risk score model which classified HCC patients into high- and low-risk groups. Survival analyses revealed that patients in the high-risk group exhibited worse OS, both in the training and validation groups. The risk score was also identified as an independent prognostic factor of OS, and a nomogram was established and verified with superior prediction capacity. Besides, the risk score significantly correlated with the expression of immune checkpoint genes and immune subtypes.

Conclusions

These findings indicated the significant role of m6A-related lncRNAs in HCC and the potential application of the novel risk score model for prognostic prediction.

Keywords

Hepatocellular carcinoma, N6-methyladenosine, Long noncoding RNAs, Prognosis, Immune checkpoints

Introduction

Liver cancer is one of the most prevalent malignancies with increasing incidence worldwide.1,2 Hepatocellular carcinoma (HCC) accounts for the majority of liver cancer, with an overall 5-year survival rate of less than 20%. Due to the insidious feature of HCC, most patients are diagnosed at advanced stages, and thus lose the opportunity for radical treatments.3 Systemic treatments with targeted therapy dramatically alter the situation of treatment for extrahepatic metastases and vascular invasion.4 Nevertheless, the improvement of survival is limited. Recent reports have suggested that immunotherapy based on immune checkpoint inhibitors is a promising therapeutic strategy for patients with advanced HCC, while its efficiency remains less than 20%.5 Therefore, advances in the treatment of HCC are urgently needed, which requires comprehensive investigation on the pathogenesis of HCC.

N6-methyladenosine (m6A), the most abundant mRNA modification in eukaryotic cells, has been found to play important roles in various pathophysiological processes, including the occurrence, development, and metastasis of cancers, including HCC.6–8 Moreover, it is a dynamic reversible modification implicated in post-transcriptional gene regulation, referring to the RNA splicing, transporting, translation, stability, and the high-level structure.9 Notably, m6A modification is highly conserved and can be catalyzed or removed by specific enzymes, including the m6A methyltransferases (“writers”) and the demethylases (“erasers”). It is recognized by m6A-binding proteins (“readers”).10 Increasing evidence reveals that m6A is a novel regulatory mechanism regulating gene expression and cell fate decisions. The potential mechanisms are briefly summarized in Figure S1.

With the advent of next-generation sequencing, the landscape of m6A in the transcriptome has been determined, and the characteristics and patterns of m6A modification on mRNAs have also been delineated.11 Interestingly, researchers have found that the m6A modification is not only ubiquitous in mRNA but also prevalent in noncoding RNAs (ncRNAs), including microRNAs (miRNAs), long noncoding RNAs (lncRNAs), and circular RNAs (circRNAs).11–13 The m6A modification of ncRNAs potentially affects their functions. As such, the interaction between m6A modification and ncRNAs provides a novel approach to explore the pathophysiological mechanisms implicated in diseases, including cancer.14

LncRNAs are a large class of ncRNAs with a length of above 200 nucleotides. Generally, they lack protein-coding capacity and primarily regulate gene expression.15 With the identification of their multiple functions, lncRNAs have been supposed as key players in various physiological processes and diseases by regulating gene expression at both transcriptional and post-transcriptional levels.16

Recent studies have revealed that lncRNAs can also be extensively m6A-modified, and the potential regulatory mechanisms for lncRNA functions have also been identified.13,17–22 First, the occurrence of m6A modifications might provide binding sites for the m6A “readers” or trigger RNA-binding protein entry by altering the structure of local RNA. Second, m6A modifications potentially modulate the RNA-DNA triple helix structure and influence the binding between lncRNAs and corresponding DNA sites. Nonetheless, the crosstalk and mechanisms between m6A modification and lncRNAs remain to be further elaborated.

In the present study, we comprehensively investigated the role of m6A-related lncRNAs in HCC. First, the m6A-related lncRNAs were identified in HCC based on the differentially expressed m6A-related genes between tumor and normal tissues in The Cancer Genome Atlas (TCGA) database. Subsequently, the m6A-related lncRNAs with prognostic significance were screened out and used for consensus clustering, survival analysis, and pathway enrichment evaluation. A novel risk score model based on the m6A-related lncRNAs was constructed by the least absolute shrinkage and selection operator (LASSO) Cox regression. Additionally, multivariate Cox analysis was used to assess the independent prognostic significance in the training and validation groups. A prognostic nomogram was established for the overall survival (OS) of HCC patients. Finally, we evaluated the correlation between the risk model and the immune checkpoint expression, immune subtype, tumor stemness, and drug susceptibility. Our findings provide novel perspectives on the interaction between m6A and lncRNAs in HCC, and a novel prognostic indicator for HCC patients.

Methods

Datasets and patients

The level 3 RNA-seq transcriptome profiles of 374 HCC samples and 50 normal samples were obtained on March 1, 2021, from the TCGA-LIHC database (https://portal.gdc.cancer.gov/ ). The corresponding clinical and pathological characteristics, including age, sex, tumor grade, TNM stage, follow-up time, and survival status, were also collected. Duplicated sample information from the same patient was combined as averages. Finally, the data from 370 HCC patients were included. The transcriptome data of TCGA-LIHC were further distinguished with mRNA and lncRNAs, and 19,604 mRNA and 14,086 lncRNAs in HCC were obtained.

TCGA is a public database, and all cases involved in the database have been consented to use for TCGA and obtained ethical approval, which is available online to be downloaded and analyzed by individual researchers. This work was based on the open-source data and associated identification information of all enrolled cases had been blinded; hence, the study protocol was exempted from ethical approval. The study was also performed and reported in accordance with the Helsinki Declaration. The data analysis procedure is shown in Figure 1.

Flowchart of the data analysis procedures.
Fig. 1  Flowchart of the data analysis procedures.

Identification of the m6A-related lncRNAs with prognostic significance in HCC

First, the m6A-related genes were summarized by reviewing the previous literature, from which 36 genes were obtained (summarized in Table S1). Then, the comparison analyses between tumors and normal controls were conducted using the Wilcoxon test with R package “limma” to identify the differentially expressed m6A-related genes (m6A-DEGs) in HCC. Subsequently, the lncRNAs associated with m6A-DEGs were screened through Pearson’s correlation analyses (coefficient >0.5, and p<0.001), through which 259 lncRNAs were obtained. Finally, univariate Cox regression analyses using the R package “survival” were conducted to identify m6A-related lncRNAs with significant prognostic values for the OS of HCC patients. A p-value <0.001 was considered a prognostic candidate for m6A-related lncRNA.

Consensus clustering analysis

Based on the consensus expression of prognostic m6A-related lncRNAs in tumors, the 370 HCC cases were clustered into two subgroups using the R package “ConsensusClusterPlus”. A survival analysis based on the Kaplan-Meier method was conducted using the R package “survival” to compare the differential OS between clusters. The clinical characteristics were also compared by the Chi-square test. Gene set enrichment analysis (GSEA) was conducted to evaluate the enriched pathways within different subgroups. Finally, the expression of several immune checkpoint-related genes was evaluated between the two subgroups.

Construction of the risk score model based on the prognostic m6A-related lncRNAs in HCC

Based on the results of univariate Cox analyses, the LASSO Cox regression was utilized to eliminate the highly correlated lncRNAs using the R package “glment”. Ultimately, only six m6A-related lncRNAs were identified and entered to construct the novel risk score model. With the coefficients and expression levels of lncRNAs, the risk score for each HCC patient was calculated by the following formula:

Risk score=i=1nCoefixi(Coefi=coefficient, xi=lncRNA expression).

Then, all patients were randomly divided into two cohorts (186 cases in the training group, and 184 cases in the validation group). In each cohort, the HCC patients were further divided into low- and high-risk groups by the median of the risk score. Thereafter, the survival analyses using the R package “survival” for OS were conducted. Receiver operating characteristic (ROC) curves were also adopted to evaluate the prognostic accuracy. Besides, the survival analyses were further evaluated in different subgroups divided by the clinical factors.

To evaluate the independent prognostic significance of the risk score model in HCC, we further conducted both univariate and multivariate Cox regression analyses. Then, a nomogram for OS prediction was established. The prediction performance was assessed by the concordance index (C index) and calibration curves.

Correlations between risk score model and immune-related genes and immune subtypes in HCC

To explore the potential mechanisms of m6A-related lncRNAs in HCC, we further analyzed the relationship between the scoring model and tumor immunity. First, we obtained the immune-related genes with negative regulation functions from the Tracking Tumor Immunophenotype (TIP) database (http://biocc.hrbmu.edu.cn/TIP/index.jsp ).23 The difference of these negative immunoregulatory genes between the high- and low-risk group were compared using the Wilcoxon test.

Besides, the profile of immune subtypes for HCC in the TGGA cohort was obtained from the UCSC-Xena (https://xenabrowser.net/datapages/ , version: 2018-04-03). The Wilcoxon test was used to compare the risk score between different immune subtypes.

Statistical analysis

All statistical analyses were conducted using R (https://www.r-project.org/ , version 4.0.3). A p-value <0.05 (two-sided) was considered statistically significant.

Results

Identification of the differentially expressed m6A-related genes in HCC

By comparing the expression of 36 m6A-related genes in 374 HCC tumors and 50 normal tissues, 24 genes with significant differential expression were obtained (Fig. 2A, details in Table S2). Out of these, 22 genes were upregulated (METTL3, METTL16, WTAP, VIRMA, RBM15B, PCIF1, YTHDF1, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, RBMX, NXF1, TRMT112, ZCCHC4, CPSF6, CBLL1, SETD2, SRSF3, SRSF10, and PRRC2A; all p<0.05), while two genes (IGF2BP1 and IGF2BP3, all p<0.05) were down-regulated in the tumor. The expression profiles, correlations, and potential regulatory function and pathways of these m6A-related genes in HCC are summarized in Figure S2.

Identification of m<sup>6</sup>A-related lncRNAs with prognostic significance in HCC.
Fig. 2  Identification of m6A-related lncRNAs with prognostic significance in HCC.

(A) Bar plot for exhibition of the m6A-DEGs expression in HCC tumors and normal controls. (B) The correlation plot of the m6A-related genes (n=24) with the lncRNAs (n=259). (C) Univariate analysis for m6A-related lncRNAs with prognostic significance in HCC (n=29). (D) Heatmap of the expression of m6A-related lncRNAs in HCC tumors and normal controls. (E) The correlation among these m6A-related lncRNAs in HCC. *p<0.05. HCC, hepatocellular carcinoma; lncRNAs, long noncoding RNAs; m6A, N6-methyladenosine; m6A-DEGs, differentially expressed m6A-related genes.

Identification of the m6A-related lncRNAs with prognostic significance in HCC

With annotation of the TCGA-LIHC transcriptome profile, 14,086 lncRNAs expressed in HCC were identified. After the correlation analyses of m6A-DEGs and lncRNAs, 259 m6A-related lncRNAs were screened out (Fig. 2B, details in Table S3). Further univariate Cox analyses revealed that 29 m6A-related lncRNAs were related to the OS of HCC (all p<0.001), and all significant m6A-related lncRNAs were the risk factors for the prognosis (all hazard ratio >1) (Fig. 2C). The expression profiles of each m6A-related lncRNAs in the HCC cohort are shown in the heatmap (Fig. 2D), and all lncRNAs appeared up-regulated in HCC tumors, which indicated the important roles of these m6A-related lncRNAs in HCC. Furthermore, correlation analyses revealed that nearly all m6A-related lncRNAs in HCC were significantly correlated (Fig. 2E).

Consensus clustering of HCC patients based on the prognostic m6A-related lncRNAs

To evaluate the significance of m6A-related lncRNAs in the development of HCC, consensus clustering analysis was performed to separate the HCC tumors into various clusters based on the expression profile characteristics of these prognostic m6A-related lncRNAs. The cumulative distribution function (CDF) with different clustering strategies from k=2–9 and the relative changes of the area under CDF curves are illustrated in Figure 3A–B. The sample distribution determined by the clustering strategies from k=2–9 is shown in Figure 3C. Considering the CDF increase and consistent expression of m6A-related lncRNAs in HCC, k=2 was adopted as the clustering criterion (Fig. 3D). There were a total of 310 and 60 cases in clusters 1 and 2, respectively.

Consensus clustering of HCC based on the m<sup>6</sup>A-related lncRNAs.
Fig. 3  Consensus clustering of HCC based on the m6A-related lncRNAs.

(A–C) The CDF, relative changes in area under the CDF curves, and tracking plots show with the index from 2 to 9. (D) The distribution of different clusters with the clustering index was 2. (E) Survival curves in different clusters. (F) The heatmap with visualization of the expression of 29 m6A-related lncRNAs in HCC patients and the correlation of the risk score with other clinical factors. (G) Pathway analyses by GSEA in different clusters. (H) The expression profiles of immune checkpoint genes (PD-1, PD-L1, CTLA-4, LAG3, TIM3, and TIGIT in the different clusters. ns, p>0.05, *p<0.05, **p<0.01, ***p<0.001. CDF, cumulative distribution function; CTLA-4, cytotoxic T lymphocyte-associated antigen-4; GSEA, gene set enrichment analysis; HCC, hepatocellular carcinoma; LAG3, lymphocyte-activation gene 3; lncRNAs, long noncoding RNAs; m6A, N6-methyladenosine; PD-1, programmed cell death 1; PD-L1, PD-1 ligand 1; TIGIT, T cell immunoreceptor with immunoglobulin and ITIM domain; TIM3, T-cell immunoglobulin and mucin domain 3.

Prognostic and functional differences between HCC clusters

Survival analysis revealed that HCC patients in cluster 2 had worse OS than patients in cluster 1 (p<0.001) (Fig. 3E). Furthermore, the Chi-square tests confirmed the significant difference of several clinical factors between the two clusters (Fig. 3F). The proportions of female patients (p<0.05) and advanced pathological grade (grade 3–4) (p<0.001), T stage (T3–4) (p<0.01), and clinical stage (stage III–IV) (p<0.001) were higher in cluster 2 than in cluster 1. Besides, the expression of all enrolled m6A-related lncRNAs was upregulated in cluster 2. These results suggested that the clustering strategy based on the m6A-related lncRNAs had important clinical significance.

Additionally, we evaluated the potential functional differences between the two clusters divided by the m6A-related lncRNAs. GSEA indicated that the pathways of the spliceosome, cell cycle, RNA degradation, nucleotide excision repair, and ubiquitin-mediated proteolysis were enriched in cluster 2, while the pathways of PPAR signaling, retinol metabolism, complement and coagulation cascades, drug metabolism cytochrome P450, and fatty acid metabolism were enriched in cluster 1 (Fig. 3G). The differentially enriched pathways in different clusters potentially account for the distinct expression patterns of m6A-related lncRNAs.

Based on the role of m6A modification in the regulation of the tumor microenvironment (TME), we attempted to investigate whether there were differences in the immune checkpoint genes’ expression between the two clusters. Results showed that the expression of programmed cell death 1 (PD-1), cytotoxic T lymphocyte-associated antigen-4 (CTLA-4), lymphocyte-activation gene 3 (LAG3), T-cell immunoglobulin and mucin domain 3 (TIM3), and T cell immunoreceptor with immunoglobulin and ITIM domain (TIGIT) were all higher in cluster 2 (all p<0.001), while no significance was only found with the PD-1 ligand 1 (PD-L1) (p>0.05) (Fig. 3H). This indicated that differentially expressed m6A-related lncRNAs might play a role in the regulation of the expression of immune checkpoint genes, which could help to guide the treatment of immune checkpoint inhibitors for HCC patients and provide potential for efficacy prediction.

Development of the novel prognostic risk score model based on the m6A-related lncRNAs

Based on the 29 prognostic m6A-related lncRNAs from univariate Cox analyses, we applied LASSO Cox regression analysis for further screening of independent prognostic factors for OS of HCC (Fig. 4A–B). Six critical lncRNAs (KDM4A-AS1, NRAV, AC074117.1, SNHG3, AC025176.1, and AL031985.3) were finally identified. Then, the risk score based on the six m6A-related lncRNAs for each HCC patient both in the training (n=186) and validation (n=184) group was calculated according to their expression levels and coefficients determined by LASSO Cox regression. In each group, HCC patients were further divided into high- and low-risk subgroups by the median of risk scores. Survival analyses revealed that patients with high risk had worse OS than those with low risk in the training group (p<0.001; Fig. 4C), and the prognostic difference was also noted in the validation group (p=0.012; Fig. 4D). The time-dependent ROCs were further plotted. In the training group, the area under the curve (AUC) for 1-, 3-, and 5-year OS was 0.794, 0.771, and 0.694, respectively (Fig. 4E). In the validation group, the AUC for 1-, 3-, and 5-year OS was 0.661, 0.630, and 0.692, respectively (Fig. 4F). Besides, 42% of high-risk patients died during the follow-up period, whereas only 24% of low-risk patients died in the training group (Fig. 5A). Similar results were found in the validation group. Further, 43% of patients died in the high-risk subgroup, while only 31% died in the low-risk subgroup (Fig. 5B). The risk plots of both the training and validation groups showed the risk score distribution, survival status, and expression of the six lncRNAs of each HCC patient (Fig. 5C–D). These findings indicated that the risk score model based on the m6A-related lncRNAs could effectively differentiate and accurately predict the OS of HCC patients.

Identification of critical m<sup>6</sup>A-related lncRNAs for the novel risk score model and prognostic evaluation in HCC.
Fig. 4  Identification of critical m6A-related lncRNAs for the novel risk score model and prognostic evaluation in HCC.

(A–B) Screening of the critical m6A-related lncRNAs by LASSO Cox regression. (C–D) Survival analyses of high- and low-risk groups in the training and validation cohorts. (E–F) Time-dependent ROC curves based on the risk score for 1-, 3-, and 5-year OS of HCC in the training and validation groups. AUC, area under the curve; HCC, hepatocellular carcinoma; LASSO, least absolute shrinkage and selection operator; lncRNAs, long noncoding RNAs; m6A, N6-methyladenosine; OS, overall survival; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Distribution of the risk score, survival status, and expression of the six m<sup>6</sup>A-related lncRNAs in the training and validation groups of HCC.
Fig. 5  Distribution of the risk score, survival status, and expression of the six m6A-related lncRNAs in the training and validation groups of HCC.

(A–B) Proportions of death in the high- and low-risk groups of the training and validation cohorts. (C–D) Risk maps of the risk score, survival status, and expression of lncRNAs in the training and validation cohorts. HCC, hepatocellular carcinoma; lncRNAs, long noncoding RNAs; m6A, N6-methyladenosine; TCGA, The Cancer Genome Atlas.

Independent prognostic significance of the novel risk score model based on the m6A-related lncRNAs

Both univariate and multivariate Cox analyses were performed to identify the independent prognostic factors for the OS of HCC patients. In the training group, only the clinical stage and the novel risk score were significant in the univariate [stage, p<0.001, HR=1.772 (1.275–2.462); risk score, p<0.001, HR=1.181 (1.113–1.254)] (Fig. 6A) and multivariate Cox analyses [stage, p=0.009, HR=1.573 (1.118–2.212); risk score, p<0.001, HR=1.166 (1.092–1.246)] (Fig. 6B). Similar results were also observed in the validation group, and only the clinical stage and risk score were found with significance both in the univariate [stage, p<0.001, HR=1.624 (1.244–2.119); risk score, p<0.001, HR=1.218 (1.091–1.361)] (Fig. 6C) and multivariate Cox analyses [stage, p<0.001, HR=1.585 (1.205–2.086); risk score, p=0.004, HR=1.190 (1.056–1.341)] (Fig. 6D). Furthermore, we conducted various subgroup analyses within different clinical factors to evaluate the prognostic significance of our risk score model in all HCC patients. We found that patients in the high-risk group showed worse OS both with age ≤60 years (p<0.001; Fig. 6E) and >60 years (p<0.001; Fig. 6F), female sex (p=0.003; Fig. 6G) and male sex (p<0.001; Fig. 6H), grade 1–2 (p<0.001; Fig. 6I) and 3–4 (p=0.010; Fig. 6J), and stage I–II (p<0.001; Fig. 6K) and III–IV (p=0.009; Fig. 6L). The findings showed that when patients were divided based on the risk score, a significant difference was noted in OS between the high- and low-risk patients regardless of the clinical factor subgroups.

Evaluation of the independent prognostic significance of the novel risk score model based on m<sup>6</sup>A-related lncRNAs in HCC.
Fig. 6  Evaluation of the independent prognostic significance of the novel risk score model based on m6A-related lncRNAs in HCC.

(A–B) Univariate and multivariate Cox analyses with risk score and other clinical factors for OS of HCC in the training group. (C–D) Univariate and multivariate Cox analyses with risk score and other clinical factors for OS of HCC in the validation group. (E–L) Survival analyses of the risk score in different subgroups of various clinical factors: age (≤/>60 years), sex (female/male), pathological grade (1–2/3–4), and clinical stage (I–II/III–IV). HCC, hepatocellular carcinoma; HR, hazard ratio; lncRNAs, long noncoding RNAs; m6A, N6-methyladenosine;

Besides, the distribution of different clinical factors was compared between high- and low-risk groups of all HCC patients. Results showed significant differences with pathological grade, clinical stage, and cluster between these two risk groups. There were more patients with advanced pathological grade (3–4) (p<0.001) and stage (III–IV) (p<0.05) in the high-risk group, and all patients in cluster 2 were classified into the high-risk group (p<0.001) (Fig. 7A). Higher average risk scores in grade 3–4, stage III–IV, and cluster 2 (all p<0.001) were also noted (Fig. 7B).

Clinical correlation of the risk score model and the prognostic nomogram for OS of HCC established based on the risk score model.
Fig. 7  Clinical correlation of the risk score model and the prognostic nomogram for OS of HCC established based on the risk score model.

(A) Heatmap with visualization of the correlations between risk score and other clinical characteristics and the expression of the six m6A-related lncRNAs in HCC. (B) Distribution of risk score in subgroups of different pathological grades, clinical stages, and clusters. (C) Prognostic nomogram established based on the risk score and clinical stage for prediction of 1-, 3-, and 5-year OS of HCC. (D) Calibration curves for evaluation the prognostic accuracy of the nomogram. *p<0.05, **p<0.01, ***p<0.001. HCC, hepatocellular carcinoma; lncRNAs, long noncoding RNAs; m6A, N6-methyladenosine; OS, overall survival.

A novel prognostic nomogram was constructed based on the clinical stage and risk score (Fig. 7C), which could provide an easy-to-use method for predicting the 1-, 3-, and 5-year OS of HCC. The nomogram exhibited good efficacy in estimating the OS with a high C-index of 0.674. The calibration curves for the 1-, 3-, and 5-year OS rates were also largely overlapped with the standard lines (Fig. 7D).

Correlations between the risk score and expression of immune-related genes and immune subtypes of HCC

Considering the important role of the immunosuppressive microenvironment in tumors, we sought to find out whether the risk score based on m6A-related lncRNAs was correlated with the expression of immune-negative regulatory genes. The immune-negative regulatory genes were first extracted from the list of immune-related genes obtained from the TIP database, and a total of 54 genes were identified. Then, the expression of these genes was compared between the high and low-risk patients in the training and validation groups, respectively. Consequently, 20 and 26 immune-negative regulatory genes displayed significant differences in the training and validation groups, respectively (Fig. S3A–B). Only 19 differentially expressed genes were found in both groups, including PDCD1, CTLA4, HAVCR2 (TIM3), ARG1, CCL28, CXCL8, DNMT1, EZH2, HAVCR1, ICAM1, LAIR1, LGALS9, MICA, MICB, SMC3, TGFB1, TNFRSF14, VEGFA, and VTCN1. PD-1, CTLA-4, and TIM3 expression in patients with different risk scores, and the correlations between these genes and risk score are illustrated both in the training (Fig. 8A) and validation (Fig. 8B) groups. Higher expression levels of PD-1, CTLA-4, and TIM3 were found in the high-risk group and positive correlations were also identified between these immune-related genes and risk scores.

Correlations between the risk score and immune-related genes, immune subtypes, tumor stemness, and drug susceptibility in HCC.
Fig. 8  Correlations between the risk score and immune-related genes, immune subtypes, tumor stemness, and drug susceptibility in HCC.

(A–B) Correlation between the risk score and the expression of PD-1, CTLA-4, and TIM3 in the training and validation groups. (C–D) Distribution of the risk score and immune subtypes in HCC. *p<0.05, **p<0.01, ***p<0.001. CTLA-4, cytotoxic T lymphocyte-associated antigen-4; HCC, hepatocellular carcinoma; lncRNAs, long noncoding RNAs; m6A, N6-methyladenosine; PD-1, programmed cell death 1; TCGA, The Cancer Genome Atlas; TIM3, T-cell immunoglobulin and mucin domain 3.

Furthermore, the distribution of risk scores was compared in different immune subtypes of HCC. The subtype of C1, C2, C3, C4, and C6 is the wound healing type, interferon-gamma dominant type, inflammatory type, lymphocyte depleted type, and transforming growth factor-beta dominant type, respectively. The results showed that the risk scores of patients belonging to C1 and C2 immune subtypes were significantly higher than those belonging to C3 and C4 immune subtypes (Fig. 8C). In the high-risk group, the proportions of C1 and C2 subtypes were higher, whereas the C3 subtype was significantly less (Fig. 8D).

Discussion

HCC is one of the leading causes of cancer-related death.1 A high recurrence rate restricts its prognosis critically, while no effective prevention methods have been reported. With the advancement of systematic therapies, increasing drugs have emerged after sorafenib, but they are still unsatisfying in improving survival.3 Immunotherapy based on immune checkpoints has also been recently used in HCC.5 However, its efficiency is limited to less than 20%. Therefore, there is a desperate and urgent need for comprehensive research on the mechanisms of development, progression, recurrence, metastasis, and resistance of HCC. This will help reveal key regulatory pathways or networks in cancer, as well as promote the development and improvement of corresponding therapies. Besides, this provides the potential to break the current treatment dilemma.

The m6A modification of RNA provides a new mechanism of gene expression regulation, which is mainly involved in the regulation of RNA splicing, stability, transporting, and translation.24,25 Dysregulation of m6A-related regulators, including different “writers”, “readers”, and “erasers”, are involved in all stages of tumor development and progression.9,10 In the present study, the expressions of 24 m6A-related genes were found dysregulated in HCC tumors, which were also evaluated by other researchers. Zhang et al.26 discovered that 11 out of 13 m6A-related genes were significantly increased in HCC. Most of the differentially expressed genes in their study were also included in our results. However, because of much more m6A-related genes enrolled by reviewing the latest literature, we obtained additional differentially expressed genes. This indicates that research on the modification of m6A is rapidly progressing and increasing discoveries are being revealed.

In addition to the modifications on mRNA, increasing evidence has confirmed a close correlation between m6A modification and ncRNAs, specifically for lncRNAs.11,20–22 Previous studies have demonstrated that m6A methylation potentially promote the lncRNA X-inactive specific transcript (XIST)-mediated transcriptional repression.27 The lncRNA of metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) could be modified by METTL3 and facilitate the aggravation of renal fibrogenesis.28 Besides, the m6A methylation could also modify lncRNAs to participate in carcinogenesis and metastasis of various tumors, including head and neck cancer, esophageal cancer, lung cancer, gastric cancer, pancreatic cancer, colorectal cancer, and HCC.29 The m6A level of XIST was found to be increased by METTL14, which could inhibit the invasion of colorectal cancer.30 However, LINC00958 was found to be upregulated by METTL3 and promoted the tumor migration and invasion in HCC, and METTL16 could interact directly with MALAT1 to facilitate the tumor proliferation.19,21 Hence, these findings demonstrated significant crosstalk between lncRNAs and m6A in tumorigenesis and development, which provided novel mechanisms and potential targets for cancer regulation and treatment.

Based on the significant roles of m6A in cancer and the close interactions between m6A and lncRNAs, we investigated the potential regulatory mechanisms and prognostic values of m6A-related lncRNAs in HCC on the whole with the transcriptome data from TCGA. After rigorous selection, six critical m6A-related lncRNAs (NRAV, SNHG3, KDM4A-AS1, AC074117.1, AC025176.1, and AL031985.3) were identified and used to construct a novel risk score model, which could distinguish the prognosis of patients with high accuracy. Notably, the risk model was also identified as an independent prognostic factor for HCC. Then, a prognostic nomogram providing high accuracy for 1-, 3-, and 5-year OS prediction of HCC was established. This extended the accessibility of m6A-related lncRNAs for clinical use.

Unexpectedly, significant correlations were observed between m6A-related lncRNAs and immune-related genes and immune subtypes of HCC. Considering the essential support of the TME in tumorigenesis and progression, the malignant behavior of the tumor is also regulated by other cells and factors in the microenvironment. The interaction between cancer cells, immune cells, and the immunosuppressive microenvironment regulates the tumor progression at all stages. Therefore, the changes of tumor biological behaviors mediated by m6A modification may play a part via immune-related mechanisms.31 In the present study, the expression of immune checkpoints, including PD-1, PD-L1, CTLA-4, LAG3, TIM3, and TIGIT, were found with significant differences between the two clusters divided based on the 29 m6A-related lncRNAs. The expression levels of all the immune checkpoint genes were higher in cluster 2, which had higher expression of m6A-related lncRNAs. When the risk score model based on the six critical m6A-related lncRNAs was further constructed, the expression profiles of immune checkpoint genes with PD-1, CTLA-4, and TIM3 were highly expressed in the high-risk group. These findings indicated the potential role of m6A in the regulation of immune checkpoints expression in the TME. Furthermore, the distribution of immune subtypes in the high- and low-risk groups was significantly different. The proportions of C1 (wound healing type) and C2 (interferon-gamma dominant type) were higher in the high-risk group, whereas the C3 proportion (inflammatory type) was higher in the low-risk group of HCC. These findings provided potential guidance for the selection of immune checkpoint inhibitor-based therapy for HCC patients.

As for the six critical m6A-related lncRNAs, the roles of NRAV and SNHG3 in cancer have been investigated in some studies. Xu et al.32 discovered that NRAV was upregulated in HCC tumors and could serve as a poor prognostic factor for HCC patients, which was consistent with our findings. However, the NRAV was identified as an m6A-related lncRNA of HCC in our study, whereas Xu et al.32 reported an immune-related lncRNA. This indicated the complex roles of lncRNA in tumors. SNHG3 has been identified as a tumor-promoting lncRNA in various cancers.33–36 Wu et al.33 found that SNHG3 was upregulated in HCC and could promote cell proliferation, migration, and invasion. Zhang et al.37 found that SNHG3 could induce epithelial-mesenchymal transition and sorafenib resistance in HCC. Notably, a recent study revealed that SNHG3 could epigenetically regulate its neighboring MED18 gene methylation by binding to EZH2 and inhibiting its expression in gastric cancer.38 Besides, Zhang et al.34 found that the increased expression of SNHG3 induced by the platinum treatment could regulate the m6A level by targeting METTL3 in esophageal cancer. These findings indicated that the important role of SNHG3 in epigenetic regulation of cancer progression and resistance. Yet, the role of SNHG3 in the regulation of m6A in HCC remains unreported. As such, additional studies on these newly discovered lncRNAs are essential.

This work has several limitations worth mentioning. First, the data used and the results obtained in this study are based on the transcriptome data from TCGA, and therefore should be further verified by local experiments. Secondly, the training and validation groups were divided from the TCGA cohort, which required an external cohort for further validation. Thirdly, the transcriptome data of HCC in TCGA are mainly derived from Americans, which will cause inevitable selection bias.

In conclusion, a novel risk score model was developed based on the m6A-related lncRNAs with significant potential for HCC prognosis prediction. Our findings provide new perspectives for the potential mechanisms of m6A-related lncRNAs in regulating the immune microenvironment of HCC.

Supporting information

Fig. S1

Schematic diagram of the potential mechanisms of RNA N6-methyladenosine (m6A) methylation.

(TIF)

Fig. S2

Expression profile and function enrichment of m6A-DEGs in HCC.

(A) Heatmap of 24 differentially expressed m6A-related genes in HCC tumors and normal controls. (B) Correlation among these m6A-related genes. (C) Function enrichment of differentially expressed m6A-related genes in HCC. (D) Pathway analysis of differentially expressed m6A-related genes in HCC. HCC, hepatocellular carcinoma; m6A, N6-methyladenosine; m6A-DEGs, differentially expressed m6A-related genes.

(TIF)

Fig. S3

Expression profile of immune-negative regulatory genes in HCC patients with different risk scores.

(A) Differentially expressed immune-negative regulatory genes (n=20) in the training group. (B) Differentially expressed immune-negative regulatory genes (n=26) in the validation group. *p<0.05, **p<0.01, ***p<0.001. HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas.

(TIF)

Table S1

m6A-related genes.

(DOCX)

Table S2

Coefficients of m6A-related lncRNAs in the risk score model.

(DOCX)

Table S3

Clinical characteristics of HCC patients in the training and validation groups.

(DOCX)

Abbreviations

AUC: 

area under the curve

C-index: 

concordance index

CDF: 

cumulative distribution function

circRNAs: 

circular RNAs

CTLA-4: 

cytotoxic T lymphocyte-associated antigen-4

GSEA: 

gene set enrichment analysis

HCC: 

hepatocellular carcinoma

HR: 

hazard ratio

LASSO: 

least absolute shrinkage and selection operator

LAG3: 

lymphocyte-activation gene 3

lncRNAs: 

long noncoding RNAs

m6A: 

N6-methyladenosine

m6A-DEGs: 

differentially expressed m6A-related genes

MALAT1: 

metastasis-associated lung adenocarcinoma transcript 1

miRNAs: 

microRNAs

ncRNAs: 

noncoding RNAs

OS: 

overall survival

PD-1: 

programmed cell death 1

PD-L1: 

PD-1 ligand 1

ROC: 

receiver operating characteristic

TCGA: 

The Cancer Genome Atlas

TIGIT: 

T cell immunoreceptor with immunoglobulin and ITIM domain

TIM3: 

T-cell immunoglobulin and mucin domain 3

TIP: 

Tracking Tumor Immunophenotype

TME: 

tumor microenvironment

XIST: 

X-inactive specific transcript

Declarations

Acknowledgement

The results shown here are in whole or part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga. We also thank Home for Researchers editorial team (www.home-for-researchers.com) for providing assistance in preparing the manuscript.

Data sharing statement

The datasets presented in this study can be found in online repositories. The names of the repositories can be found in the article.

Funding

The work was supported in part by grants from the Guangdong Natural Science Foundation (Nos. 2015A030313038, 2015A030312013), Science and Technology Program of Guangzhou City (No. 201607010024), and Guangdong Key Laboratory of Liver Disease Research (No. 2017B030314027).

Conflict of interest

The authors have no conflict of interests related to this publication.

Authors’ contributions

Study concept and design (TD, YY, GW), acquisition of data (JL, LY, HY), analysis and interpretation of data (TD, MD, WL, HL), drafting of the manuscript (TD, JL, LY), critical revision of the manuscript for important intellectual content (WL, HL, YY, GW), and administrative, technical, or material support, study supervision (YY, GW).

References

  1. Llovet JM, Kelley RK, Villanueva A, Singal AG, Pikarsky E, Roayaie S, et al. Hepatocellular carcinoma. Nat Rev Dis Primers 2021;7(1):6 View Article PubMed/NCBI
  2. Villanueva A. Hepatocellular carcinoma. N Engl J Med 2019;380(15):1450-1462 View Article PubMed/NCBI
  3. Llovet JM, Montal R, Sia D, Finn RS. Molecular therapies and precision medicine for hepatocellular carcinoma. Nat Rev Clin Oncol 2018;15(10):599-616 View Article PubMed/NCBI
  4. Liu Z, Lin Y, Zhang J, Zhang Y, Li Y, Liu Z, et al. Molecular targeted and immune checkpoint therapy for advanced hepatocellular carcinoma. J Exp Clin Cancer Res 2019;38(1):447 View Article PubMed/NCBI
  5. Kalasekar SM, Garrido-Laguna I, Evason KJ. Immune checkpoint inhibitors in novel combinations for hepatocellular carcinoma. Hepatology 2021;73(6):2591-2593 View Article PubMed/NCBI
  6. Zhou Z, Lv J, Yu H, Han J, Yang X, Feng D, et al. Mechanism of RNA modification N6-methyladenosine in human cancer. Mol Cancer 2020;19(1):104 View Article PubMed/NCBI
  7. Lu J, Qian J, Yin S, Zhou L, Zheng S, Zhang W. Mechanisms of RNA N(6)-methyladenosine in hepatocellular carcinoma: from the perspectives of etiology. Front Oncol 2020;10:1105 View Article PubMed/NCBI
  8. Pan XY, Huang C, Li J. The emerging roles of m(6)A modification in liver carcinogenesis. Int J Biol Sci 2021;17(1):271-284 View Article PubMed/NCBI
  9. Jiang X, Liu B, Nie Z, Duan L, Xiong Q, Jin Z, et al. The role of m6A modification in the biological functions and diseases. Signal Transduct Target Ther 2021;6(1):74 View Article PubMed/NCBI
  10. Zaccara S, Ries RJ, Jaffrey SR. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol 2019;20(10):608-624 View Article PubMed/NCBI
  11. 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(3):270-288 View Article PubMed/NCBI
  12. Ma S, Chen C, Ji X, Liu J, Zhou Q, Wang G, et al. The interplay between m6A RNA methylation and noncoding RNA in cancer. J Hematol Oncol 2019;12(1):121 View Article PubMed/NCBI
  13. Chen Y, Lin Y, Shu Y, He J, Gao W. Interaction between N(6)-methyladenosine (m(6)A) modification and noncoding RNAs in cancer. Mol Cancer 2020;19(1):94 View Article PubMed/NCBI
  14. Yang X, Hu X, Liu J, Wang R, Zhang C, Han F, et al. N6-methyladenine modification in noncoding RNAs and its function in cancer. Biomark Res 2020;8(1):61 View Article PubMed/NCBI
  15. He RZ, Jiang J, Luo DX. The functions of N6-methyladenosine modification in lncRNAs. Genes Dis 2020;7(4):598-605 View Article PubMed/NCBI
  16. Kopp F, Mendell JT. Functional classification and experimental dissection of long noncoding RNAs. Cell 2018;172(3):393-407 View Article PubMed/NCBI
  17. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3′ UTRs and near stop codons. Cell 2012;149(7):1635-1646 View Article PubMed/NCBI
  18. 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(7540):560-564 View Article PubMed/NCBI
  19. Brown JA, Kinzig CG, DeGregorio SJ, Steitz JA. Methyltransferase-like protein 16 binds the 3′-terminal triple helix of MALAT1 long noncoding RNA. Proc Natl Acad Sci U S A 2016;113(49):14013-14018 View Article PubMed/NCBI
  20. 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(18):4612-4626 View Article PubMed/NCBI
  21. 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(1):5 View Article PubMed/NCBI
  22. Dai F, Wu Y, Lu Y, An C, Zheng X, Dai L, et al. Crosstalk between RNA m(6)A modification and non-coding RNA contributes to cancer growth and progression. Mol Ther Nucleic Acids 2020;22:62-71 View Article PubMed/NCBI
  23. Xu L, Deng C, Pang B, Zhang X, Liu W, Liao G, et al. TIP: a web server for resolving tumor immunophenotype profiling. Cancer Res 2018;78(23):6575-6580 View Article PubMed/NCBI
  24. Lee Y, Choe J, Park OH, Kim YK. Molecular mechanisms driving mRNA degradation by m(6)A modification. Trends Genet 2020;36(3):177-188 View Article PubMed/NCBI
  25. Tian S, Lai J, Yu T, Li Q, Chen Q. Regulation of gene expression associated with the N6-methyladenosine (m6A) enzyme system and its significance in cancer. Front Oncol 2020;10:623634 View Article PubMed/NCBI
  26. Zhang L, Qiao Y, Huang J, Wan D, Zhou L, Lin S, et al. Expression pattern and prognostic value of key regulators for m6A RNA modification in hepatocellular carcinoma. Front Med (Lausanne) 2020;7:556 View Article PubMed/NCBI
  27. Patil DP, Chen CK, Pickering BF, Chow A, Jackson C, Guttman M, et al. m(6)A RNA methylation promotes XIST-mediated transcriptional repression. Nature 2016;537(7620):369-373 View Article PubMed/NCBI
  28. Liu P, Zhang B, Chen Z, He Y, Du Y, Liu Y, et al. m(6)A-induced lncRNA MALAT1 aggravates renal fibrogenesis in obstructive nephropathy through the miR-145/FAK pathway. Aging (Albany NY) 2020;12(6):5280-5299 View Article PubMed/NCBI
  29. Yi YC, Chen XY, Zhang J, Zhu JS. Novel insights into the interplay between m(6)A modification and noncoding RNAs in cancer. Mol Cancer 2020;19(1):121 View Article PubMed/NCBI
  30. Yang X, Zhang S, He C, Xue P, Zhang L, He Z, et al. METTL14 suppresses proliferation and metastasis of colorectal cancer by down-regulating oncogenic long non-coding RNA XIST. Mol Cancer 2020;19(1):46 View Article PubMed/NCBI
  31. Shulman Z, Stern-Ginossar N. The RNA modification N(6)-methyladenosine as a novel regulator of the immune system. Nat Immunol 2020;21(5):501-512 View Article PubMed/NCBI
  32. Xu Q, Wang Y, Huang W. Identification of immune-related lncRNA signature for predicting immune checkpoint blockade and prognosis in hepatocellular carcinoma. Int Immunopharmacol 2021;92:107333 View Article PubMed/NCBI
  33. Wu J, Liu L, Jin H, Li Q, Wang S, Peng B. LncSNHG3/miR-139-5p/BMI1 axis regulates proliferation, migration, and invasion in hepatocellular carcinoma. Onco Targets Ther 2019;12:6623-6638 View Article PubMed/NCBI
  34. Zhang M, Bai M, Wang L, Lu N, Wang J, Yan R, et al. Targeting SNHG3/miR-186-5p reverses the increased m6A level caused by platinum treatment through regulating METTL3 in esophageal cancer. Cancer Cell Int 2021;21(1):114 View Article PubMed/NCBI
  35. Wan Q, Tang M, Sun SL, Hu J, Sun ZJ, Fang YT, et al. SNHG3 promotes migration, invasion, and epithelial-mesenchymal transition of breast cancer cells through the miR-186-5p/ZEB1 axis. Am J Transl Res 2021;13(2):585-600 View Article PubMed/NCBI
  36. Xu B, Mei J, Ji W, Bian Z, Jiao J, Sun J, et al. LncRNA SNHG3, a potential oncogene in human cancers. Cancer Cell Int 2020;20(1):536 View Article PubMed/NCBI
  37. Zhang PF, Wang F, Wu J, Wu Y, Huang W, Liu D, et al. LncRNA SNHG3 induces EMT and sorafenib resistance by modulating the miR-128/CD151 pathway in hepatocellular carcinoma. J Cell Physiol 2019;234(3):2788-2794 View Article PubMed/NCBI
  38. Xuan Y, Wang Y. Long non-coding RNA SNHG3 promotes progression of gastric cancer by regulating neighboring MED18 gene methylation. Cell Death Dis 2019;10(10):694 View Article PubMed/NCBI