MAPK-RAP1A Signaling Enriched in Hepatocellular Carcinoma Is Associated With Favorable Tumor-Infiltrating Immune Cells and Clinical Prognosis

Background MAPK-RAP1A signaling, which is involved in cancer progression, remains to be defined. Upregulation of MAPK-RAP1A signaling accounts for most cancers that harbor high incident rate, such as non-small cell lung cancer (NSCLC) and pancreatic cancer, especially in hepatocellular carcinoma (HCC). MAPK-RAP1A signaling plays an important function as clinical diagnosis and prognostic value in cancers, and the role of MAPK-RAP1A signaling related with immune infiltration for HCC should be elucidated. Methods Microarray data and patient cohort information from The Cancer Genome Atlas (TCGA; n = 425) and International Cancer Genome Consortium (ICGC; n = 405) were selected for validation. The Cox regression and least absolute shrinkage and selection operator (LASSO) were used to construct a clinical prognostic model in this analysis and validation study. We also tested the area under the curve (AUC) of the risk signature that could reflect the status of predictive power by determining model. MAPK-RAP1A signaling is also associated with tumor-infiltrating immune cells (TICs) as well as clinical parameters in HCC. The GSEA and CIBERSORT were used to calculate the proportion of TICs, which should be beneficial for the clinical characteristics (clinical stage, distant metastasis) and positively correlated with the survival of HCC patients. Results HCC patients with enrichment of MAPK-RAP1A signaling were associated with clinical characteristics and favorable T cell gamma delta (Vδ T cells), and STMN1, RAP1A, FLT3, HSPA8, ANGPT2, and PGF were used as candidate biomarkers for risk scores of HCC. To determine the molecular mechanism of this signature gene association, Gene Set Enrichment Analysis (GSEA) was proposed. Cytokine–cytokine receptor interaction, TGF-β signaling pathway, and Intestinal immune network for IgA production gene sets were closely related in MAPK-RAP1A gene sets. Thus, we established a novel prognostic prediction of HCC to deepen learning of MAPK-RAP1A signaling pathways. Conclusion Our findings demonstrated that HCC patients with enrichment of MAPK-RAP1A signaling were associated with clinical characteristics and favorable T cell gamma delta (Vδ T cells), which may be a novel prognostic prediction of HCC.


INTRODUCTION
Hepatocellular carcinoma (HCC) is the sixth leading cancer and accounts for the second highest number of lethal malignant cancer globally (1). HCC accounts for most primary malignancy of liver cancers and is biologically heterogeneous. Despite treatment advances in surgery, liver transplantation, and radiofrequency ablation, the 5-year survival rate for HCC is still less due to heterogeneous tumor (2). At present, the development of sorafenib represented the main treatment for advanced liver cancer (3). However, a larger proportion of patients will have recurrence or metastasis (4). Thus, it is worth identify that distinct molecular genes and pathways, which accelerate molecular studies with the prognostic outcomes of HCC (5).
Among them, the MAPK and RAP1A are the mediator pathways in regulating cell biology of cancer cells (6,7). MAPKs are best-characterized regulators of extracellular signals that are transduced to the nucleus; MAPKs include three important protein kinase families: extracellular signalregulated protein kinases (ERKs), c-Jun NH2-terminal kinases (JNKs), and p38 family of kinases (8). Ras-associated protein-1 A (RAP1A), a small GTPase that belongs to the Ras-related protein family, was observed to regulate oncogenic Ras phenotype and Rho GTPase mediated actin cytoskeleton. RAP1A signaling is involved in cell proliferation, differentiation, and cell-cell junction in several cancer types (9)(10)(11). It was also reported that RAP1A regulates cancer through activating MAPK signaling in order to regulate motility and metastasis (12). The upregulation of RAP1A induced the MAPK signaling, which indicated that p38 may be identified as the counterparts of the RAP1A-dependent biological function (7,12). Therefore, the clinical involved in the expression of MAPK-RAP1A signaling should be strongly considered.
The tumor microenvironment (TME), which contains stromal and immune cells, is thought to be a determinant factor in the advance of HCC (13,14). Immune stroma plays a major role in the TME and development of tumor growth to the metastasis of HCC (15). Meanwhile, recent evidence indicates that immune cells may also play a role in cancer via ERK, another pathway of MAPK, thus probably essential for the development of therapeutic processes (16). Pancreatic cancer cells have been reported to be taken up by T lymphocytes to activate p38 MAPK via secreting exosomes, which ultimately causes immunosuppression (17). HCC-derived exosomes carried HMGB1 and activated TIM-1 + Breg cell; this delivery led to increase in angiogenesis by TLR and MAPK signaling pathways (18). To our knowledge, Ubc9 acts as a functional binding partner of ADAP and plays a selective role in TCR-mediated Rac1 activation via modulation of the membrane targeting of RAP1A and RapL (19). Therefore, the construction of MAPK-RAP1A signaling is known to involve the tumor microenvironment that can improve HCC prognosis.
In the present study, we constructed tumor-infiltrating immune cells (TICs) closely related to MAPK-RAP1A signaling through bioinformatics methods. Next, MAPK-RAP1A related signature genes significantly related to clinical features were detected to validate diagnosis and prognosis. Then we identified TICs by integrating MAPK-RAP1A related signature genes for HCC. We here aimed to provide potential marker for assessment of HCC clinical prognosis and playing a key role in TME.

HCC Patients
From January 2020 to December 2020, 20 human liver cancer tissue samples and corresponding para-carcinoma tissue samples were collected from patients with HCC in Hongqi Hospital Affiliated to Mudanjiang Medical University. The sample collection was approved by the research ethics committee of Hongqi Hospital Affiliated to Mudanjiang Medical University. Each patient provided informed consent to participate in the study. All collected samples were immediately frozen in liquid nitrogen until subsequent analysis.

TCGA and ICGC Data Source
All eligible sequencing datasets, clinical characteristics, and follow-up information in The Cancer Genome Atlas (TCGA) (https://gdc.cancer.gov/about-data/publications/pancanatlas) and International Cancer Genome Consortium (ICGC) database (https://icgc.org/) were downloaded. Gene sets with HCC were used to perform KEGG analyses. The KEGG pathway in the TCGA database was used as the training set, and data in the ICGC data sets were used for the validation set. Meanwhile, MAPK-RAP1A related genes were selected using the R package "limma" from the pre-processed data with log2(x + 1) transformation, and the adjusted P value <0.05 was considered statistically significant.

Construction of a Risk Signature Associated With MAPK-RAP1A Signaling
To screen the association between clinical stage and risk signature, multivariate cox regression models were used to assess the clinical characteristics of MAPK-RAP1A signaling in R language. The data in the TCGA database were used as the training sets, and data in the ICGC data sets were used for the validation sets to establish a nomogram. The HCC patients' information identified the relationship between clinical stage and risk scores, and Wilcoxon rank test as the significance test. LASSO regression was used to set precision power gene combination to identify key modules (20). R language survival and survminer were applied for the survival analysis and timedependent receiver operating characteristic (ROC) curve to assess the efficiency of the risk signature.

Association of Hub Genes' Expression With Tumor Purity and Tumor-Infiltrating Immune Cells
The ssGSEA and CIBERSORT tools were used to estimate tumor-infiltrating immune cells (TICs) for sample by R language (21,22); the final estimate outcome was sum up to three kinds of scores: stromal score, immune score, and ESTIMATE Score, which correlated with the larger ratio of the immune infiltration. Moreover, the infiltration of T gamma delta showed a relatively lower risk signature in MAPK-RAP1A signaling. Next, we applied Tumor Immune Estimation Resource (TIMER) to integrate the results of those hub gene expression of HCC hub genes and both tumor purity, TICs. Interestingly, STMN1, RAP1A, FLT3, HSPA8, ANGPT2, and PGF were all positively associated with tumor purity. Furthermore, a relationship was observed between the six hub genes and infiltration of immune cells (23,24). Combined with the above results, we suggested that the prognostic evaluation based on the high abundance of MAPK-RAP1A signaling assessed through ssGSEA or CIBERSORT analysis is reliable.

Functional Enrichment Analysis
Gene set enrichment analysis (GSEA) was performed to investigate the MAPK-RAP1A signature gene potential mechanisms (21). The gene set "c2.cp.kegg.v6.2.symbols.gmt", downloaded from the Molecular Signature Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp), was selected as the reference gene set. GSEA was used to elucidate biological processes associated with the significant survival difference; R package "ClusterProfiler" was utilized for GSEA analysis of the hub genes between low expression group and high expression group in R with adjusted P <0.05 and |log 2 FC| >0.5 (25).

Immunohistochemical Staining
Each resected specimen was fixed with 10% formalin, dehydrated, and embedded in paraffin. Paraffin sections were 3 mm thick and placed on glass slides. Firstly, the sections were deparaffined and incubated in an oven at 70°C for 1 h and 30 min. Then sections were carried out in xylene, and rehydration was carried out after gradient ethanol dehydration. Afterwards, the sections were incubated in 1× EDTA solution at low strength of microwave for 15 min to keep the endogenous peroxidase activity. Then, the sections were washed with PBS for three times and were then incubated with primary antibodies against STMN1 (1:200 . Finally, these sections were counterstained with hematoxylin, dehydrated by a gradient ethanol, followed by xylene, and mounted. The staining of each specimen was evaluated through two independent investigators blinded to the clinicopathological information. When it was present in the membrane or cytoplasm, the expression of proteins was considered positive. The staining score was assessed according to parameters of intensity and extension, and scored to determine the protein expression profiles. The staining-based expression levels were divided into four categories: positive, moderate, low, and negative; the scoring system used is based on the proportion of stained cells (>75, 25-75, or <25%); and the intensity of staining was also categorized into four: strong, moderate, weak, or negative. A final combined score of 0-12 was obtained by multiplying the intensity and percentage scores. Patients were classified into high or low protein expression groups based on median expression scores.

Western Blot
The total protein was extracted from four cell lines Huh7, HepG2, Hepa1-6 in ice by RIPA lysis buffer (Beyotime, Shanghai, China) with protease and phosphatase inhibitor cocktail (Roche, Basel, Switzerland). The protein concentration was detected by BCA Protein Assay Kit (Beyotime, Shanghai, China); proteins were separated by SDS-PAGE and transferred onto polyvinylidene difluoride membranes (Millipore, Bedford, MA). Then incubated in 5% Bovine Serum Albumin at room temperature for 2 h. The membranes were exposed to primary antibodies against ERK, p-ERK, JNK, p-JNK, p38, p-p38, and RAP1A (1:1,000 dilution; Santa Cruz Biotechnology, Santa Cruz, CA) or b-actin (1:1,000 dilution; Cell Signaling Technology, USA) at 4°C overnight. After that, the membranes were incubated with secondary antibody (1:5,000 dilution; Cell Signaling Technology, USA), and protein blots were detected by Automatic Chemiluminescence Imaging Analysis System (Tanon, Shanghai, China).

Statistical Analyses
All statistical tests including Cox regression models, LASSO regression, ROC curve analysis and K-M survival analyses were conducted using Rversion 3.5.1 (https://www.r-project. org/). Statistical differences between distributions were computed by independent t test between two groups and Kruskal-Wallis for multigroup comparison. P <0.05 was considered statistically significant. The charts, forest plots, and calibration plots were drawn using R language. This study finally selected to build the diagnostic prediction model or diagnosis guidelines using a logistic regression method.

Construction of A Nomogram for Predicting Prognostic Risk of MAPK-RAP1A
In order to detect the MAPK-RAP1A pathway, we detected the expression of ERK, JNK, p38, and RAP1A by tissue RT-qPCR and found that the expression of ERK, JNK, p38, and RAP1A in HCC tumor tissues was higher than that in adjacent tissues ( Figure 1A). We also detected the expression of ERK, JNK, p38, and RAP1A by western blot and found that the expression of phosphorylated ERK, JNK, p38, and RAP1A in Huh7, HepG2, and Hepa1-6 was higher than in L02 ( Figure 1B). To better confirm the MAPK-RAP1A related genes with risk signature, we evaluated HCC patients with detailed clinical sample information (gender, age, histological grade, pathologic stage, and TNM stage) and follow-up information in the TCGA platform. Firstly, multivariate cox regression analyses were conducted to assess the independent predictive power in the training cohort, revealing that the risk signature remained as an independent risk factor correlated with HCC patients' prognosis (p < 0.05) ( Figure 1C). Second, risk signature was significantly correlated with TNM stages and overall survival ( Figures 1D, E). Subsequently, multivariate Cox regression analyses was tested in the validation cohort by the ICGC datasets, indicating that gender, stage, and risk scores were independent risk factors for HCC patients (p < 0.05) ( Figure 1F). Finally, a nomogram integrating the factors was constructed for predicting clinical features of HCC. It was noted that the nomogram model demonstrated the probability accuracy for predicting HCC ( Figure 1G). We considered that a nomogram integrating MAPK-RAP1A risk signature might act as accurate predictive power.

Six-Related Genes Were Screened Out for Constructing a Risk Signature
Through the LASSO algorithm, 328 MAPK-RAP1A related gene were performed to build the prognostic risk signature in the TCGA training cohort, and 11 genes (STMN1, RAP1A, FLT3, HSPA8, FGF9, EFNA5, IRAK1, RAC3, DUSP10, ANGPT2, and PGF) were filtered out because of shrinking parameters ( Figure  2A). Subsequently, six genes (STMN1, RAP1A, FLT3, HSPA8, ANGPT2, and PGF) were selected to establish a prognostic model ( Figure 2B). Risk score = 1.5962 * (the expression level of STMN1) + 1.7459 * (the expression level of RAP1A) + 2.2693* (the expression level of FLT3) + 1.3980 * (the expression level of HSPA8) + 1.462 * (the expression level of ANGPT2) + 1.7673 * (the expression level of PGF). Finally, a nomogram integrating the six signatures was demonstrating the performance of the risk signature in predicting HCC ( Figure 2C). The mentioned nomogram revealed ta better predictive accuracy by the total points obtained by adding up ( Figure 2D).

Stratification of TCGA Training and ICGC Validation Cohorts Using the Risk Signature
According to the median risk score, we divided the cohort patients into high-risk and low-risk groups. The high-risk group exhibited a higher frequency of poorer overall survival (OS) than the low-risk group in the TCGA platform by K-M curve ( Figure 3A). To further assess the predictive accuracy of this MAPK-RAP1A risk signature, we performed a timedependent ROC curve analysis. The AUC was 0.75, 0.805, and 0.65 at 1, 3, and 5 years for survival in the training cohort of TCGA ( Figure 3C). The ICGC dataset protocols were used for the validation of the risk signature for this purpose, which indicated that higher scores exhibited markedly worse overall survival (OS) ( Figure 3B). The AUCs for OS in the validation cohort were 0.75, 0.805, and 0.65 at 1, 3, and 5 years, respectively ( Figure 3D). Thus, we considered that the clinical prognostic model successfully stratified cohort patients from TCGA into high and low-risk groups with predicting the clinical prognosis of HCC.

Validation of Signature Genes in the MAPK-RAP1A Datasets
Among the 11 genes screened above, six signature genes (STMN1, RAP1A, FLT3, HSPA8, ANGPT2, and PGF) were identified to validate prognostic and diagnostic value and clinical features. All of them, which were significantly upregulated (p < 0.05) in HCC samples, can potentially act as oncogenes compared with normal controls. In addition, STMN1, RAP1A, FLT3, HSPA8, ANGPT2, and PGF were identified as the key oncogenic components in HCC samples with different TNM stage, T classification, lymph node metastasis, and distant metastasis, with higher expression levels indicating higher stage ( Figures 4A-D). Conversely, all six signature genes in the module might act as potential prognostic values.

The Categories of Immunity and Tumor Purity of MAPK-RAP1A Datasets
The immunity of the MAPK-RAP1A signaling was used in the principal components analysis (PCA) for assigning patients to low-and high-immunity categories (Figures 5A, B). Further comparison analysis showed that tumor purity had significantly higher expression levels in low-immunity compared with highimmunity ( Figure 5C). Regarding prognosis, K-M curves showed that higher tumor purity of MAPK-RAP1A signaling was associated significantly with worse OS ( Figure 5D).

Association of Stromal Score, Immune Score, and ESTIMATE Score With Several Clinical Features
The immune components and clinicopathological characteristics were positively correlated; we analyzed the corresponding stromal score, immune score, and estimate score of MAPK-RAP1A gene sets from TCGA-LIHC database ( Figure 6A). Among them, stromal score and estimate score are key proportion of advanced TNM stages. The representative images of the immune score were only negatively correlated with T classification of TMN stages (p = 0.007), and estimate score significantly declined with T and M classification of TMN stages (p = 0.028 and p = 0.021). Whereas there was no correlation between ESTIMATE algorithm and lymph node metastasis (P = 0.776). These results suggested that the ratio of immune components was correlated with the progress of HCC, especially in metastasis ( Figures 6B-E).

Association of Risk Signature With Tumor Purity and Tumor Immune Cell Infiltration
It was demonstrated that the levels of TICs, including Tregs, T-cells gamma delta (Vd T cells), and monocytes were significantly higher in the normal patient group compared with the HCC patient group ( Figure 7A). As shown in Figure 7B, two immune infiltration cell subpopulations were significantly enriched in the low-risk patient group (Macrophage 0 and Vd T cells). To characterize the TICs, the contents of infiltrating immune cells were investigated. We utilized TIMER website source to identified potential associations between the expression of signature genes and immune purity of tumor, infiltrating immune cells ( Figures 7C-G). Interestingly, STMN1 and ANGPT2 were significantly positively correlated with tumor purity, whereas RAP1A, FLT3, HSPA8, and PGF were all significantly negatively correlated with tumor purity. Consequently, the present study reveals relationships linking these six genes and infiltration of Vd T cells (Figures 7C-G).

GSEA Reveal a Close Relationship Between Hub Genes and Tumor Proliferation
To explore the molecular mechanisms of STMN1, RAP1A, FLT3, HSPA8, ANGPT2, and PGF in HCC, we detected GSEA on the TCGA-LIHC RNA-seq data. This work confirmed that the role of hub genes in high expression groups of RAP1A, FLT3, HSPA8, ANGPT2, and PGF were all enriched in "cytokine−cytokine receptor interaction" pathways. Additionally, the "ECM −receptor interaction" and "TGF−beta signaling pathway" gene set were enriched in high expression groups of ANGPT2 and HSPA8, whereas "Intestinal immune network for IgA production" was enriched in the PGF, FLT3, and STMNI highexpression groups, respectively, which suggest that these hub gene sets were all closely involved in tumor proliferation ( Figures 8A-F).

Correlation of the Signature With B7 Therapy-Related Molecule
Moreover, B7(CD80) served as specific markers for T cell gamma delta (Vd T), respectively. As shown in Figure 9A, the hub genes of ANGPT2, FLT3, PGF, STMN1, RAP1A, and HSPA8 were significantly positively correlated with immune checkpoint B7 (CD80) in the TCGA database (p =1.41e-11 for ANGPT2, p <2.2e-16 for FLT3, p <2.2e-16 for PGF, p <1.9e-9 for STMN1, p <2.2e-16 for RAP1A, p =2.6e-15 for HSPA8) ( Figure 9A). Subsequently, the relationships between the risk signature with MAPK-RAP1A signaling and the expression levels of B7 were found. The data here indicated that higher immune checkpoint gene of B7 may be significantly observed in high-risk patients (p = 0.012) ( Figure 9B). Interestingly, the expression levels of PGF, STMN1, ANGPT2, FLT3, RAP1A, and HSPA8 were all significant positively correlated with B7 in HCC tissues ( Figure 9C) (Supplement Figure). To further access the prognostic model between risk signature and B7 on patients' OS, a survival comparison was made among the four groups based on the combination of risk signature and immune checkpoint gene of B7. Comparison results revealed that HCC patients with risk signature and B7 in the predictive models are able to distinguish the overall survival (p <0.0001) ( Figure 9D). Finally, the present study reveals relationships linking B7(CD80) and PDL1 in HCC tissues compared with adjacent tissues ( Figure 9E).

DISCUSSION
In the current study, a novel network was developed that enabled the identification of the MAPK-RAP1A risk signature for the was assessed. For example, a recent study indicated that the immune related nine-gene signature is involved in the status of tumor immune microenvironment (26). Another study also constructed a known lncRNAs prognostic signature for cancer and may aid in the field of tumor immunity and immunotherapy (27). Studies conducted so far concluded that investigating the prognostic value of tumor-infiltrating B lymphocytes lncRNA signature (TILBlncSig) in bladder cancer identified and validated biomarkers of immunotherapy response (28). However, these studies did not use the signaling pathway of samples to comprehensively explore the relationship between TICs and prognosis of HCC.
To probe the signaling events in HCC expression, this is the first study to reveal the relationships linking MAPK-RAP1A signaling pathway and TICs. We obtained a number of prognostic MAPK-RAP1A signature genes and established a novel patient prognosis and clinical feature model. Consistently, we have linked prognostic model to poor OS, which is consistent with other studies in HCC. Consequently, the MAPK-RAP1A related prognostic model provides a more reliable tool, which was a good prognostic factor with HCC. In addition, such model that consists of six signature genes was then successfully validated as a prognostic factor. Based on this overall hypothesis, we revealed that MAPK-RAP1A prognostic model could act as a more reliable tool for HCC prognosis prediction.
In addition, we identified MAPK-RAP1A signature genes (STMN1, RAP1A, FLT3, HSPA8, ANGPT2 and PGF) which were related with HCC worse clinical phenotype and prognosis. STMN1 is an oncogene that is a highly conserved cytosolic phosphoprotein, was found to be over-expressed in various types of cancer such as lung, breast, gastric cancer, and HCC. It plays an important role in cell differentiation, proliferation, drug resistance, and cancer stem cell properties and as an emerging target for tumor therapy (29). STMN1 regulated cancerassociated fibroblast (CAF) features through HSC by triggering the HGF/MET signaling pathway (30).
RAP1A is a small G protein that is similar to Ras oncogene and associated with different cellular processes. Other studies also showed that RAP1A may be a potential target for cell proliferation, adhesion, and invasion in different types of cancers (9,31,32). Recent clinical studies reported that RAP1A also correlated with the clinical characteristics of the advanced tumor stage in Oral Cavity Squamous Cell Carcinoma (OCSCC) (33). Another study showed that RAP1A disrupts aberrant tumor suppressor of EYA4 in HCC cells, which was associated with promoting growth and invasion (34).
FLT-3 belongs to the receptor tyrosine kinase family, which is encoded by the FLT3 gene. FLT3 has been found to control the function of normal and malignant hematopoiesis (34). Activating mutations of FLT3 associated with a poor prognosis in acute myeloid leukemia (AML) and FLT3 inhibitors have an important role in high-risk patients (35). Targeting a tyrosine kinase receptor, sorafenib is now routinely required to provide therapy benefits in FLT3 of HCC patients (36,37). Moreover, the association between patients treated with sorafenib and high FLT3 levels could be a novel predictor to improve OS in HCC patients (38).
The HSPA8 protein integrates compensatory mechanisms to drive cellular growth and is dysregulated in multiple chronic stress diseases, including cancer. Mechanistically, the molecular chaperone HSPA8, also known as Hsc70c, and directed to lysosomes, selectively degrades cellular proteins to sustain cellular homeostasis (39). The role of HSPA8 is expressed abnormally in early liver cancer, and its expression increases with cancer initiation and progression (40). So, HSPA8 detection may improve early detection of liver cancer sensitive indicators, but whether and how HSPA8 is involved in cancer initiation and development has been explored.
ANGPT2 is available in the extracellular signals that play a crucial role in angiogenesis and resistance to antiangiogenic therapy (41,42). ANGPT2 is an angiogenic factor that binds Tie2 receptor and emerged as an attractive vascular drug target with the vascular endothelial growth factor (VEGF) pathway (43). Currently, ANGPT2 has been playing an important role in tumor angiogenesis and might confer resistance to sorafenib therapy (44)(45)(46). Furthermore, a study reveals that HCC cell secreted exosomal ANGPT2 was recycled by recipient HUVECs that suppressed the epithelial-mesenchymal transition (EMT) activation (47). Placenta growth factor (PGF), a member of the VEGF family, also plays an important role in the current anti-angiogenesis therapy. In this way, the present study has generated that patients with higher PGF had poorer response to chemotherapy and poorer prognosis and identified biomarkers in epithelial ovarian cancer (EOC) (48). PGF has been involved in elevated NF-kB signaling pathway in cervical cancer (49), but its role in HCC remains unclear.
To characterize the tumor-infiltrating immune cells' (TICs') status, the relationships between risk signature model and immune cell were analyzed by GSEA and CIBERSORT. The high abundance of TICs for risk signature may also affect clinical features and survival analysis. In addition, the data here indicated that higher infiltration levels of T cell gamma delta (Vd T cells) were independent prognostic protective factors in HCC. Similar to previous reports (50,51), our data confirmed that the infiltration of Vd T cells was an independent prognostic factor. Vd T cells play a critical role in the solid tumor microenvironment, but the effectiveness in killing various tumor cells present has been shown to be limited (50)(51)(52). Vd2 T cells infiltrate several types of tumors, including the liver, and could serve as a prognostic factor (53). Activated Vd2 T cells can exhibit the functions and characteristics of dendritic cells (54). On the one hand, Vd2 T cells not only express the chemokine receptor CCR7 on the surface, but also upregulate the expression levels of MHCI and MHCII molecules, as well as costimulating molecules CD80 (B7) and CD86 (55). In order to elucidate the role of the signature genes' biological functions, the association of STMN1, RAP1A, FLT3, HSPA8, ANGPT2, and PGF was positively associated with tumor purity in HCC. Based on the exploration from TIMER website, of the risk signature, six signature genes (STMN1, RAP1A, FLT3, HSPA8, ANGPT2, and PGF) were associated with infiltrating immune cells. Characterization of immune regulation by the GSEA indicated significantly different immune function among different expression groups classified by the hub gene. The results of GSEA were in accordance with this speculation. Also, many immune-related KEGG pathways, such as cytokine−cytokine receptor interaction, ECM−receptor interaction and TGF−beta signaling pathway, suggested a signature contribution to immunity regulation. In summary, taken together, all the above research suggested that the MAPK-RAP1A risk signature maybe a potential prognostic molecular marker for HCC patients. Furthermore, the predictive value of the risk signature was validated in available clinical information and RNAseq data. The MAPK-RAP1A risk signature was not only shown to have prognostic value for HCC patients but was also related to immune cell infiltration (T cell gamma delta) and the immunotherapy signature. The MAPK-RAP1A risk signature added significant predictive power to CD80 (B7) expression, which led to a significant overall survival. Altogether, these results suggested that the risk signature could shed light on the mechanisms of immunotherapeutic targeting.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the research ethics committee of Hongqi Hospital Affiliated to Mudanjiang Medical University. The patients/ participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
HL, GH, XL, BL, and LW contributed to the bioinformatics analysis. BW and HJ performed the experiment. WW reviewed the article and provided comments. All authors contributed to the article and approved the submitted version.