Yu Pan1,2,3,4,5,6†
Geng-yuan Hu5,6,7†Shi Jiang1,2,3,4,5,6Shun-jie Xia1,2,3,4,5,6Hendi Maher1,2,3,4,5,6Zhong-jie Lin1,2,3,4,5,6Qi-jiang Mao1,2,3,4,5,6Jie Zhao1,2,3,4,5Liu-xin Cai1,2,3,4,5Ying-hua Xu8*Jun-jie Xu1,2,3,4,5*
Xiu-jun Cai1,2,3,4,5*- 1Department of General Surgery, Sir Run-Run Shaw Hospital, Zhejiang University, Hangzhou, China
- 2Key Laboratory of Laparoscopic Technology of Zhejiang Province, Hangzhou, China
- 3Zhejiang Minimal Invasive Diagnosis and Treatment Technology Research Center of Severe Hepatobiliary Disease, Hangzhou, China
- 4Zhejiang Research and Development Engineering Laboratory of Minimally Invasive Technology and Equipment, Hangzhou, China
- 5Zhejiang University Cancer Center, Hangzhou, China
- 6School of Medicine, Zhejiang University, Hangzhou, China
- 7Department of Gastrointestinal Surgery, Shaoxing People's Hospital, Shaoxing Hospital of Zhejiang University, Shaoxing, China
- 8Department of Oncology, Sir Run-Run Shaw Hospital, Zhejiang University, Hangzhou, China
Hepatocellular carcinoma (HCC) is a deadly tumor with high heterogeneity. Aerobic glycolysis is a common indicator of tumor growth and plays a key role in tumorigenesis. Heterogeneity in distinct metabolic pathways can be used to stratify HCC into clinically relevant subgroups, but these have not yet been well-established. In this study, we constructed a model called aerobic glycolysis index (AGI) as a marker of aerobic glycolysis using genomic data of hepatocellular carcinoma from The Cancer Genome Atlas (TCGA) project. Our results showed that this parameter inferred enhanced aerobic glycolysis activity in tumor tissues. Furthermore, high AGI is associated with poor tumor differentiation and advanced stages and could predict poor prognosis including reduced overall survival and disease-free survival. More importantly, the AGI could accurately predict tumor sensitivity to Sorafenib therapy. Therefore, the AGI may be a promising biomarker that can accurately stratify patients and improve their treatment efficacy.
Introduction
Globally, hepatocellular carcinoma (HCC) is the sixth most commonly diagnosed cancer and the fourth leading cause of cancer-related deaths (1). Despite advances in the treatment of HCC, its prognosis remains unsatisfactory, with a 5-year overall survival (OS) rate of 25–55% (2–4). Local recurrence, distal metastasis, and resistance to conventional therapy are the leading causes of HCC progression into late-stage cancer with limited treatment options. Genetic mutations, chromosomal instability, epigenetic changes, and molecular signaling pathway dysregulation are reported causes of hepatocellular carcinogenesis (5). Therefore, advances in the field of molecular oncology are urgently required to improve the prognosis of HCC.
Sorafenib is an oral multitargeted drug that inhibits the activity of several tyrosine kinases (6). Thus, Sorafenib can effectively suppress angiogenesis and cancer proliferation and induce tumor cell apoptosis. Since it was first approved by the US Food and Drug Administration in 2007 as the first-line treatment for advanced HCC, Sorafenib has shown favored clinical benefits. In the Sorafenib Hepatocellular Carcinoma Assessment Randomized Protocol (SHARP) trial, patients who received Sorafenib therapy showed significantly higher median OS compared with the control group (10.7 vs. 7.9 months, respectively) and a 31% reduction in the risk of death (7). In the Asia-Pacific trial, Sorafenib provided a clinical benefit, extending the median survival benefit by 2 months (8). Currently, Sorafenib is still applied as the first-line therapy for advanced HCC patients; even several other tyrosine kinases inhibitors have been evaluated by comparing them to Sorafenib, which did not demonstrate an improvement of prognosis (9–11). However, a large number of patients with HCC show poor response to Sorafenib due to the heterogeneity of the disease and the complex tumor-associated molecular signaling, which lacks generally accepted predictive biomarkers. Furthermore, primary and acquired resistance to Sorafenib are commonly reported and limit the clinical advantages of the drug (12–15).
The Warburg phenotype is a common hallmark of cancer cells, characterized by enhanced glycolysis, even under physiological oxygen conditions (16, 17). By shifting glucose metabolism from oxidative respiration to aerobic glycolysis, tumor cells display enhanced glucose metabolism for producing efficient energy and various metabolic intermediates, which are indispensable for the synthesis of macromolecules and new organelles. Several oncogenes, including Ras, Myc, and HIF1, were reported to drive metabolic adaptations toward aerobic glycolysis (18, 19). Enhanced aerobic glycolysis has been demonstrated to exhibit consistent prognostic patterns and is associated with Sorafenib resistance (20). Increased expression of aerobic glycolysis-related genes, including solute carrier family 2 member 1 (SLC2A1), solute carrier family 2 member 2 (SLC2A2), glucose-6-phosphate dehydrogenase (G6PD), glypican 1 (GPC1), procollagen-lysine, 2-oxoglutarate 5-dioxygenase 2 (PLOD2), and lactate dehydrogenase A (LDHA), has been reported to be associated with aggressive HCC (21–29). Accelerated glucose uptake and lactate synthesis were observed as responses to Sorafenib treatment (30–32). New therapeutic approaches have been reported to attenuate Sorafenib resistance by inhibiting key glycolytic enzymes including 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3 (PFKFB3), hexokinase 2 (HK2), and pyruvate kinase M1/2 (PKM2) (33–36). Tumor metabolic heterogeneity is reported to be relevant to tumor subtypes and prognosis (37), but whether heterogeneity in distinct metabolic pathways can be used to stratify HCC into clinically relevant subgroups has not been well-established.
Aerobic glycolysis is a complex biological process involving numerous genes. Thus, constructing a gene signature based on multiple glycolysis-related genes is supposed to be more suitable to represent the aerobic glycolysis pathway than single gene. The development of genomic techniques has unveiled extensive biological information that can be used to explore the underlying mechanisms of tumorigenesis and progression. In this study, we constructed a model named aerobic glycolysis index (AGI) to evaluate the signal of aerobic glycolysis, by utilizing genomic data of hepatocellular carcinoma from The Cancer Genome Atlas (TCGA) project. The AGI was calculated based on the expression of 14 glycolysis-related genes (SLC2A1, SLC2A2, G6PD, LDHA, GPC1, HMMR, PLOD2, GOT2, STC2, CENPA, RARS1, HOMER1, SRD5A3, and TKTL1). The abbreviations list and their expansions for these glycolytic genes are summarized in Supplementary Table 1. Then, we assessed whether the AGI was a predictive marker for the prognosis of HCC and sensitivity to Sorafenib. Finally, we used in vitro experiments to confirm that AGI was associated with Sorafenib resistance. Significantly, we established a methodology to quantify aerobic glycolysis signaling. The AGI was found to be a robust prognostic biomarker of HCC and a predictive factor of the response to Sorafenib.
Materials and Methods
Specimen Collection and RNA Sequencing
In total, 102 pairs of formalin-fixed paraffin-embedded HCC and corresponding normal tissue specimens from the Sir Run Shaw Hospital (SRRSH) were collected. This research was approved by the Institutional Review Board of the SRRSH of Zhejiang University. The patients provided informed consent for the use of their specimens. For RNA sequencing, TRIzol (Invitrogen, USA) was used to extract total RNA. After checking the RNA purity, integrity, and concentration, RNA sequencing was performed on an Illumina platform.
Data Acquisition and Processing
Fragments per kilobase of transcript per million mapped reads (FPKM) RNA-seq data and the clinical characteristics of TCGA samples were downloaded from the UCSC Cancer Browser database. Gene mutation data and copy number information of the TCGA samples were acquired from the cBioPortal database. RNA-seq and clinical data of the LIRI-JP cohort were downloaded from the HCCDB database. The GSE14520, GSE25097, GSE36376, GSE64041, GSE76427, GSE109211, and GSE73571 expression profile was obtained from the Gene Expression Omnibus (GEO) database. The proteomics data of TCGA samples were downloaded from the TCPA database. Drug sensitivity data of HCC cell lines were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database and Cancer Cell Line Encyclopedia (CCLE) database. The patient characteristics of TCGA, LIRI-JP, GSE14520, and SRRSH datasets are summarized in Supplementary Table 2.
Development of the AGI
Univariate Cox regression was applied to detect the aerobic glycolysis genes related to prognosis. The least absolute shrinkage and selection operator (LASSO) Cox regression model was performed to determine the coefficients for model construction with an optimal log λ (38). The AGI was established with the following formula: Risk score = expression of Gene 1 × β1 + expression of Gene 2 × β2+…expression of Gene n × βn (β was the weighted coefficient of each gene). For gene expression measured by quantitative real-time PCR (qPCR), the AGI was calculated as the method described by Zheng et al. (39). Total RNAs from cells were extracted using TRIzol (Invitrogen, USA) according to the manufacturer's instructions. Complementary DNA (cDNA) was synthesized using Hifair® II 1st Strand cDNA Synthesis SuperMix for qPCR (Yeasen, Shanghai, China). qPCR was performed using Hieff UNICON® qPCR SYBR Green Master Mix (Yeasen, Shanghai, China). Measurement was carried out by Roche LightCycler 480. Analysis was carried out using the ΔΔCt method. All assays were performed in triplicates, and results were plotted as the mean ± SD. Primer sequences are listed in Supplementary Table 3.
Gene Set Enrichment Analysis and Gene Set Variation Analysis
Datasets were divided into two groups according to tissue types or AGI scores. Annotated gene sets were downloaded from the MSigDB database. Gene Set Enrichment Analysis (GSEA) was performed using the R package “GSEAbase.” Annotated drug sets were downloaded from the DSigDB database. Gene Set Variation Analysis (GSVA) was performed using the R package “GSVA.”
Cell Culture of Liver Cancer Cell Lines
Liver cancer cell lines (SK-Hep-1, Huh7, HepG2, HCCLM3) were purchased from the American Tissue Culture Collection (Manassas, VA, USA) and cultured in accordance with the recommended guidelines. Sorafenib-resistant HCC cell lines were cultured with Sorafenib as previously reported (40, 41).
Cell Viability Test
Cells were seeded in 96-well plates in replicates of three. After incubation with Sorafenib for 48 h, cell viability analysis was performed using the Cell Counting Kit-8 (CCK-8) (Yeasen, Shanghai, China).
Apoptosis Assay
Cells were seeded in six-well plates and mock treated or treated with drugs [Sorafenib, 2-deoxy-D-glucose (2-DG), or a combination of Sorafenib and 2-DG] for 48 h before apoptosis assays. Cell apoptosis was determined using the PI/annexin V-FITC Apoptosis Kit (MULTI SCIENCES, Hangzhou, China).
Transwell Assay
Cells (1 × 105) in serum-free medium were seeded into the upper chambers of Transwell (Corning, Corning, NY, USA), and medium with 10% fetal bovine serum (FBS) was seeded into the lower chambers for 24 h in a humidified incubator at 37°C in 5% CO2. The cells remaining in the upper chamber were carefully removed using a cotton swab, and cells that migrated to the lower membrane surface were fixed in 4% paraformaldehyde and stained with crystal violet. The experiments were repeated three times.
Glucose Consumption and Lactic Acid Assays
Cells were seeded into six-well plates at a density of 1 × 106 cells and cultured overnight. Glucose consumption and lactic acid production were detected using glucose assay kit (Solarbio® BC2500) and LA assay kit (Solarbio® BC2230) according to the instruction of the manufacturers, respectively. The experiments were repeated three times.
Statistical Analysis
The univariate Cox regression, LASSO Cox regression model, and multivariate Cox regression model were performed. The OS and disease-free survival (DFS) were compared using the Kaplan–Meier method with the log-rank test. Non-parametric tests or Student's t-tests were used to determine the significance of the differences between the subgroups and clinicopathological characteristics. Spearman's correlation test was used to assess the relationship between the AGI and biological pathways and clinicopathological parameters. Statistical analyses were performed using R software (Version 3.6.0). A P < 0.05 was considered statistically significant, and all P-values were two-tailed.
Results
Establishment of the AGI
This study was conducted according to the flow chart shown in Figure 1. Initially, several key genes of glycolysis were observed to be commonly upregulated in tumor tissues compared with normal tissues in a series of datasets, including TCGA, LIRI-JP, GSE14520, GSE25097, GSE36376, GSE64041, and GSE76427 (Figure 2A). The GSEA results revealed that the glycolysis signaling pathway was significantly enriched in HCC tumor tissues compared with normal tissues (Figure 2B). These results suggested that the glycolysis signaling pathway may contribute to malignant tumor phenotypes. Because aerobic glycolysis is a complex biological process involving hundreds of genes, using a gene signature comprising multiple genes can predict the tumor characteristics and prognosis more accurately than a single gene. The univariate Cox regression analysis was conducted to determine the prognostically relevant genes related to the aerobic glycolysis, and 80 relevant genes were identified (Figure 2C). To further simplify the gene signature of aerobic glycolysis, LASSO regression was performed based on these prognostically relevant genes. Finally, 14 genes were selected to establish the AGI according to the partial likelihood deviance method with an optimal log λ (Supplementary Figures 1, 2). Additionally, as shown in Figure 2D, the correlations between the AGI and selected aerobic glycolysis-related genes were statistically significant. Furthermore, the correlations between the AGI and several important genes coding the rate-limiting enzymes of glucose metabolism were examined. The results demonstrated that the AGI correlated closely to these genes (HK2, r = 0.43, P < 2.2e−16; PFKP, r = 0.39, P = 1.4e−14; PFKFB3, r = 0.27, P = 2e−7; PKM2, r = 0.66, P < 2e−16), as shown in Supplementary Figure 3.
Figure 1. Flowchart presenting the establishment and validation of the gene signature.
Figure 2. Construction of the aerobic glycolysis index (AGI) model and validation of the AGI in tumor and normal tissues. (A) Heatmap of glycolysis-related gene expression in different datasets. (B) Gene Set Enrichment Analysis (GSEA) of the glycolysis pathway in GSE14520. (C) Bar plot showing the hazard ratio of glycolysis-related genes in The Cancer Genome Atlas (TCGA) cohort using the univariate Cox regression. The bars represent the 95% CI. (D) Correlation between the AGI and the selected signature genes in the TCGA cohort. (E) Boxplots showing AGI differences in normal and tumor tissues in the TCGA, GSE64041, GSE14520 and LIRI JP datasets. (F) Boxplots showing AGI differences in normal and tumor tissues in the SRRSH set. (G) Receiver operating characteristic (ROC) curves for tissue type prediction using the AGI as the predictor.
After the establishment of the AGI, we first tested the capacity of AGI as an acceptable indicator of aerobic glycolysis. The AGI could distinguish HCC tumor tissues from normal tissue samples in the TCGA dataset that significantly higher AGI was observed in tumor samples, indicating enhanced aerobic glycolysis activity in these tumor samples (P = 8.1e−13, Figure 2E). To further validate the power of the AGI, we examined this index in three other HCC datasets. As shown in Figure 2E, a significantly higher AGI was observed in HCC tumor samples compared with normal tissue samples in all of three datasets (GSE64041, P = 8.1e−10; GSE14520, P < 2e−16; LIRI-JP, P < 2e−16). Furthermore, RNA sequencing data from our center (SRRSH) were applied and confirmed the reliability of AGI in distinguishing HCC tumor tissues from normal tissue samples (P < 2e−16) (Figure 2F). Then, the receiver operating characteristic (ROC) curve and area under the curve (AUC) scores were further evaluated to quantify the accuracy of the AGI to classify tumor and normal tissues. High prediction accuracies were achieved in all datasets, ruling out the possibility of over fitting (TCGA, AUC = 0.916; GSE64041, AUC = 0.825; GSE14520, AUC = 0.945; LIRI-JP, AUC = 0.846, SRRSH set, AUC = 0.902, Figure 2G).
Association of the AGI With Genomic and Proteomic Alterations
Using the optimal cutoff value based on prognostic effects, the patients in the TCGA cohort were stratified into high and low AGI groups. The distribution of gene expression in the high and low AGI groups is represented in Figure 3A. In particular, we compared the genes encoding key proteins and enzymes of aerobic glycolysis, and most of them were highly expressed in the high AGI group, supporting the close correlation between aerobic glycolysis and the AGI (Figure 3B).
Figure 3. The landscape of biological processes and characteristics of the aerobic glycolysis index (AGI) subgroups. (A) Heatmap of common differentially expressed genes based on the expression data in the high and low AGI groups. (B) Box plots showing the expression of the selected glycolysis-related genes in The Cancer Genome Atlas (TCGA) cohort. (C) Proportion of TP53 and CTNNB1 mutations in the high and low AGI groups. (D) Box plots showing the AGI in patients with TP53 mutations and wild-type TP53 (left) and CTNNB1 mutations and wild-type CTNNB1 (right). (E) The oncoPrint of copy number variations was constructed in the high and low AGI subgroups. (F) Activated gene sets enriched in the high AGI subgroup. (G) Suppressed gene sets enriched in the high AGI group. (H,I) Proteins positively (H) and negatively (I) correlated with the AGI (P < 0.05 for all proteins) based on reverse-phase protein arrays analysis of 181 samples from the TCGA using Spearman's rank correlation.
Previous reports demonstrated that aerobic glycolysis genes could be regulated by transcription factors such as p53, c-Myc, and HIF-1α. In our study, the high AGI group showed an increased proportion of TP53 mutations and a decreased proportion of CTNNB1 mutations (Figure 3C), which were the most common mutations in HCC. Accordingly, patients with a TP53 mutation showed a higher AGI than patients with wild-type TP53. Conversely, patients with a CTNNB1 mutation showed a lower AGI than patients with wild-type CTNNB1 (Figure 3D). Next, we investigated the distribution of copy number variations between the high and low AGI groups. The high AGI group showed an increased amplification frequency of Myc, AGO2, EXT1, RAD21, EIF3E, RSPO2, RECQL4, RUNX1T1, NBN, PAG1, and HEY1 (Figure 3E). These outcomes may provide novel ideas for investigating the mechanism of tumor aerobic glycolysis and copy number variation.
GSEA of the transcripts in the two groups revealed that gene sets considered to be markers of high malignancy were enriched in the high AGI group, including those E2F targets, Myc targets, epithelial–mesenchymal transition (EMT) regulators, and G2M checkpoints (Figure 3F and Supplementary Table 4). Conversely, gene sets related to oxidative phosphorylation, peroxisomes, xenobiotic metabolism, and other metabolic processes were enriched in the low AGI group (Figure 3G and Supplementary Table 4). As the EMT pathway enriched in the high AGI group, we were interested in the relationship between the AGI and tumor invasion capability. The AGI in four liver cancer cell lines (SK-hep-1, Huh7, HepG2, and HCCLM3) was calculated. The result demonstrated that cell lines with high AGI (SK-hep-1, Huh7, and HCCLM3) had increased glucose intake and lactic acid level as compared with the cell line with low AGI (HepG2). More importantly, the Transwell assay revealed that cell lines with high AGI exhibited the greater invasion capability (Supplementary Figure 4). Previous studies also have reported a close relationship between aerobic glycolysis and angiogenesis (42–44). In the present study, we compared the expression of several genes involved in angiogenesis between the high and low AGI groups. The results demonstrated that the high AGI group had increased expression of vascular endothelial growth factor A (VEGFA) and VEGFB as compared with the low AGI group (Supplementary Figure 5).
Further analysis of proteomic data revealed that the AGI was strongly correlated with the expression of cell-cycle-related proteins, including pChk1, pChk2, and CyclinB1, and other tumor hallmark proteins, such as eIF4G, pPI3K, Src, Smad1, and Smad3 (Figure 3H). Furthermore, the AGI was negatively correlated with the expression of Rb, phosphatase and tensin homolog (PTEN), cleaved poly(ADP-ribose) polymerase (PARP) (Figure 3I). Collectively, the analysis of gene mutations, copy number variations, classical signaling pathway gene sets, and proteomic data suggests that the AGI is associated with cell proliferation, tumor progression, and the inhibition of apoptosis in HCC.
Correlation Between the AGI and Clinical Characteristics
Our data showed that the AGI was closely associated with TP53, Myc, cell cycle, and EMT pathways, and therefore, we next examined whether the AGI is associated with tumor progression and metastasis. We first analyzed the relationship between the AGI and clinical characteristics of HCC patients. A significantly increased AGI was observed in patients with higher tumor grades (Figure 4A), advanced T stages, and tumor–node–metastasis (TNM) stages (Figures 4B,C) and vascular invasion (Figure 4D). Here, we defined early tumor recurrence as a tumor recurring within 2 years after primary treatment and late recurrence as cancer recurring after 2 years. The AGI was significantly higher in patients with early recurrence than in those with late recurrence (Figure 4E).
Figure 4. Clinicopathological significance and prognosis prediction value of the AGI. (A) Tumor differentiation grade. (B) T stage. (C) Tumor–node–metastasis (TNM) stage. (D) Vascular invasion status. (E) Recurrence status. (F) Kaplan–Meier plot analysis of overall survival (OS) in the high and low AGI groups. (G) Kaplan–Meier plot analysis of disease-free survival (DFS) in the high and low AGI groups. (H,I) Forest plot showing the prognostic value of the AGI and clinical characteristics using univariate (H) and multivariate (I) analysis. (J) Time-dependent receiver operating characteristic (ROC) analysis comparing the AGI in predicting the 1-, 3-, and 5-year OS. (K) Time-dependent ROC analysis comparing the AGI and clinical characteristics in 5-year OS. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.
Patients with a higher AGI were associated with a worse prognosis in terms of OS [hazard ratio (HR), 3.43; P < 0.001; Figure 4F] and DFS (HR, 2.07; P < 0.001; Figure 4G). The univariate Cox regression analysis revealed that the AGI (HR, 4.132), T stage (HR, 1.675), and TNM stage (HR, 1.661) were risk factors of a worse HCC prognosis (Figure 4H), and the multivariate Cox regression analysis showed that the AGI was an independent risk factor of poor prognosis (HR, 4.865; Figure 4I). The ROC curve showed that the AGI could accurately predict the 1-year (AUC, 0.787), 3-year (AUC, 0.751), and 5-year OS (AUC, 0.714) of HCC (Figure 4J), which was superior to conventional clinical parameters such as the tumor grade, status of vascular invasion, and TNM stage (Figure 4K).
Prognostic Value Validation of the AGI
A significantly increased AGI was observed in patients with higher TNM stages in three independent datasets (Figures 5A–C). The prognostic value of the AGI in HCC was validated in three independent datasets, including GSE14520 from GEO, LIRC from ICGC, and the SRRSH set from our center. By dividing the datasets into two groups according to the AGI, the distribution of the gene expression profiles and status of patients were consistent with the AGI (Figures 5D–F). Similar to the TCGA dataset, the high AGI group showed worse OS (GSE14520: HR, 2.13; Figure 5G; LIRC: HR, 2.85; Figure 5H; SRRSH set: HR, 2.53; Figure 5I) and DFS (GSE14520: HR, 1.67; Figure 5K; SRRSH set: HR, 1.77; Figure 5L) compared with the low AGI group. The ROC curve demonstrated that the AUCs of the AGI in predicting the 5-year OS of GSE14520, LIRC, and SRRSH sets were 0.676, 0.630, and 0.621, respectively (Figure 5J), indicating a robust prognostic value of the AGI.
Figure 5. Validating the prognostically predictive value of the aerobic glycolysis index (AGI) in validation datasets. (A–C) Value of the AGI in different tumor stages in GSE14520 (A), LIRI-JP (B), and SRRSH set (C). (D–F) The heatmap and distribution of the 14 AGI-related gene expression profiles in GSE14520 (D), LIRI-JP (E), and SRRSH set (F). (G–I) Kaplan–Meier plot analysis of overall survival (OS) in the high and low AGI subgroups in GSE14520 (G), LIRI-JP (H), and SRRSH set (I). (J) Time-dependent ROC analysis comparing the predictive value of the AGI for 5-year OS in the three datasets. (K,L) Kaplan–Meier plot analysis of disease-free survival (DFS) in the high and low AGI subgroups in GSE14520 (K) and SRRSH set (L). *P < 0.05, **P < 0.01, ***P < 0.001.
AGI Predicts Sensitivity to Sorafenib in Both HCC Cell Lines and Patients
Sorafenib is the first-line therapy for advanced HCC. The increased expression of aerobic glycolysis-related genes has been demonstrated to promote Sorafenib resistance (45). Our results showed the AGI was closely related to advanced tumor stages and poor tumor differentiation. Moreover, a high AGI was strongly correlated with EMT, Myc, and cell cycle pathways (Figure 3B), which had been reported to be associated with impaired Sorafenib sensitivity and a worse prognosis (45). We were interested in the relationship between the AGI and Sorafenib sensitivity and wondered if the AGI could be a potential biomarker to predict drug sensitivity. Thus, the GSVA was performed, and the results revealed that the low AGI group appeared to be more sensitive to Sorafenib (Figure 6A). Then, we verified the relationship between Sorafenib sensitivity and the AGI in vitro using the data from the GDSC database, which showed that HCC cell lines with a lower AGI were more sensitive to Sorafenib (Figure 6B). Furthermore, a positive correlation between the AGI and the natural logarithm of the IC50 was observed (r = 0.61, P = 0.02, Figure 6C). Similarly, using HCC cell line data from the CCLE database, cell lines with a higher AGI showed decreased Sorafenib sensitivity (Figure 6D), and a positive correlation between the AGI and the natural logarithm of the EC50 of Sorafenib was observed (Figure 6E). Although not statistically significant due to a limited sample size, the trends suggested a potential relationship between the AGI and Sorafenib sensitivity in vitro. We next calculated the AGI in patients who received Sorafenib therapy using the transcription data from GSE109211 (the STORM trial). The patients who did not respond to Sorafenib showed a higher AGI (Figure 6F) and upregulated expression of AGI-related genes (Figure 6G). The ROC curve of the AGI in predicting the sensitivity to Sorafenib retrieved an AUC of 0.879 (Figure 6H), indicating that the AGI is a reliable biomarker in selecting suitable patients for Sorafenib treatment.
Figure 6. Association between the aerobic glycolysis index (AGI) and Sorafenib resistance. (A) The Gene Set Variation Analysis (GSVA) results showed that the Sorafenib sensitivity signature was enriched in patients with a low AGI. (B) AGI positively was correlated with the IC50 of Sorafenib in hepatocellular carcinoma (HCC) cell line data from the Genomics of Drug Sensitivity in Cancer (GDSC) database. (C) The IC50 of Sorafenib in HCC cell line data from the GDSC database with high and low AGIs. (D) AGI positively correlated with the EC50 of Sorafenib in HCC cell line data from the Cancer Cell Line Encyclopedia (CCLE) database. (E) EC50 of Sorafenib in HCC cell line data from the CCLE database with high and low AGIs. (F) Patient sensitive to Sorafenib presented significantly low AGI. (G) The heatmap and distribution of the 14 AGI-related gene expression profiles in GSE109211. (H) Receiver operating characteristic (ROC) analysis showed an area under the curve (AUC) of 0.879 for the AGI in predicting the response to Sorafenib.
Increased AGI in Sorafenib-Resistant HCC Cell Lines
To further study the relationship between the AGI and Sorafenib resistance, the expression of AGI-related genes was evaluated in different cell lines (SK-hep-1 and Huh-7) at 0, 24, 48, and 72 h after Sorafenib treatment using qPCR. The levels of AGI-related genes substantially elevated following the treatment with Sorafenib, resulting in an increased AGI (Figure 7A). These results suggested that the AGI and underlying metabolic remodeling may be closely related to Sorafenib treatment. In previous studies, we reported several HCC cell lines treated with Sorafenib for long term, which could be considered as models of Sorafenib-resistant (SR) cell lines (40, 41, 46, 47). The RNA-sequencing data revealed a high AGI in SR cell lines (SK-hep-1 and Huh7) (Figure 7B). Transcript data from GSE73571 also demonstrated an elevated AGI during the acquisition of Sorafenib resistance (Figure 7C). Thus, we proposed the following hypothesis. On the one hand, tumor cells adapted to Sorafenib therapy by shifting to aerobic glycolysis. On the other hand, cells predominantly using aerobic glycolysis were also selected by Sorafenib. Both of these processes resulted in cells with enhanced aerobic glycolysis activity. Because the AGI is applied as a marker of aerobic glycolysis signaling activity, we speculate that the inhibition of aerobic glycolysis may enhance the sensitivity of Sorafenib. A combination of Sorafenib and 2-deoxy-D-glucose (2-DG), a commonly used inhibitor of aerobic glycolysis in vitro, was added to cell cultures, resulting in the inhibition of cell proliferation (Figure 7D) and an increased fraction of apoptotic cells (Figure 7E).
Figure 7. Aerobic glycolysis index (AGI) is increased in Sorafenib-resistant hepatocellular carcinoma (HCC) cells. (A) Relative expression of AGI-related genes in HCC cell lines (left, SK-hep-1; right, Huh7) incubated with Sorafenib (5 μM) for 24, 36, and 72 h. (B) Distribution of the 14 AGI-related gene expression profiles in parental and Sorafenib-resistant HCC cell lines (SK-hep-1, Huh7). (C) AGI of Sorafenib-sensitive and Sorafenib-resistant xenografts from the GSE73571 dataset. (D) Combination of 2-DG and Sorafenib resulted in significantly decreased cell viability. (E) Combination of 2-DG and Sorafenib enhanced the apoptosis of Sorafenib-resistant cell lines (SK-hep-1 SR and Huh7 SR). *P < 0.05, **P < 0.01, ***P < 0.001.
Discussion
In the present study, we developed and validated an aerobic glycolysis-related gene signature (named the AGI) to predict the prognosis and Sorafenib sensitivity in patients with HCC. This AGI showed high accuracy in detecting HCC tumors and normal tissues. A low AGI score was significantly associated with early tumor stages, good differentiation, and better OS and DFS in a series of cohorts. Interestingly, the AGI score was correlated with the sensitivity of HCC cell lines to Sorafenib. More importantly, we demonstrated that the AGI could predict the response of patients to Sorafenib using data from clinical trials. Additionally, we observed that the AGI was elevated during the acquisition of Sorafenib resistance, which provides useful information for the development of a potential strategy to enhance Sorafenib sensitivity.
Aerobic glycolysis is a common biological process by which cancer cells tend to produce ATP by decomposing glucose or glycogen into lactic acid at a higher pace despite the presence of abundant oxygen. As an indicator of tumors, aerobic glycolysis activity is regulated by transcription factors, glucose transporters, and key enzymes of glucose metabolism. More importantly, aerobic glycolysis is related to multiple key cell signaling pathways, including PI3K/Akt, mTOR, and AMPK, and tightly associated with various cellular activities, including cell proliferation and EMT. In our study, the AGI was derived from a model consisting of genes encoding glucose transporters (SLC2A1, SLC2A2), key enzymes of glucose metabolism (G6PD, LDHA), and several other genes related to glycolysis (GPC1, HMMR, PLOD2, GOT2, STC2) (48–54), thereby supporting the use of the AGI as a marker of aerobic glycolysis activity. Several important genes related to glucose metabolism including HK2, PFK, and PKM2 are missed in the model of AGI, which had been reported to be associated with poor prognosis of HCC (55–57). The AGI was constructed using a bioinformatics method called LASSO regression. All the genes related to glycolysis were included as factors, and the LASSO regression selected factors to construct a model with minimal bias and acceptable reliability. Although HK2, PFK, and PKM2 were not in the model of AGI, we have calculated and found significant correlation between AGI and these genes. Elevated aerobic glycolysis activity was reported to result in a poor prognosis of multiple solid tumors including HCC (45, 58–60). Furthermore, several microRNAs, such as miR-383, miR-142-3p, and miR-100-5p, were reported to target LDHA, which subsequently inhibited cell proliferation, invasion, and glycolysis (61–63). Shang et al. reported that the transcription factor FOXM1 promoted glycolysis by transactivating SLC2A1 expression (64). In a respective cohort of 192 patients, the glucose transporter GLUT1 was significantly upregulated in HCC tumor tissues and was an independent risk factor of poor OS and relapse-free survival (24). Lu et al. reported that elevated G6PD expression contributed to the enhanced migration and invasion of HCC cells by inducing EMT (65), which was consistent with the correlation between the AGI and EMT signaling in the present study.
A positive correlation between aerobic glycolysis activity and Sorafenib resistance in both HCC cell lines and patients was observed in this study. Patients with a high pretreatment AGI tended to develop resistance to Sorafenib. Previously, Ma et al. and Li et al. found that increased aerobic glycolysis enhanced Sorafenib resistance in both HCC cell lines and xenografts (36, 66, 67). Key enzymes and transcription factors involved in aerobic glycolysis contributed to Sorafenib resistance, through reprogramming and redox adaptation (68), interacting with voltage-dependent anion channel (VDAC) and subsequently inhibiting mitochondrial apoptosis (69, 70) and increasing the expression of HIF-1α and c-Myc, thereby activating various cellular signals related to drug resistance (71, 72). In our previous review, we believed that Sorafenib resistance was associated with complex mechanisms, including metabolic remodeling, microenvironmental interplay, cellular signaling changes, genomic instability, and cancer stem cells (20). In our study, the levels of AGI-related genes substantially elevated following the treatment with Sorafenib, suggesting metabolic switch of glucose metabolism. Tesori et al. also reported metabolic shift toward glycolysis in HCC cells treated with Sorafenib in 48 h (31). Another study by Fiume et al. found that Sorafenib could hinder oxidative phosphorylation and stimulate aerobic glycolysis (32). Enhanced aerobic glycolysis activity was observed during the acquisition of Sorafenib resistance and reflected as an increased AGI. Thus, we speculate that the inhibition of aerobic glycolysis may enhance the sensitivity of Sorafenib. By inhibiting aerobic glycolysis activity, 2-DG resensitized HCC cells to Sorafenib therapy. Furthermore, several drugs targeting glycolysis-related factors such as Metformin, Aspirin, Genistein, Simvastatin, and Proanthocyanidin B2 have been shown to be effective in reversing Sorafenib resistance (33, 34, 36, 55, 70), indicating the aerobic glycolysis pathway as a promising target for exploring new therapies.
Notably, the present study highlighted that the AGI is a reliable biomarker in predicting the response to Sorafenib therapy. To date, several clinical and biological biomarkers have been proposed to evaluate responses to Sorafenib. The GIDEON trial revealed that patients with preserved liver function exhibited better OS after treatment with Sorafenib (73). Similarly, clinical characteristics such as Barcelona Clinic Liver Cancer (BCLC) stage, viral status, and Sorafenib-related adverse events were predictive of better survival (74, 75). As for biological biomarkers, Miyahara et al. reported that high levels of serum cytokines at baseline predicted poor outcomes in HCC patients treated with Sorafenib therapy (76). Arao et al. demonstrated that FGF3/FGF4 amplification was observed in 30% of HCC patients responding to Sorafenib (77). Recently, several exploratory studies investigated the roles of microRNAs in Sorafenib resistance and reported that several upregulated/downregulated microRNAs, long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs) were predictive biomarkers of survival outcomes following Sorafenib therapy (46, 47, 78–80). Different from previous studies, the present study directly evaluated the performance of the AGI in predicting Sorafenib responses using the ROC curve and observed an optimistic AUC of 0.88, which may be more accurate and suitable in clinical practice.
The present study has several limitations. First, although the available datasets with requisite gene transcript data and clinical and treatment outcome information were all included, the predictive effectiveness of the AGI was evaluated in only a few datasets. A more careful examination is required to further confirm the accuracy of the AGI using larger and multicenter clinical cohorts in the future. Second, the AGI was derived from a model of 14 gene transcripts including several genes with very low weight or minimal detectability. A simplified AGI with fewer key genes is required to improve the robustness and clinical utility of this model. Third, the AGI was demonstrated to be essential in Sorafenib resistance and associated with various tumor hallmarks. However, several genes used to construct the AGI have not yet been reported as prognostic factors of HCC and biomarker of Sorafenib, and their underlying mechanism remains unknown.
In conclusion, we developed a gene signature based on aerobic glycolysis-related genes by integrating several transcriptomic profiles, which showed great promise for predicting prognosis and the response of HCC to Sorafenib. The AGI described in our study can be developed as a predictive biomarker for Sorafenib therapy.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: FPKM RNA-Seq data and clinical information of TCGA LIHC was downloaded from the UCSC Cancer Browser (https://xenabrowser.net/datapages). Gene mutation data and GISTIC 2.0 segmentation scores as well as threshold copy number calls for the TCGA LIHC samples were acquired from the cBioPortal (http://www.cbioportal.org). RNA-seq and clinical data of the LIRI-JP cohort was downloaded from the HCCDB (http://lifeome.net/database/hccdb/download.html). The GSE14520, GSE25097, GSE36376, GSE64041, GSE76427, GSE109211, and GSE73571 expression profile was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The proteomics data of TCGA LIHC was downloaded from the TCPA database (https://www.tcpaportal.org/tcpa/index.html). Drug sensitivity data of HCC cell lines was obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org) and Cancer Cell Line Encyclopedia (CCLE) database (https://portals.broadinstitute.org/ccle/data).
Ethics Statement
The studies involving human participants were reviewed and approved by Institutional Review Board of the Sir Run Run Shaw Hospital of Zhejiang University. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
YP and G-yH: conceptualization. YP and SJ: data curation. YP, G-yH, SJ, and S-jX: formal analysis. X-jC, J-jX, Y-hX, and L-xC: funding acquisition. Z-jL and Q-jM: investigation. YP, Q-jM, and J-jX: methodology. X-jC and J-jX: project administration. YP, G-yH, and S-jX: software. X-jC, J-jX, and Y-hX: supervision, writing review and editing. YP, JZ, L-xC, and Z-jL: validation. YP: visualization. YP and HM: writing original draft. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the National Natural Science Foundation of China under Grant No. 81827804 (to X-jC), No. 81772546 (to X-jC), and No. 81902367 (to X-jC); Key Research and Development Project of Zhejiang Province under Grand No. 2018C03083 (to X-jC); Zhejiang Provincial Natural Science Foundation of China under Grant No. LQ19H160026 (to J-jX), No. LQ18H160010 (to JZ), and No. Y15H160052 (to L-xC); China Post-doctoral Science Foundation under Grant No. 2020M671755 (to J-jX); Zhejiang Clinical Research Center of Minimally Invasive Diagnosis and Treatment of Abdominal Diseases Grant No. 2018E50003 (to J-jX); and Hepatobiliary and Pancreatic Cancer Research of Hubei Chen Xiaoping Science and Technology Development Foundation under Grant No. CXPJJH11900001-2019308 (to J-jX).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thank Melissa Crawford, PhD, from Liwen Bianji, Edanz Editing China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.637971/full#supplementary-material
References
1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. (2018) 68:394–424. doi: 10.3322/caac.21492
2. Roayaie S, Obeidat K, Sposito C, Mariani L, Bhoori S, Pellegrinelli A, et al. Resection of hepatocellular cancer ≤ 2 cm: results from two Western centers. Hepatology. (2013) 57:1426–35. doi: 10.1002/hep.25832
3. Sapisochin G, Castells L, Dopazo C, Bilbao I, Minguez B, Lázaro JL, et al. Single HCC in cirrhotic patients: liver resection or liver transplantation? Long-term outcome according to an intention-to-treat basis. Ann Surg Oncol. (2013) 20:1194–202. doi: 10.1245/s10434-012-2655-1
4. Vitale A, Peck-Radosavljevic M, Giannini EG, Vibert E, Sieghart W, Van Poucke S, et al. Personalized treatment of patients with very early hepatocellular carcinoma. J Hepatol. (2017) 66:412–23. doi: 10.1016/j.jhep.2016.09.012
5. Villanueva A. Hepatocellular carcinoma. N Engl J Med. (2019) 380:1450–62. doi: 10.1056/NEJMra1713263
6. Méndez-Blanco C, Fondevila F, García-Palomo A, González-Gallego J, Mauriz JL. Sorafenib resistance in hepatocarcinoma: role of hypoxia-inducible factors. Exp Mol Med. (2018) 50:1–9. doi: 10.1038/s12276-018-0159-1
7. Llovet JM, Ricci S, Mazzaferro V, Hilgard P, Gane E, Blanc JF, et al. Sorafenib in advanced hepatocellular carcinoma. N Engl J Med. (2008) 359:378–90. doi: 10.1056/NEJMoa0708857
8. Cheng AL, Kang YK, Chen Z, Tsao CJ, Qin S, Kim JS, et al. Efficacy and safety of sorafenib in patients in the Asia-Pacific region with advanced hepatocellular carcinoma: a phase III randomised, double-blind, placebo-controlled trial. Lancet Oncol. (2009) 10:25–34. doi: 10.1016/S1470-2045(08)70285-7
9. Kudo M, Finn RS, Qin S, Han KH, Ikeda K, Piscaglia F, et al. Lenvatinib versus sorafenib in first-line treatment of patients with unresectable hepatocellular carcinoma: a randomised phase 3 non-inferiority trial. Lancet. (2018) 391:1163–73. doi: 10.1016/S0140-6736(18)30207-1
10. Faivre S, Raymond E, Boucher E, Douillard J, Lim HY, Kim JS, et al. Safety and efficacy of sunitinib in patients with advanced hepatocellular carcinoma: an open-label, multicentre, phase II study. Lancet Oncol. (2009) 10:794–800. doi: 10.1016/S1470-2045(09)70171-8
11. Cainap C, Qin S, Huang WT, Chung IJ, Pan H, Cheng Y, et al. Linifanib versus Sorafenib in patients with advanced hepatocellular carcinoma: results of a randomized phase III trial. J Clin Oncol. (2015) 33:172–9. doi: 10.1200/JCO.2013.54.3298
12. Xu J, Zheng L, Chen J, Sun Y, Lin H, Jin RA, et al. Increasing AR by HIF-2α inhibitor (PT-2385) overcomes the side-effects of sorafenib by suppressing hepatocellular carcinoma invasion via alteration of pSTAT3, pAKT and pERK signals. Cell Death Dis. (2017) 8:e3095. doi: 10.1038/cddis.2017.411
13. Xu J, Lin H, Li G, Sun Y, Shi L, Ma WL, et al. Sorafenib with ASC-J9(®) synergistically suppresses the HCC progression via altering the pSTAT3-CCL2/Bcl2 signals. Int J Cancer. (2017) 140:705–17. doi: 10.1002/ijc.30446
14. Xu J, Lin H, Li G, Sun Y, Chen J, Shi L, et al. The miR-367-3p increases sorafenib chemotherapy efficacy to suppress hepatocellular carcinoma metastasis through altering the androgen receptor signals. EBioMedicine. (2016) 12:55–67. doi: 10.1016/j.ebiom.2016.07.013
15. Zheng L, Xu M, Xu J, Wu K, Fang Q, Liang Y, et al. ELF3 promotes epithelial-mesenchymal transition by protecting ZEB1 from miR-141-3p-mediated silencing in hepatocellular carcinoma. Cell Death Dis. (2018) 9:387. doi: 10.1038/s41419-018-0399-y
16. Warburg O, Wind F, Negelein E. The metabolism of tumors in the body. J Gen Physiol. (1927) 8:519–30. doi: 10.1085/jgp.8.6.519
17. DeBerardinis RJ, Mancuso A, Daikhin E, Nissim I, Yudkoff M, Wehrli S, et al. Beyond aerobic glycolysis: transformed cells can engage in glutamine metabolism that exceeds the requirement for protein and nucleotide synthesis. Proc Natl Acad Sci USA. (2007) 104:19345–50. doi: 10.1073/pnas.0709747104
18. Gordan JD, Thompson CB, Simon MC. HIF and c-Myc: sibling rivals for control of cancer cell metabolism and proliferation. Cancer Cell. (2007) 12:108–13. doi: 10.1016/j.ccr.2007.07.006
19. Kroemer G, Pouyssegur J. Tumor cell metabolism: cancer's achilles' heel. Cancer Cell. (2008) 13:472–82. doi: 10.1016/j.ccr.2008.05.005
20. Xia S, Pan Y, Liang Y, Xu J, Cai X. The microenvironmental and metabolic aspects of sorafenib resistance in hepatocellular carcinoma. EBioMedicine. (2020) 51:102610. doi: 10.1016/j.ebiom.2019.102610
21. Feng J, Li J, Wu L, Yu Q, Ji J, Wu J, et al. Emerging roles and the regulation of aerobic glycolysis in hepatocellular carcinoma. J Exp Clin Cancer Res. (2020) 39:126. doi: 10.1186/s13046-020-01629-4
22. Alves AP, Mamede AC, Alves MG, Oliveira PF, Rocha SM, Botelho MF, et al. Glycolysis inhibition as a strategy for hepatocellular carcinoma treatment? Curr Cancer Drug Targets. (2019) 19:26–40. doi: 10.2174/1568009618666180430144441
23. Li X, Zhang Y, Ma W, Fu Q, Liu J, Yin G, et al. Enhanced glucose metabolism mediated by CD147 contributes to immunosuppression in hepatocellular carcinoma. Cancer Immunol Immunother. (2020) 69:535–48. doi: 10.1007/s00262-019-02457-y
24. Sun HW, Yu XJ, Wu WC, Chen J, Shi M, Zheng L, et al. GLUT1 and ASCT2 as predictors for prognosis of hepatocellular carcinoma. PLoS ONE. (2016) 11:e0168907. doi: 10.1371/journal.pone.0168907
25. Kim YH, Jeong DC, Pak K, Han ME, Kim JY, Liangwen L, et al. SLC2A2 (GLUT2) as a novel prognostic factor for hepatocellular carcinoma. Oncotarget. (2017) 8:68381–92. doi: 10.18632/oncotarget.20266
26. Yin X, Tang B, Li JH, Wang Y, Zhang L, Xie XY, et al. ID1 promotes hepatocellular carcinoma proliferation and confers chemoresistance to oxaliplatin by activating pentose phosphate pathway. J Exp Clin Cancer Res. (2017) 36:166. doi: 10.1186/s13046-017-0637-7
27. Guo Y, Li X, Sun X, Wang J, Yang X, Zhou X, et al. Combined aberrant expression of NDRG2 and LDHA predicts hepatocellular carcinoma prognosis and mediates the anti-tumor effect of gemcitabine. Int J Biol Sci. (2019) 15:1771–86. doi: 10.7150/ijbs.35094
28. Noda T, Yamamoto H, Takemasa I, Yamada D, Uemura M, Wada H, et al. PLOD2 induced under hypoxia is a novel prognostic factor for hepatocellular carcinoma after curative resection. Liver Int. (2012) 32:110–8. doi: 10.1111/j.1478-3231.2011.02619.x
29. Chen G, Wu H, Zhang L, Wei S. High glypican-1 expression is a prognostic factor for predicting a poor clinical prognosis in patients with hepatocellular carcinoma. Oncol Lett. (2020) 20:197. doi: 10.3892/ol.2020.12058
30. Huang Q, Li J, Xing J, Li W, Li H, Ke X, et al. CD147 promotes reprogramming of glucose metabolism and cell proliferation in HCC cells by inhibiting the p53-dependent signaling pathway. J Hepatol. (2014) 61:859–66. doi: 10.1016/j.jhep.2014.04.035
31. Tesori V, Piscaglia AC, Samengo D, Barba M, Bernardini C, Scatena R, et al. The multikinase inhibitor Sorafenib enhances glycolysis and synergizes with glycolysis blockade for cancer cell killing. Sci Rep. (2015) 5:9149. doi: 10.1038/srep09149
32. Fiume L, Manerba M, Vettraino M, Di Stefano G. Effect of sorafenib on the energy metabolism of hepatocellular carcinoma cells. Eur J Pharmacol. (2011) 670:39–43. doi: 10.1016/j.ejphar.2011.08.038
33. Feng J, Wu L, Ji J, Chen K, Yu Q, Zhang J, et al. PKM2 is the target of proanthocyanidin B2 during the inhibition of hepatocellular carcinoma. J Exp Clin Cancer Res. (2019) 38:204. doi: 10.1186/s13046-019-1194-z
34. Feng J, Dai W, Mao Y, Wu L, Li J, Chen K, et al. Simvastatin re-sensitizes hepatocellular carcinoma cells to sorafenib by inhibiting HIF-1α/PPAR-γ/PKM2-mediated glycolysis. J Exp Clin Cancer Res. (2020) 39:24. doi: 10.1186/s13046-020-1528-x
35. Yoo JJ, Yu SJ, Na J, Kim K, Cho YY, Lee YB, et al. Hexokinase-II inhibition synergistically augments the anti-tumor efficacy of sorafenib in hepatocellular carcinoma. Int J Mol Sci. (2019) 20:1292. doi: 10.3390/ijms20061292
36. Li S, Dai W, Mo W, Li J, Feng J, Wu L, et al. By inhibiting PFKFB3, aspirin overcomes sorafenib resistance in hepatocellular carcinoma. Int. J Cancer. (2017) 141:2571–84. doi: 10.1002/ijc.31022
37. Peng X, Chen Z, Farshidfar F, Xu X, Lorenzi PL, Wang Y, et al. Molecular characterization and clinical relevance of metabolic expression subtypes in human cancers. Cell Rep. (2018) 23:255–69.e4. doi: 10.1016/j.celrep.2018.03.077
38. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. (2010) 33:1–22. doi: 10.18637/jss.v033.i01
39. Zheng S, Zou Y, Liang JY, Xiao W, Yang A, Meng T, et al. Identification and validation of a combined hypoxia and immune index for triple-negative breast cancer. Molecular Oncol. (2020) 14:2814–33. doi: 10.1002/1878-0261.12747
40. Ji L, Lin Z, Wan Z, Xia S, Jiang S, Cen D, et al. miR-486-3p mediates hepatocellular carcinoma sorafenib resistance by targeting FGFR4 and EGFR. Cell Death Dis. (2020) 11:250. doi: 10.1038/s41419-020-2413-4
41. Lin Z, Xia S, Liang Y, Ji L, Pan Y, Jiang S, et al. LXR activation potentiates sorafenib sensitivity in HCC by activating microRNA-378a transcription. Theranostics. (2020) 10:8834–50. doi: 10.7150/thno.45158
42. Ren R, Guo J, Shi J, Tian Y, Li M, Kang H. PKM2 regulates angiogenesis of VR-EPCs through modulating glycolysis, mitochondrial fission, and fusion. J Cell Physiol. (2020) 235:6204–17. doi: 10.1002/jcp.29549
43. Yetkin-Arik B, Vogels IMC, Nowak-Sliwinska P, Weiss A, Houtkooper RH, Van Noorden CJF, et al. The role of glycolysis and mitochondrial respiration in the formation and functioning of endothelial tip cells during angiogenesis. Sci Rep. (2019) 9:12608. doi: 10.1038/s41598-019-48676-2
44. Zhang W, Zhang X, Huang S, Chen J, Ding P, Wang Q, et al. FOXM1D potentiates PKM2-mediated tumor glycolysis and angiogenesis. Mol Oncol. (2020). doi: 10.21203/rs.3.rs-40077/v1. [Epub ahead of print].
45. Wang ZH, Zhang YZ, Wang YS, Ma XX. Identification of novel cell glycolysis related gene signature predicting survival in patients with endometrial cancer. Cancer Cell Int. (2019) 19:296. doi: 10.1186/s12935-019-1001-0
46. Xu J, Wan Z, Tang M, Lin Z, Jiang S, Ji L, et al. N(6)-methyladenosine-modified CircRNA-SORE sustains sorafenib resistance in hepatocellular carcinoma by regulating β-catenin signaling. Mol Cancer. (2020) 19:163. doi: 10.1186/s12943-020-01281-8
47. Xu J, Ji L, Liang Y, Wan Z, Zheng W, Song X, et al. CircRNA-SORE mediates sorafenib resistance in hepatocellular carcinoma by stabilizing YBX1. Signal Transduct Targeted Ther. (2020) 5:298. doi: 10.1038/s41392-020-00375-5
48. Yao G, Yin J, Wang Q, Dong R, Lu J. Glypican-3 enhances reprogramming of glucose metabolism in liver cancer cells. Biomed Res Int. (2019) 2019:2560650. doi: 10.1155/2019/2560650
49. Zhou H, Wang L, Huang J, Jiang M, Zhang X, Zhang L, et al. High EGFR_1 inside-out activated inflammation-induced motility through SLC2A1-CCNB2-HMMR-KIF11-NUSAP1-PRC1-UBE2C. J Cancer. (2015) 6:519–24. doi: 10.7150/jca.11404
50. Du W, Liu N, Zhang Y, Liu X, Yang Y, Chen W, et al. PLOD2 promotes aerobic glycolysis and cell progression in colorectal cancer by upregulating HK2. Biochem Cell Biol. (2020) 98:386–95. doi: 10.1139/bcb-2019-0256
51. Yang H, Zhou L, Shi Q, Zhao Y, Lin H, Zhang M, et al. SIRT3-dependent GOT2 acetylation status affects the malate-aspartate NADH shuttle activity and pancreatic tumor growth. EMBO J. (2015) 34:1110–25. doi: 10.15252/embj.201591041
52. Wang T, Yao W, Li J, He Q, Shao Y, Huang F. Acetyl-CoA from inflammation-induced fatty acids oxidation promotes hepatic malate-aspartate shuttle activity and glycolysis. Am J Physiol Endocrinol Metab. (2018) 315:E496–510. doi: 10.1152/ajpendo.00061.2018
53. Sarapio E, De Souza SK, Model JFA, Trapp M, Da Silva RSM. Stanniocalcin-1 and−2 effects on glucose and lipid metabolism in white adipose tissue from fed and fasted rats. Can J Physiol Pharmacol. (2019) 97:916–23. doi: 10.1139/cjpp-2019-0023
54. López JJ, Jardín I, Cantonero Chamorro C, Duran ML, Tarancón Rubio MJ, Reyes Panadero M, et al. Involvement of stanniocalcins in the deregulation of glycaemia in obese mice and type 2 diabetic patients. J Cell Mol Med. (2018) 22:684–94. doi: 10.1111/jcmm.13355
55. DeWaal D, Nogueira V, Terry AR, Patra KC, Jeon SM, Guzman G, et al. Hexokinase-2 depletion inhibits glycolysis and induces oxidative phosphorylation in hepatocellular carcinoma and sensitizes to metformin. Nat Commun. (2018) 9:446. doi: 10.1038/s41467-018-04182-z
56. Martin SP, Fako V, Dang H, Dominguez DA, Khatib S, Ma L, et al. PKM2 inhibition may reverse therapeutic resistance to transarterial chemoembolization in hepatocellular carcinoma. J Exp Clin Cancer Res. (2020) 39:99. doi: 10.1186/s13046-020-01605-y
57. Feng Y, Zhang Y, Cai Y, Liu R, Lu M, Li T, et al. A20 targets PFKL and glycolysis to inhibit the progression of hepatocellular carcinoma. Cell Death Dis. (2020) 11:89. doi: 10.1038/s41419-020-2278-6
58. Zhang L, Zhang Z, Yu Z. Identification of a novel glycolysis-related gene signature for predicting metastasis and survival in patients with lung adenocarcinoma. J Transl Med. (2019) 17:423. doi: 10.1186/s12967-019-02173-2
59. Zhou Z, Huang R, Chai R, Zhou X, Hu Z, Wang W, et al. Identification of an energy metabolism-related signature associated with clinical prognosis in diffuse glioma. Aging. (2018) 10:3185–209. doi: 10.18632/aging.101625
60. Liu J, Li S, Feng G, Meng H, Nie S, Sun R, et al. Nine glycolysis-related gene signature predicting the survival of patients with endometrial adenocarcinoma. Cancer Cell Int. (2020) 20:183. doi: 10.1186/s12935-020-01264-1
61. Fang Z, He L, Jia H, Huang Q, Chen D, Zhang Z. The miR-383-LDHA axis regulates cell proliferation, invasion and glycolysis in hepatocellular cancer. Iran J Basic Med Sci. (2017) 20:187–92. doi: 10.22038/ijbms.2017.8246
62. Hua S, Liu C, Liu L, Wu D. miR-142-3p inhibits aerobic glycolysis and cell proliferation in hepatocellular carcinoma via targeting LDHA. Biochem Biophys Res Commun. (2018) 496:947–54. doi: 10.1016/j.bbrc.2018.01.112
63. Zhou Y, Huang Y, Hu K, Zhang Z, Yang J, Wang Z. HIF1A activates the transcription of lncRNA RAET1K to modulate hypoxia-induced glycolysis in hepatocellular carcinoma cells via miR-100-5p. Cell Death Dis. (2020) 11:176. doi: 10.1038/s41419-020-2366-7
64. Shang R, Pu M, Li Y, Wang D. FOXM1 regulates glycolysis in hepatocellular carcinoma by transactivating glucose transporter 1 expression. Oncol Rep. (2017) 37:2261–9. doi: 10.3892/or.2017.5472
65. Lu M, Lu L, Dong Q, Yu G, Chen J, Qin L, et al. Elevated G6PD expression contributes to migration and invasion of hepatocellular carcinoma cells by inducing epithelial-mesenchymal transition. Acta Biochim Biophys Sin. (2018) 50:370–80. doi: 10.1093/abbs/gmy009
66. Ma L, Liu W, Xu A, Ji Q, Ma Y, Tai Y, et al. Activator of thyroid and retinoid receptor increases sorafenib resistance in hepatocellular carcinoma by facilitating the Warburg effect. Cancer Sci. (2020) 111:2028–40. doi: 10.1111/cas.14412
67. Shen YC, Ou DL, Hsu C, Lin KL, Chang CY, Lin CY, et al. Activating oxidative phosphorylation by a pyruvate dehydrogenase kinase inhibitor overcomes sorafenib resistance of hepatocellular carcinoma. Br J Cancer. (2013) 108:72–81. doi: 10.1038/bjc.2012.559
68. You X, Jiang W, Lu W, Zhang H, Yu T, Tian J, et al. Metabolic reprogramming and redox adaptation in sorafenib-resistant leukemia cells: detected by untargeted metabolomics and stable isotope tracing analysis. Cancer Commun. (2019) 39:17. doi: 10.1186/s40880-019-0362-z
69. Pastorino JG, Shulga N, Hoek JB. Mitochondrial binding of hexokinase II inhibits Bax-induced cytochrome c release and apoptosis. J Biol Chem. (2002) 277:7610–8. doi: 10.1074/jbc.M109950200
70. Li S, Li J, Dai W, Zhang Q, Feng J, Wu L, et al. Genistein suppresses aerobic glycolysis and induces hepatocellular carcinoma cell death. Br J Cancer. (2017) 117:1518–28. doi: 10.1038/bjc.2017.323
71. Comerford KM, Wallace TJ, Karhausen J, Louis NA, Montalto MC, Colgan SP. Hypoxia-inducible factor-1-dependent regulation of the multidrug resistance (MDR1) gene. Cancer Res. (2002) 62:3387–94.
72. Liu R, Li Y, Tian L, Shi H, Wang J, Liang Y, et al. Gankyrin drives metabolic reprogramming to promote tumorigenesis, metastasis and drug resistance through activating β-catenin/c-Myc signaling in human hepatocellular carcinoma. Cancer Lett. (2019) 443:34–46. doi: 10.1016/j.canlet.2018.11.030
73. Marrero JA, Kudo M, Venook AP, Ye SL, Bronowicki JP, Chen XP, et al. Observational registry of sorafenib use in clinical practice across child-pugh subgroups: the GIDEON study. J Hepatol. (2016) 65:1140–7. doi: 10.1016/j.jhep.2016.07.020
74. Bruix J, Cheng AL, Meinhardt G, Nakajima K, De Sanctis Y, Llovet J. Prognostic factors and predictors of sorafenib benefit in patients with hepatocellular carcinoma: analysis of two phase III studies. J Hepatol. (2017) 67:999–1008. doi: 10.1016/j.jhep.2017.06.026
75. Casadei Gardini A, Marisi G, Scarpi E, Scartozzi M, Faloppi L, Silvestris N, et al. Effects of metformin on clinical outcome in diabetic patients with advanced HCC receiving sorafenib. Expert Opin Pharmacother. (2015) 16:2719–25. doi: 10.1517/14656566.2015.1102887
76. Miyahara K, Nouso K, Tomoda T, Kobayashi S, Hagihara H, Kuwaki K, et al. Predicting the treatment effect of sorafenib using serum angiogenesis markers in patients with hepatocellular carcinoma. J Gastroenterol Hepatol. (2011) 26:1604–11. doi: 10.1111/j.1440-1746.2011.06887.x
77. Arao T, Ueshima K, Matsumoto K, Nagai T, Kimura H, Hagiwara S, et al. FGF3/FGF4 amplification and multiple lung metastases in responders to sorafenib in hepatocellular carcinoma. Hepatology. (2013) 57:1407–15. doi: 10.1002/hep.25956
78. Vaira V, Roncalli M, Carnaghi C, Faversani A, Maggioni M, Augello C, et al. MicroRNA-425-3p predicts response to sorafenib therapy in patients with hepatocellular carcinoma. Liver Int. (2015) 35:1077–86. doi: 10.1111/liv.12636
79. Gyöngyösi B, Végh É, Járay B, Székely E, Fassan M, Bodoky G, et al. Pretreatment MicroRNA level and outcome in sorafenib-treated hepatocellular carcinoma. J Histochem Cytochem. (2014) 62:547–55. doi: 10.1369/0022155414537277
Keywords: hepatocellar carcinoma, aerobic glycolysis, Sorafenib, biomarker, prognosis
Citation: Pan Y, Hu G-y, Jiang S, Xia S-j, Maher H, Lin Z-j, Mao Q-j, Zhao J, Cai L-x, Xu Y-h, Xu J-j and Cai X-j (2021) Development of an Aerobic Glycolysis Index for Predicting the Sorafenib Sensitivity and Prognosis of Hepatocellular Carcinoma. Front. Oncol. 11:637971. doi: 10.3389/fonc.2021.637971
Received: 04 December 2020; Accepted: 15 February 2021;
Published: 18 May 2021.
Edited by:
Prathibha Ranganathan, Centre for Human Genetics (CHG), IndiaReviewed by:
Manoj Kumar Kashyap, Amity University Gurgaon, IndiaPriyanka Samji, Sri Ramachandra University, India
Copyright © 2021 Pan, Hu, Jiang, Xia, Maher, Lin, Mao, Zhao, Cai, Xu, Xu and Cai. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xiu-jun Cai, c3Jyc2hfY3hqJiN4MDAwNDA7emp1LmVkdS5jbg==; Jun-jie Xu, d2FsdGVyMjM1JiN4MDAwNDA7emp1LmVkdS5jbg==; Ying-hua Xu, Y29ubmljaCYjeDAwMDQwO3pqdS5lZHUuY24=
†These authors have contributed equally to this work