New Sequence Variants in HLA Class II/III Region Associated with Susceptibility to Knee Osteoarthritis Identified by Genome-Wide Association Study

Osteoarthritis (OA) is a common disease that has a definite genetic component. Only a few OA susceptibility genes that have definite functional evidence and replication of association have been reported, however. Through a genome-wide association study and a replication using a total of ∼4,800 Japanese subjects, we identified two single nucleotide polymorphisms (SNPs) (rs7775228 and rs10947262) associated with susceptibility to knee OA. The two SNPs were in a region containing HLA class II/III genes and their association reached genome-wide significance (combined P = 2.43×10−8 for rs7775228 and 6.73×10−8 for rs10947262). Our results suggest that immunologic mechanism is implicated in the etiology of OA.


Introduction
We are living in the ''Bone and Joint Decade'' (http://www. boneandjointdecade.org/). As the WHO initiative shows, bone and joint diseases are serious problems all over the world, putting us under severe medical, economical and social burden. Osteoarthritis (OA; MIM 165720) is one of the most common diseases among them. OA affects synovial joints of all over the body, mainly knee, hip, hand and spine. OA is characterized by progressive loss of articular cartilage and, often, proliferation of synovium and bone, which lead to pain, loss of joint function and disability. More than tens of millions patients in the world are suffering from this non-lethal, but intractable disease, and the number is relentlessly increasing; however, its etiological picture remains unclear and we have no fundamental treatment for it.
OA is a polygenic disease. Both environmental and genetic factors contribute to its etiology and pathogenesis [1]. To understand its genetic factor, identification of its susceptibility gene(s) must be the first step. Many OA susceptibility genes identified by candidate-gene association studies have been reported, but only a few have supporting functional evidence and replication of the results in different populations [1,2]. Largescale association studies including the genome-wide association study (GWAS) using high-density single nucleotide polymorphisms (SNPs) have been reported by a few groups in Asia and Europe [3][4][5][6], but only a gene fulfilled genome-wide significance level [2]. The genetic basis of OA susceptibility remains largely uncharacterized. To identify OA susceptibility gene(s), we conducted a GWAS for knee OA and identified two SNPs with genome-wide significance level.

Samples
Characteristics of each cohort group are shown in Table 1. Case samples of GWAS for the Japanese population were obtained from several medical institutes in Japan, as previously described [5,7]. Knee OA was diagnosed on the basis of clinical and radiographic findings using previously described criteria [5,7]. Rheumatoid arthritis (RA) and polyarthritis associated with autoimmune diseases were excluded, as were secondary OA due to crystal deposition (gout and pseudogout), posttraumatic OA and infection-induced OA. Patients who had clinical and radiographic findings suggestive of skeletal dysplasias, including overt short stature, multiple symmetric involvements of epiphyses and a definitely positive Mendelian family history were also excluded from the study. The control groups consisted of 3,396 individuals that were registered in the Leading Project for Personalized Medicine in the Ministry of Education, Culture, Sports, Science and Technology, Japan as the subjects with diseases unrelated to OA and the volunteers in the Osaka-Midosuji Rotary Club, Osaka, Japan [8]. For replication study, we recruited populationbased cohorts from inhabitants of Odai and Minami-ise town (previously Miyagawa village and Nansei town, respectively in the Mie prefecture in Japan) [9]. The Spanish and Greek knee OA and control populations were recruited as described previously from the Hospital Clinico de Santiago, the Departments of Biology and Genetics and of Orthopaedics, University of Thessaly and the Institute of Musculoskeletal Sciences [10]. All the participants provided written informed consent. This research project was approved by the ethical committees at Center for Genomic Medicine (formerly, SNP Research Center), RIKEN and the participating institutions.

SNP genotyping
For the GWAS, we genotyped 906 patients with OA and 3,396 controls using Illumina HumanHap550v3 Genotyping BeadChip. After excluding seven cases with call rate of ,0.98, we applied SNP QC (call rate of $0.99 in both cases and controls and P value of Hardy-Weinberg equilibrium test of $1.0610 26 in controls). Finally, 459,393 SNPs on autosomal chromosomes passed the QC filters and were further analyzed. Among the SNPs analyzed in the GWAS, we selected top 15 SNPs showing the smallest P values (P,1610 25 ) for the replication study using an independent 514 Japanese subjects from a resident cohort. SNPs with minor allele frequency of #0.1 in both case and control samples were excluded from the further analysis. In the replication analysis, we genotyped SNPs using the multiplex PCR-based invader assay (Third Wave Technologies) or by direct sequencing of PCR products using ABI 3700 DNA analyzers (Applied Biosystems), or by SNaPshot Multiplex System (Applied Biosystems) according to manufacturers' protocols.

Statistical analysis
In the GWAS and replication analyses, we applied Fisher's exact test to two-by-two contingency table in three genetic models: an allele frequency model, a dominant-effect model, and a recessive-effect model. We conducted the meta-analysis using the Mantel-Haenszel method. We examined heterogeneity among studies by using the Breslow-Day test. Significance levels after the Bonferroni correction for multiple testing were P = 1.09610 27 (0.05/459,393). Age, gender-and BMI-adjusted odds ratios were obtained by logistic regression analysis [11]. Odds ratios and confidence intervals were calculated using the risk allele as a reference. We analyzed the haplotype association using Haploview software [12]. We conducted a principal component analysis to detect population stratification [13].

Software
For general statistical analysis, we used R statistical environment version 2.6.1 or Microsoft Excel. Drawing the LD map, estimation of haplotype frequencies and analysis of haplotype association were performed by Haploview software.

Results
To identify genetic variants that determine OA susceptibility, we conducted a GWAS in Japanese knee OA. We examined 906 individuals with knee OA and 3,396 control individuals (Table 1) using Illumina HumanHap550v3 Genotyping BeadChip. After confirming the data quality, we compared the results of 459,393 SNPs between cases and controls by Fisher's exact test for three genetic models: allelic, dominant or recessive (Figure 1). Fifteen SNPs selected for the replication study that had the smallest P values (minimum P,1610 25 ) were next genotyped in an independent set of 167 Japanese knee OA individuals and 347 Japanese controls from a resident cohort study. Through these studies, only two SNPs, rs7775228 (combined P = 2.43610 28 ; OR = 1.34; 95% CI = 1.21-1.49) and rs10947262 (combined  P = 6.73610 28 ; OR = 1.32; 95% CI = 1.19-1.46) were significant even after the Bonferroni correction for multiple testing (P = 1.09610 27 ) ( Table 2). The two SNPs showing significant associations are located within a 340-kb region within the HLA locus, including BTNL2, HLA-DRA, HLA-DRB5, HLA-DRB1, HLA-DQA1 and HLA-DQB1 (Figure 2). Although the HLA region is known to show extensive linkage disequilibrium (LD) spanning over 7 Mb, only SNPs in the 340-kb region showed strong associations with OA (Figure 2), and SNPs outside of this region did not have significant association.
Application of the Cochrane-Armitage test to all the tested SNPs indicated that the genetic inflation factor lambda was 1.08 for GWAS (Figure 3), implying a low possibility of false positive associations due to population stratification. We also carried out age, gender-and BMI-adjusted analysis using a logistic regression model, and confirmed similar association after adjustment (data not shown). The principal component analysis [13] revealed that there was no evidence for population stratification between the two control groups used for the GWAS ( Figure S1).
To check the association of rs7775228 and rs10947262 in different ethnic populations, we examined the association of the SNPs with knee OA in two European Caucasian populations from Greece and Spain. We genotyped a total of 813 OA and 1,071 control subjects (Table 1). We conducted the meta-analysis using the Mantel-Haenszel method. The combined European results for rs7775228 were not significant with OR (95%CI) of 0.93 (0.76-1.13) ( Table 2), while those for rs10947262 were supportive with OR (95%CI) of 1.29 (1.03-1.61). rs10947262 showed replication in the Greek population and the same trend in the Spanish population (Table 2). A meta-analysis of the Japanese and two European studies gave more significant association (combined P = 5.10610 29 ).
We estimated the pairwise LD indexes (D' and r 2 ) between rs7775228 and rs10947262 using the genotype data of Japanese populations (GWAS and the replication study), and found that they were in strong LD with each other (D' = 0.82, r 2 = 0.56). They formed two frequent haplotypes (Haplotype I and II; Table 3) accounting for about 90% of all observed chromosomes. The haplotypes were also significantly associated with knee OA; Haplotype I, the most frequent haplotype was a risk haplotype (P = 1.48610 28 ; OR = 1.33; 95% CI = 1.20-1.46).

Discussion
We performed a GWAS followed by a replication in an independent population using a total of ,4,800 Japanese subjects, and identified two SNPs (rs7775228 and rs10947262) in the HLA class II/III locus associated with susceptibility to knee OA. To our knowledge, this study represents the first GWAS of OA with extensive coverage (,550,000 markers) and definite genome-wide significance even after Bonferroni's correction, which is very conservative. There were no effects of population stratification and confounding factors. Since two groups of controls were used in the GWAS, we evaluated the possibility of genetic heterogeneity between the two groups by a principal component analysis and found it unlikely ( Figure S1). Although there was large age difference between the case and control groups of GWAS (Table 1), significant association was observed after the age adjustment.
In the NCBI genome database, rs7775228 and rs10947262 located between upstream region of HLA-DQA2 and HLA-DQB1, and within the intron 1 of BTNL2, respectively ( Figure 2). HLA-DQA2 and HLA-DQB1 encode HLA-DQ a and b chains, which belong to the HLA class II molecules. HLA class II molecules are expressed in antigen presenting cells (B lymphocytes, dendritic cells and macrophages) and play a central role in the immune system by presenting peptides derived from extracellular proteins [14]. The HLA-DQA2 protein is expressed, but at a very low level in comparison with the HLA-DQA1 protein [15,16]. Moreover, the HLA-DQA2 a chain does not dimerize with class II b chains [16]. BTNL2 encodes butyrophilin-like 2, a member of butyrophilin family that shares sequence homology with the B7 costimulatory molecules. BTNL2 regulates T-cell activation through unknown receptor, distinct from CD28 and CTLA-4 [17]. In Japanese population, the haplotype association was more significant than those of respective SNPs (Tables 2 and 3). Therefore, there may be hidden SNP(s) with a lower P value than rs7775228 and rs10947262, or the haplotype may be implicated in the OA susceptibility. An association of sarcoidosis with rs2076530, a coding SNP on exon 5 of the BTNL2 gene has been reported [18], but the SNP was not in LD with rs10947262 (D' = 0.11, r 2 = 0).
Although OA has generally been considered a non-inflammatory disease, accumulating evidences suggest that this is not the case. Inflammation involving activated T cells in the synovial membrane of OA patients is well documented [22]. Recently, we identified a genetic variant of EDG2 gene encoding lysophosphatidic acid receptor associated with knee OA [23]. A GWAS has identified a genetic variant of the PTGS2 gene encoding cyclooxygenase-2 involved in risk for knee OA [6]. These genetic associations of genes such as EDG2 and PTGS2 underscore the potential role of inflammatory pathways in the pathogenesis of knee OA.
Several studies have suggested associations of OA with HLA class I and class II alleles. Study on generalized OA revealed association with HLA A1-B8 in Caucasian [24] and with HLA-Cw4 in Japanese [21]. An association of the HLA-DRB1*02 alleles with knee and hip OA was identified in a cohort of 106 patients [25]. Interestingly, chondrocyte, which are normally HLA-DP, DQ and DR-negative, become positive for them in OA [26,27], suggesting their function as antigen-presenting cells. Cartilage fragments are mechanically shaved from the joint surface and frequently found in the synovial membrane of OA patients [28]. So, physical interaction between chondrocytes and T cells is conceivable. Peripheral blood T cells from OA patients show significantly higher proliferative responses to autologous chondrocytes [27]. Our results further support the concept that OA is an immunologic disorder. Figure S1 Principal component analysis of GWAS samples. Samples in the GWAS and in HapMap database are analyzed by a program of Smartpca [12], and plotted for the first (X axis) and the second (Y axis) principal component (PC), respectively.