Unveiling the Therapeutic Potential of Baicalin in Intervertebral Disc Degeneration: Integrative Bulk and Single-Cell Transcriptome Analysis with Experimental Validation of PANoptosis Inhibition

Introduction

Intervertebral disc degeneration (IVDD) is a leading cause of chronic low back pain, spinal dysfunction, and reduced quality of life.¹ With the global aging population, the incidence of IVDD has significantly increased, becoming a major public health challenge and imposing substantial economic and medical burdens on society.^2,3 Therefore, understanding the pathogenesis of IVDD and developing effective treatment strategies are of critical importance.

The nucleus pulposus (NP), the central component of the intervertebral disc, plays a crucial role in maintaining the disc’s structural integrity and biomechanical function.⁴ Nucleus pulposus cells (NPCs) are responsible for synthesizing and secreting extracellular matrix components that ensure hydration, elasticity, and overall disc integrity.^5,6 However, aging and external stressors can lead to NPC loss and functional decline, contributing to the acceleration of IVDD.^7,8

NPC death is a hallmark pathological feature in the progression of IVDD.⁹ PANoptosis, a recently discovered form of programmed cell death (PCD), integrates the mechanisms of pyroptosis, apoptosis, and necroptosis.¹⁰ This process is regulated by the PANoptosome complex, which includes key mediators such as ZBP1, TAK1, AIM2, RIPK1, and NLRP12.¹¹ Although previous studies have shown that apoptosis, pyroptosis, and necroptosis are closely associated with IVDD, the interaction between these three modes of cell death in IVDD has not been extensively explored.¹² However, a recent study demonstrated that enhancing TAK1 expression in NPCs can alleviate mitochondrial dysfunction, reduce reactive oxygen species (ROS) production, and maintain redox balance, ultimately inhibiting PANoptosis.¹³ Therefore, targeting PANoptosis in NPCs may offer a more effective approach compared to conventional treatments that focus on a single form of cell death.

Baicalin(BAI), a natural flavonoid extracted from Scutellaria baicalensis, possesses various pharmacological properties, including anti-inflammatory, antioxidant, anti-cancer, and antibacterial effects.^14–16 Notably, BAI has been shown to modulate multiple forms of PCD, including apoptosis, pyroptosis, and necroptosis.^17–19 For instance, BAI attenuates hippocampal cell apoptosis by inducing autophagy in a rat status epilepticus model,²⁰ inhibits sodium monourate-induced pyroptosis in HK-2 cells,²¹ and reduces necroptosis in carp gills exposed to chlorpyrifos.²² Moreover, recent studies have found that in a murine model of hemophagocytic lymphohistiocytosis (HLH), BAI inhibits PANoptosis by blocking mitochondrial Z-DNA formation and ZBP1-PANoptosome assembly in macrophages.²³ These findings suggest that BAI may serve as a promising candidate for targeting PANoptosis. However, whether BAI can mitigate IVDD by modulating PANoptosis remains unexplored.In this study, we employed bulk and single-cell transcriptomic analyses to identify potential targets of BAI in inhibiting PANoptosis in NPCs. These findings suggest that BAI may be a promising drug candidate for targeting PANoptosis. However, its potential to alleviate IVDD by modulating PANoptosis remains unexplored.

Therefore, this study aims to explore the potential targets and therapeutic mechanisms of baicalin (BAI) in regulating PANoptosis in nucleus pulposus cells (NPCs) using integrated bulk and single-cell transcriptome analyses. The study is further supported by validation through clinical samples as well as in vitro and in vivo experiments (Figure 1).

Figure 1 Schematic overview of the study design.

Materials and Methods

Analysis of Bulk RNA-Seq Data

Acquisition and Organization of Bulk RNA-Seq Data and PANoptosis-Related Genes

High-throughput sequencing datasets of NP tissue from IVDD were retrieved from the Gene Expression Omnibus (GEO) database, including raw count data from GSE167199, GSE245147, and GSE266883 (Supplementary Tables S1–S3). Gene identifiers were converted to gene symbols using the annotation file of the GPL24676 platform (Supplementary Table S4). Low-expression genes with zero counts across all samples were excluded from further analysis.The TPM (transcripts per million) values for each gene were then calculated based on To account for batch effects, data were normalized using the ComBat function in the sva package (version 3.48.0) in R.²⁴ The datasets were subsequently integrated.PANoptosis-related genes(PRGs) were identified from several public databases: pyroptosis genes were sourced from GeneCards (https://www.genecards.org/), Reactome (https://reactome.org/), and AmiGO2 (https://amigo.geneontology.org/amigo); apoptosis genes were extracted from KEGG (https://www.kegg.jp/kegg/) and Reactome; necroptosis-related genes were obtained from KEGG and AmiGO2. These gene sets were merged and duplicates removed to compile a final list of PRGs.

Differentially Expressed Genes Analysis and Identification of PRGs

Differentially expressed genes(DEGs) analysis was conducted using the Limma package (version 3.56.2).²⁵ Genes with a log2 fold change > 0.585 and an adjusted p-value < 0.05 were categorized as upregulated, while those with a log2 fold change < −0.585 and an adjusted p-value < 0.05 were classified as downregulated. Volcano plots and heatmaps of DEGs were generated using the ggplot2 package (version 3.5.1). DEGs were cross-referenced with the PRGs set to identify potential biomarkers associated with PANoptosis.

Enrichment Analyses

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes(KEGG) pathway enrichment analyses were performed using the clusterProfiler package (version 4.8.3).²⁶ The enrichGO function was used for GO analysis, while the enrichKEGG function was applied for KEGG analysis. Visualizations were generated using the ggplot2 package to display the results of both GO and KEGG analyses.

Protein-Protein Interaction (PPI) Network Analysis

PPI networks were constructed using the STRING database (version 12.0, https://cn.string-db.org/). The intersecting genes from the DEGs analysis were input into Cytoscape software (version 3.7.2), where hub genes were identified based on their degree centrality and network importance.

Single-Gene Gene Set Enrichment Analysis (GSEA)

The correlation between each target gene and all other genes in the expression matrix was calculated using the Spearman rank correlation algorithm. Genes most strongly correlated with the target gene were identified. Subsequent KEGG pathway enrichment analysis was performed on these correlated gene groups using the gseKEGG function in the clusterProfiler package.

Immune Infiltration Analysis

Immune cell infiltration analysis was carried out using the CIBERSORT package. The correlation between immune cell types was determined using the cor function in the corrplot package (version 0.94)²⁷ The immune infiltration score was computed using the gsva function from the GSVA package (version 1.48.3). Additionally, correlation analysis between immune cells and genes was performed using the calculate_sig_score function in the IOBR package (version 0.99.0).

Single Cell RNA-Seq Data Analysis

Acquisition and Quality Control of Single-Cell Datasets

The GSE244889 dataset was retrieved from the GEO database. Single-cell RNA-seq data analysis was performed using the Seurat package (version 5.1.0).²⁸ First, low-expression genes were excluded by requiring expression in at least three cells. Additionally, low-quality cells were removed based on the criteria that each cell must express a minimum of 300 genes. The PercentageFeatureSet function was used to calculate the ratio of mitochondrial genes to erythrocyte genes. The screening criteria for subsequent analysis were: mitochondrial gene ratio < 20%, erythrocyte gene ratio < 3%, gene count between 300 and 700, and cells with gene counts ranking in the top 99%.

Cell Annotation

Data normalization was performed using the LogNormalize method in the NormalizeData function. The top 3000 highly variable genes (HVGs) were identified using the vst method from the FindVariableFeatures function. Principal component analysis (PCA) was then applied to the normalized data using the RunPCA function. Data integration was carried out using the canonical correlation analysis (CCA) method through the IntegrateLayers function. Dimensionality reduction was performed using the first 18 principal components (PCs), and a shared nearest neighbor (SNN) graph was constructed using the FindNeighbors function. Clustering was conducted with the FindClusters function at a resolution of 0.5. UMAP visualization was performed using the RunUMAP function, and marker genes for each cluster were identified with the FindAllMarkers function (min.pct = 0.25, logfc.threshold = 0.5) for cell type annotation.

Hub Genes Positioning and DEGs Analysis

First, the FeaturePlot function was employed to visualize the distribution of hub genes across different cell types.

Subsequently, the FindMarkers function was applied to identify differentially expressed genes (DEGs) between the control and IVDD groups for each cell type. Upregulated genes were defined as those with a log2 fold change > 0.585 and an adjusted p-value < 0.05, while downregulated genes had a log2 fold change < −0.585 and an adjusted p-value < 0.05. To further illustrate the DEGs, volcano plots were generated using the ggplot2 package, providing a clear visualization of the gene expression changes.

Enrichment Analysis and Pathway Scoring

GO and KEGG pathway enrichment analyses were performed using the enrichGO and enrichKEGG functions, respectively. Pathway scoring was carried out using the AUCell package (version 1.22.0), with gene sets downloaded from the MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/) for KEGG pathways, from Reactome for pyroptosis pathways, and from KEGG for necroptosis pathways. These gene sets were integrated for pathway scoring using AUCell. Additionally, pathway enrichment analysis was conducted using the gsva function in the GSVA package, applying the Gaussian algorithm to calculate pathway activity.

Pseudotime Analysis

Pseudotime analysis was performed using the Monocle2 package (version 2.28.0)²⁹ to model differentiation trajectories. DEGs between differentiation states were identified using the differentialGeneTest function. The BEAM function was employed to detect genes that were dependent on different differentiation directions, allowing for the construction of dynamic differentiation trajectories.

Mechanism of BAI in the Treatment of IVDD

Molecular Docking

To investigate the mechanism by which BAI exerts therapeutic effects on IVDD, molecular docking was performed to explore the interaction between BAI and five key genes identified through bulk and single-cell transcriptomic analysis. The three-dimensional structures of the target proteins were retrieved from the Protein Data Bank (PDB, https://www.rcsb.org) and pre-processed using PyMOL (version 4.60), where water molecules and co-crystal ligands were removed. Polar hydrogen atoms were added using AutoDockTools (version 1.5.7), and the protein structures were saved in PDBQT format for docking. The binding sites of the proteins were predicted using the proteins.plus website (https://proteins.plus/). The structure of BAI was obtained from the Traditional Chinese Medicine Systems Pharmacology Database (TCMSP, https://old.tcmsp-e.com/tcmsp.php), optimized using Chem3D, and converted to PDBQT format. Molecular docking simulations were conducted using AutoDock Vina (version 1.12), with the search space defined based on the predicted binding pockets to focus the docking process on relevant interaction sites. The best conformation, representing the strongest binding affinity, was visualized using PyMOL.

In vitro Experimental Validation

Acquisition of NP Samples

This study was approved by the Ethics Committee of the Affiliated Hospital of Traditional Chinese Medicine, Southwest Medical University, and conducted in accordance with the Declaration of Helsinki. Written informed consent was obtained from all participants. IVDD samples were obtained from patients with lumbar disc herniation confirmed by magnetic resonance imaging (MRI), with Pfirrmann grades IV or V. Non-degenerative control samples were collected from patients with vertebral fractures or adolescent idiopathic scoliosis, confirmed by MRI to have Pfirrmann grade I discs and no history or imaging evidence of disc degeneration. All intervertebral disc tissues were immediately frozen in liquid nitrogen after surgery and stored at −80°C for further analysis.

Experimental Cells and Reagents

Human immortalized NPCs (iCell-0028a) were purchased from iCell. BAI (HY-N0197, purity 98.92%) and IL-1β (HY-P7028) were sourced from MedChemExpress. MEM/F12 medium (A4192001) and fetal bovine serum (FBS, A5670701) were obtained from Gibco. Antibodies against FOS (T56596S), H2AC6 (PK02295S), CASP1 (TA5418S), H1-2 (PHT2378S), BCL2L11 (TN23627S), and GAPDH (TA7021) were purchased from Abmart, while the corresponding primers were sourced from Amresco. The CCK8 assay kit (96992) and BCA protein assay kit (71285-M) were obtained from Sigma. Reverse transcription and qPCR kits (HY-K0510A and HY-K0501) were purchased from MedChemExpress. Hematoxylin-eosin (HE) staining kit (G1120) and Annexin V-FITC/PI apoptosis detection kit (CA1020) were obtained from Solarbio.

Experimental Methods

Preparation of Experimental Solutions

BAI (10 mg) was dissolved in 224.03 μL of DMSO to create a 100 mmol/L stock solution, which was stored at −80°C for future use. During experiments, the stock solution was diluted with MEM/F12 medium containing 10% FBS to achieve the desired working concentration. IL-1β was prepared with PBS to the required working concentration prior to use.

Culture and Treatment of Human NPCs

NPCs were cultured in MEM/F12 medium supplemented with 10% FBS and maintained in a 37°C incubator with 5% CO₂. Cell growth was monitored regularly under a microscope. Upon reaching approximately 80% confluence, cells were passaged and the 3rd to 4th generation cells were used for experimentation. To induce an IVDD in vitro model, cells were exposed to 10 ng/mL IL-1β for 24 hours.^30–32 After IL-1β treatment, cells in the BAI treatment group were subsequently treated with the specified concentration of BAI for the appropriate duration.

CCK8 Assay

To assess the effect of BAI on NPC proliferation, cells were seeded at a density of 5000 cells/well in 96-well plates. The cells were treated with varying concentrations of BAI (0, 0.01, 0.1, 1, 10, and 100 μmol/L) for 24 and 48 hours. For the IVDD model, cells were first exposed to 10 ng/mL IL-1β for 24 hours, followed by treatment with different concentrations of BAI. Six replicates were included for each group. After the treatment period, culture medium was removed, and 100 μL of 10-fold diluted CCK8 solution was added to each well. After 1 hour of incubation, the optical density (OD) at 450 nm was measured using a microplate reader.

Western Blot

The Western blot(WB) process is as follows.For protein extraction from clinical samples, tissues were lysed in RIPA buffer containing protease inhibitors. After ultrasonic disruption, lysates were centrifuged at 12,000 rpm for 10 minutes at 4°C, and the supernatants were collected. Protein concentrations were quantified using a BCA protein assay kit. For the NPCs experiment, cells were divided into three groups: the normal group (culture medium only), the model group (IL-1β-treated), and the treatment group (IL-1β + BAI). Cells in each group were seeded in 6-well plates at 2 × 10₅ cells/well. The model group was treated with 10 ng/mL IL-1β for 24 hours, followed by culture in fresh MEM/F12 medium for another 24 hours. The treatment group was treated with BAI at the optimal concentration for 24 hours following IL-1β treatment. After washing the cells with PBS, proteins were extracted and quantified. Proteins were separated by SDS-PAGE, transferred to PVDF membranes, and incubated with primary and secondary antibodies. Protein bands were visualized, and grayscale values were calculated using ImageJ software.

Quantitative PCR (qPCR) Analysis

RNA was extracted from clinical samples and NPCs using TRIzol reagent following the manufacturer’s instructions. RNA purity and concentration were assessed using a spectrophotometer. For cDNA synthesis, 1 μg of RNA was reverse transcribed. qPCR was performed according to the manufacturer’s protocol, and gene expression levels were calculated using the 2^(-ΔΔCt) method. Relative expression levels between groups were statistically analyzed, and bar graphs were generated. Primer sequences are listed in Table 1.

Table 1 The Sequence of the Primers Used in the Experiment

Flow Cytometry

Flow Cytometry(FC) for detecting cell apoptosisCells were treated as described above. After treatment, culture medium was collected and cells were washed with PBS, digested, and resuspended in the medium. Cells were centrifuged and resuspended in cold PBS, followed by a second centrifugation. Apoptosis was detected using the Annexin V-FITC/PI apoptosis detection kit according to the manufacturer’s instructions, and analysis was conducted using FlowJo software.

X-Ray and Hematoxylin-Eosin (HE) Staining

All experimental rats were obtained from the Laboratory Animal Center of Southwest Medical University. All procedures were performed in accordance with the ARRIVE guidelines, the National Institutes of Health Guide for the Care and Use of Laboratory Animals, and the Regulations for the Administration of Affairs Concerning Experimental Animals of China. The protocol was approved by the Animal Ethics Committee of Southwest Medical University.

Eight-week-old male Sprague-Dawley (SD) rats were randomly assigned to three groups (n=8): sham operation, IVDD model, and BAI treatment. IVDD was induced by puncturing the L3-L5 intervertebral disc with an 18-gauge needle, while the sham group underwent the same procedure without disc puncture. The BAI treatment group was administered 40 mg/kg BAI daily for 1 month post-surgery. The sham and model groups received an equal volume of distilled water. After the intervention, rats were anesthetized, and X-ray imaging of the tail was performed to evaluate intervertebral disc height. Following euthanasia, the L3-L5 intervertebral discs were harvested, fixed, decalcified, dehydrated, embedded, sectioned, and stained with HE for histopathological examination.

Data Processing and Statistical Analysis

Experimental data were analyzed using GraphPad Prism (version 9.5). The Shapiro–Wilk test was applied to evaluate data normality. For datasets with a normal distribution, group comparisons were performed using a t-test for two groups or one-way ANOVA followed by a post-hoc test for multiple groups. For non-normally distributed data, the Mann–Whitney U-test was utilized. Quantitative results are presented as mean ± standard deviation (SD), with statistical significance defined as a p-value < 0.05. Figures and illustrations were created using GraphPad Prism and Adobe Illustrator (version 2023).

Results

Bulk Transcriptome Analysis

DEGs

The GSE167199, GSE245147, and GSE266883 datasets were retrieved from the GEO database, including 18 samples-9 normal and 9 IVDD samples. Normal NP tissues were sourced from patients with spinal cord injury, lumbar fractures, or idiopathic scoliosis (Pfirrmann grade I–II), while IVDD samples were obtained from patients with lumbar disc herniation or spondylolisthesis (Pfirrmann grade IV–V) (Supplementary Table S5).Following quality control, ID conversion, and TPM value normalization, the TPM expression matrix for each dataset was constructed (Supplementary Tables S6–S8). The integrated matrix, after batch effect removal, was generated (Supplementary Table S9). DEGs analysis using the limma package identified 1,463 DEGs, with 713 up-regulated and 750 down-regulated genes (Figure 2A, Supplementary Table S10). A heatmap of the top 20 up-regulated and down-regulated genes is shown in Figure 2B. After merging genes associated with pyroptosis, apoptosis, and necroptosis (Supplementary Table S11) and removing duplicates, 494 PRGs were identified (Supplementary Table S12). Intersection analysis of the DEGs and PRGs revealed 45 common genes (Figure 2C, Supplementary Table S13).

Figure 2 (A) Volcano plot of DEGs; (B) Heatmap of the top 20 up-regulated and down-regulated DEGs; (C) Venn diagram showing the intersection of DEGs and pan-apoptosis-related genes; (D) GO enrichment analysis results (top 10 in BP, CC, and MF); (E) KEGG enrichment analysis results (top 20); (F) PPI network constructed using the STRING database; (G) PPI network constructed using Cytoscape; (H) Single-gene GSEA analysis of FOS, CASP1, H1-2, BCL2L11, and H2AC6.

GO and KEGG Enrichment Analysis

GO and KEGG pathway enrichment analyses were performed on the 45 intersection genes. The GO analysis identified 97 biological processes (BP), 34 cellular components (CC), and 23 molecular functions (MF) associated with these genes (Supplementary Table S14). The top 10 significantly enriched pathways for BP, CC, and MF are shown in Figure 2D. KEGG analysis revealed 134 enriched pathways (Supplementary Table S15), with significant involvement in key pathways such as Apoptosis, Neutrophil extracellular trap formation, Systemic lupus erythematosus, Necroptosis, and Alcoholism (Figure 2E).

PPI Network Construction and Hub Gene Identification

The PPI network for the 45 intersection genes was constructed using the STRING database (Figure 2F). Five hub genes—FOS, CASP1, H1-2, BCL2L11, and H2AC6—were identified based on degree centrality using Cytoscape software (Figure 2G, Supplementary Table S16). These genes are likely to play critical roles in the pathogenesis of IVDD.

Single Gene GSEA Analysis

To further explore the biological functions of these hub genes in IVDD, single-gene set enrichment analysis (ssGSEA) was performed. The results showed significant enrichment in pathways related to apoptosis, pyroptosis, and necroptosis (Figure 2H). Additionally, Figure 3A and B highlight the key positions of these genes within the apoptosis and necroptosis signaling pathways. These findings suggest that the five hub genes may contribute to IVDD through the regulation of multiple cell death pathways.

Figure 3 (A) Apoptosis pathway in KEGG; (B) Necroptosis pathway in KEGG; (C) Expression of different immune cells; (D) Expression of different immune cells in control group and IVDD group (Compared to the control group, *P<0.05, **P<0.01); (E) Correlation between different immune cells; (F) Correlation between 5 hub genes and various immune cells.

Immune Infiltration Analysis

Immune infiltration analysis identified several key cell types, including Macrophages_M0, Macrophages_M2, T_cells_CD4_memory_resting, and B_cells_naive (Figure 3C). B_cells_naive was significantly increased in the IVDD group (P<0.05), while Macrophages_M2 showed a significant increase (P<0.01, Figure 3D). The intercellular correlation heatmap revealed a positive correlation between Macrophages_M2 and B_cells_naive (Figure 3E, P<0.05). Furthermore, correlation analysis showed a positive relationship between the five hub genes and various immune cells, including B_cells_naive, Tregs, Macrophages_M0, Macrophages_M2, Dendritic_cells_activated, and Mast_cells_resting (Figure 3F, P<0.05), with the strongest correlation observed with Macrophages_M2.These results suggest that the hub genes may regulate IVDD not only through the PANoptosis network but also by modulating the immune microenvironment, particularly the function of Macrophages_M2.

Single-Cell Transcriptome Analysis

Data Quality Control

The GSE244889 dataset, comprising 4 normal and 3 IVDD samples, was used for single-cell analysis. Normal samples were derived from trauma patients (Pfirrmann grade I–II), and IVDD samples from patients with degenerative spinal diseases (Pfirrmann grade III–IV). After quality control, 51,868 cells were retained for subsequent analysis.

Cell Annotation

Dimensionality reduction and clustering analysis revealed 14 clusters (Supplementary Figure 1). Based on characteristic marker genes (Supplementary Figure 2, Supplementary Table S17), these clusters were annotated as 8 distinct cell types: Chondrocyte, Neutrophil, Plasma cell, T cell, Pro-Myelocyte, Macrophage, Endothelial cell, and Smooth muscle cell (Figure 4A). Marker genes for each cell type are shown in Figure 4B.

Figure 4 (A) Cell annotation and proportion of various cells; (B) Expression of marker genes in various cells; (C) The distribution of hub genes in cells. (D) Volcano plot of DEGs between control group and IVDD group; (E) Volcano plot of DEGs in chondrocytes; (F) GO enrichment analysis results of DEGs in chondrocytes; (G) Expression of various pathways in control group and IVDD group, The red rectangles highlight the key pathways involved in PANoptosis, including apoptosis, pyroptosis, and necroptosis; (H) Expression of apoptosis, pyroptosis, and necroptosis pathways in control group and IVDD group.

Hub Genes Positioning and DEGs

Hub gene localization analysis showed that the five hub genes were predominantly expressed in chondrocytes, suggesting their primary role in modulating IVDD through regulation of chondrocyte function (Figure 4C). DEGs across cell types were identified using the FindAllMarkers function, with results summarized in Figure 4D (Supplementary Table S18). In chondrocytes, DEG analysis revealed 698 genes, of which 507 were upregulated and 191 downregulated in the IVDD group compared to normal controls (Figure 4E, Supplementary Table S19).

Enrichment Analysis and Pathway Scoring

GO and KEGG enrichment analyses were performed on the DEGs. The GO analysis revealed key biological processes involved in the DEGs (Figure 4F). Pathway activity scores, derived using GSVA, showed significant differences between IVDD and normal chondrocytes (Figure 4G). Analysis of pyroptosis, apoptosis, and necroptosis pathways revealed significantly higher activity in IVDD chondrocytes compared to normal chondrocytes (Figure 4H), suggesting that chondrocytes in degenerated discs are more prone to these forms of cell death.

Pseudotime Analysis

The differentiation trajectory of chondrocytes was analyzed using Monocle2, revealing a branching trajectory with three distinct states (Figure 5A). Figure 5B highlights key genes that determine the chondrocyte differentiation states, while Figure 5C identifies the critical genes driving the differentiation between S2 and S3 states. The expression of five hub genes progressively increased during differentiation (Figure 5D), suggesting their involvement in the pathogenesis of IVDD.

Figure 5 (A) Changes in cell state over time during differentiation; (B) Heat map of DEGs in various cell differentiation states; (C) Heat map of key genes determining cell differentiation states; (D) Expression levels of the 5 hub genes over differentiation time; (E) Pathway number and weight of interactions between different cell types; (F) Interaction strength between various cell types in the control group and IVDD group; (G) Interaction network between each cell type and other cells; (H) Ligand-receptor interactions between macrophages, neutrophils, and chondrocytes.

Cell Interaction

Cell-to-cell interaction analysis, performed using the CellChat package, revealed strong interactions between macrophages, neutrophils, and chondrocytes (Figure 5E). The interaction strength between various cell types is shown in Figure 5G, and specific interactions between macrophages and chondrocytes in both control and degeneration groups are highlighted in Figure 5F. These interactions were mediated by molecules such as TNF-α, CD44, and TGF-β1 (macrophage-chondrocyte) and RETN, NAMPT (neutrophil-chondrocyte) (Figure 5H), suggesting that macrophages and neutrophils influence chondrocyte function through these signaling molecules in IVDD.

Molecular Mechanisms of BAI in IVDD Treatment

Molecular Docking

Molecular docking analysis demonstrated that BAI strongly binds to the proteins encoded by the five hub genes (BCL2L11, CASP1, FOS, H1-2, H2AC6), with binding energies below −5 kcal/mol, indicating potential biological activity (Figure 6).

Figure 6 A–E Molecular docking results of FOS, CASP1, H1-2, BCL2L11, and H2AC6 with baicalin.

Expression of Hub Genes in Clinical Samples

Western blot analysis of clinical samples revealed significantly higher protein expression of the hub genes in IVDD samples compared to normal samples (P < 0.05, Figure 7A). Corresponding increases in mRNA expression were also observed (P < 0.05, Figure 7B), suggesting the critical involvement of these genes in IVDD progression.

Figure 7 (A) Protein expression and relative expression of the five hub genes in the control and IVDD groups in clinical samples (Compared to the control group, *P<0.05, **P<0.01); (B) mRNA expression of the five hub genes in the control and IVDD groups in clinical samples (Compared to the control group, **P<0.01); (C) Effects of different concentrations of baicalin on NPC proliferation at various intervention time points (Compared to the control group, **P<0.01); (D) Effects of different concentrations of baicalin on NPC proliferation under IL-1β intervention (Compared to the IL-1β group, **P<0.01); (E) Protein expression and relative expression of the five hub genes in the control, IL-1β, and IL-1β+baicalin groups (Compared to the IL-1β group, *P<0.05, **P<0.01); (F) mRNA expression of the five hub genes in the control, IL-1β, and IL-1β+baicalin groups (Compared to the IL-1β group, *P<0.05, **P<0.01); (G) Cell apoptosis in the control, IL-1β, and IL-1β+baicalin groups; (H) X-ray images and relative disc height of intervertebral discs in the control, model, and model+baicalin groups (Compared to the model group, **P<0.01); (I) HE staining results of intervertebral discs in the control, model, and model+baicalin groups (10x).

Effect of BAI on NPC Proliferation

CCK-8 assays indicated that low concentrations of BAI (0.01 and 0.1 μmol/L) did not significantly affect NPC proliferation. However, 1 μmol/L, 10 μmol/L, and 100 μmol/L concentrations significantly enhanced NPC proliferation, with 10 μmol/L showing the most pronounced effect after 24 hours (P < 0.05). After 48 hours, 1 μmol/L and 10 μmol/L still promoted proliferation, while 100 μmol/L inhibited NPC proliferation. However, the effects of 1 μmol/L and 10 μmol/L interventions at 48 hours were less pronounced compared to those at 24 hours (P < 0.05, Figure 7C).

Further evaluation under IL-1β (10 ng/mL) stimulation showed that IL-1β significantly reduced NPC activity, while treatment with 1 μmol/L, 10 μmol/L, and 100 μmol/L BAI effectively reversed this inhibition, with the 10 μmol/L concentration showing the most prominent effect. However, there was no statistically significant difference between the 0.01 μmol/L and 0.1 μmol/L BAI treatment groups and the IL-1β group (P < 0.05, Figure 7D).

Regulation of PANoptosis-Related Genes in NPC by Baicalin

To investigate the regulatory effects of BAI on core genes associated with PANoptosis, an in vitro IVDD model was established. The results showed that, compared to the control group, the protein and mRNA expression levels of 5 core genes were significantly elevated in the model group (P < 0.05). However, after baicalin intervention, the expression of these genes was significantly reduced (P < 0.05, Figure 7E and F). Additionally, flow cytometry analysis demonstrated that BAI significantly alleviated IL-1β-induced apoptosis in NPCs (Figure 7G). These findings suggest that BAI effectively inhibits PANoptosis in NPCs by modulating the expression of BCL2L11, CASP1, FOS, H1-2, and H2AC6.

X-Ray and HE

X-ray imaging revealed a significant reduction in intervertebral disc height in the IVDD model group compared to the sham operation group. In contrast, the BAI treatment group exhibited a substantial improvement in intervertebral disc height (Figure 7H). Histological examination through hematoxylin and eosin (HE) staining showed that, while the intervertebral disc structure was intact in the sham operation group, the model group displayed disorganized disc architecture, with damage to the annulus fibrosus and extrusion of the NP. Conversely, in the BAI treatment group, the intervertebral disc structure was relatively preserved, with reduced damage to the annulus fibrosus and NP, and a more orderly tissue structure. Notably, the degree of degeneration was significantly lower in the BAI-treated group compared to the model group (Figure 7I). These results further corroborate the beneficial effects of BAI in ameliorating IVDD and restoring disc structural integrity.

Discussion

The aging population has led to a significant increase in the incidence of IVDD, which has become a major global public health concern.³³ IVDD is a leading cause of chronic low back pain, severely affecting patients’ quality of life and imposing a substantial medical and economic burden.³⁴ Therefore, understanding the pathogenesis of IVDD and exploring novel therapeutic strategies are critical for improving patient outcomes.

NPCs death plays a central role in the pathological progression of IVDD, including apoptosis, pyroptosis, and necroptosis.³⁵ Apoptosis of NPCs has been widely recognized as a key feature of IVDD. Fas and Fas ligand (FasL) are critical regulators of the extrinsic apoptotic pathway, initiating cell death by activating caspase-8.^36,37 Park et al examined 23 cases of lumbar disc herniation and found that FasL expression was higher in intervertebral disc cells lacking NP, and the proportion of FasL-positive cells increased with age.³⁸ Similarly, Chen et al demonstrated that inducing NPCs apoptosis accelerated IVDD in rats³⁹ Pyroptosis, a recently discovered form of programmed cell death, is initiated by inflammasomes and mediated by Gasdermin proteins, which cause cell membrane rupture and release inflammatory mediators such as IL-1β and IL-18, thereby promoting inflammation.^40–42 The inflammatory microenvironment is a hallmark of IVDD,⁴³ and pyroptosis, as a pro-inflammatory programmed cell death (PCD) mechanism, plays an essential regulatory role in IVDD. Studies have shown that lactic acid stimulation activates the NLRP3 inflammasome, inducing pyroptosis in NPCs and exacerbating IVDD⁴⁴ Gong et al further demonstrated that inhibiting NLRP3-mediated pyroptosis could slow the progression of IVDD.³¹ Necroptosis, activated by death receptor ligands, involves RIPK1, RIPK3, and MLKL in necrosome formation, ultimately leading to cell death.^45–47 Research has shown that necroptosis also contributes significantly to the development of IVDD. Fan et al observed that RIPK3 and MLKL expression was significantly higher in degenerated intervertebral discs compared to normal discs.⁴⁸ Chen et al found that continuous mechanical stress induced necroptosis in NPCs, exacerbating IVDD, while knockdown of RIPK1 reduced cell death and slowed disc degeneration.⁴⁹

PANoptosis, a recently identified form of programmed cell death, is a unique form of inflammatory cell death that integrates apoptosis, pyroptosis, and necroptosis through a multifaceted complex known as the PANoptosome. This complex regulates all three forms of cell death simultaneously.⁵⁰ PANoptosis has been implicated in various diseases, including cardiac dysfunction and retinal ischemia/reperfusion injury. Bi et al showed that FUNDC1 deficiency exacerbates doxorubicin-induced cardiac dysfunction and PANoptosis, contributing to dilated cardiomyopathy.⁵¹ Furthermore, PANoptosis has been linked to tumorigenesis, with emerging evidence supporting its potential as a therapeutic target in cancer treatment.^52,53

Previous research on IVDD focused primarily on individual forms of cell death, without considering their potential co-regulation. However, recent studies have highlighted the role of PANoptosis in IVDD progression. A study by Zhang et al demonstrated that exposure to tert-butyl hydroperoxide (TBHP) altered the expression of key PANoptosis-related proteins in NPCs, leading to IVDD. TBHP treatment upregulated TAK1 expression to maintain mitochondrial redox balance, thereby inhibiting PANoptosis and delaying IVDD progression in NPCs.¹³ This underscores the potential significance of PANoptosis in IVDD and suggests that targeting this pathway may offer a promising therapeutic approach.

This study integrates multi-omics analysis and experimental validation to elucidate the pivotal role of PANoptosis in IVDD for the first time. By analyzing several publicly available datasets, we identified 45 genes associated with PANoptosis and narrowed this list to five hub genes (FOS, CASP1, H1-2, BCL2L11, and H2AC6), which are integral to the signaling pathways of pyroptosis, apoptosis, and necroptosis. These genes were notably enriched in IVDD samples. Further single-gene Gene Set Enrichment Analysis (GSEA) suggests that these hub genes may contribute to IVDD progression by modulating various cell death and inflammatory pathways. Immune infiltration analysis revealed that these genes are potentially involved in the infiltration of immune cells, such as macrophages, underscoring their role in regulating the inflammatory microenvironment associated with IVDD.

In single-cell transcriptomic data, we observed a marked enhancement of PANoptosis signaling in NPCs from IVDD samples. Pseudo-temporal analysis showed that the expression levels of the hub genes increased progressively with cell development, indicating their potential involvement in the degeneration of NPCs. Additionally, cell interaction analysis highlighted significant interactions between immune cells (eg, macrophages and neutrophils) and chondrocytes, mediated by molecules such as TNF-α, TGF-β1, and CD44. These interactions likely exacerbate inflammatory and immune responses in IVDD, reinforcing the role of the immune microenvironment in IVDD pathology. These findings align with previous studies that have explored immune changes in IVDD and suggest that targeting the interaction between immune cells and NPCs may offer a novel therapeutic approach.

BAI, a bioactive compound derived from traditional Chinese medicine, has demonstrated anti-inflammatory and anti-apoptotic properties in various degenerative diseases^54,55 IL-1β, a central pro-inflammatory cytokine, plays a critical role in the pathogenesis of IVDD, promoting inflammation, matrix degradation, cell senescence, and cell death, including apoptosis, pyroptosis, and necroptosis.^56,57 As such, IL-1β is commonly used in IVDD in vitro models. Existing literature suggests that BAI can mitigate IL-1β-induced damage. For instance, J et al demonstrated that BAI alleviated IL-1β-induced osteoarthritis chondrocyte damage through mitophagy activation,⁵⁸ while Li et al reported that BAI protected OA chondrocytes from IL-1β-induced apoptosis and extracellular matrix degradation.⁵⁹ Furthermore, BAI has been shown to inhibit pyroptosis, apoptosis, and necroptosis in a variety of disease models. In a mouse model of cardiac injury, BAI reduced Ang II–induced cardiomyocyte apoptosis.⁶⁰ Additionally, in an acute skin injury model, BAI attenuated damage by reducing oxidative stress and inhibiting pyroptosis via the NF-κB/NLRP3/GSDMD pathway.⁶¹ Zhao et al also found that BAI inhibited necroptosis in a model of unexplained recurrent abortion (URSA) by preventing mitochondrial fission.Notably, recent studies have demonstrated that BAI can inhibit PANoptosis in macrophages, blocking mitochondrial Z-DNA formation and the ZBP1-PANoptosome complex in a mouse model of hemophagocytic lymphohistiocytosis (HLH).²³ These studies collectively suggest that BAI may be a promising therapeutic agent for PANoptosis-related diseases.

To assess BAI’s effects on PANoptosis in IVDD, we conducted molecular docking and experimental verification. Our results show that BAI exhibits strong binding affinity to the hub genes and significantly downregulates their expression. Moreover, BAI promoted NPCs proliferation and alleviated IL-1β-induced cellular damage. Flow cytometry further demonstrated that BAI inhibits IL-1β-induced apoptosis in NPCs. In the IVDD rat model, BAI significantly ameliorated IVDD. These findings highlight the potential of BAI as a therapeutic agent for IVDD through modulation of the PANoptosis signaling pathway.

Despite these significant findings, this study has several limitations. The small clinical sample size may restrict the generalizability of the results. Additionally, the specific molecular targets through which BAI regulates PANoptosis and its downstream signaling pathways warrant further investigation.

Conclusion

This study integrated bulk and single-cell transcriptomic analyses to systematically elucidate the pivotal role of PANoptosis in IVDD. We identified key biomarkers associated with PANoptosis and highlighted their significant involvement in the pathogenesis of IVDD. Through experimental validation, we demonstrated for the first time that BAI effectively modulates the PANoptosis signaling pathway in NPCs, mitigating cell degeneration and the pathological progression of IVDD. These findings enhance the mechanistic understanding of IVDD and provide a novel theoretical foundation for the development of therapeutic strategies targeting PANoptosis.