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-9], and that FMF-unaffected first-degree 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−5FMF-associated MEFV-variant (rs61752717, M694V) with AS [10-Joint Bone Spine. 2010 ">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 case-control cohorts from Turkey and Iran.

Results

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).

Table 1

Stepwise conditional logistic regression analysis of HLA-B allelic associations with ankylosing spondylitis in the Turkish and Iranian cohorts.

Turkish Cohort
Round	HLA-B allele	OR (95% CI)	Stepwise conditional P-value	InitialP-value	Allele frequency(case/control)
1	*2705	34.31 (25.33–46.48)	1.31×10−88	1.31×10−88	0.311/0.026
2	*2702	44.93 (19.27–104.7)	1.23×10−18	4.48×10−11	0.049/9.003
3	*2704	15.94 (3.461–73.37)	0.00038	0.041	0.005/0.001
4	*5101	1.61 (1.241–2.089)	0.00034	0.15	0.1/0.117
5	*4006	3.06 (1.419–6.613)	0.0044	0.75	0.009/0.008
Iranian Cohort
Round	HLA-B allele	OR (95% CI)	Stepwise conditional P-value	Initial P-value
1	*2705	57.54 (38.17–86.74)	9.65×10−72	9.65×10−72	0.345/0.017
2	*2707	15.84 (4.97–50.53)	3.05×10−6	0.0059	0.013/0.003
3	*2704	21.3 (4.37–103.8)	0.00015	0.021	0.008/0.001
4	*4002	8.03 (2.14–30.06)	0.0020	0.036	0.008/0.001
5	*1503	8.22 (2.17–31.09)	0.0019	0.23	0.006/0.003
6	*2702	27.28 (3.02–246.5)	0.0032	0.077	0.005/0.0007

CI, confidence interval; OR, odds ratio.

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).

Table 2

Genome-wide significant and suggestive non-MHC loci for ankylosing spondylitis in the Turkish cohort.

Chr	Posa	SNP	P-value	OR (95%CI)	Nearby genes	Risk /Prot	MAFCase /Con
4	11,504,491	rs10000518	6.05×10−6	0.48(0.35–0.66)	HS3ST1	G/A	0.03/0.07
4	23,346,184	rs34939008	7.20×10−6	0.74(0.65–0.84)	GBA3	G/A	0.39/0.46
7	38,706,817	rs10280089	1.08×10−6	1.44(1.24–1.67)	FAM183B	A/C	0.32/0.24
8	2,139,471	rs4876208	5.56×10−6	1.37(1.20–1.57)	MYOM2	C/G	0.39/0.32
8	24,116,722	rs1013210	7.88×10−8	0.65(0.56–0.76)	ADAM28	G/A	0.20/0.28
15	50,785,016	rs148783236b	1.01×10−11	0.52(0.43–0.63)	USP8	A/G	0.12/0.20
15	70,338,485	rs11639037	9.52×10−6	1.40(1.21–1.62)	TLE3	C/T	0.31/0.24
15	101,767,822	rs8036083	7.24×10−6	1.35(1.19–1.55)	CHSY1	G/T	0.47/0.40
16	3,293,407	rs61752717	7.63×10−12	5.34(3.31–8.62)	MEFV	C/T	0.06/0.01
16	3,293,407c	rs61752717	1.75×10−11	5.26 (3.24–8.54)	MEFV	C/T	0.058/0.01
18	14,723,700	chr18:14723700:D	8.78×10−6	0.74(0.65–0.84)	ANKRD30B	Y/Z	0.36/0.43
21	46,272,105	rs235316	7.85×10−6	1.44(1.23–1.69)	PTTG1IP	T/A	0.25/0.19

Genome-wide significant loci details appear in italics. All others had suggestive significance.

Chr, chromosome; CI, confidence interval; Con, control; MAF, minor allele frequency; OR, odds ratio; Pos, position; Prot, protective.

aHuman genome build hg19

brs148783236 is exm1161045 in the raw data

cexcluded four FMF cases

Table 3

Genome-wide significant and suggestive non-MHC loci for ankylosing spondylitis in the Iranian cohort.

Chr	Posa	SNP	P-value	OR (95% CI)	Nearby genes	Risk /Prot	MAF Case/Con
2	62568221	rs13001372	5.33×10−7	1.56 (1.31–1.85)	2p15	A/G	0.49/0.40
5	96126308	rs27529	1.12×10−9	1.74 (1.46–2.08)	ERAP1	A/G	0.48/0.39
9	78503158	rs7874251	6.35×10−6	0.66 (0.55–0.79)	PCSK5	A/G	0.33/0.42
10	90754487	rs2031610	4.67×10−6	1.49 (1.26–1.77)	FAS	C/T	0.50/0.40
12	26043525	rs7298011	4.11×10−6	1.62 (1.32–1.99)	RASSF8	A/G	0.25/0.17
19	17451098	rs4808624	9.99×10−6	1.48 (1.24–1.75)	DDA1	G/C	0.47/0.44

Genome-wide significant loci details appear in italics. All others had suggestive significance.

Chr, chromosome; CI, confidence interval; Con, control; MAF, minor allele frequency; OR, odds ratio; Pos, position; Prot, protective.

aHuman genome build hg19

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 genome-wide 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 (r2>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 AS-associated 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.

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.

Chr	Posa	SNP	P-value	OR (95%CI)	Nearby genes	Risk /Prot
2	62,568,221	rs13001372	7.23×10−9	1.36 (1.19–1.56)	2p15	A/G
4	11,504,491	rs10000518	6.05×10−6	0.48 (0.35–0.66)	HS3ST1	G/A
4	23,346,184	rs34939008	7.19×10−6	0.74 (0.65–0.84)	GBA3	G/A
5	96,126,308	rs27529	1.36×10−9	1.39 (1.22–1.59)	ERAP1	A/G
7	38,706,817	rs10280089	1.1×10−6	1.44 (1.24–1.67)	AMPH	A/C
8	2,139,471	rs4876208	5.62×10−6	1.37 (1.20–1.57)	MYOM2	C/G
8	24,116,722	rs1013210	7.85×10−8	0.65 (0.56–0.76)	ADAM28	G/A
9	78,503158	rs7874251	6.34×10−6	0.66 (0.55–0.79)	PCSK5	G/A
10	90,754,487	rs2031610	4.6×10−6	1.49 (1.26–1.77)	FAS	C/T
12	26,043,525	rs7298011	4.07×10−6	1.62 (1.32–1.99)	RASSF8	A/G
15	50,785,016	rs148783236	1.44×10−13	0.58 (0.48–0.69)	USP8	A/G
15	70,338,485	rs11639037	9.71×10−6	1.40 (1.20–1.62)	TLE3	C/T
15	101,767,822	rs8036083	7.37×10−6	1.35 (1.19–1.55)	CHSY1	G/T
16	3,293,407	rs61752717	1.72×10−12	4.76 (2.94–7.68)	MEFV	C/T
17	147,35,492	rs6502398	5.18×10−7	0.76 (0.67–0.87)	17p12	C/A
17	260,98,781	rs7217335	1.34×10−6	1.34 (1.15–1.56)	NOS2	A/G
18	1,4723,700	chr18:14723700:D	8.79×10−6	0.74 (0.65–0.84)	ANKRD30B	Y/Z
19	17451098	rs4808624	9.99×10−6	1.48 (1.24–1.75)	GTPBP3	G/C
21	46,272,105	rs235316	7.87×10−6	1.44 (1.23–1.69)	PTTG1IP	T/A

Genome-wide significant loci details appear in italics. All others had suggestive significance.

Chr, chromosome; CI, confidence interval; Con, control; MAF, minor allele frequency; OR, odds ratio; Pos, position; Prot, protective.

aHuman genome build hg19

Suggestive associations (5×10−8 ADAM28 (encoding ADAM metallopeptidase domain 28; rs1013210, OR = 0.65, 95% CI = 0.56–0.76, P = 7.88×10−8). No association was seen at this locus in the Iranian cohort despite the MAF being similar in Turkish and Iranian controls (MAF = 0.28 and 0.25, respectively).

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 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).

Table 5

Comparison of MEFV M694V in Turkish, Iranian and the meta-analysis of both cohorts.

rs61752717	Case maf	Control maf	#cases	#controls	OR	P-value
Turkish	0.061	0.011	921	907	5.34	7.63×10−12
Iranian	0.012	0.0033	422	754	2.85	0.042
meta-analysis	0.047	0.0078	1344	1369	4.76	1.72×10−12

OR, odds ratio. maf, minor allele frequency.

Interaction analysis between MEFV variants and HLA-B27

Previous studies have suggested interaction between HLA-B27 and MEFV in AS (OR = 8.69 in HLA-B27-positive, and 21.8 in HLA-B27-negative cases) [11]. In the current study, rs61752717 shows stronger association in HLA-B27-negative patients (P = 8.93×10−15) than in HLA-B27-positive patients against all controls (P = 7.69×10−8, Table 6). In the Turkish cohort, in comparison of cases alone, rs61752717 risk allele carriage was higher in HLA-B27-negative than positive cases (OR = 2.67, 95% CI = 1.78–4.01, P = 1.28×10−6 in whole cases cohort, OR = 2.39, CI = 1.58–3.64, P = 9.54×10−6 excluding four cases with co-existent FMF and AS, Fisher Exact Test, Table 6), whereas no difference was observed in controls (P = 0.25). Compared with all-Turkish controls, the association of rs61752717 was much stronger in HLA-B27-negative cases (OR = 7.8, P = 8.93×10−15 in whole cohort, OR = 7.6, P = 3.92×10−14 excluding four cases with co-existent FMF and AS) than in HLA-B27-positive cases (OR = 4.3, P = 7.69×10−8), although no interaction observed in a specific logistic regression analysis (P = 0.97). Similarly, rs61752717 risk allele carriage was higher in HLA-B27-negative than positive cases in the Turkish-Iranian combined cohort (OR = 2.74, 95% CI = 1.86–4.03, P = 1.29×10−7 in entire cohort; OR = 2.47, 95% CI = 1.66–3.67, P = 5.78×10−6 excluding four cases with co-existent FMF and AS; Table 7). However, no significant difference was observed in rs61752717 risk allele carriage between HLA-B27 positive and HLA-B27-negative patients (OR = 1.44, 95% CI = 0.23–6.40, P = 0.70; Table 8) in the Iranian cohort, noting the extremely low power of this analysis given the rarity of rs61752717 in this cohort.

Table 6

Genotype counts and (%) of rs61752717 in HLA-B27-positive and HLA-B27-negative Turkish cohorts (including four HLA-B27-negative subjects with co-existent AS and FMF, three homozygous for rs61752717 ‘C’ allele, one homozygote for rs61752717 ‘T’ allele).

rs61752717	Case		Control		Total
	B27+	B27−	B27+	B27−
CC (+/+)	1 (0.2)	9 (3.4)	0 (0)	0 (0)	10
CT (+/−)	53 (8.2)	38 (14.3)	0 (0)	19 (2.3)	110
TT (−/−)	596 (91.7)	218 (82.3)	57 (100)	824 (97.7)	1695
CC or CT (+/+ or +/−)	54 (8.3)	47 (17.7)	0 (0)	19 (2.3)	120
Sum	650	265	57	843	1815

rs61752717 genotype counts differ between HLA-B27-positive and–negative cases (P = 1.22×10−6). Comparing allele frequencies between HLA-B27-positive and–negative controls, no difference was observed (CC genotype frequency comparison not possible due to zero counts of CC genotypes; P = 0.25). Risk allele carriage in HLA-B27-negative cases was higher than in positive cases (OR = 2.67, 95% CI = 1.78–4.01, P = 1.28×10−6 in whole cases cohort, OR = 2.39, 95% CI = 1.58–3.64, P = 9.54×10−6 excluding four cases with co-existent FMF and AS).

Table 7

Genotype counts and (%) of rs61752717 in HLA-B27-positive and HLA-B27-negative combined cohorts (including four HLA-B27-negative subjects with co-existent AS and FMF, three homozygous for rs61752717 ‘C’ allele, one homozygote for rs61752717 ‘T’ allele).

rs61752717	Case		Control		Total
	B27+	B27−	B27+	B27−
CC (+/+)	1 (0.1)	11 (3.0)	0 (0)	0 (0)	12
CT (+/−)	57 (6.1)	39 (10.6)	0 (0)	24 (1.6)	120
TT (−/−)	869 (97.3)	319 (86.4)	90 (100)	1479 (98.4)	2757
CC or CT (+/+ or +/−)	58 (6.3)	50 (13.6)	0 (0)	24 (1.6)	132
Sum	927	369	90	1503	2889

Risk allele carriage in HLA-B27-negative cases was higher than positive cases (OR = 2.74, 95% CI = 1.86–4.03, P = 1.29×10−7 in whole cases cohort, OR = 2.47, 95% CI = 1.66–3.67 P = 5.78×10−6 excluding four cases with co-existent FMF and AS).

Table 8

Genotype counts and (%) of rs61752717 in HLA-B27-positive and HLA-B27-negative Iranian cohorts.

rs61752717	Case		Control		Total
	B27+	B27−	B27+	B27−
CC (+/+)	0 (0)	2 (1.8)	0 (0)	0 (0)	2
CT (+/−)	5 (1.5)	1 (0.9)	0 (0)	5 (0.7)	11
TT (−/−)	299 (98.4)	111 (97.4)	35 (100)	710 (99.3)	1155
CC or CT (+/+ or +/−)	5 (1.5)	3 (2.6)	0 (0)	5 (0.7)	13
Sum	304	114	35	715	1168

No significant difference was observed in risk allele carriage between HLA-B27-negative cases and positive cases (OR = 0.37, 95% CI = 0.084–1.623, P = 0.14)

Interaction analysis between MEFV variants and HLA-B51

The MEFV nsSNP rs61752717 is also associated with Behçet’s Disease (BD) [15-18]. Joint involvement is common in BD, and about 5% of the BD patients developed sacroiliitis [19]. 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-B51-positive 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 co-existent AS and FMF.

Table 9

Serum cytokine levels by case or control status and MEFV M694V carrier status.

Cytokine	Control (mean (SE))	Case(mean (SE))	Case-controlP-value	MEFV 694M Carrier	MEFV 694V Carrier	GenotypeP-value
IL-1	0.19 (0.042)	1.98 (0.39)	0.00030	0.95 (0.21)	3.39 (0.83)	0.0017
IL-1*		2.16 (0.70)	0.00041	0.93 (0.21)	3.43 (0.86)	0.0093
IL-17	0.41 (0.086)	2.21 (0.68)	0.033	1.96 (0.90)	2.56 (1.06)	0.89
IL-17*		1.95 (0.41)	0.016	1.99 (0.92)	2.39 (1.11)	0.78
IL-23	8.20 (3.85)	26.90 (4.47)	0.0036	16.72 (3.86)	40.81 (8.63)	0.0068
IL-23*		28.09 (4.61)	0.0023	17.14 (3.93)	43.72 (9.00)	0.0038

Cytokine levels are expressed as pg/mL.

*excluding the subjects with co-existent AS and FMF.

Note that the carrier MEFV M694V status is for cases only.

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-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 gene-gene 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-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-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-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.

Materials and methods

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 (r2 < 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.cog-genomics.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 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.

Web resources

SHELLFISH; http://www.stats.ox.ac.uk/~davison/software/shellfish/shellfish.php

PLINK v1.90b3.36; https://www.cog-genomics.org/plink2

AS-associated variations in MEFV in Turkish and Iranian case-control cohorts separately.

(DOCX)

Click here for additional data file.

Genotype counts and (%) of rs61752717 in HLA-B51-positive and HLA-B51-negative cohorts in Turkish dataset (including four subjects with co-existent AS and FMF, three B51-positive and homozygote for rs61752717 ‘C’ allele, one HLA-B51-negative homozygote for rs61752717 ‘T’ allele).

(DOCX)

Click here for additional data file.

Genotype counts and (%) of rs61752717 in HLA-B51-positive and HLA-B51-negative cohorts in Iranian dataset.

(DOCX)

Click here for additional data file.

Genotype counts and (%) of rs61752717 in HLA-B51-positive and HLA-B51-negative cohorts in combined dataset (including four subjects with co-existent AS and FMF, three B51-positive and homozygote for rs61752717 ‘C’ allele, one HLA-B51-negative homozygote for rs61752717 ‘T’ allele).

(DOCX)

Click here for additional data file.

Manhattan plot for chr1-22 (genotyped SNPs) in the Iranian cohort.

(TIF)

Click here for additional data file.

Manhattan plot for chr1-22 (genotyped SNPs) in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 2p15 and linkage disequilibrium with lead SNP rs13001372 in the Iranian cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 5q15 and linkage disequilibrium with lead SNP rs27529 in the Iranian cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 9q21.13 and linkage disequilibrium with lead SNP rs7874251 in the Iranian cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 10q23.31 and linkage disequilibrium with lead SNP rs2031610 in the Iranian cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 12p12.1 and linkage disequilibrium with lead SNP rs1488826 in the Iranian cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 19p13.11 and linkage disequilibrium with lead SNP rs4808624 in the Iranian cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 4p15.33 and linkage disequilibrium with lead SNP rs10000518 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 4p15.2 and linkage disequilibrium with lead SNP rs34939008 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 7p14 and linkage disequilibrium with lead SNP rs10280089 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 8p21.2 and linkage disequilibrium with lead SNP rs1013210 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 8p23.3 and linkage disequilibrium with lead SNP rs4876208 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 11q21 and linkage disequilibrium with lead SNP rs1271188 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 15q21.2 and linkage disequilibrium with lead SNP EXM1161045 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 15q26.3 and linkage disequilibrium with lead SNP rs8036083 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 15q23 and linkage disequilibrium with lead SNP rs11639037 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 16p13.3 and linkage disequilibrium with lead SNP rs61752717 in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 18p11.21 and linkage disequilibrium with lead variant chr18:14723700:D in the Turkish cohort.

(TIF)

Click here for additional data file.

Locus association plot of chromosome 21q22.3 and linkage disequilibrium with lead SNP rs235316 in the Turkish cohort.

(TIF)

Click here for additional data file.

Plot of λ and λ1000 in Iranian and Turkish cohorts for non-MHC region.

λ and λ1000 are coloured by red and blue, respectively. The x-axis and y-axis are the 1st and corresponding λ and λ1000 values.

(TIF)

Click here for additional data file.

Plot of λ and λ1000 in the combined cohort for non-MHC region.

λ and λ1000 are coloured by red and blue, respectively. The x-axis and y-axis are the 1st and corresponding λ and λ1000 values.

(TIF)

Click here for additional data file.

Plot of principal component analysis (PCA) in the Iranian cohort.

Healthy controls (HC) are in red while ankylosing spondylitis patients (AS) are in blue. The x-axis and y-axis are the 1st and 2nd principal component from the final PCA.

(TIF)

Click here for additional data file.

Plot of principal component analysis (PCA) in the Turkish cohort.

Healthy controls (HC) are in red while ankylosing spondylitis patients (AS) are in blue. The x-axis and y-axis are the 1st and 2nd principal component from the final PCA.

(TIF)

Click here for additional data file.

Plot of principal component analysis (PCA) in the combined cohort of Turkish and Iranian samples.

Iranian samples are in red while Turkish samples are in blue. The ankylosing spondylitis (AS) and healthy control (HC) are indicated by solid dots and triangles, respectively. The x-axis and y-axis are the 1st and 2nd principal component from the final PCA.

(TIF)

Click here for additional data file.

Q-Q plot of the Iranian cohort. The λ1000 was 1.041 for the entire dataset and 1.022 excluding the MHC region.

Red dashed lines indicate 95% CI of expected χ2 values, while magenta and blue points are the observed χ2 values of the entire dataset and the dataset excluding the MHC region, respectively.

(TIF)

Click here for additional data file.

Q-Q plot of Turkish cohort, with first five PCA components as covariates.

The λ1000 was 1.051 for overall dataset and 1.024 excluding MHC region. Red dashed lines indicate 95% CI of expected χ2 values, while magenta and blue points are the overserved χ2 values of the entire dataset and dataset excluding the MHC region, respectively.

(TIF)

Click here for additional data file.

Q-Q plot of the combined cohort (with first four PCA components as covariates).

The λ1000 was 1.086 for the entire dataset and 1.021 excluding the MHC region. Red dashed lines indicate 95% CI of expected χ2 values, while magenta and black points are the overserved χ2 values of the entire dataset and dataset excluding the MHC region, respectively.

(TIF)

Click here for additional data file.