Hepatocellular carcinoma (HCC) is a highly aggressive primary liver malignancy and the third most common cause of cancer-related deaths worldwide. Although early HCC can be treated through surgical resection, liver transplantation, and radical ablation, the early recurrence rate after treatment is > 70% ^[1]. Systemic therapy, such as molecular-targeted agents, has made great advances in HCC treatment, but the optimal therapeutic approaches are still limited due to insidious onset and late diagnosis ^[2]. Therefore, efficient HCC treatment strategies should be developed urgently.

HCC is a typical immunogenic tumor. Infiltrating stromal and immune cells promote an immunosuppressive tumor microenvironment, allowing cancer cells to grow and evade immune monitoring. Recent studies have reported that immunosuppressive therapy can prolong the life of cancer patients ^[3]. Immune-related genes (IRGs) are the key drivers of HCC development and progression. However, approximately 20% cancer patients benefit from immunotherapy in the long term ^[3]. Therefore, identifying new immunotherapy markers and uncovering the underlying mechanisms of immune checkpoints are important.

Recent high-throughput sequencing technologies have suggested that RNA modification plays a key role in various physiological and pathological processes. N7-methylguanosine (m7G) is a common RNA modification form ^[4]. It recruits cellular proteins and mediates cap-related biological functions. Moreover, m7G modification is closely associated with tumor occurrence and development and involved in multiple tumor-related processes ^[4]. For example, the m7G regulator methyltransferase like-1 protein (MTTL1) participates in tumor immune infiltration by modifying m7G in neck squamous cell carcinoma ^[5]. However, the correlation between m7G-related genes and immunity in HCC remains unclear. Therefore, exploring the relationship between m7G-related genes and immunity in HCC and therapies targeting m7G-related genes may provide new ideas for the immunotherapeutic HCC treatment.

In this study, we obtained the HCC RNA-seq data from the UCSC Xena database (https://xenabrowser.net/datapages/), and GSE112790 datasets were downloaded from GENE Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). Then ESTIMATE, which is used to infer the gene expression data-based stromal and immune microenvironment infiltration in malignant tumor tissues using the "Estimate" R package of R 4.6.0, was applied to monitor the stromal content (stromal score), immune cell infiltration level (immune score), ESTIMATE score, and tumor purity for each TCGA-LIHC sample. The ESTIMATE score, which is defined as tumor purity, is the sum of individual stromal and immune scores. Then, we used the Spearman correlation test to analyze the correlation between m7G-related gene expression and stromal scores, immune scores, ESTIMATE scores, and tumor purity. Cluster analysis of EIF4E3 and LARP1 with immune scores was performed using the Pearson correlation distance-based consensus matrix using the R package ConsensusClusterPlus. According to the clustering results, the samples were divided into two subgroups.

Then, the tumor-infiltrating immune cells were analyzed using the CIBERSORT algorithm, a genetic analysis tool that counts 22 infiltrating immune cells in each sample. The ssGSEA method of the R software Gene Set Variation Analysis (GSVA) software package was used to analyze the infiltration level of 28 immune cells. Visual data analysis of immune checkpoints in cluster1 and cluster2 was performed using the R package.

Further study of immune infiltration in the two subgroups. The outlier genes and samples were removed using the goodSamplesGenes method of the R software package WGCNA, and the scale-free co-expression network was further constructed. We selected β = 5 (scale-free R² = 0.90) as the soft threshold, obtained eight co-expression modules, and calculated their correlation with cluster, stroma score, immunity score, ESTIMATE score, and tumor purity. Subsequently, 18 hub genes were obtained by calculating the gene significance (GS) and module membership (MM). Then, we used the protein−protein interaction (PPI) network and Kyto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis to analyze 18 hub genes. All statistical analyses were conducted using R statistical software version 4.0.2 and SPSS statistical software (version 25.0). The survival curves were estimated based on the Kaplan–Meier method and the difference between groups was determined using the log-rank test. The correlation of subtypes, genes, and immune infiltration levels was analyzed using a Pearson correlation coefficient. Hypothetical tests were two-sided, with P-values of 0.05 considered statistically significant.

m7G mRNA methylation is involved in various tumor-related biological activities and closely related to tumor occurrence and development ^[6]. Since tumorigenesis and cancer progression are also affected by the tumor immune microenvironment, we explored the relationship between m7G-related genes and immune infiltration in HCC. The analysis flowchart is shown in Figure 1.

Figure 1.  Flowchart of this study.

We assessed the stromal and immune cell fractions of by calculating four ESTIMATE indices in each sample. To explore the potential role of the m7G-related genes in HCC patient tumor immunity, we selected 32 m7G-associated genes reported in prior articles ^[7], reviews ^{[4, 7]}, and gene cards (https://www.genecards.org/) and evaluated the correlation between each m7G-related gene expression and the ESTIMATE results (Figure 2A). Considering that EIF4E3 had the highest absolute correlation with the immune score, we selected EIF4E3 for cluster analysis with NUDT11, NUDT10, and LARP1 with decreasing absolute values. Finally, EIF4E3 and LARP1 with the best clustering effect according to the consensus clustering results were selected to categorize the samples into two clusters for subsequent analysis (Figure 2B and Supplementary Figure S1A–C available in www.besjournal.com). The heatmap showed that cluster1 had high LARP1 expression and low EIF4E3 expression, whereas the opposite expression patterns were observed for the two genes in cluster2 (Supplementary Figure S1A available in www.besjournal.com). The Wilcoxon rank-sum test was performed for tumor and normal tissues using RNA-seq data from TCGA-LIHC. EIF4E3 and LARP1 expression in tumor tissues were significantly lower and higher than those in normal tissues, respectively (Figure 2C, D). EIF4E3 and LARP1 expression in HCC patients in our study was consistent with that reported in previous studies. In addition, EIF4E3 and LARP1 expression were inversely correlated with immune scores in the ESTIMATE analysis. Therefore, both EIF4E3 and LARP1 may be involved in regulating m7G modifications, which can affect immune infiltration and immunotherapy response in patients.

Figure 2.  Cluster analysis based on EIF4E3 and LARP1. (A) Correlation between m7G regulatory gene and results of ESTIMATE. (B) Patients with hepatocellular carcinoma are divided into two clusters based on EIF4E3 and LARP1. (C) Comparison of EIF4E3 expression between tumor and normal tissues. (D) Comparison of LARP1 expression between tumor and normal tissues. (E) Kaplan-Meier curves of overall survival in two clusters. (ns, no significance, ^*P < 0.05; ^**P < 0.01, ^***P < 0.001).

To characterize the differences in clinical characteristics, we performed survival curve for EIF4E3 and LARP1 in high- and low-risk groups and clusters (Figure 2E and Supplementary Figure S1) and associated analysis (Table 1). First, we analyzed EIF4E3 and LARP1 prognosis in high- and low-risk groups and clusters using Kaplan–Meier methods. The results showed that the survival curves of the high- and low-risk groups in EIF4E3 and LARP1 were not significantly different, while the survival curves of the cluster group were different, and cluster2 had a longer survival period than cluster1, which indicated that cluster2 had a better prognosis.

Therefore, we used KEGG and GSEA enrichment analyses to understand the functional differences between the two clusters and found some enriched immune-related pathways in cluster2 (Figure 3B and Supplementary Figure S2A available in www.besjournal.com). Next, we used the ESTIMATE algorithm to calculate immune scores, stromal scores, and tumor purity to understand the HCC microenvironment. The results showed that patients in cluster2 had high immune and stromal scores (Figure 3C–F), suggesting that they may have a large tumor immune microenvironment. The abundance of immune cell infiltration contributes to the tumor immune microenvironment growth ^{[8, 9]}. Endogenous T cells significantly delay malignant progression through tumor response and invasion and T-cell infiltration also increases tumor susceptibility to immunotherapy with PD1-blocking antibodies; particularly high CD8⁺ T-cell expression affects immunotherapy and survival outcomes in HCC patients ^[9]. High CD8⁺ T-cell expression can target tumor cells in the tumor microenvironment to exert anti-tumor immunity. In addition, DCs, which are key antigen-presenting cells, promote anti-tumor immunity by activating T cells. We used CIBERSORT to analyze HCC data from TCGA database and found that the proportion of CD8⁺ T-cell infiltration in cluster2 was higher than that in cluster1 (Figure 3G). Moreover, ssGSEA analysis indicated that the 24 immune cell subtypes (CD8⁺ T-cell, DC, macrophage, natural killer cell, and natural killer T-cell activation) showed significantly higher expression in cluster2 than those in cluster1 (Figure 3H). These results suggest that cluster2 may respond well to tumor immunotherapy and achieve a good prognosis. Therefore, improving tumor T-cell infiltration is one method for improving cancer immunotherapy.

Figure 3.  Comparison of immune characteristics between two clusters. (A) Differential analysis volcano map. (B) GSEA enrichment analysis. (C) ESTIMATE Score. (D) Immune Score. (E) Stromal Score. (F) Tumor Purity. (G) The proportion of immune cells in two clusters. (H) Comparison of the abundance of immune infiltrating cells in two clusters. (ns, no significance, ^*P < 0.05; ^**P < 0.01, ^***P < 0.001).

T-cell infiltration in tumors can be described as high T-cell infiltration (hot tumor) or low T-cell infiltration (cold tumor), with studies suggesting different responses to immune checkpoint inhibitor (ICI) treatment for hot and cold tumors. To identify the differences and gain thorough understanding of the immune microenvironment, we analyzed immune checkpoint expression. We found that immune checkpoints (PD1, PDL1, PDL2, CTLA4, LAG3, HAVCR2, TIGIT, CD86, BTLA, NRPQ, LAIR1, LGALS9, CH27, CD40, CD44, CD48, CD70, TNFRSF8, TNFRSF18, TNFRSF14, and VSIR), particularly CTLA4, LAG3, HAVCR2, and TIGIT, were highly expressed in cluster2 (Supplementary Figure S3 available in www.besjournal.com). ICIs are beneficial for patients with several types of cancer, including HCC, and anti-CTLA-4 immune checkpoint blockade therapy improves survival in advanced HCC patients ^[10]. In addition, LAG3, TIM3, and ITIM restore the responses of HCC-derived T cells to tumor antigens, and combine with PDL1 and CTLA4 to have additive effects. Activation of tumor-associated T cells by immune checkpoint blockade is one of the most successful immunotherapy approaches for various solid tumors. Therefore cluster2 might benefit more from ICI treatment than cluster1. EIF4E3 and LARP1 may promote cold tumor transformation into hot tumors in HCC.

Figure S3.  Comparison of immunomodulatory drug targets in HCC between two clusters. (ns, no significance, ^*P < 0.05; ^**P < 0.01, ^***P < 0.001).

To further analyze the relationship between the two clusters, we also studied the hub genes of the subtypes. We applied the DEGs (630 upregulated and 167 downregulated) to construct the WGCNA network (Supplementary Figure S4A–C). Subsequently, module-trait analysis revealed that the yellow module was significantly associated with m7G (R = 0.29, P < 0.01) and immunity (R = 0.90, P < 0.01). Then we obtained 18 hub genes (GZMA, CCL5, S1PR4, GZMH, HCST, CTSW, GZMK, NKG7, NCR3, CST7, CD2, CRTAM, CD27, GZMB, CD8B, LAG3, CD8A, TNFRSF13B) from yellow module based on MM > 0.8 and GS > 0.15 (Supplementary Figure S4D–E available in www.besjournal.com). Among these, CD27 co-stimulates OX40 and 4-1BB to promote the survival of activated T cells, which is the key to T-cell initiation and memory differentiation. Under pathological conditions, CD27 binds to tumor cell-expressed CD70. Chronic co-stimulation leads to immune checkpoint expression in T cells, thus depleting immune function. The LAG3 immune checkpoint also plays a role in immune function depletion in tumors. Therefore, we speculated that cluster2 responded better to immunotherapy than cluster1, probably due to EIF4E3- and LARP1-mediated m7G modification.

Figure S4.  Identification of immune-related modules by the WGCNA. (A) Analysis of network topology for soft powers. (B) The number of genes expressed in each module. (C) Gene dendrogram and module colors. (D) Heatmap between module eigengenes and cluster, ESTIMATE results. (E) Scatter plot of module eigengenes in the yellow module.

We developed a PPI network and performed correlation analysis to determine the interactions between proteins, which showed that the 18 hub genes were strongly correlated (Supplementary Figure S5A, B available in www.besjournal.com). Subsequently, KEGG enrichment analysis was performed to study the biological functions of the 18 hub genes (Supplementary Figure S5C). The results showed that these genes were mainly enriched in the T-cell receptor signaling pathway. The ESTIMATE and ssGSEA datasets showed the same trend (Supplementary Figure S5D and E).

Figure S5.  Analysis of 18 hub genes. (A) PPI network of hub genes. (B) Correlation between hub genes. (C) The KEGG analysis of hub genes. (D) Correlation between hub genes and results of ESTIMATE. (E) Correlation between hub genes and expression of immune cells.

We used 183 patients with GSE112790 as an external dataset to validate our approach, which further enhanced our knowledge on the association between EIF4E3/LARP1 and immunity (Supplementary Figure S6 available in www.besjournal.com). Further studies should explore the relationship between EIF4E3/LARP1 and immunity and its regulatory mechanism, which may provide a promising target for improving the immunotherapy response in HCC. This study still has some limitations. Due to the use of TCGA and GEO public databases, we could not explain EIF4E3 and LARP1 expression at the protein level or prove the direct role of EIF4E3 and LARP1 in anti-tumor immunity. Therefore, the direct mechanism needs to be further studied.

Figure S6.  GSE112790 validation of immune contexture between two clusters. (ns, no significance, ^*P < 0.05; ^**P < 0.01, ^***P < 0.001).