Seven-senescence-associated gene signature predicts overall survival for Asian patients with hepatocellular carcinoma

Xiao-Hong Xiang, Li Yang, Xing Zhang, Xiao-Hua Ma, Run-Chen Miao, Jing-Xian Gu, Yu-Nong Fu, Qing Yao,Jing-Yao Zhang, Chang Liu, Ting Lin, Kai Qu

Abstrac t BACKGROUND Cellular senescence is a recognized barrier for progression of chronic liver diseases to hepatocellular carcinoma (HCC). The expression of a cluster of genes is altered in response to environmental factors during senescence. However, it is questionable whether these genes could serve as biomarkers for HCC patients.AIM To develop a signature of senescence-associated genes (SAGs) that predicts patients' overall survival (OS) to improve prognosis prediction of HCC.METHODS SAGs were identified using two senescent cell models. Univariate COX regression analysis was performed to screen the candidate genes significantly associated with OS of HCC in a discovery cohort (GSE14520) for the least absolute shrinkage and selection operator modelling. Prognostic value of this seven-gene signature was evaluated using two independent cohorts retrieved from the GEO (GSE14520) and the Cancer Genome Atlas datasets, respectively.Time-dependent receiver operating characteristic (ROC) curve analysis was conducted to compare the predictive accuracy of the seven-SAG signature and serum α-fetoprotein (AFP).RESULTS A total of 42 SAGs were screened and seven of them, including KIF18B, CEP55,CIT, MCM7, CDC45, EZH2, and MCM5, were used to construct a prognostic formula. All seven genes were significantly downregulated in senescent cells and None.Song H upregulated in HCC tissues. Survival analysis indicated that our seven-SAG signature was strongly associated with OS, especially in Asian populations, both in discovery and validation cohorts. Moreover, time-dependent ROC curve analysis suggested the seven-gene signature had a better predictive accuracy than serum AFP in predicting HCC patients' 1-, 3-, and 5-year OS.CONCLUSION We developed a seven-SAG signature, which could predict OS of Asian HCC patients. This risk model provides new clinical evidence for the accurate diagnosis and targeted treatment of HCC.

Key words: Senescence-associated genes; Hepatocellular carcinoma; Overall survival;Risk model; Asian patients

INTRODUCTION

Hepatocellular carcinoma (HCC) becomes the third leading cause of cancer d eaths worldwide, with approximately 80% of mortalities occurring within 5 years[1,2]. During the past decades, great effects have been mad e to improve the management of HCC.As more and more HCCs are diagnosed at an early stage, treatment efficacy is greatly improved[3]. Whereas, deaths caused by HCC alw ays occur w hen patients undergo treatment and the clinical outcome of HCC patients is still poor[4]. Moreover, HCC is an extremely heterogeneous disease, w hich must be monitored for high-risk patients w ith p oor clinical outcomes and ad op t effective treatments to imp rove p atient survival[5]. Trad itional serum markers, especially alpha-fetoprotein (AFP), have been the most common prognostic ind icators in clinic. However, they significantly depend on tumor burd en, w hich limits their value in d iagnosing early stage tumors.Therefore, identifying novel prognostic biomarkers contributes to early diagnosis and reducing HCC mortality.

Cellular senescence is considered to be a response of a proliferating somatic cell to stress and d amage d erived from both exogenous and end ogenous sources, and p ersistent DNA d amage is the most common cause[6]. It is characterized as a permanent cell cycle arrest[7]. Cellular senescence has been deemed as a mechanism of limited cell division due to progressive telomere shortening[8]. In cancer cells, there exists telomere-ind ep end ent senescence due to their activation of telomerase. For instance, numerous stud ies found that RAS activation could ind uce cellular senescence in many cancer cell types. Both telomere-d epend ent and -independent pathw ays induced senescence by inducing a DNA damage response. Recently, it has been demonstrated that hepatocytes escaping from senescence played a key role in hepatic carcinogenesis and HCC p rogression[9,10]. During the senescence process, a large number of genes are expressed d ifferentially. Hence, it is reliable to focus on senescence associated genes (SAGs) as a novel ap proach in cancer d iagnosis and monitor.

Bearing this in mind, w e first id entified the SAGs by analyzing the genome profiling d ata derived from two types of senescent cells, replicative senescence (RS)and oncogene-induced senescence (OIS). Next, w e investigated the prognostic value of several cand idate SAGs in predicting survival of HCC by multistep comparisons.Finally, w e valid ated that the cand id ate SAGs w ere sup erior to classic serum biomarker AFP in p red icting the OS in HCC cohort retrieved from the Cancer Genome Atlas (TCGA). As far as w e know, the stud y is the first to investigate the expression patterns of SAGs between senescence and HCC and their association with clinical prognosis of Asian patients w ith HCC.

MATERIALS AND METHODS

Data sources and processing

Tw o genomic profiling datasets of RS cells (GSE19018 and GSE36640), three datasets of OIS cells (GSE19864, GSE40349, and GSE60652), and one HCC dataset (GSE14520)w er e obtain ed f r om th e Gene Exp ression Om n ibu s (GEO) d atabase(http://ww w.ncbi.nlm.nih.gov/geo/). RMA algorithm w as performed to normalize and transform all the selected d ata from GEO to exp ression values in the R environment (v3.5.1). Among them, GSE14520 (based on GPL3921 platform) enrolled 209 HCC patients w ith survival d ata. An ind ep end ent HCC cohort d erived from TCGA database (TCGA-LIHC) was used as a validation group. The level 3 RNAseqv2 data and clinical d ata were derived from the TCGA data portal. A total of 370 HCC samples were included.

Differential gene expression analysis

Differentlly expressed genes (DEGs) were calculated in RS and OIS cell models,respectively. Only fold change (FC) ≥ 1.5 and q-value ＜ 0.05 were consid ered statistically significant. Then, we picked up the overlapped significantly expressed genes in the two senescent models. Venn diagram was carried out using Venny 2.1.0.Moreover, the expression levels of candidate genes were plotted and analyzed by the t-test in GSE14520 and TCGA-LIHC cohorts.

Prognostic model development

A prognostic model was created by the seven-SAG signature according to least absolute shrinkage and selection operator (LASSO) analysis. LASSO is one of the most popular approaches for sparse linear regression[11]. “Almnet” package was carried out based on a series of λ in the R environment (v3.5.1)[12]and the coefficients of each gene in the risk score system were generated. We got a risk score for every patient based on their own expression levels of the seven genes after the LASSO regression analysis.

Survival analysis

Univariate and multivariate survival analyses w ere carried out and further multivariate COX regression analysis only included variables with P ＜ 0.05. All tests were carried out using SPSS (version 24.0; Chicago, United States). Kaplan-Meier curves were generated using GraphPad Prism 7.0. Comparisons between different subgroups were performed by the Log-Rank test. Patients are divided into high- and low- risk groups by the median.

Time-dependent receiver operating characteristic curve (ROC) analysis

ROC curve is extended to evaluate biomarker's accuracy of d iscriminating binary outcomes[13]. Ind ivid uals w ith a high risk of d evelop ing the d isease later may be disease-free in earlier life and their markers' value may change from baseline during follow-up. Therefore, time-d ep end ent ROC curve analysis is more appropriate and outperforms the conventional method adopted for handling censored biomarker data.In this study, the time-dependent ROC curve analysis w as performed w ith “survival ROC” package (R version 3.5.1). The prognostic performance was evaluated at 1, 3,and 5 years to compare the predictive accuracy and sensitivity of different prognostic models.

Statistical analysis

The chi-square test w as carried out to d iscover the relationship betw een gene expression and clinical parameters. Unpaired student's t-test w as used to analyze the d ifference of gene exp ression in HCC patients of d ifferent features. P ＜ 0.05 w as considered statistically different. Statistical analyses w ere performed using IBM SPSS Statistics softw are program version 24.0 (IBM Corp, NY, United States).

RESULTS

Identification of SAGs using different senescent models

The overall workflow of the data processing is presented in Figure 1A. To identify SAGs, we first integrated five different microarray profiles (GSE19018, GSE36640,GSE19864, GSE40349, and GSE60652). All datasets used in the present study were normalized before analysis. The relative expression of all samples pre- and postnormalization is shown in Figure 1B. Next, w e screened the DEGs, which w ere identified as FC ≥ 1.5 and q-value ＜ 0.05 (senescent vs proliferating cells), in RS and OIS models, respectively. A total of 781 up-regulated and 739 down-regulated genes were selected in the RS model, and 103 up-regulated and 288 down-regulated genes in the OIS model (Figure 2A). By overlapping the tw o DEG lists, 42 common differentially expressed genes (35 downregulated and only 7 upregulated genes) were selected as SAGs (Figure 2B) and the expression levels of these genes in RS and OIS models are presented as a heat map in Figure 2C.

Construction of a risk score system

To obtain a prognostic SAG signature for HCC survival prediction, we first performed univariate COX regression analysis to evaluate the prognostic value of each candidate gene. And we found that all seven genes (CEP55, MCM7, CDC45, MCM5, KIF18B, CIT,and EZH2) were proved to be risk factors for HCC patients. Next, we screened the expression pattern of the above 7 candidate genes in HCC cohorts. Intriguingly, seven dow nregulated genes in senescent cells, were significantly upregulated in HCC tissues in both discovery and validation groups, with a P-value ＜ 0.0001 (Figure 3).We then developed a risk score formula based on these seven SAGs using the LASSO method: Risk score = (0.243 × relative expression value of KIF18B) + (0.274 × relative expression value of CEP55) + (0.282 × relative expression value of CIT) + (0.266 ×relative expression value of MCM7) + (0.678 × relative expression value of CDC45) +(0.175 × relative expression value of EZH2) + (0.536 × relative expression value of MCM5). In this risk score system, the contribution of every gene to the risk score model was weighted by absolute value of coefficients. Every patient would get a risk score based on the expression of the seven SAGs of themselves.

Validation of the seven-SAG signature for prognosis

To confirm the potentiality of the seven-SAG prognostic model, Kaplan-Meier curve was carried out to evaluate the association between the overall survival (OS) and our gene signature in discovery (GSE14520) and validation (TCGA-LIHC) cohorts. The whole group was divided into the high- and low-risk subgroups according to the median of all patients' risk scores. In the discovery cohort, with the increase in the risk score, the expression of all the seven genes was increasing, and the death events accumulated (Figure 4A). The patients in the high-risk subgroup had a 1.92-fold higher death risk than the low subgroup [hazard ratio (HR), 95% confidence interval(CI) = 1.92, 1.16-3.19; log-rank P value = 0.011] (Figure 4B). We then tempted to test these findings in the validation cohort (TCGA-LIHC) (Figure 4C). Similar to the findings obtained from the discovery cohort, patients in the high-risk group [median survival time (MST) = 46.6 m] had significantly shorter OS time than patients with a low-risk score (MST = 70.5 m) [HR (95%CI) = 1.80 (1.27-2.54), log-rank P value = 0.001](Figure 4D). Interestingly, when we analyzed the data in the Asian population, we observed a highly significant association between the seven-SAG signature and OS.The majority of death events occurred in the high-risk group (Figure 4E). Asian HCC patients with a high-risk score were shown to have a ＞ 5-fold increased death risk than low-risk patients [HR (95%CI) = 5.81 (3.20-10.54), log-rank P value ＜ 0.0001]. The MST of the high-risk subgroup was only 60% of that of the low-risk group (MST =21.6 m vs 91.7 m) (Figure 4F). In order to investigate the prognostic value of the risk score system in different patient groups with different characteristics, we performed univariate/multivariate Cox regression analysis of clinicop athologic factors associated with OS in the discovery and validation cohorts. From the Cox regression results, both the seven-SAG signature and serum AFP level were confirmed to be independent risk factors of OS in the two cohorts (Table 1).

Comparison of the seven-SAG signature and serum AFP in predicting OS

To assess the prognostic accuracy of the seven-SAG signature, w e performed timedepend ent ROC analysis of 1-, 3-, and 5-year OS of the valid ation cohort. The area under the curve (AUC) of the seven-SAG signature model indicated an accep table p red ictive accuracy, w hich is sup erior to AFP, a w id ely used trad itional serum marker, at 1 year (our seven-gene model AUC = 0.708, serum AFP level AUC = 0.606)(Figure 5A), 3 years (our seven-gene mod el AUC = 0.699, serum AFP level AUC =0.568) (Figure 5B) as w ell as 5 years (our seven-gene model AUC = 0.678, serum AFP level AUC = 0.604) (Figure 5C). These results ind icated the valid ation of the prognostic signature.

Figure 1 Data sources and processing. A: Flowchart describing the process used to generate differentially expressed senescence-associated genes; B: The relative expression in all samples pre- and post-normalization.

Stratified analysis of the seven-SAG signature for prognosis prediction

Stratified analyses based on the clinical characteristics w ere carried out to identify the suitable Asian patient group s for the seven-SAG signature (Table 2). In the elderly population, patients with a high-risk score had a more than 3-fold increased risk of death than the low-risk group (Figure 6A-C). These results suggested that our seven-SAG signature w as more applicable to the HCC patients with older age in pred icting OS.

DISCUSSION

In this study, we first identified 42 overlapped DEGs using the RS and OIS models.Among them, seven downregulated genes in senescent cells, KIF18B, CEP55, CIT,MCM7, CDC45, EZH2, and MCM5, were shown to be upregulated in HCC tissues and selected to construct a prognostic model. The seven-SAG signature was shown to be associated w ith OS in both discovery and validation cohorts. Stratified analysis showed that our seven-SAG signature was significantly associated with OS in elderly Asian patients. Moreover, time-d epend ent ROC analysis show ed a favorable prognostic value of our seven-SAG signature when compared with serum AFP.

The cellular senescence is considered an aging hallmark. With the increase in age,the number of senescent cells is increasing. Cellular senescence is widely considered to be an anti-tumor mechanism. Studies have show n that a source of stress that triggers liver senescence is chronic inflammation, which causes damage to liver cell regeneration. Importantly, abrogation of senescence lead s to aggressive HCC development[14]. In the present study, we found that the HCC patients carrying high expression of seven SAGs had a shorter OS time. Stratified results further suggested the seven-SAG signature w as more applicable to the elderly HCC patients. The potential explanation might be that due to the increasing number of senescent cells,the expression of the seven-SAG signature is decreased, while its high expression indicates a higher proliferation rate and poorer OS.

Figure 2 Expression of differentially expressed genes in senescent cells. A: The differentially expressed genes were identified in replicative senescence (RS) and oncogene-induced senescence (OIS) models. There were 781 upregulated and 739 down-regulated genes selected using the RS model, and 103 up-regulated and 288 downregulated genes selected using the OIS model; B: Overlapping the two lists, 42 common differentially expressed genes were selected as senescence associated genes; C: The expression levels of 42 genes are presented as a heat map in RS and OIS models. DEGs: Differentially expressed genes; RS: Replicative senescence; OIS: Oncogeneinduced senescence; SAGs: Senescence associated genes; FC: Fold change; LASSO: Least absolute shrinkage and selection operator.

It has been widely accepted that senescence pathways are collectively at the level of activation of CDKIs, which play pivotal roles in regulating the cell cycle progression.Of the seven genes, KIF18B, a member of the kinesin-8 subfamily, is involved in cell cycle process[15]and acts as an oncogene in cervical cancer[16]. CEP55, also known as c10orf3 and FLJ10540, promotes tumorigenesis and regulates stemness in various cancers, such as lung adenocarcinoma[17-19]. CIT (Serine/threonine kinase 21) encoding a serine/threonine protein kinase, is a downstream effector of Rho family GTPases and participates in cell cycle regulation. Liu et al[20]and Xu et al[21]have demonstrated that CIT is up-regulated in HCC and regulates the G2/M transition in rat hepatocytes.EZH2 is a subunit of polycomb repressive complex 2 (PRC2), a protein complex that induces epigenetically silencing of genes[22]. And it has been reported that several lncRNAs are able to regulate gene transcription by binding to PRC2[23,24].

In most eukaryotes, the MCM complex consists of six highly conserved MCM proteins, namely MCM2-7, which functions as a replicative DNA helicase to unwind the DNA d up lex template d uring DNA rep lication[25]. Recent evid ence has d emonstrated that several MCM p roteins are tightly associated w ith tumorigenesis[26-28]. The MCM2-7 hexamer complexes with CDC45 and the heterotetrameric GINS complex, the Cdc45-Mcm2-7-GINS (CMG) complex, function as a p otential target for cancer treatment and CDC45 interacts w ith MCM 2[29].Furthermore, our previous study showed that MCM7 promotes cancer progression through cyclin D1-dependent signaling in HCC[30]. Three members of CMG complex,MCM 5, MCM 7, and CDC45, w ere includ ed in the risk formula. This evid ence indicated that more attention should be focused on the pro-oncogenic mechanisms of senescence escape of HCC cells.

Figure 3 Expression of the seven candidate genes in hepatocellular carcinoma and normal liver tissues. A: MCM5; B: MCM7; C: EZH2; D: CDC45; E: CIT; F:KIF18B; G: CEP55. Seven senescence associated genes, which were downregulated in senescent cells, were shown to be upregulated in hepatocellular carcinoma tissues in both discovery and validation cohorts. P ＜ 0.0001 for all.

Our study showed that AFP was an independent risk factor for HCC patients. AFP is often expressed at high levels in most HCC patients and is considered a reliable clinical tumor biomarker. As a classic serum biomarker, AFP w as found to be associated w ith p rognosis in HCC p atients[31,32]. Park et al[33]rep orted that AFP combined w ith PIVKA-II w ere useful in p red icting survival in the rad iological treatment of locally advanced HCC. Jiang et al[34]also found that preoperative AFP and fibrinogen show ed a p red ictive p ow er for recurrence of HCC after liver transplantation. Here, our study demonstrated that the seven-SAG signature model w as sup erior to serum AFP level and more ap plicable in p red icting OS of HCC patients w ith older age. The above findings suggested a potential clinical applicationof the seven-SAG signature in HCC patients.

Table 1 Univariate/multivariate Cox regression analysis of clinicopathologic factors associated with overall survival in GSE14520 and The Cancer Genome Atlas cohorts

However, there are some limitations in our study. First, the samples for screening SAG w ere small, w hich might cause false positive results. Second, we constructed the risk score system merely based on the gene expression levels, rather than the other genetic events that probably have an effect on the initiation and progression of cancer.Third, patients in the discovery cohort w ere from Asia, thus, the risk score system was established based on an Asian background. And further stratified analysis in the valid ation cohort also showed that this model w as more suitable for Asian patients.Hence, our HCC prognostic signature still need s to be validated in a larger group of patients from various populations.

In conclusion, w e constructed and confirmed a p rognostic risk score system comprised of seven SAGs. The seven-SAG signature could be a potential predictor for OS, p articularly in eld erly Asian HCC p atients. Our d ata provide new promising evidence on prediction biomarkers and targeted therapy for HCC.

Table 2 Stratified analysis of overall survival in GSE14520 and The Cancer Genome Atlas cohorts

Figure 4 The seven-senescence associated gene signature is associated with overall survival in the discovery and validation cohorts. A: The heat map of the expression of seven genes and patients' death status in the discovery cohort; B: Kaplan-Meier survival curves plotted to estimate the overall survival probabilities for the low-risk vs high-risk group in the discovery cohort; C: The heat map of the expression of seven genes and patients' death status in the validation cohort; D:Kaplan-Meier survival curves plotted to estimate the overall survival probabilities for the low-risk vs high-risk group in the validation cohort; E: The heat map of the expression of seven genes and patients' death status in the validation cohort-Asian only subgroup; F: Kaplan-Meier survival curves plotted to estimate the overall survival probabilities for the low-risk vs high-risk group in the validation cohort-Asian only subgroup. HR: Hazard ratio; CI: Confidence interval.

Figure 5 Comparison of the seven-senescence associated gene signature and α-fetoprotein in predicting overall survival of hepatocellular carcinoma patients. To assess the prognostic accuracy of the seven-senescence associated gene signature and serumα-fetoprotein level, time-dependent receiver operating characteristic curve analysis was conducted for 1-, 3-, and 5-year overall survival (OS). A: 1-year OS; B: 3-year OS; C: 5-year OS. AUC: Area under the curve; AFP:Alpha-fetoprotein; OS: Overall survival.

Figure 6 Association between the seven-senescence associated gene signature and overall survival in the elderly age subgroup. Kaplan-Meier survival curves were plotted to estimate the overall survival probabilities for the low-risk vs high-risk group. A: Discovery group; B: Validation group; C: Combination of discovery group and validation group. HR: Hazard ratio; CI: Confidence interval.

ARTICLE HIGHLIGHTS

Research background

Hepatocellular carcinoma (HCC) is a common malignancy that remains a serious cause of death worldwide. Recently, molecular markers and prognostic models have been used to improve the diagnosis and treatment of HCC, but few can be applied clinically. Currently, bioinformatics technology has been used for data mining in large public databases. The abundant sample size in the public database can make up for the shortcomings of small samples in real hospitals and help to seek for a more accurate and applicable prognostic model for HCC.

Research motivation

Researchers have been making efforts to find molecular markers or prognostic models that can effectively predict the prognosis of HCC. Senescence is a cell cycle arrest caused by stress in cells,but the cells are still alive. Studies have shown that the proportion of senescent cells in tissues of patients with cirrhosis increases, but a considerable number of patients with cirrhosis can develop liver cancer, and its specific molecular mechanism has rarely been reported.

Research objectives

By analyzing the database of two cellular senescence models from Gene Expression Omnibus,we screened for senescence-associated genes and validated these genes in the liver cancer databases (GSE14520 and TCGA-LIHC). Then, we constructed an HCC prognostic model and evaluate its prognostic accuracy.

Research methods

Senescence-associated genes (SAGs) w ere identified using R package “limma”. The latest statistical algorithm-the least absolute shrinkage and selection operator (LASSO) w as applied to create our prognostic model. Time-dependent receiver operating characteristic (ROC) curves w ere used to compare the prognostic accuracy betw een the seven-SAG signature and serum αfetoprotein.

Research results

The prognostic model for predicting the overall survival (OS) of HCC w as constructed by LASSO, consisting of the seven senescence-associated genes (SAGs) (KIF18B, CEP55, CIT, MCM7,CDC45, EZH2, and MCM5). All seven SAGs w ere highly expressed in HCC and proliferating cells, while lowly expressed in normal tissues and senescent cells. Survival analysis showed that our seven-SAG characteristics are closely related to OS, especially in Asian populations, both in the d iscovery and validation cohorts. In addition, the time-dependent ROC curve analysis indicated that the seven-gene marker is better than serum alpha-fetoprotein in predicting 1-, 3-,and 5-year OS of HCC patients.

Research conclusions

The seven-SAG signature was more applicable to evaluate OS of Asian HCC patients, which may provid e new clinical evidence for the diagnosis and treatment of HCC transformed from cirrhosis.

Research perspectives

The current study provides clues that the expression changes of senescence-associated gene are the molecular basis for the progression of cirrhosis to liver cancer. Finding effective senescenceassociated molecular biomarkers and predictive features of HCC prognosis is necessary.