Objectives Osteoarthritis (OA) is the most common form of arthritis with a clear genetic component. To identify novel loci associated with hip OA we performed a meta-analysis of genome-wide association studies (GWAS) on European subjects.
Methods We performed a two-stage meta-analysis on more than 78 000 participants. In stage 1, we synthesised data from eight GWAS whereas data from 10 centres were used for ‘in silico’ or ‘de novo’ replication. Besides the main analysis, a stratified by sex analysis was performed to detect possible sex-specific signals. Meta-analysis was performed using inverse-variance fixed effects models. A random effects approach was also used.
Results We accumulated 11 277 cases of radiographic and symptomatic hip OA. We prioritised eight single nucleotide polymorphism (SNPs) for follow-up in the discovery stage (4349 OA cases); five from the combined analysis, two male specific and one female specific. One locus, at 20q13, represented by rs6094710 (minor allele frequency (MAF) 4%) near the NCOA3 (nuclear receptor coactivator 3) gene, reached genome-wide significance level with p=7.9×10−9 and OR=1.28 (95% CI 1.18 to 1.39) in the combined analysis of discovery (p=5.6×10−8) and follow-up studies (p=7.3×10−4). We showed that this gene is expressed in articular cartilage and its expression was significantly reduced in OA-affected cartilage. Moreover, two loci remained suggestive associated; rs5009270 at 7q31 (MAF 30%, p=9.9×10−7, OR=1.10) and rs3757837 at 7p13 (MAF 6%, p=2.2×10−6, OR=1.27 in male specific analysis).
Conclusions Novel genetic loci for hip OA were found in this meta-analysis of GWAS.
- Gene Polymorphism
This is an Open Access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 3.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited and the use is non-commercial. See: http://creativecommons.org/licenses/by-nc/3.0/
Statistics from Altmetric.com
Osteoarthritis (OA) is the most common form of arthritis affecting 40% of people over the age of 70 years and is one of the most common disabling diseases observed worldwide.1 ,2 Its aetiology is multifactorial with a clear genetic component. Inheritance studies in twins and other family-based studies have assessed the estimated heritability for OA in the range of 40–65% depending on the joint site.3–5 The established OA loci have small ORs (range 1.10–1.20)6 ,7 and the genetic architecture of OA is likely to consist of multiple variants of similar magnitude.
During the last few years, extensive efforts have led to the identification of a number of OA susceptibility signals in European populations that have surpassed the genome-wide significance (GWS) level (p<5×10−8). A locus on chr7q22 near the orphan receptor GPR22 derived from genome-wide association studies (GWAS)8 ,9 and a variant in the GDF5 gene, originating from a candidate gene approach reached GWS for knee OA.10 An analysis using a 1000-genomes-project-based imputations identified a variant on chromosome 13q34 near the MCF2L gene.7 A recent GWAS in UK subjects revealed eight more loci that increase risk for OA. From those eight loci, four signals in 9q33 (ASTN2), 6q14 (FILIP1/SENP6), 12p11 (KLHDC5/PTHLH) and 12q23 (CHST11) were found to be associated with total hip replacement (THR) or hip OA in European populations.6 Finally, GWAS and functional studies revealed that the DOT1L gene on 19q13 is also associated with hip OA and cartilage thickness.11 ,12
In this study a large-scale GWAS meta-analysis for hip OA was performed under the auspices of the Translational Research in Europe Applied Technologies in Osteoarthritis (TreatOA) consortium including eight sample sets, in the discovery stage. With a total of 4349 hip OA cases and 46 903 controls in the discovery stage, and a total of 11 277 cases and 67 473 controls, this is the largest study of hip OA to date.
Study design and analysis plan
A two-stage design was used for the identification of potential associations. In the discovery stage imputed and directly genotyped autosomal single nucleotide polymorphism (SNPs) were assessed using a quality control (QC) procedure that is described in the online supplementary material (section 1.a). Briefly, we excluded SNPs based on low minor allele frequency (MAF <1%), low imputation quality, low call rate and deviation from the Hardy-Weinberg equilibrium. Genomic control was applied to each study before meta-analysis. The effect estimates of each study where synthesised using an additive model. The variants that surpassed the p<1×10−6 threshold in the meta-analysis were selected for further follow-up. Besides the main analysis including all participants a separate analysis stratified by sex was performed to detect possible sex-specific signals. In silico and de novo replication was sought for the discovery signals in 10 additional studies. All the derived effects from the discovery and the replication stage were finally synthesised using inverse variance fixed-effect models and the between-study heterogeneity was assessed using the I2metric.13 Moreover, a random-effects (RE) model was applied.14 A p-value of <5×10−8 was considered GWS. The associations of the top findings at the discovery stage were also assessed when adjusted for other risk factors such as age, height and body mass index (BMI). We also examined the association of these markers with height and BMI in the large publicly available sample set of the Genetic Investigation of Anthropometric Traits (GIANT) consortium (133 000 individuals in height analysis and 123 000 individuals in BMI analysis.15 ,16 The detailed analysis plan is presented in the online supplementary material (section 1.b).
Study populations and phenotype definition
The studies included in the discovery and replication efforts are described in table 1 and more details are given in the online supplementary material (section 1.c). All studies had standardised definitions of the phenotypes. Specifically, the definition of the hip OA in the studies was either a radiographic Kellgren and Lawrence (K/L) grade of ≥2 or history of a THR surgery because of OA. THR subjects were excluded from the study if they had: other major arthropathy (eg, rheumatoid arthritis, ankylosing spondylitis); Paget's disease affecting the pelvis or femur; THR due to hip trauma or avascular necrosis of the femoral head; or terminal illness. The control groups consisted of subjects who had no known affected joints. Population-based controls were used by the arcOGEN study.
Genotyping and imputation
Genotyping of GWAS genotyping was performed by each study following standard protocols and imputation was then carried out at the individual study level on the ∼2.5 million SNPs from HapMap Phase 2 release 22 using genome build 36 (Utah residents with ancestry from northern and western Europe (CEPH))17 on MACH or IMPUTE18 software. Imputation quality scores for each SNP were obtained from IMPUTE and MACH statistics, as appropriate. An overview of all studies and the genotyping platforms and imputation method used is given in table 1. Ten studies of European ancestral origin provided data for independent replication. Four datasets (arcOGEN stage 2, arcOGEN plus, Osteroporotic Fractures in Men study and the Study of Osteoporotic Fractures) provided ‘in silico’ replication whereas ‘de novo’ replication was performed in six other study groups (Icelandic, Swedish, Estonian Genome Center, University of Tartu (EGCUT), Paprika study, Greek, Spanish). QC criteria for deviation from Hardy-Weinberg equilibrium, MAF inconsistencies with the discovery data and outliers were applied in the replication data before including all the available data in the final analysis.
Expression was determined by Illumina HT-12 V3 microarrays using standard methods using 47 000 probes corresponding to over 25 000 well-characterised genes. DNA was available from blood and cartilage. Using the Beadstudio software the intensity values were normalised using the ‘rsn’ option in the Lumi R-package. The corresponding signals increase exponentially with relative levels and units are light intensity (Illumina provided values). The obtained raw probe-level data (overall mean normalised probe level value of measured genes in cartilage) were exported for analyses using Limma.19 As implemented in Limma, a paired t test was used on all samples. There were two probes, approximately 2.2 kb apart, on the array used for NCOA3.
Heritability of hip OA explained by genetic variants
We calculated the sibling recurrence risk and the expected genetic variance explained for hip OA hits that were identified previously and in this study as described in online supplementary material (section 1.e)
The final analysis included a total of 11 277 radiographic and symptomatic hip OA cases and 67 473 controls of European ancestry, with 4349 cases and 46 903 controls included in the discovery stage and 6928 cases and 20 570 controls in the follow-up effort (table 1). In the sex-specific analyses 2045/20 823 male cases and controls and 2689/25 384 female cases and controls were analysed in the discovery. After QC, 2 567 279 SNPs were analysed. Low genomic inflation factor was observed for the first stage with λ=1.028. The results of the discovery stage were uncorrected for the overall inflation factor. Quantile-Quantile plots showed an excess of signals compared with what was expected by chance, indicating the presence of true association signals that could confer susceptibility to hip OA (see online supplementary figures S1–S3). For male-specific analysis (see online supplementary figure S3) there is an early deviation from the neutrality line, therefore positive signals should be treated with caution.
Following analysis of the discovery stage it was found that eight independent loci reached the prespecified threshold of p<1.0×10−6 required for further replication; five from the combined analysis of sexes, two male-specific and one female-specific loci. The number of independent SNPs is larger compared with three estimated independent SNPs expected under the null for the main and the sex-specific analyses (binomial test p=0.012). Two of these signals (rs6094710 at 20q13 and rs640070 at 11q25) were imputed, had MAF <5% and moderate or large effect sizes were observed (OR >1.2). Therefore, we examined for possible imputation errors by comparing with de novo genotyping in a random sample of the RS-I and the Twins UK studies to exclude any chance of false positive findings. The obtained MAF estimates between the imputation and the genotyping efforts were not consistent for rs640070, which was therefore excluded from further consideration, minimising the chance of any false positive signals (see online supplementary material; section 1.d). Moreover, we included in the replication stage rs17610181 at 17q23, a SNP that was just below the desired threshold, but was previously shown to be associated with height.15 Therefore, eight independent SNPs were included in the replication phase. The combined effect sizes and p values for the eight signals are presented in table 2 and in online supplementary figures S4–S11. None of these variants associate with height or BMI in the large publicly available GIANT databases (see online supplementary table S2) nor did including these covariates, or age, in the analysis change the association of OA with these genetic variants (see online supplementary table S3).
Two SNPs, rs6094710 at 20q13 and rs1577792 at 6q14, were borderline GWS in the discovery stage with p=5.8×10−8 and 7.3×10−8, respectively, with no observed heterogeneity (I2=0). The SNP rs6094710 replicated with p=7.3×10−4 in the follow-up samples and reached GWS level with p=9.3×10−9 and OR=1.28 (95% CI 1.17 to 1.39) when we combined the discovery and the replication data, although large heterogeneity was observed I2=64% (figure 1). The SNP remained GWS after the second genomic control. rs6094710 is annotated near NCOA3 (nuclear receptor coactivator 3) gene. The OR of rs1577792 at 6q14 was close to unity in the replication effort with one study being significant in the opposite direction.
The rs3757837 SNP at 7p13 from the male-specific analysis replicated nominally (p=0.044, OR=1.15) in the follow-up samples, but did not reach GWS in the overall analysis with p=2.2×10−6 and OR=1.27. Moderate heterogeneity was observed for this SNP in the discovery analysis (I2=46%). Heterogeneity was increased, in the overall analysis with all the studies combined (I2=78%), reflecting the further heterogeneity introduced by the replication data. We, therefore, also applied the Han and Eskin (RE) model, an approach that allows more heterogeneity in the data compared with traditional models.14 Using this model rs3757837 showed a stronger association with p=7.7×10−10 in this analysis. The strength of the association of this approach was not substantially different from the fixed-effect model for any other SNP in our study (table 2). rs3757837 resides in intron 8 of the CAMK2B (calcium/calmodulin-dependent protein kinase II β) gene.
The rs5009270 SNP at 7q31 remained suggestively associated in the final analysis (p<5×10−6) with combined two-stage p=9.9×10−7 and summary OR 1.10. rs5009270 is located near the IFRD1 (interferon-related) gene. The p with the RE approach was 3.1×10−6.
Gene expression data
To further investigate our findings for the top GWS hit at 20q13 we explored mRNA expression profile of the NCOA3 gene in articular cartilage, a highly relevant tissue for OA. Expression of the gene was examined by exploring a micro array mRNA expression dataset generated on Illumina V3 Human-12 chips in cartilage samples of 33 patients (13 men and 20 women of European descent aged 54 years to 80 years) that underwent joint replacement due to end-stage OA disease. Expression levels of NCOA3 in cartilage displaying symptoms of OA were compared with expression levels in cartilage that appeared macroscopically normal but isolated from the same joint (preserved cartilage). A moderate level of expression, as determined by the mean normalised probe level value, was observed for NCOA3 (mean level 8.49), which was above the observed average expression of genes in the articular cartilage (mean normalised probe level value of measured genes in cartilage was 7.4; range 6.6–14.9). Expression was also high in blood (6.96; range 6.3–14.7). When we tested for differential expression of NCOA3 among the pairs of preserved and OA-affected cartilage, we observed a significantly lower expression of NCOA3 in the OA-affected cartilage (p=0.0064).
Heritability of hip OA explained by genetic variants
Table 3 summarises genome-wide significant and suggestive signals of hip OA including this study. Based on these findings and if we consider a sibling recurrence ratio λs=520 then the discovered signals of OA contribute 3% of the heritability in OA when all hits are considered and 2.1% if only GWS are evaluated.
In this report we attempted to further clarify the genetic architecture of the genetic background of hip OA by using the largest sample-size for hip OA to-date of more than 78 000 genotyped individuals under a GWAS framework. During the discovery phase, we identified eight signals that qualified for further independent replication. Of these, the signal at chromosome 20q13 near NCOA3 gene was found to be GWS and two other loci were suggestively significant at a p<5×10−6 level in the joint analysis of the discovery and replication stages. Adjusted analyses of the prioritised signals revealed that these markers were not associated with body size. These genetic risk factors contribute to our knowledge base in the field of the susceptibility for hip OA by conferring a medium OA risk.
The top signal identified in this meta-analysis was rs6094710, a variant that is annotated on chromosome 20q13 near the NCOA3 gene, increasing the risk for hip OA for the carriers of the A allele by almost 30%. Furthermore, we showed that the identified NCOA3 gene was expressed in articular cartilage and its expression was significantly reduced in OA affected cartilage, further supporting a role of the NCOA3 gene signal in OA disease process. Interestingly, rs6094710 is in complete linkage disequilibrium with rs6094752 (r2=1) which is a missense SNP leading to an amino acid change at position 218 in the protein (Arg>Cys). This amino acid change is predicted to have a benign and damaging effect on the protein by PolyPhen-2,21 dependent on the variant protein. The functional consequences of SNP rs6094752 are unknown, and further research is therefore needed to unravel the biological mechanism of this amino acid change in relation to OA.
The NCOA3 gene is a nuclear receptor coactivator that directly binds nuclear receptors and stimulates the transcriptional activities in a hormone-dependent fashion. In this signalling process, NCOA3 recruits histone acetyltransferases and methyltransferases for chromatin remodelling and facilitating downstream gene transcription. NCOA3 is involved in the coactivation of different nuclear receptors, such as for steroids, retinoids, thyroid hormone, vitamin D3 and prostanoids. Many of these hormones have been implicated in skeletal metabolism and OA, which makes NCOA3 a compelling causal candidate gene. Previously, NCOA3 knockout mice were generated through homologous recombination in embryonic stem cells.22 These mice showed growth retardation and reduced adult body size, but the molecular mechanism responsible for this growth retardation remains largely unknown. In addition, female mice exhibited abnormal development and function of their reproductive system and oestrogen levels were significantly lower in the knockout mice compared with the wild type,22 possibly indicating involvement of NCOA3 in steroid regulation.
NCOA3 could be also implicated through regulation of the target tissue responses to thyroid hormone (T3).23 Since intracellular T3 is tightly regulated by deiodinase, iodothyronine, type-2 and deiodinase, iodothyronine, type-3 encoded by the DIO2 and DIO3 genes, respectively, that were previously recognised as OA susceptibility genes,24 ,25 the current NCOA3 findings complement the previous outlined hypothesis that local T3 signalling may affect OA susceptibility.26
Another possible mechanism by which NCOA3 might be involved in cartilage homoeostasis is through transcriptional regulation in mechanotransduction. NCOA3 is upregulated by the signal transducer and activator of transcription 6 (STAT6) in naïve splenic B cells from BALB/c and serves as a positive regulator of transcriptional activation by STAT6.27 STAT6 is known to be the common signal transducer of interleukin (IL)-4 receptor α chain and mediates IL-4- and IL-13-induced responses.28 Chondrocytes from normal and OA cartilage signal through a type II IL-4R in in human articular chondrocyte mechanotransduction. This signalling is via a STAT6-independent pathway. Differences in IL-4 signalling are likely due to crosstalk between integrin and cytokine signalling pathways.29 Therefore, NCOA3 may be related to cartilage function and molecular signalling and transcriptional regulation in mechanotransduction.
The male-specific locus on 7p13, represented by rs3757837, showed considerable heterogeneity between studies. The signal was strongly supported by the RE model. Unlike the conservative traditional RE methods this new method has been shown to achieve higher statistical power when heterogeneity exists, allowing for new discoveries in the field of genetic epidemiology.14 rs3757837 is located in CAMK2B gene, which belongs to the calcium/calmodulin-regulated kinase (CaMKII) subfamily. There is evidence that CaMKII-signalling may be important in onset and progression of OA.30 ,31 This pathway has been described as central to the molecular events that regulate chondrocyte responses to mechanical stimulation and in particular, to the upstream effect of IL-4. This would be linked to the STAT6 pathway, which in turn is regulated by NCOA3. Thus these two genetic signals, NCOA3 and CAMK2B, both point to pathways related to cartilage mechanotransduction, suggesting that genetic defects in this pathway may be central to the degeneration of cartilage that takes place in hip OA. OA is a disease affecting articular cartilage and the underlying bone, resulting from many biological and mechanical interacting factors which change the extracellular matrix and cells and lead to increasing levels of cartilage degeneration. Joint tissues are exquisitely sensitive to their mechanical environment, and mechanical loading may be the most important external factor regulating the development and long-term maintenance of joint tissues.32 Finally, rs5009270 resides near IFDRD1 which codes an interferon-related developmental regulator and has been implicated in skeletal muscle regeneration.33 Reduction in muscle strength is strongly associated with functional decline, and individuals with lower quadriceps strength adjusted for body weight are more likely to develop OA.34
Our study has certain limitations. In our discovery stage, two out of the eight independent loci had no heterogeneity and in five signals the heterogeneity was moderate or low. However, when we included the replication data, large heterogeneity was introduced for all SNPs. Conflicting results in the replication data, besides chance, could be explained by inconsistent definitions of the OA phenotypes but also from the different population structure that can introduce heterogeneity. The TreatOA consortium has addressed the need of standardisation of OA phenotypes35 and this effort may have diminished the observed heterogeneity in the discovery stage. Efforts including more data using common and stringent QC criteria and standardisation methods should substantially improve the power of GWAS to identify novel findings in the near future. The main limitation of the expression study is that there are few individuals included, in particular only two carriers of the rs6094710 variant.
In conclusion, novel loci involved in hip OA were discovered through a large-scale meta-analysis of GWAS. The exact underlying mechanism leading to a higher risk of OA remains to be elucidated by functional experiments. It is evident from this work and other recent studies that deciphering the architecture of the genetics of OA requires major large-scale efforts, and in this regard calls for international worldwide collaborations are not fruitless. In the near future emphasis should be given to the enhancement of the total sample size, the adoption of stringent and standardised definitions of the phenotypes and the application of imputation-based meta-analysis using the panel of the 1000 Genomes Project.36 In addition, linking results from genetic association studies to for example genome-wide RNA expression data might further improve our understanding. Eventually, large-scale studies with whole-genome sequencing will be needed to target the heritability caused by less common variants.
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Files in this Data Supplement:
- Data supplement 1 - Online supplement
Handling editor Tore K Kvien
EE, HJK, US, EEN, EZ, IM, JPAI, TDS, JBvM and AMV contributed equally.
Correction notice This article has been corrected since it was published Online First. The affiliations for P Eline Slagboom have been corrected and affiliations 1 and 2 have been amended.
Acknowledgements The authors are very grateful to all study participants, the staff from all studies and the participating physicians and pharmacists. Acknowledgements for each study are given in online supplementary material (section 4).
Contributors EE, HJK, US, IM, JPAI, TDS, JBvM and AMV conceived and designed the study. EE, HJK, EEN and KKT analysed the data. EE, HJK, US, EEN, EZ, IM, JPAI, TDS, JBvM and AMV wrote the first draft of the manuscript. All authors interpreted the results, critically commented and approved the final version of the manuscript.
Funding This project was supported through Coordination Theme 1 (Health) of the European Community's FP7 (grant 200800) TreatOA. The work was also supported by EU FP7 Small-scale focused research collaborative projects EurHEALTHAging (grant 277849).
Competing interests US, IJ, GT, HTH, UT and KS are employed by deCODE Genetics/Amgen. HJK is currently employed by Pfizer.
Ethics approval Meta-analysis of GWAS, different review boards approved the individual studies.
Provenance and peer review Not commissioned; externally peer reviewed.
If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.