Individuals with Behçet's disease suffer from episodic inflammation often affecting the orogenital mucosa, skin and eyes. To discover new susceptibility loci for Behçet's disease, we performed a genome-wide association study (GWAS) of 779,465 SNPs with imputed genotypes in 1,209 Turkish individuals with Behçet's disease and 1,278 controls. We identified new associations at CCR1, STAT4 and KLRC4. Additionally, two SNPs in ERAP1, encoding ERAP1 p.Asp575Asn and p.Arg725Gln alterations, recessively conferred disease risk. These findings were replicated in 1,468 independent Turkish and/or 1,352 Japanese samples (combined meta-analysis P < 2 × 10(-9)). We also found evidence for interaction between HLA-B*51 and ERAP1 (P = 9 × 10(-4)). The CCR1 and STAT4 variants were associated with gene expression differences. Three risk loci shared with ankylosing spondylitis and psoriasis (the MHC class I region, ERAP1 and IL23R and the MHC class I-ERAP1 interaction), as well as two loci shared with inflammatory bowel disease (IL23R and IL10) implicate shared pathogenic pathways in the spondyloarthritides and Behçet's disease.

*These authors made equal contributions to the study.†These authors jointly directed this work.

Behçet's disease (BD) is a form of vasculitis that manifests with orogenital ulcers, uveitis, skin inflammation, arthritis, enterocolitis, and inflammation in other organs1,2. BD is relatively common in Turkey, Japan, and modern-day countries that fall on Marco Polo's ancient Silk Routes, and is an important cause of vision loss in these countries2. Genetic risk factors contribute to disease-susceptibility. HLA-B*51 is the most strongly associated risk factor for BD, confirmed in multiple populations3-5. Although its association was established more than three decades ago, the role of HLA-B*51 in disease pathogenesis remains elusive5. In addition to HLA-B*51, two recent independent GWASs identified variants in regions encompassing MHC-I, IL10, and IL23R associated with BD in both the Turkish and Japanese populations6,7. However, the combined effects of these genetic factors do not fully explain the observed disease heritability.

The pathobiology of BD is also largely unknown. In 1974, based on clinical features, Moll et al. proposed the concept of “seronegative spondylarthritides”, and included BD along with ankylosing spondylitis (AS), psoriatic arthritis, reactive arthritis, and inflammatory bowel disease (IBD)8. Since then, the inclusion of BD within the spondyloarthropathy (SpA) category has been debated as BD patients rarely exhibit sacroiliitis, and BD is associated with HLA-B*51 rather than HLA-B*279-11. On the other hand, overlapping extra-articular clinical manifestations (inflammation in the eyes, skin, and intestine), genetic associations at MHC-I and IL23R6,7,12-14, and the effectiveness of tumor necrosis factor (TNF)-α blockade15,16, suggest shared pathogenesis between BD and SpA. Furthermore, IL10 and IL23R variants were found associated with both BD and IBD (Crohn's disease and ulcerative colitis), implicating common inflammatory pathways between BD and these other members of the SpA group.17,18

We selected 21 SNPs from the novel loci identified by imputation and one SNP from the previously reported Japanese GWAS7 for genotyping in a Turkish replication collection, comprising newly collected 838 Turkish cases and 630 controls (Supplementary Table 2). Four promising loci, STAT4 (signal transducer and activator of transcription 4), KLRC4 (killer cell lectin-like receptor subfamily C, member 4), the CCR1-CCR3 locus, and IL12A (interleukin-12 alpha chain) were then selected for validation of the imputed data (by direct genotyping) and fine-mapping studies in the original Turkish GWAS samples. As shown in Figure 2a-c, we found the strongest signals at rs7616215 3′ of CCR1, at rs7574070 in intron 3 of STAT4, and at nonsynonymous SNPs rs2617170 and rs1841958, encoding KLRC4 p.Asn104Ser and p.Ile129Ser. In a meta-analysis, we combined the Turkish GWAS and replication data, and if polymorphic, the Japanese replication data (from the reported GWAS collection with 612 cases and 740 controls7) using Cochran-Mantel-Haenszel tests (Table 1 and Figure 2). Three loci (CCR1-CCR3, STAT4, and KLRK1-KLRC1) were associated with BD at genome-wide significance (p=1.34 × 10−9 to 4.30 × 10−13). Additionally, the IL12A locus exhibited suggestive association (p=6 × 10−7). Of the 612 Japanese cases ascertained using the Japanese diagnostic criteria7, 496 also fulfilled the International Study Group Criteria19. An analysis including only cases that met the International Study Group criteria revealed genome-wide significance for the same three loci despite the reduced numbers (Supplementary Table 3).

In an attempt to reduce genetic heterogeneity, we performed genome-wide association tests in the subset of GWAS discovery patients with uveitis (435 cases, 1,278 controls)6. Neither the basic allelic test nor the dominant model test showed associations outside of the MHC. However, when we applied a recessive model, one SNP in the ERAP1-ERAP2 locus exhibited association with a p-value close to genome-wide significance level (rs2927615, p=1.02 × 10−7). We performed fine-mapping of this region in the uveitis subset of the GWAS discovery collection and identified rs10050860 and rs17482078, encoding ERAP1 p.Asp575Asn and p.Arg725Gln, which conferred risk for BD with uveitis in a recessive model (Figure 2d). A meta-analysis of p.Arg725Gln combining the Turkish discovery collection and the Turkish replication collection (with 370 BD cases with uveitis and 630 controls) exceeded the three model threshold for genome-wide significance and revealed a large effect size of the homozygous p.Arg725Gln genotype on BD with uveitis (odds ratio=4.56, p=4.73 × 10−11, Table 1 and Figure 2d).

Because a recessive model was required to detect the ERAP1 association in BD patients with uveitis, we tested whether the recessive model would reveal the ERAP1 p.Arg725Gln association with BD susceptibility in the combined uveitis and non-uveitis samples. A meta-analysis of the GWAS and replication collections found significant association of the homozygous p.Arg725Gln genotype with BD susceptibility (p=4.35 × 10−8, Supplementary Table 4). The minor allele frequency of rs2927615 (a variant in strong linkage disequilibrium [LD] with p.Arg725Gln) was too low in the Japanese population (1.8% in BD cases, 2.0% in controls) to evaluate recessive effects. Furthermore, none of the Japanese GWAS SNPs from the regions encompassing ERAP1 or IL12A (rs17810546 was not polymorphic in the Japanese population) were associated with BD (Supplementary Table 5).

ERAP1 is an endoplasmic reticulum expressed amino peptidase that functions to trim peptides for loading onto MHC Class I20. Previous GWASs have established associations of ERAP1 variants in psoriasis21, 22 and AS12. ERAP1 p.Asp575Asn and p.Arg725Gln, which are in strong LD, confer protection against these diseases through reduced peptide trimming and antigen presentation by MHC-Class I12,22-24. Of note, recent reports have shown that these ERAP1 variants confer protection preferentially in HLA-B*27 positive individuals in AS23 and HLA-C*06 positive individuals in psoriasis22, suggesting that peptide processing and binding/presentation mechanisms contribute to the pathogenesis of these diseases.

We therefore tested for an interaction between HLA-B*51 and ERAP1 in BD. The ERAP1 variants preferentially conferred risk for BD in HLA-B*51 positive individuals (p-value for interaction=0.0009 from a logistic likelihood ratio test comparing the full model including a multiplicative interaction term with the reduced model without interaction term) in the combined Turkish GWAS and replication (including uveitis and non-uveitis) samples (Figure 3). Furthermore, ERAP1 p.Arg725Gln homozygosity was associated with an odds ratio for BD of 3.78 [95% CI 1.94-7.35] in the HLA-B*51 positive individuals versus an odds ratio of 1.48 [95% CI 0.78-2.80] in the HLA-B*51 negative individuals. This finding indicates that the disease-associated peptidase variant contributes to disease susceptibility through an interaction with the HLA-B*51 protein.

Homozygosity for the ERAP1 variants is associated with increased risk for BD, but decreased risk for AS23 and psoriasis22. The difference between risk and protection among these three diseases may depend on the variability and binding affinities of peptides loaded onto the respective MHC Class I molecules, which can affect their stability and function. Indeed, repertoires of MHC-bound peptides are altered in ERAP1 deficient mice25. ERAP1 p.Arg725Gln-related alterations might affect the repertoire of peptides that bind to HLA-B*51, which is known for its promiscuous peptide binding features5,26. The recessive nature of the ERAP1 effect in BD (one wild type copy in heterozygotes is sufficient to obscure the risk effect of the mutant allele) suggests that homozygotes fail to produce one or more disease-protective peptides.

The KLRC4 BD-associated SNP (Table 1) is within a haplotype block containing five natural killer (NK) cell receptor genes (KLRK1, KLRC1-4) (Figure 2c). Two non-synonymous variants in KLRC4, rs1841958 and rs2617170, encoding KLRC4 p.Ile29Ser and Asn104Ser, are found on the BD-protective haplotype of this LD block. This haplotype has been associated with reduced peripheral blood leukocyte cytotoxicity and increased incidence of cancer27. Conditional logistic regression analysis, conditioning on the KLRC4 Asn104Ser variant, showed no additional independent association signals within this LD block (Supplementary Figure 1a). KLRC4, also called NKG2F, encodes a c-type lectin receptor whose function is largely unknown. A possible clue to its function may be found in a related family member, NKG2D, encoded by KLRK1 and also located within the disease-associated haplotype block. NKG2D is expressed on NK cells and γδ T cells, and can act as a co-stimulatory molecule for CD4+ and CD8+ T cells28,29. Interestingly, a ligand of the NKG2D receptor is MICA (the MHC Class I chain-related protein A)28. The MICA gene is located within the MHC region and SNPs within the MICA locus are in linkage disequilibrium with HLA-B*516. The importance of NK receptors in BD pathogenesis is also supported by the observation that the strongest linkage peak in Turkish familial BD (LOD score of 3.94) is found at chr12p12-13, which includes the KLRK region locus30.

The disease-associated variants in the 3′ flanking region of CCR1-CCR3 (rs7616215) and within the third intron of STAT4 (rs7574070) are noncoding and are not in strong LD with any coding variants (Figure 2a-b). The CCR1-CCR3 locus contains a cluster of chemokine receptor genes within the LD block. Logistic regression analysis conditioning on covariate rs7616215 revealed only a single association signal in the region (Supplementary Figure 1b). ENCODE data indicated that rs7616215 and rs7574070 are located within DNase I hypersensitivity and histone 3 lysine 4 methylation sites, suggesting effects on transcription.

STAT4 mRNA expression was higher in individuals with the risk allele A (Figure 4e and 4f). The BD-associated variant rs7574070 (and its surrogate rs7572482) is in poor LD with the previously reported autoimmune disease-associated STAT4 variant, rs757486532; in fact it is located two LD blocks away, suggesting the associations are independent, although both variants are located within the large third intron of STAT4 (Figure 2b). Both variants are associated with increased expression of STAT4 (ref 33 and data presented here), but the genetically distinct disease-associations suggest different STAT4 regulatory mechanisms in BD compared with rheumatoid arthritis, systemic lupus erythematosus, and other autoimmune diseases. Smaller effect sizes observed for the associations of the CCR1 and STAT4 variants in the Japanese replication (Table 1) could be explained by “the winner's curse”34.

In conclusion, this study adds substantially to the understanding of genetic factors that contribute to BD susceptibility (the new loci are CCR1-CCR3, STAT4, KLRK1-KLRC4, and ERAP1). Furthermore, the results support an emerging concept delineating common pathogenic mechanisms for BD and the SpA. BD, AS, and psoriasis are inflammatory disorders affecting the skin, eyes, and joints, with significant MHC Class I associations (B*51 for BD, B*27 for AS, and C*06 for psoriasis). Recent genetic studies implicate variants of IL23R, encoding an upstream molecule in Th17 activation, in susceptibility to all three disorders. The present work adds ERAP1 to the list of shared genetic factors and furthermore, interactions between MHC Class I and ERAP1 are also found in all three of these diseases. ERAP1 trims peptides for proper loading onto Class I antigens, thus suggesting that peptide-MHC Class I interactions contribute to all three of these diseases. These data suggest the existence of shared inflammatory pathways among these diseases leading to the possibility of common therapeutic strategies, while raising questions about the specific disease characteristics, which may be related to their different MHC Class I associations.

Online methods

Patients

1,215 Turkish Behçet's disease (BD) cases and 1,278 genetically matched controls used in previous GWAS were studied6. Individuals who also met the Tel-Hashomer clinical criteria for the diagnosis of familial Mediterranean fever (FMF)35 were excluded (n=6). For replication, an additional similarly collected 838 Turkish cases (none fulfilled FMF criteria) and 630 controls, and 612 Japanese BD cases and 740 control samples enrolled in the previous GWAS7 were included. Turkish BD patients fulfilled the International Study Group diagnostic criteria for Behçet's disease19. All Japanese BD patients fulfilled the Japanese BD diagnostic criteria7 and 496 of them also fulfilled the International Study Group criteria. All study participants provided written informed consent, and the study was approved by the Ethics Committees of each investigative institution.

Genotype imputation

We imputed genotypes of 1,209 BD cases and 1,278 controls using MACH v1.0. For a reference panel, we used 96 Turkish healthy controls who participated in the previous BD GWAS using the HumanHap370CNV chip (Illumina) and additionally genotyped on the Human OMNI 1M chip (Illumina, San Diego, CA). For quality control, we excluded SNPs from the reference panel if they had a minor allele frequency less than 5%, deviated from Hardy–Weinberg equilibrium (P<0.0001), or had a call rate below 95%, yielding 814,474 SNPs for the imputation. Quality scores for the imputation are shown in Supplementary Table 6. SNPs with Rsq<0.3 were excluded from the association analysis. A total of 779,465 imputed SNPs were included in the genome-wide association analysis.

Validation and fine-mapping

The Turkish GWAS imputation provided the discovery data. Twenty-two SNPs from the novel genetic loci were selected for evaluation in the Turkish replication samples. Three of these SNPs with P < 0.001 and one SNP from the phenotypic subset analysis identified four regions for validation and fine-mapping by i-PLEX assays (TOF-MS, Sequenom, San Diego, CA) using the original Turkish GWAS samples6. The IL12A locus with two SNPs with P < 0.05 in the replication samples was also investigated, but failed to reach genome-wide significance. For variants that failed TOF-MS design or reaction, TaqMan genotyping was performed (Applied Biosystems, Foster City, CA). Genotyping was performed in an unbiased fashion by masking the phenotype of the samples. For the fine-mapping, we used the Tagger SNP selection tool from HapMap to select SNPs to augment the coverage of the GWAS SNPs with the intent to obtain 100% coverage of the HapMap Phase III SNPs with greater than 5% minor allele frequency in the CEU HapMap population with pairwise r2 > 0.8. Although already tagged, additional SNPs with r2 > 0.8 with the most significantly associated SNP of the region were also included. Genotyping of the same samples used for the Turkish GWAS discovery collection was performed6. After quality control, the resulting coverage for the CCR1 locus (chr3:45441901-46908964, hg18) was 92%, the STAT4 locus (chr2:191602386-191769025) was 84%, the KLRK1-KLRC1 locus (chr12:10329925-10557292) was 92%, and the ERAP1-2 locus (chr5:96026703-96305246) was 94%. The most significantly associated marker from each region was used for the replications. The Japanese replication collection genotypes were from the Japanese GWAS7. The Turkish replication collection samples were genotyped by TOF-MS or TaqMan assay.

mRNA expression data and migration assay

CCR1 mRNA expression data were extracted from the report by Zeller and colleagues, which includes genome-wide SNP data along with mRNA expression array data from monocytes of n = 1490 European ancestry individuals31. SNP rs7616215 showed association with CCR1 mRNA expression (p=9.54×10−6), whereas CCR3 data were not reported, indicating that the association of rs7616215 with CCR3 is less significant (p>5×10−5). STAT4 mRNA expression data in lymphoblastoid samples were obtained from Gene Expression Omnibus (GEO) datasets36, 37. Primary human monocytes from unrelated healthy volunteers were isolated from peripheral blood mononuclear cells by MACS Human Monocyte Isolation Kit (Miltenyi Biotec, Gladbach, Germany). RNA was isolated from PBMCs by RNeasy kit (Qiagen, Valencia, CA) and preserved at −80 C° until used. cDNAs were prepared from DNase I (Invitrogen, Carlsbad, CA)-pretreated-RNA using SuperScript II according to manufacturer's protocol (Invitrogen). Q-PCR gene expression assay for human CCR1 (Hs00174298_m1) was purchased from Applied Biosystems (Foster City, CA). Human GAPDH (4310884E) served as an internal control. Multiplex PCR (GAPDH with CCR1) was performed in triplicate or quadruplicate as in the manufacturer's protocol (Applied Biosystems). The ΔΔCT method was used for the analysis (n=93).

A monocyte migration assay was performed with 24-well Transwell 5 μm polycarbonate membrane chambers (Costar, Corning, NY). 1% bovine serum albumin RPMI1640 medium was used for incubation of the cells. 5×104 cells were seeded in the upper chamber and the CCR1 chemokine MIP1-α (0 or 10ng/ml) (R&D, Minneapolis, MN) was placed in lower chamber. After incubation for 2 hours, cells that migrated into the membrane were fixed and stained by DiffQuik (Siemens, Newark, DE). Cells were counted in 5 high power fields on membranes from each of two duplicate wells. The relative change in cell migration in response to the chemokine was determined by dividing the migrated cell count obtained with MIP1-α by the cell count obtained from the same sample without MIP1-α (n=61).

Statistical analysis

Genome-wide SNP association tests were performed based on allelic tests by comparing the allele frequencies between BD cases and controls, using Golden Helix SVS 7.5.2 software (Golden Helix, Bozeman, MT). P<5×10−8 was considered genome-wide significance. For the uveitis subset analysis, as well as the basic allelic test, we applied dominant and recessive genetic model tests and therefore employed 1.67×10−8 for genome-wide significance. Conditional analyses were performed by fitting the logistic regression model with SNPs rs2617170 (for KLRK1-KLRC1) or rs7616215 (for CCR1-CCR3) as covariates. For meta-analysis, Cochran-Mantel-Haenszel tests were performed. To test for the interactive effects, we fit the log likelihood of the full model, including additive terms of the main effects and a multiplicative term of the interaction effect, versus a reduced model of the additive terms of the main effects only. P value was calculated by likelihood ratio tests between full and reduced models based on 1 degree of freedom Chi-square test statistics. For comparisons in expression and chemotaxis data, the Kruskal-Wallis rank sum test as a non-parametric version of one-way ANOVA and Spearman's rank correlation test were performed.

This research was supported by the Intramural Research Programs of the National Human Genome Research Institute and the National Institute of Arthritis and Musculoskeletal and Skin Diseases, and the Center for Human Immunology, Autoimmunity and Inflammation, of the National Institutes of Health, USA, and by the Istanbul University Research Fund, and Research on Specific Disease of the Health Science Research Grants from the Japanese Ministry of Health, Labor, and Welfare, and the Japan Rheumatism Foundation. We thank Dr. Alexander Wilson (Genometrics Section, Inherited Disease Research Branch, National Human Genome Research Institute, National Institutes of Health, Baltimore, Maryland, USA) for helpful comments on this manuscript.