Pharmacodynamic Material Basis of the Components of Four Epimedium Species with Activities Against Hepatocellular Carcinoma Based on Biological Target Networks and Multi-Omics Analysis

Introduction

The global trend in liver cancer is becoming increasingly severe. Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer. It has a high mortality rate (>90%) in the advanced stage.¹ HCC is usually caused by infection with hepatitis B/C viruses, aflatoxins, alcoholic liver disease, or non-alcoholic hepatic steatohepatitis/disease.²

Sorafenib is currently a first-line therapeutic agent for hepatocellular carcinoma (HCC).³ As a multikinase inhibitor, it primarily suppresses tumor cell proliferation and angiogenesis through inhibition of key signaling pathways such as RAF/MEK/ERK and vascular endothelial growth factor receptors (VEGFR).^4,5 In addition, sorafenib can induce tumor cell death via mechanisms including apoptosis and ferroptosis.⁶ However, despite its clinical benefits, the therapeutic efficacy of sorafenib is often limited by the development of drug resistance. Increasing evidence suggests that sorafenib resistance is a multifactorial process involving complex regulatory networks, including alterations in cell survival signaling, redox homeostasis, and tumor microenvironment adaptation. Therefore, the identification of novel therapeutic agents with distinct mechanisms of action is of great importance for improving HCC treatment outcomes. In this context, natural active compounds derived from traditional Chinese medicine have attracted considerable attention due to their diverse biological activities and potential anti-tumor effects.

Traditional Chinese medicine (TCM) formulations are important complementary alternative therapies.⁷ The characteristic components of TCM formulations, such as angelica polysaccharide⁸ and andrographolide,⁹ have demonstrated anti-HCC activity in clinical studies.

Epimedium is a genus of flowering plants in the family Berberidaceae. Epimedium species have been used as medicinal products for more than 2000 years. The main medicinal species are Epimedium brevicornu Maxim, Epimedium sagittatum (Sieb. et Zucc). Maxim., Epimedium pubescens Maxim., and the dried leaves of Epimedium koreanum Nakai. These plants are used mainly to nourish the liver and kidney, expel wind, and remove dampness. Scholars have detected >260 components in Epimedium pubescens, including flavonoids, lignans, phenolic glycosides, phenylethanol glycosides, and polysaccharides.¹⁰ Among them, components such as epimedin, desmethyl epimedin, and epimedin C have shown anti-HCC activity through various mechanisms.^11–13

At the beginning of 2022, Epimedin soft capsules were approved for marketing in China. This was heralded as a major breakthrough in the modernization of TCM in the field of liver cancer. In addition to the isoprenoid flavonoids mentioned above, Baohuoside I was found to inhibit the proliferation of HCC cells through apoptosis signaling and the nuclear factor-kappa B (NF-κB) signaling pathway.¹³

Research on the components of the four common species of Epimedium and their relationship with anticancer activity is limited. We used liquid chromatography-mass spectrometry (LC-MS), network pharmacology, bioinformatics analysis, molecular docking, and molecular-dynamics simulation for molecular prediction and validation, as well as in vitro and in vivo experiments, to predict the mechanism of action of these species.

Materials and Methods

Reagents and Materials

Four basic and original Epimedium herbs were purchased from Bozhou Traditional Chinese Medicine Commodity Trading Center and were identified as authentic by experts from Jining Medical University.

Pentobarbital sodium (ShangHai EKEAR Bio Tech Co.LTD, analytical grade), kaempferol (Rhawn Reagent, 97%), and icariin (Rhawn Reagent, 95%) were used in the study, and all other reagents were chromatographically pure.

Cells and Experimental Animals

HepG2 human hepatocellular carcinoma cells and H22 mouse hepatocellular carcinoma cells (Wuhan Pricella Biotechnology Co., Ltd., CL-0341) were used in the study.

Specific-pathogen-free (SPF) grade male BALB/c mice, 17–18 g, were purchased from Jinan Ponyue Laboratory Animal Breeding Co., Ltd., license no. SCXK (LU) 2022 0006, and Laboratory Animal Quality Certificate No. 370726241101030187. Animals were housed under controlled environmental conditions (23 ± 2 °C. 55 ± 5% relative humidity, 12h light/dark cycle) with free access to standard chow and water.

All operations performed on animals during the experiment comply with the requirements and standards of the Animal Experiment Ethics Committee, and have been officially approved by the Ethics Committee of Jining Medical College (Ethics number: JNMC2023DW106).

Toxicity Test and LC-MS of Lyophilized Powder

Freeze-dried powder of Epimedium extract was weighed and the effects of four types of Epimedium on HEPG2 cells were determined using a cell counting kit-8 (CCK-8) assay kit.

Four types of Epimedium freeze-dried powder were weighed. The test material was prepared with 80% methanol and then placed in a chromatograph. Chromatographic conditions were as follows: (1) column: AQ-C18, 150×2.1 mm, 1.8 µm, Welch; (2) flow rate: 0.30 mL/min; (3) aqueous phase: 0.1% formic acid/water solution; (4) organic phase: methanol; (5) column incubator temperature: 35°C; (6) autosampler temperature: 10.0°C; and (7) autosampler injection volume: 5.00 µL. The elution gradient is shown in Table 1. Mass spectrometry parameters are shown in Table 2. The data collected by LC-MC were searched and compared with the mzCloud database. After merging and deduplication, a common compound was identified.

Table 1 LC-MS Elution Gradient

Table 2 Mass Spectrometer

Network Pharmacology

The compounds obtained from the LC-MS were identified in PubChem based on the CAS number. The compounds were then screened using the conditions of oral bioavailability (OB) ≥ 30% and drug-like properties (DL) ≥ 0.18 for each compound in TCMSP, and corresponding targets were identified.

Targets with a predicted likelihood >0 were identified using the SwissTargetPrediction (http://swisstargetprediction.ch/) database and combined and de-duplicated with the target compounds from TCMSP. We downloaded four liver cancer-related microarray datasets from GEO^14–18 (https://ncbi.nlm.nih.gov/geo/). The details are provided in Table 3, and a volcano plot and heat maps of differentially expressed genes of single microarrays were generated using the R 4.3.3 software package limma. Finally, the four datasets were analyzed using the Robust Rank Aggregation (RRA) algorithm based on the combined multiple microarrays. GeneCards (https://www.genecards.org/) and Online Mendelian Inherited Diseases in Humans (OMIM, https://www.omim.org/) databases were searched, and targets were identified, merged, and de-duplicated with the results of the RRA analysis and then intersected with the drug targets. The intersection of the key drug ingredient targets and disease targets was visualized using VENNY 2.1 (https://bioinfogp.cnb.csic.es/tools/venny/) and then imported into Cytoscape 3.9.1 to establish a drug–ingredient–target–disease visualization network diagram. The genes related to these intersected targets were imported into the STRING database (https://string-db.org/) to construct the PPI network and to screen the core targets multiple times using CytoNCA and then were imported into Cytoscape 3.9.1 for visualization. Finally, GO and KEGG enrichment analysis was performed using R 4.3.3 software.

Table 3 All GSE Series Chip Dataset Details

HPLC, Molecular Docking, and Molecular Dynamics Simulation

High-performance liquid chromatography (HPLC) was used to detect the content of the key components in each species, subject to the limit of detection, to obtain each monomer of Icariin and kaempferol as a complex.

Molecular docking was performed using AutoDockTools (V1.5.7), PyMOL (V3.1), and AutoDock Vina (V1.2.3),¹⁹ and the binding mode with the lowest binding energy was uploaded to Plip (https://plip-tool.biotec.tu-dresden.de/plip-web/plip/index). The final pse file was exported to PyMOL for structure optimization and visualization.

Molecular dynamics simulations were conducted using Gromacs v2022.3. Small molecules were pre-processed using AmberTools22 to add general AMBER force fields (GAFFs). Gaussian 16W was used to hydrogenate the small molecules and calculate restrained electrostatic potential (RESP), which were incorporated into the molecular dynamics system topology file. The simulation was performed at 300 K and 1 bar. The force field was an AMBER99SB-ILDN, and the solvent was TIP3P water. The system’s total charge was neutralized with Na+ ions. The molecular dynamics system was first minimized by the steepest descent method. Then, canonical ensemble (NVT) and isothermal-isobaric ensemble (NPT) equilibria were performed for 100,000 steps each, with a coupling constant of 0.1 ps and a duration of 100 ps. Finally, the free molecular dynamics simulation ran for 5,000,000 steps (2 fs step size, 100 ns total duration). After simulation, built-in tools of the software package were used to analyze the trajectories to calculate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), and molecular mechanics/generalized Born surface area (MMGBSA) of each amino acid trajectory.

Inhibition of HepG2 by the Complexes

Cell Viability Assay

The CCK-8 method was used to determine the effects of different concentrations of CMPLX and monomers on cells.

HepG2 cells were inoculated into a six-well plate at a density of 300 cells per well for plate cloning experiments. A concentration gradient of CMPLX solution was used to treat the cells. The cells were incubated for 14 days, and the solution was changed every 3 days. The cells were fixed with 4% paraformaldehyde, stained with 0.1% crystal violet, and washed with PBS. Photographs of the cells were then obtained.

Next, HepG2 cells (8 × 105 cells per well) were inoculated into a six-well plate for a scratch assay. After 24 h, the cells were scraped and photographs were obtained. After incubating with different concentrations of the CMPLX solution for 24 h, photographs were again obtained.

HepG2 cells at a density of 1.5×104 were seeded onto the upper chamber of a 24-well Transwell plate (pore size 8 μm), and complete culture medium was added to the lower chamber. Different concentrations of CMPLX solution were applied to the cells. After 24 h, the solution was removed, and the cells were washed with PBS. The cells were then removed from the upper surface of the upper chamber. The remaining cells in the chamber were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet, and photographs were obtained.

Annexin V-FITC/PI Double Staining for Apoptosis Detection

HepG2 cells were digested with EDTA-free trypsin, centrifuged, stained with Annexin V-FITC/PI using a double staining kit (Solarbio, CA1020), and analyzed using a flow cytometer (BD, Canto II).

Western Blot

Total protein was extracted using a whole protein extraction kit, and the protein concentration was determined using bicinchoninic acid assay (BCA) quantification to prepare samples for electrophoresis. The proteins were electrophoresed using SDS-PAGE and transferred to an nitrocellulose (NC) membrane, blocked with 5% skim milk, and incubated with a primary antibody (1:1000) overnight in a refrigerator at 4°C. The membrane was then washed three times, incubated with a secondary antibody (1:5,000) for 2 h, and the membrane was washed. An ECL (Enhanced chemiluminescence) kit was used to develop the color. A gel imaging analyzer was used to image the membrane, and ImageJ (fiji 1.54f) software was used to analyze the images.

Animal Experiments

H22 cells in logarithmic growth were inoculated into the abdominal cavity of mice at a concentration of 2×106 cells/0.2 mL. After about 7 days, ascites was collected and injected into the subcutaneous axilla of mice (0.1 mL) to construct a tumor-bearing mouse model. Tumor-bearing mice were randomly divided into a model group (MOD), a low-dose compound administration group (CMPLXL), and a high-dose compound administration group (CMPLXH), with six mice in each group. The mice were gavaged daily, during which each group was administered the following: MOD group, 0.5% CMC Na solvent; control group, 0.2 mL per mouse; CMPLXL group, KAE+ICA at 50 mg/kg+10 mg/kg; and CMPLXH group, KAE+ICA at 100 mg/kg+20 mg/kg. Starting from the third day of inoculation, the tumor volume (TV) was measured every two days using the following equation: TV (mm3)=1/2ab2, where a is the longitudinal diameter and b is the short diameter.

Twelve hours after the last administration, blood was collected from the orbit under deep anesthesia with the intraperitoneal injection of 0.3% pentobarbital sodium (30mg/kg). The serum was separated by centrifugation and stored at −80°C. Euthanasia was confirmed by verifying cessation of respiration and cardiac activity, along with dilated pupils persisting for at least 10 minutes, following cervical dislocation in mice. And the contents of the colon together with tumor tissue were excised. The colon contents were stored at −80°C, and the tumors were stored in 4% paraformaldehyde. After embedding and slicing, the tumor tissue was stained with hematoxylin and eosin (H&E), fixed, and observed under a microscope.

Non-Targeted Serum Metabolomics

Briefly, 200 μL of pre-cooled methanol/acetonitrile (1:1, v/v) was added to 50 μL of a serum sample, mixed well, concentrated, and dried. Then, 100 μL of 50% methanol (containing 5 ppm 2-chlorophenylalanine) was added to reconstitute the sample, and the sample solution was filtered through a 0.22-μm filter membrane to obtain the test sample. Then, 10–20 μL of the filtrate of each sample was used to form a QC sample, and the remainder of the test sample was used for LC-MS detection.

The sample was analyzed on an ACQUITY UPLC HSS T3 column (100 Å, 1.8 µm, 2.1 mm × 100 mm) at a flow rate of 0.4 mL/min, column temperature of 40°C, automatic sampler temperature of 8°C, and a sample volume of 2 μL. The positive and negative mobile phases were a 0.1% formic acid aqueous solution (mobile phase A) and acetonitrile containing 0.1% formic acid (mobile phase B). The elution gradient is shown in Table 4.

Table 4 Elution Gradient

A Thermo Orbitrap Exploris 120 mass spectrometer, paired with an Xcalibur (Thermo, V4.7), was used to acquire data. The HESI source settings were as follows: a spray voltage of 3.5 kV/−3.0 kV, sheath gas of 40 arb, auxiliary gas of 15 arb, capillary temperature of 325°C, and auxiliary gas temperature of 300°C. The primary resolution was 60,000 with a scan range of 100–1000 m/z, a standard automatic gain control (AGC) target, and a maximum injection time (IT) of 100 ms. The top four ions screened were secondarily fragmented. The dynamic exclusion time was 8 s, with a secondary resolution of 15,000, higher-energy collision dissociation (HCD) collision energy of 30%, standard AGC target, and automatic maximum IT. Both formal and QC samples were processed using the same chromatographic and mass spectrometry methods.

The data were imported into the commercial software Compound Discoverer™ 3.3 (V3.3.2.31, Thermo, Waltham, USA). The undetected peaks were filtered from >50% of the QC samples and filled in with missing values based on the software’s Fill Gaps algorithm. The sum of the total peak area was normalized. Metabolite identification was based on self-built libraries, the mzCloud online library (https://www.mzcloud.org/), LIPID MAPS (https://www.lipidmaps.org/), HMDB (https://hmdb.ca/), MoNA (https://mona.fiehnlab.ucdavis.edu/), and the NIST_2020_MSMS spectral library using an MS1 mass tolerance of 15 ppm and a MS2 match factor threshold of 50.

Using PMCMRplus (V1.9.6) to perform two sets of ANOVA and multiple sets of MACANOVA on the data, the p-value was calculated. Principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) were performed on the sample data using the R package Ropls, and overfitting testing of the model was performed using a permutation test. The Pheatmap package (V1.0.12) was used to perform cluster analysis on the abundance values of differential metabolites and to draw heatmaps and perform trend analysis. Venn diagrams were drawn using VennDiagram (V1.7.3). Correlation analysis on differential metabolites was performed using corrplot (V4.0.3). The correlation between differential metabolite products was analyzed using machine learning (mlr3verse, V0.2.7) and ROC curves (pROC, V1.18.2). KEGG enrichment analysis was performed on differential substances using clusterProfiler (V4.6.0).

Intestinal Flora Histology

We conducted 16S rRNA sequencing and data analysis of the gut microbiota with total DNA isolated using a fecal DNA isolation kit (D5635-02) according to the manufacturer’s directions. PCR amplification (ABI, 2720) was performed using NEB Q5 DNA high fidelity polymerase with the sequencing forward primer ACTCCTACGGGGAGGCCA and reverse primer GGACTACHVGGGWTCTAAT. The sequencing region was 16S-V3V4. The Greengenes2 database was used to classify each amplicon sequence variant (ASV) by comparing reference sequences. QIIME2 software was used to analyze the alpha and beta diversity indices. UniFrac was used to measure the beta diversity distances to study changes in the microbial community structure between samples.

Correlation Analysis of Differential Metabolites with Differential Gut Microbiota

The R package vegan (2.2.6.1) was used to calculate the Bray–Curtis distance matrix of the non-targeted serum metabolomics data and microbiota abundance data, and then a Mantel statistical analysis test was performed using QIIME software to calculate the p-value. Mothur (1.35.1) software was used to calculate the Spearman rank correlation coefficient between non-targeted serum metabolomics data and colony abundances, and the results were imported into Cytoscape software for visualization.

Data were analyzed using GraphPad Prism (V10.1.0) software. Experimental data were expressed as mean ± SD. The statistical significance of the results was analyzed using one-way analysis of variance (ANOVA) for multiple comparisons and t tests for two-group comparisons. A p-value < 0.05 indicated a significant difference. All experiments were conducted in at least triplicate.

Results

Cytotoxic Effects of the Extract Powders of Four Epimedium Species on Human Liver Cancer (HepG2) Cells and Determination of Their Shared Compounds

The survival rate of HCC cells decreased gradually, and toxicity increased with increasing concentrations of the lyophilized powder extracts of each Epimedium species within each concentration gradient (Figure 1A). Meanwhile, the half-maximal inhibitory concentration (IC50) values were calculated (Figure 1B).

Figure 1 Toxicity of the extract powders of four Epimedium species against HepG2 cells and shared compounds. (A) Toxicity of the extract powders of four Epimedium species against HepG2 cells and (B) IC50 (**P<0.01, ***P<0.001, ****P<0.0001).

The LC-MS-based identification of total ion flow is shown in Figure S1. The matching LC-MS results are shown in Table 5.

Table 5 LC-MS Results Matching Table

Mechanism of Network Pharmacology-Based Shared Compounds Against Liver Cancer

The Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCSMP) database and PubChem database, combined with LC-MS, were used to establish a database of 67 common components of the four Epimedium species. Ten common components were obtained after correcting the OB and DL values of the TCMSP database (Table 6). A total of 219 targets were jointly predicted and collected by the TCMSP and SwissTargetPrediction database. The differential genes obtained by RRA multi-chip co-analyses of the four Gene Expression Omnibus (GEO) liver cancer chips were integrated with the genes related to liver cancer in the GeneCards and Online Mendelian Inheritance in Man (OMIM) databases, totaling 17,318 targets. Heat maps and volcano maps of four GEO chips were showed in Figures S2 and S3A, and the heat map of RRA was showed in Figure S3B. A total of 215 intersecting targets were taken after intersecting the above-mentioned targets for CMPLX and liver cancer disease using Venny 2.1.0 (Figure 2A). As shown in Figure 2B, we constructed a network graph containing 227 nodes (10 components, 124 targets) and 640 lines to elucidate the relationships among them. After integrating the STRING database, the protein–protein interaction network (PPI) was imported into Cytoscape 3.9.1 for analysis. The top four targets by degree value were tumor protein 53 (TP53), heat shock protein 90 alpha family class A member 1 (HSP90AA1), estrogen receptor 1 (ESR1), and protein kinase B (AKT) (Figure 2C).

Table 6 Components

Figure 2 Mechanistic prediction of network pharmacology-based shared compounds against hepatocellular carcinoma. (A) Venn diagram showing intersections between genes. (B) Drug–component–target–disease network diagram. (C) PPI network map and screening of key targets. (D) Enrichment analyses using the GO database. (E) Enrichment analyses using the KEGG database.

Enrichment analyses of 198 potential therapeutic targets based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were done using R (R Institute for Statistical Computing, Vienna, Austria) and related software packages. Functional enrichment analyses using the GO database are shown in Figure 2D. Among cellular components, 95 entries were enriched, mainly involving “membrane rafts”, “membrane microstructural domains”, and “neuronal cell bodies”. With regard to molecular function, 208 entries were enriched, mainly involving “DNA-binding transcription factor binding”, “RNA polymerase II-specific DNA-binding transcription factor binding”, and “protein tyrosine kinase activity”. With respect to signaling-pathway enrichment using the KEGG database, 143 signaling pathways were obtained. According to the top-30 significantly enriched pathways sorted by q-value, the targets were mainly enriched in “neuroactive ligand-receptor interaction”, “PI3K-Akt signaling pathway”, “chemical carcinogenesis receptor activation”, “lipids and atherosclerosis”, and the “calcium signaling pathway” (Figure 2E).

Shared Component Complexes and Their Molecular Dynamics Simulations

As determined by high-performance liquid chromatography (HPLC), only kaempferol and icariin reached the limit of detection, and their contents are shown in Table 7. Molecular docking was undertaken between icariin, kaempferol, and the core targets of four PPI networks (AKT1, ESR1, TP53, HSP90AA1) in Figure 3A. The lowest binding energy between each component and the target was < −5, indicating “good” affinity and “strong” binding activity, which preliminarily confirmed the reliability of the network pharmacology results stated above. Among them, AKT1 interacted with icariin and kaempferol through hydrophobic interactions, hydrogen bonds, salt bridges, π–cation interactions, and π-stacking (Figure 3B).

Table 7 Compounds Arriving at the Detection Line

Figure 3 Screening and molecular docking of PPI core targets. (A) Heatmap of the mean binding energy of PPI core targets with icariin and kaempferol. (B) Molecular docking of AKT1 with icariin and kaempferol. (C) RMSD values of icarin (blue), AKT1 (red) and COMPLEX (grey). (D) RMSD values of kaempferol (blue), AKT1 (red) and COMPLEX (grey). (E) RMSF values, (F) Gyration values, (G) SASA values and (H) Number of hydrogen bonds of AKT1-icariin (blue) and AKT1-kaempferol (red). (I) AKT1–icariin residue free energy decomposition. (J) AKT1–kaempferol residue free energy decomposition.

Leveraging molecular dynamics simulations, we employed GROningen MAchine for Chemical Simulations (GROMACS) to characterize and understand the interaction mechanism between AKT1 and icariin and kaempferol. The RMSD curves of AKT1–icariin/kaempferol complexes showed equilibrium after 70 and 60 ns, with RMSD values of 0.26± 0.02 and 0.27±0.02 nm (Figure 3C and D). RMSF analysis identified key flexible regions in the AKT1 protein, particularly around amino-acid residues 100–200 versus 300 of the AKT1–icariin/kaempferol complex (Figure 3E). The radius of gyration (Rg) values of both complexes remained balanced at 2.41 ± 0.02 nm (Figure 3F). The corresponding solvent-accessible surface area (SASA) values remained stable approximately 240 and 238 nm2 (Figure 3G). During the 100-ns simulation, the number of hydrogen bonds of AKT1–icariin/kaempferol complexes was 0–8 and 0–6, respectively (Figure 3H). Moreover, icariin and kaempferol interacted with the key residues of AKT1, which are essential for catalysis, and the residue free energy decomposition is shown in Figure 3I and J. The data for binding free energy for the AKT1–icariin/kaempferol complexes are shown in Table 8.

Table 8 Combined with Free Energy Calculations

Complexes Inhibit the Proliferation and Migration of HepG2 Cells and Induce Apoptosis

Among monomer-delivered groups, the inhibition of HepG2 cells by kaempferol was more significant, and the inhibition rate was close to the plateau when the concentration reached 150 μM (Figure S4A). However, administration of the icariin monomer seemed to alter the possibility of inhibiting the proliferation of HepG2 cells (Figure S4B). By contrast, after administration of the two complexes, HepG2 cells were exposed to different CMPLX concentrations at fixed concentrations of 75 and 150 μM of kaempferol. Low concentrations of icariin counteracted the inhibitory effect of kaempferol, whereas increasing icariin concentrations led to a slight decrease in cell viability. Compared with kaempferol alone, low concentrations of icariin interfered with CMPLX formation, lacking significant anti-HCC activity. However, high concentrations of icariin in CMPLX exhibited significant anti-HCC activity. Administration of the kaempferol monomer with an icariin concentration in CMPLX up to 30 μM had a more pronounced inhibitory effect than administration of a single identical concentration of kaempferol (Figure 4A and B).

Figure 4 Complexes inhibit the proliferation and migration and induce apoptosis of HepG2 cells. (A) Toxicity of kaempferol (75 μM). (B) Toxicity of kaempferol (150 μM). (C) Wound-healing assay. (D) Cell-migration assay. (E) Plate clone-formation assay. (F) Annexin V-FITC/PI double staining for apoptosis detection. (G) Western blots. (ns, not significant, *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001 vs control, #P<0.05 was 150KAE+30ICA vs 150KAE).

Experiments based on wound healing, cell migration, and plate clone-formation assays showed that CMPLX could suppress the growth of HepG2 cells (Figure 4C–E). Even low-dose CMPLX retarded wound healing (P<0.0001), and high-dose CMPLX led to a reduction in cell adhesion. Low and high concentrations of CMPLX reduced cell migration (P<0.05, P<0.001). Low and high doses of CMPLX decreased the clone-formation ability of HepG2 cells (P<0.0001).

CMPLX treatment induced significant apoptosis of HepG2 cells, with total apoptosis percentages of 85.65% ± 2.46% (P<0.0001) and 93.15% ± 2.03% (P<0.0001) at low and high concentrations, respectively, compared with the control group (Figure 4F).

Western blots are shown in Figure 4G. Compared with the control group, the CMPLX-treated group exhibited downregulation of p-Akt protein expression (P<0.001) and p-PI3K protein expression (P<0.001) within the phosphoinositide 3-kinase protein kinase B (PI3K-Akt) signaling pathway, thereby inhibiting the growth of HepG2 cells. CMPLX significantly decreased Bcl-2 protein expression (P<0.0001) and induced the apoptosis of HepG2 cells.

Comparison of Serum Metabolite Differences in H22 Tumor-Bearing Mice with Different Syndromes

Changes in tumor volume in mice were determined (Figure S5A). Tumor growth was inhibited under CMPLX treatment. Changes in tumor volume and tumor growth were significantly different from those in the MOD group and were more significant at a high dose of CMPLX. Figure S5B shows the H&E staining of tumor tissues. In the high-dose group of CMPLX, more necrotic fragments and pathological dendrites were observed.

We wished to elucidate how kaempferol in combination with icariin induced anti-HCC effects. We undertook non-targeted serum metabolomics analysis on experimental mice. The overall principal components analysis (PCA) of the experimental quality control (QC) under positive and negative ions (Figure 5A) was good. QC correlation analyses showed the experimental reproducibility to be good (Figure S6A). To explore changes in levels of metabolites in MOD, CMPLXL, and CMPLXH groups, we used orthogonal partial least squares discriminant analysis (OPLS-DA) to undertake comprehensive classification and discriminant analysis. Multivariate statistical OPLS-DA and partial least squares discrimination analysis (PLS-DA) were used to show differences in scores in positive and negative ion modes for two-by-two comparisons of MOD and CMPLXL and CMPLXH groups, and three-by-three comparisons of MOD, CMPLXL, and CMPLXH groups, respectively (Figure 5B). A permutation test confirmed the validity of the model (Figure S6B). Identification of differential metabolites was based on significant P-values (P < 0.05) for variable importance in projection (VIP) ≥1.0, and Student’s t-test S-plots are shown in Figure S6C.

Figure 5 Differences in serum metabolites in H22 hormonal mice. (A) Experimental QC overall PCA, it includes MOD (yellow), CMPLXL (blue), CMPLXH (red) and QC (green). (B) OPLS-DA of MOD (blue) vs CMPLXL (red) and MOD (blue) vs CMPLXH (red) in positive and negative ionization modes. (C) Intersecting differential metabolites in MOD vs CMPLXL and MOD vs CMPLXH. (D) Enrichment analyses using the KEGG database (Human Disease in red, Metabolism in green and Organismal System in blue). (E) PLS-DA of MOD (yellow) vs CMPLXL (blue) vs CMPLXH (red) in positive and negative ionization modes. (F) The random forest analysis of the top-30 differential metabolites by importance.

We wished to identify HCC-related variables affected by CMPLX. We analyzed the intersections on Venn diagrams for the differential metabolites of MOD_vs_CMPLXL and MOD_vs_CMPLXH groups, and 84 proprietary differential variables were identified (Figure 5C). This exclusivity suggested that these variables were uniquely affected by the synergistic effects of kaempferol and icariin, making them potentially relevant to anti-HCC effects. We wished to reveal the metabolic pathways related to the anti-HCC effects of CMPLX. We conducted analyses using the KEGG database on the 84 intersecting differential variables mentioned above. Results showed that these metabolites were involved in “central carbon metabolism in cancer”, “renal cell carcinoma”, “signaling pathways in cancer”, “cysteine and methionine metabolism”, the “citric acid cycle” (TCA cycle), and “riboflavin metabolism”. Multiple metabolic pathways associated with disease progression were dysregulated in mice treated with CMPLX (Figure 5D).

In addition, we investigated trends in the relative content of relevant differential metabolites using the three groups’ joint comparative differential substance machine learning analysis. The PLS-DA depicted in Figure 5E, which was conducted in positive and negative ion modes, revealed that the two complex-treated groups were significantly different from the model group and from each other, as indicated by the clear separation in the score plot. The corresponding replacement test plots are shown in Figure S6D.

Under the data processing method of Merge mode with Random Forest analysis and standardized score (Z_Score) and averaging by samples within groups, the results are shown in (Figure 5F). Statistical analyses of absolute abundance values showed that 17 of the top-30 differential metabolites by importance were significantly upregulated in the CMPLXL group versus the CMPLXH group compared with the MOD group, and 13 of these differences were highly significant (P < 0.0001) in both comparison groups.

Dysbiosis of the Gut Microbiota in CMPLX Improves H22 Hormonal Mice

Upon categorization, 1130, 960, and 1498 amplicon sequence variants (ASVs) were obtained from MOD, CMPLXL, and CMPLXH mice, respectively.

We assessed the alpha diversity (Figure 6A) and beta diversity between groups. The Bray–Curtis algorithm under principal coordinate analysis (PCoA) plots showed separation for each of the two groups (Figure 6B). The beta diversity of intestinal bacterial flora was also analyzed using this algorithm under non-metric multidimensional scaling (NMDS) (Figure 6C). Then, we mapped the composition of gut microbiota species in each group using hierarchical clustering analysis (Figure 6D).

Figure 6 CMPLX ameliorates gut microbiota dysbiosis in H22 hormonal mice. (A) Individual indices of alpha diversity analysis (MOD in Orange, CMPLXL in yellow and CMPLXH in blue). (B) PCoA and (C) NMDS analysis plot of MOD (blue), CMPLXL (yellow) and CMPLXH (red). (D) Hierarchical clustering analysis in family level and genus level. (E) Venn diagram of the three groups. (F) Differential analysis of MOD vs CMPLXL (P<0.001 in orange) and MOD vs CMPLXH (P<0.05 in orange, P<0.01 in purple and P<0.001 in yellow) metabolic pathways. (H) LEfSe plot. (G) Metabolic pathway statistics (*P<0.05, **P<0.01).

The microbiomes of the three groups were dominated by the families Muribaculaceae, Lachnospiraceae, Marinifilaceae, and Bacteroidaceae. In the MOD group of mice, the abundance of Enterobacteriaceae_A was significantly increased relative to that in the CMPLXL and CMPLXH groups. The abundance of Marinifilaceae in the CMPLXH group significantly increased compared with that in the CMPLXL and MOD groups. At the genus level, the abundance of Escherichia species in the MOD group was not only higher than that in the treated groups but also accounted for most of the four samples comprising the group. In the MOD and CMPLXL groups, the abundance of Odoribacter species and Kineothrix species was significantly lower than that in the CMPLXH group. Species differences and marker species analysis of ASV Wayne plots are shown in Figure 6E.

The group with downregulated expression was designated as “MOD”. The groups with upregulated expression were designated as CMPLXL and CMPLXH. These classifications were made for pairwise comparisons of differential metabolic pathways between them according to the KEGG database. Compared with the CMPLXL group, the MOD group had only one metabolic pathway, “β-alanine metabolism” (P<0.05). Compared with the CMPLXH group, the MOD group showed higher enrichment in “bacterial invasion of epithelial cells”, “biosynthesis of penicillin and cephalosporins”, and “biosynthesis of non ribosomal peptide iron carriers” (P<0.05). The CMPLXH group had the highest enrichment in “β-lactam antibiotic resistance” (Figure 6F). The CMPLXL and CMPLXH groups were subjected to metabolic-pathway analyses using the KEGG database. The enrichment of “cell growth and death” and “cell motility in cellular processes” was significant. “Environmental information processing”, “enrichment of membrane transport in environmental information processing”, “folding, sorting and degradation”, “replication and repair”, and “translation in genetic information processing” showed high activity. “Carbohydrate metabolism”, “cofactor metabolism and vitamins”, and “amino acid metabolism” were the top-three enrichments (Figure 6G).

Linear discriminant analysis effect size (LEfSe) and random forest analyses were used to distinguish specific microorganisms in different taxonomic groups. Evolutionary maps and linear discriminant analysis (LDA) histograms were plotted (Figure 6H). In the taxonomic branching diagram analyzed by LEfSe, the MOD group exhibited unique evolutionary branching among the particularly enriched microbial communities evolving from phylum to genus levels. The main enrichments were in the abundance of Proteobacteria (phylum), Gammaproteobacteria (class), Enterobacteriales A (order), Enterobacteriaceae A (family), and Escherichia (genus), with LDA scores (log10) >4 for all five nodes. In the CMPLXL group, three unique evolutionary branches were formed at the phylum level: Actinobacteriota, Desulfobacterota_i, and Patescibacteria. In the CMPLXH group, the abundance of Firmicutes A was significantly different from that in the other two groups at the phylum level, and the abundance of Christensennelales, Lachnospirales, and Oscillospirales was also significantly increased at the order level compared with that in the other two groups. Random forest analyses also revealed this information in the form of heatmaps.

Correlation of Differential Metabolites with Differential Gut Microbiota Under CMPLX Treatment

Mantel correlation test for the two sets of data correlation data as in Table 9.

Table 9 Mantel Correlation Test

As shown in the heatmap (Figure 7), luteolin derivatives (luteolin 4′-sulfate, luteolin 7-O-diglucuronide) and kaempferol derivatives (kaempferol 3-glucuronide, kaempferol 3-glucuronide-7-sulfate) showed a significant negative correlation with Escherichia species, Enterococcus_H, and Mammaliicoccus species.

Figure 7 Heatmap of correlation between differential metabolites and differential gut microbiota (*P<0.05, **P<0.01).

Discussion

Recently, there has been a growing emphasis on the resistance and adverse effects of primary medications for liver cancer, particularly sorafenib. TCM, which serves as an auxiliary treatment, has often been overlooked. In 2022, the approval and marketing of the innovative Chinese medicine Epimedin soft capsules was undertaken.

We have shown that baohuoside 1, isolated from Korean Epimedium species, can inhibit the proliferation of HCC cells through apoptosis signaling and the NF-κB pathway. This is only one variety of Epimedium species, and for the other three species, there may be one or more types of common constituents associated with anti-HCC effects,¹³ a question that we attempted to answer in the present study.

We found that all varieties of the lyophilized powder of an Epimedium herb extract exhibited excellent anti-HCC activity. After identifying the shared components by LC-MS and using network pharmacology, 10 shared components with anti-HCC potential were identified. Only kaempferol and icariin reached the limit of detection after LC of these compounds.

Recently, kaempferol has been shown to induce apoptosis and autophagy in Hep3B cells by downregulating the expression of Bcl-2 and upregulating the expression of BCL2-associated X (Bax), BH3 interacting domain death agonist (Bid), caspase-3, beclin-1, and microtubule-associated protein light chain 3 (LC3).²⁰ Kaempferol can inhibit colorectal cancer via circ_000034-mediated jumonji domain containing 1C (JMJD2C)/β-catenin. Icariin possesses various immune functions. The potential molecular mechanisms may include T-helper 17 (Th17)/regulatory T lymphocyte (Treg) regulation and proliferation of NK cells, followed by anti-HCC activity.²¹ Icariin has shown strong synergistic effects in combination with arsenic trioxide²² and hydroxypropyl-gamma-cyclodextrin.²³ Notably, in HCC, enrichment analyses using the KEGG database revealed “lipid and atherosclerosis” and “PI3K-Akt” to be central pathways in cellular signaling. Moreover, enrichment analyses using the GO database indicated “modulation of serine/threonine kinase activity”, “phosphorylation of peptidylserine residues”, and “serine modification of peptide chains” to be enriched.

CCK-8 data have shown that the activity of HepG2 cells cannot be inhibited under the administration of the icariin monomer. We co-administered kaempferol and icariin. With the intervention of different concentrations of icariin, co-administration exhibited superior inhibitory effects on HCC compared with that employing administration of kaempferol monomer. CMPLX (icariin and kaempferol) exhibited significant anti-HCC activity according to experiments based on wound healing, cell migration, and clone formation. To investigate whether CMPLX induced apoptosis, we used Annexin V-FITC/PI double-staining to detect cells. Percent apoptosis at a 75+15 μM CMPLX concentration was 85.65%±2.46%. Percent apoptosis at a 150+30 μM CMPLX concentration reached 93.15%±2.03%. Western blotting illustrated that CMPLX blunted growth and promoted apoptosis of HepG2 cells by inhibiting the protein expression of p-Akt, p-PI3K, and Bcl-2 in the PI3K-Akt signaling pathway.

Studies in recent years have linked hepatitis, HCC, non-alcoholic fatty liver disease, and alcoholic liver disease to alterations in the structure and function of the gut microbiota. Ecological dysregulation of gut symbiotic communities is known to be associated with immune responses, including those affecting the liver.²⁴ In the present study, 16S rRNA gene sequencing was used to analyze the abundance and intergroup differences of gut microbiota in mice. Hierarchical clustering analysis showed that Muribaculaceae, Lachnospiraceae, Marinifilaceae, and Bacteroidaceae were dominant families across the three groups. Notably, the abundance of Enterobacteriaceae_A was higher in the MOD group, whereas Marinifilaceae was more abundant in the CMPLXH group compared with the CMPLXL and MOD groups.

Previous studies have reported that an increased abundance of Lachnospiraceae is associated with hematopoietic recovery and gastrointestinal repair following radiation exposure,²⁵ and may be linked to enhanced tumor immunosurveillance.²⁶ In addition, Marinifilaceae has been reported to be reduced in patients with liver fibrosis.²⁷ Consistent with these findings, our results showed a relatively higher abundance of Marinifilaceae in the high-dose CMPLX group. Enterotoxigenic Bacteroides fragilis (ETBF), a member of the Bacteroidaceae family, has been associated with gastrointestinal inflammation and colorectal carcinogenesis, and its relative abundance appeared higher in the MOD group but lower in the CMPLX-treated groups.

At the genus level, Escherichia spp. were more abundant in the MOD group, whereas Odoribacter spp. were relatively enriched in the CMPLXH group. Increased abundance of Escherichia has been reported in immunocompromised cancer patients, which may explain the higher level of Enterobacteriaceae_A observed in the MOD group.²⁸ Conversely, reduced Odoribacter abundance has been associated with non-alcoholic fatty liver disease.²⁹ In addition, Kineothrix spp. were mainly observed in the CMPLXH group and have been reported to be associated with improved lipid metabolism and intestinal microbial homeostasis.³⁰ The abundance of Alloprevotella spp. was higher in the MOD group, which has been suggested as a potential marker of cancer metastasis.³¹

Functional prediction based on KEGG pathway analysis suggested that microbiota in the MOD group may be associated with pathways related to bacterial invasion of epithelial cells and altered metabolic processes, whereas the CMPLXH group showed enrichment in pathways related to beta-lactam resistance. Importantly, it should be noted that the associations observed between gut microbiota alterations and tumor-related changes in this study are correlative in nature and do not imply a causal relationship. Further studies are required to elucidate the underlying mechanisms.

The intricate link between metabolism and disease development is crucial for maintaining normal body function. We investigated the HCC-suppressive mechanism of CMPLX in balb/c mice from an in vivo metabolic perspective. Non-targeted serum metabolomics analysis showed that 84 variables underwent significant changes. Analyses of signaling-pathway enrichment revealed specific metabolic characteristics between the MOD group and the high- and low-dose groups of CMPLX. Carbon metabolism, signaling pathways, the TCA cycle, as well as the metabolism of cysteine, methionine, and riboflavin, are important pathways in cancer. This unique metabolic reprogramming may be an important internal mechanism mediating CMPLX’s anti-tumor activity. The metabolites of flavonoid components considered to have anti-HCC effects, such as vincamine and scutellarin, showed higher abundance in the high-dose group. Studies have indicated that vincamine may inhibit activin receptor-like kinase (ALK) and erb-B2 receptor tyrosine kinase 2 (ERBB2) proteins, which are implicated in the progression of lung cancer, as evidenced by research on targeted therapies for ERBB2 mutations.³² Moreover, vincamine can quench hydroxyl radicals and deplete iron ions in cancer cells,³³ thus exerting anticancer potential. Scutellarin can inhibit the proliferation of gastric cancer cells and induce their apoptosis through the wingless/integrated (Wnt)/β-catenin signaling pathway.³⁴ It can also inhibit the migration and apoptosis of human colon cancer HCT-116 and human colon adenocarcinoma RKO cell lines.³⁵ Luteolin has anti-inflammatory, anti-allergic, and anticancer effects and acts as an antioxidant or pro-oxidant. Its anti-inflammatory activity may be related to its anticancer properties.³⁶ In our study, levels of the derivatives of luteolin (luteolin 7-sulfate, luteolin 7,3′-diglucuronide, luteolin 4′-sulfate) and flavonoid-containing components such as kaempferol derivatives (kaempferol 3-glucuronide, kaempferol 3-glucuronide-7-sulfate) were significantly higher than those of the other two. These groups were considered to be metabolically derived from the two factors, and they are suggested to form the basis of immune stress conditions, which in turn facilitate the material basis of hepatocellular carcinoma.

Microbiota and metabolic correlation analyses revealed key intestinal flora associated with the antitumor activity of CMPLX. Heatmaps showed that luteolin 4′-sulfate, luteolin 7-O-diglucuronide, kaempferol 3-glucuronide, and kaempferol 3-glucuronide-7-sulfate were significantly negatively correlated with Escherichia species. Enterococcus_H also showed a significant negative correlation with the metabolites stated above. A computational method known as PARADIGM, developed by Nguyen et al, has been employed to scrutinize the correlation between drug exposure and fluctuations within the human intestinal microbiome.³⁷ This method has unveiled a link between the prevalence of Enterococcus species and the likelihood of death based on patterns of drug exposure. Therefore, the intervention of CMPLX could lead to the suppression of the abundance of this bacterial group and a reduction in the risk of cancer.

Our finding that two types of CMPLX, consisting of key anti-HCC components from four varieties of Epimedium species, had significant anti-tumor efficacy was supported by in vitro and in vivo experiments. Our study provides a foundation for establishing a research method for the development of innovative drug resources.

Conclusion

Our study shows that the combined use of kaempferol and icariin, both in vitro and in vivo, exerted significant anticancer efficacy against hepatocellular carcinoma.

Metabolomics and 16S rRNA sequencing analyses revealed the effects of CMPLX on central carbon metabolism in cancer, amino acids, and various metabolic pathways. CMPLX intake also affected the abundance of Enterobacteriaceae, Bacteroidaceae, Lachnospiraceae, and Marinifilaceae in the gut microbiota, which was ultimately able to suppress tumor development.

Thus, our results suggest that when combined with kaempferol and icariin, the results obtained in the in vitro experiments were consistently demonstrated in the in vivo experiments by altering the microbiota and therefore is a key factor in inhibiting the development of hepatocellular carcinoma.