In this study, genetic variants associated with protein abundance are differentiated into cis-acting and trans-acting types. The quantitative trait loci that act on CS (pQTLs) were typically located near the genes that encode the corresponding proteins (within 1Mb in this work) and have direct and independent biological effects (11). Proteins are more likely than other molecular traits to be used as drug targets, and MR analysis combined with pQTL as instrumental variables (IVs) is valuable for drug development in human genetics. We integrated the cis-pQTL data derived from prior studies conducted by Sun BB et al. (12–17) as our genetic IVs source. Significance thresholds were set at P< 5e-8 (with a minor allele frequency > 0.01) and a linkage disequilibrium (LD) threshold of R ²< 0.8, kb = 10,000, utilizing the two-sample MR package (version 0.5.7) R (18). The explanatory power of IVs concerning exposure was quantified using the F-statistic as shown below:

In this formula, MAF represents the minimum allele frequency; β represents the effect size in GWAS for proteins and sd is its standard deviation (sd = se × N ); N is the total number of samples in the proteins GWAS; k is the number of IVs used (19). An F-statistic below 10 suggests a potentially weak IV strength, indicative of inadequate suitability for robust causal inference. Blood pQTLs associated with AS were identified from the latest public GWAS database as outcomes (https://storage.googleapis.com/finngen-public-data-r9/summary_stats/finngen_R9_M13_ANKYLOSPON.gz). We employed blood proteomic profiling data as SNPs to perform a two-sample MR analysis across two independent cohorts consisting of 2,860 AS cases and 270,964 controls. This methodological approach facilitates the identification of SNPs linked to specific protein expressions, advancing research into genetic variations associated with diseases. The workflow of the study is depicted in Figure 1 .

MR analysis mandates the fulfillment of three crucial assumptions: genetic instruments are strongly correlated with exposure, genetic instruments are independent of confounding factors, and genetic instruments solely influence the outcome through their impact on exposure. The Wald ratio (WR) method was employed for preliminary explorations of the top SNPs. (most significant SNPs) (20, 21). Where only a single IV was available with a P-value of<1e-05, the WR was computed. A Bonferroni correction was applied to the P-value (0.05/1,654 = 0.00003023) to mitigate the false discovery rate (FDR), and findings were articulated using statistical metrics (Odds Ratio and 95% Confidence Interval) (22). Throughout these analyses, only results with a statistically significant level (P< 0.05) were retained to ensure the reliability of the results.

Colocalization analysis serves as a methodological approach for investigating protein-protein interactions and functions (23). When comparing the subcellular localization patterns of proteins, we reveal their interrelationships. This analysis facilitates the assessment of whether disease-associated and genetic regulation of protein levels are driven by identical genetic variants. In our study, we used Bayesian colocalization analysis with default parameters from the “coloc” software package to assess the likelihood of two features sharing the same causal variation (24). We evaluated five hypotheses’ posterior probabilities (PPH) to ascertain whether two features are influenced by a common variant. The coloc.abf algorithm was applied, classifying proteins as first-tier (PPH4 > 0.8) with strong evidence of colocalization, second-tier (0.5< PPH4< 0.8) with moderate support, and others as third-tier targets based on prior research (25).

The exploration of protein-protein interactions (PPI) was facilitated using the String database (https://string-db.org/) (26). It contains known protein interaction data from various sources, including experimental evidence and computational predictions. This step involved data preparation, access to the String database, extraction of the PPI network, and subsequent visualization and analytical interpretation of these interactions. The parameters for this study were configured as follows: the network edges were specifically interpreted as indicators of interaction evidence, while the active interaction sources included a range of methodologies such as textmining, experimental assays, Database entries, co-expression patterns, proximity in genomic context, Gene Fusion events, and co-occurrence statistics. The threshold for inclusion in the network was set at a medium confidence score of 0.400 to ensure relevance and reliability of the interactions.

Regarding the identification of candidate genes implicated in AS pathogenesis, genes derived from the Bonferroni correction statistical approach were pinpointed as key targets. Further insights into the biological pathways and functional attributions of these genes were gleaned through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. These analyses were performed using Metascape (http://metascape.org) and the Oebiotech platform (https://cloud.oebiotech.cn/, last accessed on 12 January 2024) (27). These tools facilitated the visualization and interpretation of potential pathways involved in the disease’s progression and pathogenesis, enhancing our understanding of the molecular mechanisms at play (27, 28).

Phenome-wide mendelian randomization (PhenMR), integrating principles from both Mendelian Randomization (MR) and Phenome-Wide Association Studies (PheWAS), is employed to elucidate the causal relationships between genes and a spectrum of phenotypic traits (29). This method enables a broad assessment of genetic impacts across health and disease phenotypes, revealing gene-phenotype associations and facilitating the exploration of genetic functions and biological underpinnings.

To exclude potential pleiotropy of AS, we performed a PheWAS analysis using variants identified from the GWAS ATLAS and Pheno-Scanner databases (30–32). This analysis allows for the investigation of how the six primary proteins linked to AS are associated with various traits, shedding light on their roles in AS pathogenesis and their potential as therapeutic targets.

The drug gene interaction database (DGIdb) (https://dgidb.genome.wustl.edu) is a database that integrates drug-gene interaction data from 30 existing databases, helping researchers to explore existing data to generate hypotheses about how genes can be used in therapy or drug development (33). It allows users to search for drugs that target specific genes of interest, including FDA-approved, experimental, and investigational drugs, along with information on their indications, target genes, and interaction types. Our search within the DGIdb was guided by prior research, focusing on genes associated with potential treatment modalities for AS (34). Ethical approval had been obtained from their institutional review board in each study, and all participants had provided informed consent, obviating the need for additional ethical approval.

The three-dimensional configuration of the MAPK14 protein was sourced from the Protein Data Bank (https://www1.rcsb.org, PDB ID: 6zwp) and refined by removing non-essential entities using PyMOL software, version 2.5.4 (35, 36). The molecular structures of four therapeutic agents, as identified in the Drug Gene Interaction Database (DGIdb), along with a co-crystallized ligand, were retrieved from the PubChem repository, (https://pubchem.ncbi.nlm.nih.gov) and subjected to structural optimization employing OpenBabel (version 3.1.1) and ORCA (version 5.03) (37–39). The root mean square deviation (RMSD) values for co-crystallized ligand structures were calculated before and after docking using PyMOL software. An RMSD value of less than 2 Å indicates successful methodological validation. Computational assessments were conducted using the B3LYP hybrid functional combined with the def2-TZVP basis set. To prepare for molecular docking, modifications to enhance hydrogen bonding and molecular flexibility were executed using Auto Dock Tools version 1.5.6. The docking procedure itself was facilitated by Autodock Vina software (version 1.2.0), employing the Lamarckian genetic algorithm to optimize interaction predictions (40). Specific parameters for the docking grid were set with a central point at coordinates (X, Y, Z) = (6, 0, 19) and dimensions spanning (X×Y×Z) = (26×26×24), with an exhaustiveness setting of 25 to ensure thorough sampling.

This study explored how certain proteins, for which genetic markers are known, are related to the risk of AS. We focused on 1,949 proteins that have identifiable genetic signals known as pQTLs and examined their influence on the development and severity of AS. Figure 2 provides a summary of these findings. We match a single cis-SNP for each protein by integrating genetic data from pQTL, and we identify 1,654 unique circulating plasma proteins that met the criteria for having a causal relationship with AS ( Supplementary Table S1 ). Among these, 868 proteins showed a positive causal relationship with AS, while 786 proteins showed a negative causal relationship with AS. This suggested that specific genetic variations may be associated with the risk of AS by influencing the expression or function of these plasma proteins, thus affecting the occurrence of AS.

To ensure the reliability of our findings, we used the Bonferroni correction method to adjust for multiple comparisons, setting a strict significance threshold (P value = 0.05/1,654, P< 3.03E-05). Among the 1,654 plasma proteins mentioned above, we identified 18 that were significantly associated with the risk of developing AS. These proteins included AGER, AIF1, ATF6B, C4A, CFB, CLIC1, COL11A2, ERAP1, HLA-DQA2, HSPA1L, IL23R, LILRB3, MAPK14, MICA, MICB, MPIG6B, TNXB and VARS1. Specifically, seven of these proteins were associated with an increased risk of AS, including ATF6B, C4A, COL11A2, etc. Additionally, 11 proteins were found to lower the risk of AS, such as AIF1, MAPK14, MICA, MICB, etc. Interestingly, the most significant factors influencing AS were the proteins IL23R (OR: 2.49), ATF6B (OR: 2.21), and C4A (OR: 1.38). These findings provide compelling information for further research on the role of these proteins in AS, warranting further investigation ( Table 1 ).

We conducted a Gene Ontology (GO) analysis on 18 shared protein targets, covering biological processes (BP), Molecular Function (MF) and Cellular Component (CC), as seen in Figure 3A . Notable GO terms included positive regulation of interleukin-12 production (GO:0032735), regulation of T cell-mediated cytotoxicity (GO:0001914). The KEGG enrichment analysis showed that Th17 cell differentiation signaling pathway was significantly enriched in target proteins (see Figure 3B ). We also constructed a map showing interactions between these proteins, featuring 38 connection points and 82 connecting lines, which suggests a highly significant network based on a statistical measure (P-value:1.45e-12), as seen in Figure 3C using the String platform.

In our research, we performed a detailed analysis to study how causal genetic variations and proteins were connected to AS. This technique, known as co-localization analysis, helps us understand which genetic changes are likely to influence the disease and its treatment. Specifically, we found that four proteins—MAPK14, AIF1, ATF6B, and MICA— can serve as first-tier drug targets. HSPA1L and COL11A2 were second tiers, according to the results, while others were categorized as third-tier targets. Our findings suggested possible shared causal genetic variations between MAPK14 and AS (PPH4 = 0.977), while ATF6B could pose a risk factor for AS (PPH4 = 0.887). These results underscore our identification of several credible drug-related genes via MR and colocalization analysis, revealing a common genetic effect between target proteins and AS risk ( Table 2 ).

We sourced data on gene-trait associations from the GWAS Catalog accessible at: https://gwas.mrcieu.ac.uk/. Our analysis focused on examining the link between gene variants and traits. We utilized cis-pQTLs to understand these relationships and conducted thorough searches to pinpoint studies that reported statistically significant genetic associations (P< 0.001) with specific traits. Our study incorporated a detailed PheWAS analysis to investigate how the top six proteins that colocalize—HSPA1L, MICA, COL11A2, MAPK14, AIF1, and ATF6B were connected to We discovered that these proteins not only contribute to targeting AS but also significantly impact the development of other diseases. For example, AIF1, MICA, ATF6B, and HSPA1L were found to be associated with an increased risk of rheumatoid arthritis (RA), while AIF1 was also associated with an increased risk of atopic dermatitis. On the other hand, MAPK14, ATF6B, and COL11A2 were associated with a lower risk of both rheumatoid arthritis and primary sclerosing cholangitis (refer to Supplementary Table S2 ). exploration enables us to deepen our understanding of how specific genetic variations can influence disease mechanisms, providing critical insights for genetic research and biomedical applications beyond just treating AS.

Our research identified MAPK14 as a potential genetic drug target and biomarker in AS. By analyzing the DGIdb database, we discovered a link between MAPK14 and the drug pirfenidone, which is currently under clinical investigation. Interestingly, MAPK14 is also associated with three FDA-approved drugs: dasatinib, doxorubicin, and sorafenib. Although our study did not directly involve pharmaceutical action on AS, these findings provided new prospects for the treatment of AS and related diseases ( Supplementary Table S3 ).

In the realm of molecular docking (MD) investigations, the concept of binding energy serves as a fundamental metric for assessing ligand-receptor interactions. A binding energy less than 0 kcal/mol is indicative of a spontaneous binding process, confirming the thermodynamic feasibility of the ligand associating with the receptor. Specifically, a binding energy threshold of -7.2 kcal/mol is often utilized as a benchmark for strong molecular interactions, characterized by elevated affinity and specificity towards the receptor. In our analysis, we evaluated the docking profiles of pirfenidone, dasatinib, doxorubicin, and sorafenib with the kinase domain of MAPK14. The results demonstrated binding energies of -9.7 kcal/mol, -9.0 kcal/mol, -8.5 kcal/mol, and -12 kcal/mol, respectively. All values significantly surpass the -7.2 kcal/mol threshold, suggesting robust binding affinities (see Figure 4 , Supplementary Table S4 for details). Concurrently, we employed the methodology of molecular docking to validate the co-crystallized ligand with MAPK14 ( Supplementary Figure S1 ). The RMSD (Root Mean Square Deviation) value between the co-crystal ligand and MAPK14 before and after docking is 0.321Å, which is less than 2Å, indicating that the molecular docking methodology has been validated. These findings imply a high potential for these compounds to inhibit MAPK14 activity effectively, which may be therapeutically beneficial in modulating pathways regulated by this kinase in AS ( Figure 5 ).