Abstract
Ankylosing spondylitis (AS) is a chronic immune-mediated arthritis marked by persistent inflammation and progressive structural damage. Although dysregulation of programmed cell death (PCD) is increasingly recognized in AS pathogenesis, the full spectrum of transcript-level regulation remains unclear. Here, we employed Oxford Nanopore Technologies (ONT) long-read RNA sequencing to comprehensively profile peripheral blood transcriptomes from six AS patients and six matched healthy controls. Our analysis identified 1,088 differentially expressed genes (DEGs) and 1,812 differentially expressed transcripts (DETs), with upregulated transcripts enriched in apoptosis, autophagy, and transcriptional regulation. We further detected 50 transcripts with significant differential usage and 304 alternative splicing events affecting immune- and PCD-related genes, including FCGR2B, TLR2, and STAT5B. Integrative multilayered analysis revealed 26 core genes, such as NAMPT, GATA2, and DDIT3, showing consistent dysregulation at gene, isoform, and splicing levels, highlighting convergent regulatory networks underlying immune imbalance and cell death in AS. These findings provide the first isoform-resolved transcriptomic landscape of PCD regulation in AS, which unveils extensive regulatory complexity and nominates a set of core genes for future mechanistic and therapeutic exploration.
Introduction
Ankylosing spondylitis (AS) is a chronic immune-mediated arthritis primarily affecting the sacroiliac joints and axial skeleton, often resulting in persistent pain, structural damage, and functional impairment (1). The global prevalence of AS varies by population, with higher rates observed in individuals carrying human leukocyte antigen B27 (HLA-B27), and a significant disease burden remains due to early onset and lifelong morbidity (2). Despite the advent of biologic therapies, including tumor necrosis factor-alpha (TNF-α) and interleukin-17A (IL-17A) inhibitors, a considerable proportion of patients still experience suboptimal responses or relapse (3, 4). This underscores the need to better understand upstream mechanisms potentially contributing to AS pathogenesis, such as programmed cell death (PCD).
Accumulating evidence implicates dysregulation of programmed cell death (PCD) pathways, including apoptosis, autophagy, necroptosis, and pyroptosis, as central to the pathogenesis and immune dysregulation of AS (5–9). However, transcriptomic investigations in AS currently have predominantly relied on gene-level analyses with short-read sequencing platforms. Such approaches are unable to powerfully resolve isoform diversity or detect alternative splicing events, which represent essential regulatory layers modulating PCD pathways and immune cell fate (10–12). For example, alternative splicing of FAS (CD95) can yield both pro-apoptotic and dominant-negative isoforms, fundamentally altering cellular outcomes (13), yet the scope and impact of such regulatory plasticity in AS remain undefined.
Recent advances in long-read sequencing, such as those offered by Oxford Nanopore Technologies (ONT), now enable accurate profiling of full-length transcripts, thereby capturing isoform usage and splicing complexity with unprecedented resolution (14). These platforms have already revealed novel regulatory shifts in various immune-mediated diseases (15, 16). We hypothesize that AS pathogenesis involves multi-layered transcriptomic disturbances, including differential gene expression, alternative splicing, and isoform switching, that converge to disrupt PCD networks. Integrative analysis of these transcriptomic layers is likely to identify novel molecular targets and mechanistic insights.
In this study, we apply ONT-based long-read RNA sequencing to systematically profile peripheral blood transcriptomes from AS patients and matched healthy controls (HC). By combining analyses of gene expression, isoform usage, and alternative splicing, we generate the first comprehensive landscape of transcriptomic regulation of PCD in AS and identify core genes operating across both transcriptional and post-transcriptional regulatory layers.
Methods
Sample collection
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.
Table 1
| Characteristic | AS | HC | p value |
|---|---|---|---|
| Number | 6 | 6 | |
| Age, years, mean (SD) | 32.5 (2.8) | 31.7 (3.1) | 0.68 |
| Male, n (%) | 3 (50%) | 3 (50%) | > 0.99 |
| BASDAI, mean (SD) | 5.2 (0.6) | — | — |
| BASFI, mean (SD) | 4.3 (0.5) | — | — |
| CRP (mg/L), mean (SD) | 13.2 (2.9) | 2.1 (0.4) | < 0.001 |
| HLA-B27, n (%) | 6 (100%) | 0 (0%) | < 0.0001 |
| Comorbidities | None | None | — |
Clinical characteristics of AS patients and healthy controls.
RNA extraction, library construction, and ONT sequencing
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 sequencing data processing and transcript identification
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.
Expression quantification
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 and transcript usage 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.
Pathway enrichment analysis
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).
Differential alternative splicing analysis
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.
Programmed cell death 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 analysis
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.
Microarray data analysis
The gene expression profiles of microarray datasets GSE25101and GSE73754 were downloaded from GEO. Differential gene expression was performed using limma (26).
Identification of core AS-related genes
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.
Statistical analysis
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.
Results
Transcriptome landscape of AS revealed by long-read sequencing
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.
Figure 1
Differential transcript expression landscape in AS
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.
Figure 2
Differential gene expression profile in AS
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.
Figure 3
Differential expression of programmed cell death-related genes and transcripts 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.
Figure 4
Differential transcript usage 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.
Figure 5
Alternative splicing landscape 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.
Figure 6
Integration of alternative splicing and transcript expression reveals programmed cell death regulation 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.
Figure 7
Multilayer regulatory convergence identifies core genes associated with AS pathogenesis
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.
Figure 8
The up-regulated expression of NAMPT and DDIT3 in AS is confirmed by additional datasets
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).
Discussion
Our integrated long-read transcriptomic analysis in AS systematically revealed a multilayered landscape of transcriptomic dysregulation, encompassing differential gene and transcript expression, isoform usage, and alternative splicing, particularly within PCD pathways. This represents the first study in AS to comprehensively dissect transcriptional and post-transcriptional regulatory alterations at isoform resolution, and despite the exploratory nature of our cohort (n=6 per group), it provides a novel, high−resolution molecular map that highlights previously inaccessible regulatory layers in AS. We acknowledge that the sample size, though typical for initial long−read sequencing studies due to cost and analytical complexity, limits statistical power and generalizability. Nevertheless, the consistent separation of AS and HC groups in global expression and splicing profiles, together with stringent statistical thresholds and external validation of key findings (see below), supports the robustness of the reported transcriptomic landscape.
In our cohort, we identified 1,088 DEGs and 1,812 DETs distinguishing AS from HC, with a strong bias towards upregulation. Functional enrichment highlighted apoptosis, autophagy, angiogenesis, and transcriptional regulation among upregulated genes and transcripts, whereas downregulated ones were linked to chemokine signaling and mRNA splicing. This pattern reflects a transcriptional shift toward activation of cell death and cellular stress pathways, together with suppression of certain immune regulatory processes, which may contribute to chronic inflammation and tissue remodeling in AS (27, 28). Importantly, many DETs represented isoform - specific changes not captured by conventional gene - level analysis, highlighting the unique capacity of long-read sequencing to uncover previously hidden layers of transcriptomic complexity in autoimmune disease (29).
A unique advantage of our long-read approach is the identification of 50 transcripts with significant DTU between AS and HC. This highlights a layer of regulatory divergence beyond overall gene expression. For example, we observed pronounced shifts in isoform fractions of LINGO3 and SAMD9 in AS. LINGO3 is implicated in mucosal inflammation, whereas SAMD9 is involved in antiviral defense and interferon signaling (30, 31). These isoform changes may influence domain composition and downstream signaling. Previous research in autoimmune diseases has shown that isoform switching can act as a molecular switch in immune activation and tolerance, underscoring the pathogenic relevance of DTU (32). Together, these findings suggest that transcript usage modulation is an important but underappreciated mechanism underlying immune dysregulation in AS.
Integrative analysis across gene, transcript, and splicing layers pinpointed 26 core genes with consistent multi-level dysregulation. Among these, NAMPT, GATA2, and DDIT3 emerge as focal regulators within PCD networks. emerged as focal regulators within PCD networks. These genes exhibited convergent dysregulation at expression, isoform, and splicing levels, suggesting that multi−layered regulatory convergence is a hallmark of PCD and immune dysfunction in AS. NAMPT is a key NAD+ biosynthetic enzyme involved in energy metabolism, inflammation, and neutrophil extracellular trap formation, and has been targeted in rheumatoid arthritis models (33). GATA2 is essential for hematopoietic stem cell maintenance and immune lineage specification. In addition to its role in natural killer (NK) cell development. GATA2 also regulates transcriptional programs involved in immunological tolerance and leukemogenesis, with mutations linked to familial myelodysplastic syndromes and increased infection susceptibility, and its dysregulation may affect immune−cell development and function in AS (34, 35). DDIT3 (also known as CHOP) is a canonical ER stress sensor promoting apoptosis and ferroptosis under oxidative and metabolic stress (7, 36). The observation that these genes are dysregulated at both expression and splicing levels suggests that multi-layered regulatory convergence is a hallmark of PCD and immune dysfunction in AS.
Notably, the up-regulation of NAMPT and DDIT3 was also observed in two independent peripheral blood transcriptomic datasets (GSE205812, RNA-seq; GSE25101, microarray) (Figures 8E, F), reinforcing their potential relevance in AS. Furthermore, comparison of our DEG/DET set with DEGs from two publicly available RNA-seq studies (GSE205812, GSE181364) revealed a significant overlap enriched in apoptotic and transcriptional regulatory processes (Figures 8G, H), indicating that PCD-related transcriptional remodeling is a recurrent feature across AS blood transcriptomic studies (37). For instance, JUNB—a shared DEG in our and two prior datasets—has been implicated in AS pathogenesis through single-cell multi-omics studies (38).
Although our data are derived from peripheral blood, a readily accessible but indirect compartment, the observed transcriptomic shifts likely reflect systemic immune dysregulation that may contribute to or mirror processes in affected tissues. Circulating immune cells can serve as both mediators and reporters of inflammation, and shared pathways between blood and synovial/entheseal tissue have been documented in spondyloarthritis (39). Thus, the isoform-level alterations identified here provide a set of testable hypotheses for future investigations in target tissues or appropriate cellular models.
Furthermore, our analysis revealed that diverse PCD mechanisms are transcriptionally intertwined in AS pathogenesis. Prior research has largely focused on necroptosis and autophagy in synovial inflammation, but our data show co-regulation across multiple PCD subtypes, including lysosome-mediated cell death and pyroptosis. Recent studies have shown that these forms of cell death interact dynamically with inflammation, affecting IL-1β release, macrophage polarization, and tissue damage (6, 40, 41). The link between PCD dysregulation and HLA-B27 misfolding pathways may further connect our transcriptomic observations with known molecular mechanisms of AS. Importantly, recent advances in targeting cell death regulators, such as inhibitors of receptor-interacting serine/threonine-protein kinase 1 (RIPK1), ferroptosis modulators, and autophagy blockers, are beginning to show promise in preclinical rheumatic disease models, opening potential translational pathways from transcriptomic findings to therapeutic innovations (42–44). In particular, NAMPT is a druggable enzyme with existing inhibitors explored in other inflammatory conditions (45) while DDIT3 and its upstream ER−stress pathway represent emerging therapeutic targets (46). Although GATA2 deficiency is classically associated with immunodeficiency and myeloid neoplasia, recent studies highlight its role as a master regulator of endothelial function through miRNA networks such as miR−126 and miR−221, which are critical for vascular homeostasis and angiogenesis. Thus, GATA2 downregulation in AS may contribute to vascular inflammation or impaired endothelial repair —processes relevant to AS-related cardiovascular comorbidities—and could represent a novel biomarker or therapeutic target for modulating immune−vascular crosstalk in this disease (47). The convergent dysregulation of these genes across multiple regulatory layers nominates them as high−priority candidates for further validation as potential biomarkers or therapeutic targets in AS.
Despite these advances, this study has limitations. The relatively small sample size limits statistical power and may introduce cohort-specific bias. Since HLA-B27 itself could drive broad transcriptional and splicing changes in blood, without including HLA-B27 positive healthy controls in this study, we cannot exclude the possibility that some of the disease-related changes in transcript expression and alternative splicing are resulted from HLA-B27. Although our analysis revealed compelling transcript-level signatures, it remains descriptive at the transcriptomic level. Functional validation through protein-level assays, isoform−specific manipulations in immune cells, and investigations in disease−relevant tissues are essential next steps to establish causality and clarify the physiological impact of the identified isoform switches and splicing events. Isoform-specific knockdown or overexpression experiments in immune cells or fibroblasts could help clarify causal relationships.
Conclusion
In conclusion, our long-read transcriptomic profiling of AS provides novel insights into the multilayered regulation of gene activity, especially as it pertains to PCD-related pathways. Beyond revealing differential gene and transcript expression, the study uniquely captures isoform-specific switching and splicing changes that reshape immune signaling architecture in AS. Among the key findings, NAMPT, GATA2, and DDIT3 emerged as candidate regulators situated at the intersection of inflammation and cell death. While the exploratory cohort size warrants caution in over-interpretation, the integration of multi-layer analyses, external dataset validation, and the generation of specific, testable hypotheses substantially strengthen the biological plausibility of these findings. These results not only enhance our mechanistic understanding of AS pathogenesis but also nominate high-priority targets for future validation and therapeutic exploration.
Statements
Data availability statement
The datasets generated and/or analysed during the current study are available in Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE299639).
Ethics statement
The studies involving humans were approved by the Ethics Committee of Henan Provincial People’s Hospital (AF/SC-07/05.0). All participants gave written informed consent prior to enrollment. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
XC: Conceptualization, Data curation, Formal analysis, Methodology, Software, Writing – original draft, Writing – review & editing. PL: Data curation, Methodology, Supervision, Writing – original draft, Writing – review & editing. HP: Software, Writing – original draft, Writing – review & editing. QC: Writing – original draft, Writing – review & editing. LS: Supervision, Writing – original draft, Writing – review & editing. DC: Methodology, Validation, Writing – original draft, Writing – review & editing. TC: Funding acquisition, Supervision, Writing – original draft, Writing – review & editing. YC:.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This work was funded by the National Natural Science Foundation of China (NSFC) (No. 82402529, awarded to Yanwei Cheng), and by the Henan Province Medical Science and Technology Co-construction Projects (LHGJ20220016 to Yanwei Cheng; LHGJ20240049 to Xue Cao).
Acknowledgments
We sincerely thank our team members for their support throughout the research process. We also acknowledge the helpful feedback from peer reviewers and editors, which contributed to improving the quality of this manuscript. In addition, I wish to thank Ruixing Biotechnology Co., Ltd (Wuhan) for kindly providing excellent technical assistance.
Conflict of interest
Author DC was employed by the company of Wuhan Ruixing Biotechnology Co., Ltd.
The remaining author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2026.1640271/full#supplementary-material
Supplementary Figure 1Schematic diagram of the experimental and data-analysis procedures in this study.
Supplementary Figure 2Long-read sequencing reveals transcript expression. (A) Schematic diagram of the experimental and data-analysis procedures in this study. (B) Transcript classification codes based on their relationship to reference transcripts, as generated by GffCompare. (C) The bar plots showing the classification of transcripts with quantification. (C, D) The top 10 most enriched GO terms (biological process) and KEGG pathway were illustrated for genes of new transcript.
References
1
GolderVSchachnaL. Ankylosing spondylitis: an update. Aust Fam Physician. (2013) 42:780–4.
2
DeanLEet al. Global prevalence of ankylosing spondylitis. Rheumatol (Oxford). (2014) 53:650–7. doi: 10.1093/rheumatology/ket387
3
DongQet al. IL-17A and TNF-α inhibitors induce multiple molecular changes in psoriasis. Front Immunol. (2022) 13:1015182. doi: 10.3389/fimmu.2022.1015182
4
How Shing KoyEet al. Immunomodulation with IL-17 and TNF-α in spondyloarthritis: focus on the eye and the central nervous system. Ther Adv Musculoskelet Dis. (2021) 13:1759720x211025894. doi: 10.1177/1759720X211025894
5
QianSet al. Programmed cell death: molecular mechanisms, biological functions, diseases, and therapeutic targets. MedComm (2020). (2024) 5:e70024.
6
ZhaoJet al. Apoptosis, autophagy, NETosis, necroptosis, and pyroptosis mediated programmed cell death as targets for innovative therapy in rheumatoid arthritis. Front Immunol. (2021) 12:809806. doi: 10.3389/fimmu.2021.809806
7
LiQet al. Role of ferroptosis-associated genes in ankylosing spondylitis and immune cell infiltration. Front Genet. (2022) 13:948290. doi: 10.3389/fgene.2022.948290
8
MistryPKaplanMJ. Cell death in the pathogenesis of systemic lupus erythematosus and lupus nephritis. Clin Immunol. (2017) 185:59–73. doi: 10.1016/j.clim.2016.08.010
9
YangYet al. Programmed cell death and its role in inflammation. Mil Med Res. (2015) 2:12. doi: 10.1186/s40779-015-0039-0
10
JiangZHWuJY. Alternative splicing and programmed cell death. Proc Soc Exp Biol Med. (1999) 220:64–72.
11
LiuQFangLWuC. Alternative splicing and isoforms: from mechanisms to diseases. Genes (Basel). (2022) 13. doi: 10.3390/genes13030401
12
UleJBlencoweBJ. Alternative splicing regulatory networks: functions, mechanisms, and evolution. Mol Cell. (2019) 76:329–45. doi: 10.1016/j.molcel.2019.09.017
13
TejedorJRPapasaikasPValcárcelJ. Genome-wide identification of Fas/CD95 alternative splicing regulators reveals links with iron homeostasis. Mol Cell. (2015) 57:23–38. doi: 10.1016/j.molcel.2014.10.029
14
AmentIHet al. Long-read RNA sequencing: A transformative technology for exploring transcriptome complexity in human diseases. Mol Ther. (2025) 33:883–94. doi: 10.1016/j.ymthe.2024.11.025
15
SunXet al. Nanopore sequencing and its clinical applications. Methods Mol Biol. (2020) 2204:13–32.
16
LinBHuiJMaoH. Nanopore technology and its applications in gene sequencing. Biosensors (Basel). (2021) 11. doi: 10.3390/bios11070214
17
LiH. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. (2018) 34:3094–100. doi: 10.1093/bioinformatics/bty191
18
PerteaMet al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. (2015) 33:290–5. doi: 10.1038/nbt.3122
19
PerteaGPerteaM. GFF utilities: GffRead and GffCompare. F1000Res. (2020) 9. doi: 10.12688/f1000research.23297.2
20
PatroRet al. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. (2017) 14:417–9. doi: 10.1038/nmeth.4197
21
LoveMIHuberWAndersS. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8
22
LiYet al. RNA-Seq analysis of differential splice junction usage and intron retentions by DEXSeq. PLoS One. (2015) 10:e0136653. doi: 10.1371/journal.pone.0136653
23
Vitting-SeerupKSandelinA. IsoformSwitchAnalyzeR: analysis of changes in genome-wide patterns of alternative splicing and its functional consequences. Bioinformatics. (2019) 35:4469–71. doi: 10.1093/bioinformatics/btz247
24
KanehisaMet al. KEGG: biological systems database as a model of the real world. Nucleic Acids Res. (2025) 53:D672–d677. doi: 10.1093/nar/gkae909
25
TrincadoJLet al. SUPPA2: fast, accurate, and uncertainty-aware differential splicing analysis across multiple conditions. Genome Biol. (2018) 19:40. doi: 10.1186/s13059-018-1417-1
26
RitchieMEet al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47–7. doi: 10.1093/nar/gkv007
27
MauroDet al. Ankylosing spondylitis: an autoimmune or autoinflammatory disease? Nat Rev Rheumatol. (2021) 17:387–404.
28
SieperJPoddubnyyD. Axial spondyloarthritis. Lancet. (2017) 390:73–84. doi: 10.1016/S0140-6736(16)31591-4
29
LiRet al. Direct full-length RNA sequencing reveals unexpected transcriptome complexity during Caenorhabditis elegans development. Genome Res. (2020) 30:287–98. doi: 10.1101/gr.251512.119
30
ZulloKMet al. LINGO3 regulates mucosal tissue regeneration and promotes TFF2 dependent recovery from colitis. Scand J Gastroenterol. (2021) 56:791–805. doi: 10.1080/00365521.2021.1917650
31
LiuJMcFaddenG. SAMD9 is an innate antiviral host factor with stress response properties that can be antagonized by poxviruses. J Virol. (2015) 89:1925–31. doi: 10.1128/JVI.02262-14
32
KatsoulaGet al. A molecular map of long non-coding RNA expression, isoform switching and alternative splicing in osteoarthritis. Hum Mol Genet. (2022) 31:2090–105. doi: 10.1093/hmg/ddac017
33
BussoNet al. Pharmacological inhibition of nicotinamide phosphoribosyltransferase/visfatin enzymatic activity identifies a new inflammatory pathway linked to NAD. PLoS One. (2008) 3:e2267. doi: 10.1371/journal.pone.0002267
34
CollinMDickinsonRBigleyV. Haematopoietic and immune defects associated with GATA2 mutation. Br J Haematol. (2015) 169:173–87. doi: 10.1111/bjh.13317
35
MaceEMet al. Mutations in GATA2 cause human NK cell deficiency with specific loss of the CD56(bright) subset. Blood. (2013) 121:2669–77. doi: 10.1182/blood-2012-09-453969
36
MarciniakSJet al. CHOP induces death by promoting protein synthesis and oxidation in the stressed endoplasmic reticulum. Genes Dev. (2004) 18:3066–77. doi: 10.1101/gad.1250704
37
LiSet al. Sex dimorphism of IL-17-secreting peripheral blood mononuclear cells in ankylosing spondylitis based on bioinformatics analysis and machine learning. BMC Musculoskelet Disord. (2024) 25:490. doi: 10.1186/s12891-024-07589-6
38
XuHet al. Integrative single-cell RNA-Seq and ATAC-seq analysis of peripheral mononuclear cells in patients with ankylosing spondylitis. Front Immunol. (2021) 12:760381. doi: 10.3389/fimmu.2021.760381
39
SimoneDet al. Single cell analysis of spondyloarthritis regulatory T cells identifies distinct synovial gene expression patterns and clonal fates. Commun Biol. (2021) 4:1395. doi: 10.1038/s42003-021-02931-3
40
BenderECTareqHSSuggsLJ. Inflammation: a matter of immune cell life and death. NPJ Biomed Innov. (2025) 2:7. doi: 10.1038/s44385-025-00010-4
41
YangGet al. Regulation of alveolar macrophage death in pulmonary fibrosis: a review. Apoptosis. (2023) 28:1505–19. doi: 10.1007/s10495-023-01888-4
42
PiaoLet al. The Bcr-Abl inhibitor DCC-2036 inhibits necroptosis and ameliorates osteoarthritis by targeting RIPK1 and RIPK3 kinases. BioMed Pharmacother. (2023) 161:114528. doi: 10.1016/j.biopha.2023.114528
43
XiaoJet al. Icariin inhibits chondrocyte ferroptosis and alleviates osteoarthritis by enhancing the SLC7A11/GPX4 signaling. Int Immunopharmacol. (2024) 133:112010. doi: 10.1016/j.intimp.2024.112010
44
WangYGaoW. Effects of TNF-α on autophagy of rheumatoid arthritis fibroblast-like synoviocytes and regulation of the NF-κB signaling pathway. Immunobiology. (2021) 226:152059. doi: 10.1016/j.imbio.2021.152059
45
WeiYXiangHZhangW. Review of various NAMPT inhibitors for the treatment of cancer. Front Pharmacol. (2022) 13:970553. doi: 10.3389/fphar.2022.970553
46
ZhuSLet al. A novel DDIT3 activator dehydroevodiamine effectively inhibits tumor growth and tumor cell stemness in pancreatic cancer. Phytomedicine. (2024) 128:155377. doi: 10.1016/j.phymed.2024.155377
47
KotmayerLet al. GATA2 deficiency and MDS/AML: Experimental strategies for disease modelling and future therapeutic prospects. Br J Haematol. (2022) 199:482–95. doi: 10.1111/bjh.18330
Summary
Keywords
alternative splicing, ankylosing spondylitis, isoform switching, long-read sequencing, programmed cell death
Citation
Cao X, Li P, Peng H, Chen Q, Shi L, Chen D, Chu T and Cheng Y (2026) Integrated long-read transcriptomic profiling of peripheral blood from ankylosing spondylitis patients identifies regulatory shifts and core genes associated with programmed cell death. Front. Immunol. 17:1640271. doi: 10.3389/fimmu.2026.1640271
Received
03 June 2025
Revised
30 December 2025
Accepted
16 February 2026
Published
06 March 2026
Volume
17 - 2026
Edited by
Rosaria Talarico, University of Pisa, Italy
Reviewed by
Yuting Tang, Huazhong University of Science and Technology, China
Sadettin Uslu, Manisa Celal Bayar University, Türkiye
Updates
Copyright
© 2026 Cao, Li, Peng, Chen, Shi, Chen, Chu and Cheng.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yanwei Cheng, 18234069506@139.com; Tianshu Chu, 13568836206@139.com
†These authors have contributed equally to this work and share first authorship
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.