Metabolic alterations accompanying drug resistance in hepatocellular carcinoma

The development of resistance to systemic therapies in HCC patients poses a significant barrier to extending patient survival. To elucidate the mechanisms of resistance and identify potential strategies for improvement, we established lenvatinib and sorafenib-resistant HCC cell lines, Huh-7 LR and Huh-7 SR, respectively, through a prolonged drug escalation method in Huh-7 cells (Fig. 1a and Supplementary Fig. 1a, b). Compared with the parental cell line (Huh-7 P), the resistant variants presented substantially greater IC50 values for lenvatinib and sorafenib. Specifically, the IC50 value for lenvatinib in Huh-7 LRs increased 3.89-fold, and the IC50 value for sorafenib in Huh-7 SRs increased 2.46-fold, indicating stable resistance characteristics (Fig. 1a and Supplementary Fig. 1c). Further analysis revealed cross-resistance between the two cell lines, in which Huh-7 LRs exhibited resistance to sorafenib and Huh-7 SRs displayed resistance to lenvatinib (Fig. 1a and Supplementary Fig. 1c). We also assessed the sensitivity of these resistant cell lines to other clinically relevant tyrosine kinase inhibitors (regorafenib, gefitinib, and lapatinib) and chemotherapeutic agents (5-fluorouracil, irinotecan, and oxaliplatin). The results indicated a general reduction in sensitivity to these drugs (Fig. 1b, Supplementary Fig. 2, and Table 1). In 3D cultured tumor spheroids, we observed that spheroids derived from resistant cells similarly exhibited resistance to both targeted and chemotherapeutic agents, with a lower rate of cell death (Supplementary Fig. 3). Moreover, compared with spheroids derived from Huh-7 P cells, those from resistant lines demonstrated enhanced tolerance to commonly used clinical treatment regimens for HCC patients (Supplementary Fig. 4).

Fig. 1

Multidrug resistance and associated metabolic adaptations in drug-resistant cells. a Schematic representation of the drug-resistant cell model development process alongside lenvatinib and sorafenib IC50 values. The drug concentrations used for induction ranged from 1–20 μM for the Huh-7 LR cells and from 1–5 μM for the Huh-7 SR cells. Schematic figures were generated with BioRender (https://app.biorender.com/). b Drug sensitivity testing for targeted therapies and chemotherapeutic agents in drug-resistant and parental cells, presented as IC50 values. c Drug sensitivity evaluation in CDX models derived from drug-resistant and parental cells (n = 5/group). CDX model drug application concentrations: solvent: 5‰ carboxymethyl cellulose sodium; lenvatinib: 5 mg/kg/d. d Tumor growth kinetics in CDX models derived from drug-resistant and parental cells. Testing method: Unpaired Student’s t-test. e Multicolor IF staining of metabolic enzymes in tumor tissues from CDX models of drug-resistant and parental cells. GLUT1 (yellow), FABPs (purple), CD31 (red), FASN (blue), CD36 (green), and DAPI (gray). Scale bar = 100 μm. f Quantification of multicolor IF staining intensity in tumor tissues from CDX models. Testing method: Unpaired Student’s t-test. g Multicolor IF staining of metabolic enzymes in tumor tissues from HCC patients treated with or without systemic therapy. GLUT1 (yellow), FABPs (purple), CD31 (red), FASN (blue), CD36 (green), and DAPI (gray). Scale bar = 100 μm. h Quantification of multicolor IF staining intensity in patient-derived tumor tissues. Testing method: Unpaired Student’s t-test. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

Table 1 IC50 values and fold increase in drug resistance of parental vs. drug-resistant cells

When implanted into nude mice and subjected to lenvatinib intervention, tumors derived from resistant cells demonstrated heightened drug tolerance, characterized by increased tumor volumes and masses (Fig. 1c, d and Supplementary Fig. 5a, b). Consistent with the antiangiogenic properties of lenvatinib, quantitative CD31 IF demonstrated a significant reduction in microvascular density in treated parental tumors (78.87% decrease vs. solvent), whereas resistant tumors retained elevated vascularization (Fig. 1c and Supplementary Fig. 5c, d). This persistence of vascularization in resistant models aligns with their accelerated growth kinetics despite TKI administration. However, the volume of the Huh-7 LR-resistant cell-derived xenografts did not decrease, and the cells remained highly proliferative, indirectly suggesting the presence of alternative nutrient pathways supporting tumor growth (Supplementary Fig. 5e, f). Multicolor IF analysis revealed elevated expression of key enzymes associated with glucose and lipid metabolism (GLUT1 for glucose transport; FASN for fatty acid synthesis; CD36 and FABPs for fatty acid transport) in resistant tumor tissues compared with parental cells, which was unaffected by lenvatinib treatment, indicating metabolic reprogramming in resistant tumors (Fig. 1e, f). Furthermore, in tumor samples collected from certain HCC patients, metabolic enzymes related to glucose and lipid synthesis and transport were significantly upregulated in tissues from resistant patients compared with those from untreated patients and those from patients with a therapeutic response (Fig. 1g, h and Supplementary Table 1).

Single-cell transcriptome sequencing highlights increased stemness and metabolic activity in drug-resistant cells

To gain deeper insights into the molecular alterations associated with the acquisition of drug resistance in HCC cells, we performed an integrated analysis utilizing bulk RNA sequencing and single-cell transcriptomics. A total of 1,292 commonly upregulated genes were identified in two drug-resistant cell lines (Huh-7 LR and Huh-7 SR) (Fig. 2a, Supplementary Fig. 6a, and Supplementary Data 1–3). Gene Ontology (GO) biological process analysis revealed that these upregulated genes were enriched predominantly in pathways associated with cellular metabolism, exosome biogenesis, and extracellular matrix (ECM) organization (Fig. 2b). Further analysis via the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database revealed significant enrichment in metabolic pathways, ECM‒receptor interactions, drug metabolism, lysosomal processes, ATP-binding cassette (ABC) transporters, and pyruvate metabolism in drug-resistant cells (Fig. 2c and Supplementary Data 4). Gene set enrichment analysis (GSEA) additionally demonstrated that fatty acid metabolism and epithelial‒mesenchymal transition (EMT)-related processes were markedly upregulated in drug-resistant cell lines compared with their parental counterparts (Supplementary Fig. 6b).

Fig. 2

Transcriptomic analysis highlights stemness characteristics and metabolic reprogramming in drug-resistant cells. a Venn diagram depicting the upregulated genes shared between Huh-7 LR and Huh-7 SR (FC > 1.5, p < 0.05) based on bulk RNA sequencing. b PPI network analysis using ClueGo for commonly upregulated genes in the Huh-7 LR and Huh-7 SR. c KEGG analysis identified the top 20 enriched pathways of commonly upregulated genes in Huh-7 LR and Huh-7 SR cells. d Pseudotime trajectory analysis via Monocle2 and cell state clustering derived from single-cell sequencing of drug-resistant and parental cells. e Heatmap generated through GSVA via the MSigDB hallmark dataset for the 9 states derived from pseudotime analysis. f Stemness assessment of drug-resistant and parental cells from single-cell sequencing data conducted via CytoTRACE. g Clustering visualization through t-SNE dimensionality reduction from single-cell sequencing data. h Transcription factor analysis derived from single-cell sequencing of drug-resistant and parental cells. i PPI network analysis of downstream target gene clusters regulated by the transcription factors KLF4 and CEBPG. j GSEA was performed on clustering data from single-cell sequencing of drug-resistant and parental cells. k Proportional KEGG pathway enrichment analysis of upregulated genes across drug resistance-related datasets from multiple GEO database sources

Single-cell transcriptomic sequencing revealed significant heterogeneity between drug-resistant cells and their parental counterparts. The t-distributed stochastic neighbor embedding (t-SNE) dimensionality reduction analysis clearly distinguished drug-resistant cell lines from drug-sensitive ones (Supplementary Fig. 7). Pseudotime trajectory analysis, performed via the Monocle2 algorithm, classified the cells into 9 distinct states, with states 1–3 corresponding to the parental Huh-7 P cells and states 4–9 representing the drug-resistant Huh-7 LR cells (Fig. 2d and Supplementary Fig. 8). GSVA of each state revealed a gradual increase in stemness-related regulation as cells transitioned from parental to drug-resistant states. In contrast, pathways related to metabolism and cell survival initially declined, followed by peak activity at state 9 (Fig. 2e). Furthermore, CytoTRACE analysis indicated that, overall, drug-resistant cells displayed significantly higher stemness levels than did parental cells (Supplementary Fig. 9). Both the resistant and parental populations exhibited notable intratumoral heterogeneity, consisting of subpopulations with differing levels of stemness (Fig. 2f).

To characterize the molecular distinctions between drug-resistant and parental cells further, we classified the resistant cells into two clusters (clusters 0 and 1) and the parental cells into three clusters (clusters 2, 3, and 4) (Fig. 2g). Transcription factor (TF) analysis revealed that KLF4, CEBPG, SMAD3, and FOS were upregulated in Huh-7 LR cells, whereas NR1H4, CEBPD, and SOX5 were downregulated (Fig. 2h). Further investigation demonstrated that the transcription factors KLF4 and CEBPG in resistant cells primarily modulate biological processes related to metabolism, the cell cycle and differentiation, immune and inflammatory responses, and ECM remodeling (Fig. 2i, Supplementary Fig. 10, and Supplementary Data 5). Notably, KLF4 has been demonstrated to induce cancer stem cell-like phenotypes in nonstem cancer cells,27 whereas CEBPG has been implicated in promoting disease progression across various cancers.28,29,30 Conversely, TFs silenced in resistant cells, such as CEBPD and NR1H4, play critical roles in promoting immune and inflammatory gene expression while enhancing tumor cell chemosensitivity31,32,33 (Supplementary Fig. 11 and Supplementary Data 6). GSEA of the identified clusters highlighted unique biological profiles. Cluster 0 (resistant cells) displayed the upregulation of genes associated with EMT, metabolism, the cell cycle, and immune evasion pathways. Clusters 1 (resistant cells) and 2 (parental cells) were associated with pathways related to mitochondrial function and signal transduction. Clusters 3 (parental cells) and 4 (parental cells) were associated with pathways involving adherens junctions, immune and inflammatory responses, HCC-specific signaling, and ECM organization (Fig. 2j and Supplementary Fig. 12). Additionally, a multidimensional pathway analysis was conducted using nine publicly available datasets, including drug-resistant cell lines, CDX (cell line-derived xenograft) models, and HCC tissues. This analysis revealed that genes upregulated in drug-resistant cells, compared with their parental counterparts, were associated primarily with lipid metabolism, drug resistance, glutathione metabolism, and signal transduction pathways (Fig. 2k and Supplementary Data 7).

Metabolic reprogramming features of drug-resistant HCC cells and their association with resistance mechanisms

The metabolomic analysis revealed substantial variations across multiple metabolite categories, including heterocyclic compounds, lipids, nucleotides, and organic oxygen compounds, between the drug-resistant HCC cell line and its parental counterpart (Fig. 3a and Supplementary Data 8). Significant alterations in FFA metabolism, including carbon chain elongation and increased unsaturation, were evident in the drug-resistant cells (Fig. 3b and Supplementary Data 9). Pathway enrichment analysis of significantly upregulated metabolites in drug-resistant cells revealed notable alterations in pathways, including β-alanine metabolism, the Warburg effect, glutamine metabolism, de novo triglyceride synthesis, and glutathione metabolism, all of which are pivotal for cellular growth, energy production, and substrate availability (Fig. 3c and Supplementary Data 10). This study corroborated the phenomenon of aerobic glycolysis, as first described by Otto Warburg, whereby tumor cells produce high levels of lactate even under aerobic conditions, using a Seahorse extracellular flux analyzer. Compared with their parental counterparts, drug-resistant cells presented significantly increased glycolytic activity, glycolytic capacity, and glycolytic reserve, emphasizing the centrality of glycolysis as a primary energy source (Fig. 3d and Supplementary Fig. 13). The substantial increase in the levels of NAD+, a primary hydrogen acceptor involved in both glycolysis and the TCA cycle, further indicated an increased metabolic state of glycolysis and increased TCA cycle activity in the drug-resistant cells. Moreover, the observed increases in NADP+ and NADPH levels, which are linked to the pentose phosphate pathway, biosynthesis, biotransformation, and antioxidant responses, suggest enhanced energy metabolism and biosynthetic capacity in drug-resistant cells (Supplementary Fig. 14).

Fig. 3

Metabolic adaptations in drug-resistant HCC cells. a Nontargeted metabolomics profiling of parental and Huh-7 LR cells. b FFA profiling in parental and Huh-7 LR cells. c Enrichment analysis of upregulated metabolites (FC > 1.2, p < 0.05) in Huh-7 LR cells compared with parental cells. d Seahorse extracellular flux analysis of glycolytic activity (ECAR) in parental and drug-resistant cells. Testing method: Unpaired Student’s t-test. e Visualization of intracellular lipid droplet accumulation during drug resistance development. Lipid droplets (green), the cell membrane (red), and the nucleus (blue). Scale bar = 25 μm. IC50 values: Huh-7 LR_lenvatinib (S0: 13.02 μM; S1: 19.15 μM; S2: 36.40 μM; S3: 50.61 μM); Huh-7 SR_sorafenib (S0: 6.77 μM; S1: 9.46 μM; S2: 13.53 μM; S3: 19.46 μM). f Correlation analysis between the intracellular lipid droplet content and drug IC50 in drug-resistant cells. Testing method: Pearson’s correlation coefficient test. g Integrated metabolomic (four replicates per cell type, averaged in pairs), transcriptomic (two replicates per cell type), and supplementary kit-based assay profiling to map metabolic adaptations in drug-resistant cells. h High-resolution 13C-glucose metabolic flux analysis of glycolysis, the TCA cycle, and glutathione metabolism pathways. The “M+number” indicates the number of additional ¹³C atoms in the metabolite molecule. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

In the evaluation of energy storage metabolites, the triglyceride content in drug-resistant cells was significantly greater than that in parental cells, which aligns with lipid content assessments in tumor tissues (Supplementary Fig. 15a). Analysis of the lipid droplet size, number, and distribution via the neutral lipid dye Nile Red revealed significant increases in both the quantity and size of lipid droplets in drug-resistant cells compared with those in parental cells (Supplementary Fig. 15b–e). Raman spectroscopy (Raman shift: 2850 cm⁻¹), which is utilized for lipid molecule identification, further confirmed the elevated lipid content in the drug-resistant cells (Supplementary Fig. 16). Stimulated Raman scattering (SRS) vibrational imaging with a deuterium-labeled palmitic acid (D31-PA) probe was employed to observe lipid uptake, revealing a significantly stronger CDL channel signal (Raman shift: 2107 cm⁻¹) in drug-resistant cells than in parental cells (Supplementary Fig. 17). Real-time quantification of FFA uptake, conducted with an FFA tracer and high-content fluorescence imaging, corroborated the findings from SRS imaging, which revealed markedly increased fatty acid uptake in drug-resistant cells (Supplementary Fig. 18 and Supplementary Data 11). Additionally, the level of FFA β-oxidation (FAO) was significantly elevated in drug-resistant cells, indicating enhanced metabolic activity that supports survival and proliferation under drug-induced stress (Supplementary Fig. 19). Importantly, the IC50 values of lenvatinib and sorafenib enabled us to define four resistance stages (S0, S1, S2, and S3) during the establishment of the drug-resistant cell model. The quantification of lipid droplet accumulation at each stage revealed a progressive increase in lipid droplet abundance corresponding to the development of resistance. (Fig. 3e and Supplementary Fig. 20). A strong correlation was observed between the degree of resistance and the number of lipid droplets (Fig. 3f). Furthermore, exposure of parental cells to palmitic acid (PA) demonstrated that increased lipid droplet accumulation facilitates both cell proliferation and drug resistance. These findings highlight the pivotal role of lipid droplet accumulation in the progression of drug resistance (Supplementary Fig. 21).

By integrating transcriptomic data, which reflect the expression of metabolic enzymes, with metabolomic data, which capture metabolite levels, a comprehensive metabolic reprogramming map associated with hepatocellular carcinoma drug resistance was constructed. This integrative analysis revealed a high degree of concordance between the regulation of metabolic enzymes and metabolites in drug-resistant cells, with glucose, lipid, and amino acid metabolic pathways demonstrating significantly elevated activity (Fig. 3g and Supplementary Data 2, 8). Metabolic flux analysis via 13C-labeled glucose further revealed significantly increased activity across several critical metabolic pathways, including glutathione metabolism, the TCA cycle, glycolysis, the pentose phosphate pathway, and the hexosamine biosynthetic pathway (Fig. 3h, Supplementary Fig. 22–24, and Supplementary Data 12). These findings highlight the pivotal role of these metabolic pathways in meeting the elevated energy demands and supporting the survival of drug-resistant cells.

Knocking down AKR1B1 enhances drug sensitivity in drug-resistant hepatocellular carcinoma cells

Through metabolomic analysis, metabolic pathways significantly enriched during the progression of drug resistance in HCC were identified. Key rate-limiting enzymes upregulated within these pathways were further identified through the integration of transcriptomic data. As shown in Fig. 4a, AKR1B1 was identified as a prominently upregulated rate-limiting enzyme closely associated with fructose metabolism, the pentose phosphate pathway, glutathione metabolism, and triglyceride metabolism (Supplementary Data 13). Previous studies have established the critical role of the fructose metabolism pathway in tumor progression.34,35 To elucidate the underlying mechanisms, key enzymes involved in the polyol‒fructose‒lipid metabolic axis, including AKR1B1, sorbitol dehydrogenase (SORD), ketohexokinase (KHK), aldolase B (ALDOB), α-glucosidase 14 (AGPA14), and diacylglycerol O-acyltransferase 1 (DGAT1), were analyzed (Supplementary Fig. 25a). Notably, the majority of these enzymes were progressively upregulated during acquired resistance (Supplementary Fig. 25b). ALDOB, while included for completeness, presented baseline expression levels below reliable detection thresholds (FPKM < 1.0) across all the cell models. The silencing of these rate-limiting enzymes universally enhanced drug sensitivity, with AKR1B1, the most upstream enzyme in this pathway, demonstrating the most pronounced effect (Supplementary Fig. 25c, d).

Fig. 4

AKR1B1 overexpression in drug-resistant cells modulates drug sensitivity. a Integrated transcriptomic (upregulated genes, FC > 2, p < 0.05) and metabolomic profiling identified the key enzymes implicated in HCC drug resistance. b IF validation of AKR1B1 expression in parental and drug-resistant cells. AKR1B1 (red) and DAPI (blue). Scale bar = 25 μm. c WB analysis of AKR1B1 expression in cells during the development of drug resistance. d Correlation analysis of AKR1B1 expression levels with drug resistance in drug-resistant cells. Testing method: Pearson’s correlation coefficient test. e IHC staining of AKR1B1 in tumor tissues from HCC patients treated with or without systemic therapy. The results of the IHC quantification of AKR1B1 are presented in the right panel. Scale bar = 100 μm. Testing method: Unpaired Student’s t-test. f Quantification of serum AKR1B1 levels in HCC patients treated with or without systemic therapy. Testing method: Unpaired Student’s t-test. g Heatmap illustrating drug sensitivity in Huh-7 LR cells following AKR1B1 knockdown. h Drug sensitivity evaluation in CDX models derived from drug-resistant cells following AKR1B1 knockdown (n = 5/group). CDX model drug application concentrations: solvent: 5‰ carboxymethyl cellulose sodium; lenvatinib: 5 mg/kg/d. i Tumor growth kinetics in CDX models generated from drug-resistant cells following AKR1B1 knockdown. Testing method: Unpaired Student’s t-test. j Tumor weight measurements in CDX models generated from drug-resistant cells following AKR1B1 knockdown. Testing method: Unpaired Student’s t-test. k Visualization of intracellular lipid droplet levels in drug-resistant cells following AKR1B1 knockdown. Lipid droplets (green), the cell membrane (red), and the nucleus (blue). Scale bar = 25 μm. Testing method: Unpaired Student’s t-test. l WB analysis of key glutathione-regulating enzymes in drug-resistant cells following AKR1B1 knockdown via statistical analysis. Testing method: Unpaired Student’s t-test. m Quantification of intracellular GSH levels in drug-resistant cells following AKR1B1 knockdown. GSH (teal green), the cell membrane (red), and the nucleus (blue). Scale bar = 25 μm. Testing method: Unpaired Student’s t-test. n Quantification of ROS levels in drug-resistant cells following AKR1B1 knockdown. ROS (red), nuclei (blue). Lenvatinib: 50.61 μM (IC50), 12 h; sorafenib: 16.70 μM (IC50), 12 h. Scale bar = 25 μm. Testing method: Unpaired Student’s t-test. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001

As shown in Fig. 4b and Supplementary Fig. 26, AKR1B1 protein expression levels were significantly elevated in drug-resistant cells compared with their parental counterparts, with a progressive increase observed during the development of resistance (Fig. 4c). Notably, AKR1B1 expression was strongly positively correlated with the IC50 values of lenvatinib and sorafenib in Huh-7 LR and Huh-7 SR cells (Fig. 4d). To further validate these findings, tissue and blood samples from HCC patients who had undergone systemic therapy for more than 3 months and those without such treatment were analyzed. IHC and IF staining revealed substantially greater AKR1B1 expression in the treated group than in the untreated group (Fig. 4e and Supplementary Fig. 27). Moreover, the serum levels of AKR1B1 were markedly elevated in patients receiving sustained treatment compared with those in untreated individuals (Fig. 4f and Supplementary Table 2). These results suggest that AKR1B1 is upregulated during drug therapy, thereby mediating cellular survival under therapeutic stress conditions, indicating its potential utility as a biomarker for predicting therapeutic resistance in HCC patients.

Silencing AKR1B1 in drug-resistant cells substantially improved their sensitivity to the therapeutic agents responsible for inducing resistance (Supplementary Fig. 28). Heatmap analysis of multidrug resistance further demonstrated that AKR1B1 silencing significantly increased cellular sensitivity across a spectrum of targeted therapies and chemotherapeutic agents, including regorafenib, gefitinib, lapatinib, irinotecan, 5-fluorouracil, and oxaliplatin, underscoring AKR1B1’s pivotal role in mediating drug resistance phenotypes (Fig. 4g, Supplementary Fig. 29, and Supplementary Table 3). To further investigate the functional implications of AKR1B1 in vivo, tumor growth assays were performed in a nude mouse model. Compared with control tumors, tumors derived from AKR1B1-silenced cells presented significantly reduced volumes, slower growth rates, and lower tumor masses (Fig. 4h–j and Supplementary Fig. 30a). Additionally, Lenvatinib-treated tumors presented reduced vascular density, as indicated by Ki-67 staining, indicating pronounced suppression of tumor proliferation in the AKR1B1-silenced group (Supplementary Fig. 30b–d). Even in untreated tumors, proliferation was moderately inhibited, although to a comparatively lesser extent. Collectively, these findings highlight AKR1B1 as a compelling therapeutic target for enhancing the sensitivity of HCC cells to conventional treatments.

AKR1B1 enhances drug resistance in HCC cells via multiple pathways

To determine the established role of AKR1B1 in cellular drug resistance, this study further investigated the impact of altered AKR1B1 expression on the metabolic phenotype of resistant cells. Figure 4k and Supplementary Fig. 31 revealed that silencing of the AKR1B1 gene led to a marked reduction in lipid droplet accumulation and triglyceride content within resistant cells, accompanied by a significant decrease in FAO levels. Integrated transcriptomic and metabolomic profiling revealed coordinated hyperactivation of glutathione metabolism, glycolysis, and polyol/fructose metabolic pathways during acquired resistance (Fig. 3c, g, h). Among these pathways, the glutathione pathway has been shown to play a pivotal role in maintaining cell survival. Research has demonstrated that AKR1B1 modulates cystine uptake and glutathione synthesis flux in lung cancer cells via the STAT3/SLC7A11 signaling axis.26 Consistently, this study revealed similar patterns in drug-resistant HCC cells, where AKR1B1 knockdown significantly suppressed STAT3 phosphorylation and SLC7A11 protein expression (Fig. 4l). To further elucidate the role of AKR1B1 in regulating the glutathione regulatory pathway in drug-resistant HCC cells, the baseline activity of this pathway was assessed. The results revealed that the glutathione regulatory pathway was significantly upregulated in resistant cells compared with parental cells (Supplementary Fig. 32a). The quantification of glutathione levels further confirmed that compared with parental cells, drug-resistant cells presented significantly increased glutathione levels (Supplementary Fig. 32b). Concurrently, measurement of ROS levels revealed that parental cells exhibited increased sensitivity to drug treatment with increased ROS accumulation (Supplementary Fig. 32c). These findings suggest that resistant cells utilize a robust self-protective mechanism involving increased glutathione levels and decreased ROS accumulation, both of which contribute to increased cell survival. However, upon AKR1B1 knockdown, this self-protective adaptation was effectively disrupted. The glutathione regulatory pathway was significantly suppressed, as evidenced by the suppressed expression of key regulatory proteins, leading to a decrease in intracellular glutathione levels and a pronounced increase in ROS levels following drug treatment (Fig. 4l–n). Furthermore, in vivo experiments validated the in vitro findings and further confirmed that AKR1B1 plays a pivotal role in the glutathione regulatory pathway in drug-resistant HCC cells (Supplementary Fig. 33). These results highlight the pivotal regulatory role of AKR1B1 in sustaining the metabolic adaptations that underlie chemoresistance in HCC.

To further investigate the mechanisms underlying the upregulation of AKR1B1 in response to drug-induced selective pressure, we employed multiple predictive tools, including Cistrome, TCGA_LIHC, hTFtarget, and ENCODE, to identify potential transcription factors regulating AKR1B1 expression. By integrating these predictions with RNA-seq data from drug-resistant cells, we identified FOSL2, a member of the FOS transcription factor family, as a key regulator of AKR1B1 expression (Supplementary Fig. 34). This finding aligns with the results from single-cell transcriptomic sequencing (Fig. 2h). Previous studies have reported that FOSL2 can indirectly promote angiogenesis in tumor tissues, even when classical VEGF signaling is inhibited by anti-VEGF antibodies or axitinib.36 These findings suggest that FOSL2 may be adaptively upregulated in response to antiangiogenic therapies such as lenvatinib and sorafenib. In support of this hypothesis, our data demonstrated that FOSL2 expression was markedly elevated in drug-resistant cells compared with their parental counterparts. Additionally, GO enrichment analysis revealed significant upregulation of the Wnt signaling pathway, which may have contributed to the observed increase in FOSL2 expression (Supplementary Fig. 35a). We observed significant nuclear accumulation of nonphospho (active) β-catenin, a critical protein in the Wnt pathway, in drug-resistant cells (Supplementary Fig. 35b). Furthermore, the Wnt signaling inhibitor XAV9393 downregulated FOSL2 expression, underscoring the regulatory role of the Wnt pathway in modulating FOSL2 expression (Supplementary Fig. 35c). This observation is consistent with previous reports, which demonstrated the role of the Wnt pathway in promoting TKI resistance in liver cancer cells37,38,39 and enhancing FOSL2 expression in drug-resistant colorectal carcinoma cells.40 Through analysis via the JASPAR database,41 a potential binding motif for FOSL2 was predicted on the AKR1B1 promoter region: “CAGTGACTCAT” (Supplementary Fig. 36a, b). Furthermore, analysis of downstream target genes potentially regulated by FOSL2 via GO enrichment analysis revealed that these genes were associated primarily with the regulation of cellular metabolic processes (Supplementary Fig. 36c, d). These findings suggest a critical role for FOSL2 in driving the adaptive upregulation of AKR1B1 and orchestrating broader metabolic reprogramming in drug-resistant cancer cells.

AKR1B1 mediates the transmission of drug resistance and correlates with poor prognosis

These findings demonstrate that AKR1B1 expression is induced during the development of drug resistance. Moreover, patients with intrinsically elevated AKR1B1 expression tend to exhibit primary resistance to therapy. To explore this phenomenon, we collected pretreatment tumor tissues from HCC patients and stratified them into two groups on the basis of their clinical response: partial response (PR) and disease progression (DP) (Fig. 5a and Supplementary Table 1). IHC analysis of AKR1B1 expression revealed markedly higher AKR1B1 levels in tumors from DP patients than in those from PR patients (Fig. 5a, b and Supplementary Fig. 37). Consistently, proteomic data from Jiang et al. revealed markedly elevated AKR1B1 expression in the tumor tissues of HCC patients with poor prognosis compared with adjacent nontumorous tissues.42 AKR1B1 expression is progressively increased in higher-grade malignant subtypes42 (Supplementary Fig. 38). Additionally, data from David A. Wheeler’s study, which employed the iCluster algorithm to integrate information from five molecular platforms (DNA copy number variation, DNA methylation, mRNA expression, miRNA expression, and RPPA proteomics), classified HCC patients into three molecularly distinct subgroups. Among these clusters, iCluster_3 exhibited marked chromosomal instability, increased TP53 mutation rates, and widespread hypomethylation at multiple CpG loci.43 Survival analysis of TCGA data revealed that patients in this group with high AKR1B1 expression had significantly shorter survival times, highlighting the robust correlation between elevated AKR1B1 expression and unfavorable outcomes in HCC (Fig. 5c).

Fig. 5