Eight heat-processed HR products (labeled JZ, JPT, XFT, TQT, JJ, SFT, MDYY, and JD) were purchased from different manufacturers in China. Unprocessed JD HR was also obtained and designated SJD HR. Each HR product was prepared in triplicate, air-dried, flash-frozen in liquid nitrogen, and stored at −80 °C until analysis.

Targeted metabolomics analysis focusing on secondary metabolites was conducted following the standards outlined in the ConPhyMP statement (Supplementary Tables S1, S2) (Heinrich et al., 2022; Wang et al., 2022; Sun et al., 2023). Metabolites were extracted with 1,200 μL pre-cooled (−20 °C) 70% aqueous methanol containing 2-chlorophenylalanine (1 ppm; purity 98%) as an internal standard. After vortexing and centrifugation, the supernatant was filtered through a 0.22 μm microporous membrane (ANPEL, Shanghai, China) and transferred to an injection vial. Metabolomics analysis was performed using an ultra-high-performance liquid chromatography system (Nexera X2 UPLC, Shimadzu, Tokyo, Japan) coupled with an electrospray ionization-triple quadrupole-linear ion trap mass spectrometer (4500 QTRAP mass spectrometer, Applied Biosystems, Waltham, MA, USA). Chromatographic separation was achieved using an Agilent SB-C18 column (1.8 µm, 2.1 mm × 100 mm) with mobile solution A (0.1% formic acid in water) and solution B (0.1% formic acid in acetonitrile). Gradient elution was performed at a flow rate of 0.35 mL/min with an injection volume of 2 μL. The gradient program was: 95% A/5% B at 0 min, linearly ramped to 5% A/95% B at 9.0 min and held until 10.0 min, then returned to 95% A/5% B at 11.1 min and maintained until 14.0 min. The elution was directly introduced into the mass spectrometer with key operational parameters as follows: turbo spray ion source, source temperature 550 °C, ion spray voltage 5500 V (positive) and 4500 V (negative), ion source gas I 50 psi, gas II 60 psi, curtain gas 25 psi (Wang et al., 2022; Sun et al., 2023). The declustering potential and collision energy were individually optimized for each transition. Metabolite identification and quantification were conducted by searching the Metware Database (MWDB, MetWare Biological Co., Ltd., Wuhan, China). Chromatographic peaks were integrated and corrected using MultiQuant software (AB SCIEX, Framingham, MA, USA). SIMCA-P 18.0 software was used for principal component analysis (PCA), while significantly changed metabolites were determined using the orthogonal partial least squares discriminant analysis (OPLS-DA), followed by an independent t-test (SPSS 22.0 software) (Wan et al., 2022). VIP (Variable Importance in Projection) values from the OPLS-DA model reflected the contribution of each variable, with VIP values >1 indicating significance. The natural variation of the content of each secondary metabolite in HR was presented in a Z-score plot generated using R software (version 4.0.2) with three SD as the cutoff. Metabolite-metabolite correlation analysis was also performed using R software for Pearson’s product-moment correlation (Pearson’s r), p-values were calculated using the cor.test function.

The Simplified Molecular Input Line Entry System (SMILES) structure of formononetin was retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/) to predict its potential molecular targets based on its structural properties and bioactivity profiles. Target prediction was performed in SwissTargetPrediction (http://www.swisstargetprediction.ch/), TargetNet (http://targetnet.scbdd.com/), and PharmMapper (http://lilab-ecust.cn/pharmmapper/submitfile.html), which integrate multiple computational approaches to achieve comprehensive coverage of plausible targets. Duplicate entries were removed to retain unique targets. Targets related to CRC were compiled by querying DisGeNET (https://www.disgenet.org/), GeneCards (https://www.genecards.org/), and OMIM (https://www.omim.org/) databases using the keyword “Colorectal cancer” to select genes with documented relevance to CRC pathogenesis. Middle targets between formononetin and CRC were identified using Venny 2.1.0 (https://bioinfogp.cnb.csic.es/tools/venny/) to generate a Venn diagram including potential therapeutic targets of formononetin in CRC. The overlapped targets were then selected and subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses via the Bioinformatics online platform (https://www.bioinformatics.com.cn/) to elucidate the biological processes and signaling pathways involved.

Human SW620 colon cancer cells (ATCC, Manassas, VA) were cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum at 37 °C in a 5% CO₂ incubator. Cells were seeded in 96-well plates (2,000 cells/well) and treated with formononetin (catalog #F408902, Aladdin, Shanghai, China) at final concentrations of 0, 5, 10, 20, 40, and 80 μM for 48 h. Cell viability was assessed using CCK-8 reagent (catalog #DB884, Dojindo, Kumamoto, Japan). After 1 hour of incubation with CCK-8, absorbance at 450 nm was measured using a microplate spectrophotometer.

Proteomic profiling was performed using a label-free DIA approach to identify protein changes in formononetin-treated SW620 cells as previously reported (Wang et al., 2020; Wan et al., 2022). Proteins were extracted, reduced, and digested using the Filter-Aided Sample Preparation (FASP) method and analyzed on an EASY-nLC 1,200 system (Thermo Fisher Scientific, San Jose, CA, USA) coupled to a Q Exactive HF mass spectrometer (Thermo Fisher Scientific, San Jose, CA, USA). Chromatographic separation was performed on a reversed-phase C18 column with a binary gradient mobile solutions comprising 0.1% formic acid in water (A) and 0.1% formic acid in 80% acetonitrile (B). The gradient program was: 6%–10% B (0–5 min), 10%–30% B (5–47 min), 30%–45% B (47–55 min), and 95% B (55–60 min). MS/MS analysis in DIA mode was set at a resolution of 30,000, 32 isolation windows (first mass: 200 m/z), automatic maximum injection time, an AGC target of 3e6, and normalized collision energy of 28% (Wang et al., 2020). DIA raw data was processed using Spectronaut 18.0 (Biognosys) for protein identification. Differentially expressed proteins were defined by a fold change >1.2 and p value <0.05. GO and KEGG enrichment analyses of these proteins were conducted using OmicsBean (http://www.omicsbean.cn). The raw MS data have been deposited to the ProteomeXchange Consortium (https://proteomecentral.proteomexchange.org) via the iProX partner repository with the dataset identifier PXD059226 (Chen et al., 2022).

Phosphoproteomic analysis was performed on formononetin-treated SW620 cells using DIA-based methods. Protein extraction, digestion, and phosphopeptide enrichment were carried out using the High-Select™ Fe-NTA Kit (Thermo Fisher Scientific, San Jose, CA, USA) following the manufacturer’s instructions (Song et al., 2022). Briefly, lyophilized peptides were reconstituted in 200 μL of immobilized metal affinity chromatography (IMAC) binding/wash buffer. IMAC spin columns were activated by centrifugation (1,000 × g, 30 s) and equilibrated with two washes of the same buffer. After loading the samples, the resin was gently resuspended and incubated for 30 min with manual mixing at 10-min intervals. Following incubation, the columns were washed three times with binding/wash buffer and once with HPLC-grade water. Phosphopeptides were eluted twice with 100 μL elution buffer each time (1,000 × g, 30 s each), and the two elution were combined and dried under vacuum. The eluted phosphopeptides were analyzed on an EASY-nLC 1,200 system coupled to a Q Exactive HF mass spectrometer under the same parameters as in the proteomic analysis. Phosphorylation site abundances were normalized by dividing the abundance of each site by its corresponding protein abundance (Fan et al., 2020). Phosphorylation sites exhibiting significant alterations (fold change >1.2 and p < 0.05) were subjected to functional annotation through GO and KEGG pathway enrichment analyses.

Molecular docking simulations were conducted using MOE software (version 2019.0102). The molecular structures of formononetin and four reference drugs (5-Fluorouracil, Capecitabine, Oxaliplatin, and Irinotecan) were retrieved from the PubChem database. The structures of target proteins CA9, MME, and YAP1 were obtained from the Protein Data Bank (https://www.rcsb.org/). Protein structures were processed in MOE, and semi-flexible docking simulations were performed to evaluate interactions between formononetin and the selected targets.

A targeted UPLC-ESI-QTRAP-MS/MS metabolomic analysis focusing on secondary metabolites was conducted on eight heat-processed HR products from different manufacturers. We identified 1,292 metabolites (Supplementary Table S3), classified into eight major categories: flavonoids, phenolic acids, terpenoids, alkaloids, lignans and coumarins, and others. Flavonoids constituted the largest category (29.41%), followed by phenolic acids (18.89%), others (12.54%), alkaloids (12.31%), terpenoids (12.31%), and lignans and coumarins (10.91%) (Figure 1A). There were 26 metabolites that accounted for more than 1% each of the total relative content of the HR metabolome, such as 3-hydroxy-1-methylpyrrolidin-2-one, hypaphorine, formononetin, and ferulic acid. Notably, formononetin and its five derivatives including formononetin-7-O-glucoside and formononetin-7-O-(6″-malonyl)glucoside constituted approximately 4.22% of the total HR metabolome, while the most abundant single metabolite, 3-hydroxy-1-methylpyrrolidin-2-one, was 3.06% of the HR metabolome.

We further searched the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) to screen for metabolites with potential drug activities (Li et al., 2023). The selection criteria were oral bioavailability (OB) ≥30%, drug-likeness (DL) ≥0.18, molecular weight (MW) ≤500, logP ≤5, hydrogen bond donors (nOHNH) ≤5, and hydrogen bond acceptors (nOH) ≤10. The criteria were consistent with standard TCMSP thresholds. OB ≥ 30% ensures sufficient systemic exposure after oral administration, DL cutoff of 0.18 corresponds to the mean value of known drug molecules in DrugBank. Using these criteria, we identified 11 metabolites with potential drug activity (Table 1), and six of them—medicarpin, calycosin, formononetin, naringenin, quercetin, and vestitol—were detected in our metabolomic profiling of HR.

To examine metabolic alterations induced by heat-processing (called paozhi in Chinese) of HR, we performed UPLC-ESI-QTRAP-MS/MS analysis on SJD, the unprocessed form of JD HR. The PCA plot demonstrated clear separation between the three replicates of SJD HR and the other heat-processed HR samples including heat processed JD HR (Figure 1B). Processed HR samples clustered more tightly than SJD HR, indicating increased metabolic consistency post-processing. Z-score analysis identified 11 metabolites in SJD HR as outliers relative to processed HR products, including rhamnetin-3-O-rutinoside-5-O-rhamnoside, chrysoeriol-5,7-di-O-glucoside, physcion-8-O-(6-acetyl)-glucoside, and trans-5-O-(p-coumaroyl)shikimate (Figure 1C). Six metabolites had increased levels post-processing, and three of them were flavonoids, including Rhamnetin-3-O-rutinoside-5-O-rhamnoside and chrysoeriol-5,7-di-O-glucoside whose absolute Z-scores were bigger than 5.

The metabolic differences between the unprocessed and processed JD HR were also investigated using supervised OPLS-DA, and 108 differentially expressed metabolites (VIP >1, p ≤ 0.05) were identified to contribute to the separation of the two groups in OPLS-DA plot (Figure 2). Among these metabolites, 29 were significantly upregulated post-processing, approximately 34% of which were flavonoids, including acetyl wistin, bonannione A, and hesperidin (Supplementary Table S4). Although the total formononetin content in JD HR did not change significantly after heat processing, two formononetin derivatives, 8-methoxy acetyl ononin and formononetin acetyl glucoside, were significantly upregulated. These results indicated that heat processing significantly changed the HR metabolome and particularly increased the contents of flavonoid derivatives.

To reveal the metabolite regulatory network in HR, we performed correlation analysis among the 1,292 detected HR metabolites. The resulted heatmap included 833,986 correlations ranging from −0.9939 to 1 (Figure 2C). When setting the threshold of r² ≥ 0.49 (r ≥ 0.7 or r ≤ −0.7) and p ≤ 0.05, 88,955 significant correlations remained including 67,457 positive and 21,498 negative correlations. Flavonoids dominated these significant metabolite-metabolite correlations, with 41,174 significant correlations among 380 flavonoids (32,739 positive and 8,435 negative). Next to flavonoids, there were 30,908 significant correlations among 244 phenolic acids, 20,583 significant correlations among 141 lignans and coumarins, 18,625 significant correlations among 159 terpenoids, and 18,520 significant correlations among 159 alkaloids. When screening highly significant correlations with formononetin by r² ≥ 0.81, 26 correlations were identified (Figure 2D). Four of these 26 correlations were negative, including three alkaloids (hypaphorine, valerine, salicylamide) and one terpenoid (6′-O-Glucosylaucubin). Formononetin had high positive correlation with ten flavonoids, including prunetin, afrormosin, biochanin A, and genkwanin. The other 12 high correlations with formononetin were five lignans and coumarins, four phenolic acids, and three quinones.

Since it has been established via Zinc database analysis (https://zinc.docking.org/) that formononetin is not a pan-assay interference compound (Bolz et al., 2021), we performed network analysis to search formononetin’s putative molecular targets against CRC. We retrieved 220 predicted targets from SwissTargetPrediction, TargetNet, and PharmMapper databases. A total of 14,050 CRC-related genes were obtained from GeneCards, OMIM, and DisGeNET. The Venn diagram revealed 194 overlapping genes between predicted formononetin targets and CRC-related genes (Figures 3A,B). GO enrichment analysis of the 194 genes identified 517 enriched terms: 358 biological processes (BP), 43 cellular components (CC), and 116 molecular functions (MF). The top 10 enriched terms in each category were shown in Supplementary Figure S1, highlighting key processes such as peptidyl-tyrosine phosphorylation, cytosol localization, receptor complex formation, cytoplasmic signaling, steroid binding, zinc ion binding, and protein kinase activity. KEGG pathway enrichment analysis of the 194 genes identified 114 significantly enriched pathways. The top 30 pathways (Figure 3C) included several cancer-related signaling pathways such as pathways in cancer, MAPK signaling, Ras signaling, and PI3K-Akt signaling. These results suggested that formononetin may exert anti-CRC effects by modulating multiple oncogenic signaling pathways.

We evaluated the inhibitory effect of formononetin on SW620 colon cancer cells by measuring cell viability using the CCK-8 assay. Compared to the untreated control cells, formononetin treatment significantly reduced SW620 cell viability at tested concentrations (5, 10, 20, 40, and 80 μM), although the effect was not strictly dose-dependent (Figure 4A). The concentration of 10 μM was selected for further proteomics study since it gave pronounced inhibitory effect without inflicting overwhelming cytotoxicity. DIA proteomics identified 7,381 proteins, of which 291 were significantly altered upon formononetin treatment, and the fold changes ranged from 0.0019 to 8.2382 (Figure 4B; Supplementary Table S5). Specifically, 105 proteins were significantly upregulated, while 186 proteins were significantly downregulated including CCDC18, CA9, IGF1R, and MME.

GO enrichment analysis of these 291 differentially expressed proteins identified 3,457 enriched BP terms, 372 CC terms, and 404 MF terms. The top ten enriched GO terms in each category were displayed in Supplementary Figure S2, including single-organism cellular process, multicellular organismal process, cellular anatomical entity, cytoplasmic part, protein binding, and carboxylic acid binding. KEGG pathway enrichment analysis revealed 17 significantly enriched pathways, including metabolic pathways, Hippo signaling, HIF-1 signaling, nitrogen metabolism, one-carbon pool by folate, and cholesterol metabolism (Figure 4C). KEGG pathway enrichment analyses of upregulated and downregulated proteins were also conducted and shown in Supplementary Tables S6, S7. These findings indicated that formononetin caused broad proteomic changes involving multiple metabolic and oncogenic pathways in SW620 cells.

Further investigation into formononetin’s molecular mechanisms of cancer cell inhibition was conducted using DIA-based phosphoproteomic analysis. We quantified 19,038 phosphorylation sites, among which 2,587 sites in 1,535 phosphorylated proteins were significantly altered (Figure 5A; Supplementary Table S8). These phosphorylation sites comprised 79.1% phosphoserine, 16.4% phosphothreonine, and 4.5% phosphotyrosine modifications. Among them, 1,244 phosphorylation sites were significantly upregulated (involving proteins such as VAPA, CA9, YY1, FOXK2, SQSTM1), while 1,343 sites were downregulated (in proteins such as AKT2, YWHAB, MME, YAP1). GO enrichment analysis of the phosphoproteins identified 7,832 significantly enriched terms (p < 0.05): 6,104 BP, 898 MF, and 830 CC. The most enriched terms were cellular component organization or biogenesis, protein binding, poly(A) RNA binding, intracellular organelle part, and nuclear compartmentalization (Supplementary Figure S3).

KEGG pathway enrichment analysis highlighted 41 significantly enriched pathways (p < 0.05), such as cell cycle, autophagy, Hippo signaling, mTOR signaling, AMPK signaling, glycolysis/gluconeogenesis, lysine degradation, and insulin signaling (Figure 5B). While there was no commonly enriched pathways across network analysis, proteomics, and phosphoproteomics, 19 pathways were enriched by at least two analyses, including Hippo signaling and lysine degradation by both proteomics and phosphoproteomics, and proteoglycans in cancer and glycolysis/gluconeogenesis by both network analysis and phosphoproteomics (Figure 5C). These findings highlighted formononetin’s multifaceted regulatory effects on key oncogenic and metabolic pathways.

To further validate key molecular targets of formononetin, we compared proteins identified by network analysis, proteomics, and phosphoproteomics (Figure 5D). Two targets, CA9 and MME, were consistently identified across all three analyses, suggesting their central roles in formononetin-mediated anti-cancer effects. Additionally, 50 proteins were common between proteomic and phosphoproteomic datasets (Supplementary Table S9), while 11 targets, including FYN and DNMT1, were shared between network analysis and phosphoproteomics. We confirmed the binding affinity of formononetin to CA9, MME, and YAP1 (a key regulator of Hippo signaling) by performing molecular docking. Control drugs 5-Fluorouracil, Capecitabine, Oxaliplatin, and Irinotecan were included for comparison. The binding energy scores for formononetin-CA9, formononetin-MME, and formononetin-YAP1 in molecular docking analysis were −4.8162, −6.1931, and −4.3183 kcal/mol, respectively. These scores were higher than those for Capecitabine and Irinotecan, but lower than those for 5-Fluorouracil and Oxaliplatin (Supplementary Table S10). The complete docking analysis results demonstrated that formononetin had strong binding affinity to CA9, MME, and YAP1, indicating they were potential formononetin targets (Supplementary Figure S4).