Landscape of genetic and expression variation of m6A regulators in hepatocellular carcinoma.
We summarized mutations and CNVs (Copy number variations) to explore genetic characteristics of m6A regulators. Among 23 m6A regulators, 19 showed somatic mutation in TCGA-LIHC cohort (Fig. 1A). KIAA1429 had the highest mutation frequency (1.4%), followed by ZC3H13(1.1%), YTHDC2(1.1%) and HNRNPC (1.1%). METTL3, METTL14, YTHDF2 and ALKBH5 did not show any mutation in database. Figure 1B showed a prevalent CNV alteration in 23 regulators. In addition, KIAA1429, YTHDF3 and IGF2BP2 had the highest frequency of CNV amplification (over 45%). Considering that m6A regulators play a vital role in cancer development,[15, 46, 47] we compared their mRNA expression level in tumor and normal tissue. Of the 23 m6A regulators, 20 were significantly upregulated, while only 3 genes (ALKBH3, IGF2BP1 and YTHDC1) were downregulated in HCC (Fig. 1C). Based on the differential expression characteristics, tumor and normal samples were well distinguished based on the expression profiles of m6A regulators in principal component analysis (Fig. 1D). Furthermore, we investigated the association between the expression of m6A regulators and overall survival. Ten regulators were associated with worse overall survival, while five regulators were associated with better overall survival. (Fig. 1E).
Correlation among 23 m6A regulators and identification of m6Aclusters.
We obtained protein-protein interaction network from STRING,[38] and found complicated correlation among 23 m6A regulators (Fig. 2A). Then, we performed spearman correlation analysis based on mRNA expression and found the majority of m6A regulators significantly correlated with each other (Fig. 2B, Table S1). The GO analysis demonstrated that function of m6A regulators enriched in DNA and RNA modification (Figure S1A. These processes included the regulation of alternative mRNA, RNA stabilization, DNA dealkylation, DNA repair, oxidative demethylation. Importantly ,the m6A regulators not only participated in the same biological process, but also regulated each other. KIAA1429 (“writer”), showed the strongest correlation with YTHDF3 (“reader”) (Spearman r = 0.62, P < 0.0001), indicating important cross-talk among different functional categories (“writers”, “readers”, “erasers”).
Using consensus clustering analysis based on expression of 23 m6A regulators of 699 HCC patients, we identified 328 patients in cluster A and 341 patients in cluster B (Fig. 2C, Figure S1B-S1D). The expression patterns of 23 m6A regulators were significantly different between clusters (Fig. 2D). The m6Acluster A was prominent with higher expression of ELAVL1, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, KIAA1429, METTL3, RBM15B and YTHDF1 (Figure S1E). Patients in m6Acluster B had a better overall survival compared to those in m6Acluster A (P = 0.016; HR = 1.43 (1.07–1.90)) (Fig. 2E).
GSEA analysis shown that differential pathways were all enriched in metabolic processes (Figure S1F), including glucose metabolism, fatty acid metabolism, retinol metabolism and ABC transporters. In addition, pathways related to drug metabolism and steroid hormone metabolism were significantly upregulated in m6Acluster B.
Pathway enrichment of DEGs and identification of m6A-gene-clusters
A total of 147 DEGs were identified between m6Acluster A and B (P < 0.05, FDR < 0.01) (Table S2). We performed GO and KEGG analysis to enrich pathways that m6A-related genes involved in. Consistent with results of GSEA between m6Acluster A and B, 147 DEGs were significantly involved in metabolism, including steroid metabolism, retinol metabolism, glucose metabolism and Tyrosine metabolism (Fig. 3A). Enrichment in growth and stem cell differentiation was also showed in Go analysis. Then we constructed PPI (protein-protein interaction) network (Figure S4A) and identified the most significant modules (Fig. 3B). In addition, 12 hub genes ranked by degrees were obtained and were shown in Figure S4B. These hub genes were also involved in metabolic pathways (Figure S4C).
Two m6A-gene-clusters were identified based on expression of 147 DEGs by Consensus clustering analysis (Fig. 3C and Fig. 3D). The prominent differences in the expression of m6A regulators between m6A-gene-clusters was in accordance with the expected results of m6A methylation modification patterns, indicating that 147 DEGs were associated with m6A modification (Fig. 3E). M6A-gene-clusterB had a better prognosis than m6A-gene-clusterB (P = 0.0007, HR = 1.49) (Fig. 3F).
Differential metabolic characteristics in distinct m6A modification patterns
We demonstrated significantly differential pathways involved in metabolism between m6A-gene-clusters (Figure S1F). M6A-gene-cluster A exhibited a higher activity of hypoxia and glycolysis than m6A-gene-cluster B. The m6A-gene-cluster B was characteristic with higher activity of tricarboxylic acid cycle (Fig. 4A, 4B). Expression of FBP1, a key inhibitor of glycolysis,[48] was significantly lower in m6A-gene-cluster A (Figure S3A). The correlation between median glycolytic expression and DEGs was showed in Table S3. MAPK13 was most strongly positively correlated with glycolysis (Spearman r = 0.44, P < 0.0001) and AQP9 was most strongly negatively correlated with glycolysis (Spearman r=-0.45, P < 0.0001) (Fig. 4C).
Expression of cholesterol metabolism related genes was higher in m6A-gene-cluster A (Fig. 4D). In addition, SOAT1 and SREBF2, two key regulators of cholesterol metabolism, were also significantly upregulated in m6A-gene-cluster A (Fig. 4E).M6A-gene-cluster A exhibited a higher expression of fatty acid biosynthetic process and a lower expression of fatty acid oxidation (Fig. 4F, 4G), causing fatty acid accumulation and promoting proliferation and migration of HCC cell.
In conclusion, m6A-gene-cluster A was characterized with higher activity of hypoxia, glycolysis, cholesterol metabolism, and fatty acid biosynthesis, while m6A-gene-cluster B exhibited higher activity of TCA (tricarboxylic acid cycle), fatty acid oxidation, steroid hormone metabolism (Figure S3B) and retinol metabolism (Figure S3C).
Differential immune characteristics in distinct m6A modification patterns
The m6A-gene-cluster A were remarkable with abundant B cell, CD8 T, Exhausted T, nTreg, Treg1 and higher expression of all immune checkpoints (Fig. 5A, 5B). Expression of most immune related genes was higher in m6A-gene-cluster A (Fig. 5C). However, expression of TGFβ was also higher in m6A-gene-cluster A (Fig. 5D), which was consistent with higher Treg abundance and expression immune checkpoints. Furthermore, ESTIMATE score and immune score of m6A-gene-cluster A were also higher than m6A-gene-cluster B (Fig. 5E, 5F). Therefore, m6A-gene-cluster A was characterized as immune-activated and immune-suppressive simultaneously.
In correlation analysis between immune score and DEGs (Table S4), SAA1 (Serum Amyloid A), a protein associated with Amyloidosis, was found most positively correlated with immune score (Spearman r = 0.36, P < 0.0001) (Fig. 5G). In addition, HCCs with high expression of SAA1 exhibited higher immune score, stromal score, and Estimate score, indicating abundant immune infiltration and lower tumor purity (Fig. 5H, 5I). A better overall survival was showed in SAA1-High group (Fig. 5J). Above all, SAA1 might be a potential target for immune therapy.
Characteristics of m6A score in prognosis, immune microenvironment and metabolism
In univariate analysis, prognostic factors associated with m6A related genes from DEGs was showed in Table S4. We constructed m6A score using principal component analysis in a cohort of 669 HCC patients. Survival analysis demonstrated patients with high-m6A score (Median m6A score = 0.11) had a better OS (P = 0.0003, HR = 1.69) (Fig. 6A). As expected, m6A score in m6Acluster B and m6A-gene-clusterB was significantly higher than m6Acluster A and m6A-gene-clusterA (Fig. 6B, 6C). The m6A score was negatively correlated with TNM stage (AJCC, 2010) (Fig. 6D). In addition, HCCs with TP53 mutation and vascular invasion also had lower m6A score (Fig. 6E). In validated group of 225 HCC patients in GSE14520 cohort, high-m6A score was associated with better recurrence-free survival (RFS) (Fig. 6F) and overall survival (Fig. 6G). M6A score was associated with tumor stages in various staging systems, including TNM stage (AJCC, 2010), BCLC stage and CLIP stage in training cohort and validation cohort. HCCs with larger tumor size (> 5cm) and multiple nodules also had lower m6A score (Fig. 6H).
In correlation between m6A score and tumor immune microenvironment, the low-m6A score group was characteristic with higher abundance of B cell, CD8 + T, NK, Th1, Exhausted, Treg1, nTreg, iTreg. However, the abundance of Gamma-delta T, MAIT, Monocyte, Neutrophil, Tfh and Th17 were higher in the high-m6A score group (Fig. 6I). TGFβ was significantly upregulated, indicating immune-inhibition in low-m6A score group (Fig. 6J). we also compared CYT between groups and found no differences (Fig. 6K). In addition, the m6A score was negatively correlated with expression of PDCD1 (Fig. 6L) and CTLA4 (Figure S5A). The stromal score was also lower in low-m6A score group (Figure S5B). The m6A score was positively correlated with Th17 and Tfh abundance, and negatively correlated with nTreg and B cell (Figure S5C).
In addition to the immune cell infiltration, we compared expression of glucose metabolism related genes between groups. Consistent with results of GSEA, GO and KEGG analysis, low-m6A score was correlated with significantly higher activity of glycolysis and lower activity of TCA (Fig. 6M). These characteristics could promote proliferation and migration of HCC and contribute to worse prognosis.
Characteristics of m6A score in targeted therapy and immunotherapy resistance
To examine whether m6A score was associated with targeted therapy and immunotherapy resistance, we compared mRNA expression of several pathways involved in drug resistance. Hypoxia,[49] EGF-signaling,[50] FGF-signaling,[51] PI3K-AKT pathway,[52] SCD[53] were associated with sorafenib resistance. In addition, TGFβ,[54] hypoxia,[55] stemness,[56] WNT pathway,[57] ENTPD1 and NT5E[58] were associated with immunotherapy resistance. In low-m6A score group, pathways involved in sorafenib resistance (hypoxia, EGF-signaling, FGF-signaling, MEK/ERK pathway, PI3K-AKT pathway, RTK signaling, SCD) and immunotherapy resistance (TGFβ, hypoxia, stemness, WNT pathway, ENTPD1 and NT5E) were significantly upregulated (Fig. 7A). We found differences in expression of KLF4, OCT4, MYC and SOX2 between groups (Fig. 7B). In addition, m6A score was negatively correlated with mRNAsi (Fig. 7C). Above all, m6A score was negatively correlated with stemness.
HCCs with low expression of HNRNPC, IGF2BP1, METTL3 and YTHDF1 had significantly higher m6A score, and low expression of FTO, METTL14 and ZC3H13 was associated with lower m6A score (Fig. 7D). In 29 HCC patients treated with Sorafenib from TCGA-LIHC cohort, patients with low expression of HNRNPC, IGF2BP1, METTL3 and YTHDF1, namely high m6A score, had a better overall survival (Fig. 7E). So did patients with high expression of FTO, METTL14 and ZC3H13 and high m6A score also had better overall survival (Fig. 7F).