Integrated Analysis of Angiogenesis-Mediated Tumor Immune Microenvironment Pattern in Hepatocellular Carcinoma (HCC) and a Novel Prognostic Model Construction to Predict Patient Outcome

Background

As one of the most important public health concerns around the world, HCC has become the sixth most common neoplasm and the fourth primary cause of cancer-related death all over the world [1]. Because of its biology characteristics heterogeneity, complex pathogenesis, and aggressive biology features, HCC is infamous for its poor therapy responsiveness, low early detection rate, and poor patient outcomes [2]. In recent years, comprehensive treatment combined with immunotherapy, radiotherapy, chemotherapy and surgery is widely applied clinically. However, because of its low rate of early detection and high disease recurrence rate, few patients have the chance to receive comprehensive and effective treatment [3,4]. Thus, there is an urgent need for thorough research on the mechanism underlying HCC tumorigenesis and progression and related poor prognosis of certain tumor subtypes to facilitate precise treatment.

Because of the expanding tumor mass during HCC cell proliferation, its elevated high oxygen and metabolic substrate requirement have become important tumor characteristics [5–7]. Angiogenesis is the core biological feature of HCC formation. Angiogenesis begins with an imbalance between pro-angiogenic and anti-angiogenic factors [8]. Vascular endothelial growth factor (VEGFA) is the critical promoter of the tumor’s endothelial proliferation, migration, and neovascularization to support its rapid tumorigenesis and local metastasis by activating the PI3K/AKT signaling pathway [9,10]. Angiogenesis plays a crucial role in HCC progression, and drugs targeting angiogenesis, such as sorafenib and lenvatinib, are important components of first-line treatment and systemic treatment of HCC [11,12]. Predictive biomarkers such as pERK and imaging criteria have been proposed for identification of tumors with the active angiogenesis subtype, but none of them have been approved for clinical use [13,14], and this gap has been one of the obstacles to precise treatment of patients with the active angiogenesis HCC subtype. To solve this problem, we developed a machine-learning-based prognostic model to help evaluate a tumor’s angiogenesis activity and to predict patient outcome.

The present study analyzed the angiogenesis-mediated biological features and tumor microenvironment immune cell infiltration alterations to reveal its underlying pro-tumor mechanism. To promote angiogenesis HCC subtype identification and related individualized treatment in clinic practice, we constructed an angiogenesis-based prognostic model based on the clinical information and transcription data acquired from The Cancer Genome Atlas (TCGA) and assessed the model’s accuracy based on the dataset acquired from GSE14520 and ICGC-JP. According to our results, our prognostic model has satisfactory accuracy in predicting HCC patient outcomes and shows potential to be implemented for HCC diagnosis and precision treatment.

Material and Method

DATA ACQUISITION AND PROCESSING:

Transcription data and clinical information of the training set were obtained from The Cancer Genome Atlas LIHC project (TCGA-LIHC). The samples that had transcription data and follow-up information were included in analysis. The TCGA data were obtained using the R package “TCGAbiolinks” [15]. The RNA-seq data were normalized with transcripts per million (TPM) method and then transformed to log2 (TPM+1). The validation set’s clinical and transcription data were acquired from ICGC-JP (https://dcc.icgc.org) and GSE14520 [16], separately. To make the gene expression profiling suitable for comparison, the RNA-seq data of ICGC-JP were also normalized by log2 (TPM+1) format, and the microarray’s probe value of GSE14520 was also log2(X+1) transformed.

The angiogenesis-related genes markers were downloaded from the gene set “Angiogenesis” from MSigDB [17].

IDENTIFICATION OF THE SUBTYPE WITH DIFFERENT ANGIOGENESIS STATUS:

To identify the HCC subtype with different angiogenesis status, we conducted consensus clustering analysis using the R package “ConsensusClusterPlus” [18]. Consensus clustering analysis was based on the angiogenesis-related genes, with 80% of the total sample drawn 100 times.

TUMOR IMMUNE INFILTRATION ANALYSIS:

To estimate the proportion of the immune-infiltrating cells in the tumor microenvironment, we used the R package “quanTIseq” to analyze the tumor subtype’s immune status [19].

DEG’S ACQUISITION AND ENRICHMENT ANALYSIS:

To acquire the DEGs between the active angiogenesis subtype and inactive angiogenesis subtype for further function analysis, we used the R package “Limma” to obtain DEGs. DEGs were identified under the threshold of log2 (Foldchange) >2 and P<0.05. To analyze the DEGs-mediated biology function, we conducted enrichment analysis using the R package “ClusterProfiler” [20]. Function enrichment analysis was based on the Gene Ontology (GO) [21] dataset, which includes molecular function (MF), biological pathways (BP), and cellular components (CC) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [22] dataset.

CONSTRUCTION AND VALIDATION OF THE PROGNOSTIC MODEL:

To select the genes among DEGs for prognostic model construction, we used univariate Cox regression to analyze the genes’ impact on patient OS. Genes with P value <0.01 were selected and submitted for the subsequent LASSO Cox regression. Then, the R package “glmnet” was used to conduct LASSO Cox regression to determine the features and assess their contribution as an independent prognostic factor for patient OS [23]. In this way, the novel 13-genes prognostic model based on the angiogenesis features was developed according to the feature’s regression coefficients (β), which were obtained from the LASSO Cox regression model multiplied by its related mRNA expression. Then, we computed the risk score of each sample in the training set according to our prognostic model and divided them into a high-risk group and a low-risk group based on the median cut-off. Patient outcome was compared by log-rank test and the receiver operating characteristic (ROC) curve was plotted using the R package “survivalROC” [24]. The risk score of samples of the validation set were also assessed and the samples were also divided into a high-risk group and low-risk group and further compared in the same way.

CELL CULTURE AND QUANTITATIVE REAL-TIME PCR:

Human HCC cell line (HepG2) and normal liver cell line (LO2) were cultured in RPMI-1640 medium (Hyclone, Beijing, China) with 10% FBS (Hyclone, Beijing, China) and 1% streptomycin-penicillin (Hyclone, Beijing, China). The cells were maintained in a humidified chamber containing 5% CO₂ at 37°C.

Total RNA of the HCC cells was extracted using TRIzol (Invitrogen, Beijing, China). First-strand cDNA was synthesized using EasyScript First-Strand cDNA Synthesis SuperMix (TransGen, Beijing, China). Quantitative real-time PCR was conducted using SYBR Green qPCR SuperMix (Vazyme, Nanjing, China) in the ABI PRISM^® 7500 Sequence Detection System (Applied Biosystems, Shanghai, China). Relative mRNA levels were calculated using the −ΔΔCt method using 18S rRNA as control and presented as 2^(−ΔΔCt). The primers were:

STATISTICAL ANALYSIS:

During the comparison of 2 groups, we used the Wilcox test. In the comparison among multi-groups, we used the Kruskal-Wallis test. Analysis of correlations between 2 continuous variables was conducted by Spearman correlation. For the quantitative real-time PCR assay, data are presented as mean±standard deviation (SD) and the comparison between 2 groups was conducted by t test. If not specifically mentioned, P<0.05 was considered statistically significant.

Results

IDENTIFICATION OF ACTIVE ANGIOGENESIS SUBTYPE:

First, we used consensus clustering analysis to identify the active angiogenesis subtype based on the angiogenesis markers derived from MSigDB. A graphical representation of consensus clustering results is shown in a CDF curve and delta area plot (Figure 1A, 1B). According to the results, the subgroups were divided under the K value=4. To validate the robustness of the subtype identification, we test the angiogenesis marker-VEGFA’s expression in the four subtypes (Figure 1C). G4 performs the significant higher expression level of VEGFA while G1 has the significant lower expression level of VEGFA. Thus, G4 was identified as angiogenesis active subtype and G1 was identified as angiogenesis inactive subtype. Meanwhile, Kaplan-Meier plot was plotted and the result demonstrated that angiogenesis active group patients have shorter survival time comparing with the patients in angiogenesis inactive subtype. (P=0.0561; Figure 1D), which indicate that active angiogenesis can be the biology characteristics of aggressive HCC

IMMUNE INFILTRATION ANALYSIS OF ACTIVE ANGIOGENESIS SUBTYPE AND ANGIOGENESIS INACTIVE SUBTYPE:

To further analyze the immune microenvironment pattern differences between the active angiogenesis subtype and inactive angiogenesis subtype, we applied quanTIseq to estimate the infiltration level of immune cells. Results indicated that the active angiogenesis subtype has a significantly higher M2 microphage infiltration level than the inactive angiogenesis subtype (Figure 2A). To assess the M2 microphage’s activity and immune-suppressive function, M2 microphage-related genes IL4, IL13, CSF1, IL10, GFB1, and TGFB2 were enrolled in the analysis (Figure 2B) [25]. A higher M2 microphage-related genes expression level suggests that the active angiogenesis subtype’s infiltrating M2 microphage has active immune-suppressive activity.

IDENTIFICATION OF DEGS AND RELATED FUNCTION ANALYSIS:

To analyze the angiogenesis-mediated biology heterogeneity between the 2 subtypes, we acquired DEGs between the 2 subtypes. Among the 1179 DEGs, 953 are upregulated genes and 226 are downregulated genes. The result was visualized by volcano plot and heat map (Figure 3A, 3B).

To further research the biology function that these genes mediate, enrichment analysis was performed. According to Figure 3C, the PI3K-Akt signaling pathway and human papillomavirus infection are the significantly enriched terms in the KEGG pathway (up). Among the KEGG pathway (down), complement and coagulation cascades is the most enriched term. In the GO pathway (up), extracellular matrix organization is the most enriched term. In the GO pathway (down), small-molecule catabolic process is remarkably enriched (Figure 3C–3F).

CONSTRUCTION OF THE ANGIOGENESIS-BASED GENE PROGNOSTIC MODEL:

To support clinical application, we constructed a gene prognostic model to predict HCC patient outcome. Univariate Cox analysis was used to determine the genes that are suitable for LASSO regression among the DEGs. Genes with P<0.001 were considered to have significant impact on patient outcome and were submitted to subsequent LASSO Cox regression analysis. According to the LASSO Cox regression analysis results, 13 genes would be optimized to construct the gene prognostic model (Figure 4A, 4B). The risk score formula reads as follows:

To validate the efficiency of the prognostic model, we divided patients into high- and low-risk groups according to the median cut-off (Figure 4C). The Kaplan-Meier survival curve demonstrated that compared with the low-risk group, the high-risk group had a significantly worse outcome (P=2.03e-07, Figure 4C–4F). ROC curves indicated that the model has the best prediction ability for 1-year patient outcome (Figure 4G).

GENE PROGNOSTIC MODEL VALIDATION:

The external sets GSE14520 and ICGC-JP were employed to examine the prediction accuracy of the gene prognostic model. Kaplan-Meier analysis indicated that the high-risk group had significantly worse prognosis than the low-risk group (Figure 5A, 5C). ROC curves showed that AUC rate for 1-year, 3-year, and 5-year OS were 0.561, 0.626, and 0.671, respectively in GSE14520 and the AUC rates for 1-year and 3-year OS were 0.763 and 0.747, respectively, in ICGC-JP (Figure 5B, 5D). Since there were only 2 sample with OS more than 5 years, the 5-year survival is not included in ROC curve in ICGC-JP. These results proved the accuracy of our prognostic model. To further validate our prognostic model’s potential in predicting HCC’s microenvironment pattern, we analyzed the correlation between the risk score (including prognostic model component genes) and immune infiltration level and immune-related gene expression (Figure 6). Risk score and prognostic model component genes with coefficient>0 were positively correlated with M2 macrophage infiltration level and M2 macrophage-related genes expression level, while prognostic model component genes with coefficient <0 were negatively correlated. Our results show that the prognostic model accurately reflects HCC M2 macrophage infiltration and activity.

CDCA8 HAD HIGHER EXPRESSION LEVEL IN HCC CELL LINE:

Among all prognostic model component genes, CDCA8 had the highest absolute value of coefficient (Figure 7A). Thus, we further examined the expression levels of CDCA8 in normal liver cell line (LO2) and HCC cell line (HepG2), and the quantitative real-time PCR assay showed that the expression level of CDCA8 in the HCC cell line was significantly higher than in the normal liver cell line (Figure 7B).

Disscussion

As a malignant tumor with rising incidence and mortality rate in recent years, tumor heterogeneity has become one of the important obstacles for precise treatment of HCC [26]. Angiogenesis is the essential process during formation of the HCC aggressive subtype and is an important feature of the poor prognosis subtype. Thus, analyzing the tumor biology features and immune microenvironment alterations that angiogenesis confers would be essential for targeting therapy according to HCC angiogenesis subtypes. To facilitate clinical application, a robust method is needed to identify the angiogenesis HCC subtype for individualized therapy for HCC patients.

Here, we first identified the angiogenesis subtype among the TCGA-LIHC cohort. Survival analysis proved the significant impact of angiogenesis during formation of the poor prognosis HCC subtype. Angiogenesis status is an important predictor of poor outcomes for HCC patients.

The immune microenvironment plays an important role during tumor proliferation and progression [27]. Recent research demonstrated that angiogenesis can significantly influence the tumor immune microenvironment pattern and exert important effects on the anti-tumor response of immune cells [28,29]. Therefore, we assessed the immune cell infiltration level of the HCC tumor microenvironment. Our results show that in the active angiogenesis tumor subtype, M2 macrophages are associated with elevated infiltration level. M2 macrophages are immune-suppressive cell that exert a pro-tumor function by mediating CD8+T cell dysfunction in the anti-tumor response of cytotoxic killing [25]. Its pro-tumor function and division rely on the cytokines IL4, IL13, CSF1, IL10, TGFB1, and TGFB2. Thus, we analyzed the expression of IL4, IL13, CSF1, IL10, TGFB1, and TGFB2. The high M2 macrophage-related gene expression in the active angiogenesis subgroup indicates the high activity of M2 macrophage in this HCC subtype. In current clinic application, therapy combining angiogenesis inhibitor and immunotherapy has a satisfactory responsive rate [30]. Our research suggests the potential of this combined treatment strategy for use in individualized treatment for HCC patients with the active angiogenesis subtype.

To further analyze the angiogenesis-mediated biology features alterations, we used enrichment analysis to determine the DEGs between the 2 subgroups of mediated biology function. The upregulated genes are involved in the PI3K pathway, extracellular structure organization, and extracellular matrix organization. The PI3K pathway is the essential cell signal pathway for cancer cell proliferation [31]. Ying Zhou also demonstrated that PI3K/Akt pathway activation is the essential part of angiogenesis in hepatocellular carcinoma [32]. Extracellular structure organization and extracellular matrix organization are important processes during tumor expansion [33,34]. Extracellular matrix organization and scaffolding can influence the rate of angiogenesis and capillary integrity and are required for blood vessel morphogenesis [35]. The results further proved the accuracy of our subtype identification and indicate that angiogenesis may be an important factor causing the frequent progression of HCC. Downregulated genes are involved in the inherent metabolic function of a normal liver. This result indicates that the active angiogenesis tumor subtype is associated with more significant low differentiation characteristics compared to the inactive angiogenesis tumor subtype.

We also constructed a prognostic model based on angiogenesis-related genes to assist clinicians to estimate a tumor’s angiogenesis potential and predict patient outcome. In the external validation sets GSE14520 and ICGC-JP, the prognostic model was accurate. In analysis of the prognostic model’s efficiency in predicting HCC’s microenvironment pattern, the prognostic model also performed well. These results prove our prognostic model’s robustness. In the quantitative real-time PCR assay, the HCC cell line had the highest expression level of CDCA8. The process of this research is presented as a workflow diagram in Figure 8. Our prognostic model construction genes’ regulating role in HCC was also reported in previous research. As a worse prognostic marker, TG6PD is essential for tumor proliferation and angiogenesis and is correlated with HCC invasion [36,37]. It has also been reported that high MMP1 expression is important in progression of HCC angiogenesis and is a tumor invasion driver [38,39]. Psallidas et al revealed that SPP1 promotes HCC metastasis and tumor angiogenesis [40,41]. A previous study indicated that HILPDA is a marker of poor prognosis. The expression level of HILPDA was significantly correlated with TAM infiltration level in most neoplasms and was associated with tumor immunosuppression status [42].

Conclusions

To sum up, we identified the angiogenesis subtype among the liver cancer and analyzed its impact on tumor biology features, immune microenvironment pattern, and patient outcome. Finally, we developed a reliable novel gene prognostic model to predict patient outcome and tumor angiogenesis potential. Therapy targeting specific HCC characteristic has significantly improved patient survival time [43]. An efficient method which can help identify tumor subtype with specific characteristic would be meaningful for individualized precision treatment. This prognostic model can potentiality assist clinicians in decision-making.