Significance of m 6A Methyltransferase‑Related lncRNAs in the Prognosis and Immunotherapy of Hepatocellular Carcinoma

doi:10.21203/rs.3.rs-165071/v1

Download PDF

Research

Significance of m 6A Methyltransferase‑Related lncRNAs in the Prognosis and Immunotherapy of Hepatocellular Carcinoma

https://doi.org/10.21203/rs.3.rs-165071/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Hepatocellular carcinoma (HCC) is one of the leading causes of mortality worldwide. N6-methyladenosine (m6A) methyltransferase, has been proved to act as an oncogene in several human cancers. However, little is known about its relationship with the long non-coding RNAs (lncRNAs) that remains elusive in HCC.

Methods: We comprehensively integrated gene expression data acquired from 371 HCC and 50 normal tissues in The Cancer Genome Atlas (TCGA) database. Differentially expressed protein-coding genes (DE-PCGs)/lncRNAs (DE-lncRs) analysis and univariate regression & Kaplan-Meier (K-M) analysis was performed to identify m6A methyltransferase‑related lncRNAs that were related to overall survival (OS). m6A methyltransferase‑related lncRNA signature was constructed using the Least absolute shrinkage and selection operator (LASSO) Cox regression analyses. Furthermore, Cox regression analysis was applied to identify independent prognostic factors in HCC. The signature was validated in an internal validation set. Finally, the correlation analysis between gene signature and immune cells infiltration was also investigated via single-sample Gene Set Enrichment Analysis (ssGSEA) and immunotherapy response was calculated through Tumor Immune Dysfunction and Exclusion (TIDE) algorithm.

Results: A total of 21 m6A methyltransferase-related lncRNAs were screened out according to Spearman correlation analysis with the immune score (|R| > 0.3, P < 0.05). We selected 3 prognostic lncRNAs to construct m6A methyltransferase-related lncRNA signature through univariate and LASSO Cox regression analyses. The univariate and multivariate Cox regression analyses demonstrated that the lncRNAs signature was a robust independent prognostic factor in OS prediction with high accuracy. The GSEA also suggested that the m6A methyltransferase-related lncRNAs were involved in the immune-related biological processes and pathways which were very well-known in the context of HCC tumorigenesis. Besides, we found that the lncRNAs signature was strikingly correlated with the tumor microenvironment (TME) immune cells infiltration and expression of critical immune checkpoints. Finally, results from the TIDE analysis revealed that the m6A methyltransferase-related lncRNAs could efficiently predict the clinical response of immunotherapy in HCC.

Conclusion: Together, our study screened potential prognostic m6A methyltransferase related lncRNAs and established a novel m6A methyltransferase-based prognostic model of HCC, which not only provides new potential prognostic biomarkers and therapeutic targets but also deepens our understanding of tumor immune microenvironment status and lays a theoretical foundation for immunotherapy.

Surgery

Oncology

Hepatocellular carcinoma

N6-methyladenosine

LncRNAs

Immune cells infiltration

Immune checkpoint

Worldwide, hepatocellular carcinoma (HCC) has become one of the most common malignancies and the second leading cause of cancer-related death (1). HCC is the main pathological type of primary liver cancer, with a global incidence of 500,000 new cases and more than 700,000 deaths each year (2). The major risk factors for HCC are chronic hepatitis B virus (HBV) and hepatitis c virus (HCV) infections, alcoholic liver disease, diabetes, and nonalcoholic fatty liver disease (3). The current curative treatments for early-stage HCC are surgery, thermal ablation, radiofrequency ablation or liver transplantation (4–7). However, a large proportion of patients will have recurrence or distant metastasis after surgery (8). For the more than 70% of patients diagnosed with advanced stage, treatments can only provide limited therapeutic benefit for a small subset of patients. Thus, elucidating the molecular mechanisms of HCC and identifying new molecular targets are essential for its diagnosis and treatment.

N6-methylandenosine (m6A), the most popular internal modification in eukaryotic messenger RNAs (mRNAs), microRNAs (miRNAs) and long non-coding RNAs (lncRNAs), plays a critical part in RNA splicing, stability, export and translation (9, 10). The dynamic modification of m6A is regulated by methyltransferases (m6A “writers”), demethylases (m6A “erasers”), and binding proteins (m6A “readers”) (11). Numerous studies have described that m6A regulators are closely related to the occurrence and progression of a variety of cancers. It has been reported that m6A regulators are key participants in the malignant progression of gliomas and have potential roles in both prognosis and treatment strategy formulation (12). Recently, many studies have revealed that abnormal m6A modification is involved in the development and progression of HCC (13–18). METTL14 inhibits the metastatic potential of HCC by regulating N⁶-methyladenosine-dependent primary miRNA processing (13). Upregulation of YTHDF2 suppressed cell proliferation, tumor growth and activation of mitogen-activated protein kinase kinase (MAPKK/MEK) and extracellular signal-regulated kinase (ERK) in HCC cells (14). KIAA1429 facilitated migration and invasion of HCC by altering m6A modification of ID2 mRNA (16).

LncRNAs are families of non-coding molecular transcripts with over 200 nucleotides in length and have been recognized as having important functional roles in various human diseases (19, 20). The aberrant expression of lncRNAs is also closely related to the malignancy of tumors, and it is no exception in HCC (21–25). For example, lncRNA MIAT was reported to promote proliferation and invasion of HCC cells via sponging miR-214 (21); Upregulation of lncRNA SNHG16 inhibited HCC cell proliferation and chemotherapy-resistance via functionally sponging hsa-miR-93 (22); And lncRNA HULC can trigger autophagy by stabilizing Sirt1 and attenuate the sensitivity of HCC cells to chemotherapeutic agents (23). However, the role of m6A regulators in the dysregulation of lncRNAs in cancer development is still unclear, and there are few studies focusing on the relationship between m6A modification and lncRNA-dependent HCC. Therefore, understanding the relationship between the modification of lncRNAs by m6A and the progression of HCC may contribute to identify biomarkers that may serve as therapeutic targets.

In this study, we constructed an m6A-related lncRNA prognostic signature (LINC02362, SNHG20, and SNHG6) which acted as an independent prognostic factor in overall survival (OS) prediction with high accuracy. Our results showed that this model was involved in immune-related response processes and pathways that played a key role in HCC tumorigenesis, and was highly correlated with tumor microenvironment (TME) immune cells infiltration and immunotherapy responses. Our study constructed a novel m6A -based prognostic model that has potential prognostic value for HCC patients and may contribute to personalized immunotherapy counseling.

Acquisition of gene expression and clinical data

Normalized RNA-sequencing data and the associated clinical information of the HCC samples were downloaded from the TCGA database. They included 374 tumor samples and 50 normal tissue samples. Three samples with metastasis were eliminated. Thus, 371 tumor tissues were included in the study. Normalized gene expression data for the TCGA-HCC database were log₂-transformed for further analysis (Supplementary Fig. 1). Samples from the TCGA database were divided randomly into a training set (260 tumor samples, 35 normal samples) and an internal validation set (111 tumor samples, 15 normal samples) at a ratio of 7:3. Clinical data such as gender, age, pathologic stage, Tumor-Node-Metastasis (TNM) stage, treatment type, treatment or therapy, prior malignancy, and survival time were also obtained from the TCGA database.

Analysis of differentially expressed PCGs and lncRNAs

In the training dataset, the differentially expressed RNAs, including DE-PCGs and differentially expressed DE-lncRs, between tumor and normal samples, were analyzed using the EdgeR package version 3.8.5, which adopted an overdispersed Poisson model and an empirical Bayesian approach to improve the reliability of the prediction (26). The thresholds for DE-PCG and DE-lncR selection were based on a P ≤ 0.05 and |logFC | ≥ 1. Hierarchical cluster analysis served to analyze sample similarity.

Identification of m6A methyltransferase‑related lncRNAs

The m6A methyltransferase-related genes were obtained from multiple literatures (27–29). Four differentially expressed m6A methyltransferase-related genes (DE-m6A methyltransferase genes) were obtained through the intersection analysis with DE-PCGs. The ssGSEA was used to calculate the m6A methyltransferase score of each sample in the TCGA-HCC training cohort (30, 31). Then the m6A methyltransferase‑related lncRNAs were identified for a high correlation with the immune score (|R| > 0.3, P < 0.05) based on Spearman correlation analysis. K-M survival analyses were utilized to screen out the prognosis-related lncRNAs (P < 0.2). After merging the m6A methyltransferase-related and prognosis-related lncRNAs, the remained selected lncRNAs were considered to be candidate m6A methyltransferase-related lncRNAs.

Establishment and validation of the prognostic gene signature

Prognosis associated factors were selected by univariate Cox regression. Subsequently, we performed Cox regression analysis combined with LASSO regression to establish a risk model and the penalty regularization parameter lambda (λ) was chosen through the cross-validation routine with an n-fold equal to 10 by using R package glmnet (32). Meanwhile, Lambda.min was identified to pick out the variables. Finally, 3 m6A methyltransferase‑related lncRNAs were enrolled in risk cox regression and survival analysis, scatter diagram, and heatmap were performed in R software according to the risk score for each patient. The time-dependent receiver operating characteristic (ROC) curves were utilized to assess the accuracy in prognosis prediction of gene signature and the area under curve (AUC) for 1-year, 2-year, 3-year, 4-year, and 5-year OS was measured using package ‘survivalROC” in R (33). The prognostic gene signature was verified in the internal validation set. Moreover, univariate and multivariate Cox regression analyses were performed to demonstrate whether the risk score was an independent prognostic factor.

Construction of a predictive nomogram

The risk score and clinical features were merged to find independent prognostic factors through univariate and multivariate Cox regression analyses and visualized through package ‘forestplot’ in R. Then the nomogram integrating with selected independent prognostic factors were established through package ‘rms’, ‘nomogramEx’ and ‘regplot’ in R (34). Following that, calibration curves were plotted to evaluate the concordance between actual and predicted survival. Furthermore, decision curve analysis (DCA) was used to see whether our nomogram was serviceable as the ideal model.

Gene set enrichment analyses

We conducted GSEA to uncover the signaling pathways and biological processes in which differentially expressed genes (DEGs) between high- and low-risk subgroups were enriched. A total of 454 DEGs (160 upregulated and 294 downregulated) were identified as differentially expressed in high-risk compared with low-risk groups (Supplementary Fig. 2, Supplementary Table S1). We employed GSEA to obtain information for Gene Ontology (GO) including biological processes (BP), the cellular component (CC), and molecular function (MF). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis served to annotate the potential functions. A significance level of P < 0.05 was set as the cutoff criteria and the plots were constructed by the gplots package in R software.

ESTIMATE algorithm

Stromal and immune scores were calculated by the ESTIMATE package in R (version 2.15.3). ESTIMATE is an algorithmic tool for predicting tumor purity, which uses the gene expression profiles of 141 immune genes and 141 stromal genes to generate ESTIMATE score (30). The presence of infiltrating immune cells in high- and low-risk groups was calculated with ssGSEA software.

Statistical analysis

All analyses were performed using R version 3.5.1. Unless otherwise noted, P < 0.05 was considered to be significant.

Identification of DE-PCGs and -lncRNAs in HCC

With the aim of identifying PCGs and lncRNAs potentially involved in the HCC pathogenesis, in the training set of 295 (260 HCC and 35 healthy subjects), transcriptional profiles of HCC patients and healthy subjects were compared. After a robust filtering procedure (P ≤ 0.05 and |logFC| ≥ 1), 1916 PCGs were significantly modulated (Fig. 1A; Supplementary Table S2). More than one hundred lncRNAs were detectable, of which 80.6% were up-regulated and 19.4% were down-regulated (Fig. 1B; Supplementary Table S3). The DE-PCGs and -lncRNAs with similar expression levels were clustered using the systematic cluster analysis (Figs. 1C and D). As displayed in Table 1, the most up-regulated PCG was AKR1B10 with 3.84-logFC, and the most down-regulated PCG was HAMP with − 5.90-logFC. Besides, the most up- and down-regulated lncRNA were ST8SIA6-AS1 (2.60-logFC) and LINC01093 (-3.60-logFC), respectively. The expression patterns of the top 10 up-regulated and down-regulated PCGs and lncRNAs were illustrated in Figs. 1E and F (Supplementary Tables S4 and S5).

Table 1
Statistical Analysis of All of the Differentially Expressed PCGs and lncRNAs
DE RNAs	Total No.	No. of Upregulated	No. of Downregulated	The Most Upregulated (logFC)	The Most Downregulated (logFC)
PCG	1916	1500	416	AKR1B10 (3.84)	HAMP (-5.90)
LncRNA	103	83	20	ST8SIA6-AS1 (2.60)	LINC01093 (-3.59)
PCG: protein-coding gene;DE RNAs: Differentially expressed RNAs

Identification of prognostic m6A methyltransferase‑related lncRNAs

Transcription data of 260 patients with full clinical information were retrieved from the TCGA-HCC training cohort. The 45 prognostic-related lncRNAs were screened out through K–M survival analyses (Fig. 2A; Table 2; Supplementary File 1). Furthermore, the m6A methyltransferase scores in each sample were calculated via ssGSEA according to the reference of the DE-m6A methyltransferase (METTL3, VIRMA, IGF2BP1, and IGF2BP2) gene sets derived from the intersection analysis of 20 m6A methyltransferase genes and DE-PCGs (Fig. 2B). After measuring with Spearman correlation analyses, 21 lncRNAs were recognized as m6A methyltransferase-related lncRNAs (|R| > 0.3, P < 0.05; Fig. 2C). At last, we have got 8 candidate m6A methyltransferase-related lncRNAs (LINC01093, LINC02362, SNHG20, SNHG17, ZFAS1, SNHG6, SNHG7, and GAS5) by the intersection analysis for further research (Fig. 2D).

Table 2
The 45 prognostic-related lncRNAs
Prognostic-related lncRNA	p-value
CASC9	0.0148133
CH17.340M24.3	0.0268725
CRNDE	0.1003839
CYTOR	0.0001119
DUXAP8	0.0147407
ELF3.AS1	0.0517343
FAM99A	0.0046928
FAM99B	0.037008
FOXD2.AS1	0.0302527
GAS5	0.1616622
GSEC	0.0005606
HNF1A.AS1	0.1288273
HNF4A.AS1	0.0779326
LINC00511	0.0019846
LINC01093	0.0522651
LINC01370	0.1255833
LINC01419	0.1944095
LINC01554	0.0004055
LINC01703	0.0238661
LINC02037	0.0481163
LINC02055	0.1514205
LINC02163	0.1866274
LINC02362	0.1794657
LINC02377	0.1152222
LINC02506	0.1772376
LUCAT1	0.0278805
MAPKAPK5.AS1	0.0016472
MINCR	0.007126
MIR4435.2HG	0.0030908
MYLK.AS1	0.0309763
NALT1	0.0858859
PCAT6	0.0628995
PITPNA.AS1	0.0554125
PRR34.AS1	0.1776464
PVT1	0.1908601
PXN.AS1	0.0448655
RALY.AS1	0.0630094
SNHG17	0.1073125
SNHG20	0.0105894
SNHG6	0.1046746
SNHG7	0.049813
TSPEAR.AS2	0.1091516
UBR5.AS1	0.0293862
ZFAS1	0.0361043
ZFPM2.AS1	0.0013408

Table 3
The univariate cox regression analysis of eight m6A methyltransferase-related lncRNAs
Gene ID	Coefficient	HR (95% CI for HR)	wald.test	z	p.value
LINC01093	-0.14	0.87 (0.76-1)	3.7	-1.9	0.055
LINC02362	-0.16	0.85 (0.75–0.97)	6.1	-2.5	0.014
SNHG20	0.44	1.6 (1.2–2.1)	8.6	2.9	0.0033
SNHG17	0.29	1.3 (1-1.7)	5.5	2.4	0.019
ZFAS1	0.28	1.3 (1.1–1.6)	6.8	2.6	0.0093
SNHG6	0.27	1.3 (1.1–1.6)	7	2.6	0.0084
SNHG7	0.24	1.3 (1-1.6)	4.5	2.1	0.034
GAS5	0.2	1.2 (1-1.5)	4	2	0.046

Establishment of m6A methyltransferase-related lncRNAs signature

The univariate Cox regression analysis demonstrated that LINC01093 (P = 0.055) and LINC02362 (P = 0.014) were protective genes with hazard ratio (HR) less than 1, whereas SNHG20 (P = 0.0033), SNHG17 (P = 0.019), ZFAS1 (P = 0.0093), SNHG6 (P = 0.0084), SNHG7 (P = 0.034), and GAS5 (P = 0.046) were risky genes with HR greater than 1 (Fig. 3A). Following this, eight significant m6A methyltransferase‑related lncRNAs (P < 0.2) were subjected to LASSO modeling (Figs. 3B and C). Then, three m6A methyltransferase-related lncRNAs were selected based on lambda.min values, and the HR values of the three candidate lncRNAs were calculated (Fig. 3D; Supplementary Table S6).

We stratified patients into low-risk or high-risk groups by median cut-off based on the risk score. K-M survival curve indicated that the low-risk patients lived longer than high-risk patients in the training cohort (P = 0.00064) (Figs. 4A and B). Therefore, time-dependent ROC analysis showed an appropriate accuracy of the risk score in predicting OS in training cohort and area under the ROC curve (AUC) was 0.633 at 1 year, 0.636 at 2 years, 0.651 at 3 years, 0.663 at 4 years, and 0.638 at 5 years (Fig. 4C). Moreover, these prognostic factors were further validated in the testing cohort to assess the robust prediction value of the risk score. We found that results in the testing cohort were consistent with the outcome of the training cohort, indicating all the high-risk patients were associated with poorer prognosis (Fig. 4D and E). In the testing cohort, the significant prognostic value was P = 0.031 and AUC with 1-, 2-, 3-, 4-, and 5 years were 0.708, 0.674, 0.635, 0.603, and 0.611, respectively (Fig. 4F). Also, the PCA analysis indicated that there was a significant disparity in patients between high- and low-risk groups (Supplementary Fig. 3).

Prognostic risk score displayed strong correlations with clinicopathological features in HCC patients

We used the median of the risk score as the cutoff to define the groups of HCC patients with high and low scores. The heatmap showed the expression of the three m6A methyltransferase‑related lncRNAs in high- and low-risk patients of the training cohort. Significant differences between the high-risk and low-risk groups of patients concerning prior malignancy (P < 0.05), T stage (P < 0.05), and pathologic stage (P < 0.05) were identified (Fig. 5A). We then observed that the risk scores were significantly different between patients categorized by prior malignancy (Fig. 5B), T stage (Fig. 5C), and pathologic stage (Fig. 5D). These results revealed that risk scores were positively associated with the T stage and pathologic stage, and the 'no prior malignancy' subgroup exhibiting the worst prognosis had the highest risk score. Stratification survival analyses were utilized to see whether the risk score could apply in different clinicopathological characteristics. Thus, we found that the risk score could efficiently predict the OS in most subgroups from different clinical features (Supplementary Fig. 4). Similar procedures were performed in the testing set, and the results were almost consistent with the training set (Supplementary Figs. 5 and 6).

The risk score was an independent prognostic factor in HCC

Based on the univariate and multivariate Cox proportional hazards regression analyses, only one independent predictor (risk score) was identified in the HCC. Univariate Cox proportional hazards regression analysis demonstrated that age (P = 0.44), gender (P = 0.27), N stage (P = 0.31), prior malignancy (P = 0.33), treatment type (P = 0.21), and BMI (P = 0.21) had no impact upon OS (Fig. 6A; Supplementary Table S7). Multivariate Cox proportional hazards regression analysis demonstrated that the risk score (HR = 2, P = 0.0037) was associated significantly with OS in patients with HCC (Fig. 6B; Supplementary Table S8). Then, the predictor was incorporated to build the nomogram.

According to the nomogram, every patient will get a total point from the prognostic parameter to estimate the 3- and 5-year survival rates. The higher point indicated the higher mortality of the patients (Fig. 6C). We concluded that the risk score played an important role in patient prognosis. Furthermore, calibration curves demonstrated that the prediction accuracy of our model was similar to the ideal model (Figs. 6D and E). The DCA revealed that the nomogram had an advantage of the risk score alone and displaced a high potential for clinical utility (Figs. 6F and G).

Functional enrichment analyses of the DEGs between the high- and low-risk groups

To explore the underlying molecular mechanisms of the signature, we conducted GSEA to compare 454 DEGs between high- and low-risk groups (risk DEGs) in 371 TCGA-HCC patients of the whole set. We enriched the functions of these mRNA (P < 0.05). As shown in Fig. 7A, ‘negative regulation of transcription by RNA polymerase II’ and ‘microtubule cytoskeleton organization’ were significantly enriched for biological processes (Supplementary Table S9),while for the cellular component was ‘Golgi membrane’, for molecular functions were ‘transcription regulatory region sequence-specific DNA binding’ and ‘RNA polymerase II regulatory region sequence-specific DNA binding’ (Figs. 7B and C; Supplementary Tables S10 and S11). KEGG pathway analysis revealed that the DEGs participated in ‘Herpes simplex virus 1 infection’, ‘Amyotrophic lateral sclerosis’, ‘Human papillomavirus infection’, and ‘MicroRNAs in cancer’ (Fig. 7D; Supplementary Table S12). Furthermore, the DEGs participated in the ‘Axon guidance’ pathway indicating their involvement in the neurological invasion of HCC (Fig. 7E; Supplementary Table S13). Moreover, the GSEA also revealed that immune response and immune system process terms/pathways were highly enriched in the risk DEGs (Supplementary Table S14). All these demonstrated the risk score could be related to the tumor immune microenvironment, which was so important for HCC.

The landscape of TME immune cells infiltration in the two risk subgroups

The tumor immune microenvironment is considered the ‘seventh largest representative characteristic’ of tumor and is composed of innate immune cells, adaptive immune cells, cytokines, and cell surface molecules. These immune microenvironment components constitute a complex regulatory network and play a pivotal role in tumorigenesis and development (35, 36). Considering that the risk score also had a close connection with immune reactions, thus we investigated the difference in immune signatures between the two risk subgroups.

Based on the ImmuneScore, StromalScore, and ESTIMATEScore generated by the ESTIMATE algorithm, as shown in Fig. 8A, the stromal score and ESTIMATE score were significantly negatively correlated with the risk score (P < 0.01). ssGSEA analysis revealed that the ratio of pro-tumor signatures (CD56bright natural killer cell, central memory CD8 T cell, effector memory CD8 T cell, natural killer cell, type 1 T helper cell) and anti-tumor signatures (CD56dim natural killer cell, macrophage, neutrophil) was significantly elevated in the low-risk subgroup, indicating that the low-risk subgroup was interrelated with increased immune/inflammation activity. While the high-risk subgroup had relatively higher activated CD4 T cells (Figs. 8B and C).

Low-risk patients were likely to be more sensitive to immunotherapy

Pioneering investigations revealed that immunotherapy targeting immune-checkpoints and the human leukocyte antigen (HLA) offered great hope for the clinical treatment of human cancers(37). We investigated the correlation between the risk score and the expression levels of several well-known immune-checkpoints and HLA in the TCGA database. Then we found that the expression of almost all immune-checkpoints (CTLA4, HAVCR2, ICOS, PDCD1, and TIGIT) and HLA (HLA-DPB2, HLA-DQA2, HLA-DQB2, HLA-DRB6, and HLA-K) was upregulated in the high-risk patients, which indicated the poor prognosis might be owed to the inducing of the immune checkpoints (Fig. 8D and E). Here, patients in the high-risk group may be more suitable for immunotherapy. Then, we used the TIDE algorithm to predict the likelihood of the risk model for immunotherapy. Interestingly, the low-risk group was more likely to respond to immunotherapy than the high-risk group (Fig. 8F). These data further supported that the low-risk group had a better prognosis and may have a better treatment prospect for the immunotherapy.

Currently, there is lack of effective clinical interventions in HCC, which lead to the high metastasis rate and mortality rate. Thus, it is critical to understand the molecular mechanisms underlying HCC development. There is growing evidence that m6A modification participates in the development and progression of HCC (13–18). METTL3 increases the expression of hepatocellular carcinoma-derived growth factor (HDGF) by upregulating LINC00958, thereby promoting HCC progression and lipogenesis (38). Tian Lan et al demonstrated that KIAA1429 promotes the progression of HCC by modifying lncRNA GATA3 with m6A (39). However, studies on m6A-related lncRNAs in the prognosis of HCC were limited. Therefore, we focused on their interactions to identify potential prognostic biomarkers or therapeutic targets for HCC.

We identified 21 m6A-related lncRNAs from 260 HCC patients, three of which were selected into the m6A-related lncRNA prognostic signature (LINC02362, SNHG20, and SNHG6). LINC02362 was included in a prognostic model for HCC, which was significantly correlated with tumor grade, stage, and T stage (40). Wenhui Tu et al reported that hepatitis B virus x protein promoted the proliferation of HCC cell through the lncRNA SNHG20/PTEN signaling pathway (41). According to reports, upregulated SNHG6 activates the expression of SERPINH1 by competitively binding miR-139-5p, thereby promoting the development of HCC (42). These three lncRNAs have been reported to be associated with HCC, but their interactions with m6A-related genes have not been studied yet. Our results contribute to identify lncRNAs related to m6A modulators, and thus revealing their potential roles in the oncogenesis and progression of HCC.

Our study showed that the lncRNA prognostic signature was closely related to the TME, which was correlated with the prognosis of HCC patients. It has been reported that several components of the TME, including immune cells, cytokines, chemokines, inhibitory receptors, and ligands, are critical factors in tumorigenesis and progression (43, 44). We further found that almost all immune checkpoints and HLA expression were significantly up-regulated in high-risk patients, which indicated that the poor prognosis might be related to the induction of high immune checkpoint expression. In HCC, it is well known that the increase in the number of PD-1 + and CD8⁺ T cells in tumor tissues and circulation indicates a high recurrence rate and poor prognosis after surgery (45). Research also showed that the upregulation of PD-L1 in HCC cells was induced by a variety of cytokines, especially IFN-γ, which in turn impaired anti-tumor immunity and promoted apoptosis of CD8⁺ T cells (46). Although the numbers of NK cells and CD8⁺ T cells were accumulated during tumor progression, the activity of which will be greatly inhibited after binding with PD-Ls, thus weakening the anti-tumor effect of these cells (47). On the other hand, our result showed that the ratio of pro-tumor signatures and anti-tumor signatures was significantly elevated in low-risk patients, indicating that better prognosis might be attributed to the increased immune/inflammation activity. Previous studies have demonstrated that increased infiltrations of NK, T, and NKT cells are positive prognostic factors for HCC (48–50). In addition, Feng Xu et al. also showed that a higher immunoreactive environment could inhibit tumor growth, progression, invasion, and metastasis (51).

Currently, immunotherapy has been applied to the treatment of various cancers and changed the landscape of cancer care. For example, blocking the combination of PD-1 and PD-L1 can rescue the function of effector T cells to promote their function of killing tumor cells (52). The expression of PD-L1 in tumors is the main factor in determining whether a patient is eligible for PD-1/PD-L1 axis immunotherapy. In clinical practice, however, quite a few PD-L1 positive patients respond poorly to the PD-1/PD-L1 axis treatment, while some patients with negative PD-L1 have a surprisingly response to treatment (53–56). Consistently, our study observed that patients in the high-risk group with upregulation of immune-checkpoints had a weaker response to immunotherapy. On the contrary, patients in the low-risk subgroup with higher immune/inflammation activity were more likely to benefit from immunotherapy, indicating that immune cells infiltration might be a predictor of immunotherapy response.

There were still limitations in our study. Firstly, our results were only based on the TCGA database without external validation. Secondly, no experiment was used to confirm the interaction between the prognostic lncRNAs and m6A modulators in HCC. However, we will conduct further research to provide a better understanding of the above results. In conclusion, we established a novel prognostic model for HCC based on m6A-related lncRNAs. Meanwhile, this study also deepened the understanding of the immune microenvironment status of HCC, and laid a theoretical foundation for the prediction of immunotherapy for HCC patients.

Funding

None

Acknowledgments

Not applicable

Competing interests

The authors declare that they have no competing interests.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The data that generated or analyzed during the present study are openly available in TCGA database (http://www.tcga.org/). The data used to support the findings of this study are available from the corresponding author upon request.

Authors' contributions

Lili Li and Qichun Wei contributed to the study conception and design. Material preparation, data collection and analysis were performed by Rongrong Xie. The first draft of the manuscript was written by Lili Li and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Hartke J, Johnson M, Ghabril M. The diagnosis and treatment of hepatocellular carcinoma. Semin Diagn Pathol. 2017;34(2):153–9.
Torre LA, Bray F, Siegel RL, Ferlay J, Lortet-Tieulent J, Jemal A. Global cancer statistics, 2012. Cancer J Clin. 2015;65(2):87–108.
McGlynn KA, Petrick JL, London WT. Global epidemiology of hepatocellular carcinoma: an emphasis on demographic and regional variability. Clinics in liver disease. 2015;19(2):223–38.
Llovet JM, Zucman-Rossi J, Pikarsky E, Sangro B, Schwartz M, Sherman M, et al. Hepatocellular carcinoma. Nature reviews Disease primers. 2016;2:16018.
Bruix J, Reig M, Sherman M. Evidence-Based Diagnosis, Staging, and Treatment of Patients With Hepatocellular Carcinoma. Gastroenterology. 2016;150(4):835–53.
Vogel A, Cervantes A, Chau I, Daniele B, Llovet JM, Meyer T, et al. Hepatocellular carcinoma: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up. Annals of oncology: official journal of the European Society for Medical Oncology. 2018;29(Suppl 4):iv238–55.
European Association for the Study of the Liver. Electronic address eee, European Association for the Study of the L. EASL Clinical Practice Guidelines: Management of hepatocellular carcinoma. Journal of hepatology. 2018;69(1):182–236.
Bruix J, Gores GJ, Mazzaferro V. Hepatocellular carcinoma: clinical frontiers and perspectives. Gut. 2014;63(5):844–55.
Dai D, Wang H, Zhu L, Jin H, Wang X. N6-methyladenosine links RNA metabolism to cancer progression. Cell death disease. 2018;9(2):124.
Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nature reviews Molecular cell biology. 2017;18(1):31–42.
Yang Y, Hsu PJ, Chen YS, Yang YG. Dynamic transcriptomic m(6)A decoration: writers, erasers, readers and functions in RNA metabolism. Cell research. 2018;28(6):616–24.
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(4):1204–25.
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(2):529–43.
Zhong L, Liao D, Zhang M, Zeng C, Li X, Zhang R, et al. YTHDF2 suppresses cell proliferation and growth via destabilizing the EGFR mRNA in hepatocellular carcinoma. Cancer letters. 2019;442:252–61.
Rong ZX, Li Z, He JJ, Liu LY, Ren XX, Gao J, et al. Downregulation of Fat Mass and Obesity Associated (FTO) Promotes the Progression of Intrahepatic Cholangiocarcinoma. Frontiers in oncology. 2019;9:369.
Cheng X, Li M, Rao X, Zhang W, Li X, Wang L, et al. KIAA1429 regulates the migration and invasion of hepatocellular carcinoma by altering m6A modification of ID2 mRNA. OncoTargets therapy. 2019;12:3421–8.
Yang Z, Li J, Feng G, Gao S, Wang Y, Zhang S, et al. MicroRNA-145 Modulates N(6)-Methyladenosine Levels by Targeting the 3'-Untranslated mRNA Region of the N(6)-Methyladenosine Binding YTH Domain Family 2 Protein. J Biol Chem. 2017;292(9):3614–23.
Chen M, Wei L, Law CT, Tsang FH, Shen J, Cheng CL, et al. RNA N6-methyladenosine methyltransferase-like 3 promotes liver cancer progression through YTHDF2-dependent posttranscriptional silencing of SOCS2. Hepatology. 2018;67(6):2254–70.
Wapinski O, Chang HY. Long noncoding RNAs and human disease. Trends in cell biology. 2011;21(6):354–61.
Shi X, Sun M, Liu H, Yao Y, Song Y. Long non-coding RNAs: a new frontier in the study of human diseases. Cancer letters. 2013;339(2):159–66.
Huang X, Gao Y, Qin J, Lu S. lncRNA MIAT promotes proliferation and invasion of HCC cells via sponging miR-214. American journal of physiology Gastrointestinal liver physiology. 2018;314(5):G559-G65.
Xu F, Zha G, Wu Y, Cai W, Ao J. Overexpressing lncRNA SNHG16 inhibited HCC proliferation and chemoresistance by functionally sponging hsa-miR-93. OncoTargets therapy. 2018;11:8855–63.
Xiong H, Ni Z, He J, Jiang S, Li X, He J, et al. LncRNA HULC triggers autophagy via stabilizing Sirt1 and attenuates the chemosensitivity of HCC cells. Oncogene. 2017;36(25):3528–40.
Yan X, Zhang D, Wu W, Wu S, Qian J, Hao Y, et al. Mesenchymal Stem Cells Promote Hepatocarcinogenesis via lncRNA-MUF Interaction with ANXA2 and miR-34a. Cancer research. 2017;77(23):6704–16.
Li B, Mao R, Liu C, Zhang W, Tang Y, Guo Z. LncRNA FAL1 promotes cell proliferation and migration by acting as a CeRNA of miR-1236 in hepatocellular carcinoma cells. Life sciences. 2018;197:122–9.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.
He L, Li H, Wu A, Peng Y, Shu G, Yin G. Functions of N6-methyladenosine and its role in cancer. Mol Cancer. 2019;18(1):176.
Geng Y, Guan R, Hong W, Huang B, Liu P, Guo X, et al. Identification of m6A-related genes and m6A RNA methylation regulators in pancreatic cancer and their association with survival. Annals of translational medicine. 2020;8(6):387.
Dorn LE, Lasman L, Chen J, Xu X, Hund TJ, Medvedovic M, et al. The N(6)-Methyladenosine mRNA Methylase METTL3 Controls Cardiac Homeostasis and Hypertrophy. Circulation. 2019;139(4):533–45.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nature communications. 2013;4:2612.
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 research. 2012;22(9):1775–89.
Huang GZ, Wu QQ, Zheng ZN, Shao TR, Lv XZ. Identification of Candidate Biomarkers and Analysis of Prognostic Values in Oral Squamous Cell Carcinoma. Frontiers in oncology. 2019;9:1054.
Heagerty PJ, Zheng Y. Survival model predictive accuracy and ROC curves. Biometrics. 2005;61(1):92–105.
Iasonos A, Schrag D, Raj GV, Panageas KS. How to build and interpret a nomogram for cancer prognosis. Journal of clinical oncology: official journal of the American Society of Clinical Oncology. 2008;26(8):1364–70.
Junttila MR, de Sauvage FJ. Influence of tumour micro-environment heterogeneity on therapeutic response. Nature. 2013;501(7467):346–54.
Ishii G, Ochiai A, Neri S. Phenotypic and functional heterogeneity of cancer-associated fibroblast within the tumor microenvironment. Adv Drug Deliv Rev. 2016;99(Pt B):186–96.
Olson B, Li Y, Lin Y, Liu ET, Patnaik A. Mouse Models for Cancer Immunotherapy Research. Cancer discovery. 2018;8(11):1358–65.
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.
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(1):186.
Zhao JS, Li JH, Yang JH, Chen SS. Construction and analysis of prognostic model for hepatocellular carcinoma based on autophagy-related long non-coding RNAs. Chinese Journal of General Surgery. 2020;29(7):839–84.
Tu W, Yang Y, Song Y, Zhu W. Hepatitis B virus x protein accelerated the proliferation of hepatocellular carcinoma cell through lncRNA SNHG20/PTEN pathway. J BioChem. 2019;165(5):423–31.
Wu G, Ju X, Wang Y, Li Z, Gan X. Up-regulation of SNHG6 activates SERPINH1 expression by competitive binding to miR-139-5p to promote hepatocellular carcinoma progression. Cell cycle. 2019;18(16):1849–67.
Aravalli RN. Role of innate immunity in the development of hepatocellular carcinoma. World journal of gastroenterology. 2013;19(43):7500–14.
Matsui M, Machida S, Itani-Yohda T, Akatsuka T. Downregulation of the proteasome subunits, transporter, and antigen presentation in hepatocellular carcinoma, and their restoration by interferon-gamma. Journal of gastroenterology hepatology. 2002;17(8):897–907.
Liao H, Chen W, Dai Y, Richardson JJ, Guo J, Yuan K, et al. Expression of Programmed Cell Death-Ligands in Hepatocellular Carcinoma: Correlation With Immune Microenvironment and Survival Outcomes. Frontiers in oncology. 2019;9:883.
Shi F, Shi M, Zeng Z, Qi RZ, Liu ZW, Zhang JY, et al. PD-1 and PD-L1 upregulation promotes CD8(+) T-cell apoptosis and postoperative recurrence in hepatocellular carcinoma patients. International journal of cancer. 2011;128(4):887–96.
Lee YH, Martin-Orozco N, Zheng P, Li J, Zhang P, Tan H, et al. Inhibition of the B7-H3 immune checkpoint limits tumor growth by enhancing cytotoxic lymphocyte function. Cell research. 2017;27(8):1034–45.
Pang YL, Zhang HG, Peng JR, Pang XW, Yu S, Xing Q, et al. The immunosuppressive tumor microenvironment in hepatocellular carcinoma. Cancer immunology immunotherapy: CII. 2009;58(6):877–86.
Chew V, Chen J, Lee D, Loh E, Lee J, Lim KH, et al. Chemokine-driven lymphocyte infiltration: an early intratumoural event determining long-term survival in resectable hepatocellular carcinoma. Gut. 2012;61(3):427–38.
Wada Y, Nakashima O, Kutami R, Yamamoto O, Kojiro M. Clinicopathological study on hepatocellular carcinoma with lymphocytic infiltration. Hepatology. 1998;27(2):407–14.
Xu F, Chen JX, Yang XB, Hong XB, Li ZX, Lin L, et al. Analysis of Lung Adenocarcinoma Subtypes Based on Immune Signatures Identifies Clinical Implications for Cancer Therapy. Molecular therapy oncolytics. 2020;17:241–9.
Nadal E, Massuti B, Domine M, Garcia-Campelo R, Cobo M, Felip E. Immunotherapy with checkpoint inhibitors in non-small cell lung cancer: insights from long-term survivors. Cancer immunology immunotherapy: CII. 2019;68(3):341–52.
Garon EB, Rizvi NA, Hui R, Leighl N, Balmanoukian AS, Eder JP, et al. Pembrolizumab for the treatment of non-small-cell lung cancer. N Engl J Med. 2015;372(21):2018–28.
Herbst RS, Soria JC, Kowanetz M, Fine GD, Hamid O, Gordon MS, et al. Predictive correlates of response to the anti-PD-L1 antibody MPDL3280A in cancer patients. Nature. 2014;515(7528):563–7.
Robert C, Schachter J, Long GV, Arance A, Grob JJ, Mortier L, et al. Pembrolizumab versus Ipilimumab in Advanced Melanoma. N Engl J Med. 2015;372(26):2521–32.
Tumeh PC, Harview CL, Yearley JH, Shintaku IP, Taylor EJ, Robert L, et al. PD-1 blockade induces responses by inhibiting adaptive immune resistance. Nature. 2014;515(7528):568–71.

FigureS1.tif
Supplementary File 1 legend: The 45 prognostic-related lncRNAs screened out through K–M survival analyses. Supplementary Figure 1 legend: Sample quality control chart.
FigureS2.tif
Supplementary Figure 2 legend: Differentially expressed genes (DEGs) in the two risk subgroups. (A). Volcano plots of DEGs. The red dots in the plot represent upregulated genes and blue dots represent downregulated genes with statistical significance. Grey dots represent no differentially expressed genes. (B). Heatmap of significantly DEGs. Red represents higher expression while green represents lower expression.
FigureS3.tif
Supplementary Figure 3 legend: PCA analysis in the two risk subgroups. (A). PCA analysis in the training group. (B). PCA analysis in the internal validation set the testing group.
FigureS4.tif
Supplementary Figure 4 legend: Kaplan–Meier curves of clinicopathological characteristics and risk group stratification within the training set for OS. Clinicopathological features included age (A), BMI (B), gender (C), T stage (D), N stage (E), M stage (F), pathologic stage (G), prior malignancy (H), treatment or therapy (I), and treatment type (J).
FigureS5.tif
Supplementary Figure 5 legend: Relationship between the risk scores and clinicopathological factors in the testing cohort. (A). The heatmap showed the distribution of clinicopathological factors and three m6A methyltransferase‑related lncRNAs between the high- and low- risk groups. (B-K). Distribution of risk scores stratified by T stage (B), pathologic stage (C), age (D), BMI (E), gender (F), M stage (G), N stage (H), prior malignancy (I), treatment or therapy (J) and treatment type (K). *P < 0.05, **P < 0.01.
FigureS6.tif
Supplementary Figure 6 legend: Kaplan–Meier curves of clinicopathological characteristics and risk group stratification within the testing set for OS. Clinicopathological features included age (A), BMI (B), gender (C), T stage (D), N stage (E), M stage (F), pathologic stage (G), prior malignancy (H), treatment or therapy (I), and treatment type (J).
SupplementaryTableS1.xlsx
SupplementaryTableS2.xlsx
SupplementaryTableS3.xlsx
SupplementaryTableS4.xlsx
SupplementaryTableS5.xlsx
SupplementaryTableS6.xlsx
SupplementaryTableS7.xlsx
SupplementaryTableS8.xlsx
SupplementaryTableS10.xlsx
SupplementaryTableS11.xlsx
SupplementaryTableS12.xlsx
SupplementaryTableS13.xlsx
SupplementaryTableS14.xlsx
LIHCCASC90.0148SurvivalPlot.pdf
LIHCCRNDE0.1004SurvivalPlot.pdf
LIHCCYTOR1e04SurvivalPlot.pdf
LIHCDUXAP80.0147SurvivalPlot.pdf
LIHCELF3.AS10.0517SurvivalPlot.pdf
LIHCFAM99A0.0047SurvivalPlot.pdf
LIHCFAM99B0.037SurvivalPlot.pdf
LIHCFOXD2.AS10.0303SurvivalPlot.pdf
LIHCGAS50.1617SurvivalPlot.pdf
LIHCGSEC6e04SurvivalPlot.pdf
LIHCHNF1A.AS10.1288SurvivalPlot.pdf
LIHCHNF4A.AS10.0779SurvivalPlot.pdf
LIHCLINC005110.002SurvivalPlot.pdf
LIHCLINC010930.0523SurvivalPlot.pdf
LIHCLINC013700.1256SurvivalPlot.pdf
LIHCLINC014190.1944SurvivalPlot.pdf
LIHCLINC015544e04SurvivalPlot.pdf
LIHCLINC017030.0239SurvivalPlot.pdf
LIHCLINC020370.0481SurvivalPlot.pdf
LIHCCH17.340M24.30.0269SurvivalPlot.pdf
LIHCLINC020550.1514SurvivalPlot.pdf
LIHCLINC021630.1866SurvivalPlot.pdf
LIHCLINC023620.1795SurvivalPlot.pdf
LIHCLUCAT10.0279SurvivalPlot.pdf
LIHCMAPKAPK5.AS10.0016SurvivalPlot.pdf
LIHCMIR4435.2HG0.0031SurvivalPlot.pdf
LIHCMYLK.AS10.031SurvivalPlot.pdf
LIHCNALT10.0859SurvivalPlot.pdf
LIHCPCAT60.0629SurvivalPlot.pdf
LIHCPITPNA.AS10.0554SurvivalPlot.pdf
LIHCPRR34.AS10.1776SurvivalPlot.pdf
LIHCPVT10.1909SurvivalPlot.pdf
LIHCPXN.AS10.0449SurvivalPlot.pdf
LIHCRALY.AS10.063SurvivalPlot.pdf
LIHCSNHG170.1073SurvivalPlot.pdf
LIHCSNHG70.0498SurvivalPlot.pdf
LIHCTSPEAR.AS20.1092SurvivalPlot.pdf
LIHCUBR5.AS10.0294SurvivalPlot.pdf
LIHCZFAS10.0361SurvivalPlot.pdf
LIHCZFPM2.AS10.0013SurvivalPlot.pdf
LIHCLINC023770.1152SurvivalPlot.pdf

Download PDF

Version 1

posted

You are reading this latest preprint version

Significance of m 6A Methyltransferase‑Related lncRNAs in the Prognosis and Immunotherapy of Hepatocellular Carcinoma

Status:

Version 1

Abstract

Figures

Introduction

Materials And Methods

Acquisition of gene expression and clinical data

Analysis of differentially expressed PCGs and lncRNAs

Identification of m6A methyltransferase‑related lncRNAs

Establishment and validation of the prognostic gene signature

Construction of a predictive nomogram

Gene set enrichment analyses

ESTIMATE algorithm

Statistical analysis

Results

Identification of DE-PCGs and -lncRNAs in HCC

Identification of prognostic m6A methyltransferase‑related lncRNAs

Establishment of m6A methyltransferase-related lncRNAs signature

Prognostic risk score displayed strong correlations with clinicopathological features in HCC patients

The risk score was an independent prognostic factor in HCC

Functional enrichment analyses of the DEGs between the high- and low-risk groups

The landscape of TME immune cells infiltration in the two risk subgroups

Low-risk patients were likely to be more sensitive to immunotherapy

Discussion

Declarations

References

Supplementary Files

Status:

Version 1