N6-methyladenosine (m6A)-mediated messenger RNA signatures and the tumor immune microenvironment can predict the prognosis of hepatocellular carcinoma.

BACKGROUND: N6-methyladenosine (m6A)-mediated ribonucleic acid (RNA) methylation is considered to be the most significant and abundant epigenetic modification in eukaryotic cells, and plays an essential role in the carcinogenesis and molecular pathogenesis of hepatocellular carcinoma (HCC). However, the relationship between m6A regulation and immune cell infiltration of the tumor immune microenvironment (TIME) has not yet been clarified. We aimed to investigate the roles of m6A RNA gene regulators in HCC immune regulation and prognosis.
METHODS: The Cancer Genome Atlas (TCGA) database was used, and unsupervised clustering of 21 m6A regulators was performed based on differential gene expression. Gene Set Variation Analysis (GSVA), single-sample Gene Set Enrichment Analysis (ssGSEA), the empirical Bayes method, and m6A scores were used in our analyses.
RESULTS: Of 433 samples, 101 (23.22%) had m6A regulatory factor mutations. From these, we identified three m6A subtypes, which correlated with different TIME phenotypes: immune rejection, immune infiltration, and immune deficiency. Tumors with low methyltransferase-like 3 (METTL3) expression had increased infiltration of dendritic cells (DCs) in the TIME. Reduced METTL3 expression also led to an overall increase in expression of major histocompatibility complex (MHC) molecules, costimulatory molecules, and adhesion molecules. The m6A subtypes were scored and analyzed for correlations. Patients with epithelial-mesenchymal transition (EMT) subtypes had lower m6A scores than the other three molecular subtypes. Survival analysis found that patients with low m6A scores had better overall survival [hazard ratio (HR) 1.6 (1.1-2.3)] and a 1.16 times better 5-year survival rate than patients with high m6A scores (56% vs. 48%).
CONCLUSIONS: Our results demonstrated that three different m6A modification subtypes contribute to immune regulation in HCC and have potential as novel prognostic indicators and immune therapeutic targets. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Hepatocellular carcinoma (HCC) remains a high incidence malignancy with high mortality rate and poor prognosis worldwide. It has been reported that there are approximately 850,000 new HCC patients annually, with about 500,000 deaths due to HCC-associated diseases (1,2). Hepatitis B and C viral infections, excess alcohol consumption, aflatoxin B1 exposure, obesity, diabetes, and metabolic diseases are major risk factors for the development of HCC (3). Owing to the difficulty in making an early diagnosis, HCC patients are typically diagnosed at an advanced stage of the disease, and thus, treatments are generally associated with a low survival rate, high risk of recurrence, malignant metastasis, as well as a high risk of selection and spreading of drug resistance (4). Recently, immune therapies, such as immune checkpoint inhibition, have been used in advanced HCC to modify the carcinogenic process and to augment adaptive immunity (5). It has been shown that the tumor immune microenvironment (TIME), which is composed of both pro- and anti-tumor immune cells, can be reprogrammed by tumor-derived factors, which are involved in immune evasion and tumor progression (6,7). Accumulating evidence supports a role for the TIME in determining HCC progression, recurrence, metastasis, and poor outcome (8,9). Tumor-infiltrating immune cells (TIICs) are important components of the TIME and have been shown to be significant predictors of HCC survival (10). However, correlations between TIME parameters and HCC genesis and development are still poorly understood (10).

Currently, the most well studied messenger ribonucleic acid (mRNA) modification is N6-methyladenosine (m6A). It is the most pervasive internal mRNA modification (10), and correlates with tumor progression (11). Previous studies have shown that m6A demethylation of the fat mass and obesity-associated protein (FTO) is up-regulated and promotes tumor proliferation and metastasis in breast cancer (12,13). However, the role of the m6A modification in HCC tumorigenesis remains unclear.

We hypothesized that m6A may regulate HCC tumorigenesis and progression by regulating the immune system and thus may be a therapeutic drug target for HCC. To evaluate this hypothesis, data from The Cancer Genome Atlas (TCGA) were used as an effective cohort, and data from the Gene Expression Omnibus (GEO) database were used as a validation cohort. Genes associated with the immune microenvironment were downloaded from the Immunology Analysis Portal database (https//www.immport.org/shared/genelists), analyzed, and immune phenotypes were classified by immune scores. We validated the different m6A modification patterns and explored the relationship between immune infiltration and m6A regulation. We also identified three different immune phenotypes: the immune rejection phenotype, the immune inflammation phenotype, and the immune desert phenotype, and found that formation of the tumor microenvironment plays a role in the phenotype. We combined m6A modification patterns and immune infiltration phenotypes to establish a system to estimate HCC prognosis.

We present the following article in accordance with the REMARK reporting checklist (available at http://dx.doi.org/10.21037/atm-20-7396).

Methods

Data collection

The work-flow of this study is illustrated in Figure S1A. All related information was downloaded from TCGA data portal (http://cancergenome.nih.gov/). And the study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). Data cohorts with missing information were removed. Data collection, refinement, and model statistics are summarized in Table S1. RNA sequence data from the National Cancer Institute’s (NCIs) Genomic Data Commons (GDC, https://gdc.cancer.gov/about-data/publications/pancanatlas), Fragments Per Kilobase per total Million mapped reads (FPKM), and TCGA biolinks were also used (14,15). Data were analyzed using the statistical software package R, version 3.6.1. Raw data were processed into R/Bioconductor.

Unsupervised gene cluster analysis

We used clustering, a resampling-based clustering algorithm, to analyze the m6A gene profiles of HCCs. A total of 21 m6A regulators, including eight writers, two erasers, and 11 readers were analyzed. We performed unsupervised cluster analysis with a consistent clustering algorithm and the Consensus Cluster Plus package as described previously (16).

Gene Set Variation Analysis (GSVA) and functional annotation

For an in-depth analysis of the biological processes that associate with the different m6A modification patterns, the “GSVA” R package for GSVA enrichment analysis was utilized. The gene set “c2.cp.kegg. v6.2. symbols” was downloaded from the MSigDB database and used to run GSVA analysis. For functional annotation of genes, we utilized cluster analysis and a false discovery rate <0.05.

Estimation of the immune cells infiltrating the TIME

A gene set of a variety of human immune cell subtypes was obtained from Charoentong et al. (17). Single-sample Gene Set Enrichment Analysis (ssGSEA) calculates separate enrichment scores for each sample and gene pair. This analysis showed the relative abundance of each cell type that infiltrated the TIME for each sample.

Classification of m6A phenotypes based on differential gene expression

To validate and classify m6A-associated genes, differentially expressed gene (DEG) analysis (fold-changes) and statistical computations of the two conditions were conducted using the EBSeq R package (https://bioconductor.org/packages/release/bioc/html/EBSeq.html).

Construction of m6A gene markers and m6A scoring

To validate m6A regulator modifications in HCC samples, we established a scoring system to estimate m6A gene characteristics. To establish m6A gene markers, DEGs of the different m6A clusters were normalized to extract overlapping genes, the unsupervised clustering method was used to analyze the overlapping data, and patients were divided into multiple groups for in-depth analysis. The number and stability of gene clusters were defined by the consistent clustering algorithm. We then used a univariate Cox regression model to analyze the association of each gene cluster with prognosis. We selected principal component 1 and principal component 2 as signature scores. Subsequently, we defined the m6a score using a method like that used in analyzing gene-gene interactions (GGIs): m6A Score= Σ (PC1i + PC2i).

Analysis of m6A modifications and their immune-related signatures in HCC

To investigate associations between m6A modifications and the immune microenvironment, we downloaded immune-related genes from the ImmPort database (http://www.immport.org), and identified prognostic gene signatures based on genes with independent prognostic values as determined previously (18). Our m6A scoring system revealed that m6A gene characteristics are associated with particular biological pathways.

Statistical analysis

DEGs were applied in this study with the threshold of absolute log2-fold-change >1 and an adjusted P value <0.05 using the R package “limma”. The overall survival rate was evaluated using the log-rank test with Kaplan-Meier estimation. Differences between the groups were evaluated by analysis of variance (ANOVA). All P values were two-sided, and P values <0.05 were considered significant. Comparative studies were analyzed using the Kruskal-Wallis test and Wilcoxon test. *P<0.05, **P<0.01, ***P<0.001. All statistics were calculated using R language.

Results

Expression of m6A regulators in HCC

To better understand the search strategy and selection of datasets for comprehensively understand m6A regulation in HCC progression, we constituted a flow chart of this study (Figure S1A). In total, 371 patients were involved in this study and three clusters were identified. The transcription levels of 21 m6A regulators, which included 11 readers, eight writers, and two erasers, as well as the copy number variations (CNVs) and somatic mutations of these regulators, were described in a previous study (19). The chromosomal locations of the amplification mutations were distributed randomly and unequally across the chromosomes (). Of 64 HCC samples, 44 contained mutations related to m6A modifications (68.75% mutation rate). The most prevalent mutation were found in the KIAA1429 (8%) and LRPPRC (8%) genes, followed by HNRNPC (6%) and YTHDC2 (6%) (). The 21 m6A regulators exhibited widespread gene gain and loss mutations, KIAA1429 and YTHDF3 contained the most gain mutation type, while ALKBH5 and ZC3H13 had the most frequent loss mutation type (). We also explored the 21 regulators in tumor tissue compared with normal tissue, and easily identified the two subgroups; the red spot represented normal tissue and the blue spot represented tumor tissue ().

Figure 1

Overview of m6A gene locus and gene information. (A) The m6A regulators mutation location. (B) Waterfall plot of m6A regulators mutation gene and mutation type. (C) Frequency plots of copy number gains (in green) and losses (in red) defined in all m6A regulators. (D) Regulators could be divided into two subgroups (red and green) based on the gene expression level. (E) Comparison of gene expression level of 21 regulators between the normal and tumor tissue cohort. **P<0.01, ***P<0.001. m6A, N6-methyladenosine.

To further investigate the overall survival between the regulators, HCC patients were divided into mutated subgroup and wide type subgroup based on whether they were combined with m6A mutation information. Obviously, the mutated subgroups had a poor prognosis in YTHDF1 [P=0.018, hazard ratio (HR) =4.7], LRPPRC (P=0.003, HR =6.44), and FTO (P=4.12e7, HR =11.3) (Figure S1B). We also performed a meta-analysis, and the results revealed that the regulators have prognostic potential in HCC cohorts (Figure S1C). The mean expression levels of m6A regulators in HCC specimens and normal tissues are shown in ; the difference was significant and m6A regulators were highly expressed in both tumor tissue and normal tissue. In this study, we mapped the gene somatic mutation information of 21 m6A regulators between tumor tissue and normal tissue, and confirmed that m6A modification plays an important role in HCC tumorigenesis and progression.

m6A regulators have close relationship with cellular crucial functional pathways regulation

In order to comprehensively explain the crosslink among the 21 m6A regulators, we divided 21 m6A regulators into three clusters as illustrated in . Through GSVA of functional genes, we determined the biological roles of the m6A regulators. The expression of genetic information of the m6A regulators was controlled through intricate regulatory networks and a complex regulatory relationship. The overall survival curve implied a significant survival advantage for cluster 2 (*P<0.0001, ).

Figure 2

Crosslinks between 21 m6A regulators and associated biological activities. (A) Landscape and inner crosslink between 21 m6A regulators. (B) Kaplan-Meier curve of survival probability between three clusters. (C) Cluster enrichment analysis of biological activities between cluster 1 and cluster 2. (D) Cluster enrichment analysis of biological activities between cluster 2 and cluster 3. m6A, N6-methyladenosine.

In the present study, we explored the relationship between different modification patterns and biological functional activities. The results indicated that three clusters exhibit significant pathway enrichment. The cluster 1-associated pathways gather on translation functions, such as olfactory transduction, neuroactive ligand receptor interaction, arachidonic acid metabolism complements and coagulation cascades, folate biosynthesis, arginine and proline metabolism, metabolism of xenobiotics by cytochrome, glutathione metabolism and Parkinson’s disease cardiac muscle contraction, amyotrophic lateral sclerosis, prion diseases, and autoimmune diseases. Cluster 2 functional pathways were significantly enriched in the ERBB signaling pathway, cell cycle signaling pathway, cell adhesion, as well as other pathways related to cell matrix adhesion and carcinogenesis, such as colorectal cancer, chronic myeloid leukemia, adherent’s junction, the ERBB signaling pathway, endometrial cancer, the cell cycle, oocyte meiosis, RNA degeneration, ubiquitin mediated proteolysis, basal transcription factors, non-homologous end joining, and thyroid cancer (). The third cluster-associated functional pathways act in metabolism pathways (). Our results implied that m6A regulators participated in many important cellular activities, including cellular metabolism regulation, nuclear factor regulation, signal transitional regulation, and carcinogenesis.

m6A regulators have a close relationship with the infiltration of immune cells

Next, we estimated the enrichment score of immune cells among the three clusters and identified the differences in the infiltration of immune cells among clusters 1–3. The box plot exhibited significant differences in most immune cells, including activated B cells, activated cluster of differentiation (CD)4+T cells, activated CD8+T cells, ᵞᵟT cells, regulatory T cells, type 1 helper cells, type 17 helper cells, activated dendritic cells (DCs), CD56 bright natural killer (NK) cells, CD56 dim NK cells, eosinophils, macrophages, and mast cells ().

Figure 3

Immune cell infiltration and immune associated classification among m6A regulators. (A) Comparison between different immune cells infiltration in three clusters. (B) Comparison of cellular biological activities with enriched regulation pathways among three clusters. (C) Comparison of the immune score between the high and low METTL3 expression subgroups. The immune score exhibited a weak difference between high and low METTL3 expression levels. (D) Correlation between m6A regulators and biological pathways in the HCC cohort. The Spearman analysis was applied in this study. Negative correlation was marked with blue and positive correlation with orange. *P<0.05, **P<0.01, ***P<0.001. m6A, N6-methyladenosine.

To further investigate the biological behaviors among the three m6A modification clusters, pathway enrichment analysis was applied. Our results indicated that the transforming growth factor beta (TGF-β) pathway, angiogenesis, and the tumor EMT process were remarkably different among the clusters (). Our research suggested that the three types of regulators were involved in crucial functional regulations with different pathways. According to the GSVA analysis, the 21 m6A regulators were classified as either an immune rejected phenotype, an immune inflamed phenotype, or an immune desert phenotype based on the infiltration of immune cells ().

We then focused on methyltransferase like 3 (METTL3), and the immune score of different METTL3 expression subgroups showed a weak difference (). Furthermore, Spearman correlation analysis was performed for TIME infiltration and the m6A regulators (). METTL3 expression correlated with many of the TIME infiltrating immune cells and the immune score. Our study implied that the infiltration of CD4+ T cells, type 1 T helper cells, type 2 T helpers, CD56 bright NK cells, eosinophils, macrophages, mast cells, monocytes, and neutrophils exhibited a significant difference between high and low METTL3 expression (Figure S2A), indicating that low METTL3 expression resulted in increased TIME-associated immune cell infiltration.

Subsequently, we estimated the gene expression of immune checkpoint molecules and receptors, and discovered that CD80, CD86, intercellular adhesion molecule (ICAM) 1, ICAM3, and programmed cell death ligand 1 (PDL1) gene expression showed differences across different METTL3 expression levels (Figure S2B). We then evaluated the enrichment score of some crucial cell activity pathways between two clusters, and the box plot exhibited significant differences in the TGF-β pathway, T cell activation, angiogenesis, and the tumor EMT process (Figure S2C). These results indicated that the m6A regulators play key roles in immune cell infiltration and characteristic TIME formation.

The transcriptome and clinical features of m6A regulators

To further identify the cellular biological behavior and in-depth mechanisms of the m6A regulators, we performed unsupervised cluster analysis on m6A phenotype-related genes to divide patients into one of four different clusters based on the m6A modified genomic phenotypes (). Based on the m6A regulator classifications, we confirmed that cluster 1, with the abundant immune cell infiltration subgroup, exhibited a better prognosis and survival probability (). Among the three immune cell infiltration-associated subgroups, expression of regulatory factors of 21 regulators differed significantly ().

Figure 4

Clinical information in different m6A regulators. (A) m6A regulators were divided into four clusters based on gene expression. Related clinical information was involved in this cluster enrichment analysis. (B) Comparison of overall survival rate between the three different subgroups. (C) Comparison of m6A regulators gene expression levels among the three clusters. (D) Comparison of high and low immune score clusters in cellular activities and biological related pathways. m6A, N6-methyladenosine. *P<0.05, **P<0.01, ***P<0.001.

Further analysis of the matrix-related pathways revealed that the gene expression level of the m6A regulators suggested significantly different expression levels of immune related functional annotation and activation of matrix pathways, including CD8+T cells effector, angiogenesis, the cell cycle, cell cycle regulator, deoxyribonucleic acid (DNA) damage repair, DNA replication, EMT1, Fanconi anemia, homologous recombination, mismatch repair, nucleotide excision repair, and Pan-F-TBRS (). Here, we confirmed that the 21 regulators were classified into three gene clusters based on gene expression and cluster analysis. The box plots illustrated that the 21 regulators gene expression clusters were consistent with the functional classification.

The m6A regulators gene expression classification with immune infiltration

To better annotate the functions of m6A-related genes with immune cells, our analysis results implied that the infiltration of activated CD4+ T cells, regulatory T cells, type 17 T helper cells, CD56 bright NK cells, CD56 dim NK cells, and eosinophils exhibited significant differences between the three clusters (). To explore the roles of m6A-related phenotypes in the modulation of the TIME, we analyzed the expression of chemokines and cytokines in the three m6A gene clusters. Interferon gamma (IFN-γ), CD8A, T-Box transcription factor 2 (TBX2), and tumor necrosis factor (TNF) showed significant differences among the three clusters. Other cytokines exhibited no or weak differences between these clusters (). In , we estimated the gene expression of immune checkpoint molecules and receptors and found that HAVCR2, CD80, PDCD1, TIGIT, CD86, COL4A1, ZEB, TGFB2, and TWIST1 gene expression showed remarkable differences between the three clusters, which suggested that these immune checkpoints had a significant relationship with the m6A regulator classification.

Figure 5

Comparison of different kinds of immune infiltration and prognosis analysis. (A) Comparison of enrichment score with different immune cells infiltration in the three clusters. (B) Comparison of enrichment score with different immune-related cytokines expression among the three clusters. (C) Expression level comparison of different immune checkpoint targets and associated receptors in the three clusters. (D) Expression comparison of m6A regulator-related kinases in the three clusters. (E) Overall survival analysis between the high and low immune score subgroups. (F) AUC of m6A regulators’ signature validation of the survival value of the risk score. m6A, N6-methyladenosine. *P<0.05, **P<0.01, ***P<0.001.

We further investigated the prognostic value of the m6A scores for HCC. HCC patients were divided into high and low m6A score groups, and the survival of these groups was compared. The overall survival for the low score m6A group was better than the high score group [P=0.019, HR 1.6 (1.1–2.3)]. The 5-year survival rate was 56% for the high m6A score group and 48% for the low score group (). The area under the curve (AUC) in patients treated with immunotherapy was 0.768 (). Our data demonstrated that the m6A evaluation scores based on modification patterns were correlated with immune or matrix activation, which affects tumor immune infiltrates and patient prognosis.

Discussion

Globally, HCC is still a lethal malignancy with a poor prognosis. The tumorigenesis mechanism of HCC must be elucidated to improve prognosis, identify biomarkers, and develop effective therapeutic strategies. The liver is the largest immune-related organ, and the immune system plays a definitive role in oncogenesis (20,21). Recent studies have indicated that aberrant mRNA modifications contribute to HCC prognosis and overall survival. Currently, m6A is the most well studied post-transcriptional mRNA modification, and evidence suggests that this modification epigenetically affects mRNA processing, translation, and stability (22). Dysregulation of the m6A modification has been shown to be associated with the initiation and progression of various cancers. Thus, exploring the relationship between m6A regulation and the TIME could lead to the design of immune-based targeted therapies for HCC.

METTL3 is a crucial oncogene in several cancers and has been shown to be up-regulated in gastric cancer (23), prostate carcinoma (24), and osteosarcoma (25). Also, methyltransferase like 14 (METTL14) has been shown to play a crucial role in the RNA epigenetics of colorectal cancer, and thus, has prognostic value (26). It is also been shown that YTHDF2 suppresses cancer cell migration via m6A-associated modification (26). In the present study, we identified three different m6A modification patterns and their relationship with the infiltration of immune cells, such as NK cells, macrophages, eosinophils, mast cells, myeloid-derived suppressor cells (MDSCs), and plasma cell-like DCs, into the HCC microenvironment. Based on the infiltration of immune cells and the immune score, we identified three different immune-associated phenotypes of the HCC tumor environment: immune excluded, immune inflamed, and immune desert. The immune excluded phenotype is characterized by the presence of T cells at the periphery of cancer nests. In the immune excluded phenotype of bladder cancer, there was no PD-L1 expression and only a few tumor-infiltrating lymphocytes (TILs) were present (27). Mariathasan et al. found that combination therapy with TGF-blocking and anti-PD-L1 antibodies decreased TGF signaling in stromal cells, facilitated concentration of T cells in the centers of tumors, and provoked vigorous anti-tumor immunity and tumor regression (28). Zhao et al. demonstrated that the immune inflamed phenotype of triple-negative breast cancers is characterized by the infiltration of CD8+ T cells into the tumor parenchyma, indicating that these patients are more likely to benefit from immune checkpoint inhibitor therapy (29). Job et al. reported that the inflamed subtype is characterized by massive T lymphocyte infiltration as well as activation and upregulation of inflammatory and immune checkpoint pathways. This phenotype is associated with a prolonged survival prognosis (30). Thus, patients with the immune inflamed phenotype may benefit from immunotherapy. In this study, we found that the immune desert phenotype is characterized by CD8 cell deficiency and is associated with a poor prognosis.

The TIME can predict patient survival in most cancer types, and m6A modification has been shown to affect cell infiltration in the TIME and the epigenetic regulation of the immune response in gastric cancer. M6A has also been shown to regulate the innate immune response by targeting type I interferons (18,31). In the innate immune system, neoantigen-specific immunity is regulated by mRNA methylation via the m6A-binding protein YTHDF1 (32). In this study, we explored m6A regulatory factor-mediated TIME cell infiltration and found more DCs, including activated DCs, immature DCs, and plasma cell-like DCs in tumors with low METTL3 expression, indicating that METTL3-mediated m6A methylation modification activates DCs in the TIME, thereby enhancing the anti-tumor immune response.

Conclusions

Our findings demonstrated that numerous immune-related genes regulated by m6A modification contribute to the immune phenotype. Moreover, we validated a novel model of m6A modification and immune-related gene regulation to predict HCC prognosis. The immune score and immune classification systems may lead to novel therapeutic targets for HCC.

The article’s supplementary files as