Development and Validation of a Propionate Metabolism-Related Gene Signature for Prognostic Prediction of Hepatocellular Carcinoma

Introduction

Hepatocellular carcinoma (HCC) accounts for 75–85% of liver cancer cases.¹ While immune and targeted drug therapy are the first-line treatments for HCC, they are not suitable for all patients owing to the different biological characteristics of tumors.² Even though the above treatments have greatly improved survival rates, advanced HCC still has a poor prognosis. In current practice, it is widely used to predict the prognosis of HCC using the Barcelona Clinic Liver Cancer (BCLC) classification.³ However, some patients with similar tumor stages exhibit different outcomes, suggesting that the staging system is still lacking. Therefore, exploring new biomarkers for HCC patients and constructing a risk model to aid in the diagnosis and prognosis of HCC are of significant importance.

Propionate metabolism differs across species.⁴ Propionate growth triggers fatty acid breakdown, boosting three β-oxidation enzymes. Increased amino acid and protein degradation enzymes indicate their role in supporting propionate metabolism.⁵ Recent research has revealed that propionate metabolism dysfunction is a significant factor in cancer and also represents a valuable potential target in the treatment of metastatic cancer.⁶ There were several health benefits of propionic fermentation in colon, which include the prevention of fat deposition, lowering cholesterol, reducing inflammation, and preventing cancer development.⁷ Propionate metabolism imbalance heightens metastasis in breast and lung cancers via ERK2-activated SP1/EGR1 switch, downregulating methylmalonyl-CoA epimerase (MCEE). MCEE loss weakens propionate-driven anaplerosis, causing intracellular and intratumoral methylmalonic acid (MMA) buildup. MMA, a propionate metabolism byproduct, enhances cancer cell invasiveness, exacerbating metastasis.⁶ Besides inhibiting colon cancer cell proliferation and inducing apoptosis, propionate could act as an inhibitor of proliferation of HCC by activating the GPR43 signaling pathway.^8–10 Despite this, the role of the propionate metabolism-related genes (PMRGs) in the pathogenesis of HCC remains uncertain.

We constructed a novel prognostic model for HCC based on five prognostic biomarkers (ACADS, CYP19A1, FTCD, G6PD, and GOT2) that were identified in this study, in which the coefficient of ACADS, FTCD, and GOT2 were negative and that of the other were positive. The biomarkers were acquired by univariate Cox analysis and at least absolute contraction and choice of operator (lasso) analysis. Next, univariate and multivariate Cox analyses were performed to select the independent predictor to construct the nomogram. Additionally, we focused on the function and the tumor microenvironment (TME) of different risk groups. Our findings may prove useful in identifying potential targets for the clinical diagnosis of HCC and in understanding the mechanisms whereby PMRGs play a role in the development of HCC.

Materials and Methods

Data Extraction

The RNA sequencing data, survival, and clinical information of HCC patients were downloaded from TCGA database (https://portal.gdc.cancer.gov) and GEO database (http://www.ncbi.nlm.nih.gov/geo), respectively. The former contains 371 HCC and 50 normal control (HC) samples, and 365 HCC samples have both survival and clinical information. The GES14520 dataset (GPL3921, https://www.ncbi.nlm.nih.gov/search/all/?term=GES14520) contains 225 HCC and 220 HC samples, and 221 HCC samples have survival information. Furthermore, a total of 291 PMRGs (relevance score > 10) were retrieved from the GeneCards database (https://www.genecards.org; Supplementary Table 1) with “propionate metabolism” as the key work.

Functional Analysis of Differentially Expressed PMRGs (DE-PMRGs) of HCC

First, the differentially expressed genes (DEGs) between 371 HCC and 50 HC samples were compared using the “limma” R package (version 3.48.3) (|log₂FC| > 1, adj.p value < 0.05).¹¹ Then, DE-PMRGs were obtained by intersecting the DEGs and 291 PMRGs using “venn”. Furthermore, Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DE-PMRGs were conducted using the “clusterProfiler” R package (version 4.0.2) (count ≥ 1, adj.p value < 0.05).¹²

Construction of the Survival Risk Model of HCC

In this study, 365 HCC samples from TCGA database were collected for constructing the survival risk model. The survival risk model was constructed using biomarkers obtained through univariate cox analysis and least absolute shrinkage and selection operator (LASSO) analysis. The univariate cox analysis was implemented by “survival” package¹³ (version 3.2–13) and the results was displayed by “forestplot” package¹⁴ (version 2.0.1). Furthermore, “glmnet” package¹⁵ (version 4.1–3) was used to implemented LASSO regression analysis. The predict.coxph function in survival package was used in calculating the risk score, as follows: Riskscore = β₁X₁ + β₂X₂ +. + β_nX_n. β refers to the regression coefficient of the prognostic gene, and X represents the expression of the gene. The HCC patients were divided into the high- and low-risk groups on the basis of the median risk value.

To predict the accuracy of the survival risk model, the Kaplan–Meier (KM) curve and receiver operating characteristic (ROC) curve were analyzed between these two groups. In addition, 221 HCC samples from GSE14520 were collected and used as the validation dataset to verify this survival risk model. The risk curve of high- and low-risk groups and the expression heatmap of biomarkers were plotted by “ggplot2” package (version 3.3.5).¹⁶ The KM curve and ROC curve were painted by “survminer” package (version 0.4.9)¹⁷ and “survivalROC” package (version 1.0.3),¹⁸ respectively.

Independent Prognostic Analysis and Clinical Correlation Analysis

The Wilcoxon test was used to compare the correlation between the risk score and different clinical characteristics (age, stage, grade, cancer status, gender, pathologic stage M, N, and T). Then, the significant prognostic factors explored by univariate and multivariate Cox analyses were used to construct the nomogram. Afterward, the calibration curves, decision curve analysis (DCA), and ROC curves of the above model and prognostic factors were constructed to verify the validity of the nomogram.

Functional Analysis

On the one hand, the gene set enrichment analysis (GSEA) was performed to study the function of genes in different risk groups using the “clusterProfiler” R package (version 4.0.2) and the “org.Hs.eg.db” R package (version 3.13.0) (|NES| > 1, NOM P < 0.05, q < 0.25). Next, 86 invasion-related genes were obtained from cancerSEA database,¹⁹ 1184 epithelial–mesenchymal transition (EMT)-related genes were obtained from dbEMT2 database,²⁰ and 47 angiogenesis-related genes obtained from MSigDB database.²¹ The invasion score, angiogenesis score, and EMT score of each sample were calculated using the “ssGSEA” algorithm in “GSVA” package²² based on the expression of above genes. These three scores and the mRNAsi score were compared in these two risk groups. The mRNAsi stemness index of patients in the TCGA-LIHC cohort was collected from the Annex of the published literature.²³

TME Analysis

TME is a critical part of the tumor tissue.²⁴ During this study, the relative percent of 22 immune cells in the samples was calculated by CIBERSORT²⁵ algorithm (http://CIBERSORT.stanford.edu/). The score of 22 immune cells and 47 immune checkpoints²⁶ were compared between different risk groups by Wilcoxon test. Moreover, the immunophenscore (IPS), tumor immune dysfunction and exclusion (TIDE), and ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) analyses in different risk groups were calculated for assessing the tumor purity and the response to the immunotherapy. The IPS score of HCC patients was collected from The Cancer Immunome Atlas (TCIA) database (https://tcia.at/home). The score of TIDE prediction, dysfunction prediction, and exclusion prediction was retrieved from the TIDE online database (http://tide.dfci.harvard.edu/). In addition, the immune score, ESTIMATE score, stromal score, and tumor purity of samples in the high- and low-risk groups were calculated via ESTIMATE algorithm based on the expression of biomarkers.

Expression Verification of Biomarkers

The expression of each biomarker was verified in both TCGA database and GSE14520. After that, the expression of biomarkers was validated by using a quantitative real-time polymerase chain reaction (qRT-PCR). Ten pairs of HCC tumor and adjacent tissue samples were committed by the patients and the Ethics Committee of Huaihe Hospital (Supplementary Table 2). Total RNA was extracted using TRIzol (Thermo Fisher, Shanghai, CN), and mRNA was reverse-transcribed into cDNA using SureScript-First-strand-cDNA-synthesis-kit (Servicebio, WuHan, CN). The qRT-PCR reaction system was made up of 3ul of cDNA, 5ul of 2xUniversal Blue SYBR Green qPCR Master Mix and 1ul of each upstream and downstream primers. Finally, the reactions were performed on a CFX96 real-time quantitative fluorescence PCR instrument. The amplification reactions were programmed with pre-denaturation at 95 °C for 1 min, followed by 40 cycles, each cycle consisting of 95 °C for 20s, 55 °C for 20s, and 72 °C for 30s. The relative expression of genes was calculated by the 2^−ΔΔCt method using Gapdh as the internal reference gene. The primer sequences were listed in Table 1.

Table 1 The Primer Sequences for qRT-PCR Validation of Biomarkers

Results

A Total of 132 DE-PMRGs Were Associated with Steroid and Organic Acid Metabolic Processes

There were 3690 DEGs (1627 up-regulated and 2063 down-regulated) between 371 HCC and 50 HC samples (Figure 1A). Then, 132 DE-PMRGs were obtained by intersecting 3690 DEGs and 291 PMRGs (Figure 1B, Supplementary Table 3).

Figure 1 The differentially expressed propionate metabolism-related genes (DE-PMRGs) and functional enrichment analysis. (A) 3690 differentially expressed genes (DEGs) between liver Hepatocellular carcinoma (LIHC) and normal control (NC) samples. (B) The venn diagram of 132 DE-PMRGs. (C and D) The Gene Ontology (GO) functions (C) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (D) enriched by 132 DE-PMRGs. BP, biological progress; CC, cellular component; MF, molecular function.

These 132 DE-PMRGs were enriched to 867 Gene Ontology (GO) functions with respect to function, encompassing steroid metabolism, carboxylic acid, organic acid, cellular amino acid metabolism, etc. Additionally, these DE-PMRGs were linked to 49 KEGG pathways, including those for steroid hormone production, retinol metabolism, the breakdown of valine, leucine, and isoleucine, carbon metabolism, and others (Figure 1C and D).

Five Biomarkers Were Used for Constructing the Survival Risk and Prognostic Models of HCC

In this study, five biomarkers, namely, ACADS, CYP19A1, FTCD, G6PD, and GOT2 were identified. Among them, ACADS, FTCD, and GOT2 were positive factors (hazard ratio < 1), while CYP19A1 and G6PD were negative factors (hazard ratio > 1) (Figure 2A and B). The coefficient of biomarkers was shown in Table 2. Next, the risk score was calculated, and the KM curve and risk curve indicated that the survival between the two risk groups differed significantly (p < 0.0001) (Figure 2C and D). It was determined that the survival risk model could be used as a prognostic model if the area under the ROC curve (AUC value) of 1–5 years was greater than 0.6 (Figure 2E). Moreover, the results of KM curve, risk curve, and ROC curve in the validation dataset (GSE14520) were consistent with those of the training dataset (Figure 2F–H).

Table 2 The Coefficient of the Biomarkers

Figure 2 Construction and validation of prognostic model for LIHC. (A) The forest diagram of 44 biomarkers obtained by univariate cox analysis (p < 0.01). (B) The coefficients of biomarkers and the error plot for cross-validation in the least absolute shrinkage and selection operator (LASSO) analysis. The different colored lines in the graphic above represent different genes. (C) Kaplan-Meier plot of patients in a low- or high-risk group (p < 0.0001), and the number of patients in different risk groups. (D) Up: The risk score of each patient with LIHC. Medium: The patient survival based on the risk score. Down: The heat map of the 5 biomarkers in the high- and low-risk group. (E) Receiver operating characteristic (ROC) curves for the predictive value of the prognostic model for different years. AUC: area under the curve. (F–H) Validation of the prognostic model in GSE14520 dataset. (F) Kaplan-Meier plot of patients in a low- or high-risk group (p = 0.00065), and the number of patients in different risk groups. (G) Up: The risk score of each patient with LIHC. Medium: The patient survival based on the risk score. Down: The heat map of the 5 biomarkers in the high- and low-risk group. (H) ROC curves for the predictive value of the prognostic model for different years.

In addition, the risk scores were considerably different between stages I and II, stages I and III, grades 1 and 2, grades 1 and 3, grades 1 and 4, grades 2 and 3, tumor free and with tumor, pathologic M0 and MX, pathologic T1 and T2, pathologic T1 and T3, and pathologic T1 and T4 (p < 0.05) (Figure 3A). Three factors, namely, pathologic T, cancer status, and risk score were negatively associated with patient survival (Figure 3B and C). The nomogram with these three prognostic factors was constructed, and the slopes of survival probability in calibration curve were close to 1, which indicated that the prediction model could be used as an effective model (Figure 3D and E). In addition, the benefit rate and the AUC value of the nomogram model were higher than others, and the AUC value of these three prognostic factors were higher than 0.6 (Figure 3F and G). All these results suggested that the nomogram model had accurate predictive ability for HCC, and the prediction accuracy of the nomogram model was higher than that of the single factor.

Figure 3 Independent prognostic analysis and clinical correlation analysis. (A) Comparison of risk scores among different sub-types for different clinical factors. ns, not significant; *p<0.05; ** p<0.01; ***p<0.001; ****p<0.0001. (B and C) Three independent prognostic factors obtained by univariate (B) and multivariate (C) COX regression analyses. (D) The nomogram was constructed based on three independent prognostic factors to predict the probability of 1/3/5-year survival. (E and F) The calibration curve (E) and decision curve analysis (DCA) curve (F) of the nomogram. (G) ROC curves of the nomogram and independent prognostic factors.

The Genes in High-Risk Groups Were Associated with the Functions of Invasion and EMT

The GSEA results demonstrated that the functions of alpha amino acid, fatty acid catabolic and metabolic process, monocarboxylic acid, organic acid catabolic process, and others and the pathways, including complement and coagulation cascades, fatty acid metabolism, glycine serine and threonine metabolism and degradation, PPAR signaling pathway, and others were significantly highly enriched in the high-risk groups (p < 0.05) (Figure 4A and B, Supplementary Table 4). In addition, the invasion score, angiogenesis score, EMT score, and mRNAsi score were significantly different between these two risk groups. In detail, the invasion score, EMT score, and mRNAsi score were higher in high-risk groups (p < 0.05) (Figure 4C). According to these results, genes associated with invasion and EMT functions were significantly more prevalent in high-risk groups.

Figure 4 Gene set enrichment analysis (GSEA) of high-risk groups. (A and B) The GO terms (A) and KEGG pathways significantly enriched in the high-risk group. (C) Discrepancies of invasion score, angiogenesis score, mesenchymal transition (EMT) score, and mRNAsi score between high- and low-risk groups. *p < 0.05, **p < 0.01, ***p < 0.001; ****p<0.0001.

The Response to Immunotherapy Was Worse, and the Level of Immune Escape Were Higher in High-Risk Groups

There were six significantly different immune cells and 36 significantly different immune checkpoints in different risk groups. There was a significant increase in the expressions of activated memory CD4 T cells, regulatory T cells (Tregs), and macrophages M0 in high-risk groups (p < 0.05). However, there was a significant decrease in the expression of resting memory CD4 T cells, monocytes, and resting mast cells in high-risk groups (Figure 5A and B). Additionally, all 36 immune checkpoints, with the exception of IDO2 and ICOSLG, were highly expressed in high-risk groups (p < 0.05) (Figure 5C). It is worth noting that the IPS score was lower, whereas the TIDE score was higher in high-risk groups (Figure 5D and E). Based on these findings, high-risk groups experienced higher levels of immune escape as well as poorer responses to immunotherapy. Moreover, the stromal score and immune score were higher in high-risk groups, which indicated a lower tumor purity in these groups (Figure 5F).

Figure 5 Analysis of immune infiltration and immunotherapy between high- and low-risk groups. (A) The proportion of 22 immune cells in LIHC samples. (B) Discrepancies of immune cells between high- and low-risk groups. ns, not significant; *p<0.05; ** p<0.01; ***p<0.001; ****p<0.0001. (C) Comparison of 47 immune checkpoints expression between high- and low-risk groups. ns, not significant; *p<0.05; ** p<0.01; ***p<0.001; ****p<0.0001. (D and E) The discrepancies of immunophenscore (IPS) (D) and tumor immune dysfunction and exclusion (TIDE) score (E) between high- and low-risk groups. ns, not significant; *p<0.05; ** p<0.01. (F) Comparison of immune score, ESTIMATE score, stromal score, and tumor purity between high- and low-risk groups. *p < 0.05, **p < 0.01. Abbreviation: ns, not significant.

Expression Verification of Biomarkers

The expressions of ACADS, FTCD, and GOT2 were significantly lower, while the expressions of CYP19A1 and G6PD were significantly higher in HCC groups of TCGA database (p < 0.05). Furthermore, the expression trends of these biomarkers in GSE14520 were consistent with TCGA database (Figure 6A and B). Besides, the qRT-PCR results also showed that ACADS, FTCD, and GOT2 were extremely significantly lowly expressed, while CYP19A1 was extremely significantly highly expressed in HCC samples (p < 0.001). These results were consistent with the results of our previous analyses in both TCGA-HCC dataset and GSE14520 (Figure 6C).

Figure 6 The expression of five biomarkers in LIHC and normal samples in TCGA dataset (A) and GSE14520 dataset (B). ****p<0.0001. (C) The mRNA expression of biomarkers in normal and HCC samples by real-time reverse transcription PCR (qRT-PCR). ns, not significant; ****p<0.0001.

Discussion

HCC treatment strategies vary from the stage of the disease and the underlying conditions, and patients with similar disease phenotypes might have different molecular etiologies, which could lead to the uniqueness of treatment responses. The categorization of patients according to their molecular characteristics would be beneficial to the construction of the most accurate treatment predicting models.²⁷ Nevertheless, in-depth studies have enhanced our understanding of the role of propionate metabolism in cancer progression.^28–30

In this study, we found out 132 DE-PMRGs associated with steroid and organic acid metabolic process. Kang³¹ also pointed out that the DEGs in HCC mainly included steroid metabolic process, which is consistent with our study. Deng³² proposed a new approach to HCC from the standpoint of the glycolysis and cholesterol synthesis axis and predicted its prognosis. In addition, our study also pointed out that these 132 DE-PMRGs were associated with the degradation of valine, leucine, and isoleucine, and carbon metabolism. The catabolism of amino acids is primarily carried out in the liver, and the amino acid metabolism starts from the tricarboxylic acid cycle and glycolytic pathway. Additionally, cancer cells are characterized by an enhanced glycolysis rate.³³ Wang³⁴ also noted that DEGs between high- and low-risk HCC patients played a key role in cell metabolism, especially organic acids, inorganic acids, and lipid metabolism, which is consistent with our study.

ACADS, CYP19A1, FTCD, G6PD, and GOT2 were screened and identified as biomarkers of HCC. Based on the expression heatmap of biomarkers, it is evident that CYP19A1 and G6PD exhibit higher expression levels in the high-risk group, whereas ACADS, FTCD, and GOT2 show elevated expression in the low-risk group. This observation aligns with the univariate Cox analysis where CYP19A1 and G6PD exhibit hazard ratios (HR) greater than 1, while ACADS, FTCD, and GOT2 consistently display HR values less than 1. In conclusion, CYP19A1 and G6PD are considered risk factors within the context of Hepatocellular Carcinoma (HCC), whereas ACADS, FTCD, and GOT2 are regarded as protective factors.ACADS could be considered a biomarker of HCC, and it was associated with methylation and might be associated with proliferation and metastasis.³⁵ Polymorphisms in CYP19A1 were associated with HCC and promoted prostate stem cell proliferation and metastatic homing.³⁶ As a tumor suppressor gene, FTCD plays a key role in the pathogenesis and progression of HCC in tumor tissues.³⁷ Inhibition of ferroptosis by G6PD was associated with increased proliferation, migration, and invasion of HCC cells.³⁸ Research has demonstrated that GOT2 expression is down-regulated in HCC tissues, and GOT2 knockdown promotes proliferation, migration, and invasion in HCC cells.³⁹ In summary, this study has identified five biomarkers through bioinformatics methods that may be associated with the occurrence and progression of Hepatocellular Carcinoma (HCC). These findings hold the potential to offer novel targets for improving the prognosis and treatment of HCC patients. Naturally, we will continue to monitor the roles of these genes within HCC and further delve into their underlying mechanisms.Immune infiltration results in a weakened tumor immunity, EMT, stromal cell alterations, tumor angiogenesis, and metabolic changes in TME. As a result of these changes, the tumor proliferates, invades, and metastasizes, evading the immune system and resisting treatment.^40–43 An immune cell infiltrating a tumor could influence the microenvironment, and its presence is indicative of the success of immunotherapy. TIDE and ESTIMATE are two computational tools used to predict and infer immune response and potential response to immunotherapy. According to previous studies,⁴⁴ T/B cells were detected in HCC tissue, whereas immune negative regulatory cells from high-risk cohorts were found infiltrating significantly. The high- and low-risk groups exhibited differences in the expressions of 36 immune checkpoints based on a comprehensive analysis of 47 immune checkpoints. A higher level of tumor immune evasion was observed in the high-risk population, which indicated a more adverse outcome. In the following analyses, we found that the cohort that was at a heightened risk had higher TIDE scores, stromal scores, ARG scores, invasion scores, EMT scores, and mRNAsi scores. The high-risk group exhibited a greater degree of immune scores, which indicated a greater level of immune cell infiltration. A significant association has been established between elevated IS and a longer PFS and OS following HCC.⁴⁵ PD-1 (programmed cell death protein 1) is an immune checkpoint that plays a role in down-regulating the immune system by preventing the activation of T cells. CD8+ T cells are a type of immune cell that can kill cancer cells. It has been shown that a high fraction of PD-1 hi CD8+ T cells in TME is positively associated with treatment response and patient survival in cases of cancer treated with PD-1 blockade.⁴⁶ A poor prognosis has been demonstrated in HCC patients who have PD1Hi or TIM3+ PD1Hi CD8+ T cells. Based on preclinical HCC models, a personalized neoantigen vaccine combined with PD-1 blockade increased the infiltration of CD8+ memory T cells.⁴⁷ According to our hypothesis, increased tumor immune attenuation might result from propionate-mediated modifications in the phenotype of immune cells. Additionally, propionate promoted a higher expression of immune checkpoints, modified the transition of tumor epithelial mesenchymal cells, triggered changes in stromal cells, and stimulated tumor angiogenesis. These results suggested that tumors are more invasive and less responsive to treatment in the high-risk group. Nevertheless, the presence of significant T/B cell infiltration and elevated IP scores indicate a greater possibility for the collaborative amplification of antitumor treatment approaches of high-risk HCC patients compared with patients with minimal lymphoid immune cell infiltration. Future research might be able to address the findings of this study.

Several limitations should be considered in this study. TCGA and GSE14520 revealed that most cases did not have a B-related hepatitis history. Studies of this type may exclude patients with hepatitis B-related HCC from other areas with distinct genetic phenotypes and clinical characteristics, making them susceptible to bias. Although these five genes have been annotated through computational approaches, further research is required to elucidate the mechanisms whereby these genes contribute to the development of HCC.

Conclusion

Based on the TCGA and GEO databases, we have identified a gene signature associated with the metabolism of propionate that could be used to classify the prognosis of HCC patients. It might be possible to use this method for the management of HCC patients as well as for individualized treatment planning.