RNA N6-Methyladenosine Patterns in Hepatocellular Carcinoma Reveal a Distinct Immune Infiltration Landscape and Clinical Significance.

BACKGROUND RNA N6-methyladenosine (m6A) methylation, the most abundant and prominent form of epigenetic modification, is involved in hepatocellular carcinoma (HCC) initiation and progression. However, the role of m6A methylation in HCC tumor microenvironment (TME) formation is unexplored. This study aimed to reveal the TME features of HCC patients with distinct m⁶A expression patterns and establish a prognostic model based on m⁶A signatures for HCC cohorts. MATERIAL AND METHODS We classified the m⁶A methylation patterns in 365 HCC samples based on 21 m6A modulators using a consensus clustering algorithm. Single-sample gene set enrichment analysis algorithm was used to quantify the abundance of immune cell infiltration. Gene set variation analysis revealed the biological characteristics between the m⁶A modification patterns. The m6A-based prognostic model was constructed using a training set with least absolute shrinkage and selection operator regression and validated in internal and external datasets. RESULTS Two distinct m⁶A modification patterns exhibiting different TME immune-infiltrating characteristics, heterogeneity, and prognostic variations were identified in the HCC cohort. After depicting the immune landscape of TME in HCC, we found patients with high LRPPRC m⁶A modulator expression had depletion of T cells, cytotoxic cells, dendritic cells, and cytolytic activity response. A high m⁶A score, characterized by suppression of immunity, indicated an immune-excluded TME phenotype, with poor survival. A nomogram was developed to facilitate HCC clinical decision making. CONCLUSIONS Our results highlight the nonnegligible role of m6A methylation in TME formation and reveal a potential clinical application of the m⁶A-associated prognostic model for patients with HCC.

Background

Hepatocellular carcinoma (HCC) is the sixth most frequently diagnosed cancer and the fourth leading cause of cancer-related death worldwide [1]. Early-stage HCC requires curative surgical resection or liver transplantation. However, most patients are at an advanced stage when diagnosed and do not have surgical opportunities [2]. First-line treatment multi-target kinase inhibitors, including sorafenib, can only extend overall survival (OS) of patients with advanced HCC by 3 months, and tumor progression occurs in some patients because of drug resistance [3,4]. Even for those receiving surgery, there is a high incidence of tumor recurrence, and most patients with recurrence die within a year [5]. Therefore, elucidating the molecular mechanisms underlying HCC occurrence and development is necessary and urgent.

N6-methyladenosine (m6A) is the most abundant and prominent internal modification of messenger RNAs (mRNAs) in eukaryotic cells [6]. Similar to other epigenetic alterations, including DNA methylation and histone modification, m6A modification is a dynamic and reversible biological process, which is regulated by 3 types of enzymes [7]. Methyltransferases, or “writers”, catalyze the transfer of methyl groups onto the sixth position of adenosines; demethylases, or “erasers”, remove methyl groups; and RNA binding proteins, or “readers”, recognize and bind to specific m6A sites to regulate RNA metabolism. In 2012, the transcriptome-wide m6A modification landscape was identified for the first time [8,9]. High-throughput screening revealed that most m6A sites are in termination codons, 3′-untranslated regions, and long internal exons. m6A plays an important role in maintaining cellular biological functions by regulating mRNA metabolic processes, including alternative splicing, stability, translation, and localization. The in-depth excavation of these modulators would facilitate exploration of the mechanism and role of m6A modification in post-transcriptional regulation. Mounting evidence suggests that aberrant m6A modification is involved in multiple pathological processes, including dysregulated cell proliferation and death, abnormal immune response, and malignant progression of various cancers [10,11].

Immunotherapies represented by immune checkpoint inhibitors (ICIs) such as CTLA4 and PD-1/L1 have demonstrated therapeutic efficacy in a variety of cancers [12,13]. In patients with advanced HCC, the positive response rate to anti-PDL1 blockade is lower than 20%, which might result from HCC tumor heterogeneity and is far from satisfying clinical needs [14,15]. Previous studies indicated that the tumor microenvironment (TME) in which tumor cells proliferate and evade immune surveillance plays a crucial role in tumor initiation and progression [16,17]. The TME includes the extracellular matrix, neovascular and stromal cells, such as cancer-associated fibroblasts and macrophages, and recruited immune cells such as regulatory T cells (Tregs) and bone marrow-derived cells. Cancer cells, together with other TME components, reciprocally regulate the biological behaviors of cancer, including apoptosis resistance, proliferation, neovascularization, immune evasion, and tumor response to immunotherapies [18]. Therefore, comprehensively characterizing TME cell infiltration within HCC can reveal the highly heterogeneous landscape of HCC and improve the response rate to immune checkpoint blockade therapies to allow tailoring of immunotherapeutic strategies for patients [19,20]. Additionally, recent studies have shed light on the close relationship between immune-infiltrating cells within the TME and m6A modification. Tong et al demonstrated that METTL3, an m6A writer, promotes degradation of suppressor of cytokine signaling protein transcripts by catalyzing their m6A modification, promoting the differentiation of naïve T cells, and sustaining the suppressive functions of Tregs [21,22]. Deletion of METTL3 induces severe autoimmune diseases. Since Tregs within the TME are vital for the suppression of tumor-killing effector T cells, it is reasonable to speculate that the selective depletion of m6A modulators in Tregs could benefit patients with cancer.

However, because of technical limitations, the above studies have necessarily been restricted to only 1 or 2 immune cells and a single m6A modulator, while the highly coordinated interaction of multiple tumor factors has been neglected. Therefore, the comprehensive exploration of the immune infiltration characteristics mediated by m6A modulators will help to increase our understanding of TME regulation. A clinical prognostic model can also guide patient outcomes and survival, while emerging bioinformatics resources could accurately provide a broader scale of the intratumor microenvironment landscape and avoid the limitation of tissue-based methods, such as the amount of tissue and cell type. In this study, we used integrated algorithms to evaluate transcriptional information from 365 patients with HCC and identified 2 m6A modification clusters with distinct TME immune infiltration characteristics and clinicopathological features. Moreover, LRPRRC, an m6A reader, was found to potentially correlate with innate immune activation and cytolytic activities. A predictive signature was constructed using 5 screened m6A modulators, and a nomogram was constructed combining the risk scores from the predictive signatures and other clinical features. The performance of both models was well verified.

Material and Methods

Data Acquisition and Curation Process

By searching The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) and the International Cancer Genome Consortium (ICGC, https://icgc.org/) databases, 2 HCC cohorts (TCGA-LIHC and LIRI-JP) with integrated clinicopathological information were obtained and included in the study [23]. The TCGA-LIHC cohort included a total of 365 HCC and 50 normal tissues after removal of samples without overall survival (OS) information and RNA-seq data. The corresponding public high-throughput information, including level 3 RNA-seq data (fragments per kilobase million [FPKM] value), somatic mutation data, and copy number variations (CNVs) were download from Genomic Data Commons using the “TCGAbiolink” R package [24]. The Ensembl IDs were converted into gene symbols based on the GENCODE project annotation file (version 22, GRCH38), and gene expression levels were log2(FPKM+1), transformed for narrowing the numeric span. For the LIRI-JP cohort, normalized RNA-seq data from 231 liver tumors (RIKEN, Japan) and 202 normal tissues with corresponding survival data were download using the Illumina HiSeq RNA-Seq platform. Signatures of interferon (IFN)-gamma response, transforming growth factor (TGF)-beta response, proliferation, wound healing, cancer-testis antigen (CTA) score, and intratumor heterogeneity were referenced from a previous study [25]. Metabolism, glycolysis, and autophagy related gene data were obtained from the Molecular Signatures Database (MsigDB, https://www.gsea-msigdb.org/gsea/msigdb). The immunohistochemistry information of selected m6A modulators at the translational level were obtained from the Human Protein Atlas (HPA, http://www.proteinatlas.org/) database.

Unsupervised Classification of m6A Methylation Modulators

Altogether, 21 m6A modulators were extracted from the transcriptome datasets, including 2 erasers (FTO and ALKBH5), 11 readers (YTHDC2, YTHDC1, YTHDF2, YTHDF1, YTHDF3, HNRNPA2B1, IGF2BP1, HNRNPC, LRPPRC, ELAVL1, and FMR1), and 8 writers (ZC3H13, KIAA1429, CBLL1, WTAP, RBM15, RBM15B, METTL14, and METTL3). Based on the expression of the 21 m6A modulators, hierarchical clustering was applied for classification of the TCGA-LIHC cohort to identify different m6A modification patterns. To guarantee robust classification, we used an unbiased and unsupervised consensus manner implemented in the “ConsensusClusterPlus” R package with cluster algorithm=pam and correlation method=Euclidean [26]. The cumulative distribution function curve and gap statistic were used to select the optimal number of clusters (k).

Gene Set Variation Analysis and Gene Set Enrichment Analysis

To investigate the potential biological mechanisms underlying distinct m6A phenotypes, gene set variation analysis (GSVA), as implemented in the “GSVA” package, was performed. The GSVA algorithm is commonly applied for evaluating the biological processes and pathway variation in a distinct sample population, using an unsupervised and non-parametric approach [27]. The gene sets of “h.all.v6.2.symbols” utilized in the above steps were downloaded from MsigDB. Gene set enrichment analysis (GSEA) (version 3.0) under the JAVA platform was performed to reveal the pathway differences between patients with HCC and high and low m6A modulator expression in a genome-wide level [28]. The “c2.cp.kegg.v6.2.symbols” annotated gene sets were also obtained from the MsigDB database. Adjusted P<0.05 was considered statistically significant.

Estimation of Immune Microenvironment in HCC

Stromal and immune scores were calculated to quantify the proportion of infiltrating stromal and immune components in HCC by using the “ESTIMATE” R package [29]. The cytotoxic activity (CYT) response was assessed by the geometrical mean of GZMA and PRF1 [30]. Immune activity and infiltration assessment were conducted using the single-sample GSEA program, which allows robust quantification of infiltration abundance of various immune cell populations in individual samples from the transcriptomic matrix [30]. Relevant marker gene signatures of immune cells for the single-sample GSEA algorithm were retrieved from the work of Bindea et al [31]. The immune cells evaluated in this study comprised adaptive immunity and innate immunity. Adaptive immunity included T cells, B cells, effector memory T cells (Tem), central memory T cells (Tcm), cytotoxic cells, CD8 T cells, Th17 cells, Th2 cells, Th1 cells, Treg cells, and T follicular helper cells. Innate immunity included CD56dim natural killer (NK) cells, NK cells, CD56bright NK cells, dendritic cells (DCs), immature DCs, activated DCs, plasmacytoid DCs, mast cells, neutrophils, eosinophils, and macrophages.

m6A Risk Model Generation Using LASSO Regularization

The Cox regression method, with least absolute shrinkage and selection operator (LASSO) regularization, was used to penalize the weight of model parameters and select the most powerful m6A signature prognostic biomarker [32]. We used the “glmnet” R package to fit the LASSO Cox regression model. By 10-fold cross validations, the optimal penalty parameter (λ) was determined, thus generating a sparse parameter space. In this method, the characteristics of m6A-related biomarkers involved in HCC were selected by shrinking the regression coefficient using the penalty proportional to their size. Finally, the genes represented by λ were picked to establish the m6A panel in the training set. The risk score formula, based on the m6A panel, was constructed by integrating the normalized gene expression levels and their regression coefficients:

The “survivalROC” R package was used to determine the optimal cutoff value for the m6A risk score, and the patients were separated into 2 groups by high-risk score and low-risk score [33]. Univariate and multivariate analyses were used to ascertain the independent prognostic capacity of the m6A risk score in a Cox proportional hazard model with the “LR forward” method and visualization by employing the “forestplot” R package.

Construction and Evaluation of the Nomogram

The independent prognostic factors that were identified by multivariate Cox analysis were chosen to establish the nomogram for predicting the OS of patients with HCC in a quantitative way. The consistency between the frequencies of the probabilities and actual survival outcomes of the nomogram prediction was assessed by calibration plots. The C-index was used to assess the stability and discrimination of the model prediction. The nomogram construction and evaluation were produced by the “rms” R package.

Statistical Analysis

The unpaired t test (for normally distributed variables) and Mann-Whitney U test (for non-normally distributed variables) were used to compare groups and determine statistical significance. For comparison among more than 2 groups, we used non-parametric Kruskal-Wallis tests. The Benjamini-Hochberg method was used for multiple testing. Correlations between m6A modulators and immune cell infiltration levels were calculated by Spearman’s correlation and distance analyses. The principle component analysis, based on specific genes, was used to distinguish tumor tissue from normal tissue in patients with HCC. The survival curve was generated using the Kaplan-Meier program, and the log-rank test was used to determine the statistical significance of differences. We applied univariate Cox regression analysis to compute the hazard ratios for m6A modulators. A receiver operating characteristic (ROC) curve was plotted to evaluated the sensitivity and specificity of the m6A risk score, and the area under the curve (AUC) was generated using the “ROC” R package. The single nucleotide polymorphism (SNP) profile and CNV landscape in human chromosomes of the TCGA cohort for 21 m6A modulators were visualized by the R packages of maftools and RCircos, respectively. All P values were 2-sided, and P<0.05 was considered statistically significant. In this study, the primary clinical endpoint was set as OS, and the secondary endpoints were progression-free survival (PFS) and disease-free survival (DFS).

Results

Landscape of Genetic Alteration by m6A RNA Methylation Modulators in HCC

A total of 21 m6A modulators, including 2 “erasers” (m6A demethylases), 11 “readers” (m6A binding proteins), and 8 “writers” (adenosine methyltransferases), were incorporated into this study. Aberrant m6A modulators can result in tumor occurrence and progression, so we first analyzed the incidence of somatic mutations and CNVs in 21 m6A modulators in HCC. The CNV investigation revealed a more widespread alteration frequency of m6A modulators in HCC tissues than in normal tissues (adjusted P<0.05). Specifically, we observed that KIAA1429, YTHDF3, HNRNPA2B1, CBLL1, YTHDF1, IGF2BP1, YTHDC2, and FMR1 had a prevalent frequency of CNV amplification, while ZC3H13, ALKBH5, WTAP, FTO, METTL14, YTHDF2, and YTHDC1 mainly had deletions in copy number (Figure 1A). The altered locations of CNVs of m6A methylation modulators at chromosomes is demonstrated in Figure 1B. In addition, we analyzed the correlation between altered CNVs in m6A modulators and their mRNA expression levels and found that most m6A modulators had markedly higher mRNA expression, with amplification of CNVs in HCC. Regarding SNPs, we observed a low somatic mutation frequency in the 21 modulators in HCC, with mutations evident in only 44 (12.09%) of the 364 HCC samples (Figure 1C). Notably, compared with normal liver tissues, all 19 m6A modulators (except for ZC3H13 and METTL3) demonstrated significantly higher expression in HCC (Figure 1E). Using principal component analysis, HCC tissues could be completely distinguished from normal liver tissues based on the expression of these m6A modulators (Figure 1D). As a contrast, we performed principal component analysis again using 21 randomly selected genes and found the tumor sample was mixed with the normal sample (Supplementary Figure 1A). These results indicated that deregulated m6A modulators may play important roles in HCC initiation and progression and that CNV alteration might be a potential factor leading to the perturbation in, and heterogeneity of, m6A modulator expression in HCC.

Figure 1

The landscape of N6-methyladenosine (m6A) methylation modulator-related genetic aberration in hepatocellular carcinoma (HCC). (A) The frequency of copy number variation (CNV) for 21 m6A modulators in The Cancer Genome Atlas (TCGA) cohort. The red point represents amplification frequency and the blue point represents deletion frequency. (B) Circle plot of the specific location of the CNV of 21 m6A modulators in the human chromosomes. (C) The somatic mutation profile of m6A modulators in 364 patients with HCC from the TCGA cohort. (D) Principal component analysis for samples from International Cancer Genome Consortium and TCGA cohorts. Tumor samples could be well distinguished from normal samples based on the expression profile of the 21 m6A modulators. Normal samples are labeled with yellow and tumor samples are labeled with blue. (E) Expression variation of m6A modulators: comparison between normal tissues and tumor tissues from the TCGA cohort. The black lines in boxes represent the median value and black points represent the outliers. Statistical significance is represented by asterisks (*** P<0.001; ** P<0.01; * P<0.05). R (version 3.6.1) software was used to create the pictures.

Distinct m6A Modification Patterns Showed Different Prognosis Benefits in HCC

We performed univariate Cox regression analysis to explore the prognostic significance of 21 m6A RNA methylation modulators in patients with HCC. By applying hierarchical and K-means cluster analysis, a comprehensive characterization of the internal association of m6A modulators and their prognostic value for patients with HCC was demonstrated (Figure 2A). The 21 m6A modulators were positively correlated with HCC, and 13 modulators served as adverse prognostic factors for patients with HCC (HR>1, P<0.05).

Figure 2

Survival outcomes and biological characteristics of distinct N6-methyladenosine (m6A) methylation modification patterns. (A) Interaction of 21 m6A methylation modulators in hepatocellular carcinoma (HCC). The node size, measured by log10 (P value), represents the impact of each m6A modulator on prognosis. The black and green dots represent overall survival risk and protective factors, respectively. The thickness of the lines indicates the correlation strength between m6A modulators. Blue lines, positive correlation; red lines, negative correlation. (B) Consensus classification of patients with HCC for k=2. (C–E) Kaplan-Meier survival analyses for distinct m6A modification patterns in The Cancer Genome Atlas cohort. The m6A cluster B presents worse (C) overall survival, (D) progression-free survival, and (E) disease-free survival than the m6A cluster A. (F) Gene set variation analysis enrichment illustrates the activation score of biological function between 2 methylation patterns and is visualized in the heatmap. The activated pathway is marked with gold and the inhibited pathway is marked with blue. (G) Immune and stromal component differences between 2 methylation patterns. R (version 3.6.1) software was used to create the pictures.

To clarify the molecular heterogeneity of patients with HCC, we applied a consensus unsupervised approach to explore whether m6A RNA methylation modification presented discernable patterns. Based on the transcriptional expression of 21 m6A modulators, the consensus matrix revealed that k=2 was an optimal selection and reflected a balanced partition (Figure 2B). Therefore, 2 distinct m6A modification patterns from the TCGA cohort, termed m6A cluster A (n=191) and m6A cluster B (n=174), were determined. We found that m6A cluster B showed higher expression of the 21 m6A modulators. This suggested that patients in cluster B may have higher m6A-related RNA methylation status (Supplementary Figure 1B). The Kaplan-Meier curve of the relationship between the patterns and prognosis revealed that m6A cluster A presented a particularly prominent survival advantage in OS (P=0.00067) and PFS (P=0.0048) (Figure 2C, 2D). For DFS, we observed a tendency toward patients in m6A cluster B having a worse survival benefit, but it was not statistically significant (P=0.059) (Figure 2E).

The Landscape of Immune Characteristics in Patients with Different m6A Modification Patterns

To investigate the biological behaviors, heterogeneity, and different survival outcomes among the m6A modification clusters, we performed GSVA enrichment analysis based on the HCC RNA-seq profiles. We found that m6A cluster A was enriched in immunogenic related pathways, including inflammatory response, IL-6/JAK/STAT3 signaling, interferon (IFN)-α response, IFN-γ response, and TNF-α signaling via NFKB (Figure 2F), while m6A cluster B was prominently enriched in pathways associated with cell-cycle regulation and tumor promoting, including the G2M checkpoint, MYC target V1, TGF-beta signaling, WNT beta catenin signaling, PI3K/AKT/MTOR signaling, and MTORC1 signaling. Using the ESTIMATE algorithm, we further compared differences in immune and stromal components between the 2 m6A clusters. Consistent with the GSVA results, m6A cluster A led to a markedly higher immune and stromal score than did m6A cluster B (Figure 2G). These results strongly suggested that distinct m6A modification patterns may not only be associated with hepatocarcinogenesis, but also with tumor TME formation in patients with HCC.

We sought to better delineate the immune characteristics of m6A clusters. Single-sample GSEA was used to evaluate immune cell infiltration in HCC (Figure 3A). By K-means and hierarchical cluster analysis, we constructed a comprehensive HCC TME landscape, demonstrating the interaction between 24 tumor-infiltrating immune cells and their prognostic value for patients with HCC (Figure 3B). Notably, patients in m6A cluster A had a significantly higher proportions of B cells, cytotoxic cells, immature DCs, mast cells, neutrophils, plasmacytoid DCs, T cells, and Th17 cells than did those in m6A cluster B (P<0.05) (Figure 3C). Increased expression of these immune cells generally led to improved survival in patients with HCC. Additionally, we found that patients in m6A cluster A had a higher IFN-gamma response and a lower TGF-beta response, CTA score, proliferation signature, wound healing signature, and intratumor heterogeneity than did those in m6A cluster B (Figure 3D). These results suggest that patients in m6A cluster A exhibited a distinct immune phenotype, characterized by increased immune activation, cytotoxic potential, and immune infiltration.

Figure 3

Immune cell infiltration characteristics in distinct N6-methyladenosine (m6A) methylation modification patterns. (A) Unsupervised classification of patients with hepatocellular carcinoma (HCC) from The Cancer Genome Atlas cohort using normalized single-sample gene set enrichment analysis scores of 24 types of immune cells. Patients were classified as having high-, median-, and low-immune infiltration status. (B) The interaction among the 24 immune cell types in the HCC tumor microenvironment. The node size was calculated by Log10(log-rank P value) and represents the impact of each immune cell type on prognosis. (C) Differences in the abundance of immune cell infiltration between the 2 m6A modification patterns. (D) Differences in IFN-gamma response, TGF-beta response, proliferation, wound healing, cancer-testis antigen score, and intratumor heterogeneity signatures between the 2 m6A modification patterns. The P value is represented by asterisks (*** P<0.001; ** P<0.01; * P<0.05). R (version 3.6.1) software was used to create the pictures.

High LRPPRC m6A Modulator Expression is Associated with DC Depletion and Impaired Cytolytic Activity

Inspired by the intimate cross-talk of m6A modification patterns and immune features, we investigated whether m6A modulators affected the immune microenvironment in HCC. The specific relationships between the 21 m6A modulators and the infiltration level of assorted immune cells was examined (Supplementary Figure 2A). We focused on 5 modulators (CBLL1, LRPPRC, METTL3, KIAA1429, and YTHDF1) because their high expression was significantly associated with proportions of depleted B cells, T cells, neutrophils, DCs, and cytotoxic cells. CYT response, which is responsible for effective natural antitumor immunity, is characterized by dramatically enhanced T cell activation and is associated with improved patient prognosis [30]. We found that elevated CBLL1, LRPPRC, and METTL3 expression was associated with impaired CYT in patients with HCC (Supplementary Figure 2B). DCs, which function as a bridges of adaptive and innate immunity, are responsible for naïve T cell activation and antigen presentation, and their activation requires abundant major histocompatibility complex (MHC) molecules [34].

Our results indicated that high LRPPRC expression in patients was significantly related to reduced infiltration of diverse DCs, including plasmacytoid DCs and immature DCs. We also observed that elevated LRPPRC expression tended to correspond with downregulated MHC molecule expression in patients with HCC (Figure 4A), whereas no significant association was observed between CBLL1 and METTL3 expression and MHC (Supplementary Figure 2C). Patients with high LRPPRC expression also showed lower immune scores (Figure 4B), indicating that the TME in these patients had decreased immune cell infiltration. Importantly, the GSEA results also indicated that tumorigenesis and immune-associated pathways were significantly enriched in patients with high LRPPRC expression, including MAPK, MTOR, ERBB, WNT, T cell receptor, and TGF-beta signaling and pathways in cancer (Figure 4C). We further analyzed the expression pattern of altered LRPPRC genomic targets in HCC tissues and noncancerous liver tissues at the protein level by using the HPA database. Immunohistochemistry staining results showed a significant upregulation of LRPPRC expression in HCC tissue (Figure 4D), which was consistent with the corresponding transcriptional results. Collectively, it is reasonable to speculate that LRPPRC-mediated m6A methylation may promote tumor progression and immune suppression characteristic shaping, which may have been responsible for its unfavorable prognoses with patients in the TCGA-LIHC cohort (hazard ratio [HR]=1.8; log-rank P<0.001) (Figure 4E).

Figure 4

Potential roles of the LRPPRC N6-methyladenosine (m6A) modulator in immune microenvironment formation and hepatocarcinogenesis. (A) Differences in human leukocyte antigen molecule expression between LRPPRC high-expression and low-expression groups. (B) Differences in immune components between LRPPRC high-expression and low-expression groups. (C) Gene set enrichment analysis reveals the significant Kyoto Encyclopedia of Genes and Genomes pathways enriched in patients with high LRPPRC expression. (D) Representative LRPPRC immunohistochemical staining in normal liver tissues and those from patients with hepatocellular carcinoma in the Human Protein Atlas database. (E) Kaplan-Meier curve survival analyses for The Cancer Genome Atlas cohort patients with high and low LRPPRC expression. R (version 3.6.1) software was used to create the pictures.

Establishment and Validation of an m6A-based Prognostic Model

To further explore the clinical application of m6A signature in patients with HCC, we performed LASSO Cox regression on 13 prognostic m6A methylation modulators for dimension reduction. Patients in the TCGA-LIHC cohort were randomly divided into a training set (n=184) and testing set (n=181). Comparisons of patient characteristics between the 2 sets showed no significant differences. By constructing a penalty parameter, the regression coefficients were compressed to less than a fixed value in the LASSO model (Figure 5A, 5B) and cross validation was performed to avoid over-fitting. The 5 greatest prognostic features, LRPPRC, METTL3, YTHDF2, YTHDF1, and YTHDC1, were finally determined in training sets with individual non-zero coefficients (Figure 5C). Conversely, potential biomarkers with regression coefficients of zero were eliminated. Then, based on the formula described earlier, we calculated the m6A risk score of each HCC patient in the TCGA training sets. The optimal cutoff point (3.870) was determined using the “survivalROC” R package (Figure 5D), and patients in the training sets were stratified into high-risk and low-risk groups (Figure 5E). The Kaplan-Meier survival curve indicated that patients in the high-risk group showed a significantly worse survival outcome than did those in the low-risk group (P<0.0001) (Figure 5F). ROC curves were generated to assess clinical predictive performance, and the AUC for OS was 0.727, 0.709, and 0.683 at 1, 2, and 3 years, respectively, indicating a good forecasting ability (Figure 5G).

Figure 5

Establishment and validation of the prognostic panel using 5 N6-methyladenosine (m6A) modulators in The Cancer Genome Atlas (TCGA) cohort. (A, B) The least absolute shrinkage and selection operator (LASSO) Cox regression model identified 5 core prognostic m6A modulators in the TCGA training set. (C) The corresponding regression coefficients: YTHDF2, 0.6744; YTHDF1, 0.1318; YTHDC1, −0.4059; METTL3, 0.1954; and LRPPRC: 0.3962. (D) The optimal cutoff point (3.870) could distinguish patients with high and low risk. (E) Risk score distribution and survival overview for patients in the TCGA training set. (F) Prognostic analysis showed that the overall survival of patients was significantly lower in the high-risk score group than the low-risk score group in the TCGA training set. (G) Receiver operating characteristic curve was used to assess the predictive performance of m6A risk score. (H–J) The predictive performance of the m6A risk score was validated in the TCGA testing. R (version 3.6.1) software was used to create the pictures.

To confirm our TCGA training cohort findings, we validated the m6A-related prognostic models in the TCGA testing and ICGC cohorts. The same risk score calculation formula and cutoff values were applied to distinguish risk groups (Figure 5H and Supplementary Figure 3A). Survival plots indicate that, compared with those in the low-risk group, patients in high-risk group had significantly worse OS in the TCGA testing cohort (P<0.0001, Figure 5I) and the ICGC cohort (P=0.0087, Supplementary Figure 3B). The AUC in the TCGA testing cohort was 0.748, 0.701, 0.726 at 1, 2, and 3 years, respectively (Figure 5J). The AUC for the ICGC cohort was 0.737, 0.745, and 0.713 at 1, 2, and 3 years, respectively (Supplementary Figure 3C). In the meantime, we further investigated the prediction ability of the model for PFS and DFS in patients with HCC, and our results indicated that patients in the low-risk group demonstrated a particularly prominent survival advantage in both PFS (Figure 6A) and DFS (Figure 6B). Therefore, our results indicated that m6A-based prognostic models present relatively robust and pleiotropic clinical predictive efficiency. To illuminate the characteristics of the m6A risk score, we analyzed the specific correlation between TME immune infiltration and m6A risk score. The result revealed that a high m6A risk score was significantly negatively associated with tumor-killing immune cell infiltration (T cells, B cells, CD8 T cells, cytotoxic cells, DCs, and neutrophils), but also with CYT response (Figure 7A). We also found that patients with an advanced pathologic stage and grade exhibited a higher m6A risk score in the TCGA cohort (Figure 7B, 7C). The Sankey map showed that patients with high m6A risk scores were mainly linked to m6A cluster B and low immune-infiltrating subtypes, which were associated with poor survival status (Figure 7D).

Figure 6

(A, B) Kaplan-Meier curves for progression-free survival and disease-free survival of patients with hepatocellular carcinoma (HCC) indicated that those with high N6-methyladenosine (m6A) risk scores had a worse outcome than did those with low m6A risk scores. R (version 3.6.1) software was used to create the pictures.

Figure 7

The correlation between clinical characteristic and N6-methyladenosine (m6A) risk score. (A) The m6A risk score was significantly negatively correlated with infiltration of protective immune cell types in hepatocellular carcinoma. (B) The correlation between m6A risk score and pathologic stage. (C) The correlation between m6A risk score and pathologic grade. (D) Sankey diagram of m6A risk score in classification with different methylation modification patterns, immune infiltration status, and survival outcomes. R (version 3.6.1) software was used to create the pictures.

Construction of an m6A-based Nomogram in Clinical Practice

We applied univariate and multivariate Cox regression analyses to test whether the m6A risk score could serve as a biomarker to independently predict OS for patients in the TCGA-LIHC cohort (Figure 8A, 8B). The results of the multivariate Cox model, incorporating the clinical information about sex, age, pathologic grade, tumor-node-metastasis (TNM) status, alcohol consumption, immune score, m6A cluster, and CYT response, suggested that the m6A risk score was an independent prognostic factor, confirming its robust predictive efficiency for OS in patients with HCC (HR: 3.82, 95% confidence interval [CI]: 2.05–7.11, P<0.01). To develop a more sensitive quantitative method for the clinical forecast of mortality of patients with HCC, we constructed a nomogram that integrated classical TNM pathologic stages and independent prognostic factors (age, CYT response, and m6A risk score) in the TCGA cohort (Figure 8C). The C-index of the constructed nomogram was 0.705 with 1500 bootstrap iterations (95% CI: 0.6–0.76), which was better than the predictive performance of pathologic stage (C-index: 0.61). The calibration plots also indicated that our nomogram performed with good consistency when comparing predicted survival and actual observed outcomes (Figure 8D–8F). We found that the m6A risk score exhibited the greatest weight for the total points in the nomogram, which was consistent with our multivariate regression model results.

Figure 8

Construction of a nomogram integrating clinicopathological features and N6-methyladenosine (m6A) risk score. (A, B) Univariate and multivariate Cox regression analyses of the correlation between m6A risk score and clinicopathological features in terms of overall survival (OS). (C) Nomogram constructed for the prediction of OS at 1, 2, and 3 years in patients with hepatocellular carcinoma. (D–F) Calibration curve for evaluating the predictive ability of nomogram for OS. R (version 3.6.1) software was used to create the pictures.

Discussion

The functional implications of m6A on the immune system are receiving increasing attention. Aberrant m6A modulator expression and dysregulated m6A modification play critical roles in inflammatory activation, immune system imbalance, and antitumor immune responses [11]. Most biological research has concentrated on a specific type of infiltrating cell within the TME or single m6A modulator. For example, Shen S et al recently identified METTL3 to be associated with the infiltration of DCs within the TME, which might be a promising immune therapeutic target [35]. Collective studies indicate that molecular patterns allow the classification of tumors into distinct phenotypes associated with diverse prognostic and clinicopathologic traits [36]. Therefore, identifying distinct m6A methylation clusters and their relationship with particular types of infiltrating cells within the TME can facilitate our understanding of m6A-regulated antitumor TME inflammatory responses and guide clinical decisions on immunotherapies.

Based on the expression of 21 m6A modulators, 2 m6A modification clusters, with distinct survival outcomes and disease progression, were identified in HCC. The 2 clusters exhibited variations in TME components. Cluster A had higher immune and stromal scores and lower expression of m6A modulators than did cluster B. Moreover, cluster A was characterized by the activation of inflammatory response and tumor stroma, consistent with its higher immune and stromal scores. Several highly expressed signaling pathways, including KRAS, IL6/JAK/STAT3, IFN-α, IFN-γ, and TNF-α, interacted to trigger multiple inflammatory responses in cluster A, which could eventually inhibit the progression of HCC to some extent. The favorable prognosis associated with cluster A also suggested that activated antitumor inflammatory response had overridden tumor stroma-induced tumorigenesis to inhibit tumor development. Cluster B, with a higher expression of m6A modulators, was characterized by an immune suppressive state, potentially indicating a cold tumor with low immune infiltration within the TME. It was reported that Wnt-β-catenin and TGF-β signaling pathways, which were significantly activated in cluster B, are associated with reduced cytotoxic T cells and can cooperatively promote tumor growth [37,38].

To characterize the infiltrating cells within the TME, 24 cell types were compared between the 2 distinct methylation patterns. Significantly higher proportions of neutrophils, DCs and T, B, cytotoxic, and mast cells were found in the m6A cluster A, consistent with an activated immune response. Therefore, m6A modification may play an important role in the maturation of innate immune cells, such as DCs, which are responsible for tumor antigen presentation and activation of adaptive immune cells. Our signature scores showed phenotypes that were more aggressive in cluster B, including higher proliferation, wound healing, and intratumoral heterogeneity, which might lead to decreased CD8+ T cell infiltration and weakened tumor-killing effects [39,40]. These results indicated that m6A modification is crucial for regulating tumor biology by shaping the tumor immune microenvironment. Previous studies also reported that m6A modification actively participates in innate immunity by regulating immune transcript translation. Wang et al showed that METTL3-mediated m6A methylation promotes DCs maturation and stimulates T cell activation [41]. Additionally, Han et el reported that YTHDF1, an m6A reader, promotes translation of lysosomal cathepsins in DCs by binding to their m6A-modified transcripts, which then promotes tumor progression [42]. Depletion of YTHDF1 enhances the cross-presentation of DCs, strengthens the CD8+ T cell antitumor effect, and improves therapeutic efficacy of immune checkpoint blockade therapies. We identified the distinct characteristics of infiltrating cells within the TME induced by different m6A methylation patterns, and determined that patients from cluster A, with highly activated innate immunity, are potential candidates for m6A-relevant immunotherapies.

We then examined the 21 m6A modulators to identify dominant factors influencing the immune responses. Five m6A modulators, CBLL1, LRPPRC, METTL3, KIAA1429, and YTHDF1, were strongly correlated with DCs, cytotoxic cell infiltration, and CYT response. This showed the significance of m6A modulators in the activation of innate immune systems and their subsequent antigen presentation processes for effector tumor-killing cells. LRPPRC, a member of the PPR family, was identified as a reader protein capable of recognizing m6A modifications [43]. Although there is no direct evidence linking LRPPRC with tumor innate immune response, several studies have reported that LRPPRC is overexpressed in a variety of cancers and is associated with a series of malignant behaviors including apoptosis resistance, tumor invasion, and proliferation [44-46]. In our study, LRPPRC was negatively correlated with expression of the HLA family in HCC, which encode the MHC for antigen presentation, indicating its effect on innate immunity [47]. Recent studies demonstrated that DC activation by pro-inflammatory molecules causes them to switch metabolic sources from oxidative phosphorylation (OXPHOS) toward glycolysis, and LRPPRC is critically involved in tumor energy metabolic processes, including OXPHOS and FAO [48,49]. Hoss et al demonstrated that alternative splicing of leucin-rich repeat domains can regulate vertebrate innate immunity [50]. Together, these results suggest that aberrant LRPPRC and m6A modification exert inhibitory effects on the activation of innate immune systems and hinder the antitumor inflammatory effect.

Given the individual heterogeneity of m6A methylation in the HCC TME, it is urgent and necessary to develop methods to quantify m6A risks and to personalize immunotherapeutic strategies for patients with HCC. In the present study, a predictive signature based on 5 m6A modulators was constructed and then trained with the TCGA internal sets. The robustness of the survival forecast capacity of this signature was validated using the TCGA internal cohort and ICGC external cohort. This signature demonstrated significant differentiation capacity for DFS and PFS in patients with HCC. Also, the m6A risk score was closely related to the immune infiltration level and pathological stage in HCC. It can not only predict patient OS but can also assess the levels of methylation and immune infiltration in patients with HCC. Therefore, the 5 modulators used in the risk score, YTHDF1, YTHDF2, YTHDC1, LRPPRC, and METTL3, could act as individual targets or be targeted in combination to produce higher efficacy for immunotherapies. The potential mechanisms of m6A modulator involvement in HCC immune responses have been well researched. Our integrated analysis confirmed the predictive accuracy and reliability of our signature, which could be used to further determine the immune landscape of HCC. Our results also demonstrated that m6A risk score was an independent prognostic factor in HCC and was combined with other clinicopathological factors including age, pathologic stage, and CYT to build a nomogram. The nomogram showed excellent consistency between real and predicted OS at 1, 2, and 3 years. Our nomogram provides a personalized risk score for patients, which might be valuable for guiding treatments and clinical decisions.

Our results shed new light on epigenomic modification and on emerging immunotherapies in HCC. First, our study provided a comprehensive evaluation of m6A patterns and the corresponding TME landscape, revealing m6A clusters with distinct m6A and TME characteristics, which could be potential candidates for m6A-relevant immunotherapies. Second, strong correlations between specific m6A modulators and innate immune responses (such as LRPPRC and DC activation) were revealed, which could be utilized as novel immunotherapeutic targets. Third, a predictive nomogram with excellent performance was built, which could be tailored for personalized treatments and prognostic prediction in patients with HCC.

Additionally, we explored the relationship between m6A risk score and clinical characteristics. We observed that the developed risk score is proportional to patient pathologic stage and grade and has high clinical efficacy. Prospective clinical trials and in vivo studies are required to validate the prognostic nomogram and the potential relationship between m6A and tumor immunotherapies.

Conclusions

This work illustrates the regulatory effect of m6A modification on TME characterization. Distinct m6A patterns play indispensable roles in the formation of heterogenous and complex TMEs. To the best of our knowledge, this is the first report of a comprehensive analysis of m6A modification in the HCC TME. The m6A risk score panel we constructed not only predicted survival of HCC patients but also evaluated the antitumor immune infiltration level and methylation pattern.

(A) Principal component analysis for The Cancer Genome Atlas cohort based on 21 randomly selected genes as the contrast effect. (B) Heatmap shows that N6-methyladenosine (m6A) cluster B has higher expression of 21 m6A modulators than does cluster A. R (version 3.6.1) software was used to create the pictures.

(A) The correlation between N6-methyladenosine (m6A) modulator expression and infiltrating levels of immune cells in the hepatocellular carcinoma tumor microenvironment. (B) Correlation between the expression of selected m6A modulators, CBLL1, LRPPRC, METTL3, and cytolytic activity response. (C) Differences in human leukocyte antigen molecule expression between high- and low-expression groups of CBLL1 and METTL3. R (version 3.6.1) software was used to create the pictures.

External validation of the N6-methyladenosine (m6A)-based prognostic model in the International Cancer Genome Consortium (ICGC) cohort. (A) The risk score rank and survival overview of individual patients in the ICGC cohort. Patients with hepatocellular carcinoma (HCC) in the ICGC cohort were well stratified into high- and low-risk groups using the same cutoff value applied in The Cancer Genome Atlas cohort. (B) Patients with HCC in the low-risk group showed favorable survival outcomes compared with those in the high-risk group in Kaplan-Meier curves (P<0.01). (C) The area under the curve of overall survival at 1, 2, and 3 years in receiver operating characteristic curves were 0.737, 0.745, and 0.713, respectively. R (version 3.6.1) software was used to create the pictures.