Microarray analysis for the identification of specific proteins and functional modules involved in the process of hepatocellular carcinoma originating from cirrhotic liver.

In order to identify the potential pathogenesis of hepatocellular carcinoma (HCC) developing from cirrhosis, a microarray‑based transcriptome profile was analyzed. The GSE63898 expression profile was downloaded from the Gene Expression Omnibus database, which included data from 228 HCC tissue samples and 168 cirrhotic tissue samples. The Robust Multi‑array Average in the Affy package of R was used for raw data processing and Student's t‑test was used to screen differentially expressed genes (DEGs). An enrichment analysis was then conducted using the Database for Annotation, Visualization and Integrated Discovery online tool, and the protein‑protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes and Cytoscape. Furthermore, the MCODE plug‑in of Cytoscape was used to conduct a sub‑module analysis. A total of 634 DEGs were identified between HCC and cirrhosis, of which 165 were upregulated and 469 were downregulated. According to the cut‑off criteria, the PPI network was constructed and Jun proto‑oncogene, AP‑1 transcription factor subunit (degree, 39), Fos proto‑oncogene, AP‑1 transcription factor subunit (degree, 34) and v‑myc avian myelocytomatosis viral oncogene homolog (degree, 32) were identified as the hub nodes of the PPI network. Based on the sub‑module analysis, four specific modules were identified. In particular, module 1 was significantly enriched in the chemokine signaling pathway, and C‑X‑C motif chemokine ligand 12, C‑C motif chemokine receptor 7 (CCR7) and C‑C motif chemokine ligand 5 (CCL5) were three important proteins in this module. Module 4 was significantly enriched in chemical carcinogenesis, and cytochrome P450 family 2 subfamily E member 1, cytochrome P450 family 2 subfamily C member 9 (CYP2C9) and cytochrome P450 family 2 subfamily A member 6 (CYP2A6) were three important proteins in this module. In conclusion, the present study revealed that CCR7, CCL5, CYP2C9 and CYP2A6 are novel genes identified in the development of HCC; however, the actual functions of these genes require verification.

Introduction

Hepatocellular carcinoma (HCC) is one of the most common types of malignant cancer, which is associated with a poor prognosis. In addition, HCC is the fourth most common cause of cancer-associated mortality worldwide (1). It was previously estimated that 39,230 new cases of HCC and 27,170 cases of HCC-associated mortality would occur in the United States in 2016 (1). Due to the highest occurrence of HCC in East Asia, this estimation should be further increased in this area (2). Chronic liver diseases, including hepatitis, fibrosis and cirrhosis, are the primary causes of HCC. In addition, ~90% of cases of HCC develop from cirrhotic livers (3), and HCC and cirrhosis share numerous risk factors, including hepatitis, alcohol consumption, obesity, diabetes, gender and advanced age (4,5). These findings suggest that cirrhosis and HCC may share pathophysiological similarities.

Despite the aforementioned correlations between HCC and cirrhosis, the detailed mechanisms linking these two diseases remain to be fully elucidated. Hepatic stellate cells (HSCs) are primary fibrogenic cells, which are involved in cirrhosis-dependent carcinogenesis (6). Su et al (7) indicated that HSCs can be activated by inflammatory signals and macrophages, which contribute to fibrogenesis, and therefore may be considered positive risk factors for HCC. Furthermore, it has been demonstrated that activated HSCs can foster a conducive environment to directly support hepatic tumorigenesis via secreting numerous cytokines, including Wnt ligands, interleukin 6 and growth factors (8). De Minicis et al (9) reported that HSCs can be activated by lipopolysaccharide produced by the intestinal microflora via the toll-like receptor 4 signaling pathway, which can upregulate the expression of proinflammatory chemokines or cytokines, and facilitate the invasion and migration of HCC (10). Increased extracellular matrix (ECM) production is another feature of cirrhosis that contributes to tumorigenesis (11). Previous studies have reported that increased ECM production may promote the growth, survival and migration of precancerous cells via integrin signaling (6,12). Furthermore, increased ECM production may disturb cell signaling by sequestering growth factors, including interleukin families, fibroblast growth factor and transforming growth factor (13). This sequestration may serve a role in the aberration of normal liver cells. However, although these existing studies have confirmed that the development of HCC is closely associated with cirrhosis, the detailed mechanisms remain to be elucidated. Therefore, further investigations into how HCC develops from cirrhosis are required.

In order to identify the epigenetic alterations and potential role of DNA methylation markers in HCC and its prognosis, Villanueva et al (14) generated a methylation-based prognostic signature based on 304 samples from patients with HCC after surgical resection using a training-validation scheme. The GSE63898 dataset contains data from this study, which analyzed whole-genome transcriptome alterations (14). Villanueva et al (14) confirmed a high prevalence of methylation deregulated genes, including Ras association domain family member 1, APC, WNT signaling pathway regulator and insulin like growth factor 2 (IGF2), and identified potential epidrivers, such as ephrin B2 and septin 9; however, the mechanisms underlying the development of HCC from cirrhosis were not determined. In order to investigate the pathogenesis of HCC, the GSE63898 dataset was used to analyze differences between HCC and cirrhotic liver tissues. The present study aimed to identify potential target molecules, which may aid HCC clinical treatment.

Materials and methods

Data acquisition

The GSE63898 expression profile (14), which contained data from HCC and liver cirrhosis samples, was downloaded from the Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo/). The samples were sequenced using the Affymetrix Human Genome U219 Array [HG-U219] platform. In this profile, a total of 396 liver tissue samples were collected, including 228 HCC tissue samples and 168 non-tumor liver adjacent cirrhotic tissue samples. The research was authorized by the Institutional Review Boards of the participating centers.

Data preprocessing

The downloaded raw data were preprocessed using the Robust Multi-array Average method (15,16) in the Affy package of R (17), including background correction, normalization and expression calculation. Gene expression was presented as the mean value of different probes.

Differentially expressed genes (DEGs) screening

DEGs between HCC samples and cirrhotic liver samples were screened by Student's t-test in the Linear Models for Microarray Data package (18). The P-value obtained from the Student's t-test was adjusted according to the Benjamini and Hochberg method (19) and the obtained adjusted P<0.01 was set as the threshold.

Functional and pathway enrichment analyses

Gene Ontology (GO) is a free database used in biological process (BP), molecular function and cellular component analyses (20). Kyoto Encyclopedia of Genes and Genomes (KEGG) is also a free database specifically used in pathway analysis (21). The Database for Annotation, Visualization and Integrated Discovery (DAVID) is a common online tool used in large-scale gene functional analyses (22). Therefore, GO functional and KEGG pathway enrichment analyses of DEGs were conducted using DAVID with the criterion of P<0.05.

Protein-protein interaction (PPI) network construction and analysis

According to the required confidence threshold set (combined score) >0.7, PPIs were analyzed using Search Tool for the Retrieval of Interacting Genes (23) and the PPI network was constructed by Cytoscape (24). CytoNCA (25), which is a plugin of Cytoscape for PPI network evaluation, was used to perform a topological analysis, including degree centrality (DC), betweenness centrality (BC) and closeness centrality (CC), for hub gene screening (26). Proteins were defined as nodes, and PPI associations were defined as edges in the PPI network.

Sub-module analysis

Based on the PPI network, the MCODE plugin of Cytoscape was used to screen specific bio-functional sub-modules (27) in the network. KEGG pathway analysis of these modules was performed using DAVID with the criterion of P<0.05. Proteins were defined as nodes, and PPI associations were defined as edges in modules screened from the PPI network.

Results

DEGs screening

Based on data preprocessing and Student's t-test, a total of 20,568 genes were identified, and 634 DEGs were identified between HCC samples and cirrhotic liver samples, of which 165 were upregulated and 469 were downregulated (Fig. 1).

Figure 1.

Volcano plot of differentially expressed genes. Red dots represent upregulated genes and blue dots represent downregulated genes.

Functional and pathway enrichment

To investigate the biological functions of the DEGs, GO_BP functional and KEGG pathway enrichment analyses were conducted. GO_BPs of upregulated DEGs were significantly enriched in cell division (P=2.91×10−10), mitotic cell cycle (P=5.57×10−10) and mitotic nuclear division (P=1.42×10−7), etc. (Fig. 2A). Furthermore, GO_BPs of downregulated DEGs were significantly enriched in small molecule metabolic process (P=1.87×10−21), xenobiotic metabolic process (P=7.15×10−15) and immune response (P=1.83×10−14), etc. (Fig. 2B). KEGG pathways of upregulated DEGs were significantly enriched in cell cycle (P=1.67×10−5), progesterone-mediated oocyte maturation (P=0.001) and p53 signaling pathway (P=0.003), etc. (Fig. 2C). In addition, KEGG pathways of downregulated genes were significantly enriched in chemical carcinogenesis (P=1.92×10−10), retinol metabolism (P=4.53×10−9) and complement and coagulation cascades (P=8.67×10−8), etc. (Fig. 2D).

Figure 2.

GO and KEGG enrichment analyses for up- and downregulated genes. The pink/purple columns represent upregulated genes and the blue columns represent downregulated genes. The color depth is negatively related to P-value. (A) GO analysis for upregulated genes; (B) GO analysis for downregulated genes; (C) KEGG analysis for upregulated genes; (D) KEGG analysis for downregulated genes. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network analysis

With required confidence >0.7, the PPI network for DEGs was constructed with 324 nodes and 915 edges (Fig. 3, Table I), including 74 upregulated genes and 250 downregulated genes. Downregulated Jun proto-oncogene, AP-1 transcription factor subunit (JUN; degree, 39), Fos proto-oncogene, AP-1 transcription factor subunit (FOS; degree, 34) and v-myc avian myelocytomatosis viral oncogene homolog (MYC; degree, 32) had high degrees and were therefore identified as the hub nodes in the PPI network.

Figure 3.

Protein-protein interaction network of differentially expressed genes. Pink nodes represent upregulated genes and blue nodes represent downregulated genes. Node size is positively related to degree.

Table I.

Proteins with a degree ≥20, as determined by topological analysis of the protein-protein interaction network.

Gene	DC	BC	CC
JUN	39	17,219.9	0.0283
FOS	34	19,632.8	0.0283
MYC	32	8,693.0	0.0282
EGR1	31	7,531.4	0.0282
CDK1	27	6,022.2	0.0280
CDKN1A	26	4,038.7	0.0280
CXCL12	25	3,844.4	0.0279
CYP2E1	22	2,721.3	0.0278
CYP1A2	22	4,144.9	0.0276
IGF1	21	6,096.5	0.0280
CCR7	21	2,227.9	0.0277
CYP3A4	21	1550.2	0.0277
PTGS2	20	6,968.9	0.0281
THBS1	20	5,690.0	0.0281
CCL5	20	3,738.3	0.0280
ANXA1	20	634.9	0.0275
SAA1	20	1,475.6	0.0275

DC, degree centrality; BC, betweenness centrality; CC, closeness centrality.

According to the selection criteria, four sub-modules were screened and named module 1, 2, 3 and 4 (Table II). In module 1 (Fig. 4A), there were 14 nodes with 91 edges, and according to the KEGG pathway analysis, the genes were significantly enriched in chemokine signaling pathway (P=3.14×10−13) and cytokine-cytokine receptor interaction (P=5.07×10−7), etc. In particular, downregulated C-X-C motif chemokine ligand 12 (CXCL12; degree, 13), C-C motif chemokine receptor 7 (CCR7; degree, 13) and C-C motif chemokine ligand 5 (CCL5; degree, 13) were the three nodes with the highest degrees in this module. In addition, there were 7 nodes with 21 edges identified in module 2 (Fig. 4B), and according to the KEGG pathway analysis, the genes were significantly enriched in complement and coagulation cascades (P=5.85×10−4) and p53 signaling pathway (P=0.04). In particular, downregulated IGF1 (degree, 6), plasminogen (degree, 6) and IGF2 (degree, 6) were the three nodes with the highest degrees in this module. In module 3, there were 7 nodes with 20 edges (Fig. 4C), and according to the KEGG pathway analysis, the genes were significantly enriched in retinol metabolism (P=5.57×10−13) and steroid hormone biosynthesis (P=1.11×10−5), etc. In particular, downregulated cytochrome P450 family 1 subfamily A member 2 (degree, 6), cytochrome P450 family 4 subfamily A member 11 (degree, 6) and cytochrome P450 family 3 subfamily A member 4 (degree, 6) were the three nodes with the highest degrees in this module. In addition, there were 9 nodes with 24 edges in module 4 (Fig. 4D), and according to the KEGG pathway analysis, the genes were significantly enriched in chemical carcinogenesis (P=2.31×10−16) and drug metabolism cytochrome P450 (P=8.13×10−10), etc. In particular, downregulated cytochrome P450 family 2 subfamily E member 1 (CYP2E1; degree, 6), cytochrome P450 family 2 subfamily C member 9 (CYP2C9; degree, 6) and cytochrome P450 family 2 subfamily A member 6 (CYP2A6; degree, 5) were the three nodes with the highest degrees in this module.

Table II.

Top 5 enriched Kyoto Encyclopedia of Genes and Genomes terms obtained from the subnet module analysis.

Term	Count	P-value	Proteins
Cluster 1
Chemokine signaling pathway	10	3.14×10−13	CCR7, GNAI1, CXCR4, CCL21, CXCL2 …
Cytokine-cytokine receptor interaction	7	5.07×10−7	CCR7, CXCR4, CCL21, CXCL9, CCL19 …
Rheumatoid arthritis	3	0.007870	CXCL6, CCL5, CXCL12.
NF-κB signaling pathway	3	0.008048	CCL21, CCL19, CXCL12.
Leukocyte transendothelial migration	3	0.014455	GNAI1, CXCR4, CXCL12.
Cluster 2
Complement and coagulation cascades	3	5.85×10−4	A2M, SERPINE1, PLG.
p53 signaling pathway	2	0.037772	SERPINE1, IGF1.
Cluster 3
Retinol metabolism	7	5.57×10−13	CYP3A4, CYP4A11, CYP2B6, UGT2B11, CYP26A1 …
Steroid hormone biosynthesis	4	1.11×10−5	CYP3A4, UGT2B11, CYP1A2, UGT2B7.
Chemical carcinogenesis	4	2.94×10−5	CYP3A4, UGT2B11, CYP1A2, UGT2B7.
Metabolic pathways	7	3.00×10−5	CYP3A4, CYP4A11, CYP2B6, UGT2B11, CYP26A1 …
Drug metabolism cytochrome P450	3	7.30×10−4	CYP3A4, CYP2B6, CYP1A2.
Cluster 4
Chemical carcinogenesis	9	2.31×10−16	GSTA2, CYP3A7, CYP2C19, CYP2C9, CYP2C18 …
Drug metabolism-cytochrome P450	6	8.13×10−10	GSTA2, CYP2C19, CYP2C9, CYP2C8, CYP2A6 …
Metabolism of xenobiotics by	5	2.30×10−7	GSTA2, CYP2C9, EPHX1, CYP2A6, CYP2E1.
cytochrome P450
Retinol metabolism	5	4.91×10−7	CYP3A7, CYP2C9, CYP2C18, CYP2C8, CYP2A6.
Linoleic acid metabolism	4	3.70×10−6	CYP2C19, CYP2C9, CYP2C8, CYP2E1.

Figure 4.

Subnet module analyses for differentially expressed genes in the protein-protein interaction network. Pink nodes represents upregulated genes and blue nodes represents downregulated genes. The node size is negatively related to its degree. (A) Module 1; (B) module 2; (C) module 3; and (D) module 4.

Discussion

Disease-specific differential gene expression reveals potential alterations associated with disease development. According to the analysis criteria of the present study, a total of 634 DEGs were identified in HCC versus cirrhotic tissue samples, of which 165 were upregulated and 469 were downregulated. After GO functional and KEGG pathway enrichment analyses, a PPI network was constructed with 324 nodes and 915 edges. In particular, JUN, FOS and MYC were the hub nodes in this network. Based on the PPI network, four specific modules were identified. According to the KEGG analysis results, module 1 was significantly enriched in chemokine signaling pathway, module 2 was significantly enriched in complement and coagulation cascades, module 3 was significantly enriched in retinol metabolism, and module 4 was significantly enriched in drug metabolism cytochrome P450.

The majority of the nodes in module 1 were chemokines or their receptors, including CXCL12, CCR7 and CCL5. Chemokines are cytokines that specifically respond to proinflammatory stimuli, and are involved in the migration of immune cells to damaged organs and are associated with HCC development (28). Chemokines are small molecules that can be divided into four groups (C, CC, CXC, and CXC3C) by the motifs of their NH2 terminals (29). CXCL12, also known as stromal cell-derived factor 1, is produced by HSCs, biliary epithelial cells and liver sinusoidal endothelial cells (30), and has a higher expression in cirrhotic liver tissue and HCC (31). It is well known that the CXCL12-CXCR4 axis serves an important role in the pathogenesis of liver diseases, including cirrhosis (32), migration (33), invasion (34) and HCC prognosis (35). However, in the present study, a significant downregulation was identified in HCC compared with in cirrhosis. Furthermore, Neve Polimeno et al (36) and Shibuta et al (37) demonstrated that the expression of CXCL12 is significantly reduced in hepatoma carcinoma cells and HCC compared with normal controls; however, the detailed mechanism of this reduction remains unclear. Therefore, it may be hypothesized that the CXCL12-CXCR4 axis serves various roles in different conditions of HCC originating from cirrhosis. CCR7 has been reported to be associated with the clinicopathological parameters of HCC (38). An in vitro study indicated that the expression of CCR7 is reduced in HCC cell lines compared with in the normal liver cell line L-02 (39). A similar result was identified in the present study; CCR7 expression was downregulated in HCC tissue samples compared with in cirrhotic tissue samples. Shi et al (40) reported that CCL1, which is exclusively expressed in non-hematopoietic stromal cells (41), may serve as a tumor suppressor by inhibiting CCR7-associated chemotaxis in HCC. Furthermore, overexpression of CCL21, the ligand of CCR7, in HCC cell lines presented a potential antitumor effect in a model of HCC (42); however, the detailed function of CCR7-CCL21 requires further investigation. CCL5 has been reported to be associated with the inflammatory cirrhosis stages in chronic liver disease (43). It has been demonstrated that CCL5 downregulation is able to inhibit the effects of the human bone marrow stromal cell line HS-5 on Huh-7 cell migration and invasion via the phosphatidylinositol 3-kinase/Akt pathway (44). These findings suggested that CCL5 may promote the migration and invasion of Huh-7 cells. In addition, Sadeghi et al (45) reported that CCL5 is upregulated in the serum of patients with HCC compared with patients with cirrhosis. Conversely, in the present study, a significant downregulation in CCL5 was detected in HCC tissue samples compared with in cirrhotic tissue. Therefore, it may be hypothesized that CCL5 expression differs between tissue and serum and in vitro and in vivo. This may be a novel opportunity for CCL5-targeted therapy.

CYPs belong to a superfamily of enzymes that serve important roles in the metabolism of procarcinogens, carcinogens and drugs (46). CYP2E1, CYP2C9 and CYP2A6 were three important CYPs identified in module 4 in the present study. CYP2E1 is often deficient in HCC cell lines (47). A significant decrease in CYP2E1 expression was identified in HCC samples compared with cirrhotic tissue samples in the present study. In addition, Kinoshita and Miyata (48) reported that CYP2E1 is significantly downregulated in HCC liver tissue. Wu et al (49) reported that CYP2E1 can be downregulated by resveratrol to attenuate diethylnitrosamine and 2-acetylaminofluorene-induced hepatocarcinogenesis in Sprague Dawley rats, thereby suggesting that CYP2E1 may have an inhibitory effect on HCC carcinogenesis. In addition, hepatitis B virus (HBV)-x protein can inhibit the expression of CYP2E1 via downregulating hepatocyte nuclear factor 4 (HNF4), resulting in the promotion of human hepatoma cell growth (48). Taken together, these results suggested that CYP2E1 may be a promising target for the drug-targeted therapy of HCC. CYP2C9 has been reported to be downregulated in HCC compared with in peri-HCC and normal control tissues (50); this result was similar to the findings of the present study. In addition, it has been reported that CYP2C9 can be suppressed by hsa-microRNA-128-3P, resulting in the invasion of HCC (51). Myung et al (52) also reported that CYP2C9 is involved in the sensitization of liver cancer stem cells to anticancer drugs via the signal transducer and activator of transcription 3 signaling pathway, particularly in advanced stages. However, the detailed mechanisms of these results remain unclear; therefore, further study is required. CYP2A6 is another important member of the CYP family that was downregulated in module 4. Similar to CYP2C9, CYP2A6 can be regulated by HNF4 and is involved in the metabolism of nitrosamines and aflatoxin B1 (53). Fushiya et al (54) reported that CYP2A6 is implicated in the metabolism of 5-fluorouracil, which is commonly used in HCC clinical chemotherapy. In addition, a previous study indicated that CYP2A6 was downregulated in HBV- and hepatitis C virus-infected livers, and in HCC (55); the present study suggested that HCC samples had a significantly lower expression of CYP2A6 compared with cirrhotic liver samples. Furthermore, a previous study reported that CYP2A6 activity was decreased in moderate or severe alcoholic liver diseases, but not in mild severe alcoholic liver disease (56). Therefore, it may be hypothesized that CYP2A6 may be negatively correlated with the stage of HCC development. Further studies are required to investigate the CYP2A6-associated functions and pathways in HCC.

Although numerous DEGs with their potential functions were identified between HCC and cirrhosis samples in an in silico analysis, there remain limitations to the present study. For example, in silico analysis of the obtained results revealed that numerous chemokines and cytokines were involved in HCC development; however, the involvement of the identified genes in vitro remains unknown. Therefore, further experimental verifications of these genes are required. Furthermore, these analyses were based on subjective criteria; therefore, some genes, which have lower degrees but important roles in HCC carcinogenesis, may have been overlooked. Despite these limitations, these analytical results may provide novel insights into the mechanism of HCC originating from cirrhosis.

In conclusion, the present study identified a series of DEGs between HCC and cirrhotic tissue samples. Based on the GO and KEGG enrichment analyses of DEGs in the PPI network, chemokines, such as CXCL12, CCR7, CCL5, and cytokines, such as CYP2E1, CYP2C9, CYP2A6, were identified as two important components in the process of HCC developing from cirrhosis as they are associated with the regulation of inflammation, growth and the invasion of pre-cancerous cells in the liver; thus, they may serve crucial roles in HCC development. However, these results were derived from bioinformatics analysis; therefore, the effects of these genes in HCC in vivo require further verification.

Competing interests

The authors declare that they have no competing interests.