Network Pharmacology with Metabolomics Study to Reveal the Mechanisms of Bushen Huoxue Formula in Intervertebral Disc Degeneration Treatment

Introduction

Intervertebral disc degeneration (IVDD) is a multifaceted pathophysiological process involving genetic, mechanical, immune, and metabolic factors. Its pathological changes manifest through a reduction in nucleus pulposus (NP) cell quantity, diminished proliferation activity, and reduced extracellular matrix synthesis.^1,2 IVDD stands as a major contributor to lower back pain, often culminating in disc herniation, spinal stenosis, and other degenerative disc conditions, profoundly impacting patients’ quality of life and productivity and imposing a significant socioeconomic burden. Current clinical approaches for IVDD primarily focus on symptom alleviation, predominantly employing nonsteroidal anti-inflammatory drugs (NSAIDs).^3,4 However, both conservative and surgical treatments fall short of completely reversing IVDD’s pathological changes.⁵

Traditional Chinese medicine (TCM) is characterized by its “multi-component, multichannel” nature, rendering it well-suited for personalized treatment and early intervention.⁶ TCM philosophy emphasizes the enhancement of renal function to effectively treat various disorders related to the lower back. The Bushen Huoxue Formula (BSHXF) stands as a classic prescription known to enhance renal function, promote blood circulation, and strengthen bones. Modified BSHXF, a frequently prescribed variant, demonstrates exceptional clinical efficacy in managing lower back pain.⁷ Recent research validates the ability of BSHXF, including its modified forms, to decelerate the progression of musculoskeletal degenerative diseases, such as osteoarthritis and IVDD.^7,8 Besides, Zhu et al found that BSHXF achieved a good efficacy in reliving lower back pain and improving lumbar vertebral function.⁹ Zhan et al verified the efficacy and safety of BSHXF through a double-blind, randomized, placebo-controlled trial, which furtherly demonstrated that BSHXF combined with lumbar traction can significantly improve the clinical symptoms and lumbar mobility in discogenic lower back pain patients.¹⁰ However, the diverse components of TCM and the intricacies of the human body have limited in-depth systematic research in this area.

Network pharmacology represents a promising avenue for investigating the therapeutic mechanisms of TCM formulations.¹¹ Network analysis results reveal that different proteins or genes can regulate the same diseases, and certain proteins can regulate a variety of diseases.¹² Moreover, this approach allows for the analysis of the collective mechanisms of action of the multiple components present in Chinese herbal compounds.^13,14 In our study, the modified BSHXF primarily comprises six herbs: Radix Salviae (DS), Chaenomeles sinensis (Thouin) Koehne (MG), Achyranthis Bidentatae Radix (NX), Cuscutae Semen (TSZ), Radix Clematidis (WLX), and Eucommiae Cortex (DZ). Despite the previous research on the chemical constituents of BSHXF, the key active components and signaling pathways critical for its efficacy in IVDD treatment remain elusive.

Metabolomics, a sub-discipline of systems biology, examines the overall functional alterations in organisms resulting from various factors. By analyzing biological samples from model animals or human for metabolite profiling, this approach facilitates the identification of specific marker metabolites, allowing for comparisons of metabolite changes in the body pre- and post-treatment with the Chinese herbal formula.

In this study, our goal is to verify the molecular mechanisms and pathways of the main bioactive ingredients within BSHXF in the treatment of IVDD, employing a dual approach of network pharmacology and metabolomics. This strategic integration compensates for the experimental validation gaps in network pharmacology and the absence of upstream molecular mechanisms and drug-binding targets in metabolomics. Our research endeavors to shed light on the therapeutic rationale underpinning BSHXF for IVDD treatment. Figure 1 illustrates the flowchart detailing the progression of this research.

Figure 1 Flowchart to explore the mechanism of BSHXF in treating IVDD.

Materials and Methods

Network Pharmacology Analysis

Screening Active Components in BSHXF and Prediction Their Targets

We sourced herbs and chemical ingredients from the Traditional Chinese Medicine Systems Pharmacology database (TCMSP), a comprehensive herbal encyclopedia providing herbal ingredient information.^15,16 Keywords like “Dan Shen”, “Mu Gua”, “Niu Xi”, “Tu Si Zi”, “Wei Ling Xian”, and “Du Zhong” were used to identify targets related to BSHXF. Candidate compounds had to meet the criteria of oral bioavailability (OB ≥ 30%) and drug-like properties (DL ≥ 0.18).¹⁷ All target names were standardized using the UniProt database (https://sparql.uniprot.org/).¹⁸

Screening of Related Targets for IVDD and Common Gene Set

Target genes associated with IVDD were screened from Genecards (https://www.genecards.org/), Online Mendelian Inheritance in Man (OMIM, https://www.genecards.org/), and PharmGKB (https://www.pharmgkb.org/) using the search term “Intervertebral disc degeneration”.^19–21 The inclusion criterion was a relevance score of ≥17.0. After eliminating duplicates, potential IVDD targets were identified.

Construction of Ingredient-Target-Disease Network

The common targets of active ingredients from BSHXF and potential targets for IVDD were visualized in an ingredients-target-disease interaction network using Cytoscape 3.8 software (http://www.cytoscape.org/). The network integrated topological parameters, and the importance of ingredients and targets was determined by their degree, representing the total number of routes related to the node by other nodes. Higher degree values indicated higher importance.

Protein-Protein Interaction (PPI) Network Construction and Core Gene Selection

To elucidate interactions between IVDD-related targets and potential BSHXF targets, common targets were analyzed using the String database (https://string-db.org/). The PPI network was constructed via Cytoscape 3.8 software. Core proteins were identified using CytoNCA software based on betweenness centrality (BC), closeness centrality (CC), degree centrality (DC), eigenvector centrality (EC), local average connectivity-based method (LAC), and network centrality (NC), with values greater than or equal to the median value.²²

Pathway Enrichment Analysis

Major biological processes and pathways were analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) platforms, respectively. The GO and KEGG enrichment analysis results were selected based on a critical criterion of P < 0.05, and the first 30 pathways were screened out. GO analysis describes gene product functions, including cell function, molecular function, and biological function, while KEGG analysis elucidates the location and role of the target in related biological signal pathways.

Preparation of BSHXF Solution and Serum

To prepare the BSHXF solution for animal experiments, Radix Salviae (DS), Chaenomeles sinensis Koehne (MG), Achyranthis Bidentatae Radix (NX), Cuscutae Semen (TSZ), Radix Clematidis (WLX), and Eucommiae Cortex (DZ) were sourced from the pharmacy department of Shanghai Guanghua Hospital. The herbs were mixed in the proportion of 1.5:1:1:1.2:1:0.9 (w/w) and immersed in 8 volumes of distilled water (v/m) for 1 hour. They were then decocted twice for 1.5 hours each time. Subsequently, the water-based decoction was concentrated to a concentration of 0.6 g/mL and stored at −20°C until use.

For the preparation of BSHXF-containing serum, mice were obtained from Shanghai Slac Laboratory Animal Co. Ltd. These mice were housed in a controlled environment with a 12:12-hour light-dark cycle, appropriate temperature, and free access to standard rodent chow and tap water. Mice were gavaged with BSHXF solution (BSHXF containing serum, n=3) or distilled water (Control serum, n=3) for 7 days. On the 7th day, abdominal aortic blood was collected using a disposable negative pressure blood collection vessel and blood collection needle. The serum was extracted through centrifugation at 3000 rpm for 15 minutes. Mice BSHXF-containing serum and control serum were collected separately, sterilized in a 55 °C water bath for 30 minutes, and stored at −20 °C.

Animal Modeling and Experimental Grouping

This study was complied with the Declaration of Helsinki and approved by the Ethics Committee of Shanghai Guanghua Hospital of Integrated Traditional Chinese and Western Medicine [Ethics approval number: 2023-K-09]. Besides, all laboratory procedures were performed following the Regulations for the Administration of Affairs Concerning Experimental Animals approved by the State Council of the People’s Republic of China. Male mice, approximately 12 weeks old and weighing around 20 to 25 g, were used in this study. An intraperitoneal anesthetic injection of a combination of 10 mg/kg xylazine and 60 mg/kg ketamine hydrochloride was administered to anesthetize the mice. They were then placed in a prone position on the experiment table, and a sagittal small skin incision was made from Co5 to Co7 to aid in locating the disc position for needle insertion in the tail. The Co5-Co6 discs were punctured using a syringe needle according to a previously described method.²³ The needle was inserted into the Co5-Co6 discs in a vertical direction and then rotated in the axial direction by 180°, being held for 5 seconds. The puncture was made parallel to the endplates through the annulus fibrosus (AF) into the NP tissue, with a locking forceps clamped at 2 mm controlling the depth of penetration. The Co6-Co7 segment was left undisturbed and served as a contrast segment. All muscles and soft tissues under the skin were retracted to the same extent, and the coccygeal discs were exposed under aseptic conditions during the operation.

The mice were randomly divided into three groups: the sham group, the model group, and the BSHXF group. Mice in the sham group did not undergo needle puncture, while those in the other two groups received needle puncture modeling. The BSHXF group received continuous gavage after modeling for 4 weeks at a dosage of 0.7 mL/100 g body weight, determined based on the dosage conversion formula. The mice in the sham group and the model group were fed the same dosage of 0.9% normal saline as a control. The experimenter conducting the operation was blinded to the animal grouping.

In vitro NP Cell Culture and Experimental Grouping

The NP tissue of the coccygeal disc was collected and digested with 0.5% collagenase type II for 2 hours at 37°C. Afterward, the digested NP tissue was incubated using 15% fetal bovine serum (Invitrogen) and 1% penicillin/streptomycin in Dulbecco’s Modified Eagle Medium/ Nutrient Mixture F‐12 (DMEM/F‐12; Invitrogen). Then, the tissue was cultured in an incubator with 5% CO2 at 37°C. Cells were harvested with 0.25% trypsin‐EDTA (Invitrogen) upon confluent and transferred into a 10 cm culture dish at a proper density. The complete medium was changed every other day, and the first three passages of NP cells were used in the experiments.

H₂O₂ was used to induce the degeneration of NP cells.²⁴ All the NP cells were divided into 4 groups: the BSHXF group with BSHXF containing serum intervention, the control group with control serum in, H₂O₂ group with 75 μM H₂O₂ and the control serum treated for 12 hours, and BHSXF+ H₂O₂ group, which was treated with 75 μM H₂O₂ and BHSXF containing serum for 12 hours.

Evaluation of BSHXF for IVDD in Animal Model and NP Cells in Vitro

Magnetic Resonance Imaging Examination

At 4 weeks post-modeling, magnetic resonance (MR) images were obtained to evaluate the effects of needle puncture and BSHXF on IVDD. A Siemens Trio Tim 3.0T MR scanner (Siemens Medical Solutions, Erlangen, Germany) was used with a quadrature extremity coil receiver, and T2-weighted sagittal plane images covering the entire experimental disc area were acquired. Disc degeneration was graded based on signal intensity alterations using a modified Thompson classification.²⁵ Two experienced observers, blinded to the study, independently assessed all images, and the mean of their evaluations was recorded. Subsequently, disc samples were collected for further analysis after drug administration.

Histological Staining

Following a three-day fixation in 4% paraformaldehyde, the tissues underwent decalcification with a 14% EDTA solution (pH = 7.4) for 21 days. Post-decalcification, the samples were dehydrated, paraffin-embedded, and sectioned coronally at a thickness of 5 μm. These sections were then stained with hematoxylin and eosin (H&E).²⁶ Disc degeneration extent was evaluated using a semi-quantitative grading scale considering the cellularity, morphology of the AF and NP, and the border between the two structures, as described by Han et al.²⁷ Two independent, blinded observers conducted the histological assessments, and the data were presented as the mean of their evaluations.

Immunohistochemistry

Paraffin-embedded sections were deparaffinized, rehydrated, and then underwent microwave treatment in 0.01 mol/L sodium citrate for 15 minutes each. Subsequently, 3% hydrogen peroxide was used to block endogenous peroxidase activity for 10 minutes, followed by blocking nonspecific binding sites with 1% bovine serum albumin for 30 minutes at room temperature. The sections were then incubated with the primary antibody overnight at 4 °C. Following washing, the sections were incubated with a secondary antibody mixture at room temperature for 40 minutes. Visualization of type II collagen was achieved through the diaminobenzidine method, with counterstaining using hematoxylin. Two independent pathologists, blinded to the immunohistochemical details, conducted the analysis of type II collagen.

Cell Viability Assay

The CCK-8 method was used to quantitatively analyze proliferation activity of NP cells of different groups. Cells were seeded in 96-well plates and categorized into four groups. Each group was subjected to treatment with BSHXF-containing serum or mouse control serum, with or without H₂O₂, for varying durations (0, 12, 24, 48, and 72 hours) to evaluate cell viability. Subsequently, 10 μL of CCK-8 reagent and 90 μL of culture media were added to each well and incubated for 2 hours at 37°C. The absorbance values at 450 nm (OD450) were measured using an enzyme marker, with four replicate wells for each group.

RNA Isolation and Reverse Transcription-Quantitative (RT-q) PCR

Following the collection of intervertebral disc tissue from different groups of mice and NP cells cultured in vitro, total RNA was extracted using TRIzol (Invitrogen; Thermo Fisher Scientific, Inc.). RNase-free water was employed to dissolve the RNA, and its concentration was determined using a spectrophotometer. SuperScript IV Reverse Transcriptase (Invitrogen; Thermo Fisher Scientific, Inc.) was utilized to reverse transcribe the mRNA to cDNA, following the manufacturer’s protocol. Primers designed by Shanghai Jierui Biotechnology Co. Ltd. were employed to amplify cDNA. The cDNA amplification was carried out using PowerUp™ SYBR Green Master Mix (Invitrogen; Thermo Fisher Scientific, Inc.). The thermocycling conditions involved an initial step at 50°C for 2 minutes, followed by 95°C for 2 minutes, and then 30 cycles of denaturation at 95°C for 3 seconds and annealing at 60°C for 30 seconds, as per the manufacturer’s protocol. GAPDH expression served as the control, and the 2-ΔΔCq method was employed to calculate the relative RNA expression. The primer sequences for the PCR genes are provided in Table 1. SOD-1 and GPX-1 were assessed in mouse disc tissue, while collagen II, aggrecan, MMP13, SOD-1, and GPX1 were evaluated in the various groups of cultured NP cells.

Table 1 Primer Sequences of the Genes for PCR

Metabolomics Analysis

Sample Preparation

Samples of BSHXF group, model group, sham group and QC group were analyzed based on UHPLC-HRMS/MS (Ultra-High Performance Liquid Chromatography- High Resolution Tandem Mass Spectrometry) -based non-targeted metabolomics platform. To 30 mg tissue, 1.2 mL of 80% methanol was added and vortexed and homogenized with ice bath by using a TissueLyser (JX-24, Jingxin, Shanghai) with zirconia beads for 4 minutes at 40 Hz. Then allowed to stand at −40 °C for 30 minutes and 4 °C for 10 min, the mixture centrifuged at 15,000 g and 4°C for 15min. A 500 μL supernatant was evaporated to dryness, and reconstituted in 50 μL of 50% methanol (including internal standard) prior to perform UHPLC-HRMS/MS analysis. Quality control (QC) sample was obtained by isometrically pooling all the prepared samples.

UHPLC-HRMS/MS Analysis

Chromatographic separation was carried out using a ThermoFisher Ultimate 3000 UHPLC system with a Waters ACQUITY UPLC BEH Amide column (2.1mm × 100 mm, 1.7μm). The mobile phases comprised (A) water with 15mM ammonium acetate (pH=9) and (B) 90% acetonitrile with 10mM ammonium acetate. A linear gradient elution was used with the following program: 0–0.5 min, 95% B; 7 min, 80% B; 8 min, 70% B; 9.5–11 min, 20% B; 11.5–14 min, 20% B. The flow rate was set at 0.25 mL/min. The eluents were then analyzed using ThermoFisher Q Exactive™ Hybrid Quadrupole-Orbitrap™ Mass Spectrometry (QE) in both Heated Electrospray Ionization Positive (HESI+) and Negative (HESI-) modes. Spray voltage was set to 3500 V. Capillary and Probe Heater Temperature were separately 350 °C. Sheath gas flow rate was 40 (Arb, arbitrary unit), and Aux gas flow rate was 10 (Arb) for positive mode. S-Lens RF Level was 50 (Arb). The full scan was operated at a high-resolution of 70,000 FWHM (m/z=200) at a range of 60–900 m/z with AGC Target setting at 3×106. Simultaneously, the fragment ions information of top 10 precursors each scan was acquired by Data-dependent acquisition (DDA) with HCD energy at 15, 30 and 45 eV, mass resolution of 17,500 FWHM, and AGC Target of 2×105.

Data Preprocessing of Raw Mass Spectral Data

The typical base peak ion (BPI) chromatograms were achieved by UHPLC-HRMS/MS. The raw data of UHPLC-HRMS/MS were firstly transformed to mzXML format by ProteoWizard and then processed by XCMS and CAMERA packages in R software platform. In XCMS package, the peak picking (method=centWave, ppm=5, peakwidth=c (5, 20), snthresh=10), alignment (bw=6 and 3 for the first and second grouping, respectively), and retention time correction (method=obiwarp) were conducted. In CAMERA package, the annotations of isotope peak, adducts, and fragments were performed with default parameters. The final data was exported as a peak table file, including observations, variables, and peak areas. The peak areas data was normalized to internal standards before performing analysis of univariate or multivariate statistics.

Identification of Differential Metabolites

For multivariate statistical analysis, the normalized data were preprocessed by par scaling and mean centering before performing PCA, PLS-DA, and OPLS-DA. The model quality is described by the R2X or R2Y and Q2 values. R2X (PCA)or R2Y (PLS-DA and OPLS-DA) is defined as the proportion of variance in the data explained by the models and indicates the goodness of fit. Q2 is defined as the proportion of variance in the data predictable by current model and indicates the predictability of current model, calculated by cross-validation procedure. In order to avoid model over-fitting, a default 10-round cross-validation was performed throughout to determine the optimal number of principal components.

For univariate statistical analysis, the normalized data were analyzed in R platform, where parametric test was performed on the data of normal distribution by Welch’s t-test, while nonparametric test was performed on the data of abnormal distribution by Wilcoxon Mann–Whitney test. The p values of univariate statistical analysis lower than 0.05 and fold change (FC) larger than 1.2 were identified as potential differential metabolites. Fold change (Log2FC) was calculated as binary logarithm of average normalized peak intensity ratio between groups. Metabolic pathway analysis was performed MetaboAnalyst5.0.

Results

Network Pharmacological Analysis of BSHXF for IVDD

Active Components in BSHXF and Target Prediction

A total of 135 active ingredients from BSHXF were identified by screening the TCMSP database using specific criteria (OB ≥ 30% and DL ≥ 0.18). These ingredients included 65 from DS, 4 from MG, 20 from NX, 28 from DZ, 7 from WLX, and 11 from TSZ. Subsequently, the DrugBank database was utilized to acquire target genes for each active ingredient, resulting in a total of 240 unique genes after removing duplicates. The gene names were standardized using the UniProt database.

Potential Targets Prediction of IVDD

Disease targets related to IVDD were retrieved from GeneCards, PharmGKB, and OMIM databases. 856 targets of IVDD were obtained from the GeneCards database, 54 targets of IVDD were obtained from the OMIM database, 1927 targets of IVDD were retrieved from the PharmGKB database. After eliminating the repeated targets, a total of 2575 targets were stored, which were closely associated with the onset and progression of IVDD (Figure 2a).

Figure 2 Network pharmacological analysis of BSHXF in treating IVDD. (a) The number of the IVDD disease targets obtained from GeneCards, PharmGKB and OMIM database. (b) Venn diagram summarizing the intersection targets of BSHXF and IVDD. (c) BSHXF-target-IVDD network. The network was constructed by linking 101 candidate compounds from the 6 herbs, which are constituents of BSHXF to the 139 target genes. (d) Protein-protein interaction (PPI) network of BSHXF and IVDD targets. (e) The process of topological screening for the PPI network. Screen the core proteins in two steps based on the thresholds of Betweenness, Closeness, Degree, Eigenvector, LAC, and Network. (f) Gene Ontology (GO) functional enrichment analysis. The top 10 of GO enrichment analysis (BP represented biological progress of core targets. CC represented cellular components of core targets. MF represented the molecular function of core targets). (g) Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis. The top 30 signaling pathways from KEGG analysis.

Ingredient-Target-Disease Network Construction

To identify the intersection of BSHXF target genes and IVDD target genes, the Venn diagram analysis was carried out. As shown in Figure 2b, a total of 139 target genes were identified. Subsequently, an IVDD target network was constructed by using Cytoscape to identify the relationships between 6 herbs and their corresponding compounds of BSHXF and IVDD common targets. As shown in Figure 2c, the network was constructed by linking 101 candidate compounds from the 6 herbs, which are constituents of BSHXF to the 139 target genes.

PPI Network Construction and Core Genes Identification

To understand the relationship among common genes and obtain core genes, 139 common genes were first imported into the STRING database, the human genome was chosen, and the PPI network was constructed (Figure 2d). Second, we used the CytoNCA to calculate topological parameters, and screen the core proteins in two steps based on the thresholds of Betweenness, Closeness, Degree, Eigenvector, LAC and Network. Consequently, a total of 15 core proteins were identified, including AKT1, MAPK1, MAPK14, JUN, TP53, CTNNB1, CCND1, MYC, HIF1A, CDKN1A, FOS, ESR1, STAT3, RB1, RELA, which were closely related to IVDD (Figure 2e).

Functional Enrichment Analysis

To explore the potential signaling pathways and biological processes regulated by BSHXF, GO and KEGG functional enrichment analyses were performed based on its potential targets. Figure 2f showed the top 10 entries of biological progress (BP), cellular components (CC) and molecular functions (MF). BP analysis showed that potential targets were primarily focused on the response to radiation, response to oxygen levels, response to lipopolysaccharide, response to oxidative stress, etc., while CC analysis showed that potential targets were primarily focused on membrane raft, membrane microdomain, transcription regular complex, protein kinase complex, etc. In addition, MF analysis showed that potential targets were primarily focused on DNA-binding transcription factor binding, RNA polymerase II-specific DNA-binding transcription factor binding, integrin binding, G protein-coupled amine receptor activity, etc.

To analyze the representative pathways related to the core targets, the KEGG enrichment analysis was carried out. The 30 most significant KEGG signaling pathways were listed, in which PI3K-Akt signaling pathway, AGE-RAGE signaling pathway, cell senescence, IL-17 signaling pathway, TNF signaling pathway and apoptosis are the most related pathway to IVDD (Figure 2g).

Experimental Validation of BSHXF Alleviates the Progress of IVDD in Vivo

MRI evaluation of the disc degeneration obtained 4 weeks after puncture had stronger T2-weighted signal intensities in the BSHXF group than in the model group. Representative MR images revealed restoration of the signal intensities and areas of the discs in the BSHXF group compared with the placebo-treated mice (Figure 3a). The sham operation group remains the highest signal intensity and intact structure. The Thompson MRI grade score, which indicates the degree of disc degeneration, was significantly lower in the BSHXF group than the model group (BSHXF 1.375 ± 0.3536 vs Model 2.250 ± 0.2673; P < 0.05) (Figure 3b).

Figure 3 MRI, HE staining and immunohistochemical and PCR evaluation for disc degeneration and oxidative stress in vivo. (a) Representative T2-weighted magnetic resonance images of each group 4 weeks after modeling. (b) Quantitative data of MRI grade based on modified Thompson classification in 3 groups at 4 weeks. (c) Representative section HE staining pictures of each group (×40). (d) Quantitative data of histological grades in all groups. (e) Representative type II collagen immunohistochemical images of each group (×40). (f) Quantitative data of type II collagen integrated optical density scores in all groups. (g) Expression of SOD-1 mRNA in all groups. (h) Expression of GPX1 mRNA in all groups.

Histological evaluation of the disc degeneration was also managed and representative pictures were showed in Figure 3c. The histological grade of discs in BSHXF group was significantly lower than the model group (BSHXF 8.250 ± 1.022 vs Model 11.38 ± 1.544; P < 0.05). The sham operation group remained the best histological structure and obtained lowest histological grade (Figure 3d).

Immunohistochemical evaluation was also managed by testing the expression and distribution of type II collagen. Representative images for type II collagen staining are shown in Figure 3e. Type II collagen staining in the BSHXF group was stronger than the model group (BSHXF 3818± 407.6 vs Model 2954 ± 367.3; P < 0.05). The sham operation group obtained most strong type II collagen staining (Figure 3f). From the results above, BSHXF was confirmed to alleviate the process of disc degeneration in vivo.

According to the results of GO analysis of network pharmacology, biology process analysis showed that response to oxygen levels, oxidative stress, hypoxia, decreased oxygen levels and reactive oxygen species are all potential targets, so superoxide dismutase-1 (SOD-1) and glutathione peroxidase 1(GPX1) mRNA were examined by PCR for the investigation of oxidative stress condition of each group. In the BSHXF group, SOD-1was significantly higher than that of the model group (BHSXF 0.6533± 0.0503 vs Model 0.3167± 0.0764; P < 0.05). Similarly, GPX1 was higher in the BSHXF group (BHSXF 0.7267± 0.1419 vs Model 0.2502± 0.0503; P < 0.05). Therefore, BSHXF showed anti-oxidative stress effect in process of disc degeneration in vivo.

Role of BSHXF Containing Serum in NP Cell Degeneration and Oxidative Stress in Vitro

Proliferation ability of NP cells of different groups was detected by CCK-8 method (Figure 4a). No obvious difference showed up between the control group and the BHSXF group in all culture time (0, 12, 24, 48, and 72 hours). After 24 hours, cells cultured with H₂O₂ obtained lower proliferation ability than those in the control group. Besides, when cultured with H₂O₂, BSHXF containing serum was able to increase the NP cells proliferation significantly at 24, 48 and 72 hours (P < 0.05).

Figure 4 Cellular viability and the mRNA expression of type II collagen, aggrecan, MMP13, SOD-1, GPX1 in cultured nucleus pulposus cells in vitro. (a) Cell viability of nucleus pulposus cells when treated with BSHXF containing serum and/or H2O2 after 12, 24, 48 and 72 hours. (b) Expression of type II collagen mRNA in 4 groups. (c) Expression of aggrecan mRNA in 4 groups. (d) Expression of MMP13 mRNA in 4 groups. (e) Expression of SOD-1 mRNA in 4 groups. (f) Expression of GPX1 mRNA in 4 groups.

In addition, H₂O₂ also decrease the expression of collagen II and aggrecan mRNA and increase the expression of MMP13. In NP cells treated with BSHXF containing serum and H₂O₂, the expression level of type II collagen and aggrecan mRNA was higher than the H₂O₂ group (P < 0.05) (Figure 4b and c). The expression of MMP13 showed no significant difference between the H₂O₂+BSHXF group and the H₂O₂ group (Figure 4d). MMP13, aggrecan, type II collagen are the markers for IVD degeneration in terms of extracellular matrix, so BSHXF containing serum could improve the expression of NP extracellular matrix in vitro.

SOD-1 and GPX1 were also examined in the cultured NP cells. No obvious difference showed up between the control group and the BHSXF group. However, the expression of SOD-1 and GPX1 were significantly higher in the H₂O₂+BSHXF group than the H₂O₂ group (p<0.05) (Figure 4e and f). From the results above, BSHXF containing serum can alleviate the degeneration of NP cells after treated with H₂O₂ in vitro, and oxidative stress may play a role in this process.

Metabolomics Analysis

Global Metabolic Analysis

Metabolic analysis of mouse intervertebral discs from the sham group, model group, and BSHXF group was conducted using UHPLC-HRMS/MS (Supplementary Materials). Base peak ion (BPI) chromatograms for the three groups were presented (Figure 5a and b). Global features of the raw data were used for PCA and PLS-DA analyses to visualize metabolic differences and differentiate the groups. PCA analysis showed a tendency of group separation with R2X=0.658 in cation mode and R2X=0.556 in anion mode (Figure 5c and d). PLS-DA analysis demonstrated R2Y=0.658, Q2=0.613 in cation mode, and R2Y=0.955, Q2=0.707 in anion mode, indicating the reliability of the results (Figure 5e and f).

Figure 5 (a) Cation current diagram. (b) Anion current diagram. (c) Cation mode principal component analysis (PCA). (d) Anion mode principal component analysis (PCA). (e) Cation mode partial least squares discriminant analysis (PLS-DA). (f) Anion mode partial least squares discriminant analysis (PLS-DA).

To further distinguish differences between each pair of groups and understand the impact of the modeling operation and BSHXF on metabolite expression, PCA, PLS-DA, and OPLS-DA analyses were carried out for both anion and cation modes (Figure 6). The results showed significant difference of the metabolite expression between the sham group and BSHXF group (Figure 6a–f), between sham group and model group (Figure 6g–l), and between model group and BSHXF group (Figure 6m–r).

Figure 6 PCA, PLS-DA, OPLS-DA diagram of each group in both cation and anion mode. (a) PCA diagram of sham group/model group comparison in cation mode. (b) PCA diagram of sham group/ model group comparison in anion mode. (c) PLS-DA diagram of sham group/ model group comparison in cation mode. (d) PLS-DA diagram of sham group/ model group comparison in anion mode. (e) OPLS-DA diagram of sham group/ model group comparison in cation mode. (f) OPLS-DA diagram of sham group/ model group comparison in anion mode. (g) PCA diagram of sham group/BSHXF group comparison in cation mode. (h) PCA diagram of sham group/ BSHXF group comparison in anion mode. (i) PLS-DA diagram of sham group/ BSHXF l group comparison in cation mode. (j) PLS-DA diagram of sham group/ BSHXF group comparison in anion mode. (k) OPLS-DA diagram of sham group/ BSHXF group comparison in cation mode. (l) OPLS-DA diagram of sham group/ BSHXF group comparison in anion mode. (m) PCA diagram of model group/BSHXF group comparison in cation mode. (n) PCA diagram of model group/BSHXF group comparison in anion mode. (o) PLS-DA diagram of model group/BSHXF group comparison in cation mode. (p) PLS-DA diagram of model group/BSHXF group comparison in anion mode. (q) OPLS-DA diagram of model group/BSHXF group comparison in cation mode. (r) OPLS-DA diagram of model group/BSHXF group comparison in anion mode.

Differential Metabolites and Correlation Analysis

Differential metabolites were selected using volcano plots. A total of 170 metabolites were identified, with 74 upregulated and 96 downregulated when comparing the model group with the sham group (Figure 7a and b). Similarly, 65 upregulated and 77 downregulated metabolites were identified between the sham group and BSHXF group (Figure 7c and d). 8 upregulated land 50 downregulated metabolites were identified between the model group and BSHXF group (Figure 7e and f).

Figure 7 Differential metabolites selection by volcano plot. (a) Volcano plot of sham group/model group in cation mode. (b) Volcano plot of sham group/model group in anion mode. (c) Volcano plot of sham group/BSHXF group in cation mode. (d) Volcano plot of sham group/BSHXF group in anion mode. (e) Volcano plot of model group/BSHXF group in cation mode. (f) Volcano plot of model group/BSHXF group in anion mode.

Pearson correlation and matrix diagrams were utilized to investigate the correlation among the differential metabolites (Figure 8a–c). Heatmaps were constructed to display the relationship between the metabolites of each group and their aggregation relationship (Figure 9).

Figure 8 Correlation matrix plot for the differential metabolites from different groups. (a) Correlation matrix plot for the differential metabolites from the sham group and the model group. (b) Correlation matrix plot for the differential metabolites from the sham group and the BSHXF group. (c) Correlation matrix plot for the differential metabolites from the model group and the BSHXF group.

Figure 9 Heat map for the differential metabolites from different groups. (a) Heat map for the differential metabolites from the sham group and the model group. (b) Heat map for the differential metabolites from the sham group and the BSHXF group. (c) Heat map for the differential metabolites from the model group and the BSHXF group.

Metabolic Pathway Analysis

Pathway enrichment analysis using MetaboAnalyst 5.0 resulted in a topological map of the metabolic pathway enrichment network for each comparison group (Figure 10). Metabolic pathways with significant differences were identified and showed in Table 2. 15 main metabolic pathways were concluded between the sham group and the model group, including Vitamin B6 metabolism, citrate cycle, pyrimidine metabolism, arginine biosynthesis, arginine and proline metabolism, etc. Besides, 13 metabolic pathways were concluded between the model group and the BSHXF group, including glutathione metabolism, alanine, aspartate and glutamate metabolism, glycine, serine and threonine metabolism, D-glutamine and D-glutamate metabolism, etc. Among them, 8 pathways were found to be common to both comparisons, including glutathione metabolism, alanine, aspartate and glutamate metabolism, glycine, serine and threonine metabolism, D-glutamine and D-glutamate metabolism, citrate cycle, arginine and proline metabolism, cysteine and methionine metabolism, tryptophan metabolism. A potential metabolic pathway network was constructed, highlighting the interconnectedness and the complex network of metabolic pathways affected by BSHXF in treating IVDD (Figure 11). It can be seen that BSHXF not only affected a single metabolic pathway but also formed a complex network of metabolic pathways through some interrelated metabolic pathways to achieve the purpose of treating IVDD.

Table 2 Differential Metabolic Pathway Table of the Three Groups of Samples

Figure 10 Enrichment analysis of metabolic pathways in three groups of samples. (a) Disturbed metabolic pathways in the sham group vs the model group. (b) disturbed metabolic pathways in the sham group vs the BSHXF group. (c) Disturbed metabolic pathways in the BSHXF group vs the model group.

Figure 11 Network diagraph of the remarkably perturbed metabolic pathways in the process of IVDD and BSHXF intervention.

Discussion

IVDD is the pathological basis for a widespread degenerative spinal disease and its incidence has been on the rise due to changes in human work habits and lifestyles, affecting individuals at a younger age. This trend is burdening healthcare systems globally.²⁸ Currently, the main interventions for IVDD-related diseases are pain relievers and surgical treatments.²⁹ The medications used are limited to nonsteroidal anti-inflammatory drugs or muscle relaxants to alleviate symptoms.^30,31 However, effective nonsurgical treatments for IVDD are lacking, making the exploration of novel and effective nonsurgical treatments a focus of research in degenerative spinal disorders.

TCM, an important component of complementary and alternative medicine, has evolved over thousands of years with its own unique system of theories, diagnostics and therapies in Asian countries.³² TCM therapies have gained increasing use worldwide, including in Western countries, in the last few decades and are recognized for their significant role in preventing and treating various diseases, including IVDD.³³ According to TCM theory, diseases related to IVDD are known as “Bi-Zheng” and are usually caused by blood stasis and qi stagnation or deficiencies in liver and kidney function.³⁴ The herbs in BSHXF play a role in promoting blood circulation and relieving pain, nourishing the liver and kidney, and strengthening muscles and bones. These properties have been proven to be therapeutically effective in treating IVDD-related diseases. However, the molecular mechanism of BSHXF on IVDD remains largely elusive, significantly hindering its modernization and international acceptance.³⁵ Based on the concept of “Disease-Gene-Target-Medicine”, network pharmacology is an effective approach to elucidate the underlying mechanism of multi-compound Chinese medicine, aligning with the holistic view of TCM.²² Some studies have demonstrated that BSHXF could promote the proliferation of NP cells and remodel the ECM during IVDD.³⁶ However, due to the complex composition and comprehensive action of BSHXF, the specific mechanism of its treatment for IVDD is still unclear. Therefore, we applied network pharmacology combined with experimental validation and metabolomics to explore the potential mechanism of BSHXF for IVDD treatment.

We used Cytoscape software to reconstruct the PPI network of common genes obtained from the STRING database and got 15 core genes (AKT1, MAPK1, MAPK14, JUN, TP53, CTNNB1, CCND1, MYC, HIF1A, CDKN1A, FOS, ESR1, STAT3, RB1, RELA) of BSHXF against IVDD. These core target proteins are involved in multiple signaling pathways such as oxidative inflammation, apoptosis, cell senescence, metabolism, and proliferation-related signaling pathways. GO and KEGG enrichment analysis results showed that the therapeutic targets of BSHXF for diseases mainly focused on the PI3K-Akt signaling pathway, AGE-RAGE signaling pathway, cell senescence, IL-17 signaling pathway, TNF signaling pathway and apoptosis. Key biological processes induced are in response to oxygen levels, oxidative stress, hypoxia, and reactive oxygen species. Previous research has emphasized the critical roles these factors play in IVDD progression. For instance, activation of the PI3K-Akt pathway can lead to a cascade of events, including ECM degradation reduction, apoptosis inhibition, and induction or suppression of autophagy, ultimately protecting against IVDD.^37,38 AGE accumulation is linked to endochondral ossification and oxidative stress, potentially inducing alterations in the oxidative microenvironment of NP tissue and osteogenic differentiation of intervertebral disc cells through the AGE/RAGE axis.³⁹ The TNF and IL-17 pathways synergistically drive IVDD progression, predominantly by promoting inflammatory factor release, NP cell apoptosis, and extracellular matrix degradation.^40–42 Moreover, the MAPK pathway is believed to primarily mediate the inflammatory response of NP cells to TNF-α stimulation.⁴³ Overall, we hypothesize that BSHXF’s therapeutic effect in IVDD may be tied to its anti-oxidative stress effects, inhibiting apoptosis and senescence of NP cells.^44,45

Therefore, experimental validation was managed to investigate the effects and the underlying pharmacological mechanism in this study. Through needle puncture induced mice coccygeal disc degeneration model and H₂O₂ induced degenerative NP cells in vitro, BSHXF showed a significant effect in slow the process of NP cells or disc degeneration. SOD-1 and GPX1 are antioxidant enzymes crucial for scavenging reactive oxygen species. We observed a significant decrease in SOD-1 and GPX1 levels after disc degeneration. However, BSHXF and its containing serum demonstrated the ability to inhibit peroxidation and enhance the expression of antioxidant enzymes, thus alleviating IVDD.

Metabolomics is an integral part of systems biology, providing real-time insights into an organism’s metabolite levels and biometabolic network characteristics. With the rapid advancement of bioinformatics, differential metabolite enrichment analysis has become a widely utilized method to identify potential metabolic processes and signaling pathways associated with differential production. Upon analyzing metabolite differences, we found that the differential marker metabolites between the sham operation group and the model group primarily influenced 11 amino acid metabolism pathways. Additionally, key carbohydrate metabolism pathways like the citrate cycle, nucleotide metabolism pathway (pyrimidine metabolism), as well as cofactors and vitamins metabolism (including vitamin B6 metabolism and nicotinate and nicotinamide metabolism) were implicated. The main biological processes associated with IVDD were identified as energy metabolism, biological oxidation, and structural protein synthesis and decomposition, aligning with the metabolic pathway findings.^46–48

Furthermore, when comparing the model group with the BSHXF group, the differential marker metabolites were prominently linked to 8 amino acid metabolism pathways, encompassing glutathione metabolism, alanine, aspartate and glutamate metabolism, glycine, serine and threonine metabolism, D-glutamine and D-glutamate metabolism, arginine and proline metabolism, cysteine and methionine metabolism, tryptophan metabolism, and lysine degradation. Additionally, 5 carbohydrate metabolism pathways were involved, including the citrate cycle, pentose phosphate pathway, pentose and glucuronate interconversions, glyoxylate and dicarboxylate metabolism, and amino sugar and nucleotide sugar metabolism. Remarkably, over half of these metabolic pathways were closely associated with oxidative stress.

Glutathione is a vital antioxidant and oxygen free radical scavenger in the body, relies on cysteine for its synthesis regulation. Methionine is crucial for neutralizing oxidative free radicals that can damage membrane lipids, thus protecting cellular and mitochondrial structures.⁴⁹ During the degeneration process, a decrease in cysteine and methionine metabolism results in a decline in the disc’s antioxidant function, leading to peroxidation of tissue membranes and promoting disc degeneration.⁵⁰ Glutamine is essential for producing glutathione and helps maintain redox homeostasis. Depletion of glutamine can make NP cells more susceptible to oxidative stress-induced injuries. Glutamine metabolism is closely linked to glycolysis and serine metabolism, where serine can be converted to glycine to achieve redox homeostasis. Glycine, in turn, synthesizes glutathione, enhancing resistance to oxidative stress. Metabolites in tryptophan metabolism act as antioxidants, removing reactive oxygen species and improving resistance to free radical damage.

The citrate cycle serves as an important metabolic hub for major nutrients.⁵¹ Downregulation of the citrate cycle suggests metabolite depletion in the discs. This depletion is attributed to apoptosis of NP cells and a decline in their secretory function. Amino sugars and nucleotide sugars play essential roles in carbohydrate conversion and energy utilization.⁵² The pentose phosphate pathway is involved in supplying energy and reducing equivalents crucial for antioxidant activity. These metabolic shifts during degeneration indicate a transition from a multimodal coexistence of energy utilization to a single mode of proteoglycan structural metabolism, signaling material disintegration in NP tissue during the disc degeneration process.

Above all the metabolomics findings, we propose that the effect of BSHXF in reducing oxidative stress levels is closely related to that of the metabolic pathways of glutathione, methionine, glutamine, cysteine, serine, glycine, tryptophan and the citrate cycle, thereby reflecting the comprehensive regulatory role of BSHXF in the treatment of IVDD. The results conform to the multi-pathway, multi-target, and multi-component characteristics typical of TCM in disease treatment.

Conclusions

In summary, our comprehensive analysis integrating network pharmacology, metabolomics, and experimental validation in both in vivo and in vitro settings demonstrates that BSHXF intervention effectively alleviates the progression of IVDD through anti-oxidative stress effect. The key metabolic pathways involving glutathione, alanine, aspartate and glutamate, glycine, serine and threonine, D-glutamine and D-glutamate metabolism, cysteine and methionine, tryptophan, and the citrate cycle play pivotal roles in reducing oxidative stress levels in IVDD. This study advances our knowledge on the function of BSHXF in treatment of IVDD, and its intricate mechanism is worth further exploration.