Identification of key pathways and biomarkers in sorafenib-resistant hepatocellular carcinoma using bioinformatics analysis.

Hepatocellular carcinoma (HCC) is one of the most common malignant types of cancer, with a high mortality rate. Sorafenib is the sole approved oral clinical therapy against advanced HCC. However, individual patients exhibit varying responses to sorafenib and the development of sorafenib resistance has been a new challenge for its clinical efficacy. The current study identified gene biomarkers and key pathways in sorafenib-resistant HCC using bioinformatics analysis. Gene dataset GSE73571 was obtained from the Gene Expression Omnibus (GEO) database, including four sorafenib-acquired resistant and three sorafenib-sensitive HCC phenotypes. Differentially expressed genes (DEGs) were identified using the web tool GEO2R. Functional and pathway enrichment of DEGs were analyzed using the Database for Annotation, Visualization and Integrated Discovery and the protein-protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins and Cytoscape. A total of 1,319 DEGs were selected, which included 593 upregulated and 726 downregulated genes. Functional and pathway enrichment analysis revealed DEGs enriched in negative regulation of endopeptidase activity, cholesterol homeostasis, DNA replication and repair, coagulation cascades, insulin resistance, RNA transport, cell cycle and others. Eight hub genes, including kininogen 1, vascular cell adhesion molecule 1, apolipoprotein C3, alpha 2-HS glycoprotein, erb-b2 receptor tyrosine kinase 2, secreted protein acidic and cysteine rich, vitronectin and vimentin were identified from the PPI network. In conclusion, the present study identified DEGs and key genes in sorafenib-resistant HCC, which further the knowledge of potential mechanisms in the development of sorafenib resistance and may provide potential targets for early diagnosis and new treatments for sorafenib-resistant HCC.

Introduction

Hepatocellular carcinoma (HCC) is a malignant tumor that has become the third leading cause of cancer-associated cases of mortality (1). Sorafenib, an oral multitarget tyrosine kinase inhibitor, targets various molecular mechanisms, including tumor growth and angiogenesis (2). It is the only systemic therapy drug for HCC that is approved by the USA Food and Drug Administration and as such has been applied in the clinic extensively (3). The clinical efficacy of sorafenib is limited and patients face poor prognosis; no difference in recurrence-free survival between sorafenib and a placebo-controlled group has been reported (4). Furthermore, time-to-tumor progression and overall survival were not observed to be different for patients treated with transarterial chemoembolization (TACE) in conjunction with sorafenib over the TACE placebo group in a stent-protected angioplasty vs. carotid endarterectomy trial (5). Several potential mechanisms of sorafenib resistance were proposed. The processes of epithelial-mesenchymal transition (EMT) and mesenchymal-epithelial transition, along with critical growth factors and signaling pathways, exhibit an impact on sorafenib resistance, including activation of the phosphatidylinositol 3-kinase/Akt signaling pathway (6). Additionally, cancer stem cells (CSCs) and CSC-like cells, which retain strong proliferation ability, multi-directional differentiation capacity and high drug-resistance properties, may not be completely cleared by sorafenib but differentiate and develop into novel cancer tissues, resulting in the metastasis and recurrence of HCC (6). Biological processes involved in tumor microenvironment, inflammation, fibrosis, angiogenesis, hypoxia, autophagy, viral reactivation and oxidative stress may serve a pivotal role in the resistance to sorafenib (6). However, the resistance mechanisms for sorafenib remain unclear and novel research may provide insight into the discovery of an effective treatment or personalized therapy for advanced HCC.

High-throughput microarray technology has been widely used to analyze the gene expression data of various cancer types (7). It has been a promising method to screen for potential biomarkers in tumor diagnosis and pathways involved in tumorigenesis and drug resistance (7–9). In the current study, microarray data for GSE73571 facilitated the investigation of differentially expressed genes (DEGs) in sorafenib-sensitive and sorafenib-resistant tumors. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed and a protein-protein interaction (PPI) network was constructed, identifying hub genes. The bioinformatics analysis of crucial genes or pathways in sorafenib-resistant HCC revealed potential strategies for improving clinical efficacy.

Materials and methods

Microarray data

Gene expression profile data (GSE73571) were obtained from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo). Four sorafenib-acquired resistant HCC and three sorafenib-sensitive phenotypes were included. The array data were acquired from Affymetrix Human Gene 1.0 ST Array [GPL6244; transcript (gene) version].

DEG analysis

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) was used to compare two or more groups of samples in a GEO series to identify genes that were differentially expressed under the same experimental conditions. DEGs in resistant and sensitive samples were analyzed by GEO2R. |log2FC|≥0.4 and P<0.01 were used as cut-off criteria and defined a statistically significant difference (10). A heat map of DEGs was generated using HemI 1.0 (https://sourceforge.net/projects/mev-tm4/).

Functional and pathway enrichment analysis

GO enrichment and KEGG pathway analysis of the screened DEGs was performed using the Database of Annotation Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/). DAVID contains a series of functional annotation programs to explore abundant biological messages of genes and information mapped in DAVID was important for the completion of high-throughput gene functional analysis (11). P<0.05 was considered to indicate a statistically significant difference (12).

PPI network analysis

PPI network analysis was performed using the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/). STRING provides information associated with predicted and proven interactions between large numbers of proteins. Identified DEGs were imported into STRING. Genes with a combined score of >0.4 were identified as significant. The PPI network was built using Cytoscape (http://www.cytoscape.org/). In addition, higher-degree nodes were regarded as hub nodes (12). Sub-modules of the PPI network were analyzed by Molecular Complex Detection (MCODE; Version 1.4.2; by Bader Lab, department of Biochemistry, University of Toronto; Toronto, Canada) (13), with the criteria set as follows: number of nodes >4 and MCODE score >3. GO enrichment and KEGG pathway analysis of DEGs of the sub modules was finished by DAVID.

Results

According to the cut-off criteria (P<0.05 and |log2FC|≥2), a total of 1,319 DEGs from sorafenib-resistant and sorafenib-sensitive specimens were identified, with 593 up- and 726 downregulated genes. A heat map of the top 50 up- and downregulated DEGs is presented in Fig. 1.

Figure 1.

Heat map of differentially expressed genes (the top 50 up- and downregulated genes; P<0.05). Red indicates a higher and green a lower expression level.

GO enrichment analysis

GO enrichment analysis was performed by DAVID. GO biological processes analysis revealed that upregulated DEGs were associated with negative regulation of endopeptidase activity, cholesterol homeostasis and fibrinolysis and downregulated DEGs were associated with DNA replication and repair (Table I). For molecular function, upregulated genes were enriched in collagen and receptor binding and serine-type endopeptidase activity, while downregulated genes were enriched in poly(A) RNA binding, helicase activity and ATP binding. Additionally, GO cell component analysis revealed that upregulated genes were primarily located in the extracellular space, extracellular exosome and extracellular region and the location of downregulated DEGs was primarily in nucleolus, nucleoplasm and centrosome.

Table I.

Gene ontology analysis of upregulated and downregulated genes in sorafenib-resistant HCC.

A, Upregulated expression
Category	Term	Count	%	P-value
GOTERM_CC_DIRECT	GO:0005615/extracellular space	110	16.492	<0.001
GOTERM_CC_DIRECT	GO:0070062/extracellular exosome	177	26.537	<0.001
GOTERM_CC_DIRECT	GO:0005576/extracellular region	100	14.993	<0.001
GOTERM_CC_DIRECT	GO:0072562/blood microparticle	21	3.148	<0.001
GOTERM_CC_DIRECT	GO:0016021/integral component of membrane	235	35.232	<0.001
GOTERM_BP_DIRECT	GO:0010951/negative regulation of endopeptidase activity	18	2.699	<0.001
GOTERM_BP_DIRECT	GO:0042632/cholesterol homeostasis	13	1.949	<0.001
GOTERM_BP_DIRECT	GO:0042730/fibrinolysis	8	1.199	<0.001
GOTERM_BP_DIRECT	GO:0002576/platelet degranulation	15	2.249	<0.001
GOTERM_BP_DIRECT	GO:0030198/extracellular matrix organization	20	2.999	<0.001
GOTERM_MF_DIRECT	GO:0005518/collagen binding	11	1.649	<0.001
GOTERM_MF_DIRECT	GO:0005102/receptor binding	29	4.348	<0.001
GOTERM_MF_DIRECT	GO:0004252/serine-type endopeptidase activity	24	3.598	<0.001
GOTERM_MF_DIRECT	GO:0001948/glycoprotein binding	10	1.499	<0.001
GOTERM_MF_DIRECT	GO:0008131/primary amine oxidase activity	4	0.600	<0.001
Extracellular space means GO enrichment of cellular primarily enriches in the extracellular space.
B, downregulated expression
Category	Term/gene function	Count	%	P-value
GOTERM_MF_DIRECT	GO:0044822/poly(A) RNA binding	99	18.099	<0.001
GOTERM_MF_DIRECT	GO:0004386/helicase activity	20	3.656	<0.001
GOTERM_MF_DIRECT	GO:0005524/ATP binding	89	16.271	<0.001
GOTERM_MF_DIRECT	GO:0004004/ATP-dependent RNA helicase activity	14	2.559	<0.001
GOTERM_MF_DIRECT	GO:0008026/ATP-dependent helicase activity	9	1.645	<0.001
GOTERM_CC_DIRECT	GO:0005730/nucleolus	75	13.711	<0.001
GOTERM_CC_DIRECT	GO:0005654/nucleoplasm	152	27.788	<0.001
GOTERM_CC_DIRECT	GO:0005813/centrosome	43	7.861	<0.001
GOTERM_CC_DIRECT	GO:0005634/nucleus	199	36.380	<0.001
GOTERM_CC_DIRECT	GO:0005814/centriole	16	2.925	<0.001
GOTERM_BP_DIRECT	GO:0006260/DNA replication	25	4.570	<0.001
GOTERM_BP_DIRECT	GO:0006281/DNA repair	27	4.936	<0.001
GOTERM_BP_DIRECT	GO:0000732/strand displacement	10	1.828	<0.001
GOTERM_BP_DIRECT	GO:0000731/DNA synthesis involved in DNA repair	11	2.011	<0.001
GOTERM_BP_DIRECT	GO:0006974/cellular response to DNA damage stimulus	23	4.205	<0.001

BP, biological process; CC, cellular component; MF, molecular function; GO, gene ontology.

KEGG pathway analysis

To gain a deeper understanding of significant DEGs, a pathway enrichment analysis was performed. As presented in Table II, upregulated genes were enriched in pathways of coagulation cascades, insulin resistance and metabolic pathways, while downregulated DEGs were significantly associated with pathways of RNA transport and cell cycle.

Table II.

KEGG pathway analysis of upregulated and downregulated genes in sorafenib-resistant hepatocellular carcinoma.

A, Upregulated genes
Name	Gene count	%	P-value	Genes
hsa04610: Complement and coagulation cascades	20	2.999	<0.001	KNG1, F12, F10, CFB, C3, C1R, SERPING1, F7, C1S, C8B, FGG, FGB, SERPINF2, KLKB1, SERPINA5, SERPIND1, C2, CFI, PLAU, F2R
hsa04142: Lysosome	14	2.099	0.003047	TCIRG1, CTSZ, CLTB, ACP5, ACP2, CTSA, CTSS, CTSL, LAPTM5, TPP1, ARSA, CTSE, SMPD1, CTSD
hsa04931: Insulin resistance	13	1.949	0.003298	SREBF1, NR1H2, PPP1R3C, PPP1R3B, SOCS3, SLC2A2, PRKAG2, PRKAB2, NFKBIA, CREB3L3, SLC27A3, CPT1A, PCK1
hsa01100: Metabolic pathways	74	11.094	0.006009	TM7SF2, ALAD, OGDHL, ADH1C, ADH1A, AGXT, CKB, PTGIS, P4HA2, ST3GAL5, MAT1A, RGN, DAO, ITPK1, AGPAT3, AGPAT2, GATM, ACADS, SPTLC3, FAXDC2, PISD, PNPLA3, PMM1, ALDH3B1, DGAT2, CYP27A1, HAO2, EXT1, AOC1, AKR1D1, LCT, AOC3, ACAA1, GALNT2, ENPP7, ASS1, HSD17B2, ALDOC, ENPP3, CERS4, PLPP3, PIPOX, THTPA, B3GNT5, DHCR7, PLA2G12B, PEMT, HAAO, ENO3, ETNK2, FUT2, BDH1, GAL3ST1, ACSL5, CYP19A1, TCIRG1, MOGAT3, ST6GAL1, NDUFA2, SI, MAOB, HOGA1, AMPD3, IDNK, PCK1, POLD4, BAAT, HMGCS2, SMPD1, LIPC, PHYKPL, PON3
hsa04975: Fat digestion and absorption	7	1.049	0.007472	APOA4, ABCG8, MOGAT3, DGAT2, PLA2G12B, PLPP3, AGPAT2
B, Downregulated genes
Name	Gene count	%	P-value	Genes
hsa03008: Ribosome biogenesis in eukaryotes	17	3.108	<0.001	XPO1, GTPBP4, UTP6, UTP15, GNL3L, BMS1, DROSHA, WDR36, POP1, WDR3, NOP58, GNL2, MDN1, RBM28, SPATA5, WDR43, GNL3
hsa03013: RNA transport	17	3.108	<0.001	XPOT, XPO1, XPO5, ALYREF, EIF5B, NUPL2, NUP155, PNN, TRNT1, UPF3B, SEH1L, NUP205, POP1, EIF3J, RANBP2, TPR, GEMIN5
hsa03460: Fanconi anemia pathway	9	1.645	<0.001	FANCM, BLM, FANCD2, FANCI, BRCA2, BRIP1, ATR, FANCA, BRCA1
hsa04742: Taste transduction	7	1.280	<0.001	TAS2R13, TAS2R14, TAS2R19, TAS2R46, TAS2R43, TAS2R20, TAS2R10
hsa04110: Cell cycle	9	1.645	0.014438	CDC7, CDC6, RBL1, PRKDC, ATR, MYC, CDC25A, ATM, STAG1

KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network construction and module screening

The PPI network of DEGs is presented in Fig. 2. The network was composed of 279 nodes and 636 edges. Degrees >10 were set as the cutoff criterion, from which eight genes were selected as hub genes. Hub genes included kininogen 1 (KNG1), vascular cell adhesion molecule 1 (VCAM1), apolipoprotein C3 (APOC3), alpha 2-HS glycoprotein (AHSG), erb-b2 receptor tyrosine kinase 2 (ERBB2), secreted protein acidic and cysteine rich (SPARC), vitronectin (VTN) and vimentin (VIM). The significant module of DEGs with the highest score was selected using the plug-in MCODE program. This module included 8 nodes and 28 edges (Fig. 3). No enrichment of GO terms and KEGG pathways was observed in this module.

Figure 2.

Protein-protein interaction network of differentially expressed genes.

Figure 3.

Significant module from the protein-protein interaction network.

Discussion

In the present study, gene expression data of GSE73571 were extracted from the GEO database. Four sorafenib-resistant and three sorafenib-sensitive HCC samples were selected for analysis. A total of 593 downregulated and 726 upregulated DEGs were identified among the sorafenib-sensitive and sorafenib-resistant phenotypes with HCC. Function annotation revealed that these DEGs were primarily associated with complement and coagulation cascades, DNA replication, synthesis and repair. PPI network analysis suggested that eight hub genes exhibited higher degrees of interaction, which may describe new targets in sorafenib resistance.

GO term analysis reveled that upregulated DEGs were primarily associated with negative regulation of endopeptidase activity, cholesterol homeostasis and fibrinolysis. Peptidases are crucial in tumor formation and development in various ways, including regulating the process of neoplastic growth, as adhesion molecules, through involvement in intracellular signaling and extracellular matrix degradation (14). Various studies indicated that the expression of various peptidases, including circulating aminopeptidase N/CD13, varies among cancer types at different stages (15–18). Cholesterol metabolites are related to the development of various cancers (19,20). Management of cholesterol homeostasis may lead to lower risks of non-Hodgkin lymphoma, melanoma, breast and endometrial cancers (21–23). The fibrinolytic system was verified to promote tumor growth and it may be involved in apoptosis inhibition, angiogenesis, cell proliferation and extracellular matrix (24,25). Downregulated DEGs were primarily involved in DNA replication and repair. DNA replication and repair may be involved in the development of certain types of cancer, by regulating the cell cycle, angiogenesis, cell differentiation and cell signaling (26–29). Additionally, KEGG pathway analysis of upregulated genes revealed involvement in complement and coagulation cascades, insulin resistance and metabolic pathways. KEGG pathway analysis of downregulated DEGs exhibited participation in RNA transport and cell cycle. A previous study reported that constituents of the coagulation cascade could affect cancer progression (30) and the complement and coagulation cascade pathway enrichment in colorectal cancer (CRC) and HCC has been described (31,32). Poor glycemic control is a prognostic element in patients with HCC and diabetes (33) and insulin resistance is considered as a risk factor in HCC development (34). Metabolic pathways, cell cycle and RNA transport have been associated with carcinogenesis, cancer cell survival and growth (35–37). Therefore, the identification of the above mentioned pathways may lead to novel prognostic and therapeutic methods in HCC and sorafenib-resistant HCC.

In the present study, the following eight hub genes were identified through PPI network construction and analysis: KNG1, VCAM1, APOC3, AHSG, ERBB2, SPARC, VTN and VIM. KNG1, a cysteine proteinase inhibitor, exhibited the highest degree of connectivity in the PPI network. It participates in blood coagulation, inflammatory response, apoptosis regulation and the prevention process of metastasis in cancer cells (38,39). KNG1 is overexpressed in CRC and HCC (40) and is regarded as a potential prognostic marker for CRC, as patients with increased KNG1 expression exhibit decreased survival rates when compared with lower KNG1 expression patients (39,41). Additional studies reported that KNG1 may be involved in cholesterol and lipoprotein metabolism disturbance and take part in the complex interplay between the hemostatic system and immune response in hepatitis C virus (HCV)-associated HCC (39,40). Differences in KNG1 expression were observed among patients with multidrug-resistant and drug-sensitive tuberculosis and healthy control patients, and it is suspected that blood coagulation serves a role in the resistant ability of KNG1 (42).

The second hub gene is VCAM1, a mediator of angiogenesis, which serves a critical function in endothelium development during angiogenesis (43,44). It promotes cancer cell adhesion to the endothelium and is associated with immune and inflammatory responses of tumors (45). High VCAM1 serum levels have been reported for various cancer types, including HCC, chronic liver disease, breast cancer and CRC, and VCAM1 is considered a potential predictor for cancer prognosis (46–50). The amount of VCAM1 in the serum is dependent on tumor stage and neoplasm metastasis (51). VCAM1 expression decreased following carboplatin intervention in mice with platinum-sensitive ovarian cancer, while high-level expression was maintained in mice with platinum-resistant tumors (52). VCAM1 has been recognized as a response monitor in the treatment of ovarian cancer and as a molecular biomarker in chemotherapy-associated sensitivity, which allows for earlier alterations in treatment decisions (52).

APOC3, the third hub gene, is expressed in the liver and participates in very low-density lipoprotein (VLDL) triglyceride (TG) metabolism by inhibiting VLDL-TG clearance in vivo and promoting absorption of intestinal TG and VLDL-TG production (53). APOC3 polymorphism is considered to be an independent risk factor of hepatocarcinogenesis and HCC development in patients with chronic hepatitis B (54). It was suggested that APOC3 may be involved in HCC familial aggregation in China (55). There is no evidence indicating an association between APOC3 and sorafenib or other drug resistance, making it a potential target for further research.

AHSG, a serum glycoprotein, is produced by hepatocytes and is involved in several metabolic disorders, including non-alcoholic fatty liver disease (56,57), genesis of diabetes (58) and metabolic syndrome (59). It is significantly increased in HCC when compared with normal patients. AHSG is also an important biomarker of recent mortality in liver cirrhosis and liver cancer patients (60). AHSG is associated with cancer progression through regulation of the transforming growth factor-β signaling pathway (61,62). In the present study, dysregulation of AHSG was observed, which is in accordance with a previous study describing that the upregulation of AHSG in HCC drug-resistant cell lines may be a predictor for HCC with chemotherapeutic drug resistance (63). The levels of AHSG and five other selected serum biomarkers were proposed to predict resistance to neoadjuvant chemotherapy in breast cancers (64). The differences in complement system and LDL oxidation may contribute to this phenomenon (64).

Another hub gene, ERBB2, encodes for a receptor tyrosine kinase, a member of the epidermal growth factor receptor (EGFR) family, and is involved in fixation and propagation of oncogene mutations (65). Comprehensive studies have demonstrated that increased ERBB2 expression in breast and ovarian cancer is associated with poor prognosis (66,67). ERBB2 was upregulated in the current study and has been regarded as a potential critical regulator for malignant transformations in early HCC (68). In specific studies on HCC, the percentage of control cases expressing ErbB2 ranged from 0–30% (69–72). ERBB2 expression level potentially refers to different HCC stages (71,73). Cetuximab and panitumumab, EGFR-targeted antibodies, are treatments for metastatic CRC (74). The amplification of ERBB2 contributes to the primary (de novo) resistance to anti-EGFR treatment, the mechanism of which is associated with the activation of the MEK-ERK cascade (74).

The other three identified hub genes in the current study were SPARC, VTN and VIM. SPARC encodes a protein involved in numerous biological processes, which are associated with various cancer mechanisms, including development, cell apoptosis, angiogenesis, cell differentiation, cell proliferation, cell adhesion and migration (75,76). Lau et al (77) reported high SPARC expression in HCC compared with non-tumorous liver. Elevated SPARC expression promotes tumor aggressiveness of melanoma (78), glioblastoma (79) and prostate cancer (80). The association between high SPARC mRNA expression and low pathological response rate in breast tumor cells following neoadjuvant anthracycline treatment has been confirmed (81). VTN, a representative acute-phase glycoprotein, is expressed and secreted by hepatocytes (82). Interactions with integrins may enhance cell adhesion and the spread in serum and extracellular matrix (82). Increased VTN was described as a poor prognostic tool for patients with HCC as it is a primary component of the stroma and accelerates leucocyte accumulation (83) and cell migration (84). VTN has an adverse effect on HCC development in HCV-infected patients with liver cirrhosis (85). Additionally, VTN adhesion was previously suggested to contribute to drug resistance in chemotherapy-treated myeloma cells through Notch signaling activation (86). Elevated VIM expression was observed in small size HCC (≤2 cm) (87). It has been suggested that the circulating level of VIM may be more sensitive and specific compared with AFP in detecting small tumors and VIM expression has been recognized as a potential biomarker for HCC diagnosis (87). VIM is overexpressed in CRC cells with butyrate or histone deacetylase inhibitor resistance, which is likely to mediate cell signaling/gene expression cascades and to integrate the EMT process in drug-resistant CRC cells (88).

The genes identified in the current study may be associated with HCC genesis, development and prognosis and potentially contribute to sorafenib resistance in HCC. Future research into these genes may provide information on the mechanism of sorafenib-resistant HCC and hold potential for novel therapeutic methods for drug-resistant HCC.

In conclusion, a comprehensive bioinformatics analysis of DEGs was performed, revealing genes that may be involved in the biological process of sorafenib-resistant HCC. A series of potential targets has been provided, which may be considered for future exploration.