An Important Role in Novel Immune Mechanism and Diagnostic Model of Ankylosing Spondylitis: The CeRNA-ADRB2 Network

Introduction

Ankylosing spondylitis (AS) is an immune-mediated inflammatory disease.¹ It is the second most common disease in the spondyloarthritis group and is mainly characterized by inflammation of the spine and sacroiliac joints, leading to abnormal bone remodeling and ankylosis.² The pathogenesis of AS is not fully understood and may be associated with genetic susceptibility, autoimmunity, and inflammation, with an estimated heritability of over 90% and the highest genetic correlation with HLA-B27.³ The average prevalence of AS per 10,000 people was 31.9 in North America, 23.8 in Europe, 16.7 in Asia, 10.2 in Latin America, and 7.4 in Africa.⁴ AS is characterized by inflammatory back pain that may lead to structural and functional impairment, including kyphotic deformity.⁵ As a disease that may not be completely cured, AS brings continuous pain and suffering to patients, and the cost of its treatment imposes a huge financial burden on families and society. A more in-depth study of the mechanisms of AS is necessary to develop more effective treatment strategies and to perform early diagnosis of the disease. The correlation between AS and genetics has become a hot issue of research since the first confirmation of a genetic factor for AS in 1961.⁶ Studies have shown that 90% ~95% of AS patients are HLA-B27 positive and 1% ~ 2% of HLA-B27 positive people are AS, and this number increases to 15% ~ 20% for those with first-degree relatives who are infected.^7,8 However, even the most highly regarded HLA-B27 did not contribute more than 20% to the overall genetic contribution,⁹ suggesting that other factors also play a key role in the pathogenesis of AS.

The intricate pathogenesis of AS involves the dynamic interplay between cellular immunity and humoral immunity. Within the realm of AS, activated T lymphocytes, particularly the CD4+ subset, orchestrate the chronic inflammatory features of the ailment.¹⁰ These T cells are postulated to infiltrate the affected tissues, such as the spine, hip joints, and sacroiliac joints, thereby releasing inflammatory cytokines like IL-17 and TNF, perpetuating the cascade of inflammatory reactions.¹¹ Conversely, the humoral immune response, epitomized by B cells and their associated antibodies, may also play a role in the autoimmune facets of AS.¹² In the context of AS, aberrations in humoral immunity give rise to the production of autoantibodies which may assail the body’s own joint tissues.² Recent studies have suggested that the metabolic microenvironment may play an important role in the arthritic progression of AS.^13–15 As a part of the metabolic microenvironment, the emergence of autoantibodies and the continued proliferation of T cells support the idea that AS is an autoimmune disease.² Studies on human genetics combined with animal models have shown that CD8+ T cells and CD4+ T cell subsets, mediated by the IL-23/IL-17 pathway and TNF-α, participate in important stages of the pathogenesis of AS,^16,17 which may lead to kyphotic deformity. These studies suggest that immune factors may play central roles in the development of AS, which still needs to be explored.

Tremendous progress has been made in exploring immune-related genes(IRGs) of AS, thanks to advances in high-throughput sequencing technology. These discoveries include genes encoding cytokine receptors, transcription factors, signaling molecules, and transporters, which provide insights into the pathogenesis of AS.¹⁸ However, the reliability of biomarkers may be weakened by insufficient sample size or being limited to a single level of sample histological types.¹⁹ Complex genetic bases underlie the molecular mechanisms of human diseases, requiring sophisticated sequencing techniques and advanced analytical methods. Interdisciplinary approaches, such as bioinformatics, contribute to a comprehensive understanding of disease mechanisms.²⁰

In this study, we aimed to identify DEIRGs in AS by using high-throughput sequencing of both ISL and PB. To gain insights into the roles played by DEIRGs in disease progression and immune cell infiltration, bioinformatics analysis was conducted, including immune infiltration analysis and functional enrichment. Moreover, we have developed a predictive model for AS and utilized immune biomarkers to subtype AS, which helps us to further understand the regulatory role of key genes on AS. Overall, our study contributes to a comprehensive understanding of AS and provides potential therapeutic targets.

Materials and Methods

Source of ISL Samples

Between September 2018 and December 2020, lumbar posterior entry surgery was performed on two groups of patients at the Department of Spinal Orthopedics, Affiliated Hospital of Guangxi Medical University: those diagnosed with ankylosing spondylitis and spinal kyphotic deformity (AS-KD group) based on the 1984 New York Revised Diagnostic Criteria for AS, and those with lumbar disc herniation (non-AS group) confirmed by MRI and Pfirrmann grade 3–4. Patients with other rheumatic diseases or autoimmune conditions such as lupus erythematosus were excluded from both groups. ISL samples were collected from both groups using a randomization principle. The collected ISL was cut into 0.5 cm pieces and washed quickly with saline to remove residual blood components. RNA Keeper was then added in a 5:1 ratio to the tissue within 10 minutes of collection and left overnight at 4°C for adequate penetration of the solution into the tissue. The tissue samples were stored at −80°C until analysis.

Gene Sequencing

The mRNA, lncRNA, and miRNA of the ISL samples were obtained by competing endogenous RNA sequencing (ceRNA sequencing). The total RNA was amplified and labeled using the Agilent Expression Profiling Kit, Low Input Quick Amp Labeling Kit and One-Color (Cat. # 5190–2305, Agilent Technologies, Santa Clara, CA, USA). The labeling process was divided into three steps: single label spike-in preparation, reverse transcription, fluorescent labeling, and purification of the labeled RNA using the RNeasy Mini Kit (Cat.# 74106, QIAGEN, GmbH, Germany). The Gene Expression Hybridization Kit (Cat.# 5188–5242, Agilent Technologies, Santa Clara, CA, USA) was used in a rolling hybridization oven (Cat.# G2545A, Agilent Technologies, Santa Clara, CA, USA) at 65°C, 10 rpm for 17 hours. Hybridized RNA was sampled at 1.65 μg, and washed in staining dishes (Cat. # 121, Thermo Shandon, Waltham, MA, USA) using the Gene Expression Wash Buffer Kit (Cat.# 5188–5327, Agilent Technologies, Santa Clara, CA, USA). The completed microarrays were scanned using an Agilent Microarray Scanner (Cat.# G2565CA, Agilent Technologies, Santa Clara, CA, USA) with the software setting Dye channel: Green, Scan resolution=3μm. Data were read using Feature Extraction software 10.7 (Agilent Technologies, Santa Clara, CA, USA).

Gene Expression Microarrays for PB of AS

To investigate gene expression patterns associated with AS, the GEO database of NCBI was searched by using keywords: Ankylosing Spondylitis, Spondyloarthropathy and Axial Spondyloarthritis. Two datasets GSE73754 and GSE25101 were identified for further analysis. To merge the datasets and eliminate any batch effects, the R package “inSilicoMerging”²¹ and the approach described by Johnson et al²² were used, respectively. After applying these methods, we obtained a gene expression matrix that was normalized and free of batch effects, enabling us to perform downstream analyses.

Screening of DE mRNAs

The R package “limma” was utilized, which is a differential expression screening method based on a generalized linear model,²³ to identify DE mRNAs between AS-KD and control groups. To account for potential differences in protein concentrations across tissue sources, a more stringent threshold for DE mRNAs was set in the ISL samples, using a difference multiplier of 2.0 and p-value < 0.05, compared to a threshold of a difference multiplier of 1.2 and p-value < 0.05 for PB samples. T-tests were used to validate the statistical significance of the DE mRNA differences. Moreover, a list of 1902 immune-related genes (Supplementary Material 1) was obtained from the immProt database (https://www.immport.org/shared/). DEIRGs were identified by using a Venn diagram.

Genes that show differential expression across different tissues may have important biological significance. Therefore, we intersected the DEGs between ISL and PB to obtain the key genes.

Functional Enrichment Analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database for storing genome, biological pathways, disease, and drug information.²⁴ From the KEGG REST API, the latest gene annotation for KEGG Pathway was obtained. In addition, gene ontology (GO) analysis was carried out by using the GO annotation of the genes obtained from the R package org.Hs.eg.db (version 3.1.0). The R package clusterProfiler (version 3.14.3) was used for enrichment analysis.

To evaluate the relevant pathways and molecular mechanisms, gene set enrichment analysis (GSEA) was used. GSEA can be used to study the synergistic changes of genes on phenotype.³ The GSEA software (version 3.0) was obtained from the GSEA website and downloaded a subset of c2.cp.kegg.v7.4.symbols.gmt from the Molecular Signatures Database.

Immune Infiltration Analysis

IOBR²⁵ is a computational tool for immuno-tumor biology research. Based on gene expression spectrum, we used R package IOBR and selected CIBERSORT²⁶ method to calculate 22 immune infiltration cell scores per sample. The correlation between immune cells and genes was analyzed by correlation heat map.

LASSO Algorithm for Screening Biomarker

In this study, we aimed to identify biomarkers for AS by analyzing gene expression data of 26 DEIRGs in 104 human PB samples. The R package “glmnet” was utilized for regression analysis using the lasso-cox method, integrating disease status and gene expression data. A 5-fold cross-validation was employed to obtain the optimal model and set the lambda value to 0.0285, resulting in the identification of 12 biomarkers.

To evaluate the predictive significance of these biomarkers, the R package SURVIVAL was used to integrate disease status and expression data of the 12 biomarkers. We then applied the Cox method and performed receiver operating characteristic curve (ROC) analysis using the R package “pROC”, yielding the area under the curve (AUC).

Finally, a nomogram model was developed based on the selected biomarkers to predict the prevalence of AS. To screen for biomarkers with the most diagnostic value for the disease, a multivariable stepwise regression analysis was used.

Collection of External Validation Set for Diagnostic Models

A new search in the NCBI GEO database was conducted using the keywords ankylosing spondylitis, inflammatory disease, autoimmune disease, and peripheral blood. As a result, the datasets GSE11886 and GSE134290 were obtained, which were used as validation sets for the diagnostic model. It is worth noting that we excluded samples treated with IFN in GSE11886 to eliminate confounding factors.

Protein-Protein Interaction (PPI) Network

The STRING database is widely used to search for relationships between known proteins and predicted proteins. Using the STRING database, a PPI network for DE mRNAs, DEIRGs, and 12 biomarkers were constructed. Cytoscape (v3.7.2) software was used to display the PPI network models.

Consistency Clustering Analysis

ConsensusClusterPlus²⁷ is an unsupervised clustering method that utilizes resampling to draw multiple subsamples from the input dataset, specifies the number of clusters (k), and evaluates cluster quality across different values of k. The optimal number of clusters is determined by an empirical cumulative distribution function plot. In our study, we employed ConsensusClusterPlus and used the expression matrix of 12 biomarkers to cluster 68 AS patients into two distinct immune patterns: Cluster1 and Cluster2.

Construct a ceRNA Network for ADRB2

The R package “limma” was used to obtain DE miRNAs and DE lncRNAs. Subsequently, we aimed to construct a ceRNA network which regulates the down-regulated biomarker ADRB2 and involves lncRNA-miRNA-mRNA interactions. The miRNAs that target ADRB2 were predicted by utilizing the TargetScan database.²⁸ Then, we compared these predictions with the DE miRNAs in ISL and identified their intersection. This allowed us to obtain a list of DE miRNAs that were likely to target ADRB2. Furthermore, the ENCORI database²⁹ was employed for predicting the lncRNAs involved in the ceRNA network.

Immunohistochemistry (IHC)

The experimental group consisted of ISL samples from AS-KD, and the control group consisted of ISL samples from patients with lumbar disc herniation without AS. IHC was performed to compare the differences in ADRB2 expression between the experimental and control groups. Immunohistochemical sections with complete staining were obtained from paraffin sections after experimental procedures such as roasting, dewaxing, hydration, and antigen repair. The images were observed under the microscope. The criteria were scored as follows: (1) by degree of staining (0 points: negative staining, 1 point: light yellow, 2 points: light brown, 3 points: dark brown), and (2) positive area (1 point: 0–25%, 2 points: 26–50%, 3 points: 51–75%, 4 points: 76–100%). The final two scores were summed (0–3 points: negative expression, 4–7 points: positive expression). ImageJ software was used to assess the positivity rate of all IHC images, and IBM SPSS Statistics 22.0 was used to perform a t-test for ADRB2 positivity in the experimental and control groups.

Results

Collection and Processing of Samples

The study flow is shown in Figure 1. The ISL samples were obtained from three lumbar disc herniation patients (Figure 2A) and three AS-KD patients (Figure 2B), with detailed information provided in clinical data of patients (Supplementary Material 2). Through ceRNA sequencing, we obtained mRNAs, miRNAs, and lncRNAs of the ISL samples, To meet the analytical requirements for gene expression data, we performed normalization and log transformation on the data (Figure 2C and D). The mRNA expression matrix of PB samples from AS patients was obtained from the GEO database (Supplementary Material 3), and the batch effect was removed between GSE73754 and GSE25101 (Figure 2E and F).

Figure 1 Flow chart. Bioinformatics methods used: KEGG (Kyoto Encyclopedia of Genes and Genomes), GO (Gene Ontology), GSEA (Gene Set Enrichment Analysis), LASSO (Least Absolute Shrinkage and Selection Operator), Cibersort (Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts), and PPI (Protein-Protein Interaction). Abbreviations: AS (Ankylosing Spondylitis), AS-KD (Ankylosing Spondylitis with Kyphotic Deformity), DE mRNAs (Differentially Expressed messenger RNAs), IRGs (Immune-Related Genes), DEIRGs (Differentially Expressed Immune-Related Genes), PB (Peripheral blood), and ISL (Interspinous ligaments), ADRB2 (beta-2 adrenergic receptor).

Figure 2 Sample collection and processing: (A) Full spine lateral radiograph of the control group. (B) Full spine lateral radiograph of a patient with AS-KD. (C and D) Normalization and log transformation were applied to the ceRNA sequencing data of ISL samples from AS-KD patients (D: disease group; H: healthy group). (E) Gene expression data from PB samples of AS patients before batch correction. (F) Gene expression data from PB samples of AS patients after batch correction.

Identification of DEIRGs

The differential analysis of ISL samples revealed 787 DE mRNAs (Supplementary Material 4), with 264 upregulated and 523 downregulated (Figure 3A). A heat map depicting the differential expression was presented (Figure 3B). By intersecting the DE mRNAs with 1902 IRGs (Supplementary Material 1), a total of 76 DEIRGs were identified (Figure 3C).

Figure 3 Screening of DEIRGs: (A) Volcano plot of gene expression differences in ISL samples from AS-KD patients. (B) Heatmap of gene expression differences in ISL samples from AS-KD patients. (D: Disease group, H: Control group). (C) Intersection of ISL-DE mRNAs and IRGs. (D) Volcano plot of gene expression differences in PB samples from AS patients. (E) Heatmap of gene expression differences in PB samples from AS patients. (F) Intersection of PB-DE mRNAs and IRGs. (G) Intersection of ISL-DE mRNAs, PB-DE mRNAs, and IRGs. (H) Expression levels of ADRB2 in ISL. (AS: AS-KD group, NO-AS: Control group). (I) Expression levels of ADRB2 in PB. (J) Expression levels of 26 PB-DEIRGs obtained by intersecting PB-DE mRNAs and IRGs. “*” represents p < 0.05, “**” represents p < 0.01, “***” represents p < 0.001, “****” represents p < 0.0001.

Similarly, the differential analysis of PB samples identified 155 DE mRNAs (Supplementary Material 5), consisting of 97 upregulated and 58 downregulated transcripts (Figure 3D). The corresponding heat map was shown in Figure 3E. Intersection of the DE mRNAs with 1902 IRGs led to the identification of 26 DEIRGs (Figure 3F) with 13 upregulated and 13 downregulated genes (Figure 3J).

Moreover, after intersecting the ISL-DE mRNAs, PB-DE mRNAs, and 1902 IRGs, gene ADRB2 was identified (Figure 3G). Notably, the expression of ADRB2 was significantly downregulated in both ISL and PB samples from AS patients (Figure 3H and I).

Functional Enrichment Analysis of DEIRGs

The KEGG analysis for ISL-DEIRGs revealed 29 pathways with p < 0.05, including key pathways such as cytokine-cytokine receptor interaction, chemokine signaling pathway, and neuroactive ligand-receptor interaction (Figure 4A). Similarly, PB-DEIRGs showed significant involvement in 17 pathways with p < 0.05, such as natural killer cell-mediated cytotoxicity, B cell receptor signaling pathway, and complement and coagulation cascades, among others (Figure 4B).

Figure 4 Functional enrichment analysis of DEIRGs and ADRB2: (A) The KEGG analysis results of ISL-DEIRGs. (B) The KEGG analysis results of PB-DEIRGs. (C) The GO analysis results of ISL-DEIRGs. (D) The GO analysis results of PB-DEIRGs. (E) The KEGG pathways enriched by ADRB2. (F) After re-grouping peripheral blood samples based on the expression level of ADRB2, GSEA analysis was performed (L: ADRB2 low-expression group; H: ADRB2 high-expression group; ES: Enrichment Score).

The GO analysis revealed that ISL-DEIRGs were predominantly involved in regulating biological processes (GObp) related to immune system function, cell surface receptor signaling pathway, and regulation of cell communication, among others. These genes also played a crucial role in regulating cellular components (GOcc) such as extracellular region, extracellular space, and extracellular matrix, among others. Furthermore, the molecular functions (GOmf) of these genes included key signaling receptor binding, growth factor activity, and receptor ligand activity, among others (Figure 4C). The GObp associated with PB-DEIRGs were primarily focused on immune system process, immune response, and response to stress, while the GOcc involved endomembrane system, vesicle, and intrinsic component of the membrane, and GOmf were mainly focused on serine-type peptidase activity, signaling receptor activity, and molecular transducer activity (Figure 4D).

In particular, ADRB2 was actively involved in regulating several KEGG pathways such as neuroactive ligand-receptor interactions, cAMP signaling pathway, cGMP-PKG signaling pathway, and calcium signaling pathway, among others (Figure 4E). Additionally, the GSEA analysis results for ADRB2 showed significant involvement in complement and coagulation cascades, TGF-β signaling pathway, and natural killer cell-mediated cytotoxicity, among others (Figure 4F).

Common KEGG and GO Terms in Both ISL and PB

The analysis of enriched DEIRGs in both ISL and PB revealed 7 common KEGG pathways that are mainly involved in the regulation of cells such as B cells, T cells, NK cells, and osteoclasts (Figure 5A).

Figure 5 The common pathways of KEGG and GO analysis for DEIRGs in PB and ISL samples: (A) The common intersection of KEGG analysis. (B) The common intersection of GO biological process (GObp) analysis. (C) The common intersection of GO cellular component (GOcc) analysis. (D) The common intersection of GO molecular function (GOmf) analysis.

In addition, the DEIRGs of ISL and PB were found to be enriched in 3 common GObp terms, including response to chemical, regulation of response to stimulus, and immune system process (Figure 5B). Furthermore, the functional enrichment analysis of DEIRGs of ISL and PB also identified 6 common GOcc terms, such as extracellular region, vesicle, and endomembrane system, among others (Figure 5C). Finally, the DEIRGs of ISL and PB were enriched in 2 common GOmf terms, including signaling receptor activity and molecular transducer activity (Figure 5D).

Identification of Biomarkers and Construction of Disease Diagnosis Models

The LASSO algorithm was applied to identify 12 biomarkers that were most significantly associated with the development of AS (Figure 6A) and the optimal diagnostic model was obtained through 5-fold cross-validation (Figure 6B). Notably, ADRB2 was included among the identified set of biomarkers. Subsequently, multivariable survival analysis demonstrated that ADRB2, IL2RB and PRF1 was important disease prediction indicators (Figure 6C). The diagnostic model yielded a high AUC value of 0.89, indicating excellent performance in diagnosing AS (Figure 6D).

Figure 6 Diagnostic model of AS: (A) Screening of diagnostic biomarkers from 26 PB-DEIRGs using the LASSO algorithm. (B) 5-fold cross-validation was performed to obtain the optimal diagnostic model. (C) Forest plot of multivariate survival analysis in AS patients. (Hazard Ratio=1: no harm, Hazard Ratio<1: reduced harm, Hazard Ratio>1: increased harm). (D) ROC curves for risk scores in AS diagnosis. (E) ROC curves of 12 diagnostic biomarkers.

Interestingly, when we analyzed the predictive ability of the 12 biomarkers independently, half of the genes exhibited an AUC value of ≥0.7, while ADRB2, CD81, and IL2RB had an AUC value of ≥0.75 (Figure 6E), further supporting their potential as diagnostic biomarkers.

To predict the prevalence of AS, a nomogram model was constructed based on the expression of the 12 identified biomarkers (Figure 7A). The calibration curves indicated that the predictions of the nomogram model were similar to those of the ideal model (Figure 7B), and the risk scores calculated by decision curve analysis were superior to those of the random model (Figure 7C).

Figure 7 Optimization and evaluation of diagnostic models: (A) Construct a nomogram for the diagnosis of AS using 12 diagnostic biomarkers. (B) The calibration curve of the nomogram. (C) The decision analysis curve of the nomogram (“All”: Predictive analysis using all factors related to the disease; “TB genes”: Predictive analysis using 12 diagnostic biomarkers; “None”: Predictive analysis using only non-genetic factors). (D) Diagnosing model with external validation.

The constructed model equation is: RiskScore = Correspondence coefficient * gene expression (−0.0096HLA-DQA1 + 0.0183PPBP - 0.2341PTGDS + 0.1133SLC40A1 + 0.0309TNFSF10 - 0.2980CD81 + 0.0107C5AR1 + 0.0676ADIPOR1 - 0.3161ADRB2 - 0.3119IL2RB - 0.0443SH2D1B + 0.1393PRF1). (Only the first four decimal places are shown.)

External Validation of Diagnostic Models

We utilized the gene expression matrix of the validation set and combined it with the diagnostic model formula to calculate the risk score for each sample. Finally, the ROC curve of the validation set was obtained, and the results showed an AUC value of 0.70 (Figure 7D).

Correlation Analysis and Interaction Among DEIRGs

In PB samples, ADRB2 was significantly positively correlated with PRF1, IL2RB, CD81, SH2D1B, and PTGDS. Additionally, ADRB2 was significantly negatively correlated with C5AR1, TNFSF10, and SLC40A1, with the highest correlation coefficient found with CD81 (Figure 8A). In AS samples of PB, ADRB2 was significantly positively correlated with PRF1, IL2RB, and CD81, with the highest correlation coefficient found with CD81 (Figure 8B). In control samples of PB, ADRB2 was significantly positively correlated with PRF1, IL2RB, CD81, SH2D1B, and PTGDS. A significant negative correlation was found with TNFSF10, and the highest correlation coefficient was identified with CD81 (Figure 8C).

Figure 8 Correlation analysis and interaction among DEIRGs: (A) Correlation analysis of 12 diagnostic biomarkers in PB of all samples (“*” represents p < 0.05, “**” represents p < 0.01, “***” represents p < 0.001, “****” represents p < 0.0001; numbers indicate correlation coefficients). (B) Correlation analysis of 12 diagnostic biomarkers in PB of AS samples. (C) Correlation analysis of 12 diagnostic biomarkers in PB of control samples. (D) PPI network of DE mRNA in PB samples. (E) PPI network of DEIRGs in PB samples. (F) PPI network of diagnostic biomarkers in PB samples. (Node: gene; Link: interaction between genes; Node size: the larger the size, the more genes with interaction.).

The PPI network of DE mRNAs from PB samples included 80 genes. Among them, PRF1 interacted with 22 genes, IL2RB interacted with 15 genes, and PPBP interacted with 8 genes (Figure 8D). The PPI network of DEIRGs in PB samples included 17 genes, with PRF1 interacting with 8 genes, IL2RB interacting with 5 genes, and PPBP interacting with 5 genes (Figure 8E). Among the 12 biomarkers, PRF1 interacted with IL2RB, PPBP, and TNFSF10, while CD81 interacted with SH2D1B (Figure 8F).

Identifying Immune Patterns of AS and Functional Enrichment Analysis of DEGs Among Immune Patterns

Consistent clustering yielded two immunisation patterns of 68 AS patients. Cluster1 (C1) contained 29 samples and cluster2 (C2) contained 39 samples. The PCA plot demonstrated successful stratification of AS patients into two groups (Figure 9A and B). Differential analysis between C1 and C2 showed 104 Cluster-DEGs, of which 79 were up-regulated and 25 down-regulated (Figure 9C).

Figure 9 Cluster analysis results of AS patients: (A and B) Cluster analysis divided AS patients into Cluster1 (C1) and Cluster2 (C2). (C) Volcano plot of differential gene expression analysis between C1 and C2. (D) KEGG analysis of DEGs between C1 and C2 (Cluster-DEGs). (E) GO analysis of Cluster-DEGs. (F) Expression levels of ADRB2 in C1 and C2. (G) GSEA analysis of the Cluster. (H) GSEA analysis of ADRB2 in the Cluster (Cluster-ADRB2-GSEA). (I) Intersection of Cluster-KEGG, Cluster-GSEA, and Cluster-ADRB2-GSEA.

KEGG analysis of Cluster-DEGs showed enrichment in pathways such as T cell receptor signaling pathway, natural killer cell-mediated cytotoxicity, and Th1 and Th2 cell differentiation (Figure 9D). The GObp terms regulated by Cluster-DEGs included immune system process, immune response, and regulation of immune system process (Figure 9E). The GOcc terms regulated by Cluster-DEGs included intrinsic component of membrane, integral component of membrane, and plasma membrane part (Figure 9E). The Cluster-DEGs were also found to regulate GOcc terms such as signaling receptor activity, molecular transducer activity, and transmembrane signaling receptor activity (Figure 9E).

The expression levels of ADRB2 were observed to be significantly higher in the C1 group (Figure 9F). Further analysis by overlapping the results from Cluster GSEA (Figure 9G), ADRB2 GSEA (Figure 9H), and Cluster KEGG analysis identified a common pathway, T-cell receptor signaling pathway (Figure 9I). The gene set associated with this pathway was found to be significantly enriched in the group with high expression of ADRB2, strongly suggesting that high levels of ADRB2 expression may activate this pathway (Figure 9H).

The Differences in Immune Infiltration Among Immune Patterns of AS

The stacked immune cell graph demonstrated that neutrophils and CD8+ T cells were the primary components (Figure 10A). In the AS group, low infiltration was observed for CD8+ T cells, activated NK cells, and M2 macrophages, while naive CD4+ T cells and M0 macrophages showed high infiltration (Figure 10B). A significant positive correlation was found between ADRB2 and CD8+ T cells and activated NK cells in AS patients, whereas a significant negative correlation was observed between ADRB2 and naive CD4+ T cells and M0 macrophages (Figure 10C).

Figure 10 Immune infiltration analysis of different clusters in PB of AS patients: (A) Stacked plot of immune cells in 68 AS patients. (X-axis: sample ID; different colors: immune cells). (B) Infiltration levels of immune cells in PB. (D: disease group, H: control group). (C) Correlation analysis between diagnostic biomarkers and immune cells in PB of AS patients. (X-axis: diagnostic biomarkers; Y-axis: immune cells; upper left of the square: correlation coefficient; lower right of the square: -log10(p value)). (D) Infiltration levels of immune cells in different AS clusters (“*” represents p < 0.05, “**” represents p < 0.01, “***” represents p < 0.001, “****” represents p < 0.0001). (E) Correlation analysis between immune cells in C1. (F) Correlation analysis between immune cells in C2. (G) Correlation analysis between diagnostic biomarkers and immune cells in C1. (H) Correlation analysis between diagnostic biomarkers and immune cells in C2.

In the C1 group, neutrophils exhibited low infiltration, and activated CD4+ memory T cells and resting NK cells showed high infiltration (Figure 10D). Furthermore, a significant positive correlation was found between neutrophils and gamma delta T cells in the C1 group (Figure 10E). In the C2 group, a significant positive correlation was also observed between neutrophils and gamma delta T cells, but a significant negative correlation was noted between neutrophils and CD8+ T cells and Monocytes (Figure 10F).

Additionally, a significant negative correlation between neutrophils and ADRB2 was observed in the C1 group (Figure 10G), whereas in the C2 group, a significant positive correlation was found between neutrophils and ADRB2 (Figure 10H).

Constructing the ceRNA Network of ADRB2

The analysis of miRNAs in ISL showed that among the 58 DE miRNAs (Supplementary Material 6), 18 were up-regulated and 40 were down-regulated (Figure 11A), which was further confirmed by the heat map (Figure 11B). In addition, differential analysis of lncRNA microarrays in ISL revealed 651 DE lncRNAs (Supplementary Material 7), of which 322 were up-regulated and 329 were down-regulated (Figure 11C), and this was also supported by the heat map (Figure 11D). By taking the intersection, we identified hsa-miR-513a-5p and its target lncRNAs (DLX6-AS1, AL008721.2, RBM26-AS1, AC006504.5, ATP1A1-AS1, MAGI2-AS3, FLNB-AS1, AC007342.4) (Figure 11E and F). Finally, the ceRNA network of ADRB2 was constructed using hsa-miR-513a-5p and its target lncRNAs (Figure 11G). Furthermore, the PPI network (Figure 11H and I) illustrated the proteins in ISL that interacted with ADRB2.

Figure 11 The ceRNA network of ADRB2: (A and B) Volcano plot and heatmap of DE miRNAs in ISL samples of AS-KD patients. (C and D) Volcano plot and heatmap of DE lncRNAs in ISL samples of AS-KD patients. (E) Venn diagram of upregulated DE miRNAs with miRNAs targeting ADRB2. (F) Venn diagram of DE lncRNA with lncRNAs targeting hsa-miR-513a-5p. (G) The ceRNA interaction network of ADRB2 (red: upregulated in AS-KD; blue: downregulated in AS-KD). (H) PPI network of ADRB2 in ISL. (I) PPI network of DE mRNAs in ISL.

Immunohistochemistry (IHC) of ADRB2

To verify the protein expression level of ADRB2 in ISL, immunohistochemistry (IHC) experiments were conducted. The results showed that the specific expression of ADRB2 was significantly lower in the AS-KD group compared to the control group (Figure 12A and B). The positive expression rate of ADRB2 was found to be 11.37% in the AS-KD group and 35.09% in the control group, with a p-value less than 0.001 (Figure 12C). These findings confirm the accuracy of our analysis in identifying the differential expression of ADRB2 between the AS-KD and control groups.

Figure 12 Immunohistochemistry of ADRB2 in ISL: (A and B) The expression of ADRB2 in ISL. (C) The data visualization of ADRB2 positivity rate results (AS: AS-KD group, NC: control group).

Discussion

Recent advances in immunological techniques have greatly improved our understanding of the role played by different cell populations, as well as the innate and acquired immune systems, in the development of diseases.³⁰ In the specific case of AS, abnormal immune cells and cytokines have been found in patients, which strongly suggests that immune factors play a major role in its pathogenesis.³¹ Therefore, in this study, we aimed to investigate the pathogenesis of AS from both the perspectives of the ISL and PB samples, identifying common components between them, and analyzing the immune microenvironment in order to identify new targets for AS.

In our study, we screened ISL samples for 76 DEIRGs, and KEGG analysis revealed their involvement in pathways related to osteoclast differentiation, TGF-β signaling, and chemokine signaling. On the other hand, PB samples were screened for 26 DEIRGs, and KEGG analysis showed their involvement in pathways related to osteoclast differentiation, natural killer cell-mediated cytotoxicity, and Th17 cell differentiation. In both ISL and PB samples, most of the KEGG-enriched pathways were associated with inflammation and immunity, which is consistent with previous studies on AS.³² Previous research has shown that chondrocyte differentiation occurs in the ligaments of early AS patients, with calcification occurring after cartilage formation. Furthermore, a large number of giant osteoclasts have been found in the inflammatory environment of ligaments and on the bone surface of calcified cartilage. The activity of these osteoclasts produces higher levels of TGF-β, which in turn induces new bone formation in the ligament.³³ Interestingly, our findings suggest that some DEIRGs in the ISL of AS-KD patients may significantly regulate osteoclast production, leading to increased levels of TGF-β and subsequent cartilage calcification, spinal fusion, and deformity.

In this study, 1902 IRGs, 787 DE mRNAs from ISL and 155 DE mRNAs from PB were analyzed to identify immune-related key genes: ADRB2. ADRB2, which encodes a beta-2 adrenergic receptor that mediates catecholamine-induced activation of adenylate cyclase via G proteins, was found to be significantly downregulated in the ISL and PB of AS patients. Previous studies have suggested that ADRB2 may be associated with bronchial asthma,³⁴ osteoarthritis,³⁵ rheumatoid arthritis,³⁶ and cancer,³⁷ but its effect on AS remains unclear. Our results showed that ADRB2 was involved in multiple KEGG pathways, including cAMP signaling, cGMP-PKG signaling and calcium signaling. Cyclic adenosine monophosphate (cAMP) was discovered in 1958 as the first identified second messenger for intracellular or extracellular ligand action.³⁸ In the immune system, cAMP plays a crucial role in regulating both pro-inflammatory and anti-inflammatory activities. Drugs that increase intracellular cAMP levels have been shown to decrease pro-inflammatory mediators while increasing anti-inflammatory factors.³⁹ The downregulation of ADRB2 may lead to decreased activation of the cAMP signaling pathway, ultimately leading to a decrease in intracellular cAMP levels. This decrease in cAMP levels could contribute to an increase in pro-inflammatory substances released from immune cells and may ultimately exacerbate the progression of inflammation in patients with AS. Inflammatory stimulation induces activation of the NO/cGMP/PKG-I signaling pathway in the spinal cord, which promotes the phosphorylation of mitogenic proteins and contributes to inflammatory pain.⁴⁰ Prior studies have shown that activation of the peripheral or local NO/cGMP/PKG-I signaling pathway can help alleviate inflammatory pain.⁴¹ Our findings suggested that the cGMP-PKG pathway may serve as a potential therapeutic mechanism for treating AS-associated inflammatory pain. Calcium signaling is essential for the development and function of all T-cell lineage cells, and calcium-regulating neurophosphatase inhibitors, such as cyclosporine and tacrolimus, are currently used to treat autoimmune diseases.^41–44 ADRB2 may regulate T-cell activity through this pathway and potentially affect AS. Overall, our study provides new insights into the potential mechanisms underlying the development and progression of AS, and highlights the importance of ADRB2 in modulating several critical signaling pathways involved in AS pathology, which may represent a novel therapeutic target for AS.

Another important genetic factor in the pathogenesis of AS is the non-major histocompatibility complex (non-MHC), which plays a regulatory role in the activation and differentiation of CD4+ or CD8+ T cells.^45,46 Our analysis of immune infiltration levels in PB samples revealed that activated naive CD4+ T cells and M0 macrophages were significantly more abundant in AS patients, and their levels were negatively correlated with ADRB2 expression. Conversely, the infiltration levels of CD8+ T cells and NK cells were significantly reduced, and their levels were positively correlated with ADRB2 expression. These results suggested that ADRB2 may plays a key role in regulating immune cells in AS. Immune cells express adrenergic receptors, which allow the sympathetic nervous system (SNS) to directly regulate immune function.⁴⁷ It has been reported that ADRB2 signaling in immune cells can attenuate the severity of autoimmune diseases,^48–50 indicating that low expression of ADRB2 in AS patients may contribute to disease exacerbation. In an experimental autoimmune encephalomyelitis (EAE) model of multiple sclerosis (MS), cAMP early blockers (ICER) were induced to inhibit IL-2, IFN-gamma and GM-CSF production via the ADRB2 pathway, thereby limiting T-cell autoimmune responses in EAE. Deficiency of ADRB2 signalling in immune cells leads to an increase in pathogenic T cells and the progression of EAE.⁴⁹ These studies further suggest that ADRB2 deficiency may influence the infiltration levels of CD4+ T cells and CD8+ T cells. Saligrama et al reported in Nature that clonally expanded CD8+ T cells can suppress specific CD4+ T cells,⁵¹ which may help explain why CD8+ T cells exhibit a low infiltration level and CD4+ T cells exhibit a high infiltration level in our study. Moreover, interleukins produced by immune cells, such as interleukin-4 (IL-4) and interleukin-10 (IL-10), play different roles in AS. IL-4 is pro-inflammatory,⁵² whereas IL-10 has potent immunosuppressive and anti-inflammatory functions.⁵³ Previous studies have found a higher proportion of IL4+CD8+ T cells and a lower proportion of IL10+CD8+ T cells in peripheral blood mononuclear cells of AS patients.^54,55 Further studies have shown that norepinephrine can induce immune cells to secrete IL-10 by acting on ADRB2 receptors located on immune cells.⁴⁷ Thus, the low expression of ADRB2 in AS patients may lead to reduced IL-10 secretion by immune cells. This suggests that the poorly infiltrating type of CD8+ T cells in AS patients are mainly IL-10+CD8+ T cells.

Contactin-1 (Cntn1), a neuronal cell surface adhesion/recognition molecule, has been linked to reduced intestinal ADRB2 levels. The loss of ADRB2 can directly impair intestinal function and disrupt the balance of immunity.⁵⁶ Previous studies have suggested that impaired bowel function may be a mechanism of AS.² These findings suggest a potential connection between ADRB2 and AS, and indicate that the gut might be a future focus of research for further exploring this link.

Hsa-miR-513a-5p may be involved in various biological processes by binding to the 3’ untranslated region (UTR) of ADRB2.⁵⁷ In AS-KD patients, the upregulation of hsa-miR-513a-5p in the ISL may lead to the silencing of ADRB2. The DE lncRNAs (DLX6-AS1, AL008721.2, RBM26-AS1, AC006504.5, ATP1A1-AS1, MAGI2-AS3, FLNB-AS1, AC007342.4) identified in this study may act as “sponges”, competing with hsa-miR-513a-5p for binding sites, degrading the miRNA, and releasing ADRB2 from miRNA-mediated inhibition.⁵⁸ Recent studies have shown that miR-15a-5p directly targets ADRB2 and inhibits IL-3-induced inflammatory responses in the nasal epithelium.⁵⁹ The ceRNA network of ADRB2 may help discover more comprehensive regulatory mechanisms and refine new therapeutic strategies. The PPI network results showed that ADRB2 interacts with angiotensinogen (AGT) and uncoupling protein 3 (UCP3). AGT promotes IL-6-induced inflammatory responses in osteoarthritis through the JAK2/STAT3 pathway,⁶⁰ suggesting its possible association with inflammation in AS-KD. UCP3 regulates the differentiation of immune cells and is a potential therapeutic target for the treatment of Th17 cell-mediated autoimmune diseases.⁶¹ These interacting proteins further reveal ADRB2’s role in the development of AS-KD.

We constructed predictive models for AS, and our findings suggested that ADRB2, IL2RB, and CD81 have good predictive ability with ROC values of 0.75 or more. The reappearance of ADRB2 further highlights its importance in the development of AS. Interestingly, all three genes had a risk score of less than 1 and were downregulated in AS, indicating their potential role as protective factors. IL2RB has been previously identified as one of the key genes involved in AS,⁶² and our study supports this finding. We also observed a significant positive association between IL2RB and resting NK cells,⁶³ providing further evidence of its relevance to AS. Cluster of differentiation 81 (CD81), a widely expressed cell surface protein, has been shown to affect the proliferation and differentiation of B and T cells.⁶⁴

In neoplastic diseases, clustering algorithms are often used to identify different molecular subtypes that help explain the heterogeneity of the disease.⁶⁵ Not all AS patients progress to kyphotic deformity, to further investigate the pathogenesis of AS-KD, we utilized unsupervised clustering to identify two AS immune patterns labeled as C1 and C2. Our analyses revealed the T-cell receptor signaling pathway to be involved in both the C1 and ADRB2 high-expression groups, as determined by GSEA and KEGG analyses as well as ADRB2 GSEA analysis. The gene set related to this pathway had ES values of 0.4969 and 0.4410 in the C1 and ADRB2 high-expression groups, respectively, indicating that ADRB2 may have an impact on key genes within this pathway, resulting in its activation in both groups. The T cell receptor (TCR) complexes play a crucial role in mediating the activation of T cells, which is essential for maintaining normal T cell homeostasis and immune response. The strength of TCR signaling determines the differentiation of T cells into various effector cell lines.⁶⁶ Previous studies have demonstrated that the T cell receptor signaling pathway is critical in the development of AS.⁶⁷ Abnormal activation of neutrophils at the site of disease predilection has been shown to be associated with AS.⁵⁶ Our findings revealed an intriguing observation that the C1 group displayed a markedly decreased level of neutrophil infiltration compared to the C2 group. Furthermore, the C1 group showed a high expression of various protective genes, including ADRB2 and CD81. These results indicate that protective genes may play a crucial role in regulating neutrophils and potentially reducing the severity of AS.

To sum up, our study revealed that the ADRB2 gene might be implicated in the pathogenesis of AS-KD through the regulation of immune cells. However, there are some limitations to our research that need to be addressed. Firstly, further experiments are required to validate some of the findings. Secondly, the molecular mechanisms underlying the interactions between key IRGs need to be further elucidated in both in vitro and in vivo models. Moreover, we plan to conduct additional investigations in the future to address these limitations and expand our understanding of the role of IRGs in the development of AS-KD.

Conclusion

This study employed high-throughput sequencing and bioinformatics methodologies to unravel the intricate pathogenesis of AS-KD at various levels, delving into the role of DEIRGs through functional enrichment and immune infiltration analyses. The application of the LASSO algorithm led to the identification of 12 promising biomarkers, enabling the construction of a robust disease diagnostic model for AS and two immune models. Furthermore, the exploration of the ADRB2-AS-KD relationship through a constructed ceRNA network expanded our understanding. The collective insights garnered from this study could pave the way for potential future treatments for AS-KD, offering a comprehensive perspective on the disease’s underlying mechanisms.