Prognostic comparative genes predict targets for sorafenib combination therapies in hepatocellular carcinoma.

With the increasing incidence and mortality of human hepatocellular carcinoma (HCC) worldwide, revealing innovative targets to improve therapeutic strategies is crucial for prolonging the lives of patients. To identify innovative targets, we conducted a comprehensive comparative transcriptome analysis of 5,410 human HCCs and 974 mouse liver cancers to identify concordantly expressed genes associated with patient survival. Among the 664 identified prognostic comparative HCC (pcHCC) genes, upregulated pcHCC genes were associated with prognostic clinical features, including large tumor size, vascular invasion and late HCC stages. Interestingly, after validating HCC patient prognoses in multiple independent datasets, we matched the 664 aberrant pcHCC genes with the sorafenib-altered genes in TCGA_LIHC patients and found these 664 pcHCC genes were enriched in sorafenib-related functions, such as downregulated xenobiotic and lipid metabolism and upregulated cell proliferation. Therapeutic agents targeting aberrant pcHCC genes presented divergent molecular mechanisms, including suppression of sorafenib-unrelated oncogenic pathways, induction of sorafenib-unrelated ferroptosis, and modulation of sorafenib transportation and metabolism, to potentiate sorafenib therapeutic effects in HCC combination therapy. Moreover, the pcHCC genes NCAPG and CENPW, which have not been targeted in combination with sorafenib treatment, were knocked down and combined with sorafenib treatment, which reduced HCC cell viability based on disruption to the p38/STAT3 axis, thereby hypersensitizing HCC cells. Together, our results provide important resources and reveal that 664 pcHCC genes represent innovative targets suitable for developing therapeutic strategies in combination with sorafenib based on the divergent synergistic mechanisms for HCC tumor suppression.

Introduction

Hepatocellular carcinoma (HCC) is the main subtype of primary liver cancer and the sixth most common cancer and the fourth leading cause of death of cancer patients worldwide [1]. Most HCC patients are found in developing countries in Asia and Sub-Saharan Africa. However, the incidence of HCC in developed countries, such as the USA, Japan and some European countries, is rising because of increasing cases of hepatitis C virus infection and nonalcoholic steatohepatitis (NASH) associated with metabolic syndrome or diabetes mellitus [2]. Hence, the World Health Organization estimates that more than 1 million patients will die from liver cancers in 2030 [3]. Depending on the tumor burden and liver function of HCC patients, treatments with hepatic surgical resection, local ablation and liver transplantation are the major curative interventions for early HCC patients[4]. Nevertheless, with the vast majority of HCC patients being first diagnosed at the noncurable advanced stage, patients will receive systemic therapies with antiangiogenic multiple kinase inhibitors and/or metronomic chemotherapeutic drugs as well as immune checkpoint inhibitors [4], [5], [6], [7]. Sorafenib is the standard FDA-approved first-line therapeutic drug for the treatment of advanced-stage HCC, and it has a response rate of<10% and a therapeutic benefit of prolonging life by a median of 3 months for HCC patients [8], [9]. Although therapeutic response rates exceeded 20% in HCC patients after treatment with the newly developed antiangiogenic multiple kinase inhibitors (e.g., lenvatinib) and/or immune checkpoint inhibitors (e.g., nivolumab and atezolizumab) [10], [11], [12], applying innovative genomic approaches to identify additional biomarkers for guiding suitable therapies and targets for developing novel therapeutic strategies remains critical for prolonging survival for the majority of HCC patients [4], [13].

After applying comprehensive sequencing approaches to the HCC genomes, only a limited number of somatic mutations were revealed as potential therapeutic targets; moreover, frequently observed mutations are not currently applied for predicting therapeutic responses in clinical practice [14], [15]. The genetic heterogeneity of HCC genomes could be due to the accumulation of somatic mutations during long-term exposure to divergent etiological risk factors, such as viral hepatitis, alcohol abuse, aflatoxin ingestion and metabolic syndrome. In particular, intratumoral heterogeneity with frequent HCC genomic alterations has been shown to confer selective advantages of cell proliferation, and it is commonly associated with tumor recurrence, drug resistance and metastasis, thus leading to poor outcomes in HCC patients [16], [17], [18].

Comparative and integrated cancer genomic approaches are commonly applied to identify cross-species aberrant cancer genes and in best-fit mouse models of cancer subtypes of several cancers to reveal the tumorigenic mechanisms and develop therapeutic approaches to improve therapeutic interventions [19], [20], [21], [22], [23], [24]. Notably, limited efforts were applied to reveal the overall consensus of pathophysiological functions shared in human and mouse cancer genomes, reveal the common tumorigenic pathways and predict innovative targets for dissecting the molecular mechanisms underlying combination therapies for therapeutic improvement [22], [25].

In this study, we applied the most comprehensive comparative transcriptome analysis to date (to the best of our knowledge) that included large-scale liver cancer transcriptome datasets including 5410 human HCC samples and 974 mouse liver cancers to reveal the major concordant pathophysiological functions in HCC progression. As enriched HCC differentially expressed genes (DEGs) are commonly associated with drug metabolism and cell proliferation [26], we hypothesized that comparative HCC genes associated with patient prognosis might include a large number of known and unknown theranostic targets suitable for sorafenib combination therapies for developing innovative strategies to prolong the lives of HCC patients. Indeed, we identified 664 prognostic comparative HCC (pcHCC) genes, ∼44% of which have not been thoroughly studied in HCC, uncovered targets for sorafenib combination therapies with divergent molecular mechanisms to synergize HCC tumor suppression, and revealed an innovative oncogenic p38/STAT3 axis that enhances sorafenib tumor suppression for potential improvement of therapeutic interventions in HCC patients.

Results

Identification of 664 pcHCC genes associated with patient prognosis

We collected large-scale transcriptome datasets of human HCCs and mouse liver cancers from public databases, including 5410 human HCCs (52 datasets from 26 array platforms, Supplementary Table 1) and 974 mouse liver tumors (52 datasets from 23 array platforms, Supplementary Table 2), for comparative analysis. All human transcriptomes were derived from HCC tumors, with the exception of the cell line transcriptomes. The majority of mouse transcriptomes were annotated as HCC (∼70%) or liver cancer (∼30%) in genetically modified models with and without carcinogen induction (Supplementary Table 2 for details). After inferring the gene expression intensity based on the log2 tumor/normal ratio and Z score normalization, we generated expression profiles of human HCCs that consisted of 3,339 tumor/normal ratios of 26,182 differentially expressed genes (DEGs) and expression profiles of mouse liver cancer that consisted of 633 tumor/normal ratios of 18,157 DEGs (analysis pipeline in graphic abstract). We then determined the concordant aberrantly expressed comparative DEGs and revealed that 1,904 (766 upregulated and 1,138 downregulated) DEGs participated in both human HCC and mouse liver cancer. To emphasize the clinical value, we matched the 1,904 comparative HCC DEGs with prognostic HCC DEGs in the TCGA_LIHC dataset and identified 664 pcHCC genes for experimental and translational evaluation (Fig. 1A and Supplementary Table 3). Notably, we found that few HCC studies have focused on 293 out of 664 (∼44%) pcHCC genes (labeled “unknown” as HCC genes) based on PubMed searches.

Fig. 1

Distribution and validation of 664 pcHCC genes based on the TCGA_LIHC. (A) Aberrantly expressed prognostic comparative HCC genes are shown based on their proportion (% of HCC samples with upregulation less the % of HCC samples with downregulation) of aberrant expression in overall HCC samples (y-axis) and the aberrant portion of HCC samples (x-axis); and (B) validation and distribution of pcHCC genes with the TCGA_LIHC dataset. The pcHCC genes were divided into known and unknown HCC genes based on literature searches in PubMed. The top 10 genes in each of the four categories are listed between the figure panels and marked with colored dots, with known upregulated, unknown upregulated, known downregulated, and unknown downregulated pcHCC genes shown in orange, red, sky blue and blue, respectively.

Validation of the aberrant known and unknown pcHCC genes

To validate aberrantly expressed pcHCC genes in HCC, we selected two independent HCC cohorts from two different platforms, namely, iCOD (integrated clinical omics database with transcriptomes conducted using Affymetrix HG-U133 plus 2.0 arrays) [27] and TCGA_LIHC (transcriptomes conducted by next-generation sequencing technology, RNA-Seq) [14], to confirm the altered gene expression in HCCs. Consistently, our results showed concordant aberrant pcHCC gene expression between two selected transcriptomic datasets (pcHCC genes from the meta-analysis results of this study, iCOD and TCGA_LIHC) with a very significant positive correlation (Pearson’s correlation coefficient r > 0.61, P value < 2.2 × 10−16) (Fig. 1B and Supplementary Fig. 1A and B, and Supplementary Table 4).

To detail the involvement of the 664 pcHCC genes in HCC tumor progression, we searched PubMed and found that a total of 354 upregulated (53.3% of 664) pcHCC genes, including 219 (62%, orange dots) that were previously reported to be related to HCC and 135 (38%, red dots) that were not previously associated with HCC (Fig. 1A and Supplementary Table 3). The majority of the top 10 known upregulated pcHCC genes (orange dots), such as TOP2A, PRC1, and AKR1B10, were studied thoroughly, with multiple reports focusing on their roles in HCC tumor progression and poor prognosis [28], [29], [30]. In contrast, the top 10 unknown pcHCC genes (red dots) were rarely studied or even had zero HCC results (e.g., RRAGD) in PubMed searches. For the 310 (46.7% of 664) downregulated genes, 152 (49%, sky blue dots) were known to be downregulated, such as TAT, CYP2C9 and SLCO1B, and they were associated with better prognosis of HCC patients in some reports [31], [32], [33]. On the other hand, 158 unknown and downregulated pcHCC genes (51%, blue dots) have been seldom studied or even yielded no search results with HCC, such as CFHR3 and TTC36. Moreover, other known HCC genes that were upregulated (e.g., GPC3, AURKA and others) and downregulated (e.g., CLEC1B, CRHBP and others) were also identified in the majority of HCC samples (as percentages on the x-axis), with differential expression intensities between tumor and normal HCC samples (as percentages on the y-axis) (Fig. 1A and Supplementary Table 3). Interestingly, several pcHCC genes were found to have very little or no representation among HCC studies, including CENPW and CENPU (upregulated) as well as OIT3, IDO2, CNDP1 and ANGPTL6 (downregulated), according to PubMed searches. Together, our results provide evidence that the 664 pcHCC genes not only participate in HCC tumorigenesis but are also useful for the development of innovative theranostic biomarkers for therapeutic interventions of HCC.

Validation of pcHCC genes associated with clinical features and patient survival

To further confirm the aberrant expression of the 664 pcHCC genes involved in altered liver function during tumor progression, we examined their aberrant expression associated with the clinical and demographic characteristics of HCC patients from the iCOD dataset (Supplementary Fig. 2). Consistent with previous reports, aberrant expression of pcHCC genes in HCCs was associated with known independent predictors of patient prognosis, including vein invasion, tumor size, diabetes and late tumor stage [34], [35], [36]. Then, we conducted unsupervised clustering to divide the genes into two groups, namely, upregulated and downregulated pcHCC genes. As expected, according to the Kaplan–Meier plot, the cluster of HCC patients with high combined expression of upregulated pcHCC genes in the TCGA-LIHC dataset exhibited worse survival than the cluster of HCC patients with low combined expression, with the comparison of Cox coefficients showing a significant P value (log rank test, P < 0.0001) (Fig. 2A). In contrast, the cluster of HCC patients with high combined expression of downregulated pcHCC genes showed a better prognosis than the cluster of HCC patients with low combined expression in the TCGA-LIHC dataset, with the comparison of Cox coefficients showing a significant P value (log rank test, P = 0.0021) (Fig. 2B).

Fig. 2

Validation of prognostic comparative HCC genes with Kaplan–Meier plots of HCC patient survival based on the TCGA-LIHC dataset after the Cox proportional hazards regression survival analysis. (A, B) Kaplan–Meier plots of (A) upregulated and (B) downregulated pcHCC genes, with significant Cox coefficient P values separated by the combined expression levels using unsupervised clustering; (C-E) Kaplan–Meier survival curves of selected pcHCC genes that are (C) upregulated known HCC genes; (D) upregulated (upper panel) and downregulated (bottom panel) unknown HCC genes; and (E) downregulated known pcHCC genes based on patient prognosis information from the TCGA-LIHC dataset.

We also validated the aberrant expression of pcHCC genes individually in association with the prognosis of HCC patients in the TCGA-LIHC dataset, including known (NCAPG, ASPM, TOP2A, PRC1, CCNB2, AKR1B10, KIF20A, UBE2T, CDC20 and PBK) and unknown (CENPW, HMMR, CKAP2 L, TRIP13 and GINS1) upregulated pcHCC genes associated with poor survival in Kaplan–Meier plots (Fig. 2C and D). In contrast, known (SOCS2, C6, CYP8B1, TAT, LCAT, CYP2C9, CYP2C8, ADH4, SLC22A1 and NAT2) and unknown (HGFAC, TTC36, DNASE1L3, AKR1D1 and KLKB1) downregulated pcHCC genes were associated with better survival of HCC patients (Fig. 2D and E). Collectively, we concluded that the 664 pcHCC genes are important resources as theranostic biomarkers and targets; thus, their pathological roles in HCC tumor progression should be further investigated and associated therapeutic interventions should be developed to prolong the life of HCC patients.

Gene set enrichment analysis (GSEA) of 664 pcHCC genes

To reveal the core pathophysiological functions of the 664 pcHCC genes in HCC progression, we compared the functional hallmarks of the GSEA signatures of these genes with two HCC DEGs, including the HCC DEGs from the transcriptomes of human and mouse microarrays and from the TCGA-LIHC datasets. Based on the significance of the functional gene signatures, we found that genes associated with cell proliferation-related hallmarks, including genes for E2F targets, MYC targets, mitotic spindle, and G2M checkpoint, were upregulated, and genes associated with metabolism-related hallmarks, including genes for the metabolism of xenobiotics, bile acids, coagulation, fatty acids, and adipogenesis, were downregulated (Fig. 3A and Supplementary Table 5). To further validate whether aberrant hallmarks of cell proliferation and metabolism play critical roles in the prognosis of HCC patients, we collected cell proliferation-related pcHCC genes (n = 116, Supplementary Table 6) and metabolism-related pcHCC genes (n = 88, Supplementary Table 7) to understand their relationship with the prognosis of HCC patients. Indeed, our results demonstrated that upregulation of cell proliferation-related pcHCC genes and metabolism-related genes was significantly associated with poor (log rank test, P = 0.004) and better (log rank test, P < 0.001) HCC prognosis, respectively (Fig. 3B). Together, our results showed that the core pathological features in the upregulation of cell proliferation and downregulation of certain metabolic pathways, especially lipid and xenobiotic (drug) metabolism-related pathways, presented significant involvement in tumor malignancy in human HCC and mouse liver cancer models.

Fig. 3

Gene set enrichment analysis (GSEA) of pcHCC genes. (A) Heatmap of the normalized enrichment score (NES) compared with GSEAs of the pcHCC genes, differentially expressed genes (DEGs) in all liver cancer microarrays, and DEGs in TCGA_LIHC. Significance was defined at P < 0.05. (B) Kaplan–Meier plots of downregulated metabolism-related (n = 89) and upregulated cell proliferation-related (n = 116) pcHCC genes.

Predicting targets for sorafenib combination therapies in pcHCC genes

Based on the 664 pcHCC genes exhibiting aberrant proliferation as known targets of the HCC first-line therapy drug sorafenib [26] and xenobiotic pathways known to have roles in drug metabolism, we hypothesized that pcHCC genes might involve multiple pathways that are correlated with sorafenib-altered HCC expression profiles in patients who received sorafenib therapy and thus may be useful for identifying targets and therapeutic agents for sorafenib combination therapy. In the TCGA-LIHC dataset, 28 patients who received sorafenib treatment and had RNA-Seq transcriptome data were downloaded for the analysis. After matching the sorafenib-altered HCC genes from the TCGA_LIHC dataset with the 664 pcHCC genes, we divided these genes into four possible quadrants, namely, Q1: up/up (n = 205), Q2: down/up (n = 116), Q3: down/down (n = 194), and Q4: up/down (n = 149), which represented the alterations in gene expression relative to the pcHCC genes/sorafenib-altered HCC genes, and the values were calculated as log2 ratios (Fig. 4A and Supplementary Tables 8 and 9). For instance, the representative top 10 known and unknown pcHCC genes were redistributed into four quadrants, as shown with the red lines/rectangles, for scenarios with potential further upregulation (Q1: up/up) and downregulation (Q3: down/down) after sorafenib treatment (Fig. 5A). In contrast, the green lines/rectangles indicate scenarios with opposite changes in expression (Q2: down/up; Q4: up/down) of the HCC genes after sorafenib treatment.

Fig. 4

PcHCC genes predict potential targets for sorafenib combination therapies. (A) Redistribution of the 664 pcHCC genes after matching with DEGs of HCC from patients who received sorafenib treatments in TCGA_LIHC into four quadrants. Altered expression that could be further up- or downregulated is shown with red lines/rectangles, and changes in aberrant expression from up to down or down to up are shown with green lines/rectangles. (B) Selected redistributed genes in each of the four quadrants and the top 5 enriched gene signatures.

Fig. 5

Functional and biochemical examinations of CENPW and NCAPG in HCC cell lines. (A, B) NCAPG (A) and CENPW (B) were upregulated in the local HCC tumors and corresponding normal pairs (n = 21) based on RT–qPCR. (C–E) NCAPG knockdown by two shRNAs reduced cell proliferation and migration in the (C) Mahlavu and (D) HCC36 HCC cell lines, as shown by bar graphs and representative migrated cells after crystal violet staining. CENPW knockdown by two shRNAs reduced cell proliferation and migration in the (E) Mahlavu and (F) HCC36 HCC cell lines, as shown by bar graphs and representative migrated cells after crystal violet staining. (G, H) Western blotting analysis of the representative markers MCM2 and PCNA for cell proliferation and MAPK signaling for the ERK, JNK and p38, and STAT pathways after shRNA knockdown of (G) NCAPG and (H) CENPW in Mahlavu cells. Error bars are the mean ± s.d. **P < 0.001 and ***P < 0.0001, determined by two-tailed Student’s t test (95% confidence interval).

To further validate whether these sorafenib-altered pcHCC genes could be useful as targets in combination therapies, we conducted a GSEA and found that these genes participated in several signaling pathways, including the PI3K/AKT/mTOR, xenobiotic metabolism, angiogenesis, hypoxia and ferroptosis pathways, which were previously reported as sorafenib-related pathways and potential targets for sorafenib combination therapies (Fig. 4B and Supplementary Table 10) [56]. For instance, AKR1B10 (aldo–keto reductase family 1 member B10) in Q1 is known to convert carbonyl compounds on aldehydes or ketones to alcohols for detoxification in the liver, and its upregulation is associated with the prognosis of HCC patients [28]. Because AKR1B10 is upregulated in HCC after sorafenib treatment, combination therapy with an AKR1B10 inhibitor via suppression of the mTOR signaling pathway enhanced the inhibitory effect of sorafenib on HCC xenografts in a nude mouse model (Fig. 4B and Table 1) [43]. In addition, SLC27A5 (in Q3), PCK1 (in Q2), SLC22A1 (in Q2) and SLCO1B3 (OATP1B3 in Q3) are known sorafenib transporters and/or participate in sorafenib metabolism in the liver [45], [46], [50], [57], and altered or aberrant expression of these transporters might modulate the therapeutic efficacy of sorafenib and/or activate the TXNRD1 and NRF2/KEAP1 pathways to increase oxidative stress and tumorigenicity in HCC. Therefore, combination therapy with sorafenib and NRF2/TXNRD1 inhibitors could be a promising strategy to improve HCC therapy [45], [46], [50], [57]. Other known sorafenib-altered pcHCC genes, such as SLC7A11 (in Q4), are key modulators of ferroptosis, which has led to the application of another ferroptosis inducer, artesunate, in sorafenib combination therapy for HCC [58]. AURKB (Aurora kinase B, in Q4) is a known HCC oncogene, and its inhibitor was combined with sorafenib to further induce HCC tumor shrinkage [54] (Fig. 4B). Collectively, the results showing the different mechanisms to support sorafenib combination therapies (Table 1) provide a better understanding and rationale for exploring innovative mechanisms of sorafenib-altered pcHCC genes to further develop strategies for the improvement of sorafenib combination therapy.

Table 1

List of known targets and modalities of the sorafenib combination HCC therapies.

Target genes	Quadrants	Agent for sorafenib combination therapies	Mechanism of sorafenibHCC combination therapy	Ref
HK2	Q1	Dichloroacetate	Sensitizes sorafenib-res. cells	[37]
CD24	Q1	Gedatolisib (PKI-587)	Inhibits CD24+ cells	[38]
PCNA	Q1	TLR3 synergist (BM-06)	Suppresses PCNA expression	[39]
DNMT1	Q1	AuNPs-anti-miR221	Inactivates DNMT1 signaling	[40]
E2F1	Q1	S-1	Downregulates E2F1	[41]
LOX	Q1	β‑aminopropionitrile	Diminishes angiogenesis	[42]
AKR1B10	Q1	Epalrestat	Enhances sorafenib effect	[43]
FOXO1	Q3	Anti-miRNA27a	Upregulates FOXO1-related apoptosis	[44]
SLC27A5	Q3	Brusatol	Sensitizes cells to sorafenib	[45]
SLCO1B3	Q3	SLCO1B3 deficient	Reduces sorafenib clearance	[46]
AR	Q2	PT-2385	Increases AR to inhibit growth	[47]
TXNIP	Q2	MEAN	Increases TXNIP to inhibit tumor	[48]
HGF	Q2	Vitamin K	Reduces HGF-stimulated growth	[49]
PCK1	Q2	Auranofin	Sensitizes sorafenib-apoptosis	[50]
HTATIP2	Q4	Metformin	Upregulates HTATIP2/downregulates Thioredoxin	[51], [52]
CCNB1	Q4	MBRI-001	Downregulates CCNB1	[53]
TXNRD1	Q4	Brusatol	Downregulates TXNRD1/up KEAP1/NRF2	[45]
AURKB	Q4	PHA-739358	Suppresses AURKB	[54]
CDK4	Q4	Palbociclib	Suppresses CDK4/6	[55]

Combination treatment with sorafenib and shRNAs of NCAPG or CENPW synergized suppression of cell viability by diminishing p38/STAT3 signaling

To explore the potential of innovative targets among these pcHCC genes for sorafenib combination therapies in HCC, we investigated the mechanistic roles and testing therapeutic efficacy of two pcHCC genes, NCAPG (non-SMC condensin I complex subunit G, Fig. 3C) and CENPW (centromere protein W, Fig. 3D), which are involved in HCC cell proliferation but have not been studied in sorafenib combination therapy. In addition to the association of these genes with poor prognosis in HCC patients, aberrant upregulation of NCAPG and CENPW was significantly associated with extrahepatic metastasis (BCLC stage C) (Supplementary Fig. 3A and C) and tumor size but not with venous invasion (Supplementary Fig. 3B and D) in HCC patients. With these supportive results, we continued to evaluate the suppressive effects of knocking down the expression of upregulated NCAPG and CENPW with specific short hairpin RNAs (shRNAs) and then investigated the impact on downstream cell proliferation signaling in HCC cells to validate the sorafenib combination therapy (Fig. 5).

First, we confirmed the upregulated expression of NCAPG and CENPW in locally recruited HCC tumors compared to their adjacent normal tissues (n = 21 pairs) through RT–qPCR experiments (Fig. 5A and B) and knocked down NCAPG and CENPW with their corresponding shRNAs to evaluate their roles in the predicted functions of HCC cell proliferation and migration. As expected, knockdown of NCAPG (shNCAPG #2 and #3) and CENPW (shCENPW #3 and #4) with two independent shRNAs in two HCC cell lines, Mahlavu and HCC36, significantly diminished cell proliferation and decreased cell migration, as shown by transwell migration assays (Fig. 5C–F). Our results suggested that upregulated expression of NCAPG and CENPW plays important roles in HCC tumor growth and that cell migration leads to tumor progression.

Second, to better understand the molecular pathways involved in NCAPG- and CENPW-driven HCC tumor progression, we examined the protein expression in the MAPK and STAT3 pathways, which are known for their involvement in HCC cell proliferation, survival and migration, through Western blotting analysis [59], [60]. Our results not only confirmed the knockdown efficiency of shRNAs against CENPW and NCAPG but also validated the reduction in cell proliferation by decreasing the expression of the cell proliferation markers PCNA and MCM2 (Fig. 5G and 5H). Interestingly, knockdown of NCAPG and CENPW with shRNAs consistently diminished the expression of p-STAT3 (pSer727) and MAPK p-p38 (pThr180/pTyr182) signaling but not ERK/p-ERK (pThr202/Tyr204) and JNK/p-JNK (pThr183/Tyr185) MAPK signaling, as shown by the Western blotting analysis (Fig. 5G and H). Previous studies indicated that the p38/STAT3 axis (but not the ERK/JNK MAPK signaling pathway) is involved in the expression of proinflammatory cytokines (e.g., IL6) and the stress-induced acute phase response of HCC cells [61], [62], [63]. Nevertheless, the biological functions of NCAPG and CENPW in the p38/STAT3 axis in HCC might be different because knockdown of NCAPG diminished the expression of p-p38, p-STAT3 and STAT3 but not the overall protein level of p38, and knockdown of CENPW reduced the expression of p-p38, p38, and p-STAT3 but not the total protein level of STAT3.

Finally, we conducted combination treatments of sorafenib with shRNAs of NCAPG or CENPW to identify the viability of HCC cells. Our results demonstrated that knockdown of either NCAPG or CENPW in Mahlavu and HCC36 cells conferred significant hypersensitivity to sorafenib treatment in a dose-dependent manner (Fig. 6A–D). Moreover, sorafenib treatment of shRNA-knocked down NCAPG or CENPW Mahlavu cells further decreased the expression of p-p38 and p-STAT3 in a dose-dependent manner, as shown by the Western blotting analysis (Fig. 6E and F). Together, our results suggest that NCAPG and CENPW, which suppress p38/STAT3 signaling, are potential targets for sorafenib combination therapies in HCC. In addition, our results suggest that the 664 pcHCC genes might contain innovative targets and pathways for developing therapeutic agents for HCC sorafenib combination therapies.

Fig. 6

Combination treatments of sorafenib with shRNAs of NCAPG- or CENPW-synergized suppression of cell viability by diminishing p38/STAT3 signaling. (A–D) NCAPG and (A, B) CENPW (C, D) knockdown with the combined treatments at various doses of sorafenib synergized with a decrease in the viability of the HCC cell lines Mahlavu (A, C) and HCC36 (B, D) in a dose-dependent manner. (E, F) Western blotting analysis of the p38 and STAT pathways after shRNA knockdown of (E) NCAPG or (F) CENPW with combination treatments of sorafenib at various doses in Mahlavu cells. Error bars are the mean ± s.d. **P < 0.001 and ***P < 0.0001, determined by two-tailed Student’s t test (95% confidence interval).

Discussion

In this study, we conducted a comprehensive comparative genomic analysis of large-scale transcriptome datasets (5410 human HCC and 974 mouse liver cancer) to identify concordantly expressed HCC genes associated with altered clinical features in HCC progression and patient prognosis (pcHCC genes) and reveal innovative targets for sorafenib combination therapy. Although HCC transcriptome analyses have been conducted in recent years with a focus on altered cell proliferation and/or metabolic reprogramming with different approaches [18], [64], [65], [66], to our knowledge, our results are the first to indicate that the enriched signatures harbored innovative targets (∼44%) suitable for developing therapeutic agents/inhibitors in sorafenib combination therapy. We provided lines of evidence from in silico analyses (Fig. 4 and Table 1) and experimental studies of the selected pcHCC genes NCAPG and CENPW (Fig. 5, Fig. 6) to support our hypothesis. Experimentally, we revealed that the proinflammatory and stress-related p38/STAT3 axis, but not conventional ERK/JNK MAPK signaling, plays critical roles in HCC progression and is suitable as a target for sorafenib combination therapy. Moreover, we identified 664 pcHCC genes harboring numerous innovative HCC targets (∼44%) for developing HCC biomarkers associated with patient prognosis, thus revealing innovative signaling pathways in tumor progression and identifying theranostic targets to predict drugs/agents for sorafenib combination therapies to potentially prolong the life of HCC patients.

With the increasing incidence and mortality of HCC worldwide, rapid developments have been made in systemic therapies for advanced HCC and several therapeutic options have been FDA-approved as first-line therapies, including sorafenib [8], [9], lenvatinib [10], and atezolizumab plus bevacizumab [12]. Other antiangiogenic multikinase inhibitors, such as regorafenib and cabozantinib, as well as the antibody drug ramucirumab, which targets anti-vascular endothelial growth factor-2 (VEGF-R2), and immune checkpoint inhibitors (ICIs) against PD-1 and cytotoxic T lymphocyte antigen 4 (CTLA-4), were approved for the second-line treatment of HCC [11], [67], [68]. Despite these therapeutic advances, only approximately 25% of HCC patients show a therapeutic response and improved overall survival, and the current therapeutic regimens remain ineffective for prolonging the life of most HCC patients [69]. Although sorafenib therapy in HCC offers some benefits, including a 2.8-month longer median overall survival than the placebo, endeavors to improve HCC sorafenib therapy, including the discovery of therapeutic biomarkers and combination therapy with other drugs, remain critical due to the limited availability, affordability and effectiveness of sorafenib as a single-drug HCC therapy.

For the discovery of therapeutic biomarkers of sorafenib, since frequent substitutional driver mutations with therapeutic options have not been identified after large-scale sequencing efforts in human HCC genomes [14], copy number alterations, especially gene amplification, leading to increasing expression of target genes have attracted specific attention to select subgroups of HCC patients for targeted sorafenib therapy. Indeed, the FGF3/FGF4 amplicon at chromosomal region 11q13 [70] and the VEGFA amplicon at chromosomal region 6p21.1 [71] were amplified in ∼2% and 7 ∼ 10% of human HCC genomes, respectively, and they encode the potential oncogenic driver proteins FGF4 and VEGF-A, which showed better responses to sorafenib treatment. Although copy number alterations in the HCC cancer genome commonly affect a small population of cancer patients, the identification of therapeutic targets/biomarkers could stratify subgroups of HCC patients to develop more precise HCC therapies. Notably, we also explored the correlations of the 664 aberrantly expressed pcHCC genes with copy number alterations (CNAs) in the HCCs of the TCGA_LIHC dataset (Supplementary Fig. 4). The lack of statistical associations between CNAs and aberrant gene expression of the 664 pcHCC genes (r = 0.0058, P value = 0.1373) implied the low potential for developing pcHCC genes as CNA-based therapeutic biomarkers in HCC.

Although some comparative genomics approaches have been applied for decades to develop best-fit mouse models of cancer subtypes and identify driver genes in cancers, including HCC [21], [23], [72], [73], a systematic and detailed explorations of comparative HCC genes implicated in the core HCC pathophysiological and clinical features and patient prognosis and even the targets of sorafenib combination therapy are lacking. Based on the aberrant enrichment of upregulated cell proliferation and downregulated lipid and xenobiotic metabolism-related genes of the 664 pcHCC genes, we hypothesized and systematically explored the underlying molecular mechanisms of sorafenib therapies in HCC. In addition to matching the divergently acting molecular mechanisms of sorafenib combination therapies reported previously (Table 1), we also selected two pcHCC genes, NCAPG and CENPW, because the aberrant upregulation of these genes is associated with the clinical features and prognosis of HCC patients (Supplementary Fig. 3A–D); moreover, these genes have not been previously targeted in sorafenib combination therapies to explore the oncogenic mechanisms and develop therapeutic strategies. NCAPG encodes a non-SMC (structural maintenance of chromosome) condensin I complex subunit G protein that has been shown to be upregulated in HCC to activate the PI3K/AKT pathway and enhance cell proliferation, and it has been associated with poor prognosis in HCC patients [74], [75], [76], [77], [78], [79]. Upregulated CENPW, the centromere protein W, which is sometimes also called cancer upregulated gene 2 (CUG2), was found to enhance cancer epithelial to mesenchymal transition (EMT) and stemness via activation of the canonical TGF-β signaling pathway and tumor progression in other cancers [80], [81]; however, its role in HCC is unknown. Neither the upregulation of NCAPG nor CENPW or related inhibitors has been shown to play a role in sorafenib combination therapy in HCC. Interestingly, we found that a rarely studied p38/STAT3 proinflammatory/stress-related axis (not JNK/ERK signaling) was required to sustain HCC oncogenic features. The combination of shRNAs of NCAPG or CENPW with sorafenib in the treatment of HCC cells showed synergistic suppression of oncogenic features and downstream signaling, further suggesting that pcHCC genes are an important resource as theranostic targets for early biomarker detection and therapeutic agents/inhibitors for sorafenib combination therapy to improve HCC therapy.

Because of the increasing need for cancer therapy worldwide, including HCC therapy, and the superior therapeutic benefits of combination therapy [12], [82], [83], [84], systematic exploration of sorafenib action mechanisms and innovative targets for combination therapy is crucial for improving HCC therapy. Our systematic prediction of potential targets among 664 pcHCC genes for the sorafenib combination therapy provides rationales and resources for further improving HCC therapy with the following features. First, we provide a comprehensive landscape of predictive sorafenib combination therapies with divergent mechanisms based on different scenarios of altered HCC gene expression (Fig. 4 and Table 1), and the findings were validated using selected target genes (Fig. 6). Second, the innovative pcHCC targets identified in our study could also be biomarkers for predicting sorafenib therapeutic responses for tumor recurrence. Finally, a similar approach could be applied to other cancer types for the development of therapeutic agents/inhibitors to improve the combination therapy of other antiangiogenic multikinase inhibitors, and such approaches could be validated in comparative mouse cancer models. Our approach of applying concordantly expressed comparative cancer genes to forecast targets and therapeutic agents in combination cancer therapy could be difficult to implement due to the limited availability of transcriptomic datasets with clinical features derived from comparative cancer models and cancer patients treated with and without therapeutic agents. Nevertheless, with increasing cancer incidence and mortality predicted worldwide in the coming decade, the successes of combination cancer therapies reported in multiple cancers will encourage more clinical trials on combination therapies with the goal of reducing the toxicity and cost to prolong the survival of cancer patients.

Conclusions

To improve HCC therapies and meet emerging needs worldwide, we conducted a large-scale comparative transcriptome analysis and identified 664 pcHCC genes (∼44% are innovative targets) associated with prognostic clinical features and enriched driver gene signatures suitable for developing sorafenib combination therapies. We provided lines of evidence to demonstrate the divergent mechanisms of sorafenib combination therapies in HCC and experimentally validated selected pcHCC genes by reducing their expression in HCC cells via a diminished innovative p38/STAT3 axis to hypersensitize HCC cells in sorafenib-combined treatments. We provided important insights on 664 pcHCC genes, and they include innovative targets and molecular signaling factors suitable for systematic targeting in sorafenib HCC combination therapy, which has the potential to prolong the life of HCC patients.

Materials and methods

Data sources and preprocessing

The transcriptome datasets were downloaded from the NCBI Gene Expression Omnibus [85]. We selected transcriptomes derived from tumors of human hepatocellular carcinoma (HCC)-only and mouse liver cancer models but excluded transcriptomes from cell lines (Supplementary Tables 1 and 2). The majority of mouse transcriptomes were annotated as HCC (∼70%) or liver cancer (∼30%) in genetically modified models with and without carcinogen induction, as detailed in Supplementary Table 2. The ratio of tumor-versus-normal samples for each probe was calculated for each dataset. If there was no normal sample in a dataset, a pool of normal samples was constructed from normal samples from other datasets of the same platform. Two independent HCC datasets from iCOD and TCGA_LIHC (the level 3 RNA-Seq data of HCC tumors) were downloaded based on the original report [27] and The Cancer Genome Atlas [86] using the R package TCGAbiolinks (https://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html) [87].

Integrated meta-analysis of HCC and mouse genes

All probes were mapped to the gene ID in the Ensembl database (https://www.ensembl.org, build 92), and mouse Ensembl gene IDs were then mapped to human Ensembl gene IDs [88]. An average ratio was calculated for the multiple probes mapped into one Ensembl gene ID. For all tumor/normal ratios, log2 transformations were performed, followed by the Z-score transformation of each sample to standardize data across samples in all datasets [89]. Student’s t-test was performed for each Ensembl gene ID.

Determination of comparative differentially expressed genes (DEGs) and gene set enrichment analysis (GSEA) in HCC

To determine the comparative DEGs, we defined the upregulated probes when the processed log2 tumor/normal ratios were more than 1 and the downregulated probes when the ratios were less than −1. The percentage of DEGs was calculated from the sum of upregulated or downregulated labels versus the total sample numbers. An adjusted percentage for defining the upregulated and downregulated genes was calculated using a formula (percentage of upregulated labels minus percentage of downregulated labels) for each Ensembl gene ID. For ensembled gene IDs with <25% presented in the samples, those with a p value of more than 0.05 in Student’s t-test, <5% total altered sample percentage, and 5% to −5% adjusted regulated percentage were filtered out. Only the Ensembl gene IDs having the same regulatory trend in both human and mouse samples were selected as comparative DEGs. Gene set enrichment analysis (GSEA) was applied to detect enriched pathways of the DEGs with FDR < 0.05 and P value <0.05 for consideration of significant enrichment [90].

Archived HCC samples

Archived HCC tissues and the corresponding normal tissues (n = 21) were collected from National Taiwan University Hospital with IRB approval [17]. Tissue was harvested after curative HCC surgery, frozen in liquid nitrogen, and stored in a −80 °C freezer immediately after being sliced into small pieces.

RNA isolation, RT-PCR and quantification of RNA expression

Total RNA was extracted from HCC tissues using TRIzol (Invitrogen) according to the manufacturer's instructions. RNA purity was assessed with a NanoDrop spectrophotometer (Thermo Scientific). RNA was reverse-transcribed into cDNA with oligo-dT primers using a Maxima H Minus First Strand cDNA Synthesis Kit (Thermo Scientific). RT-PCR was performed with Taq DNA Polymerase Master Mix RED (Ampliqon). PCR cycling conditions were 95 °C for 30 s, 60 °C for 30 s and 72 °C for 20 s for 40 cycles. The PCR products were separated through gel electrophoresis with DNA View (Tools Biotech) in a 1.5% agarose gel and identified through UV exposure. RNA expression quantification was performed with GelQuant. NET software provided by biochemlabsolutions.com. Target gene expression of certain HCC samples was normalized to the internal control GAPDH expression. The list of primers used for RT-PCR is included in Supplementary Table 11.

Transwell migration assay

Cell migration assays were conducted as described previously [84]. Each condition was performed in triplicate wells with seeding of a total of 3 × 104 cells/well for Mahlavu or 8 × 104 cells/well for HCC36 HCC cell lines. After 24 h, the migrated cells were stained with Giemsa solution (Merck) and counted in 100 × microscopic fields with ImageJ (NIH, USA).

Cell viability assays

Cell viability was assessed in triplicate wells for each condition using the alamarBlue assay (BioSource International) [84]. A total of 1 × 104 Mahlavu or HCC36 HCC cells/well were seeded in 6-well plates and incubated in a 37 °C incubator. Cell viability at days 0, 1, 3, 5, 7 and 9 post-seeding was analyzed by detecting the absorbance using a microplate reader. Sorafenib dissolved in dimethyl sulfoxide (DMSO) at various concentrations was added to the medium of HCC cell cultures.

Short hairpin RNA (shRNA)-mediated gene silencing

The productions of shRNA lentivirus particles were performed according to previous report [91]. The following shRNA vectors were used for silencing gene expression: shCENPW #3 (TRCN0000142180), shCENPW #4 (TRCN0000143715), shNCAPG #2 (TRCN0000145160) and shNCAPG #3 (TRCN0000144644) (National RNAi Core Facility, Taiwan). Transduction of Mahlavu or HCC36 cells with shRNA lentiviral particles was performed in 6-well plates containing the transduction adjuvant polybrene (8 μg/ml). Stable shRNA-expressing HCC cells were selected using puromycin for 2 weeks. Knockdown efficiency was assessed with Western blot.

Western blot analysis

The following antibodies against specific proteins were used for Western blot analysis: Actin (A00702, GenScript), CENPW (SC-137988), NCAPG (AP19113a, Abgent), MCM2 (BS1221, Bioworld), PCNA (BS1289, Bioworld), p-ERK (9101, Cell Signaling), ERK (9102, Cell Signaling), p-JNK (9255, Cell Signaling), JNK (9252, Cell Signaling), p-p38 (1229–1, Epitomics), p38 (54419, AnaSpec), p-STAT3 (9134, Cell Signaling) and STAT3 (Sc-482, Santa Cruz, USA). Total protein was extracted from the cells using RIPA buffer as described previously [84]. Quantification of total protein lysate was analyzed using a BCA protein assay kit (Pierce Chemical). Protein from each condition (20 μg) was separated with 8 ∼ 12% SDS–PAGE and then blotted onto PVDF membranes (Millipore). The protein-blotted PVDF membranes were further probed with specific antibodies (1:5000 dilution) overnight at 4 °C. The membranes were washed and then reprobed with HRP-conjugated secondary antibodies (1:10,000) at room temperature for 1 hr. The signal on protein-blotted membranes was visualized using ECL (Enhanced Chemiluminescence) reagents (SuperSignal, Pierce Chemical) and autoradiography films (Kodak Rochester, NY).

Statistical analysis and data visualization

All in vitro conditions were performed in triplicate, and values are shown as the mean ± 95% confidence interval (CI). The differences between groups were determined using Student's t-test, and statistical significance was considered when P < 0.05. All analyses were performed in the R statistical programming environment. All plots were produced with the ggplot2 package (https://ggplot2.tidyverse.org/). All scripts for data analysis and visualization are available upon request.

Declaration of Competing Interests

The authors declare no conflicts of interest.