Introduction

Lung cancer is one of the tumors with the highest incidence and mortality in the world, with an overall survival rate of approximately 17%, and the 5-year survival rate of distant metastatic tumors is only 5%.1 The proportion of non-small cell lung cancer (NSCLC) cases among all lung cancer cases accounts for 80–85%, and the main pathological types are squamous cell carcinoma and adenocarcinoma. Currently, the prevalent treatments of lung cancer mainly include surgery, radiotherapy, chemotherapy, targeted drugs and immunotherapy. However, certain limitations and adverse effects associated with them still exist. Notably, patients’ metastasis of NSCLC is one of the primary contributors of treatment failure and death in clinical practice.2 Therefore, how to elucidate the metastasis mechanism of NSCLC and identify the targeted therapeutic drugs remains a clinical challenge. As an inseparable part of complementary medicine in the world, traditional Chinese medicine (TCM) has been recognized and adopted in the treatment of lung cancer for more than 2000 years. An increasing number of studies have exhibited that TCM is widely recognized in experimental studies and clinical treatments, providing novel ideas and inspiring references for the treatment of NSCLC.3,4

Qianjin Weijin Decoction (QJWJD) is a classical TCM prescription originated from Prescriptions Worth a Thousand Gold for Emergencies (Beiji Qianjin Yaofang) by Sun Simiao, a well-renowned medical scientist in the Tang Dynasty. The prescription consists of Phragmites australis (Cav). Trin. ex Steud., Coix lacryma-jobi L., Benincasa hispida (Thunb). Cogn and Prunus persica (L). Batsch.5 For almost 1500 years, it has been one of the most commonly-used prescriptions in the treatment of lung diseases in TCM. Fundamental researches have confirmed that QJWJD exhibits certain efficacy in the treatment of lung cancer in vitro and in vivo.6,7 Clinical studies have shown that this prescription can significantly relieve clinical symptoms and enhance immunity of patients.5 However, for lung cancer metastasis, QJWJD has not shown favorable curative effects in the current studies. According to the fundamental theory of TCM, lung cancer metastasis is tightly associated with blood stasis and heat generation.8 Taking this into consideration, a modified QJWJD, Jiawei Weijin Decoction (JWWJD) is warranted by adding Curcumae Rhizoma to activate blood circulation, Prunella vulgaris L. to clear away heat and Ophiopogon japonicus (Thunb). Ker Gawl to nourish yin, thus enhancing the inhibitory effect on the lung cancer metastasis. In the JWWJD, Curcumae Rhizoma is the monarch medicinal and it is combined with Prunus persica (L). Batsch to eliminate tumors and inhibit metastasis by promoting blood circulation and removing blood stasis. Recent studies have demonstrated that numerous active ingredients in Curcumae Rhizoma have anticancer effect.9 As the ministerial drugs, Prunella vulgaris L. and Ophiopogon japonicus (Thunb). Ker Gawl. can strengthen the role of Phragmites australis (Cav). Trin. ex Steud. in clearing away heat and removing the lung toxins.9 Previous pharmacological studies have shown that triterpenoid saponins, the main ingredient of Prunella vulgaris L., and steroid saponins, the main ingredient of Ophiopogon japonicus (Thunb). Ker Gawl., demonstrate desired efficacy in treating NSCLC.10,11 However, due to the multi-component, multi-target and multi-pathway characteristics, the working ingredients and mechanisms of compound prescriptions of TCM are difficult to be identified. Thereby, elucidating the mechanism of action remains important and challenging in modern scientific researches.

In this study, we present the first comprehensive investigation and validation of JWWJD and its active components through an integrated approach combining network pharmacology, transcriptomics-based bioinformatics analysis, and rigorous in vitro and in vivo experiments. Using LC-MS/MS technology, we identified curcumol as the key active ingredient in JWWJD that effectively targets SPP1, thereby inhibiting the proliferation, migration, and angiogenesis of NSCLC. These findings not only provide novel mechanistic insights into the anti-metastatic action of JWWJD but also highlight curcumol as a promising candidate for the development of targeted therapies against NSCLC. The detailed procedure of this study is shown in Figure 1.

Figure 1 Flow chart of the study. *P < 0.05, **P <0.01, ***P <0.001 vs Control.

Materials and MethodsMaterials Grades and Certification

All chemicals and reagents used in this study were selected based on their quality grades and suitability for the intended experimental applications. Curcumol (Cat: IC0620; purity ≥ 98%), used for in vitro and in vivo experiments, was purchased from Solarbio (Beijing, China), and was confirmed to be of high purity suitable for pharmacological studies. Solvents and cell culture reagents were of analytical or cell culture grade, as appropriate for each experiment. Detailed information on other materials, including antibodies and assay kits, is provided in the respective method sections.

Human Specimens and Ethical Statements

Tumorous and adjacent lung tissues were collected from NSCLC patients underwent surgical resection and pathological confirmation in the Hubei Provincial Hospital of Integrated Chinese and Western Medicine. Tumor differentiation grade and clinical stage were classified according to the American Joint Committee on Cancer TNM Classification (7th edition). None of the patients had distal organ metastases or had received chemotherapy or radiotherapy before surgery. Personal and clinical information for the aforementioned patients with NSCLC was obtained from medical records. All tissue specimens were collected with informed consents and in strict compliance with the Declaration of Helsinki. Sixty paraffin-embedded specimens collected from January 2014 to December 2019 were used for Immunohistochemistry (IHC) staining, and a total of 20 fresh tumor tissue samples collected from June to December 2023 were subjected to Western blotting (WB) analysis. The study was approved by the Institutional Ethics Committee of Hubei Provincial Hospital of Integrated Traditional Chinese and Western Medicine (Approval no. 2023007).

As our laboratory does not have its own animal facility to conduct animal experiments, we carried out the experiments at the animal research center of Shouzheng Pharma (Wuhan) Biotechnology Co., Ltd. This facility has an independent Institutional Animal Care and Use Committee (IACUC), which carried out the ethical review and oversight of all animal procedures. Animal studies were conducted in accordance with the ethical regulations approved by the Animal Welfare Ethics Review Committee of Shouzheng Pharma (Wuhan) Biotechnology Co., Ltd. (Approval No. 2024071601, Date of ethical approval: 14 July 2024, Wuhan, Hubei Province, China). All experimental procedures complied with the Laboratory animal — Guideline for ethical review of animal welfare (GB/T 35892–2018), as well as internationally recognized principles for the care and use of laboratory animals, including adherence to the to the Basel Declaration.12 Every effort was made to minimize animal suffering and distress. Humane endpoints were strictly applied throughout the study, and appropriate anesthesia (sodium pentobarbital, intraperitoneal injection) was administered during surgical procedures. Euthanasia was performed via cervical dislocation under deep anesthesia, in accordance with the recommendations of the ethics committee.

Preparation and Extraction of Drug-Containing Serum

The medicinal names have been checked with the plant list (http://www.theplantlist.org/) on October 28, 2024. And all the herbal granules forming JWWJD (Batch NO.20240616) were purchased from Hubei Tianji Pharmaceutical Co., Ltd (Wuhan, Hubei Province, China). The dosage of each medicinal material is shown in Table 1. The granules were fully dissolved in 150 mL of boiled distilled water under continuous stirring. The mixture was then heated at 80–90 °C for approximately 30 minutes to ensure complete extraction of active components. Subsequently, the decoction was filtered through a four-layer gauze and further concentrated using a rotary evaporator under reduced pressure at 60 °C until the final concentration reached 2 g/mL (calculated as crude herb weight per mL). The resulting concentrated solution was stored at −80 °C for subsequent use in animal experiments and serum preparation.

Table 1 Herbal Formula of JWWJD

Sprague Dawley rats (8 weeks old, weighing 300±20 g, n = 20) were purchased from Shulaibao (Wuhan) Biotechnology (Wuhan, Hubei Province, China). Animals were raised in a standard Specific Pathogen Free (SPF) animal feeding center (25°C, 12-h light/dark cycle, 55% humidity, n=3 rats of the same sex in each cage). Following one week of acclimatization, the rats were randomly divided into two groups (n = 10 per group): the JWWJD-treated group and the Control group. The JWWJD group were gavaged with 11.0 g raw drug/kg twice a day for 3 days, while the control group were gavaged with an equal volume of distilled water. On the 4th day, 1 h after the last gavage, all rats were anesthetized by intraperitoneal injection of 2% pentobarbital sodium (30 mg/kg). The rats were sacrificed by cervical dislocation after abdominal aorta blood was collected. The blood was incubated at room temperature for 1 h and then centrifuged at 1500 x g at 4°C for 15 min. The supernatant was extracted, inactivated at 56°C, filtered by 0.22-µM filters and stored at −80°C.

Data Source and Processing

The selection of gene expression profiles and data processing were based on publicly available datasets from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). Human non-small cell lung cancer (NSCLC) samples were retrieved, and two datasets, GSE101929 and GSE31210, were selected for analysis. The raw microarray data were processed using R software (version 4.3.1), including background correction, normalization, and log2 transformation, which were performed using the “limma” package (version 3.56.2), a comprehensive tool for analyzing gene expression microarray data.13

Missing values in the datasets were imputed using the K-Nearest Neighbors (KNN) algorithm implemented in the “impute” package (version 1.74.0), with k=10 neighbors selected through cross-validation to minimize imputation error.14 To correct for potential batch effects arising from differences in experimental platforms and sample handling across datasets, Surrogate Variable Analysis (SVA) was applied using the “sva” package (version 3.48.0). The batch-corrected datasets were then merged into a single training dataset for downstream analysis.15

For independent validation, two additional GEO datasets, GSE29013 and GSE42127, were selected and preprocessed following the same procedure described above, and subsequently combined to form the test set. These datasets were used to validate the stability and reproducibility of the differentially expressed genes and prognostic signatures identified in the training set.

Differential Expression Screening and Weighted Gene Co-Expression Network Analysis (WGCNA)

The WGCNA package in R software was utilized for weighted gene co-expression network analysis between NSCLC and normal tissues. After filtering low-variance genes and performing log2 transformation and quantile normalization, outlier samples were removed based on hierarchical clustering. A scale-free network was constructed by selecting an appropriate soft-thresholding power β. The power β = 7 was chosen as the minimum value that achieved a scale-free topology fit index (R2) > 0.85.

A topological overlap matrix (TOM) was derived from the adjacency matrix, and modules were identified using average linkage hierarchical clustering with a dynamic tree cut method. Modules with eigengene correlation > 0.75 were merged to reduce redundancy.

Module significance (MS), defined as the -log10(P-value) of the correlation between the module eigengene and the NSCLC phenotype, was used to assess the relevance of each module. Pearson correlation analysis was employed to select the module showing the highest correlation with NSCLC, which was designated as the key module for downstream analysis.

Prediction of JWWJD Active Ingredients and Targets for Treatment of NSCLC Plus Construction of a Protein-Protein Interaction (PPI) Network

The active ingredients and corresponding targets of JWWJD were retrieved from the TCMSP (https://tcmspw.com/tcmsp.php) and BATMAN (http://bionet.ncpsb.org/batman-tcm/) databases. Target names were standardized to official gene symbols using the UniProt database (https://www.uniprot.org/). To identify core therapeutic targets of JWWJD in NSCLC, a Venn diagram analysis was performed to obtain the intersection among three datasets: Drug targets predicted for JWWJD, Differentially expressed genes (DEGs) identified in NSCLC samples from the training set (|log2FC| > 0.485, P < 0.05), Genes belonging to the key module (turquoise module) derived from WGCNA, which showed the highest correlation with NSCLC phenotype.

Only genes shared across all three datasets were defined as candidate core targets. These candidates were prioritized based on their biological relevance, node degree values from the protein–protein interaction (PPI) network, and existing literature linking them to cancer progression or drug response. The top-ranking targets were selected for subsequent functional enrichment and experimental validation.

To explore interactions among these core targets, the overlapping genes were imported into the STRING database (https://string-db.org/), with “Homo sapiens” selected as the species. A confidence score cutoff of 0.4 was applied to filter both direct and indirect interactions. The resulting PPI network was visualized using Cytoscape, and hub genes were identified based on key topological features, including degree centrality, betweenness centrality, and closeness centrality.

Construction of a Prognostic Model of JWWJD for NSCLC Using Univariate Cox and LASSO Regression Analysis

To develop a prognostic model reflecting the therapeutic effects of JWWJD in NSCLC, survival data including overall survival (OS) time and status were obtained from clinical information of NSCLC patients downloaded from the GEO database. These survival data were combined with the 25 intersected targets derived from JWWJD drug targets, DEGs, and WGCNA key module genes.

Univariate Cox proportional hazards regression analysis was performed using the survival package in R to identify candidate genes associated with patient survival outcomes. Genes with P-values< 0.05 were considered statistically significant and selected for further modeling. Subsequently, Least Absolute Shrinkage and Selection Operator (LASSO) regression was applied to refine the gene set by reducing model complexity and eliminating redundant or collinear variables. The analysis was conducted using the glmnet package in R, with 10-fold cross-validation used to determine the optimal penalty parameter λ that minimized the partial likelihood deviance.16

A risk score formula was constructed based on the coefficients (β values) of the selected genes derived from the LASSO-Cox regression model. Each patient’s risk score was calculated according to the expression levels of these genes and their corresponding weights. Patients were stratified into high-risk and low-risk groups using the median risk score as the cutoff value.

The predictive performance of the model was evaluated using Kaplan–Meier survival analysis and receiver operating characteristic (ROC) curve analysis. The area under the ROC curve (AUC) was calculated to assess the model’s discriminatory power for predicting patient survival. Principal component analysis (PCA) was also conducted to visualize the separation between high- and low-risk groups.

Finally, a nomogram incorporating the risk score and clinicopathological features (eg, age, gender, TNM stage) was constructed using the rms package in R to predict individual survival probabilities.17 Calibration plots were used to evaluate the agreement between predicted and observed outcomes. The stepwise workflow for building the prognostic model is illustrated in Supplementary S1.

Use of Multiple Databases to Evaluate the Expression and Role of Key Gene in NSCLC

The differential expression of the prognostic gene SPP1 across various tumor types and adjacent normal tissues was analyzed using the TIMER2.0 database (http://timer.comp-genomics.org/), which enables systematic comparison of gene expression levels between tumor and matched normal samples in all TCGA cancer types. Expression distributions were visualized using box plots, and statistical significance was calculated by the Wilcoxon signed-rank test.

Genes showing consistent upregulation or downregulation across multiple cancer types were selected for further analysis. For lung cancer-specific investigation, TCGA-LUAD (lung adenocarcinoma) and TCGA-LUSC (lung squamous cell carcinoma) datasets were downloaded from The Cancer Genome Atlas program (https://portal.gdc.cancer.gov/) and used to analyze SPP1 expression in 1003 NSCLC patients. Patients were stratified into high- and low-expression groups based on the median SPP1 expression level. Survival differences between these groups were assessed using Kaplan-Meier survival analysis, and statistical significance was determined by the Log rank test. The association between SPP1 expression and clinicopathological variables—including age, gender, and TNM stage—was evaluated using the chi-square test or Fisher’s exact test, depending on the data type.

SPP1 expression and its correlation with metastasis-related genes (MMP2, MMP9) were analyzed using the GEPIA platform (http://gepia.cancer-pku.cn/). LUAD and LUSC datasets from TCGA were evaluated separately. Pearson correlation analysis was performed, with r > 0.3 and P < 0.05 considered significant.

The CancerSEA database (http://biocc.hrbmu.edu.cn/CancerSEA/) was employed for single cell analysis to further explore the prognostic gene-related specific mechanisms of NSCLC. The association between prognostic gene and lung cancer metastasis was further evaluated using the EMTome database (http://www.emtome.org).

Identification of Active and Blood-Entering Components of JWWJD

JWWJD drug and JWWJD drug-containing serum were thawed at 4°C, followed by vortexing of the sample for 1 min, and evenly mixing. Next, 100µL sample was accurately transferred to a 2-mL centrifuge tube, and 400 µL methanol was then added, followed by vortexing for 1 min. Next, samples were centrifuged for 10 min at 13,780 x g and 4°C, and the supernatant was transferred to a new 2-mL centrifuge tube, concentrated and dried. Next, 150 µL 2-chloro-l-phenylalanine (4 ppm) solution prepared with 80% methanol water was added to redissolve the sample, and the supernatant was the filtered by a 0.22-µm membrane and transferred into the detection bottle for LC-MS detection.

LC analysis was performed on a Vanquish UHPLC System (Thermo Fisher Scientific, Inc). Chromatography was carried out with an ACQUITY UPLC® HSS T3 (2.1x100 mm, 1.8 µm) (Waters Corporation). The column was maintained at 40°C. The flow rate and injection volume were set at 0.3 mL/min and 2 µL, respectively. For LC-electrospray (ESI) (+)-MS analysis, the mobile phases consisted of (B2) 0.1% formic acid in acetonitrile (v/v) and (A2) 0.1% formic acid in water (v/v). Separation was conducted under the following gradient: 0–1 min, 8% B2; 1–8 min, 8–98% B2; 8–10 min, 98% B2; 10–10.1 min, 98–8% B2; and 10.1–12 min, 8% B2. LC-ESI (-)-MS analysis was carried out with (B3) acetonitrile and (A3) ammonium formate (5 mM). Separation was conducted under the following gradient: 0–1 min, 8% B3; 1–8 min, 8–98% B3; 8–10 min, 98% B3; 10–10.1 min, 98–8% B3; and 10.1–12 min, 8% B3. MS detection of metabolites was performed on Q Exactive (Thermo Fisher Scientific, Inc.) with ESI ion source. Simultaneous MS1 and MS/MS (Full MS-ddMS2 mode, data-dependent MS/MS) acquisition was used. The parameters were as follows: Sheath gas pressure, 40 arb; aux gas flow, 10 arb; spray voltage, 3.50 and −2.50 kV for ESI(+) and ESI(-), respectively; capillary temperature, 325°C; MS1 range, mass-to-charge ratio (m/z) 100–1000; MS1 resolving power, 70,000 FWHM; number of data dependent scans per cycle, 10; MS/MS resolving power, 17,500 FWHM; normalized collision energy, 30 eV; and dynamic exclusion time, automatic.

Raw data were firstly converted to mzXML by MSConvert in ProteoWizard software package (v3.0.8789) and processed using R XCMS (v3.12.0) for feature detection, retention time correction and alignment. Key parameters settings were set as follows: Ppm=15, peakwidth=c(5,30), mzdiff=0.01 and method=centWave. Next, data were corrected by the area normalization method to eliminate systematic errors. Metabolites were identified by accuracy mass and MS/MS data, which were matched with HMDB (http://www.hmdb.ca), massbank (http://www.massbank.jp/), LipidMaps (http://www.lipidmaps.org), mzcloud (https://www.mzcloud.org) and the metabolite database built by Panomix Biomedical Tech Co., Ltd. The molecular weight of the metabolites was determined according to the m/z of parent ions in MS data. The molecular formula was predicted by adduct ion, and then matched with the database to realize MS identification of metabolites. MS/MS data from a quantitative table of MS/MS data were matched with fragment ions and other information of each metabolite in the database for MS/MS identification of the metabolites.

Representative positive and negative ion chromatograms and identification spectra of curcumol are provided in Supplementary S2–S6.

Molecular Docking and Molecular Dynamics (MD) Simulation

To further clarify the binding affinity of the active components of JWWJD and core proteins, the two-dimensional structure of the compound curcumol in JWWJD was downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), and the structure of the target protein SPP1 (PDB ID:3cxd) was obtained from the Protein Data Bank (PDB) database (https://pubchem.ncbi.nlm.nih.gov/). Subsequently, we employed AutoDock Tools to perform hydrogenation and charge calculation, and used AutoDock Vina to dock the small molecules curcumol to the SPP1 protein in triplicate. Parallel computing resources were utilized to accelerate the process, and the resulting binding affinities were averaged to generate a consensus prediction. Docking parameters and scoring functions were based on previously established protocol.18 Finally, we applied PyMOL Molecular Graphics System (PyMOL) 3.8 software to visualize their three-dimensional structure and used the optimally docked conformation as input files for subsequent MD simulations.

The GROMACS 2022 program was utilized for MD. Generalized Amber Force Field (GAFF) was used to process the small molecular ligand, while AMBER14SB force field19 and TIP3P water model were applied to process the protein, and the simulation system of the complex was constructed by merging the files of curcumol and SPP1. The analysis was performed under constant temperature, constant pressure and periodic boundary conditions. During the MD simulation process, all the hydrogen bonds were constrained by the LINCS algorithm, with an integration time step of 2 femtoseconds. The energy-minimized structures were subjected to 100 picoseconds canonical ensemble (NVT) and 100 picoseconds isothermal-isobaric ensemble (NPT) simulations at 298 K and 1 bar pressure. Following the completion of simulations, the simulated trajectories were analyzed using Variational Mode Decomposition (VMD) and PyMOL. Furthermore, the binding free energy between curcumol and SPP1 was calculated using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method via the g_mmpbsa plugin, providing a quantitative assessment of the thermodynamic stability of the complex.20

Surface Plasmon Resonance (SPR) Analysis

The binding affinity between curcumol and SPP1 (Sinobiological, Beijing, China) was performed using the Biacore 1K Cytiva system (GE Healthcare, Sweden), following standard SPR protocols.21 30–50ug/mL SPP1 protein diluted by Acetate 4.0 Cytiva (GE Healthcare, Sweden) was injected three-times and immobilized to the activated SPR chip CM5 (GE Healthcare, Sweden) for 600s+600s. Various concentrations of curcumol were run over the SPR using the running buffer containing 1.76 mM KH2PO4, 8.1 mM Na2HPO4, 137 mM NaCl, 2.68 mM KCl, and 0.05% Tween-20. Binding and dissociation rates were measured at a flow rate of 30 µL/min. The injection of curcumol was performed for 90 seconds, followed by a flow with ligand-free buffer to analyze the dissociation for 180 seconds. Curves were corrected for nonspecific ligand binding by subtracting the signal obtained for the the negative control flow cell. The equilibrium dissociation constant (Kd) was derived from a simple 1:1 interaction model using the BIAevaluation Software (Biacore AB, Uppsala, Sweden).

Cells and Plasmid Transfection

Human NSCLC cell lines NCI-A549, NCI-H23, NCI-H1299, NCI-H520, 95-D and SK-MES-1, human bronchial epithelial BEAS-2B cells, and Human Umbilical Vein Endothelial Cells (HUVECs) were purchased from Pricella (Wuhan, Hubei Province, China). SC-79 was purchased from Beyotime (Shanghai, China). Cis-platinum (CDDP) (the positive control) was purchased from Solarbio (Beijing, China). Cells were maintained in DMEM or DMEM/F12 from Servicebio (Wuhan, Hubei Province, China) containing 10% FBS (MeiSenCTCC, Shanghai, China), and were incubated at 37°C in a humidified incubator containing 5% carbon dioxide.

OptiMEM (Gibco, New York, USA) and Lipofectamine 3000 transfection reagent (Invitrogen, California, USA) were used according to the manufacturer’s instructions. NCI-A549 and NCI-H23 cells were transfected with pcDNA3.1(+)-SPP1 expression vector and pcDNA3.1(+)-control vector from Sangon Biotech (Shanghai, China). After 6 h of transfection, the fresh medium was replaced, and the cells were incubated overnight. Cells were treated with corresponding drugs the following day. Transfection efficiency was verified by reverse transcription-quantitative PCR (RT-qPCR) and WB.

Cell Viability Assay

NCI-A549, NCI-H23 and BEAS-2B cells were harvested with trypsin and seeded into 96-well plates at a density of 3.5×103 cells per well. After overnight incubation to allow cell attachment, the culture medium was replaced with fresh medium containing either JWWJD-treated rat serum (final concentrations: 0%, 5%, 10%, 15%, or 20%) or curcumol (final concentrations: 0, 50, 100, 200, 300, 400, 500, or 600 µmol/L). For the JWWJD group, cells were incubated for 24, 48, and 72 hours, while for the curcumol group, cells were incubated for 24 and 48 hours.

Following treatment, 10 µL of Cell Counting Kit-8 (CCK-8) solution (Beyotime, Shanghai, China) was added to each well, and the plate was incubated for an additional 1–2 hours at 37 °C. The absorbance at 450 nm was measured using a microplate reader (PerkinElmer, Massachusetts, USA). Cell viability was calculated as a percentage relative to the control group.

Colony Formation Assay

After digestion with trypsin, NCI-A549 and NCI-H23 cells were counted, inoculated into 12-well plates at a density of 6×103 cells/well and incubated overnight. After 24 h of treatment with different concentrations of curcumol (0, 100, 200 or 400 µM) and 15 µg/mL of CDDP, the culture medium was replaced with fresh medium for additional 2 weeks. Cells were washed with PBS, immobilized with 4% paraformaldehyde fixative, stained with 0.1% crystal violet (Servicebio, Wuhan, Hubei Province, China) for 15 min, and washed twice with PBS. Next, the number of colonies with >50 cells were photographed and counted.

EdU Proliferation Assay

NCI-A549 and NCI-H23 cells were inoculated in 48-well plates overnight and treated with different concentrations of curcumol (0, 100, 200 or 400 µM) and 15 µg/mL of CDDP for 48 h. According to the manufacturer’s instructions of BeyoClick™ EdU Cell Proliferation Kit with Alexa Fluor 594 (Beyotime, Shanghai, China), 10 µL of Edu was added to each well to continue incubation for 2 h, followed by immobilization with 4% paraformaldehyde fixative and permeabilization with 0.3% Triton X-100. Next, Click-iT® reaction mixture was added. After incubating for 30 min in the dark, the nuclei were stained with Hoechst 33258 (Servicebio, Wuhan, Hubei Province, China). The stained cells were captured by fluorescence microscopy (Olympus Corporation, Tokyo, Japan).

Wound Healing and Transwell Assay

NCI-A549 and NCI-H23 were inoculated into a 12-well plate. When the cells reached 70% confluency, a straight scratch was made in the cell monolayer with a 200-µL pipette tip. After removing the floating cells by washing with PBS, transfection and drug treatment were performed. Cell migration rate was calculated by measuring the scratch width under a inverted microscope (Olympus Corporation, Tokyo, Japan) at 0, 24 and 48 h.

In the Transwell migration assay, 600 µL DMEM containing 10% FBS was added to the lower chamber. NCI-A549 and NCI-H23 cells, resuspended in 200 µL serum-free culture medium after 48 h of plasmid transfection or drug treatment, were uniformly inoculated into the upper chamber at a density of 3.5×104 cells per well. After culturing in an incubator for 24 h, cells were washed with PBS twice, fixed with 4% paraformaldehyde for 15 min and stained with 0.1% crystal violet for 15 min. Cells on the upper surface of the membrane were carefully wiped away with sterile cotton and captured by microscopy. For the invasion assay, the serum-free medium was diluted Matrix-Gel™Basement Membrane Matrix (Beyotime, Shanghai, China) on ice in a ratio of 1:3 was seeded in 24-well plates, followed by incubation at 37°C for 3h. The subsequent procedure is the same as the Transwell cell migration experiment.

Cell Apoptosis Assay

NCI-A549 and NCI-H23 cells on 6-well plates were treated with different concentrations of curcumol (0, 100, 200 or 400 µM) and 15 µg/mL of CDDP for 48 h. Upon digestion with EDTA-free pancreatic enzymes (Servicebio, Wuhan, Hubei Province, China), cells were washed with PBS and resuspended in 500 µL 1X binding buffer. Then 5 µL annexin V-FITC and 5 µL propidium iodide (PI) (Apexbio, Texas, USA) were added to the cells and mixed. After incubation at room temperature for 15 min in the dark, the apoptosis of NSCLC cells was detected by CytoFLEX flow cytometer (FACSAria™ II, New Jersey, USA).

Hoechst/PI Staining Assay

NCI-A549 and NCI-H23 cells were inoculated on 48-well plates overnight, treated with different concentrations of curcumol (0, 100, 200 or 400 µM) and 15 µg/mL of CDDP for 48 h, and washed twice with PBS. After staining with PI and Hoechst 33342 (Beyotime, Shanghai, China) in the dark at room temperature for 15 min and then washing with PBS twice, images of stained cells were captured using fluorescence microscope.

HUVECs Tube Formation Assay

NCI-A549 cells were cultured in 6-well plates overnight, then supplemented with curcumol (0, 400 µM) or transfected with SPP1 plasmid. After 24 h, the conditioned medium was harvested by taking supernatant after centrifugation. HUVECs were suspended in 100 µL with different conditioned medium and seeded into Matrigel-coated 96-well plates at a density of 3×104 cells per well for 24 h. The angiogenesis of HUVECs was observed by obtaining images under an inverted microscope. The number of branches and junction points, and the length of the total segments were calculated.

Immunofluorescence (IF) Staining Assay

NCI-A549 and NCI-H23 cells were inoculated on a microscope slide overnight and then treated with curcumol (0, 400 µM) for 48 h. Next, cells were fixed with 4% paraformaldehyde, permeabilized with 0.3% Triton X-100, and incubated with 5% Bovine Serum Albumin (BSA) at room temperature for 60 min. Subsequently, cells were incubated with a primary antibody against SPP1 overnight at 4°C. Upon washing with PBS-Tween 20 (PBST) three times, CoraLite488-conjugated Goat Anti-Rabbit IgG (H+L) antibody was added and incubated for 1 h at room temperature in the dark. After washing with PBST three times, 100 µL Hoechst 33258 was added for nuclear staining. The stained cells were captured by fluorescence microscopy.

RNA Extraction and RT-qPCR Analysis

NCI-A549 and NCI-H23 cells were digested and inoculated on a 6-well plate, and treated with drugs and plasmids for 24 h. Total RNA was extracted by TRIzol reagent (Servicebio, Wuhan, Hubei Province, China), and its concentration was measured by NanoDrop ONE© (Thermo, Massachusetts, USA). cDNA was synthesized with PrimeScript RT Master Mix (Takara, Kyoto, JPN). A StepOnePlus System (Thermo Fisher Scientific, Waltham, Massachusetts, USA) was used for fluorescence qPCR amplification, and the relative expression level of the target gene SPP1 was calculated by relative quantification using the 2−ΔΔCq method, With β-actin as an internal standard control for mRNA expression. The primers were synthesized by Sangon Biotech (Shanghai, China) and the sequences were designed as follows. SPP1: Forward: TTGGGACAGCCGTGGGAAGG, Reverse: TCACATCGGAATGCTCATTGCTCTC; β-actin: Forward: CATGTACGTTGCTATCCAGGC, Reverse: CTCCTTAATGTCACGCACGAT.

WB Analysis

Total protein was extracted with Radio Immunoprecipitation Assay (RIPA) lysis buffer containing Phenylmethanesulfonyl fluoride (PMSF) (Servicebio, Wuhan, Hubei Province, China). Protein concentration was measured by using the Bicinchoninic Acid Assay (BCA) Kit (Beyotime, Shanghai, China). Equal quantities of proteins were separated with 10 or 15% Sodium Dodecyl Sulfate-Polyacrylamide Gel Electrophoresis (SDS-PAGE) (Epizyme, Shanghai, China) and transferred to nitrocellulose membranes (MilliporeSigma, Louis, USA). Upon blocking with 5% nonfat powder milk at room temperature for 1 h, the membranes were incubated overnight at 4°C with primary antibodies against SPP1, PI3 Kinase p85 alpha, phosphorylated (p)-PI3-kinase p85α (Tyr607), AKT1, p-AKT1 (T450), MMP2, MMP9, VEGF and β-Actin. After washing with Tris-Buffered Saline with Tween (TBST) three times, the secondary antibody Goat Anti-Rabbit IgG (H+L) Horseradish Peroxidase (HRP) was added and incubated at room temperature for 1h. After three washes with TBST, the gray value of each protein band was determined using a gel imaging system (ChemiScope 6200, Shanghai, China). β-Actin was used as the reference protein, and ImageJ software win64.exe (National Institutes of Health, USA) was used to compare the relative expression changes of target proteins. Information of the antibodies involved is listed in Table 2.

Table 2 Information on Antibodies Involved

Cellular Thermal Shift Assay (CETSA)

NCI-A549 and NCI-H23 cells were lysed with RIPA lysis buffer and aliquoted into PCR tubes. Upon treatment with curcumol (400 µM) or an equal volume of Dimethyl Sulfoxide (DMSO) and placed on a horizontal shaker for 1 h at room temperature, the solution was heated at a specified temperature (37–73°C) for 3 min using a StepOnePlus System and then cooled at 4°C for 3 min. After centrifugation for 15 min at 12,000 x g and 4°C, the supernatant was extracted, and WB was performed.

Animal Experiments

A total of 48 BALB/c-nu mice (4 weeks old, 18–22g) were provided by Shulaibao (Wuhan) BioteChinaology Co., Ltd (Wuhan, Hubei Province, China). After 7 days of adaptive feeding, 5×105 NCI-A549 cells were injected into the right flanks of the mice. After 10 days, mice were randomly assigned to six groups using a computer-generated randomization sequence: control (Con, n=12, saline 0.1 mL/kg, i.g. daily), cisplatin (CDDP, n=12, i.p. every two days), low-dose curcumol (Cur-L n=6, 50mg/kg, i.g. daily), high-dose curcumol (Cur-H n=6, 100 mg/kg, i.g. daily), low-dose JWWJD (JWWJD-L n=6, 11 g crude drug/kg, i.g. daily) and high-dose JWWJD (JWWJD-H n=6, 22 g crude drug/kg, i.g. daily). All treatments were administered continuously for 14 consecutive days. Tumor size and body weight were recorded every 2 days during the treatment period. On day 15, 12 h after the final administration, all mice were anesthetized and then sacrificed by cervical dislocation. The tumors were removed, photographed, weighed, and used for subsequent transcriptomic, IHC and WB analysis.

H&E Staining and IHC Analysis

Nude mouse lungs, livers and kidneys were fixed with 4% paraformaldehyde and paraffin embedded, and 5-µm sections were stained with H&E. The structural changes of the tissues in the sections were examined under an upright white light microscope (Eclipse Ci-L, Tokyo, Japan).

All clinical NSCLC specimens used in this study were consistent with the corresponding histopathological diagnosis. NSCLC and normal adjacent lung tissues were fixed in 10% neutral buffered formalin and mounted on paraffin sections. A primary antibody against SPP1 was added. The intensity of staining was evaluated blinded by two independent pathologists using the following criteria: 0, negative; 1, low; 2, medium; and 3, high. The staining degree was scored as 0, 0% stained; 1, 1–25% stained; 2, 26–50% stained; and 3, 51–100% stained. In total, 5 random fields were evaluated under a light microscope.

Sections were deparaffinized and incubated with 3% H2O2 for 30 min. Next, 3% BSA was added for blocking, followed by the addition of the following primary antibodies and incubated overnight at 4°C: SPP1, MMP2, MMP9 and VEGF. The following day, samples were incubated with a Biotin-labeled Goat Anti-Rabbit IgG(H+L) for 1 h at room temperature. 3,3’-Diaminobenzidine (DAB) and hematoxylin (Beyotime, Shanghai, China) were used for staining. Sections were photographed under a light microscope.

RNA Sequencing (RNA-Seq)

Total RNA was extracted from tumor tissues of the Control and High-dose JWWJD groups, and its quality and purity were detected. The complementary DNA (cDNA) libraries were sequenced on the Illumina sequencing platform by Metware Biotechnology (Wuhan, Hubei Province, China). DESeq2 was used for differential gene expression analysis between two groups, and the Benjamini-Hochberg correction was applied to P-value. Fold change≥1.4 or Fold change≤0.714 and FDR< 0.05 were used as thresholds for significant differential expression.

Statistical Analysis

Experimental data were presented as the mean ± standard deviation (SD) and analyzed using GraphPad Prism version 10.1.2 (GraphPad Software, La Jolla, CA, USA). For comparisons between two independent groups, an unpaired Student’s t-test was used; for paired samples, a paired t-test was applied. When comparing more than two groups, one-way analysis of variance (ANOVA) followed by Dunnett’s post-hoc test was performed to determine statistical significance. Normality and homogeneity of variances were assessed before selecting parametric or non-parametric tests. The association between SPP1 expression levels and clinicopathological parameters (eg, gender, age, TNM stage) was evaluated using either the chi-square test or Fisher’s exact test, statistical analyses were conducted using SPSS version 26.0. Survival analysis was conducted using the Kaplan–Meier method, and differences in survival curves were assessed with the Log rank test. A P-value of less than 0.05 was considered a statistically significant difference (*P < 0.05, **P < 0.01,***P < 0.001; ns indicates not significant). A detailed summary of the statistical methods used for each type of analysis is provided in Supplementary S7.

ResultsJWWJD Inhibits the Tumor Growth in vivo, and Inhibits the Proliferation, Migration and Invasion of NSCLC in vitro

To confirm whether JWWJD inhibits the tumor growth in vivo, tumor-bearing nude mice were gavaged with 0.2 mL JWWJD daily, while control mice were administered the same volume of saline. After 14 days, the tumor volume and weight of mice in the JWWJD group were significantly lower than those in the control group, particularly in the high-dose JWWJD group. There was no significant change in body weight among the groups (Figure 2A–D).

Figure 2 Continued.

Figure 2 Continued.

Figure 2 Effect of JWWJD against NSCLC in vitro and in vivo. (A) Subcutaneous tumor photos of control, CDDP, low-dose JWWJD, high-dose JWWJD groups. (B) Tumor volume of mice in each group was recorded. (C) Tumor tissues were harvested, and the weight were measured. (D) Body weights in each animal group were recorded. (E) Compound-target network of JWWJD. (F) Top 10 hub genes of JWWJD. (G) KEGG enrichment analysis of the detailed pathways by network pharmacology. (H–J) Transcriptomics analysis. Volcano plot of DEGs with Fold change≥1.4 or Fold change≤0.714 and FDR< 0.05. Red indicates upregulation, while green indicates downregulation and gray not significant. GO enrichment and KEGG pathway enrichment analysis of DEGs. (K) CCK8 assay was used to detect the viability of JWWJD drug containing serum-treated NSCLC cells for 24, 48 and 72 h. The effect of JWWJD drug containing serum and CDDP on (L) the ability of cells to migrate by wound healing assay, Scale bar=200 µM, (M) the ability of cells to migrate and invade by Transwell assay, Scale bar=200 µM, (N) cell apoptosis by flow cytometry in annexin V-FITC/PI staining assay. (O) WB assay evaluated the expression levels of PI3K, p-PI3K, AKT1, p-AKT1, MMP2, MMP9, SPP1 and VEGF. (P) RT-PCR and WB assays were applied to confirm that the SPP1 had been overexpressed and the effect of SC79. (Q and R) JWWJD significantly inhibited migration and invasion in NSCLC cells, and the inhibitory effect were partially reversed by SC79 (10 µM) or SPP1 overexpression, which was detected by wound healing and Transwell assays, respectively. *P < 0.05, **P <0.01, ***P <0.001 vs control; ns, not significant. All data are shown as the mean ± SD from five independent experiments.

According to OB > 30% and DL > 0.18, active compounds and potential targets of JWWJD were identified using the TCMSP and BATMAN databases. The Drug-Target network was constructed through Cytoscape3.7.2 (Figure 2E). The cytoHubba result suggested that AKT1 was the hub gene of JWWJD (Figure 2F). Network pharmacology and transcriptomic results suggested that JWWJD inhibited NSCLC cell wound healing, cell junction and other migration-related phenotypes, particularly through suppression of the PI3K/AKT signaling pathway and downregulation of SPP1 expression (Figure 2G–J).

Compared with drug-free serum of control group, the proliferation of NCI-A549 and NCI-H23 cells was significantly inhibited after treatment with JWWJD drug-containing serum for 24, 48 and 72 h. The inhibitory effect of JWWJD on NSCLC cell proliferation increased with increasing in drug-containing serum concentration, with 20% drug-containing serum exhibiting the best inhibitory effect on tumor cell proliferation. However, there was no significant effect on the proliferation of BEAS-2B cells (Figure 2K). Compared with drug-free serum of the control group, 15 and 20% JWWJD-containing serum could significantly inhibit wound healing and significantly reduce the number of cells invading the lower chamber, suggesting that JWWJD drug-containing serum inhibited the migration and invasion of NCI-A549 and NCI-H23 cells (Figure 2L and M). Compared with drug-free serum of the control group, the apoptosis rate in different concentrations of drug-containing serum groups was significantly increased, suggesting that JWWJD induced the apoptosis of NSCLC cells in a serum concentration-dependent manner (Figure 2N). However, compared with the positive control CDDP, JWWJD was less effective in inducing NSCLC apoptosis. WB experiments showed that JWWJD drug-containing serum could not only reduce the expression of MMP2, MMP9 and VEGF, but also reduce the expression of proteins associated with the PI3K/AKT1 pathway and SPP1 target (Figure 2O).

The transfection efficiency and expression of SPP1 were determined by RT-qPCR and Western blotting analysis (Figure 2P). In the rescue experiments, after adding a PI3K/AKT1 pathway agonist (SC79) or overexpressing SPP1, it was found that the inhibitory effect of JWWJD drug-containing serum on the metastasis and invasion of NCI-A549 and NCI-H23 cells could be reversed (Figure 2Q and R). This further verified that JWWJD drug-containing serum inhibits the metastasis and invasion of NCI-A549 and NCI-H23 cells by inhibiting the PI3K/AKT1 pathway and lowering the expression of SPP1.

Use of Network Pharmacology and Bioinformatics Analysis to Identify the Key Genes of JWWJD in the Treatment of NSCLC, and Construction of the Related Prognostic Model

A total of 538 NSCLC genes in the training set were screened by differential expression analysis. In the heat map, the normal and tumor groups were distinctly separated, with red indicating high gene expression and green indicating low gene expression (Figure 3A).

Figure 3 Continued.

Figure 3 Construction and validation of a prognostic model of JWWJD for NSCLC. (A) DEGs in NSCLC obtained from the training set. (B and C) WGCNA was employed to obtain genes related to the occurrence and development of NSCLC in the training set. (D) The intersection of the targets of JWWJD and DEGs. (E) A forest plot of 15 JWWJD-DEGs associated with NSCLC survival was constructed by univariate Cox analysis. (F and G) LASSO analysis was used to remove redundant factors and construct the prognostic model. (H and I) PCA plots based on the risk score in the training set (H) and the test set (I). (J and K) Kaplan–Meier survival analysis of overall survival in patients stratified into high-risk and low-risk groups based on the risk score, shown for the training set (J) and the test set (K). (L and M) ROC curves for 3-year overall survival prediction based on the risk signature in the training set (L) and the test set (M). (N and O) Distribution of risk scores and survival status in the high-risk (red) and low-risk (green) groups for the training set (N) and the test set (O). (P and Q) Heat map of the expression of risk genes associated with the risk score in the training set (P) and the test set (Q). (R–U) Forest plots of univariate and multivariate Cox regression analyses showing the association between clinical characteristics and survival outcomes in the training set (R and S) and the test set (T and U). (V and W) Heat map of clinical parameters correlated with risk genes in the training set (V) and the test set (W). (X) A calibration plot was constructed to assess the agreement between predicted and observed survival probabilities for the 1-, 3-, and 5-year overall survival in the training set. (Y) The nomogram integrates the prognostic risk score with clinical variables, including gender, age, smoking status, and stage, to predict individual survival probabilities in the training set. In the nomogram, variable weights are indicated as follows: *P < 0.05, **P < 0.01, and ***P < 0.001, reflecting the relative contribution of each factor to the prediction of survival outcomes.

In the WGCNA of the training set, six gene modules were identified through average linkage hierarchical clustering and further merged based on high topological similarity. To determine the relevance of each module to NSCLC, module-trait relationships were constructed using Pearson correlation coefficients. The module significance (MS), defined as the -log10-transformed p-value of the correlation between gene expression profiles and the NSCLC phenotype across all samples, was calculated for each module to quantify its overall association with NSCLC. Among the six modules, the turquoise module exhibited the highest correlation with NSCLC (cor = 0.53, P = 6e-24) and also showed the highest module significance, indicating a strong and statistically significant relationship with the disease state (Figure 3B and C). Therefore, the turquoise module was selected as the key module for downstream analysis.

The Venn package in R language was utilized to screen DEGs of NSCLC, and the genes related to the development of NSCLC identified by WGCNA were intersected with the target of JWWJD. A total of 25 key targets for the treatment of NSCLC by JWWJD were obtained (Figure 3D).

Univariate Cox analysis was performed on 25 key targets of JWWJD in the treatment of NSCLC. The results showed that 15 genes were significantly associated with survival (Figure 3E). Further, these 15 genes were used for the LASSO regression to eliminate redundant factors, and a prognostic model was constructed with 7 genes (Figure 3F and G). The model formula was as follows: (0.0981032825803083) × ALDH1A2 + (0.0983136632977013) × BIRC5+(−0.00839899177335228) × CAV1 + (0.101549134292197) ×IGFBP3 + (0.205372480091101) ×SPP1 + (0.0502142406197653) × TCN1 + (−0.0340760580748295) × THBD. Patients in the training and test sets were divided into high and low risk groups based on median risk scores. By principal component analysis (PCA), there was a significant difference in the distribution direction of patients in the high and low risk groups. Kaplan-Meier survival analysis revealed significant differences between high and low risk groups in both the training and test sets. To further evaluate the accuracy of the model, the model in the training and test sets was further verified by univariate analysis, multivariate analysis, survival state diagram, Area Under the Curve (AUC) and other methods. All these results suggested that this prognostic model has good predictive accuracy for the prognosis of patients with NSCLC (Figure 3H–U). Furthermore, through clinically relevant heat map and nomogram, it was visually observed that the risk score of the prognostic model correlated with the clinical stage and smoking status of patients with NSCLC (Figure 3V–Y).

Associations Between SPP1 Expression and Clinical Prognosis, Clinicopathological Characteristics and Roles of SPP1 in NSCLC

Compared with normal tissues, the expression of SPP1 was elevated in pancarcinoma tissues, particularly in NSCLC tissues (Figure 4A). Analysis of TCGA data revealed that SPP1 expression in NSCLC was significantly higher than in adjacent normal tissues and showed a significant association with the N stage (Figure 4B and C). Bioinformatics analyses demonstrated that patients with NSCLC and high SPP1 expression had significantly shorter overall survival compared to those with low SPP1 expression (Figure 4D). Moreover, SPP1 expression levels were significantly correlated with the N stage (P = 0.002), and varied across different lymph node metastasis stages (N0–N3), suggesting a potential role of SPP1 in tumor progression and lymph node metastasis (Figure 4E).

Figure 4 Continued.

Figure 4 Expression and role of SPP1 in NSCLC. (A) Pan-cancer SPP1 expression. (B and C) Differential expression of SPP1 in NSCLC and normal tissues was compared according to TCGA database. (D) Survival analyses of SPP1 regarding overall survival using the TCGA database. (E) Heat map of clinical correlations of the expression of SPP1 according to TCGA database. (F) The expression of SPP1 in 20 pairs of NSCLC tissues was determined by Western blotting. (G) SPP1 expression in 50 NSCLC tissues, as measured by IHC. Scale bar, 100 µm. (H) The Kaplan-Meier curves analysis shows the relationship between the expression of SPP1 and overall survival in 50 NSCLC patients. (I) SPP1 expression in human bronchial epithelioid cell and in human NSCLC cell lines, as determined by Western blotting. (J) T-SNE cluster analysis shows the distribution of cells, with every point representing a single cell, and the color of the point representing the expression level of SPP1 in the cell. (K) Correlations between SPP1 and functional states in different single-cell datasets in NSCLC. (L) Functional status profile showcasing the diverse activity of function states of NSCLC cells. (M) SPP1 and MMP2 or MMP9 correlation analysis. (N) SPP1 expression in primary and metastatic tumors. (O) CRISPR-Cas9 genetic perturbation reagents were employed to knockout SPP1 and to identify that it affected cell survival. *P < 0.05, **P <0.01, ***P <0.001 vs the control group. Data are presented as means ± SD.

In clinical fresh tumor samples, Western blot analysis confirmed that SPP1 expression was significantly higher in tumor tissues than in paired adjacent normal tissues (Figure 4F). Immunohistochemical (IHC) staining further revealed strong cytoplasmic expression of SPP1 in NSCLC specimens, whereas only weak or undetectable signals were observed in matched adjacent normal tissues (Figure 4G). Clinical survival analysis further validated that patients with high SPP1 expression had significantly reduced overall survival compared to those with low expression (Figure 4H). These data clearly indicated that SPP1 was upregulated in NSCLC. The clinical and pathological characteristics of patients with NSCLC are summarized in Table 3. SPP1 expression showed a significant correlation with N stage (P = 0.002) and varied among the different NSCLC N stages (N0, N1, N2 and N3). At the cellular level, Western blotting showed that SPP1 was highly expressed in NSCLC cell lines, especially NCI-A549 and NCI-H23, compared with the bronchial epithelial cell line BEAS-2B (Figure 4I).

Table 3 Correlation Between SPP1 and Clinicopathological Factors in NSCLC

To explore the biological functions of SPP1, single-cell analysis was performed using the CancerSEA database. The T-SNE plot illustrated the distribution of SPP1 expression across individual cells (Figure 4J), and the correlation heat map revealed that SPP1 expression was positively associated with key cellular processes such as epithelial-mesenchymal transition (EMT), metastasis, angiogenesis, and invasion (Figure 4K a