Multiple Biological Measures Reveal Possible Indicators of Bone Growth in Ankylosing Spondylitis

The overall study design is illustrated in Figure 1. In this study, we utilized publicly available GEO and genome-wide association study (GWAS) summary data, with all participants providing informed consent and corresponding ethical approvals from the relevant review boards. Animal experiments were approved by the Ethics Committee of Guangdong Provincial Hospital of Traditional Chinese Medicine (Approval No. 2022025).

The human transcriptome data used in this study were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/↗). On March 25, 2024, we searched the keyword "Ankylosing spondylitis," specifying the organism as "Homo sapiens" and the data type as "expression profiling by array, real-time polymerase chain reaction, and high-throughput sequencing." The sample content was RNA expression data, with no intervention involved.

The cis-eQTL data for AS-related DEGs in this study were obtained from the eQTLGen Consortium (https://eqtlgen.org/cis-eqtls.html↗). The genetic instrumental variables we selected met the following criteria: a p-value less than 5 × 10⁻⁸, minor allele frequency (MAF) greater than 1%, an F-statistic value exceeding 10, and the cis-eQTLs located within 1 Mb upstream or downstream of the gene.

The AS GWAS summary data used in this study were obtained from the R10 release of the FinnGen study (https://www.finngen.fi/en/access_results↗) and the UK Biobank (https://www.ukbiobank.ac.uk/↗). The FinnGen AS dataset includes 3162 AS patients and 294,770 controls, with cases defined using International Classification of Diseases (ICD) codes, including ICD10-M45, ICD9-7200, and ICD8-7124. The UK Biobank AS dataset contains 1296 AS patients and 461,637 controls, with cases primarily defined by self-reported diagnoses.

The phenotypic data used in this study were obtained from the UK Biobank summary statistics in the IEU (Integrative Epidemiology Unit) database, including a total of 3948 phenotypes involving various disease phenotypes, biochemical markers, and health conditions.

The drug information related to AS in this study was obtained from the DrugBank database (https://go.drugbank.com/↗). We searched for "Ankylosing spondylitis" in the disease search section, excluding all veterinary drugs, withdrawn drugs, and drugs without available three-dimensional structures. Additionally, we obtained three-dimensional structures of drugs and targets from the Zinc (https://zinc.docking.org/↗) and PDB (https://www.rcsb.org/↗) databases.

The included AS human RNA expression datasets were batch-combined using the "sva" package in R, and gene symbol annotation was performed using Perl scripts. We then calculated the mean RNA expression levels for both the control and AS groups and performed differential expression analysis. DEGs were selected using p value < 0.05 and |logFC| > 1.5. For the genes, statistical significance was assessed with a two-tailed t-test.

In this study, the "org.Hs.eg.db" package was used to convert DEGs' symbol IDs to Ensembl IDs, and cis-eQTL data for the DEGs were matched from the eQTLGen Consortium as genetic instrumental variables for exposure. The initial number of SNPs obtained from the eQTLGen Consortium was 19,942, filtered using criteria of p value < 5 × 10⁻⁸, MAF > 1%, F‐statistic > 10, and proximity within 1 Mb of the gene. Using the "TwoSampleMR" package, we performed MR analysis between these instrumental variables and AS GWAS datasets provided by FinnGen and the UK Biobank. The GWAS data from FinnGen served as the training set, while the GWAS data from the UK Biobank served as the validation set. To ensure the robustness of the results, we considered the intersection of positive results from both datasets as DEGs with causal relationships to AS. MR analysis needs to meet the following conditions: the SNP-exposure association must be strong, independent of confounders, and only influence the outcome through exposure. The latter two assumptions cannot be empirically verified.

In the MR analysis, the random effects inverse-variance weighted (IVW) method was employed. In instances where the genetic instrumental variable for exposure comprised only a single SNP, the Wald ratio method was utilized. SNP effect alleles were harmonized using the "harmonise_data' function in the TwoSampleMR package. Cochran's Q statistic was applied to assess heterogeneity within the results, with a threshold of p < 0.05 denoting significant heterogeneity. To evaluate horizontal pleiotropy, the MR Egger method was implemented, considering a p value of < 0.05 as indicative of significant horizontal pleiotropy. Results demonstrating significant heterogeneity or horizontal pleiotropy were subsequently excluded from further analysis.

The DEGs that were positive in the MR analysis intersection were then subjected to SMR analysis with the training set to further validate the causal relationship between these DEGs and AS. The HEIDI (heterogeneity in dependent instruments) test was used to detect linkage disequilibrium in the results, and when the p value of the HEIDI test was ≤ 0.01, linkage disequilibrium was considered significant, and the corresponding results were excluded. Finally, we identified a set of DEGs with a clear causal relationship with AS, defined as positive DEGs for AS. In addition, we analyzed the expression differences of these 33 Positive DEGs between the AS group and the control group to determine the expression levels of these genes in AS patients.

We visualized the expression differences of AS Positive DEGs between the control and AS groups to evaluate their regulatory relationships with AS. Subsequently, we used the CIBERSORT tool to estimate the relative proportions of immune cells, summing to one. The expression levels of positive DEGs were then correlated with the relative abundance of immune cells using the Pearson correlation coefficient to obtain correlation coefficients, with the results visualized in a heatmap.

We constructed four predictive models using the expression data of positive DEGs: a random forest (RF) model, a support vector machine (SVM) model, a generalized linear model (GLM), and an extreme gradient boosting (XGBoost) model. These models were chosen for their complementary strengths: RF for handling high-dimensional data, SVM for robustness to outliers, GLM for interpretability, and XGBoost for predictive accuracy. RF ranked genes by mean decrease in accuracy (MDA) from permutation importance; SVM used radial basis function (RBF) kernel coefficients; GLM applied LASSO regularization with λ selected via 10-fold cross-validated AUC optimization; XGBoost employed feature gain metrics. All models underwent 10-fold stratified cross-validation with AUC-ROC as the primary performance metric. After calculating the predictive functions, we selected the most suitable predictive model and identified signature genes from the DEGs through a comprehensive analysis involving residual boxplots and receiver operating characteristic (ROC) curves. Subsequently, a nomogram was created based on the expression of signature genes in both the control and AS groups. To evaluate the accuracy and generalization ability of the nomogram, we generated decision curves and calibration curves.

To evaluate the potential side effects of AS signature genes as potential drug targets, we conducted PheWAS on these signature genes across 3948 phenotypes in the UK Biobank database. This analysis is aimed at exploring the potential beneficial or harmful effects of these signature genes on other diseases. The MR analysis followed the same criteria as previously described.

We obtained AS related therapeutic drugs from DrugBank and performed molecular docking using Autodock Vina v1.2.0. A docking score of less than −5.0 kcal/mmol was considered indicative of a stable docking structure, which allowed us to preliminarily assess the interactions between AS-related therapeutic drugs and signature genes. The Top 4 docking combinations were selected and visualized using PyMOL v2.5.2.

SPF-grade BALB/c male mice, 8 weeks old (body weight 20 ± 2 g), were purchased from Guangzhou Bioenergy Co., Ltd. (Certificate No.44829700012174). The mice were acclimatized for 1 week under standard laboratory conditions. They were randomly divided into a normal group and a model group, with 6 mice per group. The model group was induced to develop an AS model. The specific modeling procedure was as follows: On the first day of the experiment, mice were injected via the tail vein with 2 mg/kg of Arthrogen-CIA 5-Clone Cocktail Kit (Chondrex, 53100). On Day 7, they received an injection of 50 μg/mouse of Biotinylated Lipopolysaccharide (Chondrex, 9028). The control group received an equivalent volume of saline via tail vein injection on the same schedule. All animal experiments complied with the Guidelines for Ethical Review of Laboratory Animal Welfare (GB/T 35892-2018) and were approved by the Ethics Committee of Guangdong Provincial Hospital of Traditional Chinese Medicine (Ethics Approval No. 2022025).

Mice were fasted for 12 h prior to the final treatment. Blood was collected from the eyeball, allowed to stand for 30 min, and centrifuged at 3500 rpm for 20 min to separate the serum. After cervical dislocation, the liver, spleen, and kidneys were weighed. The left hind limb was collected, fixed in 10% paraformaldehyde, subjected to micro-CT analysis, and then decalcified in EDTA solution for 3 weeks. The samples were embedded, sectioned, and subjected to histopathological examination. The right hind limb was evenly divided into two portions: one portion was placed in a cryovial with 1 mL RNA preservation solution, and the other was placed in a separate cryovial. Both vials were stored in liquid nitrogen for 2 h before being transferred to a −80°C freezer for subsequent gene and protein analyses.

The overall condition, fur condition, and activity level of the mice were observed. The weight and joint swelling of each group of mice were recorded weekly. Before the model induction, a mark was made 0.5 cm below the ankle joint, and the width and thickness of the mark were measured using a vernier caliper. After the model induction, the joint swelling was measured at the same location weekly, and the AI score was detected. The detailed rules of the AI score: Joint swelling was scored by two independent blinded investigators using a 0–4 scale: 0 = no redness or swelling, normal appearance of the joints, and free movement; 1 = mild swelling of the joints, possibly limited to a single toe or the little toe joint, with basically unrestricted joint movement; 2 = moderate swelling of the joints, possibly involving multiple toes or metacarpophalangeal joints, with slightly restricted joint movement; 3 = severe swelling of the joints, involving the entire sole or ankle joint, with significantly restricted joint movement; 4 = severe swelling of the joints, even involving the entire foot up to the ankle joint above, with complete loss of joint function, and possible stiffness or deformity.

Joint tissues from mice were homogenized in RIPA lysis buffer to extract proteins, and total protein concentration was determined using a BCA protein assay kit. Protein samples were separated by 8%–12% SDS-PAGE and transferred to a PVDF membrane (Bio-Rad) using wet electrotransfer. The membrane underwent a blocking process in Tris-Tween buffered saline (TBST) with 5% non-fat dry milk for 2 h at room temperature, followed by overnight incubation at 4°C with the following primary antibodies: RUNX2 (CST 12556, 1:1000), OPN (CST 88742, 1:1000), SOX (CST 82630, 1:1000), and β-actin (Proteintech 20536-1-AP, 1:4000). Following three washes with TBST, the membrane was incubated with secondary antibodies for a duration of 2 h at ambient temperature. The visualization of protein bands was achieved through an enhanced chemiluminescence detection system, with ImageJ software used to quantify the grayscale values.

Total RNA was extracted from joint tissue samples using an RNA extraction kit. After assessing the purity and concentration of the RNA, it was reverse transcribed into cDNA and subjected to real-time PCR reactions to detect the mRNA expression of β-actin, RIOK1, FUCA2, COL9A2, USP16, and TTC16. The PCR process began with an initial denaturation phase at 95°C lasting 15 min, followed by 40 cycles consisting of denaturation at 95°C for 10 s, annealing at 58°C for 30 s, and extension at 72°C for 30 s. β-Actin served as the internal control, and the relative expression levels of mRNA were determined using the 2 − ΔΔCt method for relative quantification. The sequences of the primers utilized are provided in Table 1.

Statistical analyses and data visualization were conducted using R Version 4.2. Between-group comparisons of gene expression profiles and immune cell infiltration patterns across clinical subgroups were evaluated through nonparametric Wilcoxon rank-sum tests (two-tailed). Relationship assessments between variables were performed using Spearman's rank correlation coefficient.

We retrieved five datasets (GSE11886, GSE134290, GSE25101, GSE41038, and GSE73754) from the GEO database. By performing batch effect correction, we reduced the batch effects between different datasets (Figure 2). We then conducted differential expression analysis between the AS and control groups using the integrated dataset, identifying 1607 DEGs with significant differences. The gene expression and differential analysis results are provided in Supporting Information 1: File S1.

Based on the differential analysis, 1366 cis-eQTLs were matched for the 1607 DEGs on the eQTLGen Consortium, with instrumental variables provided in Supporting Information 2: File S2. After the Bonferroni correction, the significance threshold for MR analysis was p = 1.83 × 10−⁻⁵ (0.05/2732). The MR analysis results of the DEGs with the AS training and validation sets are provided in Supporting Information 3: File S3 and Supporting Information 4: File S4. After excluding results with significant heterogeneity and horizontal pleiotropy, 248 DEGs were obtained from the training set, and 229 DEGs were obtained from the validation set, with 40 DEGs showing causal relationships with AS after taking the intersection (Supporting Information 5: File S5).

After the Bonferroni correction, the significance threshold for SMR analysis was p = 1.25 × 10⁻³ (0.05/40). The SMR analysis results are presented in Supporting Information 6: File S6, showing that all 40 DEGs had causal associations with AS, but seven DEGs did not pass the HEIDI test (p < 0.01), indicating linkage disequilibrium, and were therefore excluded. Ultimately, 33 DEGs with causal relationships with AS were identified as Positive DEGs.

The differential expression analysis results of the 33 Positive DEGs are shown in Figures 3a,b. The results indicate that ALDH9A1, DPYSL4, TTC16, RIOK1, ZNF74, SNX15, SNAPC5, PRDX3, PRPF4, CD37, COL9A2, POM121, EIF4H, PAICS, SYPL1, ZCCHC9, and RAPH1 were decreased in the AS group, while CEP63, CORO1C, CTH, NBPF3, PDZD8, THBD, WDR25, CYP4F3, FRY, ST6GALNAC4, SEC14L1, FUCA2, PPIF, EXOC6, SENP2, and USP16 were overexpressed in the AS group.

We identified DEGs with causal relationships to AS through MR analysis, followed by validation with SMR, and then analyzed the correlation between the expression of Positive DEGs and immune cell abundance to reveal the potential regulatory relationships between these genes and immune cells. As shown in Figure 4, 28 out of the 33 positive DEGs showed significant regulatory relationships with 19 immune cell types. Among these significant gene-immune cell pairs, 55 showed negative regulatory relationships and 49 showed positive regulatory relationships.

As shown in Figures 5a, 5b, and 5c, among the four machine learning models we constructed, the SVM method demonstrated the highest area under the ROC curve and the lowest residuals. Therefore, we selected this model for further analysis and determined the importance scores of the signature genes (Figure 5d). The Top 5 genes with the highest importance scores used for constructing the nomogram were RIOK1, FUCA2, COL9A2, USP16, and TTC16. Subsequently, we obtained scoring scales for these five signature genes (Figure 5e). By summarizing the expression scores of the signature genes, we assessed the risk probability of AS based on the causally related signature genes. The accuracy of the model is reflected in the decision curve by the distance between the red and gray lines (Figure 5f,g).

The results of the PheWAS are provided in Supporting Information 7: Files S7–S11. After the Bonferroni correction, the significance threshold for the PheWAS analysis was p = 1.27 × 10⁻⁵ (0.05/3948). Based on the PheWAS analysis results of the five AS signature genes, we filtered out all positive disease phenotypes and excluded results with significant heterogeneity and horizontal pleiotropy. Among them, TTC16 showed harmful effects on 27 disease phenotypes across 9 disease categories and beneficial effects on 36 disease phenotypes across 9 disease categories; USP16 had harmful effects on 40 disease phenotypes across 9 disease categories and beneficial effects on 66 disease phenotypes across 10 disease categories; RIOK1 showed harmful effects on 33 disease phenotypes across 8 disease categories and beneficial effects on 31 disease phenotypes across 8 disease categories; FUCA2 had harmful effects on 69 disease phenotypes across 10 disease categories and beneficial effects on 51 disease phenotypes across 11 disease categories; COL9A2 showed harmful effects on 22 disease phenotypes across 7 disease categories and beneficial effects on 24 disease phenotypes across 7 disease categories.

We collected 14 known therapeutic drugs for AS from DrugBank, including Aceclofenac, Acemetacin, Adalimumab, Celecoxib, Etoricoxib, Fenbufen, Flurbiprofen, Ibuprofen, Indomethacin, Ivarmacitinib, Ketorolac, Sulindac, Tofacitinib, and Upadacitinib. To verify whether there are potential interactions between AS therapeutic drugs and AS causal DEGs, we performed molecular docking. The molecular docking results showed that the binding energies of the majority (88.6%) of the docking combinations were below −5.0 kcal/mmol, indicating that stable structures could form, suggesting potential therapeutic effects (Figure 6a).

The top four docking results are presented in Figure 6b. Among them, Adalimumab primarily interacts with USP16 through hydrogen bonds with residues ASP866, ARG21, LEU243, GLN239, and PHE209, resulting in a strong binding energy of −9.0 kcal/mmol. Fenbufen interacts with USP16 via hydrogen bonds involving residues TYR870, TYR174, MET177, LEU181, and ASP866, with a binding energy of −8.5 kcal/mmol. Indomethacin docks with USP16 via a hydrophobic pocket, resulting in a binding energy of −8.5 kcal/mmol, involving residues such as LEU181, PHE209, LEU867, and LEU858. Additionally, Ketorolac can bind to TTC16 via hydrogen bonds involving residues TYR759, HIS750, PHE700, LEU581, and LEU585, with a binding energy of −8.4 kcal/mmol.

As shown in Figure 7, 3 days after modeling, the AS mice exhibited reduced body weight, paw swelling, lethargy, and decreased food and water intake compared with the control group. Subsequently, joint swelling progressively worsened, peaking on Day 9, accompanied by reduced mobility. Thereafter, body weight gradually increased, while joint swelling significantly decreased (p < 0.01). Compared with the control group, the AI score in the model group was significantly elevated. However, as inflammation naturally subsided, joint swelling gradually alleviated (p < 0.01).

As shown in Figure 8, micro-CT analysis revealed ectopic new bone formation at the tendon–bone junctions of the tarsal and phalangeal bones in the model group, along with bone erosion along the inner bone marrow cavities, compared with the control group. As shown in Figure 9, H&E staining further demonstrated thickening of the tarsal and metatarsal bones, as well as their articular surfaces. Safranin O-fast green staining indicated that the original ligament tissue attached to the bone and joint surfaces had differentiated into red-stained osseous and cartilaginous tissues, suggesting pathological ossification.

As shown in Figure 10, Western blot analysis revealed that the protein expression levels of RUNX2, OPN, and SOX9 in the left hind paw tissues of the model group were significantly higher than those in the control group (p < 0.01), indicating an active pathological ossification process.

As shown in Figure 11, compared with the control group, the mRNA levels of FUCA2 and USP16 in the model group were significantly increased (p < 0.01), and the mRNA level of TTC16 was decreased (p < 0.05) in the model group. However, there was no significant difference in the expression of COL9A2 and RIOK1. We further conducted WB analysis to detect the protein expression levels of FUCA2, TTC16, and USP16. It was found that compared with the normal group, the protein expression levels of FUCA2 and USP16 in the CAIA model group were increased (p < 0.01), while the protein expression level of TTC16 was decreased (p < 0.01).