Systems analysis of key genes and pathways in the progression of hepatocellular carcinoma.

The carcinogenesis of hepatocellular carcinoma (HCC) is a complex process, starting from a chronically altered hepatic microenvironment due to liver cirrhosis and ultimately progressing to HCC. However, the sequential molecular alterations driving the malignant transformation in liver cirrhosis are not clearly defined.In this study, we obtained gene expression profiles of HCC, including 268 tumor tissues, 243 adjacent tumor tissues, and 40 cirrhotic tissues (GSE25097) from Gene Expression Omnibus (GEO), to comprehensively define changes in the transcriptome of HCC during the sequential evolution of liver cirrhosis into HCC.We showed that changes in the molecular profiles of cirrhotic and adjacent tumor samples were small and quite uniform, whereas there was a striking increase in the heterogeneity of tumors in HCC tissues at the mRNA level. A massive deregulation of key oncogenic molecules and pathways was observed from cirrhosis to HCC tumors. In addition, we focused on FOXO1 and DCN, 2 critical tumor suppressor genes that play an important role in liver cirrhosis and HCC development. FOXO1 and DCN expression levels were significantly reduced in tumor tissues compared with adjacent tumor tissues in HCC. Kaplan-Meier analysis revealed that FOXO1 and DCN expression was positively correlated with overall survival, defining FOXO1 and DCN as adverse prognostic biomarkers for HCC.This system-level research provided new insights into the molecular mechanisms of HCC carcinogenesis. FOXO1 and DCN may be applied as potential targets for HCC treatment in the future.

Introduction

Hepatocellular carcinoma (HCC) is the fifth most common malignancy and the third leading cause of cancer-related death worldwide.[ HCC accounts for 80% to 90% of primary liver cancers, and the incidence of HCC is growing globally by 3% to 9% annually.[ HCC neoplasms detected at an early stage can be cured by mainly surgical resection. Treatment options for HCC at an advanced stage are often limited.[ The survival duration of patients with advanced liver cancer is less than 12 months.[ Early detection of HCC may help improve long-term survival rates.[ Therefore, there is an urgent need for a deeper understanding of the molecular mechanisms underlying the initiation and progression of HCC, and this information might be helpful for designing novel therapeutic strategies in the future.

Because the liver is especially susceptible to chronic and acute viral injury, alcoholic insults, and nonalcoholic fatty liver disease, it is extremely prone to fibrotic remodeling.[ Liver fibrosis usually progresses to cirrhosis, which can result in damage to the normal architecture of the liver, followed by an increased probability of the development of HCC.[ HCC occurs at a rate of 1% to 4% per year once liver cirrhosis is established, and liver cirrhosis underlies HCC in approximately 80% to 90% of cases worldwide.[ Increasing evidence has demonstrated that the carcinogenesis of HCC is a multistep process triggered by the accumulation of genetic alterations through the activation of different signaling pathways, which drives the transformation of normal cells into malignant cells.[ However, the mechanism behind the progression from liver cirrhosis to HCC remains largely unknown. To the best of our knowledge, no systematic study has been performed to investigate the molecular events leading from liver cirrhosis to HCC. A definition of the sequential molecular events leading from cirrhosis to HCC is urgently needed, and it represents a major challenge in the clinical management of at-risk patients.

In this study, we obtained a genome-wide expression profile of HCC, including 268 tumor tissues, 243 adjacent nontumor tissues, and 40 cirrhotic tissues (GSE25097),[ to comprehensively define changes in the transcriptome of HCC during the sequential evolution of cirrhosis into HCC, and we validated some of these results with 3 other HCC datasets (GSE22508,[ Oncopression database,[ and TCGA_LICH), which only included adjacent nontumor and tumor tissues. We showed that changes in the molecular profiles of cirrhotic and adjacent nontumor tissues were small and quite uniform in contrast to the striking increase in heterogeneity of HCC tissues at the mRNA level. A massive deregulation of key oncogenic molecules and pathways was observed from cirrhosis to HCC. In addition, we focused on FOXO1 and DCN, 2 critical tumor suppressor genes that play an important role in liver cirrhosis and HCC development. We detected the expression of FOXO1 and DCN in HCC and analyzed their correlation with clinical pathological features. Our data indicated that low FOXO1 or DCN expression was associated with poor prognosis of HCC.

Materials and methods

HCC datasets

The discovery dataset GSE25097 was obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/)[. The validation datasets were extracted from the following 3 datasets: GSE22058, which includes 96 paired adjacent nontumor and tumor samples of HCC from the GEO database[; Oncopression datasets (http://www.oncopression.com),[ including 524 tumor samples and 322 adjacent nontumor samples of HCC that integrate several gene expression datasets based on microarrays from different platforms into 1 large dataset; and the TCGA _LIHC dataset (http://tcga-data.nci. nih.gov, as of January 28, 2016), including 371 tumor samples and 50 adjacent nontumor samples of HCC with both mRNA expression data based on RNA-Seq and clinical feature information, which was used to perform the correlation analysis and survival analysis. All of the data in this study were based on previous published studies, and thus, no ethical approval and patient consent are required.

Functional enrichment analysis

Pathway analysis of different patterns of gene expression was performed using the Sigora R package version 2.0.1, which identified pathway enrichment based on statistically over-represented Pathway Gene-Pair Signatures.[ Signalling Pathway Impact Analysis (SPIA) was used to assess the importance of enriched pathways in terms of their impact and ability to activate or inhibit a pathway.[SPIA analysis was accomplished using the R Bioconductor package SPIA (version 2.18.0). Entrez IDs, log2-fold changes, and Q-values of all genes were compiled. SPIA produces a P value, which represents the significance level at which a pathway is found to be perturbed, and a false discovery rate (FDR). We ran SPIA using the recommended value of 2000 bootstrap iterations, and all parameters were set to their default values. A pathway was significant if the FDR was less than 0.1.

Statistical analysis

A gene was considered differentially expressed when it was significant at 5% FDR (q-value method) and showed an absolute log2 mean difference higher than 1 (double expression). Single comparisons between the 2 groups were determined by a Student t test. Survival analysis was performed with the Kaplan–Meier method, and the log-rank test was used to evaluate the statistical significance of the differences. Differences were considered to be statistically significant when P < .05.

Results

Group comparison of different stages of hepatocarcinogenesis

We computed the Pearson correlation coefficient of each sample at the corresponding stage based on the mRNA expression values from GSE25097, GSE22058, Oncopression database, and TCGA_LIHC datasets. The gene expression profiles of adjacent nontumor samples were quite homogeneous (the mean coefficient of the Pearson correlation was 0.94 in GSE25097, 0.97 in GSE22058, 0.89 in Oncopression, and 0.89 in TCGA_LIHC). As anticipated, high homogeneity was also observed for cirrhotic samples (the mean coefficient of the Pearson correlation was 0.93 in GSE25097). However, the homogeneity dramatically decreased upon progression to HCC (the mean coefficient of the Pearson correlation was 0.84 in GSE25097, 0.87 in GSE22058, 0.80 in Oncopression, and 0.66 in TCGA_LIHC, P < .0001), reflecting the well-recognized prototypical heterogeneity of HCC (Fig. 1).

Figure 1

The Pearson correlation coefficients of each sample at the corresponding stage based on the mRNA expression values from GSE25097 (A), GSE22058 (B), Oncopression database (C), and TCGA_LIHC datasets (D). ∗P < .0001.

Identification of differentially expressed genes (DEGs) among cirrhotic, adjacent nontumor, and tumor tissues of HCC

To investigate the gene expression alterations associated with HCC progression, GSE25097 was used as the discovery dataset for the identification of DEGs among cirrhotic, adjacent nontumor, and HCC tumor samples. This discovery dataset included the gene expression profiles of 268 tumor, 243 adjacent nontumor, and 40 cirrhotic samples. A total of 1920 genes (961 upregulated genes and 959 downregulated genes) were differentially expressed with a FDR < 5% and a log2 mean difference > 1 between cirrhotic and adjacent nontumor samples, and 2007 genes (966 upregulated genes and 1041 downregulated genes) were differentially expressed between adjacent nontumor and tumor samples (Fig. 2A, B). As shown in Fig. 2C, 1047 genes were significantly differentially expressed among cirrhotic, adjacent nontumor, and tumor samples. The DEGs between cirrhotic and adjacent non-tumor samples were classified as “Tumor-like,” “Trend,” or “Adjacent-specific” patterns based on their level of expression in tumor tissues, as defined by Sanz-Pamplona et al.[ There were 873 Tumor-like genes, 275 Trend genes, and 772 Adjacent-specific genes (Fig. 2D). Pathway analysis of different patterns of gene expression was conducted using Sigora, which limits the repetitive assignment of the same genes to multiple overlapping pathways. The enriched KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways for these gene patterns were Tumor-like pathways of Lysosome (hsa04142), Ribosome (hsa03010), Oxidative phosphorylation (hsa00190), ECM-receptor interaction (hsa04512), and NOD-like receptor signaling pathway (hsa04621). Trend pathways were enriched in Complement and coagulation cascades (hsa04610). Adjacent-specific pathways were most enriched in metabolism related pathways, such as glycine, serine, and threonine metabolism (hsa00260). These results suggest different functions for each gene expression pattern (see complete list in supplemental Tables S1, S2, and S3).

Figure 2

Differentially expressed genes (DEGs) among cirrhotic, adjacent nontumor, and tumor tissues of HCC. (A, B) Volcano plot of differential mRNA expression analysis. X-axis: log2-fold change; Y-axis: - log10 (FDR P value) of each probe; Vertical dotted lines: fold change ≥ 2 or ≤ 2; Horizontal dotted line: the significance cut-off (FDR = 5%). (A) There were 1920 genes identified as differentially expressed between cirrhotic and adjacent nontumor in GSE25097, including 961 upregulated genes and 959 downregulated genes. (B) Two thousand seven genes (1041 upregulated genes and 961 downregulated genes) were differentially expressed between adjacent nontumor and tumor in GSE25097. (C) Venn diagram showing the overlap of DEGs between cirrhotic and adjacent nontumor and tumor tissues of HCC. DEGs_CA were differentially expressed genes between cirrhotic and adjacent nontumor, and DEGs_AT were differentially expressed genes between adjacent nontumor and tumor of HCC. (D) Representative DEG patterns are displayed. DEG between adjacent and cirrhotic samples were classified as tumor-like, trend, and adjacent-specific genes.

Pathway enrichment analysis of DEGs among cirrhotic, adjacent nontumor, and tumor samples of HCC

Many of the existing pathway analysis methods are focused on either the number of DEGs in a pathway or on the correlation of genes in the pathway. Thus, information about complex gene interactions is disregarded. However, SPIA considers whether the DEGs found in a pathway have a meaningful impact within that pathway; thus, it addresses the topology of DEGs in pathways. Thus, in this study, we used SPIA to analyze the differences between aberrant pathways among cirrhotic, adjacent nontumor, and tumor tissues of HCC using the DEGs described above. A total of 59 KEGG pathways were identified as significantly perturbed in the progression from cirrhotic to adjacent nontumor (Table 1), and 40 KEGG pathways were significantly changed in the progression from adjacent nontumor to tumor (Table 2). Interestingly, most of the significantly perturbed pathways (50/59) during the transition from cirrhotic to adjacent nontumor were inhibited, whereas during the transition from adjacent nontumor to tumor, most of the pathways (29/39) were activated (P < .0001). Complement and coagulation cascades (hsa04610) and Antigen processing and presentation (hsa4612) were inhibited in the transitions from both cirrhotic to adjacent nontumor and adjacent nontumor to tumor; PPAR signaling pathway (hsa3320) was activated from both cirrhotic to adjacent nontumor and adjacent nontumor to tumor. ECM-receptor interaction (hsa4512), pathways in cancer (hsa5200), and insulin signaling pathway (hsa4910) were inhibited from cirrhotic to adjacent nontumor, but they were activated from adjacent nontumor to tumor. Focal adhesion (hsa4510) was activated from cirrhotic to adjacent nontumor, but it was inhibited from adjacent nontumor to tumor. The analysis also revealed that pathways involved in immune response were deregulated in the progression from cirrhotic to adjacent nontumor, including chemokine signaling pathway (hsa4062), natural killer cell mediated cytotoxicity (hsa4650), cytokine-cytokine receptor interaction (hsa4060), and Toll-like receptor signaling pathway (hsa4620). The tumor evolvement process was accompanied by an increase in the number of key oncogenic pathways associated with malignancy and metastatic spread, including p53 signaling pathway (hsa4155, activated), cell cycle (hsa4110, activated), transforming growth factor (TGF)-beta signaling pathway (hsa4350, activated), ErbB signaling pathway (hsa4012, activated), and NF-kappa B signaling pathway (hsa4064, inhibited). This functional analysis also suggested an active reaction of the adjacent nontumor related to the presence of the tumor or a more passive reaction induced by factors released from the tumor.

Table 1

Significantly impacted pathways between cirrhotic and adjacent nontumor as determined by SPIA.

Table 2

Significantly impacted pathways between adjacent nontumor and tumor as determined by SPIA.

FOXO1 and DCN (decorin) were underexpressed in HCC

Given that FOXO1 was one of the top-ranked dysregulation transcription factors (TFs), playing a crucial role in the dynamic regulation of gene expression programs in tumors, and that DCN was the most downregulated gene in HCC progression from cirrhosis to adjacent nontumor and tumor, we focused on these 2 genes for further analysis in this study. FOXO1 belonged to the adjacent-specific pattern. This classification highlighted its idiosyncratic role in the progression of HCC (Fig. 3A). There was a stepwise decrease of DCN mRNA expression levels in hepatocarcinogenesis from cirrhosis to adjacent nontumor and tumor tissues of HCC in microarray dataset GSE25097 (Fig. 3B). Furthermore, we validated their underexpression between adjacent nontumor and tumor using 3 other datasets (GSE22058, Oncopression, and TCGA_LIHC) from different platforms. Consistent with these findings, lower expression levels of FOXO1 and DCN were seen in malignant samples of the validation cohort (Fig. 4A–F). All these results suggested that the underexpression of FOXO1 and DCN was a common feature in HCC and that their dysregulation may be associated with tumorigenesis in HCC.

Figure 3

Gene expression levels of FOXO1 (A) and DCN (B) in cirrhotic (green), adjacent nontumor (pink), and tumor tissue (red) of HCC from GSE25097. ∗P < .0001.

Figure 4

Validation of differential expression of FOXO1 and DCN in validation datasets. (A–C) Reduced expression of FOXO1 in HCC. (D–F). Reduced expression of DCN in HCC. ∗P < .0001, †P = .001.

FOXO1 and DCN expression were prognostic indicators for patients with HCC

We also wanted to determine whether downregulated FOXO1 and DCN levels could be used as prognostic indicators for HCC. Gene expression and clinical information from The Cancer Genome Atlas (TCGA) was collected for further investigation. A total of 371 HCC patients were included, and their clinical characteristics are summarized in Table 3. FOXO1 and DCN expression was remarkably negatively associated with pathologic T stage and tumor grade (P < .05). Both the low FOXO1 expression group and the low DCN expression group had significantly poorer overall survival [FOXO1: P = .0072, hazard ratio: 0.6007, 95% confidence interval (95% CI): 0.4328–0.8514, Fig. 5A; DCN: P = .0326, hazard ratio: 0.8247, 95% CI: 0.5842–1.164, Fig. 5B]. The median survival period was 2116 days for the FOXO1 high expression group, whereas it dropped to 1271 days in the FOXO1 low expression group. The DCN low expression group had a reduced median survival period of 1397 days compared with a median survival period of 1694 days in the high expression group. These results indicated that FOXO1 and DCN were beneficial factors for survival in HCC patients.

Table 3

Clinical characteristics and correlations with the expression of FOXO1 and DCN.

Figure 5

Overall survival curves based on FOXO1 and DCN expression. (A) The survival curves of the FOXO1 high expression group and the FOXO1 low expression group. FOXO1 low expression was associated with poor overall survival (P = .0072). (B) The survival curves of the DCN high expression group and the DCN low expression group. DCN low expression was associated with poor overall survival (P = .0326). Patients were divided into low and high expression groups according to the gene's median mRNA expression value.

Discussion

The risk of HCC development is significantly increased among patients with advanced liver fibrosis caused by viral and nonviral etiologies, which then progressively evolves to cirrhosis. Although it is well-known that the prognosis of HCC patients is closely linked to levels of liver cirrhosis, the molecular mechanisms underlying the progression of liver cirrhosis to HCC remain unclear.

In this study, we found that the gene expression of cirrhotic samples and adjacent nontumor samples were surprisingly homogeneous; however, upon progression to HCC, the correlation between patients dramatically decreased, reflecting the well-recognized phenotypic heterogeneity of HCC. The high molecular heterogeneity of HCC may underscore the poor response to standard therapies in current clinical trials and the need for individualized treatment at progressed stages of HCC.

Previous studies have typically compared paired tumor and adjacent nontumor tissues, which can result in misleading interpretations. In this study, the inclusion of samples from cirrhotic tissues has allowed us to assess whether adjacent nontumor tissue from HCC patients differs from cirrhotic tissues due to tumor presence.

Changes in gene expression and cellular signaling pathways cause significant changes in the transition from liver cirrhosis to HCC. A total of 1920 and 2007 genes were specifically associated with progression from cirrhotic to adjacent nontumor and from adjacent nontumor to tumor, respectively. DEGs between cirrhotic and adjacent nontumor tissues can be grouped into 3 altered patterns based on their level of expression in tumor tissues: “Tumor-like,” “Trend,” and “Adjacent-specific.” Our results showed that each gene expression pattern has different functions. Moreover, SPIA was used for functional enrichment analysis of DEGs among cirrhotic, adjacent nontumor, and tumor samples of HCC. Most of the significantly perturbed pathways between cirrhotic to adjacent nontumor were inhibited, but between adjacent nontumor to tumor, most of these pathways were activated. Pathway analysis of cirrhotic to adjacent nontumor also revealed the deregulation of signaling pathways involved the immune response. Consistently, the majority of pathways activated during the conversion from adjacent nontumor to tumor conferred malignant and invasive properties, including p53 signaling pathway, cell cycle, TGF-beta signaling pathway, ErbB signaling pathway, and NF-kappa B signaling pathway. This analysis emphasized that the acquisition of malignant traits is a relatively late event. Our results also suggested that adjacent tumor tissue is abnormal. In fact, studies that only compare tumor and adjacent nontumor may miss good cancer biomarker candidates because many genes are deregulated in adjacent tumor tissue, mimicking tumor expression.

Another important finding of this study is identifying FOXO1 and DCN preferentially upregulation in adjacent nontumor tissue and downregulation in tumor tissue of HCC. FOXO1 belonged to the adjacent-specific gene pattern, and DCN belonged to the Trend gene pattern. The TF FOXO1 is characterized by the forkhead DNA-binding domain, and its aberration influences multiple cellular functions, including apoptosis, cell cycle control, DNA damage repair, glucose metabolism, carcinogenesis, and tumor immunity.[ In addition, as a critical modulator of many important stress pathways, the FOXO1 TF, may regulate adaptation of the liver to stress.[ Although FOXO1 dysregulation has been observed in several human cancers,[ the study of its expression in liver cirrhosis has been limited. It is well-known that hepatic stellate cells (HSCs) play a crucial role in the liver fibrotic response, as their activation, transdifferentiation, and proliferation are key steps in liver fibrosis. Intriguingly, FOXO1 was reported to participate in the proliferation and transdifferentiation process of HSCs, which were enhanced by transcriptionally inactive FOXO1. In contrast, active FOXO1 inhibited HSC proliferation by inducing cell cycle arrest to accumulate cells in the G0/G1 phase.[ It is consistent with our results that FOXO1 gene expression was downregulated in patients with cirrhosis. Our finding shows that the functional role of FOXO1 may be linked to hepatocarcinogenesis at the stage in which cirrhosis progresses into tumor in HCC.

Another potential prognostic factor for HCC is DCN. The DCN gene encodes a member of the small leucine-rich proteoglycan family of proteins regulating collagen fibrillogenesis during liver disease development. Previous studies showed that a low quantity of decorin is present in normal healthy liver, whereas the amount of decorin significantly increases during fibrogenesis.[ Consistent with these findings, our results indicated that higher expression levels of DCN were seen in the progression of liver cirrhosis. Furthermore, the evidence suggests that decorin could act as a tumor repressor in a variety of cancers. An early study has shown that extopic expression of decorin could suppress the generalized growth of various neoplastic cells.[ Regarding liver tumors, decorin inhibits the proliferation of HepG2 and Huh-7 hepatoma cell lines.[ In addition, it has been shown that although decorin-deficient mice are not associated with the development of spontaneous tumors, their tissues are permissive for tumorigenesis.[ In line with these reports, this study indicated that DCN mRNA expression is downregulated in HCC, as its lack is accompanied with significantly higher HCC prevalence.

Conclusion

We have provided the first evidence of the molecular mechanism in the transition from liver cirrhosis to HCC using bioinformatics techniques. In addition, we also found that FOXO1 and DCN are underexpressed in HCC tissues and that their downregulation may be indicative of poor survival rates; furthermore, they could have a potential role as prognostic markers in HCC patients. Functional studies are needed to reveal the molecular mechanisms of FOXO1 and DCN in HCC and their role in prognosis and therapeutic targets.

Author contributions

Data curation: Yu-Kui Shang, Fanni Li.

Formal analysis: Zi-Ling Wang.

Funding acquisition: Fanni Li, Huijie Bian, Zhi-Nan Chen.

Investigation: Zi-Ling Wang, Huijie Bian.

Methodology: Yu-Kui Shang, Fanni Li, Zekun Liu, Zi-Ling Wang.

Project administration: Zhi-Nan Chen.

Resources: Yi Zhang, Zekun Liu.

Software: Yi Zhang.

Validation: Yi Zhang, Huijie Bian, Zhi-Nan Chen.

Visualization: Yi Zhang, Huijie Bian.

Writing – original draft: Yu-Kui Shang, Fanni Li, Huijie Bian.

Writing – review & editing: Huijie Bian, Zhi-Nan Chen.