Article Text
PDF
Abstract
Objectives Osteoarthritis (OA) is the most prevalent form of arthritis and accounts for substantial morbidity and disability, particularly in older people. It is characterised by changes in joint structure, including degeneration of the articular cartilage, and its aetiology is multifactorial with a strong postulated genetic component.
Methods A meta-analysis was performed of four genome-wide association (GWA) studies of 2371 cases of knee OA and 35 909 controls in Caucasian populations. Replication of the top hits was attempted with data from 10 additional replication datasets.
Results With a cumulative sample size of 6709 cases and 44 439 controls, one genome-wide significant locus was identified on chromosome 7q22 for knee OA (rs4730250, p=9.2×10−9), thereby confirming its role as a susceptibility locus for OA.
Conclusion The associated signal is located within a large (500 kb) linkage disequilibrium block that contains six genes: PRKAR2B (protein kinase, cAMP-dependent, regulatory, type II, β), HPB1 (HMG-box transcription factor 1), COG5 (component of oligomeric golgi complex 5), GPR22 (G protein-coupled receptor 22), DUS4L (dihydrouridine synthase 4-like) and BCAP29 (B cell receptor-associated protein 29). Gene expression analyses of the (six) genes in primary cells derived from different joint tissues confirmed expression of all the genes in the joint environment.
Statistics from Altmetric.com
Introduction
Osteoarthritis (OA) is the most prevalent form of chronic joint disease and accounts for substantial morbidity and disability, particularly among older people. It is characterised by loss of joint homeostasis. The articular cartilage cannot maintain its integrity and is progressively damaged, the subchondral bone envelope is thickened changing loads in the bone-cartilage biomechanical unit, the synovium shows signs of inflammation and bony spurs (osteophytes) appear at the edges of the bone. Its aetiology is multifactorial with a significant genetic component as shown by twin and family studies.1 2
Many genetic variants have been considered as potential risk factors for OA, but most of the reported associations are inconclusive or not replicated. A recent large-scale meta-analysis found evidence that the GDF5 locus on chromosome 20 was associated with the increased risk of knee OA in Caucasians.3,–,6 Other genome-wide data have reported an association with the DVWA gene in Asians but not Caucasians7 and a PTGS2 variant that replicated but did not reach genome-wide significance (GWS).8 Recently, a genome-wide association (GWA) study identified a locus on chromosome 7q22 which has an association with combined knee OA and/or hand OA phenotype.9
In this study we have synthesised available data from four GWA studies under the auspices of the Translational Research in Europe Applied Technologies for Osteoarthritis (TreatOA) consortium (www.treatoa.eu). A total of 2371 cases of knee OA were available for this first stage of the analysis. The most significant signals were further investigated in additional samples of European descent and single nucleotide polymorphisms (SNPs) that reached GWS were further evaluated in Asian samples.
Methods
Study design
A detailed description of all samples used in this study is provided in the online supplement. A three-stage design was used for the identification of any potential associations between sequence variants and knee OA in populations of European ancestry. We first synthesised the available data from four GWA studies (deCODE, Rotterdam Study, Framingham, Twins UK) using inverse variance fixed effects models. The variants that reached the 2×10−5 level of significance were selected for further replication. These SNPs were followed up in eight additional European cohorts (arcOGEN, Greek, Spanish, Finnish, Nottingham, Chingford study, GARP, Estonian and Swedish). The SNPs that replicated in the follow-up samples were genotyped in two additional European samples (deCODE (Icelandic) and Swedish). One cohort provided computer-generated replication from an ongoing GWA study (arcOGEN, 12 SNPs were directly genotyped and 6 were imputed) while de novo replication was performed in the other cohorts. Furthermore, the top hits were followed up in Asian populations (Chinese and Japanese samples). The effect sizes from the meta-analysis of the GWA studies and the effect sizes from the replication effort were all combined to provide an overall estimate. We also synthesised the effect estimates of the European and Asian samples to provide a global summary effect estimate.
Phenotype definitions
Study subjects with a radiographic Kellgren and Lawrence (K/L) grade ≥210 or total knee replacement were included as cases in the analysis. When clinical criteria were considered (Greek, Spanish and GARP study groups), the American College of Rheumatology classification criteria were used.11 Subjects who had no known affected joints among those assessed acted as controls. For example, in a cohort that assesses knee, hip and hand OA, controls were participants with no affected hip or hand joints for the knee OA analysis. Population-based controls were used for the arcOGEN study.
Genotyping and imputation
Samples from the GWA studies were genotyped using the Infinium HumanHap300 (Illumina) for deCODE and Twins UK samples, HumanHap550v3 Genotyping BeadCHip (Illumina) for the Rotterdam Study and the Affymetrix GeneChip Human Mapping 500K for the Framingham cohort. The number of SNPs genotyped ranged from 314 075 to 500 510. Imputations were performed to increase the coverage. All the top SNPs studied had acceptable imputation quality. The genotyped and imputed SNPs that successfully passed the quality control criteria (n=2 335 627) were considered for the analyses. Detailed information on genotyping platform, quality control and imputation methods for each cohort are shown in table S1 in the online supplement.
The replication samples for the Greek, Spanish, Finnish, Chingford and GARP studies were genotyped using the MassArray iPlex Gold from Sequenom. Replication genotyping was carried out by a genotyping contractor (Kbiosciences Ltd, Hertfordshire, UK) using a competitive allele-specific PCR SNP genotyping system for the Nottingham and the Estonian cohort. The additional 622 Icelandic cases and the samples from the Swedish cohort were genotyped by deCODE genetics using the Centaurus (Nanogen) platform.12 Detailed information on genotyping is provided in the online supplement.
Statistical analysis
Association analysis
Each team performed an association test per gender for knee OA under a per-allele model. The λ inflation factor was calculated per gender-specific effect size using the genomic control method13 and the standard errors were corrected by the square root of the λ inflation factor was calculated per gender-specific effect size using the genomic control method13 and the standard errors were corrected by the square root of the λ inflation factor
Meta-analysis
The effect size for each SNP (OR per copy of minor allele as per HapMap) was calculated using inverse variance fixed effects models,14 synthesising all the sex-specific effect sizes and the corrected standard errors. Analyses combining men and women were also performed. In family studies the results from men and women combined were used to account for relatedness between women and men within families. Meta-analyses of the GWA studies were performed using the METAL software (www.sph.umich.edu/csq/abecasis/metal). Between-study heterogeneity was tested using the Cochran Q statistic, which is considered significant at p<0.1. The extent of inconsistency across studies was quantified using the I2 metric which ranges from 0 to 100%.15 Heterogeneity is considered low, moderate, high and very high for 0–24%, 25–49%, 50–74% and >75%, respectively.16 We also computed the 95% CI for the I2.17 The calculation was repeated with random effects models for all SNPs that were further evaluated in replication datasets. Meta-analyses of the 18 top hits were performed using Stata Version 10.1.
Assessment of credibility
In order to assess the credibility of the top hit, we calculated the Bayes factor under a spike and smear prior to using as an alternative an average genetic effect corresponding to an OR of 1.2 and a conservative agnostic prior of 0.0001%.18
Functional analysis
Two methodological approaches were used to investigate the functional role of genes identified by GWA studies: (1) by assessing their expression in primary human joint cells (synovial fibroblasts, chondrocytes and meniscal cells) and its change in response to the proinflammatory cytokines tumour necrosis factor α and interleukin 1β as well as comparing their gene expression profiles during chondrocyte dedifferentiation (3D pellet cultures vs monolayer culture); and (2) by assessing their expression dynamics by whole mount in situ hybridisation using zebrafish (Danio rerio) embryos aged 6 h (shield), 10 h (bud), 13 h (5–9 somites) and 1, 2, 3 and 4 days to explore their role during embryogenesis.
Results
Meta-analysis of GWA studies and replication of top findings
The descriptive characteristics of the GWA studies used for the meta-analyses are from Iceland (deCODE), the Netherlands (Rotterdam study), USA (Framingham) and the UK (Twins UK). The characteristics of these studies are shown in table 1 and in the online supplement. The four GWA datasets included a total of 2371 cases and 35 909 controls. A quantile-quantile plot comparing the meta-analysis association results of the four studies with those expected by chance showed an excess of SNP associations indicating a likely true association signal (figure 1). Data analysis showed the strongest association on chromosome 7q22 with a p value of 5.06×10−8 for rs4730250 localised in dihydrouridine synthase 4-like gene (DUS4L) (figure 2). Other associated signals in the 7q22 gene cluster were in high linkage disequilibrium (LD) (r2>0.8) with the top signal (figure 2).
Quantile–quantile plot of the expected versus observed distribution of p values.
Regional association plot of rs4730250. Statistical significance of the associated SNPs are illustrated on –log10 scale. The p value of the rs4730250 and the other 10 selected SNPs are based on the meta-analysis of all datasets (both genome-wide association (GWA) studies and replication studies); p values for the other SNPs are based on the meta-analysis of the GWA studies. The sentinel single nucleotide polymorphism (SNP) is shown in blue. The correlation of the sentinel SNP is shown on a scale from minimal (gray) to maximal (red). SNPs in red have r2≥0.8 with the sentinel SNP and SNPs in orange have r2≥0.5. Chromosome positions are based on HapMap release 22 build 36.
Characteristics of the studies included in the analysis
We selected for follow-up in replication samples all SNPs with a p value <2×10−5 in the meta-analysis association results. A total of 18 SNPs from 10 chromosomal loci satisfied this criterion (see table S2 in online supplement). However, as some of those SNPs were fully equivalent in the HapMap-CEU dataset, a total of 11 non-identical SNPs were tested for replication in 3326 cases and 7691 controls from eight European studies (see table 1 and online supplement). Two SNPs (rs4730250 and rs10953541), both located at 7q22, replicated nominally (p<0.05) in the combined analysis of the follow-up samples with p values of 6.3×10−4 and 8.3×10−3, respectively. The two SNPs rs4730250 and rs10953541 were then further genotyped in two additional replication sets.
Both SNPs reached GWS in a meta-analysis of all European sample sets (GWA datasets and replication cohorts, table 2). A total of 6709 cases of knee OA cases and 44 439 controls were analysed. SNP rs4730250 was genome-wide significant with a per-allele summary OR of 1.17 (95% CI 1.11 to 1.24) and a p value of 9.2×10−9. The minor allele frequency was 0.17 in the combined dataset. Low heterogeneity was observed (I2=15%, 95% CI 0% to 48%) which was not statistically significant (p=0.26 for Cochran Q statistic, figure 3). No gender-specific effects were seen. The summary estimates did not differ significantly in men and women (p=0.74, test of homogeneity, figure 3). Analysis of both sexes together in all cohorts did not alter the results (OR 1.17, 95% CI 1.07 to 1.27, p=4.1×10−8). The summary effect sizes of all loci under study are shown in table 2 and the results from the random effects analysis for the top hits are shown in table S3 in the online supplement.
Forest plot of study-specific estimates (black boxes) and summary OR estimates and 95% CIs (diamonds) for the association between the rs4730250 single nucleotide polymorphism and osteoarthritis of the knee.
Summary OR and 95% CI of SNPs in the analysis including all European descent data
The two significant SNPs at 7q22, rs4730250 and rs10953541, are highly correlated (D′=1, r2=0.63 in HapMap-CEU) and are likely to represent the same underlying association signal as shown by conditional association analysis (see table S4 in online supplement). Age and body mass index are considered to be significant risk factors for the development of knee OA.19,–,25 We performed an analysis where the top hit was adjusted for these risk factors in deCODE samples and the Rotterdam study. The association of the top hit remained largely unchanged in analyses adjusted for body mass index and age.
In order to assess the credibility of the associations of the two SNPs, we calculated the Bayes factor18 under a spike and smear prior using an average genetic effect corresponding to an OR of 1.2 and a conservative agnostic prior (assuming no prior knowledge of the association) of 0.0001%. The posterior credibility of these associations was 98% and remained similarly high even with a small alternative effect size of 1.1.
We also tested if the observed signal at the 7q22 region was replicated in East Asian samples (Japanese and Chinese cohorts). The total numbers of cases of knee OA and controls assessed were 1183 and 1245, respectively. rs12535761 was used as a proxy for rs4730250. The two SNPs are in strong LD (r2=1, D′=1 in HapMap Asian samples). The finding was not replicated in the Asian samples with a summary effect size of 1.03 (95% CI 0.85 to 1.25). A meta-analysis including both European and Asian samples with 7892 cases and 45 684 controls yielded a global summary effect of 1.15 (95% CI 1.10 to 1.22) with a p value of 5.7×10−8 for rs4730250 with low heterogeneity (I2=19%).
Expression patterns of genes in 7q22 cluster
The associated signal at 7q22 is located within a large (500 kb) LD block which contains six genes: PRKAR2B (protein kinase, cAMP-dependent, regulatory, type II, β), HPB1 (HMG-box transcription factor 1), COG5 (component of oligomeric golgi complex 5), GPR22 (G protein-coupled receptor 22), DUS4L (dihydrouridine synthase 4-like) and BCAP29 (B cell receptor-associated protein 29).
We performed additional experiments to get more information about the genes in the cluster and their potential role in joint biology and pathology. Analysis of mRNA expression data in a chondrocyte pellet indicates that BCAP29, COG5, DUS4L and HBP1 expression levels were higher than in monolayer cultures, suggesting that they are expressed in an environment that more accurately recapitulates articular cartilage (see figure S1 in online supplement). In contrast, no difference was seen for GPR22 and PRKAR2B mRNA expression. In a zebrafish model, the expression of all genes was detectable from the shield stage onwards (see detailed results and figures S2 and S3 in the online supplement).
Discussion
This study provides further evidence for a knee OA signal localising to the 7q22 cluster region and associated with knee OA. The statistical credibility and confidence of this evidence is very high, based on the calculations of the Bayes factor. The same locus has been identified and proposed as an OA susceptibility locus from the Rotterdam study for the prevalence and progression of OA.9 Our study and the earlier Rotterdam study do include overlapping populations. However, our study was specifically targeting the knee OA phenotype. An additional three European cohorts and two Asian populations were used for further replication. Our study uses the largest sample size in the genetics of knee OA research to date with almost 8000 cases of knee OA analysed.
The most significant hits identified by our study are located within a large (500 kb) LD block that contains six genes: PRKAR2B, HPB1, COG5, GPR22, DUS4L and BCAP29. The top hit rs4730250 is annotated in intron 3 of the DUS4L gene. Any of the genes at the 7q22 region may confer risk for knee OA as the LD pattern across the region is high.
The gene expression data support the epidemiological findings but do not exclude any of the six candidate genes. Specifically, the zebrafish experiments show that both COG5 and DUS4L are expressed in developing cartilage, supporting the notion that either of these genes could have a biological function during chondrogenesis. The studies in the dedifferentiation model of human chondrocytes (3D vs 2D culture) show that BCAP29, COG5, DUS4L and HBP1 all have different expression patterns in 3D culture (chondro-like cells) from 2D culture (dedifferentiated cells), suggesting that these four genes may play a role in cartilage metabolism.
A major issue in the field of OA is the definition of the disease phenotypes.4 26 Different criteria may introduce bias and dilute the effect. The cases in our study were defined either clinically by the presence of a knee replacement or radiographically using the K/L system. The K/L system is, however, far from perfect and can be affected by differences in the position of the knee in which the x-rays were obtained, observer biases, interpretation of grading criteria and random error.27 28 Similarly, there are no standard criteria for replacing knee joints. This may introduce heterogeneity and move the observed effects towards the unity and so underestimate the true strength of an association. In our study we synthesised data with a standardised definition of the phenotype; however, small individual locus effects with ORs in the range of 1.1–1.2 as for other chronic diseases may well be plausible for knee OA, explaining the paucity of other significant hits despite the reasonable large-scale effort. These findings highlight that even larger collaborative studies and improved standardisation of the phenotypes are needed to better understand and identify further genetic variants of OA.
Moreover, even though we were able to accumulate a large sample size, the power of the study to detect very small effect sizes in the range of 1.05–1.15 is inadequate. For example, identification of a GWS signal with an effect size of 1.15 and minor allele frequency of 20% with 80% power would require almost 7000 additional cases of knee OA.
Our results confirm that the 7q22 chromosomal region confers risk for knee OA which, along with our functional work, implicates six possible genes. Further in-depth genetic analysis of the locus including deep sequencing of the region and functional work including in vitro assays and animal models will be required to deepen our understanding of the underlying molecular pathways associated with the disease.
Acknowledgments
The authors thank all the treatOA participants for their contribution in the study. TreatOA is funded by the European Commission framework 7 programme (grant 200800). The authors thank all arcOGEN participants for their contribution to this manuscript. arcOGEN is funded by a special purpose grant from the Arthritis Research Campaign (arc, grant 18030). This study used genotype data from population controls that was generated by the Wellcome Trust Case Control Consortium 2 (http://www.wtccc.org.uk) funded by The Wellcome Trust (grant 083948). The population controls were from the 1958 British Birth Cohort collection funded by the Medical Research Council (grant G0000934) and The Wellcome Trust (grant 068545) and from the UK Blood Services Collection of Common Controls funded by The Wellcome Trust. The samples used in arcOGEN derive from five centres in the UK: Nottingham, London, Oxford, Sheffield and Southampton. For Nottingham we acknowledge arc for funding the collection of the majority of cases. For London we thank the staff from the TwinsUK unit and the Chingford Study for patient ascertainment, we acknowledge financial support from arc, from the Wellcome Trust and from the Department of Health via the National Institute for Health Research (NIHR) comprehensive Biomedical Research Centre (BRC) award to Guy's & St Thomas' NHS Foundation Trust in partnership with King's College London. For Oxford we acknowledge funding support from the Collisson Foundation, the Botnar Foundation and the Jean Shanks Foundation for patient ascertainment, we acknowledge the NIHR for supporting the Biomedical Research Unit (BRU) at the University of Oxford, and we thank Bridget Watkins and Kim Clipsham for assistance in patient ascertainment. For Sheffield we acknowledge the NIHR for supporting the Sheffield Bone BRU, the South Yorkshire Clinical Research Network for part funding the Sheffield research nurse and for clerical support, the Royal College of Surgeons of England and the Cavendish Foundation. For Southampton we acknowledge the Wellcome Trust Clinical Research Facility at Southampton General Hospital and we thank Phillippa-Kate Battley and Elizabeth Arden for assistance with patient ascertainment, and Richard Keen and Anna Bara, principal investigator and trial manager for the arc-funded VIDEO study, respectively. We acknowledge the support of the UK NIHR BRC for Ageing and Age-related disease award to the Newcastle upon Tyne Hospitals NHS Foundation Trust (JL and AMcC). We acknowledge sample management undertaken by the UK DNA Banking Network funded by the Medical Research Council at the Centre for Integrated Genomic Medical Research (CIGMR), University of Manchester and we thank Kate Dixon, Kate Sherburn and Debbie Payne for their assistance. Genotyping was performed at the Wellcome Trust Sanger Institute and we thank Emma Gray, Sarah Edkins, Rhian Gwilliam, Suzannah Bumpstead and Cordelia Langford for their assistance. Analysis of the arcOGEN data was performed at the Wellcome Trust Centre for Human Genetics and at the Wellcome Trust Sanger Institute and we acknowledge the work of the arcOGEN analysis team members Nigel W Rayner, Lorraine Southam, Guangju Zhai, Katherine S Elliott, Sarah E Hunt, Hannah Blackburn, Simon C Potter, Aaron Garth Day-Williams and Claude Beazley. EZ is supported by the Wellcome Trust (WT088885/Z/09/Z), LS is supported by the European Community Framework 7 large collaborative project grant TREAT-OA, KC is supported by a Botnar Fellowship and by the Wellcome Trust (WT079557MA), NWR is supported by the Wellcome Trust (WT079557MA), JMW is supported by the Higher Education Funding Council for England. RJL is the recipient of a postdoctoral fellowship from the Flanders Research Foundation (FWO Vlaanderen). ROAD (TA, HK, AM, NN, NY) acknowledge Katsushi Tokunaga, Shigeyuki Muraki, Hiroyuki Oka and Kozo Nakamura for scientific advice and data collection. We acknowledge funding support by Grants-in-Aid for Scientific Research (S19109007, B21390417) from the Japanese Ministry of Education, Culture, Sports, Science and Technology, H17-Men-eki-009 from the Ministry of Health, Labor and Welfare, and JOA-Subsidized Science Project Research 2006-1 from the Japanese Orthopaedic Association.
References
Supplementary materials
Web Only Data
Files in this Data Supplement:
Footnotes
-
Funding European Commission framework 7 programme grant 200800 TREAT-OA, NWO Investments (175.010.2005.011).
-
Ethics approval Each study participating in this meta-analysis has obtained approval from respective ethics committe.
-
Provenance and peer review Not commissioned; externally peer reviewed.