Whole blood samples were obtained from 6 patients with AS (3 males and 3 females) who met the modified New York classification criteria, had comparable disease duration and inflammatory levels, and a Bath Ankylosing Spondylitis Disease Activity Index (BASDAI) score ≥ 4. Age- and sex-matched healthy controls (HC) (n = 6) were included. All AS patients had no history of other autoimmune or systemic diseases such as rheumatoid arthritis, systemic lupus erythematosus, inflammatory bowel disease, psoriatic arthritis, psoriasis, or major cardiovascular, hepatic, renal, or hematologic disorders. For all subjects, clinical data including age, sex, disease duration, use of biologics, BASDAI, Bath Ankylosing Spondylitis Functional Index (BASFI), C-reactive protein (CRP), and HLA-B27 status were collected. Detailed participant characteristics are summarized in Table 1. Approximately 3~4 mL of whole blood was collected from each individual, and total RNA was extracted for subsequent long-read RNA sequencing using the ONT platform.

Total RNA was extracted from whole blood using the Total RNA Extraction Kit (DP431, TIANGEN, China), and its purity and concentration were assessed by measuring the absorbance ratio at 260/280 nm using an N50 Touch spectrophotometer (IMPLEN, Germany). RNA integrity was verified by 1% agarose gel electrophoresis. For each sample, 500 ng of total RNA was adjusted to 9 μL with nuclease-free water and mixed with 1 μL of 10 μM VNP primer and 1 μL of 10 mM dNTPs. The mixture was incubated at 65 °C for 5 minutes and then immediately snap-cooled. Strand-switching buffer containing 5 × RT buffer, RNaseOUT, nuclease-free water, and a 10 μM strand-switching primer was added, followed by incubation at 42 °C for 2 minutes. Subsequently, 1 μL of Maxima H Minus Reverse Transcriptase was added, and the reaction proceeded at 42 °C for 90 minutes, followed by heat inactivation at 85 °C for 5 minutes and cooling to 4 °C. For cDNA amplification, 20 μL of first-strand product was combined with 2 μL each of PR1 and PR2 primers, 25 μL of LongAmp Taq 2 × Master Mix, and 1 μL of water. The PCR conditions were: 94 °C for 3 minutes; 12 cycles of 94 °C for 15 seconds, 50 °C for 15 seconds, and 65 °C for 2 minutes; and a final extension at 65 °C for 10 minutes. Barcoding was performed using the ONT Native Barcoding Expansion Kit (SQK-NBD114.96) according to the manufacturer’s protocol. Barcoded DNA was purified using 1 × AMPure XP beads, eluted in 20 μL of nuclease-free water, and quantified using a Qubit fluorometer. Finally, 300 fmol of adapter-ligated library was loaded onto a PromethION R10.4.4 flow cell and sequenced on the ONT platform.

Long-read RNA sequencing was performed on AS and HC samples using ONT. Raw off-machine data were obtained in POD5 format containing signal-level information. Basecalling was carried out using Dorado (v0.8.2), which converted POD5 files into FASTQ sequences with associated quality scores. Full-length reads were identified and trimmed using Pychopper (v2.7.9), and low-quality reads (mean Q score <10 or length <50 bp) were filtered with Chopper (v0.7.0). Clean reads were aligned to the Ensembl GRCh38_v45 reference genome using Minimap2 (v2.27-r1193) (17) in spliced alignment mode with parameters -ax splice -uf -k14. Isoform-level transcript detection and correction were performed with the Pinfish pipeline, and the resulting polished transcripts were output in GTF format. To create a comprehensive, non-redundant transcriptome across all samples, individual transcript assemblies were merged using StringTie (v2.1.6) (18). Annotation comparison against reference transcripts was conducted using GffCompare (v0.12.6) with parameters -G and -R, allowing classification of novel transcripts and genes (19). Transcript class codes [u, x, i, j, o] were considered novel, while [=, c] indicated known isoforms. All transcripts were merged and retained for downstream quantification and analysis.

To quantify expression at both transcript and gene levels, we used Salmon (v1.7.0) (20) in quasi-mapping mode. Transcripts per million (TPM) values were calculated for each transcript and gene across all samples to provide normalized abundance estimates for downstream differential expression analysis.

Differential expression analysis at both gene and transcript levels was conducted using DESeq2 (v1.30.1) (21), applying thresholds of fold change ≥ 2 and false discovery rate (FDR) ≤ 0.05 to identify differentially expressed genes (DEGs) and differentially expressed transcripts (DETs). Differential transcript usage (DTU) analysis was performed using DEXSeq (22), as implemented in the IsoformSwitchAnalyzeR package (v2.2.0) (23). Transcript annotations generated by GffCompare (v0.12.6) were used as input. For each isoform, the difference in isoform fraction (dIF) between AS and HC groups and its FDR-adjusted p-value were calculated. Isoforms with adjusted p < 0.05 were considered to exhibit significant DTU regardless of dIF value. At the gene level, the lowest adjusted p-value among isoforms was used to assign gene-level significance, and genes with at least one significantly altered isoform were defined as DTU-positive.

Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (24) were annotated and enriched using KOBAS (v2.0). Statistical significance was assessed by hypergeometric testing, and multiple testing correction was performed using the Benjamini-Hochberg procedure to control the false discovery rate (FDR).

Alternative splicing events were identified and quantified using SUPPA2 (v2.3) (25). Splicing levels were measured as percent spliced-in (PSI) values for each event. Differential splicing between AS and HC groups was assessed by calculating the difference in PSI (ΔPSI) and performing independent t-tests. AS events with ΔPSI ≥ 0.05 and p-value ≤ 0.05 were considered significantly differentially spliced and retained for further analysis.

A curated list of genes involved in 12 distinct types of PCD was obtained from Supplementary Table 1. These PCD-related genes were cross-referenced with the differential expression, transcript usage, and splicing datasets to identify potential regulatory patterns in AS.

RNA-seq datasets (GSE205812 and GSE181364) were downloaded from GEO. The raw reads were trimmed of low-quality bases using a FASTX-Toolkit (v.0.0.13; http://hannonlab.cshl.edu/fastx_toolkit/). Then the clean reads were evaluated using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc).

The retrieved clean reads were aligned onto the human genome version GRCh38_v45 using HISAT2(v.2.2.1). Uniquely mapped reads were selected for further analysis, and the number of reads located on each gene was calculated. The expression levels of genes were evaluated using FPKM (fragments per kilobase of exon per million fragments mapped). DEseq2 (v.1.30.1) software was used to perform differential gene expression analysis using the reads count file (21). DEseq2 was also used to analyze the differential expression between two or more samples and thus determine whether a gene was differentially expressed by calculating the fold change (FC) and false discovery rate (FDR), FC ≥ 2 or ≤ 0.5, FDR ≤ 0.05.

The gene expression profiles of microarray datasets GSE25101and GSE73754 were downloaded from GEO. Differential gene expression was performed using limma (26).

We defined genes showing AS-regulated expression at both gene and transcript levels, as well as AS-regulated alternative splicing as core AS pathogenesis genes. The differential expression of some of these genes was further validated in public datasets.

All statistical analyses and visualizations, including pattern diagrams and stacked bar charts, were performed in R (v4.2.3) using RStudio. Data are presented as mean ± standard error of the mean (SEM). Comparisons between two groups were conducted using Student’s t-test, with p-values < 0.05 considered statistically significant unless otherwise specified.

Long-read sequencing using ONT was performed on blood samples from 6 AS patients and 6 matched health control (HC) to investigate transcriptomic alterations. The data-analysis strategy in this study is illustrated (Supplementary Figure 1A). PCA based on transcript-level TPM values clearly separated AS from HC, indicating distinct global expression profiles (Figure 1A). A total of 81,940 transcripts corresponding to 25,343 genes were identified, including 1,902 novel transcripts from 543 novel genes, a slightly larger total number of transcripts (both known and novel) was detected in AS samples compared to HC samples (Figure 1B). Many genes exhibited high isoform diversity, with a substantial proportion expressing more than 10 isoforms (Figure 1C). Among the novel transcripts, GffCompare classification showed that most were assigned to class codes representing novel isoforms (j, i, o), antisense transcripts (x), or intergenic loci (u), suggesting the presence of previously unannotated transcriptional activity (Supplementary Figures 2A–C). The majority of quantified transcripts were protein-coding, and novel transcripts were generally shorter than known ones, as shown by the log-transformed length comparisons (Figures 1D, E). We performed functional enrichment analysis with genes associated with novel transcripts in AS and HC groups, demonstrating that the significantly enriched GO terms in both groups were similarly related to AS pathogenesis including inflammatory response, innate immune response, and actin cytoskeloton organization (Figure 1F). GO terms and KEGG pathways enriched by all transcripts detected in this study were also shown. AS related KEGG pathways included osteoclast differentiation, autophagy, and antigen processing (Supplementary Figures 2C, D). These findings highlight the extensive transcriptomic complexity in AS and suggest a functional relevance of novel isoforms in disease-related pathways.

We further characterized transcript-level expression differences between AS and HC, identifying 1,812 DETs, including 1,286 upregulated and 526 downregulated transcripts (Figure 2A). The broad distribution of DETs underscores extensive isoform-level alterations in AS. GO enrichment analysis revealed that upregulated DETs were significantly associated with transcriptional regulation via RNA polymerase II, apoptotic signaling, inflammation, autophagy, and innate immune activation (Figure 2B). In contrast, downregulated DETs were enriched in biological processes including humoral immune response, mRNA splicing and processing, antiviral defense, and regulation of cytoplasmic calcium ion levels (Figure 2C). These findings suggest that AS involves coordinated remodeling at the transcript level, characterized by activation of pro-inflammatory and cell death pathways and suppression of immune regulatory and RNA processing functions.

To investigate gene-level alterations in AS, we performed PCA based on gene-level TPM values, which revealed distinct separation between AS and HC samples, indicating overall differences in gene expression profiles (Figure 3A). We subsequently identified 1,088 DEGs between the two groups, including 724 upregulated and 364 downregulated genes, with a notable shift toward upregulation in AS (Figure 3B). GO enrichment analysis showed that upregulated DEGs were predominantly involved in transcriptional activation via RNA polymerase II, apoptosis, autophagy, and angiogenesis, suggesting increased cellular activity and tissue remodeling in AS (Figure 3C). In contrast, downregulated genes were enriched in chemokine signaling, calcium-mediated pathways, innate immune suppression, and host defense processes (Figure 3D). Together, these findings indicate a global shift in gene expression that promotes inflammation, stress response, and immune modulation in AS.

To explore the involvement of PCD pathways in AS, we integrated DETs and DEGs, identifying 508 overlapping genes between the two datasets (Figure 4A). GO enrichment analysis of these overlapping genes revealed their significant association with biological processes including apoptosis, calcium signaling regulation, innate and adaptive immune responses, and cytoskeletal organization (Figure 4B). To further examine their relevance to PCD, we cross-referenced the 508 overlapping genes with curated databases covering 12 types of cell death. A substantial number were associated with apoptosis (29 genes), autophagy (19 genes), lysosome-dependent cell death (7 genes), ferroptosis (2 genes), and necroptosis (1 gene) (Figure 4C), indicating broad transcriptional involvement of PCD-related mechanisms in AS. Among them, FAS (Fas cell surface death receptor), FASLG (Fas ligand), and ULK1 (unc-51 like autophagy activating kinase 1), key regulators of apoptosis and autophagy, exhibited consistent and significant differential expression at both gene and transcript levels. Specifically, FAS and ULK1 were upregulated, while FASLG was downregulated in AS compared to HC (Figures 4D, E). These findings suggest that dysregulated expression of core PCD components may contribute to immune and inflammatory imbalances in AS.

To further investigate transcript-level regulatory dynamics in AS, we analyzed DTU between AS and HC using DEXSeq implemented in IsoformSwitchAnalyzeR. A total of 50 transcripts exhibited significant differences in isoform usage, independent of overall gene expression changes (Figure 5A). GO enrichment analysis revealed that these DTU events were predominantly associated with innate immune response and regulation of RNA polymerase II transcription (Figure 5B). Integration of DTU with DET and DEG datasets identified 4 overlapping genes, including LINGO3 (leucine rich repeat and Ig domain containing 3) and SAMD9 (sterile alpha motif domain containing 9) (Figure 5C), both of which displayed distinct isoform usage alterations between AS and HC (Figures 5D, E). In LINGO3, although total gene expression was decreased in AS, the dominant isoform was preferentially used, while the alternative isoform showed reduced usage. In SAMD9, gene and isoform expression levels were relatively similar between groups, but the two isoforms exhibited opposite usage patterns between AS and HC, indicating a pronounced switch in isoform preference. These isoform changes were also associated with protein domain differences, suggesting potential functional divergence. Notably, these usage alterations occurred without proportional changes in isoform expression, reinforcing their regulatory specificity. Together, these findings support isoform-level regulation as an additional layer contributing to transcriptomic remodeling in AS.

Building on the observed isoform-level changes, we next explored alternative splicing regulation to better understand the mechanisms underlying transcript diversity in AS. A total of 7 types of alternative splicing were identified, including exon skipping (SE), retained introns (RI), mutually exclusive exons (MX), alternative 5′ and 3′ splice site usage (A5, A3), as well as alternative first exons (AF) and alternative last exons (AL), with AF being the most prevalent AS type, followed by SE (Figure 6A). Distributions of PSI values and corresponding ΔPSI metrics showed global variability between AS and HC (Figures 6B, C). Following statistical filtering, 304 significantly altered splicing events were retained, which clearly separated AS from HC samples in the PSI-based splicing heatmaps (Figure 6D). GO enrichment analysis of alternatively spliced genes indicated functional involvement in mRNA splicing, transcriptional regulation, cellular stress responses, autophagy, and apoptotic signaling (Figure 6E), suggesting that AS-related splicing alterations may affect key immune and cell death pathways. Representative examples of disease-relevant splicing events included reduced exon inclusion in FCGR2B (Fc gamma receptor IIb), TLR2 (Toll-like receptor 2), and STAT5B (Signal transducer and activator of transcription 5B) in AS compared to HC (Figure 6F). These splicing changes may alter isoform structure and function, potentially impacting immune regulation and inflammatory signaling in AS.

To investigate whether alternative splicing contributes to the regulation of PCD-related transcript expression in AS, we integrated DETs with genes exhibiting significant AS events and identified 111 overlapping genes affected at both the splicing and transcript levels (Figure 7A). GO enrichment analysis of these genes revealed strong associations with apoptotic signaling, innate immune responses, mRNA processing, DNA damage repair, and transcriptional regulation (Figure 7B), suggesting their functional relevance in AS pathology. Cross-referencing with curated PCD gene sets showed that these dual-regulated genes included 9 related to apoptosis, 5 to autophagy, and others involved in lysosome-dependent and pyroptotic pathways (Figure 7C). These findings indicate that alternative splicing may cooperate with transcript-level expression changes to modulate cell death-associated genes, contributing to immune dysregulation and disease progression in AS.

To pinpoint key genes under coordinated transcriptional and post-transcriptional regulation in AS, we performed an integrative analysis combining DEGs, DETs and alternative splicing events. This approach identified 26 genes exhibiting concurrent dysregulation at the levels of gene expression, transcript expression and splicing regulation, highlighting candidates potentially central to AS pathogenesis (Figure 8A). Among them, NAMPT (Nicotinamide Phosphoribosyltransferase), GATA2 (GATA Binding Protein 2), and DDIT3 (DNA Damage Inducible Transcript 3) were notable for their consistent and significant dysregulation across all three regulatory layers. These genes also showed significant differences in PSI values (Figure 8B), suggesting altered splicing regulation. Gene-level expression analysis revealed upregulation of NAMPT and DDIT3 and downregulation of GATA2 in AS (Figure 8C), while transcript-level data confirmed isoform-specific expression changes (Figure 8D). Functionally, each of these genes is intimately linked to PCD pathways. These results identify a subset of genes dysregulated across multiple transcriptomic layers in AS, many of which are functionally linked to immune regulation and programmed cell death pathways.

We then downloaded four other available PBMC transcriptomic datasets to investigate the AS-deregulated core gene expression. We found that the expression of both NAMPT and DDIT3 was increased in one RNA-seq dataset (GSE205812) and one microarray dataset (GSE25101), although the increase did not reach statistical significance (Figures 8E, F), while the down-regulated expression of GATA2 was not evident (data not shown). This finding supports the potential importance of NAMPT and DDIT3 in AS pathobiology.

We noted that differential was marginal in microarray datasets on AS (GSE25101and GSE73754), with only a few genes identified as differentially expressed (data not shown). We performed GO analysis of the overlapped genes between the 508 AS-regulated DEGs/DETs identified in this study (Figure 4A) and AS-regulated DEGs identified from two previously published RNA-seq datasets (GSE205812 and GSE181364) (Figure 8G). These overlapping genes were enriched in transcriptional regulation and the apoptotic process (Figure 8H).