Genome-wide association study in Turkish and Iranian populations identify rare familial Mediterranean fever gene (MEFV) polymorphisms associated with ankylosing spondylitis

Ankylosing spondylitis (AS) is a highly heritable immune-mediated arthritis common in Turkish and Iranian populations. Familial Mediterranean Fever (FMF) is an autosomal recessive autoinflammatory disease most common in people of Mediterranean origin. MEFV, an FMF-associated gene, is also a candidate gene for AS. We aimed to identify AS susceptibility loci and also examine the association between MEFV and AS in Turkish and Iranian cohorts. We performed genome-wide association studies in 1001 Turkish AS patients and 1011 Turkish controls, and 479 Iranian AS patients and 830 Iranian controls. Serum IL-1β, IL-17 and IL-23 cytokine levels were quantified in Turkish samples. An association of major effect was observed with a novel rare coding variant in MEFV in the Turkish cohort (rs61752717, M694V, OR = 5.3, P = 7.63×10−12), Iranian cohort (OR = 2.9, P = 0.042), and combined dataset (OR = 5.1, P = 1.65×10−13). 99.6% of Turkish AS cases, and 96% of those carrying MEFV rs61752717 variants, did not have FMF. In Turkish subjects, the association of rs61752717 was particularly strong in HLA-B27-negative cases (OR = 7.8, P = 8.93×10−15), but also positive in HLA-B27-positive cases (OR = 4.3, P = 7.69×10−8). Serum IL-1β, IL-17 and IL-23 levels were higher in AS cases than controls. Among AS cases, serum IL-1β and IL-23 levels were increased in MEFV 694V carriers compared with non-carriers. Our data suggest that FMF and AS have overlapping aetiopathogenic mechanisms. Functionally important MEFV mutations, such as M694V, lead to dysregulated inflammasome function and excessive IL-1β function. As IL-1 inhibition is effective in FMF, AS cases carrying FMF-associated MEFV variants may benefit from such therapy.


Introduction
Ankylosing spondylitis (AS) is a common form of arthritis affecting primarily the spine and pelvic sacroiliac joints. Twin studies indicate that the disease is highly heritable, with over 90% of the risk of developing the condition being genetically determined [1,2]. To date 114 loci have been identified as being associated with the disease, contributing roughly 30% of the overall risk [3]. The most strongly associated variants at these loci are all common (minor allele frequency (MAF)>5%) or low-frequency (MAF 1-5%) variants, and no rare variant (MAF<1%) has yet been demonstrated to be AS-associated at genome-wide significance (P<5×10 −8 ).
AS is found in most ethnic groups with the notable exceptions of Africans and Australian Aboriginals, in whom the prevalence of HLA-B27 is very low, or as yet unidentified environmental factors protect against disease development. No GWAS has been performed to date in AS cases of Turkish or Iranian descent. These groups are of particular interest because of evidence that patients with the monogenic autoinflammatory disease, Familial Mediterranean Fever (FMF), have an increased prevalence of sacroiliitis [4][5][6][7][8][9], and that FMF-unaffected firstdegree relatives of FMF patients have an increased frequency of AS [5], suggesting an aetiopathogenic link with AS. There have also been four candidate gene case-control association studies of MEFV variants in AS, three of which have demonstrated nominal (10 −5 <P<0.05) association of the main FMF-associated MEFV-variant (rs61752717, M694V) with AS [10][11][12][13], with one marginally negative study (P = 0.065) [14]. All of the participants in these studies were from Turkey. However, the sample sizes of these studies were relatively small (number of patients ranging from 62 to 193) and, therefore, lacked power to achieve definitive significance.
In order to identify new genetic variants associated with AS, and to investigate the potential association of MEFV variants in AS, we performed a genome-wide association study in casecontrol cohorts from Turkey and Iran.

AS-Associated loci in Turkish and Iranian GWAS cohorts
As expected, SNPs in the major histocompatibility complex on chromosome 6p21, and imputed HLA-alleles, were strongly associated with AS in both the Turkish and Iranian cohorts ( Table 1). In the Turkish cohort the strongest SNP associations were with genotyped SNP rs17192932 (OR = 22.84, 95% CI = 17.55-29.72, P = 5.64×10 −120 ), and in the Iranian cohort with rs117486637 (OR = 31.08, 95% CI = 22.27-43.38, P = 8.62×10 −91 ). The strength of association of HLA-B27 with AS was similar in both Turkish and Iranian cohorts, but the proportion of HLA-B27 carriers was lower than in most studies of European descent AS cases (71% and 74%, respectively, compared with 80-95% in most European descent cohorts), consistent with previous reports. Stepwise conditional analyses confirmed risk associations with HLA-B40 subtypes in both cohorts, with HLA-B � 5101 in the Turkish cohort, and with HLA-B � 1503 in the Iranian cohort. HLA-B � 5101 also achieved nominal significance in the Iranian cohort (OR = 1.51, 95%CI = 1.055-2.15, P = 0.024).
Association at genome-wide significance (P<5×10 −8 ) was observed with two non-MHC loci in the Turkish cohort (Table 2), including the known AS-associated locus USP8, and the MEFV locus, and with one known AS-associated locus in the Iranian cohort (ERAP1; Table 3).
No additional loci were identified at genome-wide significance in the meta-analysis combining both cohorts. In the meta-analysis the strongest non-MHC association was with SNPs in the USP8 locus, which has a strong protective effect in AS (exm1161045/rs148783236, OR = 0.58, P = 1.43×10 −13 ). The novel MEFV variant, rs61752717, is also identified at genomewide significance in the meta-analysis. Genome-wide significant association was also replicated with SNPs in chromosome 2p15 (2p15) and ERAP1 (Table 4). In each of these cases, the most strongly associated SNP is in strong linkage disequilibrium (r 2 >0.8) with the previously reported AS-associated variant in European descent populations. Only nominal significance was seen with SNPs at IL23R, a known strongly AS-associated locus, in either dataset or in the meta-analysis analysis. There were 105 out of 113 non-MHC tag SNPs for known ASassociated loci in the combined cohort, while eight of them were missing in the dataset. Two of them (2p15 and ERAP1) were GWS in the combined cohort. One was at suggestive significance (NOS2). Twenty-one of them were at nominal significance. All these 24 tag SNPs were in the same direction of effect as they were in the cross-disease study. Suggestive associations (5×10 −8 <P<1×10 −5 ) were observed at nine novel loci in the Turkish cohort, and three novel and two known AS-associated loci (2p15 and FAS) in the Iranian cohort (Tables 2 and 3, S1-S20 Figs). In the Turkish cohort, near genome-wide significant association was seen with SNPs in ADAM28 (encoding ADAM metallopeptidase domain 28; In the Turkish cohort, the strongest MEFV association observed was with the non-synonymous SNP rs61752717 (OR = 5.34, 95% CI = 3.31-8.62, P = 7.63×10 −12 ), which encodes the M694V protein variant which is the most strongly associated FMF allele. Excluding four cases who have co-existent AS and FMF, this association remains genome-wide significant (OR = 5.26, 95% CI = 3.24-8.54, P = 1.75×10 −11 ). This polymorphism was also nominally associated with AS in the Iranian cohort (OR = 2.85, 95% CI = 1.039-7.83, P = 0.042), its MAF being 3.3-times lower than in the Turkish controls (MAF = 0.0033 in Iranian controls, 0.011 in Turkish controls; Table 5). In the meta-analysis, association was seen with OR = 4.76, 95% CI = 2.94-7.68, P = 1.72×10 −12 . Nominal association was seen with multiple other MEFV SNPs Table 4. Genome-wide significant and suggestive non-MHC loci for AS in the meta-analysis of both Iranian and Turkish cohorts. The loci in italics were in the known 113 loci. in both the Turkish and Iranian cohorts, but none of these are considered FMF-associated, and none were shared between the two cohorts (S1 Table).
The major genetic association with BD is with HLA-B51 [20], which is also AS-associated [21]. rs61752717 was also more strongly associated with AS in HLA-B51-negative than in HLA-B51positive restricted analyses (OR = 5.88, 95% CI = 3.31-10.42, P = 1.36×10 −9 vs OR = 4.82, 95% CI = 1.72-11.66, P = 2.11×10 −3 , or OR = 5.88, 95% CI 3.32-10.43, P = 1.33×10 −9 vs OR = 4.24, 95% CI 1.58-11.42, P = 4.22×10 −3 if excluding four cases with co-existent FMF and AS; S2 Table) but the difference was not statistically significant, nor was interaction observed in a specific logistic regression analysis (P = 0.69). No significant difference was observed in rs61752717 risk allele carriage between HLA-B51-positive and HLA-B51-negative patients (OR = 1.44, 95% CI = 0.23-6.40, P = 0.70) in the Iranian cohort. Similarly, there was no significant interaction identified in the Iranian cohort (P = 0.975; S3 Table) nor in the combined cohort (P = 0.994; S4 Table). Serum IL-1β and IL-23 levels were significantly higher in MEFV 694V variant carriers. MEFV encodes the protein pyrin which functions to control caspase-1-induced IL-1β production. FMF-associated MEFV variants are known to lead to excess IL-1β production, and IL-1 inhibition is effective in treatment of FMF and in a subset of AS patients [22]. Previous studies have also suggested the involvement of Th17 lymphocytes in FMF [23,24]. We therefore compared serum IL-1β, IL-17 and IL-23 cytokine levels in MEFV risk variant carriers in Turkish AS cases and healthy controls. Levels of all three cytokines were significantly elevated in cases compared with controls (Table 9). Among AS cases, serum IL-1β and IL-23 levels were significantly increased in MEFV 694V variant carriers compared with non-carriers (P = 0.0017 and 0.0068, respectively). The result is similar after removing the subjects with coexistent AS and FMF.

Discussion
This study confirms that MEFV is associated with both HLA-B27-positive and-negative AS, even in the absence of clinical evidence of FMF. The association represents the first rare variant association with AS, and has the highest odds ratio for disease of any non-MHC reported locus to date, indicating a major effect on disease pathogenesis. The strength of association observed in the Turkish and Iranian populations was, as expected, closely related to the MAF of the polymorphisms involved, and therefore the absence of association of MEFV with AS in previous studies of other populations may simply due to ethnic differences in gene frequencies rather than differences in the underlying mechanisms driving disease. Further studies examining, for example, gene-expression in patients and healthy controls and in relation to disease activity, will be required to determine whether MEFV-driven autoinflammatory processes are a factor in AS in patients in populations where MEFV functional variants are much rarer. The strength of the association of the M694V MEFV polymorphism with AS is much stronger in HLA-B27-negative than-positive AS. Whilst a formal test of interaction between HLA-B27 and rs61752717 was negative, the stronger association in HLA-B27-negative cases is consistent with previous studies [11], and suggests that this is the third example of epistasis found in AS genetics, following the previously demonstrated examples of ERAP1 variants and HLA-B27 and HLA-B40. The strength of association in HLA-B27-negative subjects is considerable and is to our knowledge the strongest effect in terms of odds ratio of any non-MHC variant in a common immune-mediated disease. In the absence of a definite explanation as to how HLA-B27 induces AS, it is not possible to provide a functional explanation for to this finding. However, one possible hypothesis is that HLA-B27-positive disease is primarily driven by different immunological pathways than HLA-B27-negative MEFV-positive disease.
The replication of the association of HLA-B51 with AS, previously observed in western European-descent AS [21], in the Turkish and Iranian cohorts in the current study also provides further support that AS and Behcet's disease have overlapping aetiopathogenesis. Behçet's disease is also associated with MEFV polymorphisms, and Behçet's disease cases carrying MEFV variants have more severe disease [15][16][17], suggesting that autoinflammation may also contribute to its development. This study also provides further confirmation of the association of HLA-B40 alleles with AS [25,26]. The study sample size was too small to test whether genegene interaction with ERAP1 variants was observed, as has previously been demonstrated in AS cases and controls of western European-descent with both HLA-B27 and HLA-B40 [21].
We also demonstrated that serum IL-1β levels are higher in Turkish AS cases than healthy controls, and are higher in carriers of the MEFV M694V polymorphism. Pyrin, encoded by MEFV, has been shown to negatively or positively regulate caspase-1 and IL-1β activation, depending on the experimental system used [27][28][29]. In a loss of function model, MEFV variants operate by attenuating the apoptosis-associated speck-like protein containing a caspase recruitment domain (ASC, encoded by Pycard) -dependent and ASC independent inhibitory effects of pyrin on caspase activation and subsequent IL-β production [27,28], whereas in a gain-of-function model they lead to caspase-1 activation through an ASC-dependent, NLRP3-independent, pyrin inflammasome [29]. In either model MEFV variants lead to excessive IL-1β production. Pyrin is also a member of the tripartite motif (TRIM) family proteins, which are critically involved in regulation of autophagy and innate immunity. Pyrin/TRIM20, not only acts as a specific receptor for NLRP1, NLRP3 and pro-caspase 1, but also serves as a platform for assembly of key regulators (ULK1, Beclin1, and ATG16L1) and effector factors (mAtg8s) of autophagy initiation machinery [30]. Through these abilities, pyrin leads to degradation of key inflammasome targets in a highly selective manner, resulting in suppression of caspase-1 activation and IL-1β production. Pyrin variants harboring FMF-associated B30.2 mutations, including the M694V variant, have been shown to be associated with a deficiency in the autophagic degradation of NLRP3 [30]. This finding adds more to the complexity and multiplicity of possible mechanisms underlying the excessive IL-1β production seen in patients with FMF and many other autoinflammatory diseases. Pyrin/TRIM20, whose expression is significantly up-regulated by IFN-gamma [29], is one of the required mediators of IFN-gamma-induced autophagy [30], providing further support for pyrin's regulating role in restraining excessive inflammation induced by innate immunity.
In the light of above-mentioned findings, increased M694V prevalence in AS patients observed in our study adds to the mounting evidence that AS is an autoinflammatory disease [31], and also provides support to the recent hypothesis that autophagy may be involved in AS pathogenesis [32]. In the latter context, it is worth mentioning a recent study that showed decreased expressions of autophagy-related genes (MAP1LC3A, BECN1, and ATG5) in the peripheral blood mononuclear cells of AS patients as compared to controls [33]. Notably, even further decreased levels of expression of MAP1LC3A and BECN1 were found in patients with more advanced spinal damage in the same study [33]. Further, there is substantial evidence of activation of the innate immune system in AS, with studies demonstrating that ILC3, MAIT and γδ cells are major sources the pro-inflammatory cytokines IL-17 and IL-22 in the disease [34][35][36], consistent with a role for autoinflammation in AS pathogenesis.
Previous AS genetic studies have demonstrated genome-wide significant associations of the IL-1 receptor genes, IL1R2 and IL1R1 [37], further supporting involvement of this cytokine pathway in the disease. Serum IL-1β levels have been shown previously not to be elevated in AS patients of western European descent [38,39]. Similarly, no difference in serum IL-1β levels were seen in a previous study of Turkish AS cases and controls [40], although in this study no information is reported on the age and gender matching of the cases and controls. Put together, this evidence suggests that IL-1 related processes are important in AS in different ethnic populations, rather than being restricted to MEFV M694V carriers, but that in this group IL-1 plays a major immunopathogenic role. This may explain why minimal association was seen with IL23R variants in this population, perhaps suggesting that the IL-23 cytokine pathway is less important in this group of patients, where IL-1 effects may predominate.
The low efficacy of IL-1 blockade with the IL-1 receptor antagonist anakinra in AS suggests that the disease, at least in the western European population studies, is driven by factors largely independent of IL-1 [22,41]. In contrast, IL-1 blockade is highly effective for FMF caused by the same MEFV polymorphism we demonstrate here to be associated with AS. Our genetic and cytokine level findings support the hypothesis that, at least in patients with the M694V MEFV variant, IL-1 blockade is a potential worthwhile therapeutic option. Anakinra has been shown to be effective in relieving signs and symptoms of spondyloarthritis associated with FMF [42][43][44]. Of particular interest, one of these cases had severe axial pain and elevated acute phase markers persisting despite treatment courses with adalimumab and etanercept, but both the clinical symptoms and acute phase reactant levels improved dramatically with IL-1 blocking therapy [44].
Bacterial pathogens are thought to also be drivers of inflammation in AS, as suggested by microbiome studies demonstrating differences in the gut microbiome in AS cases and controls [45], and the association of the bacterial lipopolysaccharide innate immune receptor TLR4 variants with AS [3]. MEFV variants, and particularly the M694V polymorphism, are also associated with both ulcerative colitis and Crohn's disease in Turkish patients, suggesting a particular role for MEFV in gut mucosal inflammation [46,47]. FMF patients have been shown to have intestinal microbiome dysbiosis, both whilst in remission and during attacks, consistent with models of FMF whereby gut bacterial interactions with a genetically primed host immune system drive disease. At genus level resolution, similarities between what was seen in FMF patients and AS cases are striking, with increases in carriage of Ruminococcus, Porphyromonadaceae, and reductions in Prevotella seen in both diseases [45,48]. Further studies in AS and FMF cases matched for ethnicity and other relevant covariates would appear indicated to determine if specific bacterial species or components are common drivers or have protective effects in these two diseases.
In conclusion, this study demonstrates that the MEFV M694V (rs61752717) SNP markedly increases the risk of AS even in patients not suffering from FMF, and is associated with increased serum IL-1β levels in those patients, suggesting that IL-1 inhibition may be effective in that subset of AS patients, and that MEFV-driven autoinflammation is a factor in the aetiopathogenesis of AS.

Ethics statement
All case and control participants gave written, informed consent. The study was approved by the relevant research ethics authorities at each participating centre, and overall approval was given by Queensland University of Technology Health Research Ethics Committee (approval number HREC/05/QPAH/221).

Study participants
AS was defined according to the modified New York criteria. In total, the study involved 1001 Turkish AS patients, 1011 Turkish controls, 479 Iranian AS patients and 830 Iranian controls. All cases were specifically screened by recruiting clinicians for a history of FMF.

Genotyping
DNA from all subjects were genotyped using the Illumina CoreExome chip following standard protocols at the Australian Translational Genomics Centre, Princess Alexandra Hospital, Brisbane. All Turkish samples were genotyped by CoreExome chips v1.0. 479 Iranian AS cases and 197 Iranian controls were genotyped by CoreExome chips v1.0, while the remaining 633 Iranian controls were genotyped by CoreExome chips v1.1. Bead intensity data was processed and normalized for each sample and genotypes called using the Illumina Genome Studio software.
A subset of 20 carriers of the MEFV M694V had the microarray genotypes tested by Sanger sequencing, with complete concordance of results.

Data quality control
SNP genotypes from CoreExome v1.0 and v1.1 chips from the Iranian cohort were merged based on the manifest files and plink files. 349 SNPs were excluded due to uncertain strandedness or large difference (> 10%) in control allele frequencies between chip versions. We also removed the redundant SNPs with the same minor and major alleles in the same position but with different SNP IDs by keeping the one with lower sample missing rate. Quality control (QC) was completed separately on individual cohorts, and included assessment of missingness by individual (threshold <5%), missingness by genotype (threshold<5%), Hardy-Weinberg equilibrium in controls (P < 1×10 −6 ), extreme heterozygosity (threshold > 3 standard deviations from mean) and identity by descent threshold at the inflation point of PI_HAT (0.2013 and 0.0981 for Turkish cohort and Iranian cohort, respectively). For each pair of related samples (PI_HAT > threshold), the sample with the higher missing rate was removed, and where the pair involved a case and a control with similar call rates (absolute difference in missing rate < 0.0005), cases were preferentially selected for inclusion.
Genotyped SNPs with MAF > 0.05 were then used to perform principal component analysis for ethnicity identification using SHELLFISH (http://www.stats.ox.ac.uk/~davison/ software/shellfish/shellfish.php). Principal components analysis (PCA) was performed with 0-10 eigenvectors. Ethnic and ancestry outliers (more than 6 standard deviations from the mean on any principal component) were excluded. Ski-plots of λ and λ 1000 for different numbers of eigenvectors as covariates for the Iranian, Turkish and combined cohorts are shown in S21 and S22 Figs, and plots of the principal components analysis (first vs second principal component) in S23-S25 Figs. For the Iranian cohort, adding any principal components increased the λ nor λ 1000 . The best λ was produced using no principal components (1.011, λ 1000 = 1.022). In the Turkish cohort, use of the first five principal components produced λ at 1.022 (λ 1000 = 1.024). Use of additional components did not reduce the λ nor λ 1000 further. The quantile-quantile plots for the Iranian, Turkish and combined cohorts are shown in S26-S28 Figs, respectively.
After QC, there were 1,828 Turkish samples (921 cases and 907 controls) and 1,176 Iranian samples (422 cases and 754 controls). After SNP QC, there were total of 294,091 and 293,283 SNPs in the Turkish and Iranian cohorts, respectively.

Imputation
Genotypes were imputed in candidate regions for both cohorts using the 1000 Genomes Project (Phase 1 integrated v3). Genotype data were phased with SHAPEIT [49], and genotypes were imputed with Impute2 [50]. SNPs with low imputation quality (r 2 < 0.6) were excluded. HLA alleles and Amino Acid Polymorphisms in HLA region were imputed by SNP2HLA with T1DGC reference (collected by Type 1 Diabetes Genetics) [51].

Association analysis
Logistic association analyses were performed using PLINK (v1.90b3.36, https://www.coggenomics.org/plink2) [52] to perform with the first five principal components as covariates for population stratification control in the Turkish cohort. Covariates were not added in the Iranian cohort as there was little population stratification (with any covariates the λ increases). Genome-wide significance was accepted at P<5×10 −8 , and genome-wide suggestive at P<1×10 −5 . Association analysis in imputed genotypes was assessed with PLINK best guess genotypes.
Cluster plots for reported SNPs were checked manually in the cases and controls in the respective cohorts and with each different CoreExome chip version separately.

Meta-analysis
The meta-analysis was performed by inverse-variance method implemented in the software package METAL [53].

HLA-B27, HLA-B51, and rs61752717 interactions
Testing for interaction between HLA-B51 and HLA-B27 with the MEFV lead SNP rs61752717 in the corresponding cohort was performed by logistic regression fitting a dominant term for the respective HLA-B allele status (positive or negative) and an additive term for SNP rs61752717, including a multiplicative interaction term, and the N corresponding principal components for population stratification correction (N = 5, 0 and 4 for Turkish, Iranian and combined cohort, respectively): where Phenotype was coded as 2 (patient) and 1 (healthy control), rs61752717 genotype was coded as 0 (genotype CC), 1 (genotype CT) and 2 (genotype TT) to reflect additive effect, HLA-B allele status was coded as 0 or 1 to reflect dominant effect; HLA − Bstatus � rs61752717 genotype was the interaction term for HLA-B and rs61752717 genotype; PCi codes was the ith principal component from PCA from the corresponding cohort.

Serum IL-1β, IL-17 and IL-23 cytokine levels
Concentrations of IL-1β and IL-17A in serum were measured using the Human IL-1β/IL-1F2 Quantikine HS ELISA kit (R&D Systems) and Human IL-17A High Sensitivity ELISA kit (Life Technologies) according to manufacturer's instructions. Statistical significance of differences was evaluated using the Student's t-test. Supporting information S1 Table.