Introduction

The most common primary malignant tumor of the liver is hepatocellular carcinoma (HCC)1. Based on the GLOBOCAN 2020 figures, hepatocellular carcinoma (HCC) was the sixth most often diagnosed cancer and the third most significant contributor to cancer-related fatalities worldwide in 2020. It was reported that there were around 906,000 new cases of HCC and 830,000 deaths2. Furthermore, there is evidence to suggest that liver cancer was the second leading cause of cancer-related deaths in China in 2022, with 316,500 deaths reported3. Despite continuous improvements in treatment modalities, including surgery, targeted therapies, and immune checkpoint inhibitors (ICI), the objective response rate in HCC patients remains at approximately 30%, and the three-year overall survival rate (OS) is still below 50%4. Hence, it is imperative to investigate novel treatment agents and mechanisms in order to enhance the prognosis of patients with hepatocellular carcinoma (HCC).

In recent decades, Phytochemicals are attracting increasing attention and focus in cancer treatment and chemoprevention due to the characteristics of easy accessibility, high safety margins, and relatively low costs associated with natural dietary compounds5. According to reports, there are more than 3000 plant species that have been found to possess anticancer properties. Thus far, around 30 substances originating from plants have been extracted and evaluated in clinical studies for the purpose of treating cancer6.

Naringenin, a flavonoid compound, is predominantly present in the form of glycosidic elements derived from the hydrolysis of naringin or narirutin. This bioactive phytochemical is abundantly found in grapefruits, lemons, oranges, and tomatoes. It is acclaimed for its pharmacological benefits, including anti-inflammatory, antioxidant, neuroprotective, hepatoprotective, and anticancer properties7. Substantial evidence underscores naringenin’s potent anticancer activity; for instance, it can significantly inhibit the migration and lung metastasis of breast cancer cells by reducing the release of TGF-β18. Additionally, naringenin impedes the progression of epithelial ovarian cancer by suppressing the expression of the PI3K pathway and fostering a healthy gut microbiota9. However, the intricate molecular mechanisms underpinning naringenin’s antitumor effects remain to be fully elucidated, particularly in the context of hepatocellular carcinoma.

The advent of high-throughput mRNA sequencing (mRNA-seq) has revolutionized our capacity to analyze gene expression changes across the entire genome, offering an unprecedented view of the transcriptomic landscape in response to therapeutic agents10. This technique has yielded valuable insights into the molecular mechanisms of drug action and has facilitated the identification of potential biomarkers and therapeutic targets. The present study aims to employ mRNA sequencing coupled with experimental validation to pinpoint the targets of naringenin in the treatment of hepatocellular carcinoma (HCC), systematically evaluating the mechanisms by which naringenin exerts its effects on HCC. Moreover, we assessed the clinical relevance and distribution of central targets in HCC and conducted cellular experiments on MHCC-97 H and HUH7 cell lines to confirm their inhibitory actions and the expression of central target genes. Ultimately, our in vitro experiments substantiated that naringenin can inhibit the proliferation and metastasis of hepatocellular carcinoma cells. This research lays the groundwork for a more profound understanding of the molecular underpinnings of liver cancer pathogenesis and the identification of novel drug targets, potentially leading to more effective treatments for patients with liver cancer.

Materials and methods

Cell culture and reagent

The human hepatocellular carcinoma (HCC) cell lines of HUH7 and HMCC-97 H were acquired from Wuhan Punosai Life Science and Technology Co., Ltd (Priscilla). To provide ideal growth conditions, the cell lines were grown in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% Fetal Bovine Serum (FBS) and 1% penicillin-streptomycin. A temperature of 37ºC was maintained for the cells in an environment consisting of 5% CO2 and 95% air. The culture media underwent systematically replaced on an alternate day.

The structure and Pharmacological properties data of naringenin

The information regarding the structure (Fig. 1) and pharmacological properties of Naringenin was acquired from the DrugBank Database and Analysis Platform11using the search term “Naringenin”. The naringin utilized in these experiments was procured from MedChemExpress, bearing the respective product identification number HY-N0100.

Fig. 1
figure 1

Molecular structure of Naringenin (molecular mass :272.25).

CCK-8 assay

The cytotoxic impact of naringenin on HUH7 and MHCC-97 H hepatic carcinoma cell lines was ascertained via the colorimetric CCK-8 assay. The study was divided into three groups: (1) experimental groups treated with naringenin at concentrations of 80 µM, 160 µM, 240 µM, and 320 µM; (2) a negative control group containing cell suspension in complete medium; and (3) a blank control group with complete medium alone (no cells or drugs). Each group comprised six replicate wells. After 24 h of initial incubation, the old medium was aspirated using a pipette. For experimental groups, 100 µL of fresh medium containing the designated naringenin concentrations was added to each well. The 96-well plate was then incubated at 37 °C under 5% CO₂ for 24, 48, or 72 h. Following incubation, the medium was removed, and 100 µL of basal medium supplemented with 10% CCK-8 reagent (Cell Counting Kit-8, Dojindo Laboratories) was carefully added to each well to avoid bubble formation. The plate was further incubated at 37 °C for 1 h. Post-incubation, the optical density at 450 nm was quantitatively assessed utilizing an ultra-microplate reader (SpectraMax; Molecular Devices Corp, Sunnyvale, CA, USA). To affirm the precision and the reproducibility of the experimental results, the assay was performed in triplicate and iterated on three separate occasions. The resultant data were normalized and presented as a percentage relative to the untreated control group.

Transwell cell migration and invasion assay

The Boyden chamber was employed for the transwell migration and invasion studies. In the upper chamber, a total of 5 × 104-1 × 105 cells were cultured in 200 µL of serum-free DMEM. Naringenin was added to the upper chamber in a gradient manner(80 µM, 160 µM, 240 µM, and 320 µM), while the lower wells were supplemented with 500 µL of medium containing 10% FBS. In the case of invasion experiments, Matrigel was applied to the bottom of the small chamber using a procedure similar to that employed for the migration experiments. After an incubation period of 20–24 h, the cells were fixed, washed with PBS, and stained with crystal violet. The assessment of hepatocellular carcinoma cell migration and invasion was conducted utilizing an inverted light microscope. The quantification of migrating or infiltrated cells was conducted by enumerating the cells within 10 individually chosen fields at a magnification of ×100.

Hoechst

Initially, cells were cultured in a 24-well plate until they attained the requisite confluency. Following the initial treatment, the cells were exposed to naringenin at graded concentrations of 80 µM, 160 µM, 240 µM, and 320 µM, and subsequently incubated for 24 h under standard culture conditions (37 °C, 5% CO₂).After incubation, cells were subjected to a phosphate-buffered saline (PBS) wash. To facilitate effective staining, 500 µL of Hoechst 33,342 solution was dispensed into each well, ensuring uniform coverage of the cellular monolayer. The cells were then incubated at ambient temperature and shielded from light for 30 min. Post-incubation, the staining solution was carefully aspirated, and the cells were washed thrice with PBS, with each wash extending for 5–10 min. The cells were subsequently filtered and imaged at 100x magnification in three arbitrarily selected fields. The process of capturing images was then followed by analyzing them using ImageJ software for further assessment and interpretation of the data.

Flow cytometry

The experimental protocol involved culturing cells until they reached 70–80% confluency. Following the initial treatment, cells were treated with naringenin at graded concentrations (80 µM, 160 µM, 240 µM, and 320 µM) and incubated under standard conditions (37 °C, 5% CO₂) for 24 h. Cells were harvested by trypsinization (using trypsin without EDTA), followed by centrifugation (300 × g, 5 min) to remove the supernatant. To ensure complete supernatant removal, cells were washed three times with ice-cold PBS and subjected to additional centrifugation (300 × g, 5 min). The cell pellet was resuspended in 500 µl of binding buffer (Annexin V Binding Buffer, BD Biosciences) for subsequent apoptosis analysis. Apoptosis was assessed using an Annexin V-FITC/PI Apoptosis Detection Kit (Sigma-Aldrich, Cat# APOAF) according to the manufacturer’s protocol. Briefly, 5 µl of Annexin V-FITC and 5 µl of propidium iodide (PI) were added to 500 µl of cell suspension, followed by gentle mixing. After incubation at room temperature for 15 min in the dark, samples were analyzed by flow cytometry to distinguish viable, early apoptotic, and late apoptotic/necrotic cell populations.

Transcriptome analysis of MHCC-97 H cells

Cells were cultivated to an optimal density and treated with naringenin at an IC50 of 182.2 µM for 24 h, designated as the NG group. Cells that did not undergo any treatment were referred to as the NC group. Following treatment, the culture medium was aspirated, and the cells were washed one to two times with 1X PBS to ensure complete removal of the medium. TRIzol reagent was then applied to lyse the cells, with pipetting performed repeatedly to achieve a clear, non-viscous lysate. The cell lysate was subsequently transferred into RNase-free cryovials with screw caps. After a 5-minute rest at room temperature, the samples were stored at −80 °C and transported on dry ice. Independent observations were conducted on all samples, with sextuplicate determinations. A differential gene expression analysis was performed to investigate the probable activities and underlying biological mechanisms of differently expressed genes. The objective was to identify candidate genes that are linked with cell proliferation, migration, and apoptosis, which are impacted by NG therapy.

Screening of differential genes

We conducted a comparative analysis of the expression profiles between the treatment group and the control group using R software (version 4.3.1). For each gene, we performed differential expression analysis using the edgeR Bioconductor software package12. We identified differentially expressed genes (DEGs) using the thresholds of |log2(fold change [FC]) | > 1 and a p-value < 0.05. In total, we identified 4498 DEGs. Heatmaps and volcano plots were generated using ggplot2 packages in R software. The same methodology was applied for the differential expression analysis of the LIHC expression profiles in the TCGA database. (https://portal.gdc.cancer.gov/)

Enrichment analysis for differential genes

Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed to elucidate the functional pathways implicated by the genes differentially expressed after pharmacological intervention using the clusterProfiler package R13. The threshold for significance in the enrichment analyses was established at a q-value and a p-value of 0.05. The results were then represented graphically as a bar plot to facilitate the visualization of the analyzed data.

Sequencing differential genes pathway map analysis

We performed visual analysis using the KEGG Mapper tool14 on the KEGG15 website (https://www.kegg.jp/) to demonstrate the complex relationships of genes that are expressed differently within key pathways. This methodology enables us to visually depict genes’ spatial arrangement and interconnectedness throughout different pathways, thus clarifying the intricacy and importance of gene regulatory networks.

Collection of Naringenin-related target genes

After conducting a search for “Naringenin” on The Chinese Medicine System Pharmacology Database and Analysis Platform(TCMSP)16we collected the corresponding targets. Due to the limited number of targets available in the TCMSP database, we supplemented our search with SwissTarget Prediction12 and DrugBank databases. The results from the three databases were combined after removing duplicates and standardized using the UniProt database17for gene names. Our analysis identified a total of 131 Naringenin-related targets.

Analysis of clinical value

We harnessed the data available in The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) to evaluate the expression levels of pivotal genes. We generated Kaplan-Meier survival curves through the Xiantao Academic platform, which interfaces with TCGA data. These curves facilitated the comparison of overall survival (OS) differences. Our analysis was thorough, with the exclusion of extraneous data. Based on the expression profiles of the core genes of interest, patients were categorized into two cohorts: high-expression and low-expression. To validate the prognostic utility of these genes, The R packages pROC and ggplot2 were utilized to generate Receiver Operating Characteristic (ROC) curves, facilitating patient diagnosis evaluation. Univariate and multivariate Cox regression analyses were applied to the TCGA data for these genes, and a prognostic nomogram was established for those genes that were achieved. The multivariate analysis yielded a statistically significant result (P-value < 0.05).

RNA extraction and quantitative real-time PCR (qRT-PCR)

TransZol Up (Cat#ET111, Beijing Full-Style Gold Bio-Technology Co., Ltd., Beijing, China) was used to extract total RNA from MHCC-97 H cells within the negative control (NC) and treatment groups (NG), following the manufacturer’s directions. The concentration and purity of the isolated RNA were quantified using a BioSpec-nano spectrophotometer (Shimadzu Co., Ltd., Japan), accepting only samples with an optical density (OD) ratio of A260/A280 between 1.8 and 2.0, which were subsequently preserved at −80 °C for later analysis. Reverse transcription was carried out with commercially sourced cDNA synthesis kits, following a protocol of incubation at 37 °C for 15 min and a subsequent inactivation at 85 °C for 5 s, with the reaction mixture thereafter being chilled to 4 °C on a specialized reverse transcription apparatus. The primer sequences employed are delineated in Table 1. The experiment required the CFX-96 Real-Time PCR Detection System (Bio-Rad Labs, Hercules, California, United States) for quantitative real-time PCR (qRT-PCR). The protocol employed a 2* SYBR Green qPCR Master Mix (Cat#SM143, Severn Innovation Biotechnology Ltd., Beijing, China). The experimental procedure involved an initial denaturation stage at 95 °C for of 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 10 s and subsequent annealing at 60 °C for 30 s. Relative quantification of gene expression was computed via the 2^-ΔΔCT method, and to ensure the reproducibility of results, all samples were analyzed in triplicate.

Table 1 Univariate and multivariate Cox regression analyses are used to assess the impact of key targets on patient prognosis.

Western blot analysis

We extracted proteins from HUH7 and MHCC-97 H cells using cold radioimmunoprecipitation assay lysis buffer (RIPA, Pierce). Total protein was quantified using the BCA assay kit (P0012, Beyotime Bio, Beijing, China). Proteins (20 µg per well) underwent separation through the utilization of sodium dodecyl sulfate-polyacrylamide gel electrophoresis, followed by their transfer onto polyvinylidene difluoride membranes. The membrane-bound proteins were blocked with rapid blocking solution for 1 h and then incubated overnight at 4 °C with primary antibodies against CHRNA7, IGFBP3, PGF, CA9, KLK1, AKR1C3 (all diluted 1:1000), and β-Actin (diluted 1:5000). The appropriate secondary antibodies (diluted 1:10000) were then added and incubated for 2 h. Membranes were imaged using a Tanon 5200 Chemiluminescent Imaging System (Tanon Science & Technology Co., Ltd., Shanghai, China) with auto-exposure settings under enhanced chemiluminescence (ECL) mode. PGL (25 kDa) and CHRNA7 (55 kDa) were analyzed on separate membranes due to size disparity, preventing signal overlap and ensuring quantification accuracy. All samples originated from the same biological source, loaded uniformly.The protein bands were quantified using the BeyoECL Plus Moon system (P0018, Beyotime Bio, Beijing, China) and analyzed with the software ImageJ (Bethesda, MD, USA) after washing the membrane twice with TBST. Western blot results represent three independent biological replicate experiments, with control and treatment samples from each replicate analyzed within the same gel, and data normalized by internal reference and combined.

Molecular docking of key targets genes

AutoDock Vina utilizes a multivariate gradient-based optimization algorithm during local refinement to enhance computational precision in ligand-receptor docking. Systematic virtual screening was performed through molecular docking simulations between six target proteins (retrieved from PDB) and ligands, with binding affinities (kcal/mol) quantifying interaction strengths - values ≤−5.0 kcal/mol indicating high-affinity conformations. Hesperidin’s 2D structure was acquired from PubChem18. The six most thermodynamically stable complexes (lowest ΔG values) underwent PyMOL-based structural analysis to map hydrogen bonds and hydrophobic interactions at binding interfaces19. This integrated approach combining quantitative scoring with structural visualization enables atomic-level characterization of molecular recognition mechanisms.

Animal experiments

All murine-oriented investigations were carried out in adherence to the guidelines set forth for the care and use of laboratory animals at Harbin Medical University, China, and subsequently sanctioned by the Institutional Animal Care and Use Committee (IACUC) of Harbin Medical University. The present research report is compliant with the stipulations proffered by the ARRIVE guidelines. Female athymic nude mice, aged four weeks, were bought from Beijing Weitonglihua Experimental Animal Technology Co., Ltd. The Harbin Medical University Ethics Committee granted approval for all experimental protocols (Permission Number: YJSDW2023-282). Each mouse received an injection into the right forelimb with 8 × 10^6 MHCC-97 H cells. Upon the tumors reaching a volume of 30 mm^3, the mice were allocated into two cohorts, each comprising four animals. The treatment group was administered naringenin intraperitoneally (i.p.) at a dosage of 50 mg/kg for a duration of 20 days, while the control group received an equivalent volume of phosphate-buffered saline (PBS) via the same route. Body weight and tumor volume were monitored biweekly post-tumor detection. Tumor volume (V) was gauged with calipers and calculated employing the formula: V = (length × width^2)/220. The study was concluded on day 35 with the humane euthanasia(Cervical dislocation method) of the mice and subsequent tumor excision for histopathological analysis.

Statistical analysis

All quantitative statistical analyses in this study were conducted using GraphPad Prism 9 (GraphPad Software, La Jolla, CA, USA). To assess diagnostic potential, receiver operating characteristic (ROC) curves were constructed, and the area under the curve (AUC) along with the 95% confidence interval (CI) were calculated. The optimal cutoff was determined by the Youden Index, where the sum of sensitivity and specificity reached its maximum. Survival analysis was performed using the Kaplan–Meier method to plot the overall survival curve, and both univariate and multivariate Cox proportional hazards models were utilized to identify prognostic factors for survival time. Moreover, to evaluate cell proliferation, migration, and invasion, three independent experiments were conducted, with the results reported as mean ± standard deviation (S.D.). For all comparative analyses, a p-value of less than 0.05 was considered statistically significant. Quantitative analysis of protein levels was carried out using ImageJ software (National Institutes of Health, Bethesda, MD, USA). Comparisons between two groups were made using an unpaired Student’s t-test. In this study, “ns” denotes a lack of statistical significance. Based on the p-value, levels of significance were deemed slight, moderate, and strong at p < 0.05, p < 0.01, and p < 0.001, respectively.

Results

Naringenin significantly inhibits the proliferation of cancer cells in vitro

To determine the inhibitory effect of naringin on hepatocellular carcinoma (HCC) and normal hepatic cells in vitro, we evaluated the cytotoxic effects of naringin on one normal cell line (WRL) and two tumor cell lines (HUH7, MHCC-97 H). For the normal cell line, low concentrations (80 µM, 160 µM) and short-term treatment (24 H) did not significantly affect cell viability. Naringin’s growth inhibition in the cancer cell lines (HUH7, MHCC-97 H) exhibited a distinct reliance on dose and time. Also, a comparison indicated that the adverse effects of naringin were more prominent in cancer cells than in normal cells (Fig. 2a-c). To evaluate the effects of naringin on the three cell lines, we constructed time-dependent and dose-dependent curves. The IC50 values at 24 h were determined by analyzing the dose-response viability curves. The IC50 values for MHCC-97 H and HUH7 were 182.2 µM and 165.7 µM, respectively. Based on the results shown in Figs. 2d and e, a concentration of 180 µM naringin induced a 24-hour treatment of the MHCC-97 H cell line in subsequent sequencing experiments.

Fig. 2
figure 2

The effects of naringenin on cell proliferation. (A) 24 h, (B) 48 h, and (C) 72 h cell viability of MHCC-97 H HUH7 WRL treated with a serial concentration of Naringenin. (C-D) Dose-dependent activity curves at 24 h. The information is presented as the average ± standard deviation of values obtained from experiments conducted in triplicate. A significance level of 0.05, 0.01, and 0.001 was observed when comparing the ctrl groups.

Naringenin inhibits cancer cell migration and invasion in vitro

The process of tumor cell migration is of paramount importance in the development and advancement of hepatocellular carcinoma. To investigate the effects of naringenin on cell migration and invasion, we utilized the scratch test and transwell experiment. Our study revealed a significant inhibition of cell migration and invasion with varying concentrations of naringenin (Fig. 3.a-h). In MHCC-97 H cells, the migration inhibitory rates were 47.4% and 57.9% for 80µM and 160µM naringenin, respectively. Similarly, in HUH7 cells, the migration inhibitory rates were 43.6% and 55% for 80µM and 160µM naringenin, respectively. The results obtained from the transwell assay were consistent with these findings, which demonstrated migration inhibitory rates of 28% and 44.5% for 80µM and 160µM naringenin in MHCC-97 H cells and 40% and 52% in HUH7 cells, respectively.

Fig. 3
figure 3

vitro naringenin has been identified to impede cell migration and invasion. A visual depiction of the outcomes obtained from the wound healing assay (A-D) and Transwell assay (E-H) was employed to assess cellular migration. The mean ± standard deviation of experimental values obtained in triplicate expresses all data. The p-value is 0.001 when compared to the control.

Naringenin induces apoptosis in cancer cells in vitro

Hoechst 33,258 staining

Hoechst 33,258 is a cell membrane-permeable fluorescent dye that can infiltrate living cells. Upon interaction with double-helical DNA, it emits a characteristic blue fluorescence. During apoptosis, nuclei stained with Hoechst 33,258 exhibit intensified fluorescence alongside morphological alterations indicative of apoptosis, including chromatin condensation and nuclear fragmentation. Our findings reveal a dose-dependent augmentation in apoptotic rates correlating with increasing concentrations of naringenin(Fig. 4a). Specifically, after 24-hour exposure to 160 µM naringenin, MHCC-97 H, and HUH7 cell lines displayed apoptotic rates of 10% and 8%, respectively, marking a statistically significant elevation from the untreated control group. The ratios of naringenin concentration increased to 30% for MHCC-97 H cells and 27% for HUH7 cells at a concentration of 320 µM (Fig. 4b-c). The data presented in this study highlight a consistent upward trend in apoptosis as naringenin concentrations increase.

Fig. 4
figure 4

Apoptosis was observed through Hoechst 33,258 staining. (A) The number of apoptotic cells increased progressively with the drug dosage. Cells treated with varying dosages demonstrated notable morphological alterations indicative of apoptosis, such as chromatin condensation and nuclear disintegration. (B, C) Statistical results of apoptosis in 97 H and HUH7 cells. Error bars represent standard deviation. Compared to the control group, *** p < 0.001.

Flow cytometry

To ascertain cellular apoptosis with higher precision and to distinguish apoptotic from necrotic cells, flow cytometry was employed. The analysis delineated a dose-dependent decrement in the population of viable cells concomitant with an increment in both the early and late stages of apoptosis (Fig. 5a-b). The incidence of mechanically induced necrosis did not exceed 4%. Upon quantification of the total apoptotic fraction, it was discerned that post 24-hour treatment with 160 µM naringenin, approximately 8.9% of the MHCC-97 H cohort underwent apoptosis, a figure that stands in significant contrast to the 3.25% observed in the control group. In parallel, HUH7 cells manifested an apoptotic rate of 6.48%, again significantly divergent from the control value of 3.7%. For cells exposed to 320 µM naringenin, the apoptotic rates were approximately 24.7% and 22.28%, respectively. Comparative analysis of the results obtained through flow cytometry and morphological examination revealed that under identical conditions, the apoptotic rate of cells treated with 160 µM naringenin for 24 h was consistent.

Fig. 5
figure 5

Apoptosis in MHCC-97 H and HUH7 cells after treatment with naringenin. (A) The flow cytometry technique was employed to detect and evaluate the various stages of apoptosis in MHCC-97 H and HUH7 cells. Necrototic, late apoptotic, early apoptotic, and healthy cells are represented by cells located in the Q1, Q2, Q3, and Q4 quadrants, respectively. (B) The ratios among healthy and apoptotic cells were evaluated using the outcomes of the apoptosis stage detection. The results are displayed as the average ± standard deviation, with each sample being measured three times separately. The standard deviation is represented by error bars. The indicated significant differences are *p < 0.05, **p < 0.01, ***p < 0.001, and #p < 0.0001 (or ****p < 0.0001).

Differential expression analysis and identification of key genes

In our quest to identify potential novel targets of naringenin in hepatocellular carcinoma, we initially conducted transcriptome sequencing and identified 4390 differentially expressed genes. Subsequently, we performed differential analysis on the LIHC dataset from TCGA, revealing 4498 differentially expressed genes (Fig. 6a-b). Following this, we employed screening techniques using The Chinese Medicine System Pharmacology Database and Analysis Platform(TCMSP), SwissTarget, and DrugBank databases, identifying 131 naringenin-related target gene. Finally, utilizing a Venn diagram21we pinpointed overlapping genes as our key targets for further investigation(Fig. 6c).

Fig. 6
figure 6

Selection of differentially expressed genes and target genes. (A) Volcano plot of sequencing results. (B) TCGA volcano plot. (C)TCGA Differentially Expressed Genes, Sequencing Differentially Expressed Genes, and Naringenin Drug Target Venn Diagram.

Functional annotation and the anti-HCC comprehensive pathway analysis of naringenin

To elucidate the regulatory mechanisms underlying the differential expression of genes and their impact on biological metabolic processes, we conducted functional enrichment analyses on the annotated genes. This approach enabled us to extract more profound functional insights. The analyses culminated in the generation of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)22 enrichment diagram(Fig. 7a-b), highlighting the top 30 enriched terms and pathways, respectively. Within the GO categories, the foremost 10 terms included response to fatty acid, positive regulation of DNA-binding transcription factor activity, cellular response to fatty acid, organic anion transport, carboxylic acid transport, organic acid transport, Wnt signaling pathway, cell-cell signaling via Wnt, regionalization, and epithelial tube morphogenesis. The leading 10 KEGG pathways were the Wnt signaling pathway, Phosphatidylinositol signaling system, Inositol phosphate metabolism, MAPK signaling pathway, Sphingolipid metabolism, Amino sugar and nucleotide sugar metabolism, Hepatocellular carcinoma, Neurotrophin signaling pathway, Human papillomavirus infection, and Axon guidance. The functional enrichment analyses have uncovered potential therapeutic targets, delineating the pathways through which naringenin may modulate hepatocellular carcinoma cells. Notably, the Wnt and MAPK signaling pathways emerged as the most prominent pathways implicated in the therapeutic effects of naringenin (Fig. 7c-d).

Fig. 7
figure 7

Enrichment Analysis of Related Genes (A, B) The top 30 significantly enriched GO and KEGG pathway terms for the 4390 intersecting target differentially expressed genes (P < 0.05). The integrated pathways map of Naringenin in an OS setting. (C) WNT signaling pathway14 (D) MAPK signaling pathway14.

Clinical value of key genes in HCC

To verify the role of key targets in hepatocellular carcinoma, based on the data in the TCGA database, we compared the expression levels of key targets in normal liver tissue and hepatocellular carcinoma. As shown in the Figure.8a, AKR1C3, PGF, CA9, and KLK1 expression levels in hepatocellular carcinoma tissue were significantly higher than in normal liver tissue. In the sequencing results, the expression levels of CA9, KLK1, and AKR1C3 were significantly reduced, further indicating that these genes are the key targets of naringenin. To assess the relationship between clinical prognosis and these key target genes, we utilized Xiantao to construct Kaplan-Meier survival curves for survival analysis. The results indicated that patients with high expression of key targets had significantly lower overall survival (OS) Compared to those with low expression of key targets (p < 0.01), suggesting a close correlation between these key targets and patient survival rate(Fig. 8b). Furthermore, to predict the 1-year, 2-year, and 3-year survival rates, we plotted time-dependent survival ROC curves for all key targets and found that all the AUC values were greater than 0.5(Fig. 8c). These findings demonstrate that the identified key targets possess effective prognostic value.

Fig. 8
figure 8

The expression levels and diagnostic value of key targets of Naringenin. (A) The expression levels of central targets in normal bone and liver tissues compared to HCC tissues. (B) Survival curves for OS (overall survival) in HCC patients with high (red) and low (blue) expression of key target genes (P < 0.05). (C) ROC (Receiver Operating Characteristic) curves confirming the accuracy of key targets expression in diagnosis and prognosis prediction (AUC = 1.000).

Fig. 9
figure 9

Nomogram for assessing the probability of HCC based on the TCGA cohort. (A) Nomograms for predicting the probabilities of 1-year, 3-year, and 5-year overall survival (OS). (B, C) Time-dependent AUC curves. (D–F) Calibration curves for the nomograms predicting the probabilities of OS at (D) 1 year, (E) 3 years, and (F) 5 years.

Cox regression analysis and nomogram

A univariate Cox regression analysis was first conducted to determine the potential independence of the six chosen main targets as predictive variables for overall survival (OS) in hepatocellular carcinoma (Table 1). The study found a notable association between IGFBP3, PGF, CA9, AKR1C3, and OS. Afterward, we conducted multivariate Cox regression analysis for the important targets that had p-values more than 0.05 in the univariate Cox regression study. Among them, CA9 and AKR1C3 were identified as independent prognostic factors for HCC. A nomogram is useful for generating individual probabilities of clinical events by integrating various prognostic and determinant variables. In our study, we constructed a nomogram based on age, gender, tumor stage, pathological stage, CA9, and AKR1C3. The results, as shown in the Figure.9, demonstrated that the nomogram had time-dependent AUC values consistently exceeding 0.6 for predicting OS within 5 years. Additionally, the nomogram exhibited relatively good consistency with observed values in predicting survival outcomes at 1, 3, and 5 years.

In conclusion, the nomogram based on CA9 and AKR1C3 can effectively predict both short-term and long-term OS in HCC. This prediction model provides valuable support for clinical diagnosis, treatment, and management.

Validation of key target gene expression by PCR

Validation of gene expression levels for 6 key target genes was performed using qRT-PCR(Fig. 10). The results revealed a general consistency between the expression trends observed in qRT-PCR and RNA-seq data, indicating a high reliability of these key target genes in transcriptional regulation within signal transduction pathways. However, it is important to note that one gene (CA9) exhibited a discordant expression trend compared to the RNA-seq results. This discrepancy may be attributed to the inherent nature of RNA-seq as a tool for large-scale screening of overall gene expression trends, which cannot guarantee a perfect concordance with qRT-PCR results for every individual gene.

Fig. 10
figure 10

The qRT-PCR technique was applied to evaluate the quantitative expression of specific important targets about the internal reference gene. The qRT-PCR data is reported as the mean ± standard deviation, with a sample size of n = 3. The statistical analysis reveals significant differences at the following levels: *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001.

Validation of key target gene expression by WB

To validate the consistency of protein expression levels of key target genes with RNA-seq results, Western blot (WB) experiments were conducted. The results demonstrated that the protein expression trends of key target genes, as verified by WB, largely agreed with the RNA-seq data, as shown in the Fig. 11. This observed concordance further supports the pivotal roles of these key target genes in signal transduction pathways, providing crucial insights for a deeper understanding of their functions and regulatory mechanisms. These findings offer valuable clues for investigating this signaling pathway’s regulatory mechanisms and potential therapeutic targets.

Fig. 11
figure 11

Representative images of the expression of central target proteins detected by Western blotting (A). Data from 3 independent experiments (B) are presented as the mean ± SD of triplicate experiments. *p < 0.05; **p < 0.01; ***p < 0.001 compared to the control group.

Molecular Docking

To further investigate the interaction between naringenin and the key targets, we utilized the CB-Dock2 online analysis tool for molecular docking. According to relevant literature, a binding energy of less than − 5 kcal/mol between a compound and a protein is considered indicative of a stable binding. The results revealed that naringenin exhibited a binding energy of less than − 5 kcal/mol with all the proteins, suggesting a strong binding affinity between naringenin and the key targets. Interestingly, the best binding effects were observed for AKR1C3 and CA9, with binding energies of −9.5 kcal/mol and − 8.5 kcal/mol, respectively, which is consistent with the results of clinical prognostic analysis. These findings further support the notion that naringenin has a strong binding affinity with the key targets, particularly AKR1C3 and CA9. The binding affinity and binding details of hesperidin with these target genes are presented in Fig. 12. The key interacting parameters between hesperidin and target amino acid residues are summarized in Table 2.

Fig. 12
figure 12

Molecular docking results for Naringenin and key targets. (A) Naringenin-AKR1C3(PDBID:1s1p). (B) Naringenin-CA9(PDB-ID:3iai). (C) Naringenin-IGFBP3(PDB-ID:7wqr). (D) Naringenin-PGF(PDB-ID:1fzv). (E) Naringenin-KLK1(PDB-ID:1spj) (F) Naringenin-CHRNA7(PDB-ID:3sq6).

Table 2 The BE of 6 key targets and Naringenin.

Validation of naringenin against HCC in animal models

To assess the in vivo anticancer activity of naringenin, a robust subcutaneous xenograft tumor model was created using the MHCC-97 H cell line (Fig. 13a). The tumor was visible ten days after inoculation. Naringenin was injected intraperitoneally once a day from the 15th day until the 35th day. During this period, there were no differences in body weight between the treated group and the untreated group. It was clearly observed that the tumor volume of the untreated group increased significantly (Fig. 13b). On the 36th day, the mice were put to death, and the tumors were removed for morphological and histological investigations. The treated group observed reduced tumor weight and volume compared to the untreated group. HE pictures revealed a lower amount of necrotic tissue in the treated group compared to the untreated group. In addition, they exhibited normal characteristics. The mice’s liver, spleen, and kidneys in both the treated and untreated groups were normal (Fig. 13c-e).

Fig. 13
figure 13

Naringenin inhibited tumor growth in vivo. Diagram of the Experimental Workflow for In Vivo Studies (A). Nude Mouse and Tumor Imaging (B). Tumor size (C). Tumor weight (D). Representative H&E is staining images of the tumor, liver, kidney, and spleen. Scale bar = 100 μm. The data are presented as the mean ± ST. The control group exhibited p-values of *P < 0.05, **P < 0.01, and ***P 0.001, respectively.

Discussion

Hepatocellular carcinoma (HCC) is a prevalent form of malignancy on a global scale, exhibiting an annual mortality rate increase of roughly 2–3%23. This alarming trend necessitates the exploration of innovative and more efficacious therapeutic approaches. Natural products represent a vast reservoir of compounds for drug discovery24among which flavonoids have garnered particular attention since the 1990 s for their anticancer properties25. Naringenin, a prevalent flavonoid, possesses substantial anticancer potential. Our study is designed to identify potential targets and pathways of naringenin in the treatment of HCC through mRNA sequencing, experimental validation, and molecular docking.

Naringenin is a naturally occurring citrus flavonoid known for its significant anticancer activity and minimal toxicity25. It exhibits less cytotoxicity against the normal liver cell line WRL compared to other cancer cell lines, a finding corroborated in ovarian cancer models9. Moreover, comparative analyses of liver, spleen, and kidney tissues from mice treated with naringenin showed no statistically significant differences, further affirming its safety profile.

Previous studies have demonstrated the efficacy of naringenin in inhibiting proliferation and inducing apoptosis in human liver cancer cell lines such as HepG226. Furthermore, it has been demonstrated to demonstrate anticancer properties by suppressing the expression of PKC and EGF, thereby inhibiting the MAPK and PI3K/Akt signaling pathways in liver cancer cells triggered by TPA27. However, the precise anticancer targets of naringenin and its clinical value remain unclear. In this work, for the first time, we utilized sequencing data combined with the extensive gene information from the TCGA database to identify a set of core target genes, including IGFBP3, PGF, CA9, KLK1, CHRNA7, and AKR1C3. These genes not only respond to naringenin but also show clinical relevance to HCC. Subsequent KEGG pathway analysis suggested that the WNT, MAPK, and P53 signaling pathways might be integral to naringenin’s anti-HCC mechanism. Our findings indicate that the expression of target genes in tumor tissues differs significantly from that in normal tissues. Moreover, PGF, CA9, KLK1, and IGFBP3 demonstrated diagnostic value for overall survival, with CA9 and AKR1C3 identified as independent prognostic factors for HCC. In addition, molecular docking validation confirmed the strong interaction between naringenin and these target genes.

Cellular experiments further validated that naringenin inhibits the viability and invasiveness of HCC cells in a time- and dose-dependent manner. Mechanistically, these effects may be mediated by the modulation of CA9, IGFBP3, AKR1C3, CHRNA7, and KLK1 by naringenin. Moreover, the suppression of CHRNA7 expression has been proven to increase the sensitivity of cancer cells to chemotherapeutic drugs28. A study conducted by Mina Afrashteh Nour et al. revealed that reducing the expression of α7nAChR greatly improves the ability of resistant HepG2 and Hepa 1–6 cells to respond to SOR treatment. This demonstrates the crucial role of α7-nAChR in determining the drug sensitivity of HCC cells29. p53 upregulates the transcription of insulin-like growth factor-binding protein 3 (IGFBP3), which acts as the main transporter of insulin-like growth factors (IGF) in the bloodstream30. Our Hoechst staining and flow cytometry results revealed a marked anti-apoptotic effect of naringenin, likely related to the upregulation of IGFBP3 and activation of the P53 signaling pathway. Aldehyde ketone reductase family member AKR1C3 plays a pivotal role in tumor progression and metastasis. Research findings have shown a notable increase in the expression of AKR1C3 in hepatocellular carcinoma (HCC), which subsequently facilitates the proliferation and spread of tumors by activating the NF-κB signaling pathway31. Intriguingly, our sequencing results revealed a marked downregulation of AKR1C3, suggesting that naringenin may exert its inhibitory effects on the NF-κB pathway through the modulation of AKR1C3 expression—a potentially novel finding that warrants further investigation.

The discoveries made in this study pave the way for the development of naringenin as a targeted therapeutic agent and provide new insights into its mechanisms of anticancer activity. Nonetheless, this research has its limitations. Although the clinical significance analysis is derived from a robust database, it remains retrospective and hypothesis-generating. Prospective clinical studies are required to confirm the therapeutic potential of naringenin and establish its efficacy and safety in the treatment of liver cancer. While our molecular docking analysis provided valuable insights into the potential interactions between naringenin and its putative targets (AKR1C3 and CA9), several limitations inherent to computational approaches must be acknowledged. First, docking simulations rely on static protein structures (e.g., PDB: 1s1p for AKR1C3), which may not fully capture the dynamic conformational changes occurring in vivo during ligand binding. Second, scoring functions used to estimate binding affinities (e.g., binding energies of −9.5 kcal/mol for AKR1C3) are approximations and may not reflect true thermodynamic stability. Furthermore, the transition from in vitro discoveries to clinical application is fraught with complexity. While our study delineates the effects of naringenin at the cellular level, in vivo dynamics, including metabolism, bioavailability, and interaction with the tumor microenvironment, require further clarification.

Conclusions

The comprehensive research findings indicate that naringenin effectively inhibits the proliferation of hepatocellular carcinoma cells in vitro and significantly impedes the growth of HCC tumors in vivo. Utilizing an array of scientific methodologies, including mRNA sequencing, experimental validation, and molecular docking, this study has substantiated the potential of naringenin as an efficacious multi-targeted agent against hepatocellular carcinoma. Notably, no adverse effects were observed in either in vitro or in vivo experiments, corroborating the inherent safety profile of naringenin as a natural constituent of dietary plants. This research provides a solid scientific rationale for the pharmacological mechanism of naringenin in the treatment of hepatocellular carcinoma and lays a foundational framework for future clinical translational studies.