Introduction
Ankylosing spondylitis (AS) is an autoimmune disease characterized mainly by chronic progressive inflammation of the axial joints. Its core pathological changes are sacroiliitis and spondylitis [1]. AS a typical seronegative spondyloarthropathy, the core clinical manifestations of AS include latent low back pain accompanied by morning stiffness, progressive limitation of spinal movement, and in the advanced stage, it can lead to spinal ankylosis and kyphosis deformity, significantly damaging the patient’s function and quality of life [2, 3].
In addition to the axial joints, extra-articular involvement of peripheral joints, attachment points (attachment points of tendons, ligaments and bones), eyes (uveitis), intestines, cardiovascular systems, etc. is also not uncommon [4]. Epidemiological data show that AS has significant regional and ethnic differences. The prevalence rate in the Chinese population is approximately 0.3%, and it is more common among young and middle-aged men, with a notable familial aggregation [5].
The etiology and pathogenesis of it have not been fully clarified so far. Currently, it is generally believed that it is the result of the combined effect of multiple factors such as genetic susceptibility (such as strong association with HLA-B27), environmental factors (such as intestinal flora imbalance), and abnormal immune responses (including innate immunity and adaptive immunity imbalance) [6]. Although biological agents represented by tumor necrosis factor -α (TNF-α) inhibitors and interleukin-17 (IL-17) inhibitors have greatly improved the treatment pattern of AS in recent years, significantly alleviating inflammation and clinical symptoms, However, the insufficient response rate of drugs, potential side effects of long-term use, high costs, and the inability to completely prevent the aggravation of spinal ankylosis caused by new bone formation (osteophytes, ligament ossification) in some patients remain severe challenges currently faced in clinical practice, highlighting the urgency of in-depth exploration of new mechanisms of AS pathogenesis and the search for new therapeutic targets [7].
In the cutting-edge field of exploring the pathogenesis of immune inflammatory diseases, programmed cell death, especially PANoptosis which has been deeply recognized in recent years, is receiving extensive attention [8]. PANoptosis does not refer to a single, independent pathway of cell death, but rather a highly coordinated inflammatory programmed cell death pattern that integrates the key molecular features of Apoptosis, Pyroptosis and Necroptosis. Its core molecular mechanism involves the assembly of a large multi-protein complex called “PANoptosome” [9–11].
This complex can simultaneously recruit and activate apoptosis-related caspase-enzymes (such as Caspase-3 and Caspase-8), pyroptosis related Caspase-1/GSDMD pathways, and necrosis related RIPK1/RIPK3/MLKL pathways [12]. Once triggered, PANoptosis not only leads to cell membrane rupture and contents release, but also drives the explosive production of a large number of pro-inflammatory factors (such as IL-1β, IL-18 and TNF-α) and damage-associated molecular patterns (DAMPs) [13, 14]. This intense “inflammatory death” pattern makes PANoptosis, while eliminating infection sources or abnormal cells, highly prone to triggering severe local or systemic inflammatory storms.
It has been confirmed to play a key role in the pathological processes of various infectious diseases, autoinflammatory diseases and autoimmune diseases. Its unique integration and powerful pro-inflammatory ability make it a new window for understanding the mechanisms of chronic inflammation and tissue damage.
Given that the core pathological features of AS are uncontrollable chronic inflammation, persistent immune cell infiltration, and the subsequent abnormal bone destruction and new bone formation, exploring the role of PANoptosis in AS has a solid theoretical basis and potential significant implications. Excessive activation and dysfunction of innate immune cells and adaptive immune cells exist both locally in the joints and in the circulation of AS patients [15, 16]. Although scattered studies have suggested abnormal expression of molecules related to apoptosis, pyroptosis or necrosis in AS and other types of osteoarthritis, systematic research on PANoptosis, an integrated death pattern, is still a blank in the field of AS [17].
Its exact contribution to the occurrence and development of the disease, regulatory network and its connection with key clinical phenotypes need to be clarified urgently.
This study aims to systematically screen and identify the key PANoptosis-Related Genes (PRGs) in the pathogenesis of AS by integrating bioinformatics analysis and experimental verification. This research not only holds the promise of revealing a new mechanism of PANoptosis in the occurrence of inflammation and ossification in AS, providing a fresh perspective for a deeper understanding of the immunopathology of AS, but more importantly, the core PRGs identified and the regulatory network they constitute It has the potential to become a novel biomarker for the diagnosis of AS, a predictive indicator for assessing disease activity or ossification risk, and a novel therapeutic target for developing precise targeted intervention strategies. This will provide important scientific basis and innovative directions for overcoming the current treatment predicament of AS, achieving more effective disease control and improving the long-term prognosis of patients (Fig.
1).
Fig. 1The workflow of this study
Materials and methods
Screening of PANoptosis genes in AS
Data collection
“Ankylosing spondylitis” was searched as a keyword in the comprehensive gene expression database (https://www.ncbi.nlm.nih.gov/geo/). Four datasets were ultimately selected for the analysis of this study, including GSE25101 (from the GPL6947 platform, with 16 samples of AS patients and 16 samples of healthy controls), and GSE73754 (from the GPL10558 platform, with 52 samples of AS patients and 20 samples of healthy controls). GSE11886 (from the GPL570 platform, including 15 samples of AS patients and 18 samples of healthy controls), GSE134290 (from the GPL570 platform, including 3 samples of AS patients and 3 samples of healthy controls).
Among them, GSE25101 and GSE73754 were used as training datasets, while GSE11886 and GSE134290 were employed as external validation datasets. The collection of PANoptosis-related genes (PRGs) mainly originated from a self-built database (see Appendix 1). Specifically, these PRGs mainly comes from Genecard (https://www.genecards.org/), NCBI (https://www.ncbi.nlm.nih.gov/) and uniprotKB (https://www.uniprot.org/) retrieval, finally received a total of 157 CRGs.
Differentially expressed genes (DEG) analysis
The GSE25101 and GSE73754 datasets were used by the limma package of R language (version 4.4.1) for DEG filtering. Based on the grouping information of the samples, a matrix is constructed and designed. After fitting and expressing the data with lmFit(), empirical Bayesian correction is performed through eBayes().
The screening of significantly different genes was conducted using dual thresholds. After adjustment, the p-value (adj.P.Val) was less than 0.05 and the absolute log₂ folding change (logFC) was greater than 1. The visualization of results includes the use of “ggplot2” to draw volcano maps and heat maps.
Weighted gene co-expression network (WGCNA)Analysis
The GSE25101 and GSE73754 datasets were used by the R package WGCNA (v1.72) for WGCNA analysis, and weighted gene co-expression networks were respectively constructed. Firstly, genes with expression variance values higher than the 75% quantile were screened for subsequent analysis. The optimal soft threshold β was determined through scale-free topological fitting (fitting index threshold > 0.8), and co-expression modules were identified using the dynamic tree cutting algorithm (minimum module size = 20, cutting sensitivity = 2).
Merge redundant modules with feature vector similarity greater than 0.75; Furthermore, the association between the module characteristic gene (ME) and the target trait was calculated. A Pearson correlation coefficient (|r| > 0.3 and p < 0.05) was used as the condition for screening significantly correlated modules [18].
Differential PANoptosis-related gene screening of AS (APRGs)
DEG and the most significant genes related modules and extensive apoptosis related genes have been uploaded to the bioinformatics and evolutionary genomics platform screened APRGs, and draw the Venn diagram.
Protein-protein interaction (PPI) and enrichment analysis
APRGs was imported into the STRING database (https://cn.string-db.org), with the minimum interaction threshold set to “Highest confidence > 0.9”, and the tsv format file of the default interaction node was downloaded. The tsv file was imported into Cytoscape to remove free nodes. The CytoHubba algorithm was used for the screening of Hub genes and to construct a visual protein-protein interaction (PPI) network.
In addition, the APRGs were uploaded to the DAVID database (https://david.ncifcrf.gov/tools.jsp) for gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. Based on the condition of P-value < 0.05, the biological processes (BP), molecular functions (MF), cellular components (CC), and pathway enrichment results were selected, and the enrichment analysis maps were drawn using R language (https://david.ncifcrf.gov/tools.jsp). Based on miRwalk (http://mirwalk.umm.uni-heidelberg.de/), TRRUST database (https://www.grnpedia.org/) screening Hub - gene transcription regulatory factors (TF) and associated micrornas.
Machine learning analysis
For pre-screened APRGs, multiple machine learning is adopted to screen key biomarkers. Firstly, through random forest (RF) analysis, the Gini coefficient importance of genes was iteratively calculated using 500 decision trees, and the top 10 genes in terms of importance scores were retained. A stratified sampling strategy was adopted in the least absolute shrinkage and selection operator (LASSO) logistic regression analysis.
The internal dataset was randomly divided into a training set and a test set in a ratio of 7:3 to ensure consistent proportions of cases and controls in both groups. On the training set, a 10-fold cross-validation procedure combined with three repetitions was employed for model training and hyperparameter optimization. Finally, gene sequencing is carried out using the Support Vector Machine (SVM) with radial basis kernel function recursive feature elimination.
By integrating the results of the three algorithms, the top 10 hub genes that repeatedly occur and rank high in the three methods are extracted AS the core candidate markers for PANoptosis regulation in AS. To evaluate the statistical significance of the model performance, a permutation test was conducted. The specific procedure was as follows: the same 10-fold cross-validation process was applied on the training set, but the outcome variables (such as AS/Control labels) were randomly shuffled 1,000 times.
After each shuffle, the model was retrained and its performance metrics (AUC-ROC) were recorded.
Receiver operating characteristic (ROC) analysis
The core APRGs screened by machine learning were subjected to ROC analysis by the pROC package (v1.18.0) in the GSE11886 and GSE134290 datasets to verify their diagnostic efficacy. Firstly, determine the cut-off value of the optimal expression level of each gene through the maximum Youden index; Then calculate the area under the curve (AUC) and its 95% confidence interval (with 500 repeated samples using the Bootstrap method). The evaluation criteria are AUC > 0.6 (good discriminative ability) and p value < 0.05.
The collaborative diagnostic performance of the hub gene set was ultimately verified through multi-gene combined ROC analysis. To eliminate batch effects, the expression matrices of the internal training set and external validation sets were merged and corrected using the ComBat algorithm, with the data source specified as the batch variable and disease status (AS/Control) retained as the biological variable. Based on the corrected data, features were re-extracted, and the final model was applied to the validation sets for performance evaluation.
The effectiveness of batch effect correction was validated through principal component analysis (PCA).
Gene set enrichment analysis (GSEA)
Convert the probe names in the series expression matrix of core APRGs to Gene names, eliminate abnormal and non-gene name-free data, and then upload to gene set enrichment analysis (GSEA) software (v3.0) for GSEA analysis of core APRGs. Based on the gene expression profile and phenotypic grouping, the minimum gene set was set to 5, the maximum gene set was set to 500, and 1000 resampling was conducted. The P value < 0.05 and FDR of < 0.25 were used as limiting conditions to screen the related pathways of core APRGs.
Immune infiltration analysis
The immune microenvironment of the pre-screened APRGs is analyzed by using the deconvolution algorithm (CIBERSORT). Based on the standard expression profile of APRGs, the proportion of various immune cells in the sample was calculated through deconvolution using the LM22 feature matrix (including 22 immune cell subtypes). Data preprocessing included normalization of expression quantiles and batch effect correction.
Furthermore, for each differentially expressed PANoptosis gene, calculate the rank correlation coefficient between its expression level and the abundance of immune cells to identify the significantly associated immune cell subpopulations. Ultimately, the key gene-immune cell interaction network is visualized through heat maps [19]. The tissue samples involved in this study may contain different proportions of target and immune cells.
We adopted the following steps to ensure the quality of the data input to CIBERSORT and the purity of the filter. Firstly, we assess the quality of RNA based on the expression levels of the internal reference genes. Secondly, we examined the P-value output of CIBERSORT for each sample and retained only the samples with a P-value < 0.05 to ensure that the estimation of the proportion of immune cells was statistically significant.
A total of 104 samples that passed quality control were included in the final analysis.
Verification of genes related to PANoptosis
Animal modeling and grouping
Ten 6-week-old male Wistar rats, weighing 200 ± 20 g, were provided by the Experimental Animal Center of Heilongjiang University of Chinese Medicine and placed in an SPF-grade animal room for adaptive rearing for one week. Based on the power analysis and the actual usage of animal tissues, we finally determined 5 rats in each group. The rats were divided into the Sham group and the AS group according to the random number table method, with 5 rats in each group.
Referring to previous reports [20], after anesthesia, the bilateral sacroiliac joint areas of the AS group rats were accurately located, and 0.05 mL of well-mixed complete freund’s adjuvant (CFA) suspension was injected on each side. The rats in the Sham group were injected with the same amount of normal saline for 4 consecutive weeks. Postoperative analgesia treatment must be given and the condition of the animals must be closely observed.
This animal experiment was approved by the Animal Experiment Ethics Committee of Heilongjiang University of Chinese Medicine (DWSYXK-20250710). All procedures involving the use of experimental animals comply with ethical requirements and follow the ARRIVE principle. After the modeling was completed, the abdominal aortic blood of the rats was collected after anesthesia with an excess of pentobarbital sodium.
Subsequently, the neck was severed and the rats were sacrificed, and the spinal tissues were rapidly dissected. Some of them were fixed in 4% paraformaldehyde neutral buffer for 24 h and then used for histopathological examination, while the other part was stored at -80℃ and used for the detection of other items.
Tail suspension test and arthritis score
To evaluate the spinal stiffness of rats in each group, the tail suspension test was used to assess the joint range of motion of rats. Fix the tail end of the rat on a horizontal bar, allowing its body to be naturally suspended in the air. Keep it suspended for 90 s and record the percentage of the time the body remains completely stiff to the total suspension time.
An increase in the proportion of rigidity time is regarded as a sign of reduced spinal range of motion and aggravated stiffness. Before and after the injection of CFA in rats, the manifestations of arthritis such as redness and swelling of peripheral joints in each group of mice were observed respectively. The severity of arthritis in the limbs of mice was scored using the internationally common arthritis index (AI).
No inflammatory swelling is scored as 0 points, swelling of one toe as 1 point, swelling of more than one toe as 2 points, joint stiffness and loss of function as 3 points, ankle swelling and stiffness as 4 points. The total score ranges from 0 to 16 points. All evaluations of behavioral and histological outcomes were conducted by researchers who were unaware of the groupings.
Hematoxylin-eosin (HE) staining
The fixed spinal samples were decalcified with 10% EDTA solution, and then embedded after rinsing with running water, gradient dehydration, transparency and paraffin impregnation. Prepare 5 μm continuous sections and place them in an oven at 60℃. Then, successively immerse them in xylene solution for dewaxing (twice, 10 min each) and gradient ethanol rehydration (100%, 95%, 90%, 80%, 70%).
The staining was carried out for 5 min each, followed by hematoxylin staining for 5 min, differentiation in 1% hydrochloric acid ethanol solution for 3 s, 0.5% ammonia water for 30 s, and 0.5% eosin staining for 1 min. After gradient ethanol dehydration and xylene transparency, the tissue was finally sealed with neutral gum and the pathological features of the spinal tissue were observed under an optical microscope.
Immunohistochemical (IHC) staining
Based on the operations in Section "Tail suspension test and arthritis score", paraffin sections of the paraspinal muscle tissue were dewaxed and treated with gradient ethanol, and then repaired with sodium citrate. Incubate at room temperature with 3% H2O2 for 15 min, block the non-specific site with 5% BSA for 30 min, then add RIPK1 primary antibody (1:800, beyotime (Shanghai) Trading Co., Ltd, Cat: AF7896) and incubate overnight in a 4 ° C wet chamber; Subsequently, the secondary antibody was added and incubated at room temperature for 1 h. After adding DAB chromogenic solution, it was rinsed with distilled water.
Hematoxylin staining was performed for 1 min, followed by differentiation in 0.1% ethanol hydrochloric acid for 5 s. After dehydration and transparency were completed, it was sealed with neutral gum, and the positive area was calculated using Image J.
Immunofluorescence (IF) staining
The sections of spinal tissue, after dewaxing and gradient rehydration, were placed in an EDTA microwave oven and heated for 10 min to complete the antigen retrieval. After blocking with 5% donkey serum for 30 min, ZBP1 (1:500, Wuhan Sanying Biotechnology Co., Ltd, Cat: 13285-1-AP) and AIM2 (1:500, Abcam (Shanghai) Trading Co., Ltd, Cat༚ab204995) primary antibodies were added and incubated in a 4 °C wet box away from light for 24 h. The next day, the fluorescently labeled secondary antibody was incubated at room temperature in the dark for 1 h, and DAPI (Shanghai Beyotime Biotechnology Co., Ltd, Cat༚C1002) staining was performed for 5 min.
Anti-fluorescence quenching mounting agents are used for sealing. The co-expression of AIM2 and ZBP1 was observed under a fluorescence microscope (Zeiss LSM 900), and the fluorescence intensity was analyzed using Image J.
Safranin O and fast green (SO/FG) staining
The paraffin-embedded spinal tissue sections were baked, followed by routine deparaffinization and hydration through a graded ethanol series. Staining was then performed using Safranin O and Fast Green (Shanghai Beyotime Biotechnology Co., Ltd, Cat: C0621S): the sections were first stained with fast green solution for 5 min, rinsed briefly in tap water until the cartilage became colorless, and subsequently differentiated in 1% acetic acid solution for 30 s. After a quick wash in tap water, the sections were counterstained with sSafranin O solution for 8 min.
Following staining, the sections were rapidly dehydrated through three changes of absolute ethanol, cleared in xylene three times, and finally mounted with neutral resin.
ELISA
The blood samples obtained in animal modeling and grouping section were centrifuged at 15,000 r/min for 10 min, and the supernatant was collected. The contents of TNF-α, IL-6, IL-10 and IL-4, procollagen I n-terminal Propeptide (PINP), type I collagen c-terminal peptide (CTX-I) (Wuhan Saipei Biotechnology Co., Ltd, Cat: SP12250, SP12279, SP12294, SP12281, SP13375, SP12422) in the serum of rats in each group were detected according to the operation instructions of the kit. The absorbance was detected by an enzyme-linked immunosorbent assay (ELISA) reader, and the expression levels of inflammatory factors in the serum were calculated.
qPCR
100 mg of mouse spinal tissue was rapidly ground into powder in a mortar pre-cooled with liquid nitrogen, and TRIzol was added to extract total RNA. Take 1 µg of total RNA for reverse transcription using the SuperScrip IV kit (Thermo Fisher Scientific Inc., Cat: 12594100) (37 °C, 15 min; Synthesize cDNA at 85℃ for 5 s). Specific primers were designed for the target genes AIM2, CASP8, ADAR, ALKBH5, IFNG and TNF using Primer-BLAST.
Amplification was carried out in the CFX96 Touch system (Bio-Rad) using SYBR Green Master Mix (Thermo Fisher Scientific Inc., Cat: A46110), with pre-modification at 95 °C (3 min), followed by 40 cycles (95℃, 10s; 60℃, 30s); The amplification specificity was finally verified through melting curve analysis at 60℃ for 30 s. The relative expression levels of genes were calculated by the 2−ΔΔCt method. Primer design and synthesis are provided by Nanjing Genscript Biotechnology Co., LTD (Table 1).
Table 1Primer sequenceGeneSequenceSIZE(bp)AIM2F: TGGAGACCTTCGTGACCCTTAR: GGTCTTCGGGTTCTTGTTGTC196CASP8F: CCTGGGACTGTGGAGGTCTTR: GCTCACAAACGTGCTGGTAA154ADARF: GGGAGGAAGCTGAGGTGATGR: TCTGGCTGCTCTTCATCTGG132ALKBH5F: CCTGGGATCTGCCTCTACCTR: GGTCTTGGCATAGGCGTAGA119IFNGF: AGCAACAGCAAGGCGAAAAR: CTGGACCTGTGGGTTGTTGA148TNFF: CTCCAGGCGGTGCCTATGR: GAGGCCATTTGGGAACTTCT173GAPDHF: ACCACAGTCCATGCCATCACR: TCCACCACCCTGTTGCTGTA138
Western blot
Weigh 100 mg of mouse spinal tissue and cut it into small pieces. Add RIPA lysis (Shanghai Beyotime Biotechnology Co., Ltd, Cat: P0013C) buffer for homogenization and lysis, and collect the supernatant. The protein concentration was measured by the BCA method.
After adding the sample loading buffer, it was mixed evenly and heated for denaturation. SDS gel constant pressure electrophoresis, wet transfer membrane onto polyvinylidene fluoride membrane. Seal with 5% skimmed milk powder solution for 2.5 h, and rinse with TBST twice, each time for 5 min.
After adding AIM2, CASP8, ADAR, ALKBH5, IFNG and TNF primary antibodies respectively, incubate at 4 °C overnight, with the antibody dilution ratio being 1:1000 (Abcam (Shanghai) Trading Co., Ltd, Cat: ab204995, ab264958, ab307585, ab195377, ab267369, ab183218). Subsequently, PVDF membranes were incubated with HRP-labeled secondary antibodies at room temperature for 1 h. Finally, ECL chromogenic solution was added for development, and images were collected using a gel imager.
Taking β-actin as the internal reference control, the ratio of the gray values of the target protein to those of the internal reference protein was calculated using Image J software.
Statistical analysis
The experimental data were all expressed as mean ± standard deviation (mean ± SD), and statistical analysis was conducted using SPSS 23.0 and GraphPad Prism. Repeated measures analysis of variance was used to compare data such as the tail suspension test and arthritis score between groups. For the comparison among multiple groups of other single measurement data, One-way ANOVA analysis was used, and P < 0.05 was considered statistically significant.
Results
Results of DEGs screening in AS
Firstly, we conducted DEGs analysis on the GSE25101 and GSE73754 datasets with p-value < 0.05 and log2FC > ± 1 as the screening conditions to clarify the changes of differentially expressed genes in AS (Fig. 2). The results showed that a total of 530 DEGs were identified in the GSE25101 dataset, among which 202 genes were up-regulated and 328 genes were down-regulated.
A total of 1555 DEGs were identified in the GSE73754 dataset, among which 766 genes were up-regulated and 789 genes were down-regulated. Furthermore, we found that even after merging the two datasets, the batch effect correction results were excellent. Secondly, we integrated the DEGs from the two datasets and demonstrated the biological functions and pathways of these DEGs through enrichment analysis (Fig.
3). There are a total of 500 items in the GO analysis, among which 224 items are BP, mainly related to aerobic electron transport chain. positive regulation of NF-κB transcription factor activity, ATP synthesis coupled electron transport; CC has a total of 49 entries, mainly involving mitochondrial respirasome, respiratory chain complex, respirasome; MF has a total of 38 entries, mainly involving NF-kappaB binding, oxidoreduction-driven active transmembrane transporter activity, and methyl-CpG binding. KEGG analysis revealed that these DEGs were enriched in 20 pathways, mainly involved in Huntington disease, Oxidative phosphorylation, and Thermogenesis.
Fig. 2DEGs Analysis of GSE25101 and GSE73754. A DEG volcano plot of the GSE25101 dataset; (B) DEG heat map of the GSE25101 dataset; (C) DEG volcano map of the GSE73754 dataset; (D) DEG heat map of the GSE73754 dataset; (E) PCA analysis after merging the GSE25101 and GSE73754 datasets
Fig. 3Enrichment Analysis of DEGs. A GO analysis BP entries of DEGs; (B) DEGs’ GO analysis of CC entries; (C) GO analysis of MF entries in DEGs; (D) KEGG analysis of DEG
The results of screening the most significantly relevant modules in AS
WGCNA is an analytical method for analyzing the gene expression patterns of multiple samples. It can cluster genes with similar expression patterns and analyze the association between modules and specific traits or phenotypes. Therefore, after obtaining the DEGs, we conducted a WGCNA analysis to further explore and identify the characteristic gene modules of the AS datasets GSE25101 and GSE73754, and to explore the association between the gene network and the phenotype (Fig.
4). The results show that in the GSE25101 dataset, the first value to reach the scale-free topological fitting index (0.8) is 12, thus the power value is 12 (Fig. 4A), and 14 modules are identified (Fig.
4B), among which module green contains 53 genes and has the highest correlation with module traits (0.49). And there was a significant positive correlation (Fig. 4C).
In the GSE73754 dataset, the first value to reach the scale-free topological fitting index (0.8) was also 4, so the power value was selected as 4 (Fig. 4E), and 12 modules were identified (Fig. 4F), among which module tan contains 142 genes and has the highest correlation with module traits (0.56).
And there was a significant positive correlation (Fig. 4G). Furthermore, by integrating DEGs and the most significantly related module genes, we ultimately obtained a total of 31 core DEGs (Fig.
4D), and these core DEGs were highly expressed in AS patients (Fig. 4H).
Fig. 4WGCNA analysis of the AS dataset. A The power value of the GSE25101 dataset is 12; (B) Based on the selected power values, a weighted co-expression network model was established, and ultimately 846 genes were divided into 14 modules; (C) The pearson correlation algorithm was used to calculate the correlation coefficient and p-value between the characteristic genes and traits of the module, among which the correlation of the MEgreen module was the highest; (D) DEGs of the AS dataset and venn plots of the most significantly related module genes; (E) The power value of the GSE73754 dataset is 4; (F) Based on the selected power values, a weighted co-expression network model was established, and ultimately 10,555 genes were divided into 12 modules; (G) The pearson correlation algorithm was used to calculate the correlation coefficient and pValue between the characteristic genes and traits of the module, among which the MEtan module had the highest correlation; (H) Heat map of the expression of 31 core DEGs in the AS dataset
Results of APRGs screening
After obtaining 31 core DEGs, we compared them with our self-built database of PANoptosis-related genes and ultimately identified 23 APRGs (Fig. 5A). 23 APRGs were uploaded to the STRING database for PPI analysis.
The results showed that a total of 23 nodes and 99 edges were obtained, and the average node degree was 8.61 (Fig. 5B). We selected the top 10 values through the degree, MCC, Betweenness and Closeness algorithms of cytohubba as the screening conditions for Hub genes.
Eventually, a total of 6 Hub genes were screened out, including AIM2, CASP8, NLRP3 and IRF1. TNF and IFNG (Fig. 5C), and these six Hub genes are highly expressed in patients with AS (Fig.
5D). Subsequently, based on TRRUST, miRTarbase and Mirwalk, we predicted the regulatory relationship between TF and miRNA of six Hub genes. The results showed that the TFS of these Hub genes were STAT4, STAT1, CREB1, STAT3, RELA and NFKB1.
And 163 mirnas were matched (Fig. 5E).
Fig. 5Identification of APRGs. A After screening through the PANoptosis-related gene database, a total of 23 APRGs were identified in the core DEGs of AS; (B) PPI maps of 23 APRGs after free nodes were removed; (C) The MCC, degree, closeness, and betweenness algorithms based on cytohubba ultimately obtained a total of 6 Hub genes.
D Heat map of the expression of 6 Hub genes in AS; (E) The “TF-mRNA-miRNA” network of six Hub genes
Results of APRGs enrichment analysis
To further clarify the biological functions and pathways of these 23 APRGs, we conducted GO and KEGG analyses on these APRGs through the DAVID database. Under the condition of P < 0.05, a total of 134 items were obtained in the GO analysis. We found that these APRGs mainly involved 107 items in the BP analysis of GO, including positive regulation of inflammatory response, positive regulation of interleukin-1 beta production and defense response to virus (Fig.
6A); In the CC analysis, 12 items were mainly involved, including cytosol, cytoplasm and nucleoplasm (Fig. 6B); The MF analysis mainly involves 15 items, including identical protein binding, cysteine-type endopeptidase activator activity and signaling adaptor activity (Fig. 6C).
In the KEGG analysis, it was found that these APRGs mainly involved 58 signaling pathways, mainly NOD-like receptor signaling pathways, Cytosolic DNA-sensing pathway and TNF signaling pathway (Fig. 6D). Furthermore, we also presented six Hub gene-related signaling pathways (Fig.
6E).
Fig. 6Enrichment analysis of APRGs. A BP entries in GO analysis; (B) CC entries in GO analysis; (C) MF entries in GO analysis; (D) KEGG analysis of APRGs; (E) Sankey plots of six Hub gene-related pathways
Results of co-APRGs screening
In this study, we applied three machine learning methods, RF, SVM and LASSO-logistic regression, aiming to identify the key diagnostic markers of APRGs to reveal the biological significance of PANoptosis in AS. The results showed that a total of 10 key APRGS were screened out in LASSO-logistic regression, including AIM2, TNF, IFNG, CASP8, MEFV, ADAR, NLRP12, ALKBH5, ZBP1 and RIPK1 (Fig. 7A-B); A total of 10 key APRGs were screened out in SVM, including TNF, ALKBH5, CASP8, IFNG, AIM2, IRF1, ADAR, XIAP, TMEM173 and NLRP3 (Fig.
7C); We screened out the top 10 key APRGs in RF, including TNF, ALKBH5, ADAR, DNM1L, AIM2, IFNG, CASP8, NLRC5, HIF1A and MAPK1 (Fig. 7D-E). By integrating the results obtained from three machine learning methods, we ultimately identified six key APRGs, namely AIM2, TNF, IFNG, CASP8, ADAR and ALKBH5 (Fig.
7F). To evaluate the model prediction accuracy of the three machine learning methods, we conducted ROC analysis. The results showed that the ROC values of LASSO, SVM, and RF were 0.765, 0.760, and 0.770 respectively, indicating that machine learning is accurate in screening key APRGs (Fig.
7G).
Fig. 7Machine Learning. (A-B) Screening out the top 10 APRGs in LASSO-logistic regression analysis; (C) Screen out 10 APRGs in the SVM analysis; (D-E) Screen out the top 10 APRGs in RF analysis; (F) A total of six key PDCRGs were screened out by three machine learning analysis methods, including AIM2, TNF, IFNG, CASP8, ADAR and ALKBH5; (G) The ROC curves of the three machine learning methods, LASSO, SVM and RF, have ROC values of 0.765, 0.760 and 0.770 respectively
Results of external validation and pathway analysis of co-APRGs
To determine whether the six key APRGs have significance in predicting the diagnosis of AS, ROC analysis was conducted on AIM2, TNF, IFNG, CASP8, ADAR and ALKBH5 using the external datasets GSE11886 and GSE134290, in order to clarify the accuracy of these key APRGs AS characteristic genes of AS disease. By calculating the area under the ROC curve (AUC), the closer the AUC value is to 1, the better the predictive performance. The results show that we detected the correction of batch effects after the combination of the datasets GSE11886 and GSE134290 through PCA analysis.
The results indicate that the aggregation effect of the samples in the same group is good, while the samples in the disease/control group show a significant separation trend, indicating that the quality of the combined data is good. The results showed that in the GSE11886 dataset, the expressions of these six key APRGs were significantly increased compared with the control group, and the AUC values were all higher than 0.5 (Fig. 8A, D).
Additionally, in the GSE134290 dataset, The expression of APRGs in these pipe fittings was also significantly increased, and the AUC values were all higher than 0.5 (Fig. 8B, D), indicating that the six key APRGs can serve AS diagnostic markers for the occurrence of PANoptosis in AS (Fig. 8C).
To further explore the pathways involved in AIM2, TNF, IFNG, CASP8, ADAR and ALKBH5 in the AS dataset, we conducted GSEA analysis of individual genes for these six key APRGs. The results showed that the ADAR in the key APRGs mainly involved 22 signaling pathways (Fig. 9A), among which the adipocytokine signaling pathway had the highest significance (P = 0.0041); AIM2 mainly involves six signaling pathways (Fig.
9B), among which linoleic acid metabolism has the highest significance (P = 0.005); ALKBH5 is mainly involved in 7 signaling pathways (Fig. 9C), among which amyotrophic lateral sclerosis has the highest significance (P = 0.0101); CASP8 mainly involves four signaling pathways (Fig. 9D), with the pentose phosphate pathway having the highest significance (P-value = 0.0041); TNF is mainly involved in three signaling pathways (Fig.
9E), among which the p53 signaling pathway has the highest significance (P = 0.0016); IFNG mainly involves 7 signaling pathways (Fig. 9F), among which basal transcription factors (P = 0.0241).
Fig. 8ROC analysis. (A-B) The top 10 APRGs were screened out in LASSO-logistic regression analysis; (C) PCA analysis after merging the datasets GSE11886 and GSE134290; (D) The expression of AIM2, TNF, IFNG, CASP8, ADAR and ALKBH5 in two datasets
Fig. 9GSEA analysis of 6 key APRGs. (A) The Adipocytokine signaling pathway was the most significant pathway for ADAR (P-value = 0.0059); (B) Linoleic acid metabolism is the pathway with the highest significance of AIM2 (P-value = 0.005); (C) amyotrophic lateral sclerosis is the pathway with the highest significance of ALKBH5 (P-value = 0.0101); (D) The pentose phosphate pathway is the most significant pathway for CASP8 (P-value = 0.0041); (E) The p53 signaling pathway is the most significant pathway for TNF (P-value = 0.0016); (F) basal transcription factors were the most significant pathway for IFNG (P-value = 0.0241)
Results of APRGs immune infiltration
In this study, the expression matrix of 23 APRGs was thoroughly analyzed using the CBERSORT algorithm, with the aim of precisely quantifying the abundance and spatial distribution characteristics of various immune cells in samples from AS patients. Compared with healthy patient samples, significantly decreased expression of Macrophages M2, Neutrophils, activated memory CD4 + T cells, and Monocytes was observed in AS patient samples, while significantly increased expression of activated Mast cells, activated and resting NK cells, naive CD4 + T cells, and regulatory T cells (Tregs) was demonstrated (Fig. 10A, C).
Additionally, an increase in the proportion of activated memory CD4 + T cells and activated dendritic cells was found in AS patient samples compared to normal healthy samples (Fig. 10B). Furthermore, significant positive correlations were observed between five key APRGs (AIM2, TNF, IFNG, ADAR, and ALKBH5) and Neutrophils, Monocytes, Macrophages M0, regulatory T cells (Tregs), CD8 + T cells, while CASP8 was found to be significantly negatively correlated with these cells (Fig.
10D).
Fig. 10Immune infiltration analysis. A Box plots of 22 types of immune cells in each sample of the AS dataset; (B) A stacked bar chart of the proportion of 22 types of immune cells in each sample of the AS dataset; (C) Heat maps of the expression of 22 types of immune cells in AS patients and normal healthy individuals; (D) Expression of AIM2, TNF, IFNG, CASP8, ADAR and ALKBH5 in 22 types of immune cells
Results of in vivo validation of APRGs in AS
In the bioinformatics analysis section, we identified six key APRGs. To verify the obtained results, we constructed an AS animal model by injecting CFA into the bilateral sacroiliac joint regions of rats and detected the expression levels of inflammatory factors and key APRGs in the spinal tissues of rats. After 4 weeks of modeling, compared with the Sham group, the proportion of resting time during tail suspension and the arthritis score of rats in the AS group were significantly increased (P < 0.01).
In the spinal tissue, it was observed that the Safranine O staining in the nucleus pulposus and cartilage endplate was significantly weakened, along with abnormal bone formation and ligament ossification, accompanied by local congestion and swelling (Fig. 11A-C). In addition, the expressions of PINP and CTX-I in the serum of rats in the AS group were significantly increased compared with those in the Sham group (P < 0.01), indicating that bone formation activity and bone resorption activity were enhanced in AS (Fig.
11D). We speculated that this pathological and behavioral change might be related to the occurrence of PANoptosis. Therefore, we examined the co-expression of AIM2 and ZBP1 (Fig.
11E), and the results showed that the co-expression of AIM2 and ZBP1 in the spinal tissue of rats in the AS group was enhanced (P < 0.01). The HE staining results showed that compared with the Sham group, the muscle tissue of rats in the AS group had uneven muscle fiber sizes, some fibers atrophied and were disordered in arrangement, and inflammatory cell aggregation could be seen around the spine. However, no such changes occurred in the Sham group rats (Fig.
12A-B). Moreover, we also found that the levels of inflammatory factors TNF-α, IL-6, IL-10, and IL-4 in the serum of rats in the AS group were significantly increased (P < 0.01) (Fig. 12D).
It is worth noting that the expression of the PANoptosis marker protein RIPK1 in the muscle tissues of the two groups of rats was detected by IHC. We found that the expression of RIPK1 in the AS model rats was significantly higher than that in the Sham group (P < 0.01) (Fig. 12C).
In addition, we also examined the mRNA and protein expression of six key APRGs in the spinal tissues of rats in each group (Fig. 12E). The results showed that compared with the Sham group rats, the mRNA and protein expressions of AIM2, TNF, IFNG, CASP8 and ALKBH5 in the spinal tissues of the AS group rats were significantly increased, while the mRNA and protein expressions of ADAR were significantly decreased (P < 0.01).
These results confirmed our conjecture. In the CFA-induced AS model rats, we found that the level of PANoptosis was elevated and activated the inflammatory response, causing massive cell death and subsequently leading to the occurrence of AS.
Fig. 11Construction and evaluation of AS rat models. A; The proportion of resting time in the tail suspension test of each group of rats, n = 5, compared with the Sham group, **P < 0.01; (B) Statistics of the arthritis score scale scores of each group of rats, n = 5, compared with the Sham group, **P < 0.01; (C) SO/FG staining of spinal tissues in each group of rats, 1:50 μm; (D) The expressions of PINP and CTX-I in the serum of rats in each group, n = 5, compared with the Sham group, **P < 0.01; (E) The expression of ZBP1 and AIM2 in the spine of rats in each group, n = 3, compared with the Sham group, **P < 0.01
Fig. 12The expression of key PDCRGs in AS model rats. A HE staining of spine, ×200; (B) HE staining of paravertebral muscles, ×200, 1:50 μm; (C) The expression level of RIPK1 in the muscle tissues of the two groups of rats, ×200, n=3, 1:50 μm.
Compared with the Sham group, **P<0.01; (D) The expression levels of inflammatory factors in the serum of the two groups of rats, n=3, compared with the Sham group, **P<0.01; (E) The mRNA expression levels of AIM2, TNF, IFNG, CASP8, ADAR and ALKBH5 in the spinal tissues of the two groups of rats, n=3, compared with the Sham group, **P<0.01; (F) The protein expression levels of AIM2, TNF, IFNG, CASP8, ADAR and ALKBH5 in the spinal tissues of the two groups of rats were n=3. Compared with the Sham group, **P<0.01
Discussion
AS, as a chronic autoimmune inflammatory disease, is characterized by progressive bone erosion and ectopic ossification of the axial skeleton as its core pathological features, significantly affecting the quality of life and social burden of young people. Recent studies have found that PANoptosis An inflammatory programmed cell death pattern integrating pyroptosis, apoptosis and necroptosis - coordinating multiple death effector molecules through PANoptosome complexes (such as ZBP1-PANoptosome, AIM2-PANoptosome). It drives irreversible inflammatory cascades in infections, cancers and autoimmune diseases [21, 22].
In AS, the chronic inflammatory microenvironment can activate PANoptosis, promoting the release of inflammatory factors and tissue damage, and forming a positive feedback loop of “inflammation - cell death - bone destruction”. Therefore, analyzing the regulatory network of PANoptosis-related genes in AS provides a new perspective for breaking through the therapeutic bottleneck.
Through PPI analysis of the AS dataset, we screened out six core genes related to PANoptosis, including AIM2, CASP8, NLRP3, IRF1, TNF, and IFNG, which constitute the core hubs for the regulation of PANoptosis in AS. Melanoma deletion factor 2 (AIM2) contains the PYD domain and the HIN200 DNA binding domain, which can recognize cytoplasmic dsDNA and recruit ASC/caspase-1 to form inflammasomes, driving GSDMD-mediated pyroptosis [23–25]. In AS, the locally released self-DNA of the sacroiliac joint may activate AIM2, promote the secretion of IL-1β and the differentiation of Th17, and aggravate spinal inflammation [26, 27].
CASP8, as a death effector domain (DED) protein, integrates death receptor signals through the FADD adaptor and activates apoptosis executive protein (caspase-3/7). It simultaneously inhibits the RIPK3-MLKL necrotic pathway and promotes the assembly of NLRP3 inflammasomes under specific conditions [28, 29]. The activity of caspase-8 in the joint fluid of AS patients is elevated, suggesting that it expands inflammatory damage by coordinating multiple modes of death [30, 31].
The LRR domain of NoD-like receptor heat protein domain-associated protein 3 (NLRP3) senses mitochondrial reactive oxygen species (ROS) or potassium ion efflux, recruits caspase-1 through ASC, and induces IL-1β maturation and pyroptosis [32–34]. The expression of NLRP3 in monocytes of AS patients is upregulated and positively correlated with disease activity, which may amplify inflammation through the IL-23/IL-17 axis [35]. Interferon regulatory factor 1 (IRF1), as a transcription factor, can induce the expression of PANoptosome components such as ZBP1 and AIM2.
In AS, IRF1 may respond to IFN-γ, promoting PANoptosis in joint macrophages and releasing DAMP molecules (such as HMGB1) to aggravate bone erosion. TNF and IFNG, as key inflammatory cytokines, induce RIPK1-PANoptosome assembly by synergistically activating the JAK-STAT-IRF1 axis. The levels of both in the serum of AS patients are significantly elevated, directly associated with the process of spinal fusion [36–38].
Subsequently, through machine learning analysis, we further screened out six key APRGs, including AIM2, TNF, IFNG, CASP8, ADAR, and ALKBH5. Among them, the newly added RNA modify-related genes (ADAR, ALKBH5) have expanded the epigenetic regulatory dimension of PANoptosis. ADAR (adenosine deaminase acting on RNA) contains the Zα domain and inhibits endogenous Z-RNA accumulation through A-to-I RNA editing, blocking ZBP1-mediated PANoptosome activation [39].
In AS, ADAR dysfunction may lead to abnormal activation of ZBP1, triggering the release of joint synovial cells PANoptosis and IL-18, and promoting chronic inflammation. ALKBH5 (AlkB homolog 5) acts as an m6A demethylase, and its DSBH domain mediates RNA substrate recognition [40–42]. It enhances tumor immune escape by stabilizing PD-L1 mRNA, and may regulate the m6A modification of innate immune genes such AS TRAF6 in AS, affecting the activation of the NF-κB pathway and inflammasome activation [43].
It is worth noting that AIM2, TNF, IFNG, and CASP8 repeatedly appear in PPI and machine learning, confirming their core position in the AS pan-apoptotic network. The inclusion of ADAR and ALKBH5 revealed a new mechanism by which RNA epigenetic modifications affect the course of AS by regulating the expression of ZBP1-PANoptosome (ADAR) and immune genes (ALKBH5).
The main enriched pathways of the above-mentioned APRGs include the NOD-like receptor signaling pathway, the cytoplasmic DNA sensing pathway, the TNF signaling pathway and the p53 signaling pathway. These pathways form cross-pathway collaborative networks in AS and PANoptosis. The NOD-like receptor pathway, with NLRP3 at its core, assembles inflammasomes in response to mitochondrial damage (such as ROS) and activates the caspase-1/GSDMD axis to drive pyroptosis [44].
In AS, HLA-B27 misfolding may induce mitochondrial stress and activate this pathway. The cytoplasmic DNA sensing pathway relies on AIM2 to recognize its own DNA in the joint microenvironment (such as dsDNA released by NETs), promotes the maturation of IL-1β through the ASC-caspase-1 consortium, and collaborates with ZBP1 to form AIm2-panoptosome, inducing PANoptosis [45–47]. In the TNF signaling pathway, TNF-α and IFN-γ activate IRF1 through JAK-STAT1 and upregulate the expression of ZBP1 [48]; Meanwhile, TNF-α induces the ubiquitination state transition of RIPK1 and promotes the assembly of RIPK1-panoptosome (containing RIPK1/RIPK3/caspase-8) [49].
This pathway is directly associated with inflammatory cell infiltration and bone destruction in the sacroiliac joint in AS. The p53 signaling pathway affects the cell cycle and DNA repair by regulating genes such as ALKBH5. ALKBH5-mediated m6A demethylation can stabilize p53 target genes and promote apoptosis [50]; In AS, p53 may amplify inflammarome activation through metabolic reprogramming, forming the “p53-metabolism-PANoptosis” axis.
IL-23 drives CD4⁺ T cells to differentiate into Th17 cells and secrete a large amount of IL-17, thereby activating fibroblast-like synovial cells and macrophages, releasing TNF and IFNG, and forming an early inflammatory microenvironment. As a cytoplasmic DNA sensor, AIM2 can be activated by its own DNA released by inflammation or metabolic stress, assembling into inflammasomes. Through CASP8, it synergically activates CASP1, inducing the maturation of IL-1β and IL-18 and triggering pyroptosis, further amplifying the IL-17-mediated inflammatory cascade.
Meanwhile, CASP8, as a core executor of PANoptosis, integrates signals from death receptors such as TNF receptors and AIM2 inflammasome signals, coordinating the cross-activation of pyroptosis, apoptosis and necroptosis, leading to uncontrollable death of stromal cells and immune cells as well as tissue damage [51]. ADAR regulates the stability and translation efficiency of key gene mRNAs such as IL23R and IL17A through RNA editing and ALKBH5 through m6A demethylation, thereby positively enhancing the activity and persistence of the IL-23/IL-17 signaling pathway (Fig. 13).
Fig. 13The mechanism by which key APRGs are involved in the course of AS through the IL-23/IL-17 axis
Based on the results of GSEA analysis, it was identified that key APRGs (ADAR and AIM2) may drive the pathological process of ankylosing spondylitis through the adipocytokine signaling pathway and linoleic acid metabolism, respectively. The adipocytokine signaling pathway is modulated by ADAR, which is posited to function as an “immunometabolic homeostasis regulator.” Dysregulation of ADAR may alter the signaling balance of pro-/anti-inflammatory adipokines, including leptin and adiponectin and their receptors, through RNA editing. Concurrently, a failure to effectively suppress endogenous nucleic acid-triggered type I interferon responses collectively shapes and amplifies the chronic pro-inflammatory microenvironment required for AS.
Furthermore, cell stress induced by abnormal linoleic acid metabolites leads to self-DNA leakage into the cytoplasm, which specifically activates the AIM2 inflammasome, triggering caspase-1-mediated pyroptosis and the maturation and release of IL-1β and IL-18. More critically, the potent inflammatory signals driven by AIM2, when superimposed upon the inflammatory background shaped by ADAR, transcend the boundaries of individual cell death modalities and synergistically initiate a broader PANoptosis program. This cascade ultimately results in uncontrolled inflammatory propagation, tissue destruction, and aberrant bone remodeling in AS.
The central role of key APRG components, particularly CASP8, in human arthritic diseases is strongly supported by existing clinical research. As evidenced in the provided literature, CASP8 is not only a differentially expressed gene in rheumatoid arthritis (RA) fibroblast-like synoviocytes (FLS) governing their aggressive phenotype [52], but it is also identified as a core target in traditional medicine formulations for RA treatment [53]. Most compellingly, a high-throughput proteomic analysis of psoriatic arthritis (PsA) patient synovial fluid specifically highlighted CASP8 as one of the most significantly altered inflammatory proteins [54].
This independent clinical evidence from a related spondyloarthritis condition strongly corroborates our bioinformatics predictions and underscores the necessity of verifying the expression patterns of our identified APRGs (CASP8, TNF, IFNG) in human AS synovium or joint fluid. Our research results show that in the spinal tissues of CFA-induced AS model rats, the mRNA and protein expressions of AIM2, TNF, IFNG, CASP8 and ALKBH5 were significantly increased, while the mRNA and protein expressions of ADAR were significantly decreased. We also observed an upregulation of inflammatory factors (IL-6, IL-10, IL-4, TNF-α) and the bone resorption marker (CTX-I, PINP).
In addition, we also observed a significant increase in the expression of the pan-apoptotic marker protein RIPK1, as well as an enhanced co-expression of ZBP1 and AIM2. These results indicate that PANoptosis-mediated inflammatory infiltration in AS is strongly associated with the aggravation of the pathological process of AS. However, our research also has limitations.
Our findings are based on transcriptome data analysis and related research. Although we verified the differential expression of key APRGs in animal models, this study lacks direct in vivo or in vitro mechanism experiments, such as gene knockdown or overexpression, to establish a causal relationship between the identified APRGs and the execution of PANoptosis in AS. The inferred interactions in the PANoptosis complex remain a predictive model and require further experimental confirmation.
Conclusion
In summary, we identified six PANoptotic markers in AS through bioinformatics, including AIM2, TNF, IFNG, CASP8, ADAR, ALKBH5, NLRP3 and IRF1. Further analysis through PPI and machine learning identified four key PDCRGs. They are AIM2, TNF, CASP8 and IFNG respectively.
These APRGs are mainly related to the P53 pathway. In the CFA-induced AS model rats, the expression of these APRGs increased, which was significantly associated with the occurrence of AS.
Supplementary Information
Supplementary Material 1.