Article Text
PDF
Abstract
Objective Osteoarthritis (OA) is the most common form of arthritis and the leading cause of disability in the elderly. Of all the joints, genetic predisposition is strongest for OA of the hand; however, only few genetic risk loci for hand OA have been identified. Our aim was to identify novel genes associated with hand OA and examine the underlying mechanism.
Methods We performed a genome-wide association study of a quantitative measure of hand OA in 12 784 individuals (discovery: 8743, replication: 4011). Genome-wide significant signals were followed up by analysing gene and allele-specific expression in a RNA sequencing dataset (n=96) of human articular cartilage.
Results We found two significantly associated loci in the discovery set: at chr12 (p=3.5 × 10−10) near the matrix Gla protein (MGP) gene and at chr12 (p=6.1×10−9) near the CCDC91 gene. The DNA variant near the MGP gene was validated in three additional studies, which resulted in a highly significant association between the MGP variant and hand OA (rs4764133, Betameta=0.83, Pmeta=1.8*10−15). This variant is high linkage disequilibrium with a coding variant in MGP, a vitamin K-dependent inhibitor of cartilage calcification. Using RNA sequencing data from human primary cartilage tissue (n=96), we observed that the MGP RNA expression of the hand OA risk allele was significantly lowercompared with the MGP RNA expression of the reference allele (40.7%, p<5*10−16).
Conclusions Our results indicate that the association between the MGP variant and increased risk for hand OA is caused by a lower expression of MGP, which may increase the burden of hand OA by decreased inhibition of cartilage calcification.
- genetics
- hand osteoarthritis
- genome-wide association study
- functional study
- Matrix-Gla protein
Statistics from Altmetric.com
Introduction
Osteoarthritis (OA) is the most frequent joint disorder worldwide. An estimated 22% of the adult population has a joint affected by OA and this incidence increases to 49% in individuals over 65 years of age.1 All synovial joints can be affected by OA, with hand OA as one of the most common forms of OA. Hand OA is characterised by osteophyte formation, bony enlargements of finger joints and cartilage degradation in the joints. One of the factors contributing to cartilage degradation is the increase of calcified cartilage in the joint.2 3 In addition, hand OA is related to the occurrence of OA at other sites, most notably with knee OA.4 5 Patients affected by hand OA suffer from pain and disability, impacting their quality of life. OA is a leading cause of chronic disability,6 yet currently no effective therapeutic treatments against OA are known. It is therefore imperative to dissect the underlying mechanism of disease aetiology as this may enhance effective and targeted drug development.
OA has a strong genetic component. Depending on the joint affected, the heritability of OA is estimated in the range of 40%–60%,7 8 with hand OA having the largest heritability, that is, ~60%.9 10 Therefore, in recent years, several large-scale genetic studies have been performed to identify the underlying genes and pathways leading to OA. Multiple significantly associated loci for OA of the hip and knee have been identified through genome-wide association studies (GWAS).11–18 However, thus far, only one report has described a robust association with OA of the hand.19 In this previous report, common variants in the ALDH1A2 and rare variants in chromosome 1p31 were genome-wide significantly associated with hand OA using a discovery cohort of 837 cases and 77 325 controls.
In this study, we aimed to identify novel genes and pathways involved in the aetiology of OA of the hand by performing a large-scale genome-wide association study (GWAS). We used a semiquantitative measure for OA of the hand in order to increase statistical power. We gathered a large sample size of 12 754 individuals for analysis, by combining data from three studies in the discovery phase and an additional three cohorts for replication. Next, we conducted functional follow-up of our top finding to investigate the underlying mechanism.
Methods
Discovery GWAS, replication and meta-analysis
For a detailed description on the GWAS methods, participating studies, quality control procedures for genotyping and imputation, see online supplementary text S1 and table S1.
Detailed phenotype description of Kellgren and Lawrence sum score
We have used a semi-quantitative bilateral measure of OA of the hand based on the radiographic Kellgren and Lawrence (KL) score.20 Using radiographs of both hands, the KL score was determined for each joint in the hand. Using these KL scores we defined the KL sum score: the total KL score, the sum, of the following hand joints for both hands (left and right): all distal interphalangeal joints, all proximial interphalangeal joints, all metacarpophalangeal (MCP) joints, the interphalangeal joint and the first carpometacarpal joint, which gives the sum of 15 joints on each hand, and in total 30 joints for both hands together, resulting in a minimum score of 0 and a maximum score of 120. In the Leiden Studies (LS) cohort, no KL scoring was done of the MCP joints, resulting in a KL sum score of maximum 88. Individuals lacking KL grading for both hands or one hand and individuals with missing age or gender information were excluded from all analyses in all cohorts. As the KL sum score has a skewed distribution, the top finding of the meta-analysis was repeated in the discovery cohorts using a Poisson regression.
Visualisation of the associated loci and the regulatory landscape
For the top GWAS-associated single nucleotide polymorphism (SNP), the linkage disequilibrium (LD) region (r2 >0.8) was determined using the 1000G Phase-1 population using the HaploReg V3 tool.21 Using the ROADMAP-generated reference epigenomes, we determined if any of the variants in high LD were located in potential gene regulatory regions in primary osteoblasts (generated by ENCODE) and bone marrow-derived chondrocytes (ROADMAP).22 23 The 18-state chromatin reference epigenomes were downloaded from the ROADMAP epigenomes data portal.23 SNPs and regulatory annotations were visualised using the UCSC Genome Browser on GRCh37/hg19.24 For each variant, it was also determined if the alternative allele would disrupt a protein binding motif; this was done using the HaploReg V3 tool.21
RNA sequencing data
Post-RNA isolation (Qiagen RNeasy Mini Kit, RIN >7) of 40 knee (15 paired preserved (P) and OA lesioned (OAL), 7 P only and 3 OAL only) and 28 hip (six paired P and OAL, 14 P only and 2 OAL only) cartilage samples (online supplementary table S2), paired-end 2×100 bp RNA library sequencing (Illumina TruSeq RNA-Library Prep Kit, Illumina HiSeq2000) resulted in an average of 10 million fragments per sample. Reads were aligned using GSNAP against GRCh37/hg19, in which SNPs from the Genome of the Netherlands consortium with a minor allele frequency (MAF) >1% were masked to prevent alignment bias. Number of fragments per gene were used to assess quantile-adjusted conditional maximum likelihood (edgeR, R-package). Subsequently, differential gene expression analysis was performed pairwise between P and OAL samples for which we had RNA of both (n=21). Allele-specific expression (ASE) was assessed using SNVMix225 with default settings(min coverage=25, 10 reads per allele). The extent of ASE was defined as the fraction of risk allele among all counts at the respective location. Meta-analysis was done only across P samples or OAL when no P counterpart sample was present. p Values were calculated using canonical binominal test (metagen R-package).
Supplementary file 2
TaqMan assay
Conventional TaqMan genotyping was performed on both genomic DNA (gDNA), articular cartilage and subchondral bone cDNA. An allele-specific custom TaqMan assay for rs1800801 (Thermo Fisher Scientific) was used to quantify the allele ratio in cDNA samples and were normalised against the gDNA ratio, which was used as an 1:1 allele ratio reference. Each sample has been measured in four (cartilage) or eight (subchondral bone) times, while calculations and statistics were performed as described previously.19 26 Cartilage samples that yielded fewer than four measurements (n=2) were discarded prior to further analyses. All subchondral bone samples were assessed by eight technical replicates.
Results
GWAS of KL sum score
We conducted a GWAS of a semiquantitative measure of hand OA, a bilateral summed score of KL scores,20 which grades radiographic OA severity, across all hand joints (KL sum score, range of 0–120). The discovery set consisted of three Rotterdam Study cohorts (RSI, RSII and RSIII) and included 8743 participants with KL sum scores. Replication was done in another 4011 individuals from three different cohorts; LS, Framingham Heart Study and Twins UK (TUK). General characteristics of the discovery cohorts and replication cohorts can be found in online supplementary table S3 and text S1.
The discovery analysis yielded two novel independent genome-wide significant loci (p≤5*10−8) on chromosome 12, an intergenic region between matrix Gla protein (MGP) and ERP27 and an intronic region in CCDC91. We also identified seven other novel loci with suggestive significance (p<5*10−6) (figure 1). In total, nine loci were selected for replication in 4011 individuals from three different cohorts (LS, FHS and TUK). Using a Bonferroni corrected p value <5.56*10−3, we significantly replicated one of nine loci, rs4764133 (Betameta=0.83, SEmeta=0.10, p valuereplication=3.4*10−07, p-valuemeta=1.8*10−15) with the same direction of effect as identified in the discovery analysis (table 1 and online supplementary figure S1 Forrest plot). This locus maintained genome-wide significance and another locus near ENPP3 reached near genome-wide significance (chr6: 132063842:D, Betameta=0.58, SEmeta=0.11, p valuemeta=3.8*10−7) in the combined discovery and replication joint meta-analysis (table 1). Since the KL sum score has a skewed distribution, the top hit was also reanalysed in the discovery set using a Poisson regression (rs4764133, Betapoisson=0.12, SEpoisson=0.02, p valuepoisson=1.98*10−11).
GWAS results for association with the KL sum score in the discovery phase. Manhattan plot for association with the KL sum score, adjusted for age and sex, in the discovery cohorts of RSI, RSII and RSIII. The −log10 p values, for each of the ~11 million SNPs analysed (remaining after EASYQC quality control) as part of the genome-wide association with the KL sum score, plotted against their position per chromosome. The red dotted horizontal line corresponds to the genome-wide significant threshold (p=5×10−8). The dotted grey line corresponds to the selection for replication threshold (p=5×10−6). SNP location represented by [], if the SNP is localised intergenic the dashes denotes the distance, -≤10 kb, --≤100 kb, ---≤1000 kb, ----≤1 Mkb, -----≥1 Mkb. GWAS, genome-wide association studies; KL, Kellgren and Lawrence score.
Results GWAS quantitative bilateral phenotype of osteoarthritis of the hand (KL sum score), discovery, replication and meta-analysis
Our top replicated and genome-wide significant finding, rs4764133 [T] (Pmeta=1.80*10−15, Beta=0.83, MAF=0.39) is located in a non-coding intergenic region between MGP and endoplasmic reticulum protein 27 (ERP27). However, variants in high LD with rs4764133 (r2≥0.8) span a ~80 Kb region encompassing multiple genes, including MGP and an open-reading frame C12orf60 (figure 2A). Moreover, several of these variants are located in an mRNA transcript, including a non-synonymous variant in MGP, and variants in 3′ and 5′UTR of MGP and C12orf60 (table 2, figure 2B). The non-synonymous variant in MGP, rs4236, is predicted to be non-damaging (STIFT=1, tolerated; polyPhen=0, benign) causing a threonine to alanine amino acid substitution. Two variants are located in predicted active promoter region of MGP (rs1800801) and C12orf60 (rs9668569) in chondrogenic cells and primary osteoblasts (table 2).
Locus zoom plot for rs4764133, 150 kb upstream and downstream of rs4764133 has been taken as plotted region (A). Zoom in on MGP and three SNPs in high LD with top SNP that are located in the MGP mRNA transcript (B). Also represented is ROADMAP chromatin 18-state data of two tissue types: human mesenchymal stem cell (hMSC)-derived cultured chondrogenic cells and primary osteoblasts. In both these cell types, the chromatin contains active marks surrounding the MGP promoter. LD, linkage disequilibrium; MGP, matrix Gla protein.
rs4764133 LD block (r²>0.8) annotation of potential functional elements in osteoblasts and chondrogenic cells, X marks no potential functional annotation, that is, enhancer region, promoter region or altered protein binding motifs
Next, we investigated the association of rs4764133 with bilateral severe hand OA and bilateral finger OA using the discovery set (RSI, RSII and RSIII). We found a strong association with finger OA (p value=3.09*10–08, OR=1.25) and nominal significant association with severe hand OA (p value=2.80*10–2, OR=1.36), which has a low frequency in the population (online supplementary table S4). To see if rs4764133 also confers risk for other forms of OA, that is, OA of the hip and knee, we used the GWAS summary data of the treat OA consortium27 and the recently published minimal joint space width of the hip (mJSW) meta-analysis.18 No association was found between rs4764133 and hip or knee OA (online supplementary table S4). However, we did find a nominal significant association between rs1049897 and cartilage thickness in the hip joint(mJSW) (r2=0.98 with rs4764133) (p value=1.28*10–2, Beta=−0.398).
Gene expression analyses
In order to identify potential causal genes located in the LD block surrounding rs4764133, we assessed gene expression of MGP, ERP27, ART4, SMOC3 (C12orf69) and C12orf60 in articular cartilage, the primary OA affected tissue. RNA sequencing was obtained on articular cartilage from patients with primary OA who had total joint replacement surgeries of either the knee (n=25) or hip (n=22) joint. Expression levels of ERP27, C12orf60, ART4 and SMOC3 were substantially lower than MGP expression levels in articular cartilage (online supplementary figure S2A). Nonetheless, neither MGP, ERP27, ART4,SMOC3 nor C12orf60 showed significant difference in gene expression between paired P and OAL articular cartilage. However, while these genes are not differentially expressed in OA affected cartilage, it is possible that the identified GWAS SNPs affect gene transcription. When we analysed the relationship between the top SNP and expression analysis in a classical expression quantitative trait loci (eQTL) analysis, we did not to detect significant correlations between rs1049897, rs4236 or rs1800801 and absolute MGP, ERP27, ART4, SMOC3 or C12orf60 expression levels (online supplementary figure S2B). However, we did observe several variants in high LD located in the mRNA transcript of MGP and C12orf60, allowing us to assess allele specific expression (ASE) for these genes. We were unable to study ASE for ART4, SMOC3 and ERP27, since no SNP in high LD with rs4764133 is present in the coding region. In ASE, the influence of exonic alleles on gene expression in -cis is measured within heterozygote subjects, circumventing strong effects from environmental or trans-acting influences. This property results in ASE analysis to be a more statistically powerful approach, when compared with classical eQTL analysis.28 Subsequently, we found that the OA risk alleles for three coding variants in high LD with the lead variant, rs4236 (figure 3SA, 39.6% C allele, p<5*10–16), rs1049897 (online supplementary figure S3B, 44.4% A allele, p<5*10–10) and rs1800801 (figure 3A, 40.7% T allele, p<5*10–16), were significantly correlated with lower expression of MGP, marking imbalanced expression among heterozygotes, independent of the disease status of the articular cartilage. No ASE was observed between SNPs rs11276, rs3088189 and rs1861698 (residing in C12orf60 and in high LD with the lead SNP, r2>0.8, table 2). Technical and biological replication was performed using a custom allele-specific TaqMan assay for rs1800801 in eight additional heterozygous individuals for which we isolated RNA from P cartilage (n=2), OAL (n=2) or both (n=4) from patients with primary knee OA and confirmed the observed imbalance in preserved articular cartilage (figure 3B, relative allelic difference=0.92, p<1*10−6), as well as in eight knee subchondral bone samples (figure 3C, relative allelic difference=0.78, p<1*10−4).
Supplementary file 1
Allelic imbalanced expression of MGP marked by the alleles among heterozygotes of rs1800801 (A), in the assessed cartilage RNA sequencing dataset. Validation of selected rs1800801 using a custom TaqMan assay confirmed the imbalance (B). Allelic imbalance was also assessed with a custom TaqMan assay in subchondral bone samples (C). Preserved (P) and OA lesioned (OAL) samples are shown respectively in blue and red, and genomic DNA (TaqMan control) in black (G). For ASE results for rs4236 and rs1049879, see online supplementary figure S3 and for information on the samples, see online supplementary table S2. ASE, allele-specific expression; MGP, matrix Gla protein; OA, osteoarthritis.
Discussion
Here, we show for the first time, that there is a robust genome-wide significant association between rs4764133, located near MGP, and hand OA. Furthermore, we performed functional validation showing that MGP coding variants in LD with rs4764133 are associated with ASE of MGP, which may increase risk of hand OA by lowering inhibition of articular cartilage calcification, since MGP is an essential inhibitor of cartilage calcification.29 30 These findings suggest that MGP could be considered a prioritised drug target for hand OA, since genetically supported drug targets double the success rate of therapeutics in clinical development.31
MGP is an essential inhibitor of cartilage calcification, and genetic deficiencies of MGP in humans and mice have been linked to abnormal mineralisation of soft tissues, including cartilaginous tissue.29 32 Furthermore, MGP has been previously implicated in relation to OA. A small candidate study reported marginally significant association between hand OA and genetic variants in MGP (rs1800802 and rs4236).33 This is consistent with our findings that the minor allele for rs4764133 and related coding variants in high LD (r2>0.8), rs1800802 and rs4236, increase the risk of hand OA and that we found high expression of MGP in both P and OAL articular cartilage. In contrast, another study showed that an MGP protein complex is excreted by healthy articular chondrocytes, but not by OA-affected chondrocytes,34 although we only assessed MGP expression and not MGP protein complex excretion.
Although the loci with ASE are known to be enriched for eQTLs,35 we were unable to detect an association between the MGP genotype and MGP RNA expression levels in cartilage. This could have been due to our modest sample size (knee joint, n=25 and or hip joint, n=22) in combination with large heterogeneity of the tissue. Notably, the available cartilage samples originated from different joint sites (knee and hip) and different disease stage (preserved versus affected) and had large age range of the individuals. Also, it is known that ASE is a more powerful technique than classical eQTL analysis to identify functional SNPs influencing expression of genes.28 While the extent of imbalance could be considered relatively modest, an increasing number of OA associated SNP alleles appear to mark ASE by comparable amount.19 36–38 From a more biological perspective, one could consider a prolonged, although slight, deviation from homeostasis due to modest ASE of cartilage relevant genes to be of substantial influence over time. This latter hypothesis could contain the molecular basis for increased risk towards developing OA among the ageing population. Additionally, we observed that the rs1800801 alleles also affected expression of MGP in subchondral bone samples. This could imply that, in parallel to an effect in cartilage, the presumed disturbed cartilage homeostasis is further affected by the underlying bone, further enabling the view that OA is a pathology of the entire joint.
Our findings may give an explanation for the known vitamin K association with OA: MGP-mediated calcification inhibition is dependent on γ-carboxylation by vitamin K.39 It has been shown that low vitamin K intake is correlated with OA.40 Thus, vitamin K intake may be a potential therapeutic treatment in OA. Recently, a first randomised control trial testing the effects of vitamin K on OA was published, which reported no overall effect of vitamin K on hand OA.41 Despite the low power of the trial, there was a significant beneficial effect on joint space narrowing (cartilage degradation) among those individuals that were vitamin K deficient at the start of the trial.41 Thus, an adequately powered study of vitamin K may be justified based on the found MGP association. Furthermore, genetic predisposition for hand OA was not taken into account in the trial. Perhaps, genetic predisposition for hand OA (MGP-risk variants) in combination with insufficient vitamin K intake might potentiate cartilage calcification and subsequent risk for developing hand OA. Therefore, future OA trails, therapeutic and preventive treatments might benefit from taking a personalised medicine approach since genetically supported drug targets double the success rate of therapeutics in clinical development.31
Styrkarsdottir et al19 reported on common genetic variants that associate with severe hand OA, among the replication cohorts were the Leiden and Rotterdam cohorts.19 Although we observe suggestive signals at the reported locus (ALDH1A2 gene, 1p31), the respective variants did not meet the genome-wide significance threshold in our analyses (online supplementary table S5). This difference is likely caused by the markedly different phenotypes that were used for either analyses. Where Styrkarsdottir et al studied a dichotomous severe hand OA phenotype, our phenotype was semiquantitatively phenotype.
To conclude, we here present coding variants in MGP that are associated with radiographic hand OA, and the hand OA risk allele marks lower expression of MGP in articular cartilage. Our findings suggest that MGP might play an important role in hand OA pathogenesis through pathways related to articular cartilage calcification and vitamin K. Better understanding of MGP gene and protein regulation and its relation to vitamin K intake and OA may reveal novel therapeutic drug targets for hand OA.
Acknowledgments
The authors are grateful to the study participants, the staff from the Rotterdam Study and the participating general practitioners and pharmacists. We thank Pascal Arp, Mila Jhamai, Marijn Verkerk, Lizbeth Herrera and Marjolein Peters, MSc, and Carolina Medina-Gomez, MSc, for their help in creating the GWAS database, and Karol Estrada, PhD, Yurii Aulchenko, PhD, and Carolina Medina-Gomez, MSc, for the creation and analysis of imputed data. We thank Nico Lakenberg, Ruud van der Breggen and Eka Suchiman for their help in preparing DNA and RNA samples.
References
Footnotes
IM and JJBM contributed equally,
WH and CGB contributed equally.
Handling editor Tore K Kvien
Contributors WdH and CGB contributed equally to this work. DJH, MSY, YFMR and SM performed replication analysis for this work, and LB provided analysis help. LAC and FR provided data. MK provided phenotypic contribution to the GARP study. MP provided data and analyses. TDS contributed data for replication. AH contributed data of the RS cohorts. JD, and PES contributed to genotyping data and analyses of LLS cohort. RGHHN provided contribution to the RAAK study. AGU contributed genotype data of RS cohorts. DTF and AMV contributed replication data for this work. IM and JJBvM jointly supervised this work.
Funding This study was funded by The Netherlands Society for Scientific Research (NWO) VIDI Grant 917103521. The Rotterdam Study is funded by Erasmus Medical Center and Erasmus University, Rotterdam, Netherlands Organization for the Health Research and Development (ZonMw), the Research Institute for Diseases in the Elderly (RIDE), the Ministry of Education, Culture and Science, the Ministry for Health, Welfare and Sports, the European Commission (DG XII) and the Municipality of Rotterdam. The generation and management of GWAS genotype data for the Rotterdam Study (RS I, RS II and RS III) was executed by the Human Genotyping Facility of the Genetic Laboratory of the Department of Internal Medicine, Erasmus MC, Rotterdam, The Netherlands. The GWAS datasets are supported by the Netherlands Organisation of Scientific Research NWO Investments (nr. 175.010.2005.011, 911-03-012), the Genetic Laboratory of the Department of Internal Medicine, Erasmus MC, the Research Institute for Diseases in the Elderly (014-93-015; RIDE2), the Netherlands Genomics Initiative (NGI)/Netherlands Organisation for Scientific Research (NWO) Netherlands Consortium for HealthyAging (NCHA), project nr. 050-060-810. The Leiden University Medical Centre, the Dutch Arthritis Association and Pfizer Inc., Groton, CT, USA, support the GARP study, while the LLS was supported by the Netherlands Organization of Scientific Research (MW 904-61-095, 911-03-016, 917-66-344 and 911-03-012), Leiden University Medical Centre and by the ’Centre of Medical System Biology' and the ’Netherlands Consortium of Healthy Aging' in the framework of the Netherlands Genomics Initiative (NGI). Furthermore, the research leading to the RAAK biobank and the current results has received funding from the Dutch Arthritis Association (DAA 2010_017) and the European Union’s Seventh Framework Programme (FP7/2007-2011) under grant agreement no. 259679. TwinsUK is funded by the Wellcome Trust, Medical Research Council, European Union, the National Institute for Health Research (NIHR)-funded BioResource, Clinical Research Facility and Biomedical Research Centre based at Guy’s and St Thomas’ NHS Foundation Trust in partnership with King’s College London. The Framingham Heart Study of the National Heart, Lung, and Blood Institute of the National Institutes of Health and Boston University School of Medicine was supported by the National Institutes of Health (contract no. HHSN268201500001I, N01-HC-25195, AG18393, AR47785) and its contract with Affymetrix, Inc. for genotyping services (N02‐HL‐6‐4278). Analyses reflect intellectual input and resource development from the Framingham Heart Study investigators participating in the SNP Health Association Resource (SHARe) project. MSY is supported by the National Institutes of Aging (T32AG023480).
Competing interests None declared.
Patient consent Detail has been removed from this case description/these case descriptions to ensure anonymity. The editors and reviewers have seen the detailed information available and are satisfied that the information backs up the case the authors are making.
Ethics approval Ethics committees of the participating studies.
Provenance and peer review Not commissioned; externally peer reviewed.
Author note IM and JJBvM: these authors jointly supervised this work.
Correction notice This article has been corrected since it published Online First. The equal contribution statement has been added.