Role of CD5L and SRD5A2 as Prognostic Biomarkers for Hepatocellular Carcinoma

Introduction

Hepatocellular carcinoma (HCC) is the sixth most common malignant disease and the second leading cause of cancer-related mortality in males.¹ Despite the improvements in the clinical diagnosis and treatment strategies, the 5-year overall survival (OS) rate of HCC patients remains less than 40%.^1,2 Intrahepatic recurrence and extrahepatic metastasis are the major causes of death.² The identification of specific biomarkers and responsive therapeutic targets could significantly help to predict early recurrence and subsequently improve the prognosis of HCC patients.³ High throughput gene expression profiling allows identification of potential biomarkers for constructing prognosis model and prediction of latent molecular mechanisms. In the last decade, there have been few reports of the combined use of multiple genes for prognostication of HCC.^3–5 However, there is a clear need to identify a more sensitive and specific prognostic biomarker for clinical practice.

CD5L is highly expressed in macrophages and acts as regulator of lipid synthesis, immune homeostasis and inflammatory response.^6–10 CD5L could inhibit chronic liver injury by controlling SMAD7 expression and TGFß signal pathway. But, representative analysis found CD5L protein had higher expression in HCC tissues than that in adjacent liver tissues, and the patients with high expression of CD5L had terrible prognosis comparing to those with under expression.¹¹ In vitro study indicated CD5L prompted the proliferation and colony formation of HCC cells and helped them to escape from cisplatin induced apoptosis.¹¹ The aim of the current study was to investigate the expression and clinical significance of CD5L and other differentially expressed genes (DEGs) in HCC using public databases.

Materials and Methods

Patients and Gene Expression Data

Gene expression and clinical data were obtained from the gene expression omnibus (GEO) database of the National Center for Biotechnology Information (NCBI) (accession numbers GSE45267, ¹² GSE62232,⁴ GSE121248,¹³ and GSE14520-GPL571).¹⁴ The data on the effect of gene expression on prognosis in liver cancer were acquired from the Kaplan-Meier plotter (KM plotter, http://www.kmplot.com/). Kaplan-Meier plotter is a public database containing 54,675 genes of 18,674 cancer patients, including 364 liver cancer cases with relapse-free and overall survival data derived from TCGA database.

Analysis of Differentially Expressed Genes (DEGs)

DEGs in HCC and precancerous samples were screened using the GEO2R online tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/) through which we compared two or more groups of samples in a GEO dataset in order to identify genes that are differentially expressed across experimental conditions. The GEO2R online analysis tool was used to screen the DEGs in four datasets which detected the expression profile of paired samples from HCC and adjacent non-cancerous tissues. To take care of false-positive results, an adjusted p-value was calculated using the Benjamini and Hochberg false discovery rate. An adjusted P-value <0.01, P-value <0.01, and a log fold change (log FC) ≥2 were considered as the thresholds for DEGs screening. Gene with more than one probe set was average. The common DEGs among the four datasets were analyzed using a Venn diagram. The data processing has been shown in Figure S1.

Gene Function Annotation and PPI Network Assay

The database for annotation, visualization and integrated discovery (DAVID) (https://david.ncifcrf.gov/), Omicsbean (http://www.omicsbean.cn/) and KEGG Orthology-Based Annotation System (KOBAS) (http://kobas.cbi.pku.edu.cn/kobas3/) were used to determine the function enrichment for biological process (BP), molecular function (MF), and cellular component (CC), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of DEGs. P<0.05 was considered to be statistically significant. The protein-protein interaction network analysis of DEGs was conducted using the String online database (http://www.string-db.org/, version 11.1). A combined score of more than 0.4 was set as the cut-off value.

Screening for Survival-Related DEGs

The univariate survival analysis of overlapped DEGs in four datasets was performed according to the K-M method using the KM plotter (http://kmplot.com/). The screening filters for gene selection were as follows: multiple genes (gene symbols), split patients by median value, threshold overall survival (OS) of 60 months, median survival, censored at threshold, the pathology and clinical features of the patients were set as default all (stage, grade, AJCC T, vascular invasion, gender, race, sorafenib treatment, alcohol consumption, hepatitis virus were set all). P< 0.05 was used for the purpose of analyzing significantly different. The data processing has been shown in Figure S1.

Statistical Analysis

Statistical analysis was performed using the EmpowerStats, Graph Pad Prism 7.0, IBM SPSS 23.0, and R 3.5.3 software packages. Continuous variables were assayed by two-sided Student’s t-test, or one-way analysis of variance (ANOVA), while nonparametric or χ² test (Fisher exact probability method) were used for categorical variables. Normally distributed data are presented as means ± standard deviation, and the data with skewed distribution are presented as median. Factors that appeared statistically significant on univariate analysis were subjected to multivariate Cox regression analysis. P <0.05 was considered to be statistically significant. The R 3.5.3 software was used to construct the nomogram and the calibration curve.

Results

Identification of DEGs with GEO2R and Analysis of Genetic Characteristics

On analysis, 241, 346, 200 and 124 DEGs were identified in GSE14520, GSE45267, GSE62232, and GSE121248 datasets, respectively (Figure 1). The four datasets contained 47 overlapped DEGs. Out of these, 81% of DEGs (38/47) were downregulated, and 19% of DEGs (9/47) were upregulated. The 47 DEGs were subsequently analyzed. The TOP2A, CDKN3, and IGFBP3 genes were over-expressed, and CXCR12, CXCR14, CD5L, LCAT were under-expressed (Table 1). The analysis of genetic features of the DEGs indicated that most of them were located in the chromosome regions 5q, 6p26, 11p, and 5p (Figure 2A) and distributed in the liver and CD105-positive endothelial cells (Figure 2B). Downregulated DEGs were strongly linked to HCC (Figure 2C), while the upregulated DEGs were likely associated with the pathogenesis of non-alcoholic steatohepatitis (NASH) and breast cancer (Figure 2D).

Table 1 Forty-Seven DEGs

Figure 1 Venn diagram analysis of DEGs in different GEO datasets. Individual studies are indicated in different colors. The overlapping parts indicate DEGs common to different GSE datasets. Together, 47 DEGs were screened from four studies. Abbreviation: DEGs, differentially expressed genes.

Figure 2 The GO analysis of the 47 identified DEGs performed with the enrichR software. (A) chromosome distribution of DEGs; (B) anatomic distribution of DEGs; (C) down-regulated expression of DEGs in HCC; (D) up-regulated expression of DEGs in HCC. Abbreviation: DEG, differentially expressed gene; GO gene ontology.

Functional Classification of DEGs and Protein-Protein Interaction Network

Using the Gene Ontology (GO) annotation, the number of DEGs in each GO term, including CC, BP, MF and corresponding KEGG pathways were also identified. The DEGs were predominantly involved in protein binding, monooxygenase activity, oxidoreductase activity, and the reduction of molecular oxygen (Figure S2). The BP included steroid metabolism, lipid metabolism, monocarboxylic acid metabolism, and response to chemicals (Figure S2). The KEGG pathways included retinol metabolism, bile secretion, tryptophan metabolism, cytokine-cytokine receptor interaction, and the p53 signaling pathway, as shown in Table S1. To understand the interactions involving DEGs, PPI, (Figure 3), network was generated using String 11.1. The detected interactions indicated that TOP2A and CDKN3 were the key nodes in the network. They participated in the regulation of cell cycle and were related to carcinogenesis. TOP2A was also involved in metabolism of protein and CDKN3 was associated with HCC.^15–18

Figure 3 Analysis of PPI network of the 47 identified DEGs using Sting 11.1. Abbreviations: DEG, differentially expressed gene; PPI, protein-protein interaction.

Identification of the Prognosis-Related DEGs Through KM Plotter Analysis (TCGA Cohort)

The K-M method was used in the present data to gain insight into the association between the DEGs and the prognosis (5-year OS, RFS) of patients with HCC. A total of 47 DEGs were used to analyze their association with 5-year OS and RFS in 364 and 313 patients from TCGA database, respectively. Twenty-six genes (BBOX1, ZG16, ASPM, CD5L, TOP2A, DTL, ESR1, FCN2, NAT2, LCAT, GLS2, GHR, as seen in Table 2) were associated with 5-year OS, and 40 genes were linked to RFS in HCC patients (Table 2). Twenty-two genes related to both OS and RFS were analyzed further using a regression model, as detailed below.

Table 2 Validation of DEGs Correlated to RFS and 5-Year OS with K-M Plotter

Validation of Candidate Genes and Identification of Independent Prognostic Factors (GSE14520 Cohort)

To confirm the relationship between the DEGs and prognosis, validation EGs was performed using different cases from GSE14520 cohort. The clinicopathologic features and 22 genes linked to 5-year OS were analyzed by univariate Cox regression, and the results are reported in Table 3. Tumor stage (AJCC TNM stage, BCLC stage, CLIP stage), risk signature for predicted metastasis (RSPM), tumor size, multiple nodular, cirrhosis, AFP, ZG16, ASPM, ACSM3, CD5L, PRC1, DNASE1L3, NAT2, RACGAP1, TOP2A, GLYAT, SRD5A2, SLC22A1, LCAT, ESR1, and GHR were associated with 5-year OS. Moreover, tumor stage (AJCC TNM stage, BCLC stage), RSPM, tumor size, gender, TOP2A, SLC22A1, LCAT, SRD5A2, CRHBP, GLYAT, NAT2, DNASE1L3, ZG16, and LPA were linked to RFS.

Table 3 Clinicopathologic Features and Their Correlation with 5-Year OS and RFS on Univariate Analysis

Based on statistical significance, the BCLC stage was selected for OS and RFS analysis. The multivariate COX regression hazard analysis demonstrated that RSPM, cirrhosis, BCLC stage, and CD5L were independent prognostic factors for OS, while gender, BCLC stage, and SRD5A2 were risk factors for RFS, as shown in Figure 4.

Figure 4 Forest plots summarizing the analysis of OS and RFS. Multivariate analysis of OS (A) and RFS (B) in HCC patients. The tangle squares on the transverse lines indicate the HR, and the transverse lines represent 95% CI. Abbreviations: OS, overall survival; RFS, relapse-free survival; HCC, hepatocellular carcinoma; CI, confidence interval.

Construction and Validation of a DEGs-Based Predictive Nomogram Model Using GSE14520 Cohort

The cases of GSE14520 were divided into the training (n=144) and validation cohort (n=98). The validation cohort was used for the external validation of the nomogram model. The stage of cirrhosis (hazard ratio [HR] 4.97, 95% confidence interval [95% CI] 1.17–21.18, P=0.028), RSPM (HR 11.73, 95% CI 1.63–84.58, P=0.015), BCLC (HR 0.55, 95% CI 0.27–1.13, P<0.001), and CD5L (HR 0.712748, 95% CI 0.54–0.93, P=0.036) were found to be independent risk factors for 5-year OS (Figure 4A). BCLC stage (HR 4.48, 95% CI 1.66–7.04, P=0.005), gender (HR 0.42, 95% CI 0.22–0.81, P=0.031), and SRD5A2 (HR 0.87, 95% CI 0.77–0.97, P=0.044) were risk factors for RFS (Figure 4B). The independent prognostic factors were incorporated into a nomogram to estimate the OS and RFS, respectively (Figure 5A and D). The performance of the nomogram was assessed by the index of concordance (C-index). The calibration plot displays ideal consistency for 5-y OS using bootstrap sampling, with a C-index of 0.73 (95% CI 0.71–0.90) in the training cohort and 0.72 (95% CI:0.70–0.89) in the validation cohort, respectively (Figure 5B and C). The C-index for the nomogram predicting 1-yr RFS was 0.67 (95% CI 0.68–0.87) in the training cohort and 0.68 (95% CI 0.68–0.89) in the validation cohort, respectively (Figures 5E and F).

Figure 5 Nomogram to estimate the prognosis with external /internal validation. To use the nomogram, find the position of each variable on the corresponding axis and determine the number of points for each variable. The corresponding points of all the variables are added, and the probabilities of the survival are determined using the lines at the bottom of the nomogram. (A) nomogram for 5-year OS. (B) Internal calibration plot showing the performance of the proposed nomogram in predicting the 5-year OS in the training cohort (n=144). (C) External validation plot showing the predictive performance of the nomogram in estimating the 5-year OS in the validation cohort (n =99). The closer the blue curve is to the red line, the better is the performance. (D) nomogram for 1-year RFS. (E) the internal calibration plot shows the performance of the proposed nomogram in predicting 1-year RFS in the training cohort (n=144). (F) the external validation plot showing the predictive performance of the nomogram in estimating the 1-year RFS in the validation cohort (n=99). Abbreviations: OS, overall survival; RFS, relapse-free survival.

Correlations Between CD5L, SRD5A2 and Clinicopathological Features (GSE14520)

The correlation between CD5L and SRD5A2 with the clinicopathological features was investigated, and the results indicated that CD5L was related to the risk index for predicted metastasis (RIPM, P=0.015), TNM stage (P=0.05), but not to gender (P=0.336), cirrhosis (P=0.811), multiple nodularity (P=0.06), ALT (P=0.296), HBV (P=0.329), tumor size (P=0.827), BCLC stage (P=0.186), CLIP stage (P=0.633), and AFP (P=0.197). On the other hand, SRD5A2 expression was associated with RIPM (P<0.001), TNM stage (P=0.012), BCLC stage (P=0.011), CLIP stage (P<0.001), and AFP (P<0.001), but not with gender (P=0.178), cirrhosis (P=0.473), multiple nodularity (P=0.118), ALT (P=0.117), HBV (P=0.308), and tumor size (P=0.068). Lastly, all 22 genes correlated with CD5L and SRD5A2 (Table 4).

Table 4 Correlation Between CD5L and SRD5A2 Expression and Clinicopathologic Features

Discussion

The pathogenesis of HCC is related to the dysregulation of oncogenes and tumor suppressor genes, which are also responsible for the tumor progression and metastasis. In the current study, we identified 47 DEGs in HCC tissues, most of which were located on chromosomes 5q and 6q, and expressed mostly in the liver and CD105-positive endothelial cells. Previous studies also had similar observations with the tumor suppressor gene of non-cirrhotic HCC located on chromosome 5q and the gene related to recurrence located on chromosome 6q.^19–22 In this study, most of DEGs with low expression were located on chromosome 5q, suggesting that the decreased expression of these genes promoted hepatocarcinogenesis.²¹ The upregulated genes were related to NASH and breast cancer, while the downregulated genes were closely associated with HCC. The GO and KEGG pathway analysis revealed that the 47 aberrant genes were involved in pathways of steroid, lipid, and retinol metabolism, bile secretion, cytokine-cytokine interaction, and p53 signalling, which are implicated in the cell cycle progression, proliferation, and invasion of cancer cells. The above biological phenotypes were in conformity to the role of the hub genes (TOP2A, ESR1, CDKN3, and PRC1).

In this study, 26 and 40 genes related to 5-year OS and RFS in HCC were identified using K-M plotter through the third-party database (TCGA cohort database), respectively. Subsequent analysis demonstrated that CD5L and SRD5A2 were independent risk factors for 5-year OS and RFS, respectively. Importantly, all independent prognostic factors were combined to build a predictive nomogram model, which was effective in predicting the prognosis of HCC patients.

CD5L is mostly expressed in the macrophages, lymphoid and inflamed tissues and regulates inflammatory responses and lipid synthesis.^11,23,26 It is also known as the inhibitor of apoptosis in macrophages. It promotes macrophage survival by protecting them from the apoptotic effects of oxidized lipids in atherosclerosis.¹⁰ Moreover, CD5L is involved in the early response to the infection by bacteria and other pathogens, where it acts as a pattern recognition receptor and activates autophagy.^10,24 CD5L also controls the metabolic switch in T-helper Th17 cells and regulates their expression of pro-inflammatory genes.^9,23 CD5L accumulation on the hepatic surface could inhibit chronic liver injury by attenuating CCL4-induced injury and fibrosis, repressing TGF-ß signal and immune cell infiltration.⁶ Circulating CD5L potentially protects from the development of fatty liver and HCC.²³ The current assay also identified CD5L was associated with improved OS in HCC. Contrastingly, Aran et al observed that CD5L was upregulated in HCC and it enhanced HCC cell growth and antiapoptotic responses by binding to HSPA5 (GRP78).¹¹ It was also argued that CD5L was more suitable for HCC or cirrhosis accompanied by a viral infection than in the absence of an inflammatory response because CD5L itself was a factor associated with immune regulation.^9,23

Previous studies pointed to the CD5L-related modulation of immune responses in malignancies, such as lung adenocarcinoma and HCC. Interestingly, CD5L has opposite effects on these tumors, promoting lung cancer but inhibiting liver cancer. It was established that CD5L accumulates on the surface of transformed hepatocytes and induces necrosis of the tumor. CD5L-deficient mice were susceptible to HCC and formed multiple liver tumors after feeding with a high-fat diet for one year. It was found that mouse CD5L was internalized together with CD36 by normal hepatocytes and modulated intracellular lipid metabolism.²⁴ Therefore, CD5L may serve as a potential target for the treatment of HCC. In humans, CD5L protein is present at a high concentration in the serum, especially in women. However, CD5L peaks in women in their 20s and decreases with age. Several proteomic assays pointed to the CD5L protein as a putative biomarker for inflammatory conditions as well as liver diseases.²⁵ The continuous inflammation caused by liver damage due to hepatitis virus infection, alcohol abuse, and NASH leads to hepatic fibrosis, which frequently triggers cirrhosis and, ultimately leads to HCC. Based on these findings, hCD5L can also be considered as a plasma biomarker for early detection of liver fibrosis and HCC, and the ratio of hCD5L-to- liver marker score might discriminate between HCC and non-HCC patients. HCC. Thus, CD5L may serve as a potential target in the development of HCC treatment. Future studies are required to verify the role of CD5L in HCC and develop strategies to alter its expression.

SRD5A2 is expressed in androgen-dependent tissues and responsible for converting testosterone to the more metabolically active dihydrotestosterone. Therefore, SRD5A2 V89L gene polymorphism has been associated with breast cancer and prostate cancer in previous studies.^26–28 In the present study, we found that SRD5A2 was an independent risk factor for RFS in HCC. Liver regulates the metabolism and activity of sex hormones. Hence, SRD5A2 might be aberrantly expressed in HCC and serve as an useful biomarker for early diagnosis of HCC.^29–31 However, in the current study, the C-index of SRD5A2 was 0.67 and 0.68 indicating that SRD5A2 may not have an adequate performance for the prognostic predictive nomogram. Future studies are required to delineate the role of SRD5A2 in predicting the recurrence and survival of HCC patients.

Conclusion

In this study, 47 genes associated with HCC were identified, with most of them being located on chromosomes 5q and 6q. The potential pathways involving these genes were steroid metabolism, lipid metabolism, retinol metabolism, bile secretion, and p53 signalling pathway. There were 26 and 40 genes associated with the 5-year OS and RFS of HCC patients, respectively. Among them, CD5L and SRD5A2 were independent risk factors for 5-year OS and RFS. A nomogram model combining CD5L, cirrhosis, RSPM and BCLC stage was constructed for accurate prognostication of patients with HCC. CD5L might be useful as a potential biomarker for HCC.