The primary cohort of this study included 33,998 participants enrolled from five genome-wide association study (GWAS) datasets on ankylosing spondylitis, involving individual-level data from 4,913 participants of Chinese descent and summary-level data from 2,012 participants of a Turkish cohort, 1309 of an Iranian cohort and a collection consisted of 22,647 individuals of European ancestry and 3,117 of East Asian ancestry enrolled in International Genetics of Ankylosing Spondylitis Consortium (Table 1) (International Genetics of Ankylosing Spondylitis Consortium et al., 2013; Li et al., 2019). The 1984 modified New York criteria were used to determine the outcome of ankylosing spondylitis in all the above studies (Linden et al., 1984). We also used the datasets of GWAS meta-analyses on coronary artery disease involving 86,995 individuals and on type II diabetes involving 659,316 individuals both of European ancestry in exploratory analyses (Schunkert et al., 2011; Xue et al., 2018). We obtained summary datasets reported by the Global Lipids Genetics Consortium on LDL cholesterol levels based on 188,578 participants for the exposure assessment in this study (Willer et al., 2013). The covariate and multivariate analyses in our study additionally included four GWAS datasets involving a total of 2,509,298 individuals, that included data on one of the following traits: 25-hydroxyvitamin D levels (Jiang et al., 2018), smoking initiation (Liu et al., 2019; Matoba et al., 2019), drinking per week (Liu et al., 2019), and overall physical activity time (Doherty et al., 2018) (Supplementary Table 1). Additional details for each included study are provided in the Supplementary Appendix. All contributing studies were approved by their local research ethics authorities, and written informed consent was obtained from all participants.

To measure the magnitude of genetically proxied inhibition of HMG-CoA reductase (HMGCR), we constructed a genetic score including single nucleotide polymorphisms (SNPs) within a ± 100 kb window from HMGCR that were significantly associated with the LDL cholesterol level (P value < 5.0 × 10^–8). A second inclusion criterion included the fact that all SNPs of the score were only in a weak pairwise linkage disequilibrium (r² ≤ 0.2) in both European as well as East Asian populations. Four SNPs (rs12916, rs17648288, rs3857388 and rs10064936) that met these criteria were included in the HMGCR score as the primary exposure in this study to proxy for the inhibition of HMG-CoA reductase. The exposure alleles for each SNP were as those in association with a lower LDL cholesterol level, and their per allele effect on LDL cholesterol decrement was obtained from the GWAS meta-analysis as reported by the Global Lipids Genetics Consortium (Willer et al., 2013). For each individual, exposure of the HMGCR score was equivalent to the value by summing the number of exposure alleles inherited at each SNP multiplied by its corresponding conditional effect on the LDL cholesterol level. Using the similar criteria, we constructed the PCSK9 and NPC1L1 genetic scores composed of weakly correlated SNPs (r² ≤ 0.2) that were associated with LDL cholesterol levels (P value < 5.0 × 10^–8) and located ± 100 kb from PCSK9 and NPC1L1, respectively. In addition, we constructed the LDL polygenetic score by combining independently inherited SNPs (r² < 0.001) from different gene regions in association with LDL cholesterol levels (P value < 5.0 × 10^–8), as a natural indicator to proxy for lifelong stable LDL cholesterol levels. Differences in the genetically proxied inhibition of PCSK9 and NPC1L1 as evaluated by PCSK9 and NPC1L1 scores and genetically determined serum LDL cholesterol levels by the LDL score were defined as secondary exposure in this study.

The objective of this study is to evaluate the association of genetic scores with the risk of ankylosing spondylitis. Because inheritance of an exposure allele at each genetic variant is approximately random at birth and can exert a long-term effect throughout the lifetime of an individual, exposure differences in a combined genetic score by virtue of low linkage disequilibrium among each variant would also be naturally randomized (Palmer et al., 2012). Therefore, according to a paradigm termed as Mendelian randomization (Supplementary Figure 1), the genetic score variable might have an independent effect on events that occurred incident to the genetic exposure from the perspective of the time of birth, which is analogous to the intention-to-treat effect in a randomized controlled trial (Burgess et al., 2018). In this study, the primary outcome was the odds ratio (OR) for ankylosing spondylitis, estimated by an association with the HMGCR score. Secondary outcomes were the OR for ankylosing spondylitis associated with the PCSK9, NPC1L1 and LDL genetic scores. Exploratory outcomes included the OR for coronary artery disease and type 2 diabetes as the positive control outcomes, for which a causal effect had already been established in previous studies explained by statin therapy or genetically proxied inhibition of HMG-CoA reductase (Cholesterol Treatment Trialists’ Collaboration et al., 2010; Ference et al., 2015; Swerdlow et al., 2015). We obtained effect estimates for individual SNPs in each genetic score by performing a meta-analysis of summary data across available GWAS datasets. Data for indirectly genotyped SNPs were obtained from imputation or from a suitable proxy SNP in strong linkage disequilibrium (r² ≥ 0.9) with reference to 1000 Genomes Phase 3 panel in the same population (1000 Genomes Project Consortium et al., 2015). The overall estimates for genetic scores were produced by combining the effect of individual SNPs using summary-level data in an inverse-variance weighted model that had taken into account weak correlation caused by the low pairwise linkage disequilibrium among variants included in genetic scores, which can be largely equivalent to effect estimates using individual-level data (Burgess et al., 2016).

F statistics were calculated to examine the strength of association between the genetic variants and the exposure factor (Palmer et al., 2012). An F statistic greater than 10 indicates a valid genetic instrument having sufficient relevance with the exposure factor (Palmer et al., 2012). The statistical power was calculated as previously described (Burgess, 2014).

With the summary-level data, we obtained the proportional effects of variants on the outcome relative to their effects on LDL cholesterol reductions using the Wald type ratio method (Burgess et al., 2013). Given the fact that a few variants were selected from the same gene region in the HMGCR, PCSK9, and NPC1L1 scores, a fixed-effects inverse-variance weighted model was utilized to combine individual SNP effects into these weighted scores, with adjustment for the weak correlations between variants (Burgess et al., 2016). In terms of potential heterogeneity from multiple LDL-cholesterol-lowering variants, a multiplicative random-effects inverse-variance weighted model was utilized for the LDL polygenetic score (Burgess et al., 2013). With the individual-level data, we evaluated the dose-response relationship between the HMGCR score and the odds ratio for ankylosing spondylitis using restricted cubic splines by remodeling the logistic regression equation. All estimates were described as odds ratios (ORs) or mean effect sizes (log OR) with 95% confidence intervals (CIs) standardized for per mmol/L decrement in LDL cholesterol levels.

For inverse-variance-weighted estimates, heterogeneity was quantified using Cochran’s Q and I² statistics. Several sensitivity analyses were performed to examine the presence of pleiotropy meaning that genetic variants may influence the outcome via other pathways and not only dependent on the drug target or exposure factor of interest. For the genetic score of drug targets, we examined whether variants included in the score were associated with potential factors that had been known to affect the risk or disease activity of ankylosing spondylitis. The major risk factors previously reported for ankylosing spondylitis include male gender (Lin and Gong, 2017), HLA-B27 antigen (Lin and Gong, 2017), serum 25-hydroxyvitamin D levels (Cai et al., 2015), smoking (Videm et al., 2014), drinking (Min et al., 2019), and physical activity (O’Dwyer et al., 2017). A detectable significant association (P value < 0.05) between the genetic score and any of these risk factors indicated the presence of an underlying horizontal pleiotropy. Therefore, a multivariable analysis was performed to determine the independent effect of the genetic score on ankylosing spondylitis with adjustment for the effects of those significantly associated traits as covariates. For the LDL score composed of polygenetic variants, the MR Egger regression was performed to detect the unbalanced pleiotropy among multiple variants (Bowden et al., 2015). The presence of directional pleiotropy would be inferred, provided that an intercept in the regression equation was statistically different from null (zero) (P value < 0.05). In addition, sensitivity analyses including maximum likelihood, weighted median and MR-PRESSO outlier-corrected models were also conducted to test the robustness of the results, under different statistical assumptions addressing the bias and pleiotropy (Burgess et al., 2013; Bowden et al., 2016; Verbanck et al., 2018).

To provide complementary evidence supportive of the linkage of two traits through shared causal variants, we performed the Bayesian statistical framework colocalization analyses to prioritize genetic variants located ± 100 kb from HMGCR that may influence both the molecular trait of HMG-CoA reductase and the risk of ankylosing spondylitis (Supplementary Figure 1) (Giambartolomei et al., 2018). The per SNP posterior probability (PPA) as a shared variant between the two traits was calculated to assess whether the association between two traits were driven by different causal variants correlated due to linkage disequilibrium. In addition, given the fact that the GWAS meta-analysis by the Global Lipids Genetics Consortium was generated primarily in a European population, we performed subgroup analyses for association with ankylosing spondylitis only involving participants of European ancestry. This analysis was used to examine potential bias caused by population stratification, in which unevenly distributed genotypes across geographical regions may cause spurious association with both the ancestry and the outcome. Given the fact that the functional evidence relating genetic variants to the expression of HMGCR is currently lacking with the exception of rs12916 (Swerdlow et al., 2015; Würtz et al., 2016), we performed a sensitivity analysis only including rs12916 and removing the other variants from the HMGCR score to test its association with the risk of ankylosing spondylitis.

As the primary outcome for the HMGCR score, a significance threshold was set at P value < 0.05. For secondary outcomes, Bonferroni correction was applied for multiple testing, whereby a P value < 0.017 [0.05/3 tests (PCSK9, NPC1L1 and LDL scores)] was considered to indicate statistical significance. There was no plan on adjustment for multiple comparisons for other outcomes or in sensitivity models, and therefore these results should be interpreted as exploratory. All statistical tests were conducted as two-sided with the use of Stata version 15.0 (StataCorp, College Station, TX, United States) and R version 3.6.3 (R Foundation, St. Louis, MO, United States). For extended methods, detailed descriptions are provided in the Supplementary Appendix.

Of 33,998 participants in the primary cohort, 12,596 individuals had been diagnosed with ankylosing spondylitis and 21,402 were unaffected participants (Table 1). Up to 22,233 cases with coronary artery disease and 64,762 controls as well as 62,892 type II diabetes cases and 596,424 controls from two meta-analyses of GWAS were included in the exploratory analysis for the outcome on coronary artery disease and type II diabetes, respectively. The GWAS datasets and sample size used for the covariate and multivariate analyses are summarized in Supplementary Table 1.

Four SNPs were included in the HMGCR score, eight SNPs in the PCSK9 score, two SNPs in the NPC1L1 score, and 67 SNPs in the LDL polygenetic score (Table 2 and Supplementary Table 2). Each of the chosen SNPs as well as the combined scores were strongly associated with the trait of serum LDL cholesterol levels with all F statistics greater than 10.0, respectively (Supplementary Table 3). This indicated that these scores could be used as an appropriate instrument as the genetically proxied drug target inhibition or lifelong LDL cholesterol levels.

For the primary outcome, assuming that the per mmol/L decrement in LDL cholesterol by the HMGCR score resulted in a twofold decrease in the odds of ankylosing spondylitis, the current sample size of the primary cohort would have a 97.4% power to obtain an estimate with a two-sided significance level of 0.05 (Supplementary Figure 2).

In the inverse-variance-weighted model, genetically proxied inhibition of HMG-CoA reductase equivalent to a 1 mmol/L decrement in LDL cholesterol levels by the HMGCR score was associated with a decreased risk of ankylosing spondylitis (OR, 0.57; 95% CI, 0.38–0.85; P value = 5.7×10^–3) (Figure 1). No statistically significant heterogeneity (P value for Q statistic = 0.121) was detected in the estimate (Supplementary Table 4). In a dose-response analysis using the individual-level data, an overall trend for a progressive reduction in risk of ankylosing spondylitis was observed when given a stepwise increase in the number of inherited exposure alleles (P value for trend = 0.020) or in the value of the HMGCR score (P value for trend = 0.046) (Figure 2). This finding suggested a dose-dependent association between genetically proxied inhibition of HMG-CoA reductase and the risk of ankylosing spondylitis.

No significant associations with ankylosing spondylitis were observed for the PCSK9 score (OR, 0.89; 95% CI, 0.68–1.16; P value = 0.378) and the NPC1L1 score (OR, 1.50; 95% CI, 0.39–5.77; P value = 0.556) (Figure 1). For the LDL score, the genetically determined lower LDL cholesterol level was associated with a reduced risk of ankylosing spondylitis, whereby the odds ratio was 0.64 (95% CI, 0.43–0.94; P value = 0.024) standardized for per mmol/L decrease in LDL cholesterol levels (Figure 3). However, when accounting for multiple testing across secondary outcomes by the Bonferroni correction, the statistical significance of the association was lost [P value > 0.017 (0.05/3 tests)]. In addition, a detectable heterogeneity was presented across estimates for multiple SNPs included in the LDL score (P value for Q statistic = 3.7 × 10^–8) (Supplementary Table 4), indicative of the underlying presence of directional pleiotropy (Supplementary Figure 3). The MR Egger intercept test also demonstrated an unbalanced pleiotropy in the effect estimate for the LDL polygenetic score (MR Egger intercept = −0.035; P value = 0.036) (Supplementary Table 4).

The male gender related to Chromosome Y and the HLA-B27 antigen encoded on Chromosome 6 have been reported as strong risk factors for ankylosing spondylitis, but are independent of HMGCR genetic variants located on Chromosome 5. Analyses of a subgroup of HLA-B27-antigen-tested participants also showed no association of HMGCR score with HLA-B27 antigen positivity or male gender (Supplementary Table 5). When examining the association with other known risk factors for ankylosing spondylitis, per change in HMGCR score was found to be significantly associated with a decrease in serum 25-hydroxyvitamin D levels and an increase incidence in smoking, indicating that genetic variants in the HMGCR score were likely to affect the outcome of ankylosing spondylitis via other pathways by influencing these two traits (Table 3). After adjusting the observed association for 25-hydroxyvitamin D level and smoking in the multivariate analyses, the odds ratio for ankylosing spondylitis per change in HMGCR score scaled to a 1 mmol/L decrement in LDL cholesterol levels was 0.12 (95% CI, 0.04–0.39) (Table 3). Similar results between the primary analysis and this sensitivity analysis demonstrated the absence of substantial bias caused by the putative horizontal pleiotropy.

For the LDL polygenetic score, several different statistical models were used to adjust for potential bias and pleiotropy (Supplementary Table 6). Directionally consistent estimates on the risk of ankylosing spondylitis were obtained by maximum likelihood (OR, 0.63; 95% CI, 0.49–0.83), MR Egger regression bootstrap (OR, 0.78; 95% CI, 0.47–1.29), weighted median (OR, 0.98; 95% CI, 0.64–1.49) and MR-PRESSO outlier-corrected (OR, 0.75; 95% CI, 0.54–1.03) models, but with wider confidence intervals.

Evidence of colocalization across variants nearby HMGCR proxying HMG-CoA reductase inhibition and ankylosing spondylitis with a posterior probability (PPA) > 0.01 was detected for three SNPs, rs7703051, rs11749783 and rs3846663, of which all are in strong linkage disequilibrium with each other (r² = 1.0) (Supplementary Figure 4). The top-colocalized SNP rs7703051 (PPA = 0.515) is also in high linkage disequilibrium with a variant (rs12916) included in the HMGCR score (r² = 0.98). The effect estimate for rs7703051 was consistent with the primary analysis in an odds ratio of 0.62 (95% CI, 0.40–0.98) for ankylosing spondylitis proportional to per mmol/L decrement in LDL cholesterol levels. These findings implied that there might be a biological mechanism linking HMG-CoA reductase inhibition and ankylosing spondylitis via a shared causal locus within HMGCR. The observed association between the two traits was thus less likely to be driven by different gene regions.

Furthermore, subgroup analyses only including the European population provided a similar estimate for the HMGCR score (OR, 0.56; 95% CI, 0.36–0.85) as well as other genetic scores, as compared with the results from the full primary cohort (Supplementary Table 7). Therefore, these data suggested that when combining different populations to maximize the sample size and statistical power, the primary analysis did not suffer from substantial bias influenced by population stratification.

There was functional evidence relating rs12916 to the mRNA expression of HMGCR, supporting its validity as genetic instruments in our analysis (Swerdlow et al., 2015; Würtz et al., 2016). However, despite the strong association of the selected SNPs with the LDL cholesterol level (P value<5.0×10^–8), the biological function of the other genetic variants included in HMGCR score is currently unclear. A sensitivity analysis leaving out the other SNPs and only including rs12916 in the HMGCR score showed a consistent association with the risk of ankylosing spondylitis (OR, 0.58; 95% CI, 0.38–0.88) as compared with the primary analysis.

For exploratory outcomes, inverse-variance-weighted analyses showed a decreased risk for coronary artery disease as well as an increased risk for type 2 diabetes especially associated with genetic variants related to genetically proxied inhibition of HMG-CoA reductase (Supplementary Table 8). These results were in line with expectation based on the protective role of statins on coronary artery disease and a higher incidence of type II diabetes in statin-treated patients (Cholesterol Treatment Trialists’ Collaboration et al., 2010; Ference et al., 2015; Swerdlow et al., 2015), which suggested that genetic variants included in the HMGCR score can reliably assess the effect of HMG-CoA reductase inhibition.