To ensure highly reliable target spectrum data for silymarin, 10 target databases were accessed to gather both known and predicted target profiles, as shown in Table I. The UniProt database (http://www.uniprot.org) (16) was employed to standardize the names of the targets, and only genes corresponding to ‘Homo sapiens’ were retained for further analysis.

‘MM’ was used as the search term to search for pathogenic genes associated with MM from five sources: Open Target Platform (17), DisGeNET (18), Comparative Toxicogenomics Database (19), GeneCards (20) and text mining. Gene names were standardized using DAVID v6.8 (https://david.abcc.ncifcrf.gov/) (21). To enhance data reliability, only genes that were confirmed or listed in ≥3 database sources were retained for further analysis.

To identify potential candidate genes for silymarin in treating MM, the targets of silymarin were cross-referenced with genes associated with MM. The genes found at the intersection were regarded as prospective candidates for action of silymarin on MM. A network model was then created using Cytoscape software v3.7.1 (22) to depict the interactions between silymarin, its targets and MM.

To further identify the key targets of silymarin in treating MM, a protein-protein interaction (PPI) network was constructed through the online STRING database version 11.5 (https://string-db.org/) (23). The PPI networks were visualized using Cytoscape v3.7.1 (22). The central nodes in the network, or hub nodes, were assessed based on key topological parameters, specifically the degree (24). The NetworkAnalyzer plugin (25) was employed to calculate the degree (24).

GO and KEGG enrichment analyses for the candidate genes of silymarin in MM were performed using clusterProfiler version 4.0.3, an R Bioconductor package (26). During enrichment analysis, the thresholds for the P-value and q-value were set at 0.05. This criterion was used to determine the significantly enriched GO terms and KEGG pathways associated with the candidate genes.

The molecular docking method was employed to virtually assess the binding affinity between the compound and target proteins, providing a basis for subsequent experimental validation. AutoDock Vina 1.1.2, an open-source software for molecular docking and virtual screening, offers markedly improved accuracy in predicting binding modes compared with AutoDock 4 (27). A spatial data file format of silymarin was downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) (28) and converted into Mol2 format using Open Babel 2.4.1 (29). The crystal structure of JAK2 [Protein Data Bank (PDB) identification (ID), 3UGC] and STAT3 (PDB ID, 6NJS) were downloaded from the PDB database (https://www.rcsb.org/) (30) NVP-BBT594 (PubChem CID, 59596344) and SD36 (PubChem CID, 139600321), co-crystalline ligands of JAK2 and STAT3 respectively, were selected as positive controls. Small molecules were energy minimized using PyRx-0.8 software (Dr. Sargis Dallakyan; The Scripps Research Institute). Both ligands and receptors were prepared following the AutoDock Vina 1.1.2 tutorial. During rigid docking, grid properties were determined based on the spatial position of the co-crystallized ligands within the proteins. Each structure was prepared by removing water molecules, adding non-polar hydrogen, calculating Gasteiger charges and saving in PDBQT format. Typically, a lower Vina score indicates a higher affinity between the ligand and receptor, with a general threshold for binding energy being ≤-5 kcal/mol (31,32). The interaction modes between the compound and target proteins were visualized using PyMOL software version 3.1 (Schrödinger, Inc.).

The RPMI 8226 and H929 cell lines were obtained from American Type Culture Collection. These cell lines were cultured in RPMI 1640 medium (Gibco; Thermo Fisher Scientific, Inc.), supplemented with 10% fetal bovine serum (Biological Industries), and maintained without antibiotics at 37°C in an incubator with 5% CO₂. The culture medium was refreshed every 2 days.

A total of 1×10⁵ cells/well were seeded in 96-well plates. The cells were treated with silymarin (cat. no. HY-N7073; MedChemExpress) dissolved in dimethyl sulfoxide (DMSO) at gradient concentrations from 0 to 100 µM (in 10 µM increments), or with vehicle control (0.1% DMSO, equivalent to the highest solvent concentration in treatment groups) for 24 h at 37°C prior to the addition of Cell Counting Kit-8 (CCK-8) reagent. A total of 100 ml serum-free medium supplemented with 10 ml CCK-8 reagent (MilliporeSigma) was then added to each well and incubated for 2 h at 37°C in the dark. Cell viability was assessed by measuring the absorbance at 450 nm using a spectrophotometric plate reader (BioTek; Agilent Technologies, Inc.).

In a preliminary experiment, cells began to die within 24 h of silymarin administration when examined under a microscope. Therefore, 24 h was selected as the treatment duration. After the treatment of RPMI 8226 or H929 cells, which were seeded into 6-well plates at a density of 1×10⁶ cells/well, with 50 µM silymarin for 24 h at 37°C, the cells were collected and centrifuged at 1,000 × g for 5 min at 4°C. The cell suspensions were then washed twice with chilled phosphate-buffered saline (PBS) and incubated with 10 µl Annexin-V FITC (100 mg/ml; BD Biosciences) for 30 min in the dark. Subsequently, 5 µl PI (100 mg/ml; BD Biosciences) was added, and the cells were incubated for an additional 30 min in the dark on ice. Apoptotic rates were analyzed using flow cytometry (FACSCalibur; BD Biosciences). Data acquisition was performed with BD CellQuest™ Pro software (version 5.2; BD Biosciences), and quantitative analysis was conducted using FlowJo software (version 10.8.1; FlowJo LLC).

Total RNA from RPMI 8226 and H929 cells was isolated using the TransZol Up Plus RNA Kit (cat. no. ER501-01; TransGen Biotech Co., Ltd.) and subsequently reverse transcribed using the GoScript™ Reverse Transcription System (cat. no. A5001; Promega Corporation). The resulting cDNA was used as a template for qPCR, which was performed using the MonAmp™ ChemoHS qPCR Mix kit (cat. no. MQ00401S; Monad Biotech Co., Ltd.). The primer sequences used are detailed in Table II. The following thermocycling protocol was used: Stage 1, activation at 50°C for 2 min; stage 2, pre-soak at 95°C for 10 min; stage 3, 40 cycles of denaturation at 95°C for 15 sec, annealing at 60°C for 1 min; stage 4, melting curve, 95°C for 15 sec, 60°C for 15 sec, 95°C for 15 sec. All procedures were performed according to the manufacturer's guidelines.

Cells were fixed with 4% paraformaldehyde for 10 min at room temperature, followed by permeabilization with 0.1% Triton X-100 for 15 min. Subsequently, 5% bovine serum albumin (BSA; cat. no. A9418; MilliporeSigma) was applied to block nonspecific binding for 15 min. Primary antibodies targeting JAK2 (1:500; cat. no. PA5-11267; Thermo Fisher Scientific, Inc.), phosphorylated-JAK2 (1:500; cat. no. PA5-99341; Thermo Fisher Scientific, Inc.), STAT3 (1:100; cat. no. MA1-13042; Thermo Fisher Scientific, Inc.) and phosphorylated-STAT3 (1:100; cat. no. PA5-17876; Thermo Fisher Scientific, Inc.) were added to the cells and incubated overnight at 4°C. After washing three times with PBS, the cells were incubated with Alexa Fluor™ 488-conjugated goat anti-rabbit IgG (1:500; cat. no. SA00013-2; Proteintech Group, Inc.) for JAK2, p-JAK2, and p-STAT3 detection, and Alexa Fluor™ 594-conjugated goat anti-mouse IgG (1:500; cat. no. SA00013-3; Proteintech Group, Inc.) for STAT3 at 37°C for 1 h. Nuclei were stained with DAPI (1 µg/ml) for 15 min in the dark prior to imaging. Imaging was performed using a ZEISS laser scanning confocal microscope (ZENA GmbH). Quantitative analysis of fluorescence intensity was conducted using ImageJ software (version 1.53q; National Institutes of Health) by defining regions of interest (ROIs) around nuclei (DAPI channel) and measuring mean pixel intensity in corresponding JAK2/STAT3 channels. Data were normalized to background fluorescence (unstained controls) and expressed as fold change relative to control groups.

All results are presented as the mean ± standard deviation for normally distributed data or median (interquartile range) for non-normally distributed data. Statistical analyses were performed using SPSS 22 (IBM Corp.). Initially, data were tested for normal distribution using the Shapiro-Wilk test and for homogeneity of variances using Levene's test. Based on these results, appropriate statistical tests were selected. One-way analysis of variance (ANOVA) was utilized to assess the significance of differences in normally distributed data, followed by Dunnett's multiple comparisons test. For data that were not normally distributed, the Kruskal-Wallis test was applied followed by Dunn's test. For comparisons between two independent groups, an independent Student's t-test was applied. P<0.05 was considered to indicate a statistically significant difference for all analyses.

From 10 target databases, a comprehensive set of 594 targets associated with silymarin, encompassing both known and predicted targets, was compiled. Additionally, by querying six databases and selecting only those genes appearing in at ≥3 sources, a total of 2,178 genes strongly associated with MM were identified. Intersection analysis of the silymarin targets and MM genes revealed 221 common genes. These genes were considered potential targets for silymarin in the treatment of MM and were selected for further analysis (Fig. 1A).

The 221 overlapping targets identified were subsequently uploaded into the STRING database to construct a PPI network. The resulting PPI network data were exported as a TSV file and subsequently analyzed using Cytoscape, revealing 209 connected targets and 1,359 interactions between them. Analysis of the network based on degree centrality highlighted several key targets including TP53, HSP90AA1, AKT1, BCL2, STAT3, GAPDH, EGFR, HSP90AB1, SRC, ESR1, IL6, CCND1, NFKB1, CASP3, IL1B, BRCA1, MAPK1, MAPK3, BCL2L1, MDM2, RELA, PPARG, SIRT1, GSK3B, STAT1, AR, MMP9, NR3C1, JAK2 and CASP9 (Fig. 1B).

The candidate genes of silymarin in MM were enriched in 159 pathways, including ‘Apoptosis’, ‘Prostate cancer’, ‘Chemical carcinogenesis-receptor activation’, ‘JAK-STAT signaling pathway’, ‘Endocrine resistance’, ‘Pancreatic cancer’, ‘Chronic myeloid leukemia’, ‘FoxO signaling pathway’, ‘MicroRNAs in cancer’, ‘Prolactin signaling pathway’ and ‘PI3K-Akt signaling pathway’ (Fig. 2A). Silymarin was also suggested to exert therapeutic effects on MM by engaging in 2,287 biological processes, such as ‘cellular response to peptide’, ‘response to oxidative stress’, ‘positive regulation of proteolysis’, ‘cellular response to chemical stress’ and ‘regulation of apoptotic signaling pathway’ (Fig. 2B). Additionally, silymarin in MM was noted to be involved in 191 molecular functions, including ‘ubiquitin-like protein ligase binding’, ‘ubiquitin protein ligase binding’, ‘DNA-binding transcription factor binding’, ‘RNA polymerase II-specific DNA-binding transcription factor binding’ and ‘transcription coregulator binding’ (Fig. 2C). Furthermore, it was demonstrated to interact with 89 cellular components, including ‘membrane raft’, ‘membrane microdomain’, ‘protein kinase complex’, ‘transferase complex’ and ‘serine/threonine protein kinase complex’ (Fig. 2D).

To evaluate the therapeutic effects of silymarin on MM, its impact on cell viability was investigated using the CCK-8 assay. The dose-response curves demonstrated that silymarin has a dose-dependent inhibitory effect on cell viability in both MM cell lines (Fig. 3A and B). Furthermore, flow cytometric analysis revealed that silymarin promoted early and late apoptosis in MM cells compared with that in the control groups (Fig. 3C-F).

Network pharmacology analysis identified JAK2 and STAT3 as central nodes within the PPI network (Fig. 1B), and KEGG enrichment analysis pinpointed the ‘JAK-STAT signaling pathway’ as a crucial mechanism of silymarin in MM (Fig. 2A). Previous research has reported that inhibiting the JAK2/STAT3 pathway can induce apoptosis in MM cells (30,31). These results suggested that the JAK2/STAT3 pathway may be notable in apoptosis induced by silymarin. To further evaluate the molecular mechanisms by which silymarin can induce MM apoptosis, the expression of the anti-apoptotic genes Bcl-2 and Bcl-xL, which are transcriptional targets of JAK2/STAT3 signaling, were assessed in silymarin-treated MM cells. The results demonstrated that silymarin reduced the expression levels of Bcl-2 and Bcl-xL in the tested MM cell lines at both 10 and 50 µM compared with those in the control groups (Fig. 4A-D). Moreover, the phosphorylation and nuclear localization of JAK2 and STAT3 were reduced with silymarin treatment compared with those in the control groups (Fig. 4E-K). Collectively, these findings indicated that silymarin may inhibit the JAK2/STAT3 signaling pathway.

At the molecular level, silymarin was demonstrated to potentially interact with JAK2 and STAT3. Based on empirical rules, the effective binding energy between small molecules and target protein should be <-5 kcal/mol. The present findings of molecular docking analysis indicated that silymarin has a particularly high affinity for JAK2, with a binding energy of −9.4 kcal/mol (Fig. 5A). The docking analysis of silymarin with JAK2 suggested that silymarin fits well into the active site of the JAK2 protein. It was shown to form hydrogen bonds with the amino acids Glu930, Asp994 and Phe995 and to engage in a pi-sigma interaction with Leu902. Additionally, pi-alkyl interactions were noted with Leu901, Leu983, Val911, Ala880 and Met929. In the silymarin-STAT3 protein docking model (Fig. 5B), silymarin was shown to be located within the binding site of the STAT3 protein, where it can form hydrogen bonds with Ile653 and pi-sulfur interactions with Met648, along with pi-alkyl interactions with Ile653.