Identification of VEGFs-related gene signature for predicting microangiogenesis and hepatocellular carcinoma prognosis

Microangiogenesis is an important prognostic factor in various cancers, including hepatocellular carcinoma (HCC). The Vascular Endothelial Growth Factor (VEGF) has been shown to contribute to tumor angiogenesis. Recently, several studies have investigated the regulation of VEGF production by a single gene, with few researchers exploring all genes that affect VEGF production. In this study, we comprehensively analyzed all genes affecting VEGF production in HCC and developed a risk model and gene-based risk score based on VEGF production. Moreover, the model’s predictive capacity on prognosis of HCCs was verified using training and validation datasets. The developed model showed good prediction of the overall survival rate. Patients with a higher risk score experienced poor outcomes compared to those with a lower risk score. Furthermore, we identified the immunological causes of the poor prognosis of patients with high-risk scores comparing with those with low-risk scores.


INTRODUCTION
Primary hepatocellular carcinoma (HCC) is a highly invasive malignant tumor posing a significant burden on healthcare globally.Its 5-year survival rate is less than 21% [1].HCC is the main type of liver cancer, accounting for 90% of all reported liver cancer cases and the fourth leading cause of cancer deaths worldwide [2].The development of HCC is driven by several factors, among which angiogenesis is the most important biological feature considered to be a key treatment target [3].Angiogenesis is a process in which new blood vessels are formed from existing ones to build new vascular networks.Angiogenesis have been implicated in numerous physiological processes, like the tissue repair, embryonic development and inflammatory response.However, in pathological situations, such as tumors, excessive angiogenesis provides nutrients and oxygen to tumor cells thereby promoting tumor growth, invasion, and metastasis [4][5][6][7].HCC is a highly angiogenic tumor, and the level of its angiogenesis correlates with the degree of the tumor stage, grade, prognosis, and therapeutic efficacy [8].
Angiogenesis is a complex process involving multiple cellular and molecular interactions.In this process, vascular endothelial growth factor (VEGF) functions as a crucial role, which can bind to its receptor (VEGFR), thereby activating multiple signaling pathways to regulate the functions of vascular endothelial cells, such as proliferation, migration, lumen formation, and vascular permeability, etc. [9][10][11][12].VEGF expression in HCC is modulated by various factors, such as hypoxia, inflammation, tumor microenvironment, and gene mutations [13][14][15].VEGF is the dominant factor influencing HCC angiogenesis, and its level is closely associated with the degree of malignancy, metastatic ability and prognosis of HCC [3,16,17].In addition, VEGF not only regulates HCC angiogenesis, but also affects the biological characteristics of HCC cells, like the proliferation, apoptosis, invasion and metastasis [18,19].Therefore, VEGF is an important molecular marker and therapeutic target for HCC, and drugs targeting VEGF or VEGFR have shown promising potential in the clinical treatment of HCC [20][21][22][23].
To date, several prognostic prediction models based on angiogenesis-related genes have been developed for HCC, but no study has investigated this key factor of HCC angiogenesis in detail.Moreover, the role of VEGF in hepatocellular carcinoma angiogenesis is not fully understood [24][25][26][27].Moreover, considering the heterogeneity and complexity of hepatocellular carcinoma, the detection and inhibition of VEGF or VEGFR alone may not be sufficient to reflect and predict the angiogenic status and therapeutic efficacy of HCC [28].Therefore, it is necessary to explore the gene expression profiles of VEGF from a global perspective to build a gene signature that might comprehensively assess and predict microangiogenesis and prognosis of HCC patients.For our design, we constructed a gene signature capable of predicting microangiogenesis and prognosis of HCC by screening genes associated with VEGF production.Its clinical significance and biological function in HCC were validated.

Public data and collection of samples
Clinical, molecular, and whole-genome RNA-seq expression data were downloaded from TCGA database (https://portal.gdc.cancer.gov/).After exclusion of samples without prognostic information or incomplete clinical information, 374 HCC samples were included.In addition, the TCGA database was used to obtain whole-genome RNA-seq expression data of 50 normal liver samples.A total of 232 HCC patient samples were retrieved from International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/).Moreover, the GSE76427, which includes the expression and corresponding clinical data of another 114 HCC patient samples, was gained from the Gene Expression Omnibus (GEO) database (https:// www.ncbi.nlm.nih.gov/geo/).The samples were merged and standardized to create validation sets.We excluded cases with less than or exactly 10 days of survival or patients without survival data since they could have died from unrelated complications.In such cases, 50 normal liver samples and 310 HCC samples from TCGA, 231 HCC samples from ICGC, and 114 HCC samples from GSE76427, were selected for subsequent analyses.Furthermore, we considered the possibility of batch effects present among different databases as well as within the same database.To address this issue, "normalize between arrays" function [29,30] of R package "limma" was utilized to remove multiple batch influences when merging the mRNA_seq data of ICGC, TCGA and GSE76427.

Patient samples
The tissues were obtained with informed consent from all participating patients.Six control samples were collected between September 2019 and June 2022, five samples from patients with HCC and five samples from non-tumor liver tissues.Preoperative chemotherapy or radiotherapy was not used to treat any of the HCCs.Independent samples from our hospital were verified by the hub genes of GSRS mRNA levels in HCC.

Identifying differentially expressed genes between normal liver tissues and hepatocellular carcinoma
The training set was established using HCC and correlated normal liver samples from TCGA databases.DEGs were identified from Vascular endothelial growth factor Production-RGlated genes (VPRGs) using the R package "limma" [30].The criterion set for this study was a p-value less than 0.05, with an abs value of logFC greater than 0.585.

Protein-protein interaction network analysis
The STRING database (https://cn.string-db.org/) was used to construct the PPI network containing 34 VPRGs.Nodes with interaction confidences greater than 0.4 were displayed.AGING

Genomic alterations of 34 vascular endothelial growth factor production-related genes
The cBioPortal dataset (http://www.cbioportal.org/)was utilized to explore the CNV amplification, truncating mutation, in-frame mutation, missense mutation, CNV deep deletion, and fusions involving the 34 genes.

Construction of vascular endothelial growth factor production-related risk signature
To assess the prognostic value of VPRGs in HCC, we conducted univariate Cox regression analysis with the "survival" R package.Genes with a p-value < 0.05 were selected for analysis.Related risk signature was then developed for HCC by incorporating survival status, survival time, and expression levels of prognosis-related genes.The Least Absolute Shrinkage and Selection Operator (LASSO) regression algorithm [31] was utilized, with the penalty parameter λ determined through 10-fold cross-validation.By identifying the genes and their respective regression coefficients, we determined the most appropriate λ value.A formula applied for this purpose was as follows: Riskscore = exprgene (1) × coefficient gene (1) here, n represents the prognostic genes, exprgene represents the coefficient gene.The coefficient of the gene in the risk signature represents the expression value of the gene.

Principal components analysis
The samples were divided into low and high-risk groups based on the median risk scores.The Principal Component Analysis (PCA) was employed to clarify the variances among these two groups.The variation in mRNA expression data among the high and low-risk groups in TCGA-, ICGC-and GSE76427-HCC datasets were observed using dimensionality reduction.

Prognostic analysis of vascular endothelial growth factor production-related risk signature
Both Kaplan-Meier survival curves and Cox regression analysis were applied to evaluate the prognostic significance of vascular endothelial growth factor production-related risk signature (VPRS) in HCC in training and verification datasets.The area under the ROC curves (AUC) was calculated for the predictive value of VPRS for the 1-, 3-, and 5-year Overall Survival (OS) in HCC patients.The diagnostic procedure's precision was evaluated with the AUC-ROC (Area Under the Receiver Operating Characteristic curve), employing three defined categories based on [32]: low accuracy (0.5 < AUC-ROC ≤ 0.7), moderate accuracy (0.7 < AUC-ROC ≤ 0.9), and high accuracy (0.9 < AUC-ROC ≤ 1).

Clinicopathological importance of the vascular endothelial growth factor production-related risk signature
Patients in both training and verification sets were classified as low and high-risk sets.The chi-square test was performed to investigate the correlation among the clinicopathological features and risk score, among them, age, gender, N classification, T classification, M classification, Eastern Cancer Oncology Group (ECOG), vascular invasion, tumor grade, and stage.

Tumor-infiltrating immune cells profiles
The population of immune cells in low-risk and highrisk groups was estimated using the CIBERSORT computational method.In both the TCGA and ICGC-HCC cohorts, we utilized Pearson correlation analysis and the Wilcoxon test to examine the association between the fraction of Tumor-Infiltrating Immune Cells (TIICs) and the risk score.We then utilized the ESTIMATE algorithm, specifically the "estimate" package, to calculate the immune score, stromal score and tumor purity for each HCC sample [33].

Single-sample gene sets enrichment analysis
All significant genes of the 29 immune-related pathways were obtained from a previous study [34].The TIICs level was estimated by the single-sample GSEA (ssGSEA; [35]) using melanoma mRNA Tumor Mutation Burden (TMB) data.Moreover, the enrichment level of gene hallmarks for 29 hallmarks related to immune was analyzed in the training and verification cohorts, among the low and high-risk groups, with a significance threshold of p-value < 0.05.Additionally, given the significance of Immune Checkpoints (ICPs) and Immunogenic Cell Death (ICD) modulators in cancer immunity, the expression levels of these proteins were compared among the high and lowrisk groups.

Analysis of tumor mutational load
In the TCGA-HCC dataset, total tumor somatic mutations detected in each sample was used for computation.The prognostic value of TMB in HCC was also analyzed.The R package "maftools" was applied to calculate and compare the mutational status in high and low-risk groups.

Analysis of hallmark gene sets
Further analysis focused on hallmark gene sets derived from the Molecular Signatures Database.These sets represent and summarize explicit, well-defined biological states or displays and processes that are coherently expressed.The R package "GSVA" was applied on low and high-risk groups in TCGA-HCC to implement GSVA of hallmark gene sets [36].

Verification of hub genes of vascular endothelial growth factor production-related risk signature
In both training and validation datasets, the genes correlated with the prognosis of VPRS were evaluated with the Kaplan-Meier survival curves.For further investigation, the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) was adopted to assess the gene protein level identified in both normal liver and HCC tissues.

Quantitative real-time PCR (qRT-PCR)
RNA was extracted from liver tissues using the TRIzol reagent (as described in [37]) for qRT-PCR analysis.SYBR Green assays and specific primers were used to quantify gene expression.Relative Ct values, normalized to GAPDH, were employed to compare gene expression between the experimental and control groups.

Genetic alterations of vascular endothelial growth factor production-related genes in hepatocellular carcinoma
We performed a differential analysis of 34 VPRGs in the training group and identified 22 DEGs between normal and HCC cases (Figure 1A, 1B).Supplementary Table 1 shows the upregulated and downregulated VPRGs with their analogous logFC values.In addition, Spearman's correlation analysis was performed to characterize the expression correlations among the VPRGs and the results were presented in Figure 1C.A PPI network analysis also revealed a substantial co-expression correlation among the VPRGs, as demonstrated in Supplementary Figure 1A.Additionally, for further understanding of the genomic characteristics of VPRGs in HCC, we analyzed somatic mutational status and the copy number variation of GSRGs using the cBioPortal database (Supplementary Figure 1B).

Construction and verification of vascular endothelial growth factor production-related risk signature
From 34 VPRGs, univariate Cox analysis identified 7 genes significantly correlated with prognosis (p < 0.05) (Figure 2A).LASSO regression further narrowed this down to 4 key genes: ADAMTS3, CCR2, NDRG2, and NODAL (Supplementary Figure 2A, 2B).A risk score was then calculated for each patient based on the expression levels and coefficients of these 4 genes (Figure 2B).The PCA results further revealed significant separation between high-and low-risk groups based on the median VPRS in both training and validation sets (Figure 2C, 2D and Supplementary Figure 2C).
Survival analysis demonstrated significantly different clinical outcomes between the high-and low-risk groups in both training and validation sets (Figure 2E, 2F and Supplementary Figure 2D).Notably, the VPRS across TCGA, ICGC, and GSE76427-HCC datasets exhibited distinct patterns in risk gene expression, risk scores, and survival states, highlighting the VPRS's potential as a valuable prognostic predictor compared to other clinical factors (Supplementary Figure 2E-2G).Supplementary Table 2 summarizes the 4 genes and their corresponding coefficients in the optimal model.Altogether, above results demonstrated that the risk score based on VPRS could be a better indicator for predicting the prognosis of HCC compared with other clinical factors.

The risk score acts as an independent factor for predicting overall survival of hepatocellular carcinoma patients
We conducted both univariate and multivariate COX analysis to assess the independent prognostic value of the VPRS to predict the outcomes of HCC patients.In the TCGA-HCC group, the risk score independently predicted the OS (overall survival) (HR = 2.137, p = 0.002) (Figure 3A, 3B).Additionally, the risk score exhibited a higher AUC-ROC unlike additional clinical features, including age, gender, ECOG, vascular invasion, tumor grade, and stage.The AUC of the risk score in the training group was 0.679, 0.753, and 0.684 for 1-, 3-, and 5-year OS of patients, respectively (Figure 3C-3E).The findings were also validated in the ICGC-and GSE76427-HCC cohort.In multivariate COX regression, the risk score had an HR of 1.320 in ICGC (Figure 3F, 3G, p = 0.031), while it was 1.226 in GSE76427 (Supplementary Figure 3A, 3B, p = 0.039).Moreover, the AUC values of 1-, 3-, and 5-year OS of patients in ICGC-HCC were 0.692, 0.591, and 0.539, AGING respectively (Figure 3H-3J), while there were 0.474, 0.419, and 0.473 in GSE76427-HCC, respectively (Supplementary Figure 3C-3E).

Relationship between vascular endothelial growth factor production-related risk signature and the features of clinicopathology
During the training phase, we identified a total of 310 cases that included clinical data such as age, gender, T classification, N classification, M classification, ECOG, vascular invasion, tumor grade, and stage (Figure 4A).In the ICGC sets, we narrowed down the cases to 231 by considering only those with clinical features of age, gender, and tumor stage (Figure 4B).To compare the distribution of clinical characteristics among different risk groups, we utilized the "chisq.test"function in R to conduct Chi-square tests.The results of these tests for both the TCGA and ICGC cohorts can be found in Table 1.Based on VPRS, the risk scores substantially correlated with age, tumor grade, and T classification in TCGA cohorts (Figure 4C-4F).Furthermore, we noted a major correlation between the risk scores and tumor stage in ICGC groups (Figure 4G).HCC samples with higher age, T classification, tumor grade, and stage specifically exhibited markedly higher risk scores than those with lower values, with risk score and vascular invasion also showing a significant correlation.Thus, the risk score values were significantly correlated with age, T classification, vascular invasion, tumor grade, and stage of HCC.

Tumor-infiltrating immune cells profiles
TIICs are essential components of the TME.This study examined the relative fraction of 22 immune cell types before computing using the "CIBERSORT" algorithm.The boxplots display the differences in fractions among low and high-risk sets (Figure 5A, 5D).Notably, the highrisk group had substantially higher proportions of M2,  M1, and M0 macrophages, and regulatory T cells (Tregs) than the low-risk group.Additionally, M1 macrophages and risk scores had a significant negative correlation (Figure 5B, 5E), whereas M2 macrophages positively correlated to risk scores (Figure 5C, 5F) in TCGA and ICGC cohorts, respectively.Consequently, it is suggested that the risk scores associated with the decreased M1 macrophage and increased M2 macrophages may alter the prognosis of HCC patients.The ssGSEA scores were used to quantify the abundances and activities of various pathways, immunocytes, or functions.Higher ssGSEA scores revealed significant infiltration of immune cells and more active immune-related pathways.The heatmaps (Figure 6A, 6C) and boxplots (Figure 6B, 6D) revealed that the high-risk group had higher ssGSEA scores for various immune cell types, such as checkpoint, macrophages, and Tregs.Furthermore, we evaluated the stromal and immunological scores, which represent the presence of stromal cells and infiltration of immune cells in the tumor microenvironment, respectively.The analysis revealed that samples in the high-risk group exhibited significantly higher immune and stromal scores unlike those in the low-risk set (Figure 6B, 6D).The HCC samples with elevated risk scores exhibited higher levels of infiltrating stromal and immune cells.The TIICs and immune-related pathways were not substantially distinct among high and low-risk groups.In the high-risk group, the levels of immunosuppressive TIICs like M2 and Tregs were substantially higher than the baseline.

Correlation between immune modulators and risk score
The levels of genes related to ICD and ICPs were analyzed to explore the potential impact of these modulators on anticancer immunity.All 20 genes, related with ICPs, showed differential expression among the different risk groups in the TCGA and ICGC cohort (Figure 7A, 7C).The expressions of key ICPs including CTLA4, CD274 (PD-L1), and PDCD1 (PD-1) were significantly elevated in the high-risk group.Moreover, both in TCGA and ICGC cohorts, all 12 ICD genes had differential expression among the low and high-risk groups (Figure 7B, 7D).These findings imply that the risk score value, which reflects the expression levels of ICD and ICPs modulators, could be a promising immune treatment biomarker.

The correlation of risk score with somatic mutation rates and gene set variation analysis
The mutation landscape of both groups was analyzed.
Results showed that the top 20 frequently mutated genes in each subtype included PCLO, ALB, MUC16, TTN, TP53, and CTNNB1 (Figure 8A, 8B).These suggests that the VPRS-based risk score in HCC patients can predict the somatic mutation rate.Moreover, patients with higher risk scores showed good response to anticancer immunity.Using the GSVA, the pathway activities were compared between the low and high-risk groups.Signaling pathways associated with cell cycle, tumor proliferation, tumor invasion, angiogenesis, and tumorigenesis were significantly enriched in the high-risk groups (Figure 8C).The pathways included MYC target v1 and v2, E2F target, G2M checkpoint, MTORC1 signaling, epithelial-mesenchymal transition, and others.Collectively, these findings showed that the VPRSbased risk score can be used as a novel biomarker of HCC, associated with the signaling pathways involved in angiogenesis and tumor invasion.

Identification of hub genes of vascular endothelial growth factor production-related risk signature
The prognostic value of the 4-gene VPRS was evaluated in the training and validation groups (Figure 9A-9H).Analysis of the effect of the 4 genes on the survival outcomes (Table 2) revealed that CCR2, NDRG2, and NODAL levels were positively correlated to the OS of HCC patients.On the other hand, those with higher ADAMTS3 levels predicted poor OS.Given the significant correlation between these four genes and the prognosis of HCC patients, ADAMTS3, CCR2, NDRG2 and NODAL were considered as hub genes.Next, we explored the protein expression of related hub genes utilizing the HPA database.Consequently, the expression of ADAMTS3 (liver#img/ENSG00000156140-ADAMTS3/HPA021369 available from v22.0.proteinatlas.org;liver+cancer#img/ENSG00000156140-ADAMTS3/HPA021369 available from v22.0.proteinatlas.org)and NODAL (liver#img/ENSG00000156 574-NODAL/HPA045201 available from v22.0.proteinatlas.org;liver+cancer#img/ENSG00000156574-NODAL/HPA045201 available from v22.0.proteinatlas .org)were significantly higher in HCC tissues compared with the normal liver tissues, whereas NDRG2 (liver#img/ ENSG00000165795-NDRG2/HPA002896 available from v22.0.proteinatlas.org;liver+cancer#img/ENSG0000016 5795-NDRG2/HPA002896 available from v22.0.proteinatlas.org)expression was lower in HCC tissues.The CCR2 protein (liver#img/ENSG00000121807-CCR2/CAB003793 available from v22.0.proteinatlas.org;liver+cancer#img/ENSG00000121807-CCR2/CAB00379 3 available from v22.0.proteinatlas.org)expression was not detected neither in HCC nor normal tissue (Figure 9I-9L).Furthermore, the expression levels of the four hub genes in our hospital samples were measured using the real-time PCR technique (Figure 9M).Supplementary Table 3 shows the primer sequences of 4 hub genes.The mRNA expression of ADAMTS3, CCR2, NDRG2, and NODAL was significantly higher in HCC than in normal liver tissues, similar to the documented results in the public database.

Establishing a prognostic nomogram
A new prognosis nomogram was developed by combining risk score and tumor stage to predict the survival of HCC patients (Figure 10A).The nomogram predicted the OS of patients at 1, 3, and 5 years, with higher scores indicating lower survival probability.The results were comparable with the tumor stage, as higher stage status corresponded to a worse prognosis.To assess the effectiveness of the nomogram in predicting patient survival rates, we performed ROC curve analyses for 1-, 3-, and 5-year survival predictions in both training and validation sets.In the training set, the AUC values indicating the model's discriminatory power was 0.619, 0.693, and 0.698 for 1-, 3-, and 5-year survival outcomes, respectively (Figure 10B).In the validation set, AUC values were 0.844, 0.638, and 0.634 for the same prediction periods (Figure 10C).

DISCUSSION
HCC is characterized by high vascularity and degree of malignancy due to impaired angiogenesis [38,39].This is because the tumor vascular network determines the    supply of oxygen and nutrients to the growing tumor cells [40,41].Existing blood vessels are transformed into new ones through a process referred to as angiogenesis [42].The process of angiogenesis has been shown to modulate cancer progression [43].Antiangiogenic therapy is now recognized as an effective strategy for treating cancers.Numerous retrospective studies have shown that microangiogenesis can contribute to the early recurrence and prognosis of HCC [44,45].Microangiogenesis levels are often linked to a higher probability of recurrence even poorer survival outcomes [46].Given the negative prognostic implications of microangiogenesis in recurrent HCC, clinicians should consider whether more aggressive treatment options are necessary for such patients, particularly those experiencing recurrence within two years of their initial curative hepatectomy.This is because HCCs that recur within this time frame are assumed to arise from residual microscopic lesions.However, this issue has not been investigated.Limited research has explored this topic [45].VEGF functions as an irreplaceable role in the angiogenesis of HCC tissues which is bound up with its tumorigenesis.
Although some studies have constructed predictive models of angiogenesis in HCC, none of them has involved VEGF as an important factor in promoting angiogenesis in HCC.
For our study, we first obtained the VPRGs from the Molecular Signatures Database, and found the expression levels of most VPRGs substantially fluctuated among HCC and normal tissues.In addition, these genes had lower CNV and SNP in HCC tissues, implying that these genes are associated with HCC progression.Similarly, it suggested that they had high genetic stability and could inhibit the gene mutations.Through single-factor analysis, we identified prognostic genes in HCC.Using Lasso regression, we built a risk model based on four vascular pattern regulatory genes (VPRGs): ADAMTS3, CCR2, NDRG2, and NODAL.This model outperformed other clinical factors in predicting patient prognosis significantly.Furthermore, this signature showed better predictive performance for HCC prognosis compared to other clinical factors.A positive correlation was observed between the constructed risk score and the grade, stage, and T stage of the tumor, demonstrating a substantial positive correlation between the risk model we constructed and the HCC progression.Moreover, there was a significant association of our risk score with the clinical characteristic of vascular invasion.Specifically, the risk of significant vascular invasion increases with score, confirming the correlation between our constructed risk model and angiogenesis.We further investigated the validity of our risk model by examining its correlations with clinical features of HCC patients.Through GSVA analysis, the potential mechanisms driving the model's predictive power were also explored.Furthermore, experiments were conducted to validate the aberrant expression of four key model-related genes identified as VPRGs in HCC tissues.The results showed significant upregulation of ADAMTS3 and NODAL, while CCR2 and NDRG2 were downregulated.A Disintegrin and Metalloproteinase with Thrombospondin Motifs 3 (ADAMTS3) is a metalloproteinase, belonging to the ADAMTS family, which is capable of degrading components of the extracellular matrix, like the elastic fibers, collagen and proteoglycans [47].ADAMTS3 participates in normal physiological processes, like the embryonic development, tissue remodeling and wound healing [48].A study by Kim et al. reported that ADAMTS3 was highly expressed in glioma stem cells and significantly correlated with its tumorigenic and proliferative activities [49].However, at present, no study has reported the correlation between ADAMTS3 and HCC.We found, for the first time, that ADAMTS3 was highly expressed in HCC samples in our study.
Based on these findings, future mechanistic studies should investigate the role of ADAMTS3 in hepatocarcinogenesis. C-C motif chemokine receptor 2 (CCR2), a C-C chemokine receptor, is mainly expressed on immune cells, like the monocytes, macrophages, and T-cells, and can recognize and bind to chemokines, such as CCL2, as well as modulate the migration of immune cells to inflammatory and tumor sites [50].In the advanced stage of HCC, CCR2 inhibits anti-tumor immune responses by recruiting and activating immunosuppressive cells, such as the macrophages, myeloid-derived suppressor cells, and regulatory T cells, thereby promoting immune escape of HCC cells [51].In addition, CCR2 promotes HCC development by modulating fibrosis and inflammatory responses in the liver [52].In contrast, in the early stages of HCC, immune cells such as macrophages recruited by CCR2 can prevent early tumor formation by eliminating senescent hepatocytes [53].This intricate observation suggests that CCR2 plays a multifaceted role in HCC progression.It can promote and suppress cancer, which need to be further investigated to clarify the mechanisms governing its divergent effects [54,55].NDRG2, a tumor suppressor gene, has been shown to inhibit the growth and metastasis of many malignant tumors, including HCC, and its expression level is positively associated with the prognosis and survival outcomes of cancer patients [56][57][58].In HCC, the expression of NDRG2 is downregulated, which is consistent with our findings, and down-regulation of NDRG2 expression can significantly increase tumor angiogenesis via the VEGFA pathway [59].Nodal Growth Differentiation Factor (NODAL) is a secreted protein belonging to the transforming growth factor β (TGF-β) superfamily.It plays a role in physiological AGING processes such as embryonic development, tissue remodeling and stem cell maintenance [60].However, NODAL are highly expressed in several malignant tumors, such as melanoma, breast cancer, gastric cancer, HCC, etc., and negatively correlates with the degree of differentiation, clinical stage, metastasis, and prognosis of the tumors [61][62][63].NODAL is an important component of the TGF-β signaling pathway which enhances the proliferation, migration, and invasion of HCC cells by promoting the phosphorylation of Smad3 and the expression of Snail, drug resistance, and epithelial-mesenchymal transition (EMT) [64].In addition, NODAL has been shown to modulate tumor cell plasticity by promoting the formation of angiogenic mimetic structures [7].These results suggest that these signature-associated VPRGs contribute to HCC development and progression, and further studies are advocated to explore their mechanisms and identify potential therapeutic targets.
Immune cells determine response to immunotherapy and tumor survival as the primary component of TIME.
For example, clinical research confirms that the immune cell composition of tumors in HCC is strongly associated with patient prognosis and affects therapeutic response [65][66][67].Individuals in the low-risk group experienced intense immune cell infiltration compared to their counterparts in the high-risk group.Highdensity Tumor-Associated Macrophages (TAM) infiltration in HCC is a poor prognosis marker, and TAM is primarily M2 macrophages [68,69].M2 macrophages can up-regulate the PDL1 expression in HCC, limiting the activity of CD8+ T cells [70].In our study, high-risk patients exhibited higher levels of M2 macrophage and Treg cell infiltration, but M1 macrophage infiltration was minimal.Our results strongly demonstrated that the risk scores were negatively correlated with the degree of M1 macrophage infiltration and positively associated with the degree of M2 macrophage infiltration.In our analysis, the high-risk group exhibited higher immune and interstitial scores unlike the low-risk group, although the tumor purity was lower, suggesting better response to immunotherapy.Therefore, boosting antitumor immune responses is critical for effective clinical treatment.
Cancer immunotherapy is a new treatment approach for tumors which aims to restart the tumor immune cycle and reestablish the normal anti-tumor immune response [71,72].Several studies have investigated the effect of tumor immunotherapy on the TIME [73,74].The ICP is an effective approach for stimulating the anti-tumor immune response.The PD-1-and CTLA-4-activated immune T cells are the reliable point for the treatment.Substances that suppress PD-1 and CTLA-4 could significantly enhance the treatment for advanced cancers [75].Chemokines can also influence cancer migration and progression [76,77].Specifically, several key immune checkpoints (e.g., PD-1, CTLA-4) and genes involved in immunogenic cell death pathways (e.g., ERp5, CRTL) were upregulated in the high-risk group.Therefore, higher gene expression levels are more sensitive to ICIs and immunotherapy.Thus, our risk-scoring model is expected to improve individualized treatment of HCC patients.
The intricacy of HCC pathogenesis pathways leads to variation across diverse populations in various countries and regions, posing a significant challenge to customized therapy [78,79].Therefore, the use of a single DEG as a biomarker for customized HCC patients is not reliable.Moreover, it is difficult to identify and characterize the immune properties of HCC.In our study, we developed a new prognostic biomarker of HCC that can differentiate the immune status and the malignancy of HCC as well as predict the efficacy of immune checkpoint blockade (ICB) treatment in HCC patients.In addition, a microangiogenesis prediction model for HCC prognosis was developed and validated.Nevertheless, we acknowledge that our study has certain limitations.Firstly, the retrospective nature of our study, although validated with diverse datasets to assess performance, limits the application of our model.Secondly, the TIME exhibited significant heterogeneity, making it difficult to perform precise evaluation.Despite efforts to minimize biases by leveraging the relative order of immune cell infiltration, residual biases may persist.Future biological analyses are advocated to expand our understanding of survival outcomes.

Figure 2 .Figure 3 .
Figure 2. Construction of 4-genes VPRS.(A) Forest plot for the survival analysis of HCC patients with a univariate Cox model after adjustment for VPRGs; red color represents p < 0.05.(B) Radar diagram of efficiency of the 4 genes in VPRS; the closer the red dot is to the outside, the greater the value it represents.(C) PCA of HCC samples in TCGA; dots in red and green represent samples in high-risk and low-risk groups respectively.(D) Overall survival analysis of risk score for HCC patients in TCGA.(E) PCA in ICGC-HCC.(F) Survival analysis in ICGC-HCC.

Figure 4 .
Figure 4.The correlation between clinicopathological factors and risk score.Heatmap of the correlations between clinicopathological characteristics of HCC and risk score in the TCGA (A) and ICGC (B) cohorts.Distribution of vascular endothelial growth factor production-related risk signature within HCC patients stratified by age, tumor grade, T classification, and vascular invasion in TCGA (C-F) and tumor stage in ICGC (G) cohorts.

Figure 5 .
Figure 5.The association between TIICs and VPRS.(A) Differential analysis of 22 kinds of TIICs in the low and high-risk groups in training sets.(B, C) Spearman's correlation analysis for risk score and M1 and M2 macrophages in TCGA cohorts, each dot plot represents a subject, and the correlation is fitted into a straight blue line.(D) Differential analysis of 22 kinds of TIICs in the low and high-risk groups invalidation sets.(E, F) Spearman's correlation analysis for risk score and M1 and M2 macrophages in ICGC cohorts, each dot plot represents a subject, and the correlation is fitted into a straight blue line.R, rho; ***p < 0.001, **p < 0.001, *p < 0.05.

Figure 7 .
Figure 7. Correlation between risk subtypes and ICPs and ICD modulators.(A) Differential expression of ICP genes among the risk subtypes in TCGA cohorts.(B) Differential expression of ICD modulator genes among the risk subtypes in TCGA cohorts.(C) Differential expression of ICP genes among the risk subtypes in ICGC cohorts.(D) Differential expression of ICD modulator genes among the risk subtypes in ICGC cohorts.

Figure 8 .
Figure 8. Gene set variation analysis and correlation between mutation and risk subtype.(A) Top 20 highly mutated genes in HCC high-risk group.(B) Top 20 highly mutated genes in HCC low-risk group.(C) Heatmap for the contribution of GSVA scores of hallmarks in highand low-risk groups.The red color represents up-regulated terms in the high-risk group, green color shows the down-regulated terms in the low-risk group.

Figure 9 .
Figure 9. Verification of the prognostic value and expression of hub genes of VPRS.Survival analysis of ADAMTS3, CCR2, NDRG2, and NODAL for patients in TCGA (A-D) and ICGC (E-H) cohorts.The protein expression level of ADAMTS3 (I), CCR2 (J), NDRG2 (K), and NODAL (L) in normal and HCC tissues based on the HPA database.(M) The relative mRNA expression levels of ADAMTS3, CCR2, NDRG2, and NODAL are compared between HCC and non-tumor tissues with real-time PCR results.***p < 0.001, **p < 0.001, *p < 0.05.

Figure 10 .
Figure 10.Prognostic nomogram was established by combining risk score and tumor stage characters.(A) Nomogram for assessing the 1-, 3-, and 5-year OS for HCC patients in the TCGA dataset.ROC curves of 1-, 3-, and 5-year in the TCGA (B) and ICGC (C) datasets.