Exploring the mechanisms of PANoptosis in osteoarthritis and the therapeutic potential of andrographolide through bioinformatics and single-cell analysis

Abstract

Background

Osteoarthritis (OA) is a degenerative joint disease marked by the breakdown of cartilage, where apoptosis plays a key role. Although apoptosis-related genes in OA have been studied, a detailed analysis of PANoptosis-related genes and the search for therapeutic drugs remains limited.

Methods

We performed a bioinformatics analysis combined with single-cell RNA sequencing to examine PANoptosis-related gene expression in OA cartilage. Key PANoptosis genes and critical cell populations involved in OA progression were identified. Drug prediction led to the selection of Andrographolide (AG), whose effects were validated through molecular docking, Western blotting, and qRT-PCR in chondrocyte models.

Results

Several PANoptosis-related genes, including CASP8, TLR3, CASP1, and IL18, were significantly differentially expressed in OA. These genes are linked to processes such as apoptosis, pyroptosis, and the inflammasome complex. Pathway analysis revealed necroptosis, Toll-like receptor, and apoptosis signaling pathways as important in OA pathology. Single-cell analysis identified HomC, EC, and preHTC as key cell populations. AG was predicted to regulate PANoptosis genes, which was confirmed experimentally, demonstrating AG's potential to modulate key genes involved in cartilage degeneration.

Conclusion

This study highlights PANoptosis-related genes in OA and identifies Andrographolide as a promising therapeutic drug, offering new insights into OA treatment strategies.

Introduction

Osteoarthritis (OA) is a common degenerative joint disease characterized by the gradual destruction of articular cartilage, the formation of osteophytes, and inflammation within the joint cavity [1, 2]. The major risk factors include age, hormone levels, obesity, genetics, and joint injury [3]. This disease not only severely affects the quality of life of patients but also imposes a significant economic burden on society and families [4]. Although the pathogenesis of OA is complex and not yet fully understood [5]. There is growing evidence that apoptosis plays a crucial role in the development and progression of OA [6,7,8].

Apoptosis is a form of programmed cell death that plays a crucial role in maintaining tissue homeostasis and regulating development [9]. In the pathological process of OA, the apoptosis of chondrocytes is considered one of the important mechanisms leading to cartilage degradation and destruction [10]. Previous studies have shown that apoptosis-related genes are abnormally expressed in OA cartilage and are closely related to the survival and function of chondrocytes [11, 12]. However, most existing studies focus on the role of individual apoptosis genes, lacking a comprehensive analysis of the expression and function of PANoptosis-related genes in OA. Recent research has revealed that PANoptosis is an inflammatory programmed cell death regulated by the PANoptosome complex [13]. It encompasses key features of pyroptosis, apoptosis, and necroptosis [14]. PANoptosis cannot be solely characterized by any single pathway of pyroptosis, apoptosis, or necroptosis; it is a process integrating multiple cell death features. PANoptosis has emerged as a significant factor contributing to the development of various diseases [15]. However, there is a lack of research on the specific mechanisms of its role in OA progression.

With the advancement of high-throughput sequencing technologies, bioinformatics analysis has become a vital tool for studying the molecular mechanisms of complex diseases [16]. Single-cell RNA sequencing (scRNA-seq) technology can reveal gene expression heterogeneity at the single-cell level, providing us with new insights into the disease microenvironment and cellular functions [17]. In this study, we systematically analyzed the expression patterns of PANoptosis-related genes in OA cartilage by integrating bioinformatics analysis and single-cell sequencing technology. We also identified key PANoptosis genes and crucial cell populations associated with OA progression. Despite the complexity of OA's pathological mechanisms, finding effective therapeutic drugs remains one of the key focuses of clinical research. Andrographolide (AG) is an active compound extracted from the traditional Chinese herb Andrographis paniculata, and has been demonstrated to have various pharmacological effects, including anti-inflammatory, anti-cancer, and immune-modulatory activities [18]. In this study, we selected AG as a potential therapeutic drug for OA through drug prediction and validated its regulatory effects on key hub genes and therapeutic efficacy in chondrocyte models via molecular docking, WB, and qRT-PCR experiments. This study aims to explore the involvement of PANoptosis in OA and to identify potential therapeutic drugs through an integrated approach combining bioinformatics, single-cell sequencing, and experimental validation.

Materials and methods

Microarray data source and Limma differential expression analysis

The GSE55235 dataset was obtained from the Gene Expression Omnibus (GEO) database, comprising 20 osteoarthritis (OA) samples and 10 control samples, all derived from human tissue. The data in GSE55235 were sequenced on the GPL96 platform. To identify differentially expressed genes (DEGs), we employed the "Limma" method, based on a generalized linear model. Specifically, we conducted the analysis using the limma package in R (version 3.40.6) [19]. Initially, the expression profile dataset was acquired, followed by multiple linear regression analysis using the lmFit function. Subsequently, the eBayes function was utilized to calculate the standard errors adjusted by empirical Bayes, yielding moderated t-statistics, moderated F-statistics, and log-odds of differential expression. Ultimately, we determined the significance of each gene's differential expression.

Weighted gene co-expression network analysis (WGCNA)

We applied the Weighted Gene Co-Expression Network Analysis (WGCNA) approach to investigate the correlations between genes [20]. Initially, the median absolute deviation (MAD) was computed for each gene, and the least variable 50% of genes were excluded based on their MAD values. Subsequently, the samplegenes function filtered the differentially expressed genes (DEGs) expression data, removing unsuitable genes and samples, and facilitating the creation of a scale-free co-expression network. In the third step, an optimal "soft" threshold power (β) was chosen to estimate the adjacency among genes, calculated from their co-expression similarities. Following this, the adjacency matrix was converted into a topological overlap matrix (TOM), which allowed us to determine gene connectivity and differences, thus enhancing the representation of the gene co-expression network. Hierarchical clustering combined with the dynamic tree cut method was then employed to detect gene modules, organizing genes with analogous expression profiles into various clusters through average linkage hierarchical clustering. In the final step, module eigengenes were compared for dissimilarity, a threshold was set for the module dendrogram, and similar modules were merged for comprehensive analysis. Each identified module consists of genes with correlated expression patterns. To discern the intersection between core module genes and DEGs, the VennDiagram tool was utilized, leading to the identification of key genes potentially involved in OA pathology. This process enabled the construction of a gene co-expression network and highlighted gene modules likely playing crucial roles in osteoarthritis.

Acquisition of PANoptosis-related genes

We obtained a total of 70 PANoptosis-associated genes from the GeneCards database (https://www.genecards.org/) and relevant literature. [15, 21,22,23,24,25,26]. All gene names were converted from symbols to Gene IDs using the "bitr" function in the ClusterProfiler (4.7.2) package.

Osteoarthritis disease dataset acquisition

To investigate the target proteins associated with the disease, we searched the GeneCards database (https://www.genecards.org/) using "Osteoarthritis" as the keyword. The search results from this database were then consolidated and deduplicated to obtain a comprehensive gene set related to osteoarthritis [27].

Protein–protein interaction (PPI) network analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) (https://cn.string-db.org/) is an online platform designed to analyze networks of functional protein associations [28,29,30]. We input the intersected set of WGCNA core module genes, DEGs, GeneCards, and PANoptosis-related genes into the STRING database. Genes with an experimental interaction score greater than 0.4 were considered significant. The PPI network was subsequently visualized using Cytoscape software (version 3.9.1) (www.cytoscape.org/) [31]. To evaluate the significance of nodes within the network, CytoHubba was employed to calculate topological features and centrality metrics. The nodes were color-coded to reflect their importance, with a higher score or index correlating with increased significance. Nodes of high importance were highlighted in red, indicating their potential key roles in the network. From this analysis, the top four red-highlighted genes were identified as PANoptosis hub genes in OA.

Functional enrichment analysis

To uncover the biological functions of the identified genes, we conducted GO and KEGG enrichment analyses on the DEGs. The GO analysis encompassed biological processes (BP), cellular components (CC), and molecular functions (MF) [32]. KEGG analysis was employed to identify signaling pathways associated with OA [33]. The R package ClusterProfiler was used for the functional enrichment analysis, and significantly enriched GO terms and KEGG pathways were determined with adjusted p-values (FDR < 0.05).

Single-gene gene set enrichment analysis (GSEA)

To investigate the key signaling pathways related to PANoptosis hub genes, we conducted a gene set enrichment analysis using GSEA software (version 3.0), which was downloaded from the GSEA website (http://software.broadinstitute.org/gsea/ index.jsp). Samples were categorized into high-expression (≥ 50%) and low-expression (< 50%) groups based on the expression levels of significant hub genes. The c2.cp.kegg.v7.4.symbols.gmt subset from the Molecular Signatures Database (https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/btr260, http://www.gsea-msigdb.org/gsea/downloads.jsp) was used for analysis. We set the minimum gene set size to 5 and the maximum to 5000, with 1000 resamplings. A p-value < 0.05 and FDR < 0.25 were considered statistically significant [34].

Identification of immune-infiltrating cells indisease

The immune microenvironment typically comprises immune cells, inflammatory cells, fibroblasts, and various cytokines and chemokines. Analyzing immune cell infiltration is crucial for predicting disease progression and patient response to treatment. We deconvolved gene expression profiles using linear support vector regression via the CIBERSORT method. Using RNA sequencing data, we estimated the immune cell composition in the samples by performing immune cell infiltration analysis with the CIBERSORT (0.1.0) software tool [35, 36].

Single-cell analysis based on scRNAseq

Single-cell data processing

The single-cell sample data GSE169454 was downloaded from the Gene Expression Omnibus (GEO) database [37]. After importing the data in standard format, the PercentageFeatureSet function in the Seurat package of R was used to calculate the mitochondrial content and rRNA content. The correlation between mitochondrial content and nCount (UMI) and nFeature (gene count) was analyzed. Following data quality control, a significant number of debris and non-target cells remained; therefore, cell filtering criteria were established based on the quality control data: UMI counts were set between 100 and 15,000, and the percentage of mitochondrial content was limited to no more than 10%. After cell filtering, the UMAP distribution map of each sample was examined. If cells were found to cluster within the tissue, batch effect removal was necessary. The CCA method was employed to remove batch effects until the data were free from such effects.

Cell subpopulation clustering and annotation

Typically, different cell types possess distinct cellular structures, allowing for the clustering of similar cells to identify their common features. For dimensionality reduction analysis, 2000 highly variable genes were selected, and principal component analysis (PCA) was applied. Subsequently, the t-distributed stochastic neighbor embedding (t-SNE) algorithm, along with the FindNeighbors and FindClusters (resolution = 0.3) functions in the Seurat package, was used for optimal clustering visualization. Marker genes were defined as those that are highly expressed in most cells within a specific cluster, lowly expressed in other clusters, and significantly upregulated in the target cluster compared to others. The top five most highly expressed genes were used as references, with higher ranking indicating greater expression levels. The bimod test was utilized to perform differential expression analysis between the designated cell clusters and all other clusters, identifying specific marker genes for each cell population. Finally, the SingleR function in the Seurat package, in conjunction with relevant literature, was used to annotate the different cell subpopulations to determine their cell types.

Cell trajectory analysis based on cell type using Monocle2

To explore the developmental trajectories and transitional relationships between different cell types, we performed cell trajectory analysis using Monocle2 [38]. First, the single-cell RNA sequencing data processed and annotated through Seurat were extracted and converted into the format required by Monocle2, ensuring that the data included high-quality cells and representative highly variable genes. Next, the expression data were converted into a CellDataSet object using the newCellDataSet function in the Monocle2 package, with input parameters including the expression matrix, cell information, and gene annotation. After constructing the CellDataSet object, data normalization and the removal of low-quality genes were performed, using the estimateSizeFactors and estimateDispersions functions for normalization and dispersion estimation. Dimensionality reduction analysis was then conducted using the DDRTree algorithm in Monocle2, achieved through the reduceDimension function. Subsequently, the orderCells function was used to construct cell trajectories based on the pseudotemporal order of the cells, revealing dynamic changes in cells during development. The constructed cell trajectories were visualized using the plot_cell_trajectory function in Monocle2, and a detailed analysis was conducted in conjunction with cell type information. By comparing gene expression changes at different pseudotime points, we explored potential developmental trajectories and transitional relationships. Key nodes and branch points were carefully annotated based on the cell type annotations.

Proportion of cell subpopulation abundance

To assess the abundance of different cell subpopulations within the samples, we analyzed the annotated single-cell RNA sequencing data. First, the Idents function in the Seurat package was used to annotate the cells, ensuring that each cell was assigned to a specific subpopulation. The abundance of each subpopulation was then evaluated by calculating its proportion within the sample. The table function was used to count the number of cells in each subpopulation, and this number was divided by the total cell count of the sample to obtain the abundance proportion of each subpopulation. To visually represent this information, a bar plot was generated using the ggplot2 package, allowing for a clear comparison of the relative abundance of different cell subpopulations within the sample. This approach provides valuable insights for further biological analysis and research into disease mechanisms.

Differential expression analysis using the scde package

To identify differentially expressed genes between different cell types, we performed differential expression analysis using the scde package. First, we extracted the quality-controlled and annotated single-cell RNA sequencing data and converted the Seurat-processed data into the format required by scde, ensuring the inclusion of high-quality cells and highly variable genes. Next, the scde.error.models function was used to fit error models for each cell, estimating sequencing errors and enhancing the accuracy of model fitting through default parameters and recommended preprocessing steps. The scde.expression.prior function was then employed to estimate the prior distribution of gene expression, providing a reference for subsequent differential analysis. Finally, the scde.expression.difference function was used to perform statistical analysis of gene expression differences between different cell populations, identifying significantly differentially expressed genes. Through these steps, we gained insights into the significant gene expression differences between cell types, further elucidating their potential functions and mechanisms in disease.

Expression of specific genes in cell populations

To investigate the expression of specific genes across different cell populations, we conducted a detailed analysis of the single-cell RNA sequencing data. First, the genes of interest were extracted from the preprocessed single-cell data. Then, the FeaturePlot and VlnPlot functions in the Seurat package were used to visualize the expression of these genes in different cell populations. The FeaturePlot function allows for the visualization of specific gene expression on UMAP or t-SNE dimensionality reduction plots, showing the spatial distribution of the gene across different cell populations.

Small molecule agents screening and molecular docking analysis

The primary aim of the connectivity map database (CMap, https://clue.io/) is to explore the functional relationships between genes, small molecule drugs, and diseases, serving as a drug prediction platform based on differential gene expression. To perform molecular docking, the 3D structure of the target molecule is obtained in mol2 format from the PubChem database. Using AutodockTools 1.5.6, the ligand molecule is prepared by adding hydrogen atoms, calculating charges, locating the ligand root, identifying and specifying rotatable bonds, and then saving it as a pdbqt file. The RCSB Protein Data Bank (www.rcsb.org/) provides the core 3D structure of the target protein, which is downloaded and prepared as a docking protein. This preparation involves opening the structure, calculating Gasteiger charges, adding nonpolar hydrogens, defining it as a receptor, and saving it as a pdbqt file with all hydrogen atoms added using AutodockTools 1.5.6. For Vina molecular docking, parameters such as exhaustiveness are set to 15, along with the box size and coordinates, while other settings remain at their defaults. Semiflexible docking is performed using Autodockvina 1.1.2, and the conformation with the highest binding affinity is selected as the final docked conformation [39].

Experimental consumables

The immortalized chondrocytes were purchased from Chongqing Bopei Biotechnology Co., Ltd. (iCell-0092a). We bought IL-1β and Andrographolide (A407762) from Aladdin Company. Fetal bovine serum (FBS) was obtained from Gibco Inc.; the PC-1iCell Primary Chondrocyte Cell Culture System was acquired from HyClone Inc. with an HPLC grade of ≥ 94%. We bought CASP8, CASP1, TLR3, IL18, GAPDH antibody from Biobay Bio. Beijing Solabio Technology Co. was the supplier of the CCK-8 kit (cell proliferation and toxicity detection assay) that was acquired.

Effect of CCK8 assay intervention on cell proliferation

Chondrocyte were grown until they reached 80% confluence and then digested into cell suspensions using a 0.25% trypsin solution. These cells were then seeded into 96-well plates at a density of 3000 cells per well and incubated at 37 °C with 5% CO2 for 24 h to allow for cell attachment. Following this, 200 μL of medium containing varying concentrations of IL-1β (0, 1, 10, 20, 50 ng/mL) was added to each well, and the cells were cultured for an additional 24 h. To assess cell proliferation, 10 μL of CCK-8 solution was added to each well and incubated for 2 h. The optical density (OD) at 450 nm was measured using a multifunctional microplate reader.

Western blot analysis

Proteins (10 μg of total protein per sample) were separated according to their molecular weight using SDS-PAGE gel electrophoresis and subsequently transferred onto a PVDF membrane. This membrane was blocked with 5% skimmed milk for 1 h at room temperature. Following the blocking step, the membrane was incubated overnight at 4 °C with the specific primary antibody. The next day, the membrane was incubated with the secondary antibody for 1 h at room temperature. GAPDH was used as the reference protein for normalization.

Quantitative real-time reverse transcription PCR (qRT-PCR)

To evaluate the mRNA levels of CASP8, CASP1, TLR3, and IL18 in normal, model, and AG treatment groups of chondrocytes, qRT-PCR was performed. Total RNA was isolated from both control and treated groups using a Total RNA Extraction Kit. The extracted RNA was reverse-transcribed into cDNA using the iScript cDNA Synthesis Kit. Relative mRNA levels were determined and normalized to GAPDH using the 2^(-ΔΔCt) method. Quantification was conducted through real-time fluorescence detection on the Bio-Rad CFX96 platform [40]. The following is the list of primer sequences:

Gene	Species	Primer sequence (5′ to 3′)
CASP8	Human	Forward: TTTCTGCCTACAGGGTCATGC
		Reverse: GCTGCTTCTCTCTTTGCTGAA
CASP1	Human	Forward: GCCTGTTCCTGTGATGTGGA
		Reverse: ACTCTTTCAGTGGTGGGCAT
TLR3	Human	Forward: TTGCCTTGTATCTACTTTTGGGG
		Reverse: TCAACACTGTTATGTTTGTGGGT
IL18	Human	Forward: GGTATGGCTGTAACTATCTCTGT
		Reverse: ATGTCACTTTTTGTATCCTTGATGT
GAPDH	Human	Forward: CCAGCAAGAGCACAAGAGGA
		Reverse: TGAGGAGGGGAGATTCAGTGT

In vivo validation in rats

To further validate the efficacy and safety of AG in the treatment of OA, we used adult male SD rats (8–10 weeks old) and established the OA model by injecting sodium iodoacetate. The animal experimental procedures were approved by the Ethics Committee of Southwest Medical University. The specific methods are as follows:

A unilateral intra-articular injection of 50 μL saline containing 2 mg of sodium iodoacetate (8% solution) was administered into the rat knee joint to induce the OA model [41]. Experimental Groups: After the OA model was established, the rats were randomly divided into three groups, with 6 rats in each group: (1) Control group, (2) Model group (received only the OA model), (3) AG treatment group (received Andrographolide Sulfate (10 mg/kg) by intra-articular injection with saline, twice a week for 4 weeks) [42].

Tissue Collection and Histological Analysis: After the experiment, all rats were euthanized, and the knee joint cartilage was collected. The cartilage tissue was fixed in 10% formalin, sectioned, and stained with HE staining to evaluate the morphological changes of the cartilage. Safranin O/fast green staining was used to assess the degradation of cartilage matrix. Safranin O stains the cartilage matrix, and Fast Green stains the connective tissue.

Statistical analysis

GraphPad Prism 9.0 was employed for statistical analysis of the study data. The t-test and Kruskal–Wallis test were utilized to determine statistical significance. Significance levels were denoted as follows: ns for no significance, * for P < 0.05, ** for P < 0.01, *** for P < 0.001, and **** for P < 0.0001.

Results

Identification of differentially expressed genes

Using the Limma rapid differential analysis method, we identified a total of 2291 differentially expressed genes (DEGs). Of these, 918 genes exhibited down-regulated expression, while 1373 genes were up-regulated. To better visualize and analyze the expression patterns of these DEGs, volcano plots and heatmaps were created. These visual representations of OA DEGs, shown in (Fig. 1, Supplement Table.S1), highlight the distinct gene expression patterns and the statistical significance of the differences. These results provide a robust foundation and clear direction for our subsequent research efforts.

Fig. 1

Heatmap and volcano plot for the DEGs identified from the integrated IVDD dataset. A Red and green plot triangles represent DEGs with upregulated and downregulated gene expression. B Each row displays DEGs and each column refers to one of the OA case or control group samples. Red and blue colors represent DEGs with up- and down-regulated gene expression

Weighted gene co-expression network analysis and key module identification

For the weighted gene co-expression network analysis, we determined the appropriate parameter values, setting the β-value to 14 (scale-free R2 = 0.86) as the "soft" threshold based on scale-independence and average connectivity metrics (Fig. 2A, B). Figure 2C displays the gene clustering dendrograms of both the normal and OA groups. From this clustering, we identified 9 distinct gene co-expression modules (GCMs), each represented by a different color (Fig. 2D). Genes that did not fit into any specific module were assigned to the grey module and excluded from further analysis.

Fig. 2

Weighted gene co-expression network (WGCNA) analysis. A, B Screening for soft thresholds based on scale independence and average connectivity. C Gene co-expression modules indicated by different colors under the gene tree. D Heatmap of gene adjacency. E Heatmap of the association between modules and OA. The turquoise module is shown to be correlated significantly with OA. Numbers at the top and bottom brackets represent the correlation coefficient and p-value. F Correlation plot between module membership and gene significance of genes included in the turquoise module

To pinpoint the module with the highest correlation to OA, we assessed the correlation between each gene co-expression module and OA (Fig. 2E). The turquoise module exhibited the strongest correlation with OA (correlation coefficient 0.93, p < 0.01). Consequently, this module, comprising 103 hub genes, was selected for further analysis (Supplement Table.S2).

Subsequently, we evaluated the module member associations between genes within the turquoise module and OA. Our findings indicated a significant positive correlation between the module members and the gene significance (Fig. 2F, r= 0.86, p < 0.01). Therefore, it can be concluded that the turquoise module is significantly associated with OA.

Screening of PANoptosis hub genes and GSEA pathway analysis

To build the protein–protein interaction (PPI) network and identify PANoptosis hub genes in OA (Fig. 3A, Supplement Table.S3), we combined key modular genes from the WGCNA analysis with those from the PANoptosis gene database, GeneCards disease database(Supplement Table.S4) and DEGs, resulting in thirteen intersecting genes (Fig. 3B). These genes were visualized using the STRING (v11.5) database and Cytoscape software. The analysis revealed interactions among twelve genes (Fig. 3C), which were subsequently identified as PANoptosis hub genes. CytoHubba was then employed to evaluate the significance of nodes within the network, utilizing topological characteristics and centrality metrics; node shading indicated their relative importance (Fig. 3D). The expression levels of these twelve genes across different groups are displayed in Fig. 3E.

Fig. 3

PANoptosis hub gene screening. A PANoptosis gene interaction. B Venn diagram showing the intersection between PANoptosis gene, DEGs, GeneCards and key modular genes obtained by WGCNA analysis. C, D PANoptosis hub gene interaction in OA. E Box plot of PANoptosis hub gene expression

To explore the biological pathways involving these PANoptosis hub genes, we conducted single gene set enrichment analysis (ssGSEA) on the four most significant hub genes. The findings indicated that these four hub genes were significantly enriched in pathways related to apoptosis, pyroptosis, necroptosis, and other signaling processes (Fig. 4A). These results suggest that apoptosis, necroptosis, and pyroptosis are likely crucial to the pathogenesis of OA.

Fig. 4

PANoptosis hub gene GSEA pathway analysis and analysis of immune infiltration. A GSEA Pathway Analysis for CASP8, CASP1, TLR3, IL18. B Heat map of immune-infiltrating cells correlation of hub genes. C Immune infiltration stacking plot of hub genes. D The statistical significance of the differences between the pivotal gene immune-infiltrating cell groups is shown

Analysis of immune infiltration results

We assessed the immune infiltration scores of the PANoptosis hub genes using the CIBERSORT algorithm. Figure 4B illustrates the correlations among various immune cell types, while the immune infiltration scores are depicted in a stacked bar plot (Fig. 4C). This plot revealed that the primary immune cells associated with PANoptosis hub genes included Plasma cells, T cells, Mast cells, Macrophages (M0, M1, and M2), and B cells. The statistical significance of these immune cells across different subgroups is shown in Fig. 4D. The analysis indicated that Plasma cells, T cells, and Macrophages (M0, M1, M2) exhibited significant differences (*P < 0.05) among the subgroups. (Supplement Table.S5).

Enrichment analysis

We performed functional enrichment analysis on the 12 intersecting genes to uncover their biological roles and pathways in OA. KEGG pathway analysis revealed that these genes are involved in 85 signaling pathways (Supplement Table S6), with significant enrichment observed in pathways such as Apoptosis, Necroptosis, Cytosolic DNA-sensing pathway, Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, and C-type lectin receptor signaling pathway (Fig. 5A, B). These pathways are crucial in the development and progression of OA.

Fig. 5

Functional enrichment analysis. A, B KEGG pathway analysis of 12 hub genes. C, D GO analysis of 12 pivotal genes for biological processes, cellular components and molecular functions, respectively

GO enrichment analysis indicated that these genes were enriched in 346 GO terms (Supplement Table S7), including 245 biological processes (BP), 50 cellular components (CC), and 51 molecular functions (MF). In the category of biological processes, the top 10 significantly enriched terms included positive regulation of canonical NF-kappaB signal transduction, apoptotic process, pyroptosis, and others. For cellular components, significant terms included mitochondrion, cell body, extracellular region, etc. In terms of molecular functions, notable terms were protease binding, cysteine-type endopeptidase activity, lipopolysaccharide binding, among others (Fig. 5C, D).

Figure 6A, B presents detailed enrichment results of the intersecting genes in the Apoptosis and Necroptosis signaling pathways. These pathways play pivotal roles in the development and progression of OA. This visualization aids in understanding the specific roles and mechanisms of these genes and pathways in OA pathology. The enrichment findings offer valuable insights into the molecular mechanisms of OA and provide direction for future research.

Fig. 6

A Apoptosis signaling pathway in KEGG. B Necroptosis signaling pathway in KEGG

Batch effect correction and definition of cell subpopulations in single-cell data

To ensure accurate analysis, we corrected for batch effects in the six cartilage tissue samples, including three control samples (Normal 1, 2, 3) and three OA samples (OA 1, 2, 3). Figure 7A and B illustrate the sample distribution before and after batch effect correction, respectively. Post-correction, the single-cell data were classified into seven distinct subpopulations (Fig. 7C). We identified the top 5 marker genes for each of the seven subpopulations and visualized them using a bubble plot (Fig. 7D). The key marker genes for these subpopulations are highlighted in Fig. 7E. Based on cell marker analysis(Supplement Table S8) and literature references, we defined the subpopulations as follows: 0—EC (Epithelial Cells), 1—preHTC (Pre-hypertrophic Chondrocytes), 2—RepC (Resting Chondrocytes), 3—HomC (Homeostatic Chondrocytes), 4—ProC (Proliferative Chondrocytes), 5—FC (Fibrocartilage Cells), and 6—RBC (Red Blood Cells) (Fig. 7F).

Fig. 7

Single-cell Data Batch Correction and Definition of Cell Subpopulations. A, B Batch correction analysis of cartilage tissue cells from control groups (Control 1, 2, 3) and osteoarthritis (OA 1, 2, 3) groups. C Seven distinct cell subpopulations identified after batch correction, showing their distribution within the cartilage tissue. D Bubble plot of the top 5 marker genes for each cell subpopulation, displaying gene expression across different subpopulations. E Heatmap showing the expression of key marker genes within each subpopulation, further validating the characteristics of the subpopulations. F Definition of the seven cell subpopulations based on cell markers and literature review, identified as EC (Endothelial Cells), preHTC (Pre-Hypertrophic Chondrocytes), RepC (Reparative Chondrocytes), HomC (Homeostatic Chondrocytes), ProC (Progenitor Chondrocytes), FC (Fibrocartilaginous Cells), and RBC (Red Blood Cells)

Single-cell trajectory analysis and cell subpopulation characteristics

Cell trajectory analysis based on cell type using Monocle2

Cell trajectory analysis was conducted using Monocle2, and the results are shown in Fig. 8A and B. The trajectory analysis revealed that preHTC cells are positioned at a critical node of differentiation, likely playing a significant role in osteoarthritis (OA). As pre-hypertrophic chondrocytes, preHTC cells are crucial in cartilage formation and remodeling. In the pathological environment of OA, these cells may accelerate abnormal proliferation and hypertrophy of the cartilage, leading to structural damage and loss of function in the cartilage tissue. Therefore, the abnormal differentiation of preHTC cells may be a key factor in the progression of OA.

Fig. 8

Analysis of PANoptosis Hub Genes in Osteoarthritis Cartilage Cells. A, B Monocle2 trajectory analysis showing the pseudotime ordering of cells in both control and osteoarthritis (OA) cartilage samples. C Proportional abundance of different cell subpopulations within the OA cartilage tissue, highlighting shifts in cell population distribution compared to control samples. D Differential expression analysis using the SCDE package. The analysis identifies significantly upregulated PANoptosis hub genes, including CASP8, CASP1, TLR3, and IL18, in OA cartilage cells, indicating their potential involvement in disease progression. E Expression patterns of specific genes within the identified cell subpopulations, emphasizing the localization and levels of PANoptosis-related genes across different cell clusters

Analysis of cell subpopulation abundance

Figure 8C illustrates the abundance proportions of different cell subpopulations in both control and OA groups. We observed a significant increase in the proportion of preHTC, FC (fibrocartilage cells), EC (epithelial cells), ProC (proliferating chondrocytes), and RepC (resting chondrocytes) in the OA group, while the proportion of HomC (homeostatic chondrocytes) was markedly reduced. These changes may reflect the remodeling of cell composition during the pathology of OA. For instance, the increase in preHTC cells may be associated with their important role in the progression of OA, whereas the decrease in HomC cells may indicate a disruption in cartilage homeostasis and a decline in regenerative capacity.

Differential expression analysis using the scde package

Differential expression analysis using the scde package revealed a significant increase in the expression of PANoptosis hub genes, such as CASP8, CASP1, TLR3, and IL18, in OA chondrocytes (Fig. 8D, Supplement Table S9). The high expression of these genes may be related to the involvement of OA chondrocytes in pan-apoptosis, further exacerbating cartilage degeneration and pathology. The specific mechanisms underlying this phenomenon require further investigation.

Expression of specific genes in cell populations

Figure 8E shows the expression of specific genes across different cell populations. The results indicate that CASP8 and CASP1 are predominantly highly expressed in the preHTC cell population. This further supports the critical role of preHTC cells in OA, where their high expression may lead to abnormal apoptosis and differentiation, thereby accelerating the progression of OA. TLR3 expression is more dispersed, suggesting that it functions across multiple cell types, while IL18 is highly expressed in various cell populations, indicating its significant role in the global inflammatory response in OA.

In summary, through single-cell trajectory analysis and the study of cell subpopulation characteristics, we have revealed the potential roles of different cell types in OA pathology and identified key pan-apoptosis hub genes. These findings provide important insights for further exploration of the molecular mechanisms and therapeutic strategies for OA.

Results of molecular docking

Based on the screening results, molecular docking was conducted for four PANoptosis hub genes (CASP8, CASP1, TLR3, IL18) with Andrographolide (AG). Binding energy values less than 0 kJ/mol indicate that the small molecule ligand can spontaneously bind to the protein receptor, while values below − 5.0 kJ/mol suggest a stronger binding affinity. The docking simulations produced four results, each with binding energies lower than − 5 kJ/mol, demonstrating favorable binding interactions. The successful docking conformations were visualized using Pymol software (Fig. 9).

Fig. 9

Molecular Docking Analysis of Andrographolide (AG) with PANoptosis Hub Genes. A–D Molecular docking results of Andrographolide (AG) with four PANoptosis hub genes: CASP8 (A), CASP1 (B), TLR3 (C), and IL18 (D). The images depict the binding interactions between AG and the active sites of each protein, demonstrating the potential of AG to modulate these key genes involved in the PANoptosis pathway in osteoarthritis

Therapeutic effects of andrographolide on PANoptosis hub genes

Through extensive bioinformatics analysis, it was revealed that PANoptosis hub genes play a significant role in OA progression. To validate these findings, we conducted in vitro experiments. Figure 10A shows the impact of varying concentrations of IL-1β on chondrocyte viability, assessed using the CCK8 assay. Among the concentrations tested, 10 ng/ml concentration was determined to be the most effective for intervention. Figure 10B displays the expression levels of the four PANoptosis hub genes (CASP8, CASP1, TLR3, IL18) using Western blot analysis. In the model group, these genes were upregulated compared to the normal group. However, in the AG treatment group, their expression levels were significantly reduced compared to the model group. Figure 10C presents qRT-PCR results, which align with the Western blot data, demonstrating that AG treatment effectively downregulates the expression of PANoptosis hub genes. These findings suggest that CASP8, CASP1, TLR3, and IL18 are critical targets for regulating PANoptosis in OA. AG shows potential as a therapeutic agent by ameliorating PANoptosis in chondrocytes, offering a promising avenue for OA treatment.

Fig. 10

Experimental verification. A Effect of IL-1β on the viability of chondrocyte. B Western Blot Detection of Protein Expression in Control, Model and AG Treatment Groups. C qRT-PCR for mRNA expression in control, model and AG treatment groups. *p < 0.05, **p < 0.01, ***p < 0.001. D Histological analysis of cartilage repair in OA rats treated with AG

In vivo validation results in rats

To further validate the efficacy and safety of AG in the treatment of OA, we performed in vivo experiments using a rat OA model. Histological analysis with HE staining and Safranin O/fast green staining showed the following results:

HE staining results showed that, compared to the control group, the cartilage in the knee joints of the model group exhibited significant degeneration and damage, with an irregular cartilage surface, disorganized cell arrangement, noticeable thinning of the cartilage, and the presence of cartilage fissures. In contrast, the cartilage structure in the AG treatment group was relatively intact, with a smoother surface, organized cell arrangement, and partial restoration of cartilage thickness. The degree of cartilage degeneration was significantly reduced, and the fissures were markedly decreased, indicating that AG could effectively slow down cartilage damage and promote repair.

Safranin O/Fast Green staining results showed that, compared to the control group, the cartilage matrix in the model group exhibited significant decolorization, indicating severe degradation of the cartilage matrix, with a noticeable reduction in the red staining intensity. In contrast, the cartilage matrix in the AG treatment group was better preserved, with stronger Safranin O staining, indicating that the cartilage matrix was relatively intact. The preservation of the matrix was significantly better, suggesting that AG can slow down the degradation of the cartilage matrix and promote its repair (Fig. 10D).

Discussion

Osteoarthritis (OA) is a degenerative joint disease that significantly impacts the quality of life of affected individuals, particularly in the elderly, with millions of people worldwide suffering from it [43]. OA is characterized by the progressive degradation of articular cartilage, leading to pain, stiffness, and loss of joint function. The pathogenesis of OA is complex, involving mechanical stress, inflammation, and metabolic factors [44]. A key aspect of OA pathology is cell death, including various forms such as apoptosis, pyroptosis, and necroptosis. These forms of cell death are collectively referred to as PANoptosis, encompassing critical processes that lead to chondrocyte loss and cartilage degeneration [45].

PANoptosis may play a crucial role in OA by mediating both cell death and inflammatory responses. In chondrocytes, the occurrence of PANoptosis exacerbates cartilage damage and promotes the progression of OA. Although the significance of PANoptosis has been recognized in other diseases [46, 47]. Research specifically targeting its role in OA remains limited. Understanding the mechanisms of PANoptosis in OA could provide novel therapeutic targets for this disease, highlighting the need for further investigation in this area.

This study aims to elucidate the role of PANoptosis in OA through integrated bioinformatics combined with single-cell sequencing analysis. We identified 2,291 differentially expressed genes (DEGs) and conducted a cross-analysis of these DEGs with the PANoptosis gene database and key module genes obtained from WGCNA analysis. From this, we identified four hub genes: CASP8, CASP1, TLR3, and IL18. Functional enrichment analysis, including GSEA, KEGG, and GO, revealed that these hub genes are primarily enriched in apoptosis, pyroptosis, necroptosis, and related signaling pathways. Further studies indicated that OA patients exhibit higher levels of immune cell infiltration, including plasma cells, T cells, and macrophages, which release large amounts of inflammatory factors such as TNF-α and IL-1β. This exacerbates the inflammatory cascade and extracellular matrix (ECM) degradation, promoting the progression of OA. We hypothesize that these hub genes play a significant role in the immune-inflammatory response in OA.

The four PANoptosis hub genes identified in this study play crucial roles in regulating the progression of PANoptosis. CASP8 is not only involved in the initiation of the apoptotic pathway but also plays a key role in pyroptosis and necroptosis. The activation of CASP8 can induce apoptosis through caspase-3 and caspase-7 while also regulating necroptosis through interactions with RIPK1 and RIPK3. Additionally, CASP8 inhibits the cleavage of GSDMD, the main executor of pyroptosis, thereby indirectly participating in the pyroptotic pathway. STRING analysis results indicate that CASP8 is central to pyroptosis, apoptosis, and necroptosis [14]. CASP1, as the core molecule of pyroptosis, leads to the formation of pores by cleaving GSDMD, causing cell membrane rupture and the occurrence of pyroptosis. CASP1 also activates IL-1β and IL-18, further driving the inflammatory response [48]. TLR3, an important pattern recognition receptor, induces intracellular inflammation and antiviral responses by activating the NF-κB and IRF3 pathways, and it can indirectly promote apoptosis and necroptosis [49]. IL-18, cleaved by caspase-1 during pyroptosis, is a potent pro-inflammatory cytokine that exacerbates local inflammation by promoting the activation of T cells and NK cells [50]. Collectively, the multifaceted roles of these genes in PANoptosis suggest that targeting these key molecules and simultaneously blocking multiple pathways may offer a more effective treatment for OA than approaches targeting a single pathway. This provides new therapeutic insights for the treatment of OA.

Using the Connectivity Map (CMap) database, we predicted potential therapeutic drugs for OA and identified andrographolide (AG) as having the highest correlation and binding affinity with the hub genes. Studies have shown that treating chondrocytes with 0–20 μM of AG for 24 h is the ideal dose and duration for in vitro anti-OA research. Consequently, we performed Western blot (WB) and quantitative real-time PCR (qRT-PCR) experiments using 20 μM AG. The results revealed that AG significantly reduced the expression levels of the four PANoptosis hub genes—CASP8, CASP1, TLR3, and IL18—in the model group. These findings suggest that AG effectively ameliorates PANoptosis in chondrocytes, making it a promising therapeutic agent for OA.

In comparison with current OA therapeutics, such as NSAIDs, corticosteroids, and biologics, Andrographolide (AG) demonstrates several advantages. First, AG has multi-target effects, unlike NSAIDs that primarily inhibit COX enzymes or corticosteroids with broad anti-inflammatory actions. AG targets multiple pathways (e.g., NF-κB, inflammasome) involved in PANoptosis, offering broader mechanistic coverage. Additionally, AG has shown a favorable safety profile compared to NSAIDs, exhibiting lower gastrointestinal toxicity, and it lacks the immunosuppressive risks seen with biologics (e.g., TNF-α inhibitors). AG’s anti-inflammatory and anti-apoptotic properties could complement biologics or reduce NSAID dosage, thus minimizing side effects. This highlights the translational relevance of AG as a potential OA therapeutic.

A potential limitation of AG is its low bioavailability and rapid metabolism, which may limit its efficacy in clinical applications. However, several strategies can be employed to overcome this challenge. For instance, nanoformulations using lipid-based nanoparticles or polymeric carriers like PLGA could enhance AG’s solubility and stability, improving its bioavailability. Additionally, prodrug strategies (e.g., esterification) could prolong AG’s half-life, and combination therapies with bioavailability enhancers, such as piperine, may further improve its absorption. These strategies suggest that, despite its pharmacokinetic challenges, AG has significant potential for clinical application, and ongoing research is focused on optimizing its delivery.

In conclusion, our study underscores the critical role of PANoptosis in OA and identifies key hub genes that may serve as novel therapeutic targets. The results highlight the importance of targeting the PANoptosis pathway and the potential of using AG as a therapeutic approach for OA. This research provides foundational insights that may facilitate the development of more effective OA treatment strategies, ultimately improving patient outcomes.

Conclusion

In this study, we explored the role of PANoptosis in OA and evaluated the potential of AG as a therapeutic agent for OA. Our findings highlight the significant role of PANoptosis in the pathogenesis of OA, with key genes such as CASP8, CASP1, TLR3, and IL18 playing central roles in mediating cell death and inflammatory responses. These genes are primarily enriched in key pathways related to apoptosis, pyroptosis, and necroptosis, suggesting that targeting PANoptosis could provide a promising therapeutic strategy for OA.

Through in vitro experiments, we demonstrated that AG significantly reduced the expression of these PANoptosis hub genes, effectively alleviating PANoptosis in chondrocytes. Additionally, in vivo experiments further confirmed AG's therapeutic potential. Histological analysis showed that AG treatment promoted cartilage repair, reduced cartilage degeneration, and protected the cartilage matrix, supporting its effectiveness in slowing OA progression.

Compared to traditional OA therapeutics, such as NSAIDs and corticosteroids, AG also exhibited a favorable safety profile. AG has multi-target effects, targeting pathways associated with inflammation and cell death, with fewer side effects, particularly in terms of gastrointestinal toxicity. Although AG's low bioavailability remains a challenge, strategies such as nanoformulations, prodrug approaches, and the use of bioavailability enhancers may improve its clinical applicability.

Overall, our study provides insights into the role of PANoptosis in OA and highlights AG as a potential multi-target therapeutic agent. The findings offer new perspectives for developing more effective and safer OA treatment strategies, with significant potential to improve patient outcomes.