Abstract

Hepatocellular carcinoma (HCC) represents a common malignancy, and mechanisms of acquired sorafenib resistance during the treatment of HCC patients remain elusive. The present study performed integrated bioinformatics analysis and explored the potential action of heme oxygenase 1 (HMOX1) in sorafenib-resistant HCC cells. Differentially expressed genes (DEGs) of the sorafenib-resistant group as compared to the sorafenib-sensitive group from GSE140202 and GSE143233 were extracted. Fifty common DEGs between GSE140202 and GSE143233 were extracted. Ten hub genes were identified from the protein-protein interaction network based on common DEGs. Experimental results revealed the upregulation of HMOX1 in sorafenib-resistant HCC cells. HMOX1 silence promoted the sensitivity to sorafenib in sorafenib-resistant HCC cells; overexpression of HMOX1 attenuated the sensitivity. In addition, HMOX1 silence downregulated the mRNA expression of ABC transporters in sorafenib-resistant HCC cells, while HMOX1 overexpression upregulated mRNA expression of ABC transporter expression in HCC cells. Further analysis also revealed that high expression of HMOX1 was associated with shorter OS and DSS in HCC patients. In conclusion, our analysis identified ten hub genes associated with sorafenib resistance in HCC. Further validation studies demonstrated that HMOX1 promoted sorafenib resistance of HCC cells via modulating ABC transporter expression.

1. Introduction

Hepatocellular carcinoma (HCC) represents a common tumor malignancy, and HCC incidence is expected to increase annually [1, 2]. Among all the cancer types, HCC is one of the most frequently diagnosed malignancies around the world [1, 2]. The HCC patients’ overall survival (OS) remains low due to insufficient early diagnosis and the recurrence and metastasis of advanced HCC [35]. In terms of early diagnosis of HCC, imaging has limitations in diagnostic accuracy and sensitivity, while common serum markers show poor diagnostic performance [6]. With the advent of high-throughput sequencing technology, genetic biomarkers such as cell-free DNA (cfDNA), epigenetic changes, and circulating RNA (microRNAs (miRNAs), long noncoding RNAs (lncRNA), and circular RNAs (circRNA)) from peripheral blood have become the focus of the early diagnosis of HCC [7]. However, there are still limitations in the early diagnosis of HCC. The main limitation of most studies using nucleic acid molecules as biomarkers of HCC is the limited size of the cohort, and the heterogeneity of HCC should be also considered. Sorafenib was first approved by the FDA for HCC therapy at advanced stages. Sorafenib treatment could significantly improve the OS of patients with HCC by 2-3 months. Unfortunately, many patients with HCC had a poor response to sorafenib or exhibited sorafenib resistance after prolonged use of sorafenib [810]. The molecular mechanisms underlying the acquired sorafenib resistance are complex which may involve epithelial-mesenchymal transition, tumor microenvironment, autophagy, and cancer stem cells [1012]. Besides, studies also proposed that acquired resistance to sorafenib may be associated with the dysregulated signaling pathways including JAK/STAT, PI3K/Akt, and TNFα/NF-κB [1012]. Thus, it is urgent to fully decipher mechanisms associated with sorafenib resistance, which may provide a better therapeutic strategy for treating advanced HCC.

With the rapid development and progress in the high-throughput technologies, the investigation of tumor biology has been focusing on the genomic scale. RNA sequencing has been widely used in identifying novel targets in cancer pathophysiology, including HCC. For instance, Weng et al. performed the RNA-sequencing analysis in HCC and sorafenib-resistant HCC tissues and identified novel cricFOXM1 as crucial modulator of sorafenib resistance in HCC cells [13]. Wu et al. also carried out global transcriptomic sequencing in sorafenib-resistant HCC cells and emphasized the importance of circRNAs in mediating sorafenib resistance [14]. Wu et al. performed the RNA-sequencing analysis and revealed that mitophagy promoted sorafenib resistance through hypoxia-inducible ATAD3A dependent axis [15]. RNA-sequencing studies also proposed that epigenetically activated ADAMTSL5 is a key player in HCC drug resistance [16]. In addition, the integrated bioinformatics analysis of publicly available microarray datasets also provided another strategy for scientists to discover novel mediators associated with the acquired sorafenib resistance in HCC. Liu et al. analyzed dataset GSE109211 and proposed that HCC sorafenib resistance was correlated with identified hub genes and pathways [17]. Jiang et al. analyzed GSE62813, GSE73571, GSE151412, and GSE140202 and found other key mediators in HCC sorafenib resistance [18].

This study employed a strategy to integrate two datasets (GSE140202 and GSE14322) from the GEO database and extracted the DEGs between sorafenib-sensitive and sorafenib-resistant groups. The extracted DEGs were processed for PPI network construction and functional analysis. Furthermore, HMOX1, one of the hub genes, was further validated in sorafenib-resistant HCC cells, and the potential action of HMOX1 in sorafenib-resistant HCC cells was also preliminarily explored.

2. Materials and Methods

2.1. Extraction of GEO RNA-Sequencing Datasets

RNA-sequencing datasets were retrieved from the public GEO database. In the analysis, we searched the RNA-sequencing datasets profiling the mRNA expression between sorafenib-sensitive and sorafenib-resistant HCC cells, and in each group were excluded in the analysis. After screening the processed GEO RNA-sequencing datasets using GREIN-iLINCS tool [19], two RNA-sequencing datasets, including GSE140202 and GSE143233, were finally included in our analysis. GPL20795 HiSeq X Ten and GPL16791 Illumina HiSeq 2500 platforms were used in GSE140202 and GSE143233, respectively. In GSE140202, 6 sorafenib-sensitive HCC cells () and sorafenib-resistant HCC cells () were included for analysis. In GSE143233, three samples of HCC tissues and three samples of sorafenib-resistant HCC tissues were included for analysis.

2.2. Extraction of DEGs in the Collected RNA-Sequencing Datasets

The extraction of DEGs in GSE140202 and GSE143233 was performed using the GREIN-iLINCS tool [19]. The DEGs between sorafenib-sensitive and sorafenib-resistant HCC cells were extracted in GSE140202; the DEGs between HCC and sorafenib-resistant HCC tissues were extracted for GSE143233. The significant DEGs were extracted using the following criteria: and values < 0.05. The heatmaps of DEGs in GSE140202 and GSE143233 were plotted by using the GREIN-iLINCS tool, and the top 200 DEGs were illustrated in the heatmaps. The R software plotted the volcano plots of the DEGs in GSE140202 and GSE143233 with ggplot function. The Venn diagrams showing the common DEGs between GSE140202 and GSE143233 were also generated using R software.

2.3. Functional Analysis of DEGs between GSE140202 and GSE143233

The GO functional enrichment and KEGG pathway enrichment analyses of common DEGs between GSE140202 and GSE143233 were carried out using EnrichR, an interactive and collaborative HTML5 gene list enrichment analysis tool.

2.4. PPI Network Construction Based on Common DEGs Extracted from GSE140202 and GSE143233

PPI network of the common DEGs between GSE140202 and GSE143233 was generated by STRING [20]. The PPI network was visualized by using Cytoscape software. To further extract the key submodules in the PPI network, the cytoHubba application from Cytoscape software was applied to visualize submodules derived from the PPI network.

2.5. Cell Lines and Cell Culture

The Huh7 and HepG2 cells were obtained from ATCC (Manassas, USA). These cell lines were kept in Dulbecco’s modified Eagle’s medium (DMEM; Sigma) supplemented with 10% fetal bovine serum (FBS: Gibco) at 37°C in a humidified incubator containing 5% CO2. The generation of sorafenib-resistant HCC cells was performed by treating cells with 0.5 μM sorafenib (Sigma) followed by escalating the dose to 10 μM. Sorafenib-resistant Huh7 (Huh7-SR) and HepG2 (HepG2-SR) cells were maintained in DMEM containing sorafenib (10 μM).

2.6. siRNAs and Plasmids

The siRNA that targets HMOX1 was obtained from RiboBio (Guangzhou, China). The scrambled sequence of HMOX1 siRNA was selected as the corresponding negative control (NC). The sequences for the corresponding siRNAs are shown in Supplementary Table S1. For HMOX1-overexpressing vectors, HMOX1 was ligated into the pcDNA3.1 vector to generate the HMOX1-expressing vector (RiboBio), and the empty pcDNA3.1 vector was used as NC.

2.7. Cell Transfections

Cells (Huh7, Huh7-SR, HepG2, and HepG2-SR cells) were seeded at cells per well in 6-well plates in DMEM containing 10% FBS overnight. Transfection experiments were performed with OPTI-MEM serum-free medium and Lipofectamine 2000 reagent with respective siRNA (scrambled siRNA or HMOX1 siRNA) or vectors (pcDNA3.1 or pcDNA3.1-HMOX1). At 24 h post-siRNA or vector transfection, these cells were harvested for other experimental procedures.

2.8. Cell Viability Assay

The cell viability of respective HCC cells was measured by the Cell Counting Kit-8 (CCK-8) assay (Beyotime, Beijing, China). Briefly, respective HCC cells after siRNA/vector transfections and/or treatment with incremental concentrations of sorafenib were incubated with 10% CCK-8 in DMEM for 4 h. The cell viability was measured by calculating absorbance at 450 nm with a plate reader.

2.9. Quantitative Real-Time PCR

Total RNA was isolated using Trizol reagent (Qiagen, Hilden, Germany). A NanoDrop spectrophotometer was applied to measure RNA concentration (Thermo Fisher Scientific). Five μg RNA was transcribed into cDNA with a reverse transcription kit (Applied Biosystems). StepOnePlus Real-Time PCR System was used to perform real-time PCR (Applied Biosystems). Thermal conditions of amplification on cycles were as follows: 95°C for 2 min, 40 cycles at 95°C for 30 s, 56°C for 30 s, and 72°C for 1 min, then 72°C for 5 min. GAPDH was used as an internal control. Relative expression of the corresponding mRNA was calculated using the 2-ΔΔCt method. The primer sequences are shown in Supplementary Table 2.

2.10. Survival Analysis of Patients with HCC

The impact of HMOX1 on OS and disease-specific survival (DSS) of HCC patients was evaluated by the Kaplan Meier plotter [21]. Log-rank was statistically significant.

2.11. Statistical Analysis

The statistical analysis was carried out by GraphPad Prism 5 software (GraphPad Software). All the data were displayed as . Unpaired Student’s -test determined significant differences between treatment groups. was statistically significant.

3. Results

3.1. Analysis of DEGs between Sorafenib-Sensitive and Sorafenib-Resistant Groups in GSE140202 and GSE143233 Datasets

The DEGs between two groups in GSE140202 and GSE143233 were plotted as heatmaps (Figures 1(a) and 1(b)) and volcano plots (Figures 1(c) and 1(d)). In GSE140202, 542 upregulated and 651 downregulated DEGs were extracted between two groups in GSE140202; in GSE143233, 372 upregulated DEGs and 430 downregulated DEGs were extracted between two groups in GSE143233. The Venn diagram showed that 22 upregulated common DEGs were found between GSE140202 and GSE143233 (Figure 2(a)); 26 downregulated common DEGs were found between GSE140202 and GSE143233 (Figure 2(b)).

3.2. Functional Analysis of Common DEGs

Common DEGs in GSE140202 and GSE143233 were further processed for functional enrichment analysis. In the GO_biological process, common DEGs were enriched in GO terms such as “protection from natural killer cell mediated cytotoxicity”; “antigen processing and presentation of endogenous peptide antigen via MHC class I via ER pathway”; “antigen processing and presentation of endogenous peptide antigen via MHC class I via ER pathway, TAP-independent”; “antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP-independent”; and “myelin assembly” (Figure 3(a)). In the GO_cellular component, common DEGs were enriched in GO terms such as “MHC class I protein complex,” “COPII-coated ER to Golgi transport vesicle,” “MHC protein complex,” and “lumenal side of endoplasmic reticulum membrane” (Figure 3(b)). In GO_molecular function, common DEGs were significantly enriched in the GO terms such as “metal ion binding,” “serine-type endopeptidase inhibitor activity,” “C3HC4-type RING finger domain binding,” and “apolipoprotein receptor binding” (Figure 3(c)). In KEGG pathway analysis, common DEGs were enriched in “Antigen processing and presentation,” “Allograft rejection,” “Natural killer cell mediated cytotoxicity,” and “Graft-versus-host disease” pathways (Figure 3(d)).

3.3. PPI Network Analysis from Common DEGs in GSE140202 and GSE143233

PPI network was established based on 48 common DEGs. The 12 nodes and 18 edges were detected in constructed PPI network (Figure 4(a)). Furthermore, the cytoHubba analysis extracted four submodules. The most significant submodule contains four nodes and four edges (Figure 4(b)).

3.4. HMOX1 Exhibited an Upregulation in the Sorafenib-Resistant HCC Cells

Based on the literature analysis, we proposed that HMOX1 may be a novel gene associated with sorafenib’s sensitivity to the HCC cells. Thus, we compared HMOX1 mRNA expression in the normal HCC cell lines (Huh7 and HepG2) and the sorafenib-resistant HCC cell lines (Huh7-SR and HepG2-SR). HMOX1 was significantly upregulated in the sorafenib-resistant HCC cell lines (Huh7-SR and HepG2-SR) compared to their corresponding parental HCC cells (Figures 5(a) and 5(b)). Furthermore, the MTT assay showed that the suppressive actions of sorafenib on the cell viability were significantly attenuated in Huh7-SR and HepG2-SR cells compared to the corresponding parental cells (Figures 5(c) and 5(d)).

3.5. Actions of HMOX1 on Sorafenib Sensitivity of HCC Cells

HMOX1 siRNA transfection significantly downregulated HMOX1 expression in Huh7-SR and HepG2-SR cells (Figure 6(a)). The MTT assay showed that silencing HMOX1 markedly promoted the sensitivity of Huh7-SR (Figure 6(b)) and HepG2-SR (Figure 6(c)) cells to sorafenib. Gain-of-function results revealed that pcDNA3.1-HMOX1 transfection induced an apparent increase in HMOX1 mRNA expression of Huh7 and HepG2 cells (Figure 6(d)). Functional assay results demonstrated that HMOX1 overexpression promoted sorafenib-resistance of Huh7 and HepG2 cells (Figures 6(e) and 6(f)).

To further explore mechanisms associated with HMOX1-mediated sorafenib resistance of HCC cells, we examined actions of HMOX1 overexpression or silence on ABC transporter mRNA expression in HCC cells. The ABC transporter mRNA expression levels of ABCA6, ABCB1, ABCC1, and ABCG2 in these cells were upregulated compared to their corresponding parental HCC cells (Figure 7). The loss-of-function results demonstrated that HMOX1 silence reduced expression levels of ABCA6, ABCB1, ABCC1, and ABCG2 in these cells (Figure 8). Consistently, the gain-of-function results revealed that HMOX1 overexpression elevated mRNA levels of ABCA6, ABCB1, ABCC1, and ABCG2 in Huh7 and HepG2 cells (Figure 9).

3.6. The Prognostic Role of HMOX1 Expression in HCC Patients

The effects of HMOX1 expression on the prognosis of HCC patients were determined by an online KM-plotter tool, and high expression of HMOX1 was associated with shorter OS and DSS in HCC patients (Figure 10).

4. Discussion

The acquired sorafenib resistance has significantly limited the therapeutic potential of sorafenib in advanced HCC, while mechanisms associated with acquired sorafenib resistance remain to be further clarified [10]. In this study, two RNA-sequencing datasets (GSE140202 and GSE143233) from GEO were downloaded for analysis. The DEGs between sorafenib-sensitive and the sorafenib-resistant group were extracted in these two datasets, and a total of 48 common DEGs between GSE140202 and GSE143233 were extracted. Functional enrichment analysis revealed that the common DEGs might be associated with “antigen processing and presentation,” “natural killer cell mediated cytotoxicity,” “cell adhesion molecules,” and so on. Ten hub genes were identified from the protein-protein interaction network based on common DEGs. Experimental results revealed the upregulation of HMOX1 in sorafenib-resistant HCC cells. HMOX1 silence promoted the sensitivity to sorafenib in sorafenib-resistant HCC cells; overexpression of HMOX1 attenuated the sensitivity. In addition, HMOX1 silence downregulated the mRNA expression of ABC transporters in sorafenib-resistant HCC cells, while HMOX1 overexpression upregulated mRNA expression of ABC transporter expression in HCC cells. Further analysis also revealed that high expression of HMOX1 was associated with shorter OS and DSS in HCC patients. Collectively, our results demonstrated that HMXO1 might be associated with sorafenib resistance in HCC.

Analysis of the microarray datasets has become a powerful tool for exploring novel genes that may be associated with cancer progression. In the GSE140202 dataset, TCNOS_00284048 and TCONS_00006019 were highly expressed in the sorafenib-resistant HCC cells compared with parental HCC cells. Knockdown of TCNOS_00284048 and TCONS_00006019 promoted sorafenib-sensitivity of Huh7-SR and HepG2-SR cells [22]. In GSE143233, Lin et al. demonstrated that METTL3 was underexpressed in human sorafenib-resistant HCC and revealed that RNA m6A methylation mediated sorafenib resistance via FOXO3-mediated autophagy [23]. Jiang et al. carried out integrated transcriptomic sequencing of GSE140202 and other datasets and identified 13 hub genes and seven promising therapeutic agents for HCC [18]. Li et al. also performed bioinformatics analysis for GSE140202 and found that GINS1 was highly expressed in sorafenib-resistant HCC cells [24]. In our results, we identified 50 common DEGs between GSE140202 and GSE143233, and these DEGs may be “antigen processing and presentation,” “natural killer cell mediated cytotoxicity,” “cell adhesion molecules,” and so on. Ten hub genes were detected from the PPI network. Among these hub genes, HMOX1 was highly expressed in sorafenib-resistant HCC cells as determined by qRT-PCR assay.

HMOX1 is well known for its enzymatic role in regulating cellular homeostasis under stress [25]. Increasing evidence has demonstrated the regulatory role of HMOX1 in cancer progression. For example, Yim et al. showed that HMOX1 is a prognostic marker for bladder cancer recurrence and progression [26]. Park et al. showed that the HMOX1/carbon monoxide axis inhibited transforming growth factor-β1-induced growth inhibition in HCC cells [27]. Inhibiting HMOX1 expression could retard HCC progression via distinct mechanisms [2830]. However, the action of HMXO1 in HCC sorafenib resistance has not been reported yet. Gao et al. showed that inhibition of HMXO1 could sensitize clear-cell renal cell carcinoma to sorafenib [31]. Our study showed that HMOX1 was upregulated in the sorafenib-resistant HCC cells. HMOX1 could promote the resistance of HCC cells to sorafenib. Multidrug resistance is closely regulated by ABC transporters [32]. Studies have demonstrated that sorafenib can interact with ABC transporters such as ABCB1, ABCC1, ABCG2, and ABCC10 [33, 34]. In our study, sorafenib-resistant HCC cells exhibited high expression of ABC transporters, including ABCA6, ABCB1, ABCC1, and ABCG2 HCC cells. In addition, HMOX1 silence downregulated the mRNA expression of ABCA6, ABCB1, ABCC1, and ABCG2 in sorafenib-resistant HCC cells, while HMOX1 overexpression upregulated these transporters’ expression in HCC cells. These results indicated that HMOX1-mediated sorafenib resistance might be associated with the modulation of the expression of ABC transporters.

5. Conclusions

In conclusion, the present study identified ten hub genes linked to the sorafenib resistance in HCC according to bioinformatics analysis. Further validation studies demonstrated that HMOX1 promoted the sorafenib resistance of HCC cells via modulating ABC transporter expression. However, the bioinformatic analysis was limited to two RNA-sequencing datasets, and future studies should explore more relevant datasets to identify more novel potential mediators in the regulating sorafenib sensitivity of HCC cells.

Data Availability

All the data are available from the corresponding author.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

CW and QL designed the overall study design. CW drafted the article. XZ and YZ performed the data analysis. WD participated in the figure preparation and revised the paper. YW, WD, and GD read literatures and assist in experiments. Xian Zhu and Yinfang Zhang equally contributed to this work.

Acknowledgments

This study was funded by the Key Science and Technology Project of Shanghai Songjiang Health Bureau (2020SJ293).

Supplementary Materials

Supplementary 1. Supplementary Table 1: siRNA sequences used in the study.

Supplementary 2. Supplementary Table 2: primer sequences for qRT-PCR analysis.