Zheng Yu1, Qifeng Ou1, Fan Chen2, Jiong Bi1, Wen Li1, Jieyi Ma1, Rongchang Wang3, Xiaohui Huang4. 1. Laboratory of Surgery, The First Affiliated Hospital, Sun Yat-Sen University, Guangzhou, 510080, China. 2. Department of Orthopedics Surgery, The First Affiliated Hospital, Sun Yat-Sen University, Guangzhou, 510080, China. 3. Department of Gastrointestinal Surgery, The First Affiliated Hospital, Guangzhou Medical University, Guangzhou, 510120, China. wrc1985wrc@163.com. 4. Laboratory of Surgery, The First Affiliated Hospital, Sun Yat-Sen University, Guangzhou, 510080, China. 15626020913@163.com.
Abstract
BACKGROUND: Hepatocellular carcinoma is a malignant tumor with a highly invasive and metastatic phenotype, and the detection of potential indicators associated with its recurrence and metastasis after surgical resection is critical for patient survival. METHODS: Transcriptome data for large cohorts (n = 1432) from multicenter sources were comprehensively analyzed to explore such potential signatures. The prognostic value of the selected indicators was investigated and discussed, and a comparison with conventional clinicopathological features was performed. A survival predictive nomogram for 5-year survival was established with the selected indicator using the Cox proportional hazards regression. To validate the indicator at the protein level, we performed immunohistochemical staining with paraffin-embedded slides of hepatocellular carcinoma samples (n = 67 patients) from our hospital. Finally, a gene set enrichment analysis (GSEA) was performed to detect the underlying biological processes and internal mechanisms. RESULTS: The liver-specific protein paraoxonase 1 (PON1) was found to be the most relevant indicator of tumor recurrence, invasiveness, and metastasis in the present study, and the downregulation of PON1 might reveal poor survival for patients with hepatocellular carcinoma. The C-index of the PON1-related nomogram was 0.714, thus indicating a more effective predictive performance than the 7th American Joint Committee on Cancer (AJCC) tumor stage (0.534), AJCC T stage (0.565), or alpha-fetoprotein (0.488). The GSEA revealed that PON1 was associated with several hepatocellular carcinoma-related pathways, including the cell cycle, DNA replication, gap junction and p53 downstream pathways. CONCLUSIONS: The downregulation of paraoxonase 1 may suggest worse outcomes and a higher recurrence rate. Thus, paraoxonase 1 might represent an indicator for predicting the survival of patients with hepatocellular carcinoma.
BACKGROUND:Hepatocellular carcinoma is a malignant tumor with a highly invasive and metastatic phenotype, and the detection of potential indicators associated with its recurrence and metastasis after surgical resection is critical for patient survival. METHODS: Transcriptome data for large cohorts (n = 1432) from multicenter sources were comprehensively analyzed to explore such potential signatures. The prognostic value of the selected indicators was investigated and discussed, and a comparison with conventional clinicopathological features was performed. A survival predictive nomogram for 5-year survival was established with the selected indicator using the Cox proportional hazards regression. To validate the indicator at the protein level, we performed immunohistochemical staining with paraffin-embedded slides of hepatocellular carcinoma samples (n = 67 patients) from our hospital. Finally, a gene set enrichment analysis (GSEA) was performed to detect the underlying biological processes and internal mechanisms. RESULTS: The liver-specific protein paraoxonase 1 (PON1) was found to be the most relevant indicator of tumor recurrence, invasiveness, and metastasis in the present study, and the downregulation of PON1 might reveal poor survival for patients with hepatocellular carcinoma. The C-index of the PON1-related nomogram was 0.714, thus indicating a more effective predictive performance than the 7th American Joint Committee on Cancer (AJCC) tumor stage (0.534), AJCC T stage (0.565), or alpha-fetoprotein (0.488). The GSEA revealed that PON1 was associated with several hepatocellular carcinoma-related pathways, including the cell cycle, DNA replication, gap junction and p53 downstream pathways. CONCLUSIONS: The downregulation of paraoxonase 1 may suggest worse outcomes and a higher recurrence rate. Thus, paraoxonase 1 might represent an indicator for predicting the survival of patients with hepatocellular carcinoma.
Recurrence and metastasis after hepatic resection of hepatocellular carcinoma (HCC) usually contributes to poor long-term patient survival [1]. HCC has strong invasiveness and metastasis abilities, thus enhancing its recurrence rate [2]. Although novel therapies have been developed in recent years, the mortality of HCCpatients has not decreased [3]. The identification of potential indicators associated with recurrence and metastasis can improve treatment response in clinical trials and the quality of life of HCCpatients. Transcriptome data from very large cohorts were used in the present study to show that paraoxonase 1 (PON1) might represent a potential signature that is strongly correlated with HCC recurrence and metastasis.Evidence has shown that PON1 is a member of the paraoxonase family, and it encodes a protein of an enzyme with lactonase and ester hydrolase activity [4]. It is an antioxidant defensive factor that is relevant in the pathogenesis of several inflammatory diseases [5]. Increasing evidence has demonstrated that PON1 plays a significant role in atherosclerosis as its downregulation may lead to the disease [6, 7]. In addition, its clinical application value in cancer studies has been gradually discovered [8]. Increasing evidence has demonstrated that PON1 could serve as a significant clinical indicator for breast cancer and lung cancer [9, 10]. Meanwhile, several studies have also found that serum PON1 is highly fucosylated (Fuc-PON1) and could be utilized as a novel diagnostic biomarker of early-stage HCC [11]. However, the number of studies on PON1 in HCC is limited, and its latent prognostic value and potential in clinical applications, especially its correlation with metastasis, recurrence, overall survival (OS) and other clinical risk factors, has not been observed.In the present study, we performed a bioinformatics analysis based on sequencing data of 1432 patients’ tissues to find reliable prognostic signatures. A weighted correlation network analysis (WGCNA) was utilized using The Cancer Genome Atlas (TCGA) for a blind selection in our research. Latent prognostic indicator was selected for further analysis. Expression difference of selected indicator in tumor tissue and nontumor tissue was validated in 7 Gene Expression Omnibus (GEO) datasets from multicenters. We also established a 5-year predictive model that was visualized with a nomogram for predicting 5-year survival. Subsequently, we validated our results on protein levels via experimental tools based on clinical data and our own samples. Finally, we detected and discussed related biological pathways and how they are affected in HCC with PON1 downregulation.
Materials and methods
Transcriptome data of HCC patients in the present study
Data retrieved from multiple research centers were used for integrated analysis in this study, including data from The Cancer Genome Atlas (TCGA) project and Gene Expression Omnibus (GEO) databases. We systematically analyzed the expression profiles of transcriptomes from the following datasets to ensure the credibility of the current study: TCGA [n = 423, National Cancer Institute (NCI) and the National Human Genome Research Institute (NHGRI), USA], GSE14323 (n = 124, Virginia Commonwealth University), GSE14520 (n = 481, National Cancer Institute, Laboratory of Human Carcinogenesis, USA), GSE6764 (n = 75, Mount Sinai School of Medicine), GSE51401 (n = 64, Zhongshan Hospital affiliated with Fudan University), GSE41804 (n = 40, Kanazawa University Graduate School of Medical Sciences), GSE45436 (n = 134, National Yang-Ming University) and GSE62232 (n = 91, INSERM, UMR U-1162, Université Paris Descartes). The sequencing data used in the present study were collected on Affymetrix microarray and Illumina platforms by different researchers. Complete data on the clinicopathologic characteristics of HCCpatients in TCGA were obtained from Cbioportal (http://www.cbioportal.org/), including 7th American Joint Committee on Cancer (AJCC) staging rules, serum alpha-fetoprotein (AFP) level, disease-free time, OS time, living status and recurrence time.
Identification of differentially expressed genes
Increasing evidence has demonstrated that the expression levels of cancer-related genes are abnormally changed during the initiation of HCC. Thus, discovering these genes for further analysis is critical. Based on 373 tumor tissues and 50 adjacent nontumor tissues in TCGA datasets, we identified key differentially expressed genes (DEGs) via R/language (edgeR package, R version 3.34). Varying degrees of gene expression between tumor tissues and adjacent non-tumor tissues were evaluated by log2(fold-change). Genes with abs(log2(fold-change)) > 1 (absolute value) and P < 0.05 were selected as candidate signatures for further analysis.
Weighted correlation network analysis for discovering recurrence-related gene modules
To uncover the recurrence-related DEGs, we conducted a WGCNA using DEGs [12, 13]. In general, genes with similar expression patterns are likely to exhibit coexpressed relationships and similar molecular functions. According to the connectivity of the coexpressed genes, WGCNA clustered thousands of genes into several gene modules using a soft thresholding power [14, 15]. The relations between different gene modules and clinicopathologic characteristics of prognosis were revealed by Pearson correlation. In the present study, we investigated gene modules that were the most relevant to recurrence and metastasis. We performed the WGCNA analysis based on clinical data and sequencing data of 373 HCC samples from TCGA dataset. Samples with the duplicated TCGA barcodes and uncomplete survival data were excluded. In total, 368 HCCpatients were involved in the WGCNA, and the clinical information of these patients was also available for the WGCNA. The soft thresholding power was identified using the WGCNA algorithm and sequencing data of DEGs. Then, DEGs were clustered in different gene modules using the soft thresholding power, and the recurrence-related gene module was investigated.
Indicator selection procedure
The Kaplan–Meier method was performed for every gene in the selected gene module from the WGCNA. The gene with the most significantly different log-rank p value was retained as our final selected indicator, and its expression difference between tumors and nontumors was validated using the GEO datasets. The prognostic value of our candidate indicator was also investigated. The relationship between the clinicopathologic characteristics and expression of the candidate indicator was also investigated. Student’s t-test was utilized for continuous variables, and Fisher’s exact test was used for categorical variables. In total, 369 patients with complete clinicopathologic information were evaluated. Density plots were illustrated for clear views.
Prognostic nomogram with selected indicator
Patients with HCC usually have a poor prognosis due to tumor recurrence. Thus, establishing a recurrence-related model for predicting survival is significant. We established a nomogram to predict 5-year survival. Multivariate Cox proportional hazards regression was used for the nomogram. TCGA datasets of 368 HCCpatients were used [16]. Compared with the traditional method, PON1 expression (log2 transformed, Htseq-counts) was used as a novel variable for the establishment of a nomogram. For other variables, we selected several survival-related indicators and basic information, including age, sex, AJCC staging indicators, tumor differentiation, vascular invasion, and Child–Pugh classification. Five-year survival prediction performance was examined with the C-index. To avoid the potential bias caused by tumor heterogeneity, we performed internal validation by extracting 60% of samples randomly as the validation dataset for 10 times. C-index was calculated for testing the robustness. Calibration curves were also illustrated to test the predicting accuracy. Univariate and multivariate COX proportional hazards tests were performed for conventional clinical features, and the receiver operating characteristic (ROC) curves of conventional features for predicting 5-year survival were also displayed for comparison.
Immunohistochemical staining for validation at protein level
To validate the selected indicator at protein level, we performed immunohistochemical staining in the samples of human tissues. The slides of paraffin-embedded HCC tissues from 67 patients were obtained after surgical resection. All patients received radical resection without preoperative chemotherapy or radiotherapy in the First Affiliate Hospital of Sun-yat Sen University. All diagnoses were confirmed by pathology. The slides were incubated for 2 h at 65 °C, deparaffinized, and rehydrated. Retrieval of the heat-mediated antigen was conducted in 10 mmol/L Tris-citrate buffer (pH 7.0) with a pressure cooker. Blocking of endogenous peroxidase activity was performed by incubating sections with 3% hydrogen peroxide for 10 min at room temperature. Each section was then incubated with 5% normal goat serum in phosphate-buffered saline containing 0.1% Tween 20 for 30 min at room temperature to block nonspecific binding of the primary antibody. The slides were incubated with primary antibodies (diluted 1:250) against PON1 (Abcam, ab92466) overnight at 4 °C. After washing, each slide was incubated with the appropriate horseradish peroxidase (HRP)-labeled secondary antibody and then developed with DAB solution (DAKO, Agilent) before counterstaining with hematoxylin. Staining intensity was scored as 0, 1, 2, or 3 for negative, weak, moderate, or strong, respectively, and the staining percentage was given a score of 0 (absent) for < 5% positive staining, 1 (focal) for 5% to < 50% positive staining, or 2 (diffuse) for ≥ 50% positive staining. The sum of the intensity and distribution scores was then used to determine PON1 immunoreactivity. A score of 1 or 0 was considered to show low expression, whereas higher scores were considered to indicate high expression. Two pathologists assessed the specimens independently. Images were obtained using an Olympus BX63 microscope (Olympus, Japan).
Biological pathway and internal mechanism detection
To investigate the pathways and biological parameters correlated with PON1, we conducted a gene set enrichment analysis (GSEA) [17, 18]. The 368 HCCpatients in the TCGA dataset were separated into high-expression and low-expression groups according to the median value of PON1 and retained as the phenotypes. GSEA software was obtained from GSEA website (http://software.broadinstitute.org/gsea/index.jsp). One thousand permutations were performed to determine the statistically significant and ensure the credibility of the results. Pathways with a FDR < 0.05 and P-value < 0.05 were selected as the enriched terms. The molecular signatures database from the GSEA were utilized as the annotation file. Pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome and Pathway Interaction Database (PID) were used. A Gene Ontology enrichment analysis was also performed to discover biological functions in HCC. Correlated biological processes, molecular functions, and cellular components were identified.
Statistical software and data format
DEGs were identified using edgeR (Bioconductor/R version 3.34) [19, 20]. Htseq-counts (TCGA) and RSEM (GEO, log2 scaled) were used as the sequencing data formats. WGCNA was performed by the WGCNA package (R version 3.3.4). Nomogram was established using Regression Modeling Strategies (RMS, R version 3.34). Student’s t-test and Fisher’s exact test were performed in R (version 3.34).
Results
A total of 985 DEGs were discovered by bioinformatics analysis, including 512 upregulated and 473 downregulated genes (Additional file 1: Table S1). The variation degree of DEGs ranged from 8 to − 5 as evaluated by log2(fold-change). A volcano plot was created to illustrate the distribution of DEGs (Fig. 1a).
Fig. 1
Identifying recurrence-related prognostic indicators. a Volcano plot exhibiting the distribution of DEGs via log2(fold-change) and P values. b Sample clusters showing basic clinical information on HCC patients. c Soft threshold power selection for gene clustering; a value of 7 was selected as the threshold. d In total, 13 gene modules were discovered via gene clustering. e Heatmap exhibiting the relationship between different gene modules and clinical risk factors (Pearson correlation). Blue modules were the most relevant to recurrence and OS. f Plot showing 34 genes in blue modules and their coexpression relationships
Identifying recurrence-related prognostic indicators. a Volcano plot exhibiting the distribution of DEGs via log2(fold-change) and P values. b Sample clusters showing basic clinical information on HCCpatients. c Soft threshold power selection for gene clustering; a value of 7 was selected as the threshold. d In total, 13 gene modules were discovered via gene clustering. e Heatmap exhibiting the relationship between different gene modules and clinical risk factors (Pearson correlation). Blue modules were the most relevant to recurrence and OS. f Plot showing 34 genes in blue modules and their coexpression relationships
WGCNA analysis revealed blue gene module was related to tumor recurrence
Based on the complete clinical information and RNA-seq data of 985 DEGs (Htseq-count), we performed WGCNA analysis to cluster DEGs into different modules according to their coexpression relationships. We also revealed different gene modules and their relevant clinical information. Patients were first clustered to show basic clinical information (Fig. 1b). Then, a soft threshold power was calculated for gene clustering, and 7 was selected as the power in the present study (Fig. 1c). In total, 13 gene modules were discovered (Fig. 1d). The relationships between gene modules and primary clinical indicators were discovered with Pearson correlation. Correlation coefficients (r) were illustrated as a heatmap (Fig. 1e). Blue modules were positively correlated with OS time (r = 0.18) and recurrence time (r = 0.35). In addition, they were also negatively correlated with AJCC tumor staging rules, including T stage (r = − 0.20) and tumor stage (r = − 0.20). The results demonstrated that genes in the blue modules could be viewed as tumor suppressors. Patients with higher expression of genes in blue modules had longer survival times. A total of 34 genes were found in blue modules (Fig. 1f). Survival analysis was performed with genes in blue module. PON1 was the most relevant gene for OS using the Kaplan–Meier (KM) method according to log-rank P values (Additional file 2: Table S2).
Reduction in PON1 expression indicated greater invasiveness and a poor prognosis
We detected the expression of PON1 in other GEO datasets and found that PON1 expression in HCC was downregulated in nearly all datasets (Fig. 2a, Additional file 3: Table S3). Patients in the TCGA dataset were separated into two groups according the median value of PON1, including a high-expression group and a low-expression group. Differences between the two groups were investigated with Student’s t-test and Fisher’s exact test (Table 1). The results revealed that significant differences existed in age, sex, AJCC TNM staging rules, tumor differentiation, and vascular invasion. We compared the diagnosis performance of key clinical indicators in HCC, including marker Ki-67 (MKI67), AFP, and PON1 [21, 22]. The results indicated that the area under curve (AUC) of PON1 was 0.8119, which was higher than that of AFP (0.6857) but lower than that of MKI67 (0.9515) (Fig. 2b) [23]. Based on a density plot of PON1 expression, we clearly observed that low PON1 expression in HCC was correlated with stronger invasiveness and metastasis, including tumor T stage, AJCC tumor stage, tumor differentiation and vascular invasion (Fig. 2c–f) [24]. The KM curves of Disease-free survival and overall survival was present for visualizing survival data, whilst the Log Rank test was utilized for determining differences. Results revealed Log Rank P-values were 0.013 (DFS) and P = 0.0014 (OS) separately (Fig. 2g, h). HCCpatients with low PON1 expression were found to have a poor prognosis in long-term survival.
Fig. 2
Prognostic value of PON1 in clinical applications. a Validation of PON1 expression differences in tumor tissues and adjacent non-tumor tissues by analyzing multicenter data sources that included 1432 samples. PON1 was downregulated in HCC tissues. b Comparison of the diagnostic ability of PON1, MKI67, and AFP via ROC curves. The diagnostic accuracy of PON1 (0.8119) was higher than that of AFP (0.6857) and lower than that of MKI67 (0.9515). c–f Density plot exhibiting the relationship between PON1 expression and clinicopathologic characteristics, including AJCC tumor stage, tumor T stage, tumor differentiation, and vascular invasion. The trend in the peaks indicated that low tumor expression was likely to be classified as a poor prognosis, including stage III, G4, macro vascular invasion, and T4. g Disease-free survival analysis of patients with low PON1 expression and high PON1 expression in the TCGA dataset (P-value = 0.013). h OS analysis of patients with low PON1 expression and high PON1 expression in TCGA dataset (P = 0.0014)
Table 1
Association of PON1 with clinicopathological characteristics of HCC
Clinical factors
PON1 level
Total (n = 369)
95% CI
P-value
Group
High (n = 184)
Low (n = 185)
Gender
F
40
81
121
0.2198819
0.575337
8.10E−06
M
144
104
248
Age (mean, SD)
61.17
57.64
0.7774031
6.287814
0.01212
Tumor pathologic_stage
Stage I–II
145
129
256
1.140819
3.267362
0.009873
Stage III–IV
36
56
89
Tumor size
T1–T2
145
129
274
1.052772
2.921888
0.02976
T3–T4
36
56
92
Tumor metastasis
M0
129
136
265
0.2244094
150.4063
0.623
M1
1
3
4
Tumor nodes
N0
120
130
250
0.2181723
146.4568
0.6237
N1
1
3
4
Child–Pugh classification
A
115
101
216
0.1671
B
15
6
21
C
1
0
1
Grade
G1
36
19
55
9.62E−05
G2
99
77
176
G3
42
79
121
G4
4
8
12
Race
American indian or alaska tive
1
1
2
0.1536
Asian
68
90
158
Black or african american
10
7
17
White
98
84
182
Tumor status
Tumor free
119
114
233
0.7075
1.8518
0.6426
With tumor
52
57
109
Vascular invasion
Macro
6
9
15
0.00552
Micro
35
57
92
Non-vascular invasion
118
88
206
Hepatic_infalmmation_adj_tissue
Mild
44
55
99
0.02331
Severe
10
8
18
None
73
43
116
AFP (mean, SD)
16352.55
11481.7
0.7408
Prognostic value of PON1 in clinical applications. a Validation of PON1 expression differences in tumor tissues and adjacent non-tumor tissues by analyzing multicenter data sources that included 1432 samples. PON1 was downregulated in HCC tissues. b Comparison of the diagnostic ability of PON1, MKI67, and AFP via ROC curves. The diagnostic accuracy of PON1 (0.8119) was higher than that of AFP (0.6857) and lower than that of MKI67 (0.9515). c–f Density plot exhibiting the relationship between PON1 expression and clinicopathologic characteristics, including AJCC tumor stage, tumor T stage, tumor differentiation, and vascular invasion. The trend in the peaks indicated that low tumor expression was likely to be classified as a poor prognosis, including stage III, G4, macro vascular invasion, and T4. g Disease-free survival analysis of patients with low PON1 expression and high PON1 expression in the TCGA dataset (P-value = 0.013). h OS analysis of patients with low PON1 expression and high PON1 expression in TCGA dataset (P = 0.0014)Association of PON1 with clinicopathological characteristics of HCC
PON1-related prognostic nomogram
A prognostic nomogram was established between PON1 and several significant clinical factors, including age, sex, vascular invasion, Child–Pugh classification, AJCC staging rules, and tumor differentiation (Fig. 3a). The calibration curves indicated that the predictive performance of the model was excellent (Fig. 3b). In our nomogram, all variables were in accordance with clinical logic. PON1 was utilized as a new variable and promoted model accuracy. The expression of PON1 was negatively correlated with risk score. The C-index (an evaluating marker similar to the ROC curve) of the model was 0.714 (0.6753–0.7627, C-index in validation datasets). The results of the univariate Cox hazards analysis revealed that the AJCC tumor pathological stage, tumor size, tumor metastasis, and AFP level were related to 5-year survival. Hazard ratios were listed (Additional file 4: Table S4). As a comparison with the nomogram, we illustrated the ROC curves using the AJCC tumor stage, tumor size, and AFP in 5-year survival prediction (Fig. 3c). The results indicated that our nomogram provided a more accurate prediction than using conventional clinical features.
Fig. 3
Establishment of the prognostic nomogram for 5-year survival. a Nomogram for predicting survival probabilities. Each clinical characteristic is shown with its corresponding score. PON1 was used as a novel liver-specific variable in the nomogram. The C-index of the nomogram was 0.714. b Calibration plot of OS at 5 years to visualize the difference between true values and predicted values. c ROC curves of AJCC tumor stage (0.5349), serum AFP (0.4888), and tumor T stage (0.5656) for predicting 5-year survival are shown for comparison
Establishment of the prognostic nomogram for 5-year survival. a Nomogram for predicting survival probabilities. Each clinical characteristic is shown with its corresponding score. PON1 was used as a novel liver-specific variable in the nomogram. The C-index of the nomogram was 0.714. b Calibration plot of OS at 5 years to visualize the difference between true values and predicted values. c ROC curves of AJCC tumor stage (0.5349), serum AFP (0.4888), and tumor T stage (0.5656) for predicting 5-year survival are shown for comparison
Experimental validation of PON1 at the protein level by immunohistochemistry
To validate these results, we conducted an immunohistochemical analysis to reveal the protein level of PON1 in 67 HCCpatients to explore the relationship between PON1 and HCC clinicopathological characteristics. PON1 was found to be significantly located in the cytoplasm, and it was downregulated in HCC tissues. The junction area was illustrated for the comparison. The statistical results indicated that 42 (62.7%) of the 67 samples showed high expression levels according to staining intensity (Fig. 4a). Survival data were also available for 67 patients. KM curves of the OS and disease-free survival was presented for visualizing survival data, whilst the Log Rank test was performed for determining differences. Result revealed that patients with low expression of PON1 in HCC tissue had a poor prognosis for overall survival (Fig. 4b, c).
Fig. 4
Protein level validation using immunohistochemistry (IHC). PON1 protein expression in adjacent non-tumor tissues and tumor tissues. a PON1 location in HCC was demonstrated, particularly in the cytoplasm. Low and high expression in tumor tissues was also illustrated. The junction area is shown as a comparison. b Disease-free survival based on the IHC of 67 patients using the KM method (P = 0.14). c OS based on the IHC of 67 patients using the KM method (P = 0.043)
Protein level validation using immunohistochemistry (IHC). PON1 protein expression in adjacent non-tumor tissues and tumor tissues. a PON1 location in HCC was demonstrated, particularly in the cytoplasm. Low and high expression in tumor tissues was also illustrated. The junction area is shown as a comparison. b Disease-free survival based on the IHC of 67 patients using the KM method (P = 0.14). c OS based on the IHC of 67 patients using the KM method (P = 0.043)
Correlated biological pathways of PON1 in HCC
The GSEA was used to detect pathways that were correlated with PON1. HCC-related pathways were found to be related to PON1, including the cell cycle (P = 0.00122), DNA replication (P = 0.002093), gap junction (P = 0.01286), and p53 downstream pathways (P = 0.00252) (Fig. 5a, Additional file 5: Table S5). Researchers have found that the cell cycle pathway is abnormally changed from normal liver functions to chronic hepatitis as well as during the transition into HCC [25]. Increasing evidence has revealed that gap junction pathways could affect HCC invasion and metastasis [26]. Apoptosis was also discovered to be related to PON1 (P = 0.008673). As a reference, peroxisome (P = 0.0001) and biological oxidation pathways (P = 0.0001) were also identified (Fig. 5b, c) [27]. These findings demonstrated that the results of the GSEA were reliable. Enrichment of Gene Ontology terms was also conducted and illustrated. The results indicated that PON1 variation in HCC leads to changes in oxidation reduction processes (P < 0.001), oxygen binding (P < 0.001), extracellular exosomes (P < 0.001), and blood microparticles (P < 0.001) (Additional file 6: Table S6).
Fig. 5
Detection of biological pathways and internal mechanisms. a Plot of the GSEA; several key pathways are visualized, including the cell cycle, DNA replication, gap junction, and p53 downstream pathways. b Enriched pathways found by the GSEA using MsigDB. c Plot of the Gene Ontology enrichment analysis. Biological processes (BP), molecular functions (MF) and cellular components (CC) are illustrated
Detection of biological pathways and internal mechanisms. a Plot of the GSEA; several key pathways are visualized, including the cell cycle, DNA replication, gap junction, and p53 downstream pathways. b Enriched pathways found by the GSEA using MsigDB. c Plot of the Gene Ontology enrichment analysis. Biological processes (BP), molecular functions (MF) and cellular components (CC) are illustrated
Discussion
Although many HCC-related prognostic biomarkers have been identified in recent years, most do not exhibit tissue specificity [28]. Thus, these biomarkers might be affected by various factors, and they lack significant value in clinical applications. We performed a systematic analysis that discovered the prognostic value of PON1 in patients with HCC. We found that PON1 downregulation in HCC suggests worse tumor differentiation, higher recurrence rate, stronger invasiveness, and poorer outcomes. Researchers in recent years have demonstrated that serum PON1 levels can assist in diagnosing AFP-negative HCC at an early stage [29, 30]. Moreover, PON1 played an important role in the initiation of non-alcoholic fatty liver disease (NAFLD) [31, 32]. Meanwhile, NAFLD was found to be a potential risk factor for HCC, especially in patients without hepatitis virus infection [33, 34]. However, discoveries about PON1 in HCC were mostly based on small cohorts from a single data source in retrospective studies, which ignored diversity in terms of race, age, and hepatitis virus infection. The prognostic application of PON1 was also rarely investigated or discussed. Finding crucial cancer biomarkers from thousands of genes is usually difficult. In the present study, we introduced WGCNA algorithm for discovering recurrence-related indicators. Compared with other algorithms, WGCNA systematically combined the sequencing data and clinical information. Using of WGCNA could assist in discovering underlying clinical significance of some key genes. WGCNA also illustrated the connection of genes with similar expression patterns. Moreover, WGCNA algorithm was gradually used in analyzing the single-cell transcriptome data by researchers [35]. Meanwhile, our research was conducted on 1432 samples from multicenter data sources, and the prognostic and predictive values were validated at both the gene expression and protein level. We also detected possible internal mechanisms and biological processes related to HCC. Therefore, we promoted other researchers’ studies and provided more credible results.Even though the expression of PON1 is negatively correlated with tumor recurrence, metastasis, and invasion, we did not define it as a tumor suppressor. The PON1 gene encodes an enzyme that can be released from normal liver cells into the blood circulation and is of great antioxidant significance [36, 37]. According to our pathway detection results, PON1 downregulation may not directly affect the invasiveness and metastasis of tumor cells as PON1 is located downstream of cancer-related biological processes. In our view, PON1 downregulation in HCC might be induced by losing the original function of normal liver cells in PON1 secretion. In GSEA analysis, the results indicated that several carcinogenesis-related pathways were enriched by varying degrees by PON1, including the p53 downstream pathway, gap junction, cell cycle, apoptosis, and DNA replication. Therefore, PON1 downregulation might be increased by complicated internal mechanisms related to cell division, proliferation, and migration, which could also explain why the expression of PON1 was related to tumor recurrence and clinical outcomes. However, a decrease in PON1 might still cause changes in tumor cells, especially tumor-derived inflammation, autophagy, and apoptosis [38, 39]. However, we still need more experimental evidence to prove our conclusions from sequencing data and pathological results. In addition, further detection of these mechanisms will be very necessary.In recent years, next-generation sequencing (NGS) use has gradually increased [40, 41]. However, NGS in tumors is still based on tissue obtained from surgical resection, which limits the clinical application of this method. Although some novel test methods, including sequencing of tumor-derived circulating DNA and exosomes in the blood, have been discovered, their practicability still needs to be validated [42-44]. Noninvasive methods, including blood tests, are still very important. Published research has shown that PON1 is an enzyme that is mainly synthesized in the liver and released into the circulatory system [45, 46]. The expression of PON1 is mainly affected by liver cells, and the use of PON1 in a predictive model may reduce other possible interference factors. Therefore, we built a PON1-related nomogram with a Cox proportional hazards regression to investigate 5-year patient survival. In our model, PON1 gene expression was quantified and served as an indicator to promote its specificity. Compared with conventional clinical indicators, our nomogram exhibited excellent prediction accuracy and effectiveness. However, the limitation of our model is that we did not examine serum PON1 levels or use them as a variable. The model will be more useful if serum PON1 levels are tested and used as a variable.
Conclusions
In summary, we performed a comprehensive analysis of the prognostic value of PON1. We discovered that PON1 downregulation indicates a high recurrence rate and poor outcomes. We also provided a nomogram to use PON1 in clinical applications. We supplied a more accurate plan than conventional methods for predicting prognosis.Additional file 1. 985 differentially expressed genes in HCC.Additional file 2. Log-rank P values of genes in blue module.Additional file 3. Expression values of PON1 in 8 HCC datasets from TCGA and GEO databases.Additional file 4. Univariate and multivariate Cox proportional-hazard regression analysis for overall survivalin TCGA dataset (HCC).Additional file 5. GSEA for MSigDB pathways (KEGG, PID, REACTOME)Additional file 6. Gene ontology (GO) enrichment analysis.
Authors: Aravind Subramanian; Pablo Tamayo; Vamsi K Mootha; Sayan Mukherjee; Benjamin L Ebert; Michael A Gillette; Amanda Paulovich; Scott L Pomeroy; Todd R Golub; Eric S Lander; Jill P Mesirov Journal: Proc Natl Acad Sci U S A Date: 2005-09-30 Impact factor: 11.205
Authors: Jeremy A Miller; Chaochao Cai; Peter Langfelder; Daniel H Geschwind; Sunil M Kurian; Daniel R Salomon; Steve Horvath Journal: BMC Bioinformatics Date: 2011-08-04 Impact factor: 3.169
Authors: Chi Ma; Aparna H Kesarwala; Tobias Eggert; José Medina-Echeverz; David E Kleiner; Ping Jin; David F Stroncek; Masaki Terabe; Veena Kapoor; Mei ElGindi; Miaojun Han; Angela M Thornton; Haibo Zhang; Michèle Egger; Ji Luo; Dean W Felsher; Daniel W McVicar; Achim Weber; Mathias Heikenwalder; Tim F Greten Journal: Nature Date: 2016-03-02 Impact factor: 49.962
Authors: Sara Cuylen; Claudia Blaukopf; Antonio Z Politi; Thomas Müller-Reichert; Beate Neumann; Ina Poser; Jan Ellenberg; Anthony A Hyman; Daniel W Gerlich Journal: Nature Date: 2016-06-29 Impact factor: 49.962
Authors: Alicja E Grzegorzewska; Adrianna Mostowska; Wojciech Warchoł; Paweł P Jagodziński Journal: BMC Infect Dis Date: 2021-08-26 Impact factor: 3.090