M6a Regulator-Associated Methylation Modification Patterns Shape Tumor Microenvironment Characteristics in Hepatocellular Carcinoma

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

Download PDF

Research Article

M6a Regulator-Associated Methylation Modification Patterns Shape Tumor Microenvironment Characteristics in Hepatocellular Carcinoma

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

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Although substantial achievements in the tumor microenvironment (TME) of hepatocellular carcinoma (HCC) have led to fundamental improvements both in the basic research and clinical management, the potential mechanisms and regulatory relationships between m6A regulators and the TME are still unknown. We first conducted unsupervised clustering on the samples according to the core m6A expression, and then compared the signaling pathways, differential genes (DEGs), and TME between the m6A phenotypes, and re-validated the relationship between m6A regulators and TME by single cell sequencing. Then, the geneCluster was obtained by another unsupervised clustering of the DEGs, and the clinical as well as TME traits were evaluated among the geneClusters. Finally, the m6A scores of individual patients were calculated by principal component analysis (PCA) to verify the correlation from multiple perspectives, including survivals, clinical characters, mutations, TME, immunotherapy, and chemotherapy. Through a comprehensive analysis of 729 samples, we classified HCC patients into three m6A clusters and three geneClusters. Each group exhibited remarkable variations in terms of signaling pathways, clinical traits, and survival expectations. Notably, the m6A phenotypes corresponded to three different types of TME, namely immune-inflamed, immune-excluded, and immune-desert, respectively. In addition, the m6A regulator can accurately reflect the individualized microenvironment in HCC, and present supreme expression levels in the stromal microenvironment. However, the m6A score system is able to make accurate predictions not only in terms of clinical traits, survival prediction, and TME mentioned above, but also in the sensitivity of HCC patients to immunotherapy and chemotherapy. This study revealed the uniqueness and pluripotency of m6A regulators in the TME of HCC by combining single-cell sequencing and bulk sequencing. The quantified m6A modification indices were able to accurately predict patient survival expectations, clinical traits, TME, and sensitivity to immunotherapy and chemotherapy.

Scientific Communication

General Surgery

Gastroenterology & Hepatology

tumor microenviorment

m6A

hepatocellular carcinoma

immunotherapy

chemotherapy

Hepatocellular carcinoma is the most frequent primary liver cancer, accounting for 75%-85% of cases (1). The incidence of hepatocellular carcinoma ranks sixth worldwide and is consistently rising over time (2). Estimates from the World Health Organization suggested that liver cancer will take a roll of over one million in 2030 (3). The dismal 5-year survival rate of HCC was less than 18%, which remains the second malignant carcinoma-causing death across the human solid tumors, merely secondary to pancreatic cancer (4). This is in part results from the fact that most patients receive a diagnosis in their terminal stage, and are not able for curative-intent treatment like surgical resection and orthotopic liver transplantation(5). With the prosperity and development of immunotherapy, it has become indispensable in individualized comprehensive treatment regimens for cancer patients (6–8), multi-clinical trials have consistently provided proof that immunological therapeutic approaches achieve noticeable success in tackling malignant carcinomas (9, 10). However, the clinical benefits and sensitivities of patients to immunotherapy have always plagued clinicians' treatment decisions when determining treatment options (11, 12). This means the significance of in-depth evaluation of immune characteristics in the context of the tumor.

m6-methyladenosine (m6A) has been recognized as the most prevalent internal modification of eukaryotic ribonucleotides (RNA), comprising nearly 50% of the total methylated RNA (13). It has been well corroborated that m6A participates in a variety of post-transcriptomic activities, including the translation, splicing, degradation of both coding and non-coding RNA (14, 15). It has been proposed that m6A modification is regulated by a model of “Writer-Eraser-Reader” which is a dynamic process and can be reversed by a series of coordinated catalytic reactions (16, 17). m6A modification can be achieved by methyltransferases complex, known as writers. Demethylases, also termed erasers, are responsible for the removal of m6A. The biological effects of m6A modification are delicately controlled by various readers, which recognize and bind m6A sites, leading to different targets (18, 19). Accumulating evidence reveals the crucial role of m6A in developing a wide range of diseases, especially in malignant carcinomas (20–22). However, the understanding of the biological function of m6A is just in its fancy.

However, the potential mechanism of m6A genes in HCC is continuously explored. Li et al. reported that Fat Mass and Obesity-associated protein (FTO) drive tumorigenesis and growth of hepatocellular carcinoma through inducing the demethylation of PKM2 messenger RNA (mRNA) (23). Liu et al. discovered that YTHDF1 increases FZD5 mRNA translation on the grounds of m6A modification and, subsequently, promotes the proliferation and metastatic formation of HCC through the WNT signaling pathway (24). Similarly, a recent study from Lan and colleagues indicated that KIAA1429 promotes HCC growth and metastasis via m6A methylating 3′ UTR of GATA3 pre-mRNA (25). Besides, emerging studies concerning the effect on antitumor therapy revealed that m6A regulators may play a crucial role in the therapeutic resistance of HCC (26–29). Furthermore, certain researchers even proposed that m6A modification can be served as a biomarker to HCC diagnosing and outcomes predicting (30–32).

Numerous proof-of-concept studies reveal that the TME is essential to the initiation and progression of malignant carcinomas. As critical components, immune cells, whose infiltration is delicately regulated by the perturbated expression of immune-related genes in the context of tumor, can facilitate a niche defining the fate of neoplastic cells. Emerging experiments have been suggested close crosstalk between m6A modification regulators and the biological function of intra-tumor immune infiltrates. Mettl3, an m6A methyltransferase, facilitates the activation and functional maturation of dendritic cells via increasing transcription of CD40, CD80, and TLR4 through m6A modification, thus strengthening the antitumor activity of CD8 + T cells (33). Dampening m6A-binding protein YTHDF1 in dendritic cells strikingly potentiated its capacity in cross-presentation of tumor antigens and cross-priming CD8 + T cells, and ythdf1-defective mice experienced an enhanced therapeutic efficacy of anti-PD-L1 treatment (34). A study focusing on the m6A demethylation revealed that depletion of FTO significantly increases the expression of critical pro-tumorigenic genes, such as PD-1 (PDCD1), CXCR4, and SOX10 in melanoma and sensitizes melanoma mice to PD-L1 checkpoint blockade simultaneously (35). Due to the complicated nature of tumorigenesis, we speculated that the development and progression of malignancy are not simply driven by a single type of m6A modification, but also coordinated with m6A pattern orchestrated by a group of m6A regulators instead. Thus, we conducted the present study to explore the effect of the comprehensive m6A modification pattern on the landscape of TME in HCC.

In this study, we deeply excavate the potential mechanisms and links between m6A modification and TME by comprehensively combining the results of bulk sequencing and single-cell sequencing. The bulk sequencing data together with clinical materials were retrospectively collected from TCGA, ICGC, and GEO databases. The results of bulk sequencing and single-cell sequencing both confirmed that m6A regulators were widely distributed in the TME and can clearly distinguish the special microenvironment characteristics of different patients. This means that the m6A modification pattern plays a crucial role in shaping immune characteristics in HCC. Furthermore, a m6A score system by PCA analysis exhibited an extremely powerful function in predicting clinical characteristics, TME, and even sensitivity to immunotherapy and chemotherapy among the HCC patients

Data collection and panorama of m6A regulators in HCC

The overall flowchart of the present study is shown in Figure 1A. The HCC meta-cohort was constructed by a total of 850 samples’ RNA-seq data from three public databases, including TCGA (424) (https://portal.gdc.cancer.gov/), GEO (166) (www.ncbi.nlm.nih.gov/geo/), and ICGC (260) (https://dcc.icgc.org/releases). Ultimately 729 HCC patients with clinical information were enrolling for further analysis (Table S1). Firstly, the transcriptomic profile, clinical information, and somatic mutation data of HCC patients from TCGA were extracted. To make the TCGA values more compatible with the GEO database, we converted the TCGA matrix from FPKM to TPM format. Then, differential expression analysis of 23 m6A regulators was performed between the normal and HCC patients using R packages including "limma", "reshape2" and "ggpubr". After accessing the somatic mutation patterns, the tumor mutation burden (TMB) of individual patients was identified by perl language. The mutations of m6A regulators were further retrieved and a waterfall graph was mapped using the "maftools" package. Additionally, the copy number variation (CNV) files of aforesaid patients were obtained from the XENA website (https://xena.ucsc.edu/). We then plotted the m6A-related CNV over the 23 human chromosomes via R package "RCircos". Following the acquisition of the R-seq data from the GEO and ICGC databases respectively, we merge the three datasets using the Combat method by R package “SVA”. Furthermore, we conducted a Mutation-Expression analysis to reveal the m6A expression variance between the mutant and wild groups. Through the Timer database (http://timer.cistrome.org/) we obtained the relevance of m6A regulators to immune cells in HCC. Finally, we demonstrated the expression distribution profile of m6A regulators in the TME by single-cell sequencing through the TISCH website (http://tisch.comp-genomics.org/).

Prognostic analysis and unsupervised clustering of m6A regulators

Cox and K-M survival analysis of HCC meta-cohort were conducted to identify the m6A regulators with prognostic significance using the R packages “limma”, “survivor”, and “survminer”. Taken together, we thoroughly revealed the characters, co-expression relationship, and prognostic predictions by prognosis network diagram. The R packages "igraph", "psych", "reshape2", and "RColorBrewer" were utilized for this operation. After the corroboration, an unsupervised cluster analysis was performed based on the m6A regulators. Above that, we investigated the survival differences, m6A expressed discrepancy, and clinical traits across the m6A clusters. The packages of “ConsensusClusterPlus”, “survival”, “survminer” and “pheatmap” were allocated.

Gene set variation analysis (GSVA) and KEGG annotation

To explore the changes in signaling pathways and the molecular biological processes across m6A clusters, the gene set "c2.cp.kegg.v6.2. symbols" was downloaded from the MSigDB database and used to run the GSVA analysis by the R packages "GSEABase", "GSVA" and "pheatmap". FDR < 0.05 was considered statistically significant.

Integrated Analysis of m6A clusters and TME in HCC

With the help of the "limma", "GSEABase", "GSVA" and "ggpubr" packages, the variations of the infiltrated immunocytes in the individual patient from TCGA were quantified through the ssGSEA method among m6A clusters. Then, we considered the core m6A regulators as a signature and validated them in two separate single-cell sequencing datasets of HCC. In GSE140228, Kruskal-Wallis test was conducted among the TME to confirm the difference of m6A distribution among the patients. Whereas, in GSE125449, we compared the m6A signature across the infiltrated cells.

Distinctly expressed genes among m6A clusters

Across the three m6A clusters, the DEGs were identified for deeper insight into the mechanism regulated by m6A genes. The threshold of FDR<0.01 was considered statistical significance. After intersection of the DEGs, we performed GO and KEGG enrichment analysis of these genes by "clusterProfiler", "org.Hs.eg.db", "enrichplot", "DOSE", and "ggplot2" packages. Thus, the biological functions and molecular signal pathways in which m6A-related genes are involved were determined.

Identification of prognosis-related genes and m6A geneCluster

Another univariate cox survival analysis was performed to identify the prognostic DEGs with a threshold of P value<0.05. Then, we conducted unsupervised clustering again base on the core m6A-related DEGs. The HCC meta-cohort was classified into three geneClusters. Sequentially, survival differences, clinical traits, and m6A regulators expression were examined. To reveal the TME characteristics among the subgroups, we further identified the expression of core genes in the Wnt, EMT, TGF-β, immune response, and immune system process pathways among geneClusters. To quantify the level of m6A modification, we constructed a scoring system by PCA. According to the m6A score, we selected an optimal cutoff value and divided the HCC meta-cohort into high- and low- m6A score groups. The survival differences between the two groups were exhibited by the K-M survival curves. The Sankey diagram was drawn by the "ggalluvial" package provided an overview of the entire process of completion. m6A clusters, m6A geneClusters, high- and low- m6A score groups, and the survival status were systematically demonstrated.

Multi-omics analysis between m6A score and TME, TMB, and clinical characteristics

The "ggpubr" package was used to determine m6A score distinctions across the m6A clusters and m6A geneClusters. Moreover, the relations between the m6A score and clinical information, such as survival status, T-stage, and Stage were firstly determined by "plyr", "ggplot2" and "ggpubr" packages. Additionally, we evaluated the relationship between the m6A score and TMB through "reshape2" package. We collated the TMB of HCC patients from TCGA and classified them into high- and low- TMB groups according to median values. The survival difference between the two groups was first determined. Furthermore, the survival analysis was conducted under the comprehensive perception of TMB and m6A score. Meanwhile, the waterfall diagrams of the TMB and mutation characteristics and the immune checkpoints gene, PD-L1, between the high- and low- m6A score groups were also presented. In the clinical subgroup analysis, we further divided the patients into low-grade groups (T1, T2) and high-grade groups (T3, T4) according to the T stage and separately constructed survival analysis to explore the differences in clinical subgroups.

Correlation analysis of m6A score with TME and m6A efficiency on immunotherapy and chemotherapy

The R packages “scales”, “ggplot2”, and “ggtext” were utilized to anticipate the relations between m6A score and TME by 6 separate algorithms. Then, ssGSEA analysis was further conducted to estimate the connection between the m6A score and TME infiltration cells. Information on the patients from TCGA was obtained from the TCIA (https://tcia.at/). In accordance to the given data, the patients were divided into four groups. And the immunotherapeutic response between the high- and low- m6A score groups was operated. Finally, the R package “pRRophetic” was applied to resolve the chemotherapeutic efficacy between the high- and low- m6A groups.

Landscape of genetic variation of m6A regulators in HCC

Start of the article, we scrutinized the landscape of m6A regulators between normal and HCC samples and observed that almost all of those gene expressions were upregulated in the tumor. Inversely, the expression of IGFBP1 and IGFBP3 were decreased in the tumor (Figure 1B). Throughout the entire HCC patients from TCGA, we observed that 31 patients experienced somatic mutations in the set of m6A regulators (8.52%) and HNRNPC was the most frequently mutated m6A gene. Our analyses revealed that the most common form of mutation was Missense Mutation, however, the overall mutation frequency of m6A regulators in HCC was significantly lower than that in gastric and colorectal cancer (Figure 1C). Moreover, we confirmed that significant mutational co-occurrence was existence between LRPPRC and METTL14; YTHDC1 and IGFBP1, ZC3H13; YTHDC2, and ALKBH5, FTO, LRPPRC, METTL3, METTL16 (FigureS1). CNV alternation analysis indicated that VIRMA, HNRNPC, and METTL3 dominated the list of most prevalent amplification, whereas ZC3H13, YTHDF2, and WTAP were found the most prevalent deletion (Figure 1D). We further located the 23 m6A regulators on the chromosomes along with their CNV alteration in Figure 1E.

After the cox univariate analysis, 12 m6A regulators were identified prognostic significant (Table S2). Eleven of them were considered as risk factors for the progression, except for ZC3H13. K-M survival analysis further revealed the survival status in HCC between the high- and low- m6A expression groups (Figure S2). Almost all the m6A regulators were positively correlated with immune cells (Figure S3).

To compensate the omission of the sparse cells infiltrated in the TME due to the limits of bulk sequencing, the distribution of m6A regulators in TME was depicted in Figure S4 under the single-cell sequencing data. Through the violin plot in Figure 1F, the expression and comparison of m6A genes in the TME can be intuitively visualized.

Together, our analyses suggest that expression alterations and genetic variation of m6A genes play an indispensable role in the development and progression of HCC.

Identification of m6A methylation modification patterns in HCC mediated by the 23 regulators

The m6A regulator network has systematically exhibited the interaction among the 23 m6A regulators as well as their prognostic significance in HCC (Figure 2A). To excavate tumorigenesis underpinning of m6A modification, we performed unsupervised clustering analysis based on the 12 core m6A regulators. 729 patients from the HCC meta-cohort were stratified into three categories, with 230 samples classified into cluster A, 211 samples into cluster B, and 288 samples into cluster C (Figure 2B, 2C). Meanwhile, we found that different m6A clusters were closely associated with survival outcomes of HCC patients. As the Log-rank test showed, cluster B had a worse survival rate than clusters A and C (Figure 2D). Additionally, cluster B had a stronger association with the progressive tumor. In Figure 2E, the m6A regulators were remarkably different among the three clusters and IGFBP3 was the most evident. Thus, it is reasonable to propose that the context of the m6A modification pattern mediates tumor aggressiveness and survival outcomes of HCC. Then, the core m6A regulators were regarded as a signature and imported into a single cell sequencing dataset consisted of five HCC patients. Although they were all in stage II, the distribution of the m6A signature among the TME was significantly different among those patients. This inspired us to assume that m6A modification not simply plays a crucial role in personalized tumorigenesis of HCC, but has a huge impact on the formation of the TME (Figure 2F).

The m6A clusters mediates distinct immune landscapes

Our observations suggested that cancer-related signaling pathways, such as Bladder Cancer, Renal Cell Carcinoma, and Pancreatic Cancer were significantly activated in cluster B, whereas the metabolism-related pathways were predominant in cluster A (Figure 3A). Consistent with the trend in prognosis, cell cycle signaling pathways, as well as metabolic pathways, were more active in group B than in group C (Figure 3B). Besides the cancer pathways and metabolic pathways, immune signaling pathways, such as B Cell Receptor Signaling Pathway and T Cell Receptor Signaling Pathway were significantly upgraded in cluster C which was the promising group in terms of both short- and long- survival (Figure 3C). The results of ssGSEA revealed the distribution of immune cell compartments among the three m6A clusters (Figure 3D). Almost all of the immune and stromal cells were significantly enriched in cluster B that associated with the worst survival. However, the m6A signature was abundantly expressed in stromal cells and immune cells also confirmed from the perspective of single-cell sequencing. Those findings proved that TME modulated by m6A might be essential to mediate the prognostic outcomes of HCC patients.

m6A phenotype-regulated DEGs and geneClusters construction in HCC

Our findings suggest that HCC patients can be stratified into three categories based on distinct m6A modification patterns, with each cluster epitomizing different clinical characteristics and survival outcomes. A total of 512 key DEGs that represented the genes subject to m6A regulator modifications were obtained for further analysis (Figure 4A). Growth factor regulation, cell membrane protein binding, and extracellular matrix on functions were mainly embroiled by Go enrichment analysis (Figure 4B). Additionally, KEGG enrichment analysis revealed that not only the classical cancer pathways eg. PI3K-Akt signaling pathway, MAPK signaling pathway, Ras signaling pathway, HIF-1 signaling pathway were concentrated, but also cancer-related macromolecular organic pathways, eg. EGFR tyrosine kinase inhibitor resistance, Proteoglycans in cancer, Central carbon metabolism in cancer (Figure 4C). The above results additionally proved that the genetic profile and molecular biological process under the m6A genes control may perform their functions in the progression of HCC by regulating TME, particularly stromal cells. 172 prognostic m6A-regulated genes were screened out for the geneCluster analysis (Figure 4D & Table S3). The unsupervised consensus clustering analysis classified the meta-cohort into three geneClusters (Figure 4E, F). As expected, geneCluster B manifested the worst clinical outcomes, no matter the aspect on clinical stage, tumor grade, or survival status. Intriguingly, while the survival rates were not statistically different in the first 4 years, the survival outcomes of geneCluster C dramatically deteriorated compared with Genecluster A (Figure 4G). In Figure 4H, geneCluster B exhibited the highest level of DEGs, followed by Group C, and Group A. Noteworthily, a significantly distinct expression pattern of the 23 m6A regulators was observed among the three subgroups (Figure 4I). According to Figure 4J~M, the pathway kernel genes in both immune and stromal aspects showed the same trend, which were negatively correlated with the survival status in each geneCluster.

Construction of the m6A score system and exploration of its clinical relevance

Although m6A modification exhibited its crucial role in modulating HCC prognosis and TME, the population-based analyses cannot accurately quantify the individual m6A phenotype. Therefore, we constructed a score scheme to reflect their m6A modification pattern by PCA and the Sankey diagram was drawn (Figure 5A). Firstly, we explored the ability of the m6A score system to predict the survival outcomes of HCC patients. After dividing into high- and low- m6A score groups, the survival outcomes were significantly different between the two groups (Figure5B). Meanwhile, either in the m6A cluster or geneCluster, the overall m6A scores were dramatically different among the subgroups (Figure 5C, D). Additionally, the m6A score can accurately reflect the clinical traits in terms of the survival status, T-stage as well as clinical tumor stage (Figure5E~J). Taken together, the m6A score system was able to distinguish distinct m6A isoforms properly and associated with the prognosis and the clinical traits in HCC.

Mounting evidence indicates that oncogenic somatic mutations relate to survival expectancy and the efficacy of immunotherapy. Therefore, we divided the patients from TCGA based on the best cutoff and performed a survival analysis between the two groups. As expected, the survival outcomes in the high-TMB group deteriorated with a p value <0.001 (Figure 5K). Subsequently, we integrated the m6A score and TMB to predict the survival status, by which the patients were classified into four subgroups. Astonishingly, the patients in the H-TMB+H-m6Ascore group suffered the worst survival expectations, only with a 2-year life expectancy. However, patients in the H-TMB+L-m6Ascore and L-TMB+H-m6Ascore groups appeared to have similar survival expectations. Inevitably, patients in the L-TMB+L-m6Ascore group had the best survival outcomes (Figure 5M). Moreover, we also revealed that TMB has a negative correlation with the m6A score in HCC (R=-0.13, P=0.016) by Pearson correlation analysis (Figure 5L). Additional findings indicated that the mutation frequency was higher in the low-m6A group, which is consistent with the results of gastric and colorectal cancers (Figure 5N, O). In further details, the significant mutation gene (SMG) mutation landscape revealed that CTNNB1 (32% vs. 18%) and TTN (25% vs. 22%) had higher somatic mutation rates in the low-m6A group, while TP53 (20% vs. 36%) had higher somatic mutation rates in the high-m6A subtype. The immune checkpoint biomarker PD-L1 was selected and significantly enriched in the high-m6A score group (Figure 5P). Moreover, to validate the pluripotency and stability of the m6A score in predicting survival, we categorized the patients into the early stage group (T1 and T2) and the advanced stage group (T3 and T4). The survival curves of the two groups separately revealed that patients with a high-m6A score appeared to have a shorter life expectancy regardless of their tumor grade (Figure 5Q, R).

The role of m6A score system in predicting immunotherapeutic and chemotherapeutic benefits

In Figure 6A, we demonstrated the correlation between m6A scores and tumor-infiltrating cells by six methods of TME calculation (Figure 6A). We next evaluated the correlation between the m6A score and 23 types of immune cells by ssGSEA (Figure 6B). Both analyses revealed that the vast majority of immune cells, as well as stromal cells, were significantly and positively correlated with the m6A score, which was consistent with the previous results of single-cell sequencing and geneClusters. In immunotherapeutic analysis, IPS was significantly higher in the low-m6A score group regardless of whether the patient has a positive response to CLTA4 or PD1 immunotherapy (Figure 6 C~F). Nowadays, chemotherapy remains an essential component of first- and second-line treatment options in individualized comprehensive treatment for HCC patients. Therefore, the analysis between chemotherapeutic agents and the m6A score indicated that patients with higher m6A scores appeared to have lower sensitivity to chemotherapeutic drugs (Figure6 G~O). In summary, the m6A regulators in HCC were deeply involved in sculpting the TME, and the scoring system can sensibly provide instructions on the clinical management for the medics.

The role of RNA m6A modification has been increasingly recognized in various biological activities, such as stem cell development, immune response, and malignant progression, through affecting gene expression (36–40). Unexceptionally, evidence accumulated confirms that m6A modification involves every stage of HCC, ranging from tumor initiation to tumor growth and, finally, to metastatic dissemination (41–43). Our works further confirmed the findings of mounting studies concluding the reliable role of m6A genes in predicting survival outcomes of cancer patients (44–47).

In the present study, we integrated the results of K-M analysis and the impact of m6A modification on immune microenvironment. In HCC, risky m6A regulators were positively associated with the immune infiltration. Conversely, protective m6A regulators were negatively correlated to the immune infiltration. When comprehensively considering the whole of m6A regulators at single-cell resolution, we found that m6A regulators were remarkably appearing in the stromal compartments, particularly in fibroblasts and endothelial cells (Fig. 3E, F). We further noticed the striking activation of stroma-related signaling pathways among the m6A clusters (Fig. 3A-C). Additionally, all the kernel genes involved in the EMT, TGF-β, and Wnt signaling pathways were activated in geneCluster B. Although the degree of immune signaling pathways and immune cell infiltration was also highest in group B, its patients had the worst survival and clinical traits. Moreover, six distinct measurements were conducted to leverage the TME. And the m6A score was positively correlated with stromal cellular compartments (Fig. 6A). Therefore, we speculated that stromal cellular compartments critically modulate the evolution and reshape the TME through m6A modification in HCC. Indeed, previous studies had demonstrated that the m6A modification in medicating the TME during cancer evolution (48, 49). Additionally, the genomic foundations confirmed distinct m6A modification patterns could result in corresponding categories of TME and discrepant evolutionary outcomes (50, 51). Furthermore, m6A modification can also be used to evaluate clinical stage, survival outcomes, and individual treatment (50, 52).

However, we identified three distinct immune microenvironment landscapes of HCC according to different m6A modification patterns. As expected, the stromal microenvironment landscapes of patients were regarded as factors of tumor deterioration consistent with previous reports (53, 54). GeneCluster A with the lowest immune microenvironment can be defined as immune-desert, while the stromal infiltration is poor and the patients had the best prognosis. GeneClusterC was characterized as an immune-inflamed phenotype, which possessed a moderate infiltration of both immune and stromal cells. Surprisingly, patients in genecluster B, the immune-excluded phenotype, owned the highest immune cell infiltration and immune signaling activation, but with the worst survival outcome (Fig. 4G). The unique TME of HCC promotes tumorigenesis and tumor immune tolerance under conditions of fibrosis and chronic inflammation (55–58). This can be attributed to the fact that oncogenic stromal cells isolated the effective immune cells from the cancer cells. Although a high level of immune infiltration was observed in HCC patients, the density of antitumor-related immune cells was relatively low. As a recent study reported, the density rather than the absolute number of lymphocytes significantly associates with the survival and prognostic outcomes of HCC patients (59). Moreover, it is well established that as cancer progresses, cancer cells tend to evolve mechanisms underpinning their resistance against the killing of immune cells, like impeding antigen presentation, inhibiting immune response, and leading to immune cells dysfunction and exhaustion (60–63).

PD-1 expression was regarded as an immune checkpoint to represent immunotherapeutic responsibility (64–66). We found the high m6A score group appeared to have a higher level of PD-L1 expression (Fig. 5P). Intriguingly, patients with low m6A scores are insensitive to immune therapies (Fig. 6C-F). This effect probably results from the defect of antigen-presenting dynamics. Distinguishing from the early stage, the advanced tumors tend to upregulate PD-L1 and CTLA4 expression. PD-L1 molecules expressed on progressive carcinomas can competitively bind CD80 presented on dendritic cells and, impede the costimulatory interaction between dendritic cells and CD8 + T cells, leaving the dysfunction of CD8 + T cells (67, 68). However, the dysfunction of CD8 + T cell regulated by m6A modification in HCC remains to be validated in cell culture systems and animal models.

Taken together, in the present study, we identified m6A regulators that were integrated to shape the evolutionary course of HCC through mediating the TME. Moreover, we constructed a m6A score system based on the differential expression of m6A related genes, by which HCC patients can be assigned to three groups with distinct tumor grades and survival outcomes. Finally, we estimated that this m6A scoring system is reliable in predicting the responsibility of HCC to both immunotherapies and chemotherapies. Our work reconfirmed the essential role of m6A modification in the development of HCC and provided novel insight when evaluating the tumor aggressiveness and survival information of the disease.

Conflict of Interest

The authors declare no potential conflicts of interest.

Sung, H. et al. (2021) Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries.CA Cancer J Clin
Forner, A., Reig, M. & Bruix, J. Hepatocellular carcinoma. The Lancet, 391, 1301–1314 (2018).
Villanueva, A. Hepatocellular Carcinoma. N Engl J Med, 380, 1450–1462 (2019).
Jemal, A. et al. (2017) Annual Report to the Nation on the Status of Cancer, 1975–2014, Featuring Survival. J Natl Cancer Inst 109
von Felden, J., Garcia-Lezana, T., Schulze, K., Losic, B. & Villanueva, A. Liquid biopsy in the clinical management of hepatocellular carcinoma., 69, 2025–2034 (2020).
Greten, T. F., Lai, C. W., Li, G. & Staveley-O'Carroll, K. F. Targeted and Immune-Based Therapies for Hepatocellular Carcinoma., 156, 510–524 (2019).
Cariani, E. & Missale, G. Immune landscape of hepatocellular carcinoma microenvironment: Implications for prognosis and therapeutic applications. Liver Int, 39, 1608–1621 (2019).
Locy, H. et al. Immunomodulation of the Tumor Microenvironment: Turn Foe Into Friend. Front Immunol, 9, 2909 (2018).
El Dika, I., Khalil, D. N. & Abou-Alfa, G. K. Immune checkpoint inhibitors for hepatocellular carcinoma., 125, 3312–3319 (2019).
Khemlina, G., Ikeda, S. & Kurzrock, R. The biology of Hepatocellular carcinoma: implications for genomic and immune therapies. Mol Cancer, 16, 149 (2017).
Llovet, J. M., Montal, R., Sia, D. & Finn, R. S. Molecular therapies and precision medicine for hepatocellular carcinoma. Nat Rev Clin Oncol, 15, 599–616 (2018).
Iñarrairaegui, M., Melero, I. & Sangro, B. Immunotherapy of Hepatocellular Carcinoma: Facts and Hopes. Clin Cancer Res, 24, 1518–1524 (2018).
Dominissini, D. et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-sEq. Nature, 485, 201–206 (2012).
Liu, H. et al. A novel N6-methyladenosine (m6A)-dependent fate decision for the lncRNA THOR. Cell Death Dis, 11, 613 (2020).
Tang, Y. et al. m6A-Atlas: a comprehensive knowledgebase for unraveling the N6-methyladenosine (m6A) epitranscriptome. Nucleic Acids Res, 49, D134–d143 (2021).
Zhang, H. et al. Dynamic landscape and evolution of m6A methylation in human. Nucleic Acids Res, 48, 6251–6264 (2020).
Reichel, M., Köster, T. & Staiger, D. Marking RNA: m6A writers, readers, and functions in Arabidopsis. J Mol Cell Biol, 11, 899–910 (2019).
An, S. et al. Integrative network analysis identifies cell-specific trans regulators of m6A. Nucleic Acids Res, 48, 1715–1729 (2020).
Zaccara, S., Ries, R. J. & Jaffrey, S. R. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol, 20, 608–624 (2019).
Qin, Y. et al. Role of m6A RNA methylation in cardiovascular disease (Review). Int J Mol Med, 46, 1958–1972 (2020).
Chen, M. & Wong, C. M. The emerging roles of N6-methyladenosine (m6A) deregulation in liver carcinogenesis. Mol Cancer, 19, 44 (2020).
Widagdo, J. & Anggono, V. The m6A-epitranscriptomic signature in neurobiology: from neurodevelopment to brain plasticity. J Neurochem, 147, 137–152 (2018).
Li, J. et al. m6A demethylase FTO promotes hepatocellular carcinoma tumorigenesis via mediating PKM2 demethylation. Am J Transl Res, 11, 6084–6092 (2019).
Liu, X. et al. YTHDF1 Facilitates the Progression of Hepatocellular Carcinoma by Promoting FZD5 mRNA Translation in an m6A-Dependent Manner. Mol Ther Nucleic Acids, 22, 750–765 (2020).
Lan, T. et al. KIAA1429 contributes to liver cancer progression through N6-methyladenosine-dependent post-transcriptional modification of GATA3. Mol Cancer, 18, 186 (2019).
Zhou, T. et al. m6A RNA methylation-mediated HNF3γ reduction renders hepatocellular carcinoma dedifferentiation and sorafenib resistance. Signal Transduct Target Ther, 5, 296 (2020).
Xu, J. et al. N(6)-methyladenosine-modified CircRNA-SORE sustains sorafenib resistance in hepatocellular carcinoma by regulating β-catenin signaling. Mol Cancer, 19, 163 (2020).
Zuo, X. et al. M6A-mediated upregulation of LINC00958 increases lipogenesis and acts as a nanotherapeutic target in hepatocellular carcinoma. J Hematol Oncol, 13, 5 (2020).
Lin, Z. et al. RNA m(6) A methylation regulates sorafenib resistance in liver cancer through FOXO3-mediated autophagy. Embo j, 39, e103181 (2020).
Qu, N. et al. Multiple m(6)A RNA methylation modulators promote the malignant progression of hepatocellular carcinoma and affect its clinical prognosis. BMC Cancer, 20, 165 (2020).
Zhao, Z. et al. The Effect of m6A Methylation Regulatory Factors on the Malignant Progression and Clinical Prognosis of Hepatocellular Carcinoma. Front Oncol, 10, 1435 (2020).
Zhang, X. et al. circRNA_104075 stimulates YAP-dependent tumorigenesis through the regulation of HNF4a and may serve as a diagnostic marker in hepatocellular carcinoma. Cell Death Dis, 9, 1091 (2018).
Wang, H. et al. Mettl3-mediated mRNA m(6)A methylation promotes dendritic cell activation. Nat Commun, 10, 1898 (2019).
Han, D. et al. Anti-tumour immunity controlled through mRNA m(6)A methylation and YTHDF1 in dendritic cells. Nature, 566, 270–274 (2019).
Yang, S. et al. m(6)A mRNA demethylase FTO regulates melanoma tumorigenicity and response to anti-PD-1 blockade. Nat Commun, 10, 2782 (2019).
Tan, L. et al. Positive natural selection of N6-methyladenosine on the RNAs of processed pseudogenes. Genome Biol, 22, 180 (2021).
Fang, Y. et al. (2021) The H3K36me2 methyltransferase NSD1 modulates H3K27ac at active enhancers to safeguard gene expression.Nucleic Acids Res
Tong, J. et al. (2021) Pooled CRISPR screening identifies m(6)A as a positive regulator of macrophage activation.Sci Adv7
Zhou, J. et al. (2021) m(6)A demethylase ALKBH5 controls CD4(+) T cell pathogenicity and promotes autoimmunity.Sci Adv7
Hou, Y. et al. YTHDC1-mediated augmentation of miR-30d in repressing pancreatic tumorigenesis via attenuation of RUNX1-induced transcriptional activation of Warburg effect (Cell Death Differ, 2021).
Liu, X. et al. SIRT1 Regulates N(6) -Methyladenosine RNA Modification in Hepatocarcinogenesis by Inducing RANBP2-Dependent FTO SUMOylation., 72, 2029–2050 (2020).
Gao, Y. et al. m(6)A Modification Prevents Formation of Endogenous Double-Stranded RNAs and Deleterious Innate Immune Responses during Hematopoietic Development., 52, 1007–10211008 (2020).
Su, T. et al. (2021) Insufficient Radiofrequency Ablation Promotes Hepatocellular Carcinoma Metastasis through m(6) A mRNA Methylation Dependent Mechanism. Hepatology
Wang, Q. et al. METTL3-mediated m(6)A modification of HDGF mRNA promotes gastric cancer progression and has prognostic significance., 69, 1193–1205 (2020).
Anita, R., Paramasivam, A., Priyadharsini, J. V. & Chitra, S. The m6A readers YTHDF1 and YTHDF3 aberrations associated with metastasis and predict poor prognosis in breast cancer patients. Am J Cancer Res, 10, 2546–2554 (2020).
Li, T. et al. METTL3 facilitates tumor progression via an m(6)A-IGF2BP2-dependent mechanism in colorectal carcinoma. Mol Cancer, 18, 112 (2019).
Li, Z., Li, F., Peng, Y., Fang, J. & Zhou, J. Identification of three m6A-related mRNAs signature and risk score for the prognostication of hepatocellular carcinoma. Cancer Med, 9, 1877–1889 (2020).
Jin, Y. et al. Analysis of m6A-Related Signatures in the Tumor Immune Microenvironment and Identification of Clinical Prognostic Regulators in Adrenocortical Carcinoma. Front Immunol, 12, 637933 (2021).
Shen, S. et al. N6-methyladenosine (m6A)-mediated messenger RNA signatures and the tumor immune microenvironment can predict the prognosis of hepatocellular carcinoma. Ann Transl Med, 9, 59 (2021).
Wang, T., Kong, S., Tao, M. & Ju, S. The potential role of RNA N6-methyladenosine in Cancer progression. Mol Cancer, 19, 88 (2020).
Lan, Q. et al. The Critical Role of RNA m(6)A Methylation in Cancer. Cancer Res, 79, 1285–1292 (2019).
Sun, T., Wu, R. & Ming, L. The role of m6A RNA methylation in cancer. Biomed Pharmacother, 112, 108613 (2019).
Zhang, B. et al. m(6)A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer, 19, 53 (2020).
Chong, W. et al. m(6)A regulator-based methylation modification patterns characterized by distinct tumor microenvironment immune profiles in colon cancer. Theranostics, 11, 2201–2217 (2021).
Affo, S., Yu, L. X. & Schwabe, R. F. The Role of Cancer-Associated Fibroblasts and Fibrosis in Liver Cancer. Annu Rev Pathol, 12, 153–186 (2017).
Xiong, S. et al. Cancer-associated fibroblasts promote stem cell-like properties of hepatocellular carcinoma cells through IL-6/STAT3/Notch signaling. Am J Cancer Res, 8, 302–316 (2018).
Deng, Y. et al. Hepatic carcinoma-associated fibroblasts enhance immune suppression by facilitating the generation of myeloid-derived suppressor cells., 36, 1090–1101 (2017).
Luo, Q. et al. FOXQ1/NDRG1 axis exacerbates hepatocellular carcinoma initiation via enhancing crosstalk between fibroblasts and tumor cells. Cancer Lett, 417, 21–34 (2018).
Park, J. et al. Impact of cell density in lymphocyte-rich areas in the tumor microenvironment on prognosis and gene expression landscape in hepatocellular carcinoma. Journal of Clinical Oncology, 39, 4107–4107 (2021).
Takeda, K. et al. IFN-γ is required for cytotoxic T cell-dependent cancer genome immunoediting. Nat Commun, 8, 14607 (2017).
Shi, J. Y. et al. Margin-infiltrating CD20(+) B cells display an atypical memory phenotype and correlate with favorable prognosis in hepatocellular carcinoma. Clin Cancer Res, 19, 5994–6005 (2013).
Shao, Y. et al. Regulatory B cells accelerate hepatocellular carcinoma progression via CD40/CD154 signaling pathway. Cancer Lett, 355, 264–272 (2014).
Sung, P. S. & Jang, J. W. (2018) Natural Killer Cell Dysfunction in Hepatocellular Carcinoma: Pathogenesis and Clinical Implications.Int J Mol Sci19
Li, H. et al. Programmed cell death-1 (PD-1) checkpoint blockade in combination with a mammalian target of rapamycin inhibitor restrains hepatocellular carcinoma growth induced by hepatoma cell-intrinsic PD-1., 66, 1920–1933 (2017).
Shigeta, K. et al. Dual Programmed Death Receptor-1 and Vascular Endothelial Growth Factor Receptor-2 Blockade Promotes Vascular Normalization and Enhances Antitumor Immune Responses in Hepatocellular Carcinoma., 71, 1247–1261 (2020).
Chen, D. S. & Mellman, I. Elements of cancer immunity and the cancer-immune set point. Nature, 541, 321–330 (2017).
Cheng, Y. et al. Cancer-associated fibroblasts induce PDL1 + neutrophils through the IL6-STAT3 pathway that foster immune suppression in hepatocellular carcinoma. Cell Death Dis, 9, 422 (2018).
Wei, Y. et al. IL-17A secreted from lymphatic endothelial cells promotes tumorigenesis by upregulation of PD-L1 in hepatoma stem cells. J Hepatol, 71, 1206–1215 (2019).

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

M6a Regulator-Associated Methylation Modification Patterns Shape Tumor Microenvironment Characteristics in Hepatocellular Carcinoma

Status:

Version 1

Abstract

Figures

Introduction

Methods

Data collection and panorama of m6A regulators in HCC

Prognostic analysis and unsupervised clustering of m6A regulators

Gene set variation analysis (GSVA) and KEGG annotation

Integrated Analysis of m6A clusters and TME in HCC

Distinctly expressed genes among m6A clusters

Identification of prognosis-related genes and m6A geneCluster

Multi-omics analysis between m6A score and TME, TMB, and clinical characteristics

Correlation analysis of m6A score with TME and m6A efficiency on immunotherapy and chemotherapy

Results

Landscape of genetic variation of m6A regulators in HCC

Identification of m6A methylation modification patterns in HCC mediated by the 23 regulators

The m6A clusters mediates distinct immune landscapes

m6A phenotype-regulated DEGs and geneClusters construction in HCC

Construction of the m6A score system and exploration of its clinical relevance

The role of m6A score system in predicting immunotherapeutic and chemotherapeutic benefits

Discussion

Conclusion

Declarations

Conflict of Interest

References

Additional Declarations

Supplementary Files

Status:

Version 1