Comprehensive Network Analysis of the Molecular Mechanisms Associated with Sorafenib Resistance in Hepatocellular Carcinoma


 Background: This research aimed to investigate the potential molecular mechanism of sorafenib resistance to hepatocellular carcinoma (HCC). Methods: Differential expression analysis were performed to identified differentially expressed genes (DEGs) in sorafenib resistant HCC. Then, a series of bioinformatic analysis were performed to explore the potential crucial molecules in sorafenib resistant HCC. For example, gene function annotation, pivot regulators prediction, ROC analysis and survival analysis. Results: There were 827 differentially expressed genes were identified. Moreover, most of the differentially expressed genes are involved in immune and inflammatory-related functions and signaling pathways. Also, 18 transcription factors were predicted to regulate the transcription factors of differentially expressed genes, which play an essential role in the regulation of dysfunctional gene networks. In target genes of transcription factors, CDK1 and CDKN1A have high diagnostic value in the resistance of hepatocellular carcinoma to sorafenib. Conclusions: TAPBP has the strongest correlation with drug resistance of hepatocellular carcinoma and the highest diagnostic efficiency. In addition, CDK1 and CDKN1A have high diagnostic value in the resistance of hepatocellular carcinoma to sorafenib. Overall, our analysis shows that a large number of gene disorders occur during the development of resistance to sorafenib in hepatocellular carcinoma, and they are associated with immune and inflammatory reactions in the body. These results provide critical theoretical references for the pathogenesis and diagnosis of sorafenib resistance.

transcription factors of differentially expressed genes, which play an essential role in the regulation of dysfunctional gene networks. In target genes of transcription factors, CDK1 and CDKN1A have high diagnostic value in the resistance of hepatocellular carcinoma to sorafenib. Conclusions: TAPBP has the strongest correlation with drug resistance of hepatocellular carcinoma and the highest diagnostic efficiency. In addition, CDK1 and CDKN1A have high diagnostic value in the resistance of hepatocellular carcinoma to sorafenib. Overall, our analysis shows that a large number of gene disorders occur during the development of resistance to sorafenib in hepatocellular carcinoma, and they are associated with immune and inflammatory reactions in the body. These results provide critical theoretical references for the pathogenesis and diagnosis of sorafenib resistance.

Background
As one of the most common malignant tumors in the world, hepatocellular carcinoma (HCC) has an inferior prognosis, leading to its second leading cause of cancer death [1].
Surgical resection is the primary treatment for early hepatocellular carcinoma [2]. 3 Because most HCC patients do not qualify for curative surgery in the late stage of the disease at diagnosis, the availability and efficacy of treatment options are limited [3]. At present, the 5-year survival rate of HCC patients is less than 20% [4]. Phase III clinical trials have shown that sorafenib, a multitarget tyrosine kinase inhibitor that reduces tumor cell proliferation and angiogenesis, can improve overall survival in patients with advanced HCC [5,6]. Sorafenib is a chemotherapeutic drug approved by the Food and Drug Administration (FDA) [7]. Although sorafenib is currently the only drug approved for the treatment of hepatocellular carcinoma, its therapeutic effect is affected by a variety of regulatory mechanisms [8,9].
Sorafenib is an oral multikinase inhibitor that blocks the proliferation of cancer cells by inhibiting the serine/threonine kinase isoforms of Raf, Raf-1, and B-Raf. Thus, it inhibits the mitogen-activated protein kinase/extracellular signal-regulated kinase (ERK) signaling pathway, and decreases the expression of cyclin D1 and stagnates the cell cycle [10]. The role of sorafenib in HCC cell death is attributed to the down-regulation of Mcl-1 (myeloid leukemia-1) [11], a member of the Bcl-2 family. Sorafenib blocked Erk-mediated phosphorylation of Mcl-1 on Thr92, making Mcl-1 more stable [12]. On the other hand, the inhibition of Snail expression by sorafenib via MAPK signaling pathway in HCC cells showed a productive inhibitory activity on EMT [13]. However, it has also been reported that in sorafenib-resistant HCC cells, EMT is associated with activation of phosphoinositol 3-kinase (PI3K) / AKT pathway [14]. In recent years, various clinical studies have shown that sorafenib has anti-tumor effects not only in HCC but also in other types of cancer, including thyroid cancer, renal cell cancer and prostate cancer [15,16].
Although sorafenib therapy has been shown to prolong the overall survival of HCC patients, only 2% of patients showed partial response to treatment based on the RECIST criteria (substantial tumor response assessment criteria) [17]. The low response rate was 4 attributed to the intrinsic resistance of HCC to the toxicity of sorafenib [18]. Owing to the heterogeneity of HCC, some patients initially developed resistance to sorafenib, but in most cases, resistance was acquired through long-term exposure to drugs [19]. So far, however, there has been the no more systematic report on sorafenib resistance.
With the development of microarray technology, it is easier to identify global genetic changes and their functions in many tumors. Therefore, in this study, we combined several sets of sequencing data of hepatocellular carcinoma samples with sorafenib resistance and explored the molecular mechanisms related to drug resistance. Through a series of bioinformatics analysis, key regulators and potential signaling pathways related to sorafenib resistance were explored.

Data resources
Gene Expression Omnibus (GEO)is an international public knowledge database, which included high-throughput gene expression data submitted by research institutions around the world. The gene expression profiles in HCC xenograft mice treated with sorafenib (GSE73571) were downloaded from GEO database for differential expression analysis, including 4 sorafenib-acquired resistant and 3 sorafenib sensitive [20]. In addition, GSE62813) and TCGA HCC dataset were used to performed ROC analysis and survival analysis [14].

Differential expression analysis
The difference analysis of gene expression profile data is mainly realized by using the limma package of R language [21][22][23], which includes preliminary screening and correction of data. Firstly, the correct background function was used to calibrate and standardize the data. Then the quantile normalization method in the normalize Between Arrays function is used to filter out the control probes and the low expression probes in the probe data.

5
Finally, eBayes and lmFit functions were used to identify and screen differentially expressed genes in the dataset (p<0.05, |logFC|>0.5). In addition, In order as following formula:

GO function and KEGG pathway enrichment analysis
In this study, we analyzed the differentially expressed genes between sorafenib resistance and sorafenib-sensitive hepatocellular carcinomas, and enriched the functions and pathways by using Bioconductor's Cluster profiler package [24] in R language, including Go function (p-value cutoff = 0.01, q-valueCutoff = 0.01) and KEGG pathway (p-value cutoff = 0.05, q-valueCutoff = 0.2). We construct corresponding functional and pathway networks to identify the functions and pathways of gene participation.

Regulator prediction
Transcription factor (TF), as the core driving force of transcriptional and posttranscriptional regulation, has potential effects on sorafenib resistance. Therefore, we predict pivot regulators of differentially expressed genes. A pivot node is a node that not only interacts with two modules but also has at least two pairs of interactions with each module. The hypergeometric test significance analysis of the interaction between the node and each module was p < 0.01. Python program is compiled to find pivot nodes of interaction module genes for further analysis.

Correlation analysis
Pearson correlation analysis was used to identify the linear correlation between the expression of differentially expressed genes. The correlation coefficients of two continuous variables (X, Y) are generally expressed by r, where n is the sample size. The greater the absolute value of r, the stronger the correlation. 6

Dysregulated genes in sorafenib-resistant hepatocellular carcinoma
Resistance to sorafenib in patients with hepatocellular carcinoma reduces their prognostic survival time. Therefore, exploring the molecular mechanism of drug resistance is of great significance for the treatment of hepatocellular carcinoma. First of all, Principal Component Analysis (PCA) was performed on samples to detect abnormal samples (Fig.   1a). Then, differential expression analysis was performed to identified differentially expressed genes (DEGs) in sorafenib-resistant samples compared with sorafenib-sensitive samples. A total of 827 significantly DEGs were identified (p<0.05, |logFC|>0.5) (Fig. 1b ,   Table S1), including 475 up-regulated and 352 down-regulated. It was considered that these DEGs are partly related to sorafenib resistance in HCC. And the top 5 DEGs were further filter by PI value, namely FXYD3, FXYD1, TBX4, NPNT, and ADH1C (Fig. 1c).

Gene function annotation
Gene-set enrichment analysis (GESA) for all genes and GO function and KEGG pathway enrichment analysis for DEGs were performed, respectively. It has been showed that genes in sorafenib-resistant HCC were significantly participated in various signaling pathways associated with HCC and drug metabolism pathways (Fig. 2a), such as PPAR signaling pathway, TGF-β signaling pathway, JAK-STAT signaling pathway, MAPK signaling pathway and drug metabolism cytochrome p450, etc.

7
The differentially expressed genes were significantly enriched in 1137 biological process entries, 173 cell component entries, 189 molecular function entries, and 28 KEGG pathway entries (Table S2). It was found that these genes were significantly involved in neutrophil associated biological processes, dephosphorylation, regulation of vesicle-mediated transport, etc. (Fig. 2b). The results of KEGG pathway enrichment analysis showed that DEGs mainly participated in MAPK signaling pathway, sphingolipid signaling pathway, phagosome, etc. (Fig. 2c). Further, pathways which enriched by DEGs with PI value > 2 were focused (Fig. 2d), and for subsequent analysis.

Transcription factors significantly regulating dysregulated genes
It is well known that the expression, transcription, and post-transcription of genes in vivo are strictly regulated. The transcription and post-transcriptional regulation of dysregulated genes may be closely related to the occurrence and development of various diseases. Here, TFs significantly regulating the drug resistance-related genes were predicted to help us to further explore the transcriptional regulation mechanism of sorafenib resistance in hepatocellular carcinoma. There were 18 TFs were identified (Fig.   3a) and differentially expressed in sorafenib-resistant HCC. It was observed that NFKB1 regulated the most DEGs, suggesting NFKB1 may affect the resistance of hepatocellular carcinoma to sorafenib and play an essential regulatory role in the process of drug resistance in hepatocellular carcinoma. Further, the strong correlation between TFs and DEGs indicated that the prediction results are more authentic (Fig. 3b). Moreover, these TFs (correlation coefficient >0.85) were identified as critical regulatory factors associated with drug resistance in hepatocellular carcinoma and played an essential role in the whole regulatory network (Fig. 3c).

Effect of crucial disorders on drug resistant hepatocellular carcinoma
To confirm potential prognostic efficacy of critical molecules for drug resistant HCC, ROC 8 analysis and survive analysis were performed on two GEO datasets (GSE73571 and GSE62813) and TCGA dataset of HCC, respectively. ROC curve showed that DEG CDK1, CDKN1A, TAPBP and MYC with high AUC (AUC>=0.8, Fig. 4a), suggesting the correlation between these molecules and disease. Survive analysis also showed the four DEGs tightly associated with HCC survival and prognosis (p<0.05, Fig. 4b).

Discussion
Since the clinical approval of sorafenib, sorafenib has become a new standard therapy for advanced HCC patients [25]. Compared with untreated patients, sorafenib prolonged the overall survival of HCC patients by only 2.8 months [26,27]. Because of its side effects of drug resistance, there is a lack of clinical efficacy [28]. Therefore, it is of considerable significance to explore the resistance of sorafenib in the treatment of hepatocellular carcinoma. In our analysis, 827 DEGs in sorafenib-resistant hepatocellular carcinoma were found. These disordered genes may be the essential molecules that lead to or mediate the resistance of hepatocellular carcinoma to sorafenib. By ranking the differentially expressed genes, FXYD3, FXYD1, TBX4, NPNT, and ADH1C were found to be the most significant differentially expressed genes. Among them, the clinical significance of the sodium-potassium ATPase regulator FXYD domain-containing transport regulator 3 (FXYD3) has been proved to have potential regulatory effects in a variety of cancers. In hepatocellular carcinoma, it can be used as an independent prognostic factor [29]. Studies have shown that FXYD3 is a multidrug resistance gene and plays an essential role in the multidrug resistance mechanism of hepatocellular carcinoma [30]. FXYD1 was differentially expressed in most tumors and overexpressed in hepatocellular carcinoma [31]. In hepatocellular carcinoma, TBX4 is methylated to regulate its expression and function [32].
Meanwhile, methylation of TBX4 may play an essential role in drug transport in 9 hepatocellular carcinoma [33]. The extracellular matrix protein nephronectin (NPNT) regulates cell adhesion, proliferation, and survival in the kidney. NPNT localizes in inflammatory foci and plays a role in the inflammatory response [34]. NPNT is imbalanced in the early stage of hepatocellular carcinoma and is mainly involved in the process of the cell cycle [35]. Alcohol dehydrogenase 1C (ADH1C) is located in the TSGS region of chromosome 4q. Its inactivation may be a late event after malignant transformation [36].
Besides, ADH1C affects the functional metabolism of hepatocytes and is also a target gene for the treatment of hepatocellular carcinoma with various Chinese herbal medicines.
Therefore, ADH1C may also be a target gene for drug resistance in sorafenib therapy [37,38].
On the other hand, these disordered genes are mainly involved in immune and inflammatory-related functions and signaling pathways. In the process of resistance to sorafenib in hepatocellular carcinoma (HCC), changes in immune and inflammatory responses are essential, as well as some epigenetic modifications. First of all, when RNA interference inhibited the phosphorylation of EGFR/her-3 and combined with sorafenib, the anti-proliferation and apoptotic ability of sorafenib was improved [39]. Apoptosis, autophagy, epithelial-mesenchymal transition (EMT) and other related pathways and reactions are also one of the mechanisms of sorafenib resistance [40]. Besides, the expression and phosphorylation of c-Jun in human hepatocellular carcinoma cells were enhanced after treatment with sorafenib. Inhibiting c-Jun promotes apoptosis induced by sorafenib [41]. Sphingolipids are often disturbed in sorafenib-resistant cancer cells.
Sphingosylphorylcholine can regulate various biological processes, such as growth, proliferation, migration, invasion, and metastasis [42][43][44][45], by controlling the signaling function in the signal transduction network of cancer cells.
NFKB1 is a transcription factor predicted to regulate disordered genes. NFKB1 gene can inhibit the expression of neutrophil chemokine network by expressing anti-inflammatory p50 nuclear factor-kappa B (NF-kappa B) dimer [46]. NFKB 1-94 INS promoter polymorphism increases the risk of HCC and can be used as a predictor of clinical stage and tumor size in women with HCC [47]. Finally, CDK1 is the most useful diagnostic factor in the diagnosis of drug resistance in hepatocellular carcinoma. Study had shown that CDK1 is a gene significantly associated with the prognosis of hepatocellular carcinoma, which has also been confirmed in other studies [48]. However, the clinical diagnostic value needs further verification.