Advertisement

Characterization of dynamic changes in Matrix Gla Protein (MGP) gene expression as function of genetic risk alleles, osteoarthritis relevant stimuli, and the vitamin K inhibitor warfarin

Open AccessPublished:May 09, 2021DOI:https://doi.org/10.1016/j.joca.2021.05.001

      Summary

      Objective

      We here aimed to characterize changes of Matrix Gla Protein (MGP) expression in relation to its recently identified OA risk allele rs1800801-T in OA cartilage, subchondral bone and human ex vivo osteochondral explants subjected to OA related stimuli. Given that MGP function depends on vitamin K bioavailability, we studied the effect of frequently prescribed vitamin K antagonist warfarin.

      Methods

      Differential (allelic) mRNA expression of MGP was analyzed using RNA-sequencing data of human OA cartilage and subchondral bone. Human osteochondral explants were used to study exposures to interleukin one beta (IL-1β; inflammation), triiodothyronine (T3; Hypertrophy), warfarin, or 65% mechanical stress (65%MS) as function of rs1800801 genotypes.

      Results

      We confirmed that the MGP risk allele rs1800801-T was associated with lower expression and that MGP was significantly upregulated in lesioned as compared to preserved OA tissues, mainly in risk allele carriers, in both cartilage and subchondral bone. Moreover, MGP expression was downregulated in response to OA like triggers in cartilage and subchondral bone and this effect might be reduced in carriers of the rs1800801-T risk allele. Finally, warfarin treatment in cartilage increased COL10A1 and reduced SOX9 and MMP3 expression and in subchondral bone reduced COL1A1 and POSTN expression.

      Discussion & conclusions

      Our data highlights that the genetic risk allele lowers MGP expression and upon OA relevant triggers may hamper adequate dynamic changes in MGP expression, mainly in cartilage. The determined direct negative effect of warfarin on human explant cultures functionally underscores the previously found association between vitamin K deficiency and OA.

      Keywords

      Introduction

      Osteoarthritis (OA) is the most common degenerative disease of joints and its incidence is rising with increasing obesity and age, resulting in a high social and economic burden on society. Interacting risk factors for OA include obesity, age, sex, abnormal loading and genetic factors. The genetic component of OA is estimated to be in the range of 40–60%
      • Spector T.D.
      • Cicuttini F.
      • Baker J.
      • Loughlin J.
      • Hart D.
      Genetic influences on osteoarthritis in women: a twin study.
      ,
      • Felson D.T.
      Risk factors for osteoarthritis: understanding joint vulnerability.
      . For that matter, large-scale genome wide association studies (GWAS) have identified strong, in other words highly significant and reproducible, OA risk genes involved in the aetiology of OA, whereas follow-up studies have shown that risk single nucleotide polymorphisms (SNPs) frequently modulate pathology due to altering transcription of the genes in cis both in bone and cartilage
      • Bos S.D.
      • Bovee J.V.
      • Duijnisveld B.J.
      • Raine E.V.
      • van Dalen W.J.
      • Ramos Y.F.
      • et al.
      Increased type II deiodinase protein in OA-affected cartilage and allelic imbalance of OA risk polymorphism rs225014 at DIO2 in human OA joint tissues.
      • Reynard L.N.
      • Bui C.
      • Canty-Laird E.G.
      • Young D.A.
      • Loughlin J.
      Expression of the osteoarthritis-associated gene GDF5 is modulated epigenetically by DNA methylation.
      • Loughlin J.
      Genetic contribution to osteoarthritis development: current state of evidence.
      • Shepherd C.
      • Zhu D.
      • Skelton A.J.
      • Combe J.
      • Threadgold H.
      • Zhu L.
      • et al.
      Functional characterisation of the osteoarthritis genetic risk residing at ALDH1A2 identifies rs12915901 as a key target variant.
      .
      In this regard Matrix Gla protein (MGP) via rs4764133
      • den Hollander W.
      • Boer C.G.
      • Hart D.J.
      • Yau M.S.
      • Ramos Y.F.M.
      • Metrustry S.
      • et al.
      Genome-wide association and functional studies identify a role for matrix Gla protein in osteoarthritis of the hand.
      with proxy SNPs rs1800801 and rs4236
      • Misra D.
      • Booth S.L.
      • Crosier M.D.
      • Ordovas J.M.
      • Felson D.T.
      • Neogi T.
      Matrix Gla protein polymorphism, but not concentrations, is associated with radiographic hand osteoarthritis.
      , was previously identified as strong OA risk gene for hand OA with the OA conferring allele associated with lower expression of MGP relative to the non-risk allele
      • den Hollander W.
      • Boer C.G.
      • Hart D.J.
      • Yau M.S.
      • Ramos Y.F.M.
      • Metrustry S.
      • et al.
      Genome-wide association and functional studies identify a role for matrix Gla protein in osteoarthritis of the hand.
      in a range of joint tissues but its effect was most profound in cartilage and subchondral bone
      • den Hollander W.
      • Boer C.G.
      • Hart D.J.
      • Yau M.S.
      • Ramos Y.F.M.
      • Metrustry S.
      • et al.
      Genome-wide association and functional studies identify a role for matrix Gla protein in osteoarthritis of the hand.
      ,
      • Shepherd C.
      • Reese A.E.
      • Reynard L.N.
      • Loughlin J.
      Expression analysis of the osteoarthritis genetic susceptibility mapping to the matrix Gla protein gene MGP.
      . On the other hand, these studies could not identify significant differential expression of MGP in OA pathophysiology in macroscopically lesioned OA compared to preserved cartilage
      • den Hollander W.
      • Boer C.G.
      • Hart D.J.
      • Yau M.S.
      • Ramos Y.F.M.
      • Metrustry S.
      • et al.
      Genome-wide association and functional studies identify a role for matrix Gla protein in osteoarthritis of the hand.
      nor in macroscopically preserved cartilage compared to healthy cartilage
      • Shepherd C.
      • Reese A.E.
      • Reynard L.N.
      • Loughlin J.
      Expression analysis of the osteoarthritis genetic susceptibility mapping to the matrix Gla protein gene MGP.
      . These differential expression analyses were, however, determined in a relatively small sample size.
      MGP regulates extracellular calcium levels via high affinity to its γ-carboxyglutamic acid (Gla) residues. Low MGP levels results in higher calcification of cartilage tissue and a reduced bone mineral density
      • Luo G.
      • Ducy P.
      • McKee M.D.
      • Pinero G.J.
      • Loyer E.
      • Behringer R.R.
      • et al.
      Spontaneous calcification of arteries and cartilage in mice lacking matrix GLA protein.
      • Yagami K.
      • Suh J.Y.
      • Enomoto-Iwamoto M.
      • Koyama E.
      • Abrams W.R.
      • Shapiro I.M.
      • et al.
      Matrix GLA protein is a developmental regulator of chondrocyte mineralization and, when constitutively expressed, blocks endochondral and intramembranous ossification in the limb.
      • Tunon-Le Poultel D.
      • Cannata-Andia J.B.
      • Roman-Garcia P.
      • Diaz-Lopez J.B.
      • Coto E.
      • Gomez C.
      • et al.
      Association of matrix Gla protein gene functional polymorphisms with loss of bone mineral density and progression of aortic calcification.
      . As the OA risk allele (rs1800801) has been associated with a reduced MGP gene expression
      • Tunon-Le Poultel D.
      • Cannata-Andia J.B.
      • Roman-Garcia P.
      • Diaz-Lopez J.B.
      • Coto E.
      • Gomez C.
      • et al.
      Association of matrix Gla protein gene functional polymorphisms with loss of bone mineral density and progression of aortic calcification.
      and with increased vascular calcification
      • Sheng K.
      • Zhang P.
      • Lin W.
      • Cheng J.
      • Li J.
      • Chen J.
      Association of Matrix Gla protein gene (rs1800801, rs1800802, rs4236) polymorphism with vascular calcification and atherosclerotic disease: a meta-analysis.
      , this would suggest increased cartilage calcification in carriers of the OA risk allele. The latter was further justified by recapitulating downregulation of MGP in cartilage chondrocytes resulting in pro-catabolic (ADAMTS4, MMP13), as well as pro-hypertrophic (COL10A1, VEGFA) mRNA signalling
      • Shepherd C.
      • Reese A.E.
      • Reynard L.N.
      • Loughlin J.
      Expression analysis of the osteoarthritis genetic susceptibility mapping to the matrix Gla protein gene MGP.
      . The MGP protein is produced by the cell in inactive form and is dependent on vitamin K for activation, via carboxylation (c-MGP). As such, low vitamin K levels have been hypothesized to play a role in OA pathogenesis
      • Neogi T.
      • Booth S.L.
      • Zhang Y.Q.
      • Jacques P.F.
      • Terkeltaub R.
      • Aliabadi P.
      • et al.
      Low vitamin K status is associated with osteoarthritis in the hand and knee.
      ,
      • Shea M.K.
      • Kritchevsky S.B.
      • Hsu F.C.
      • Nevitt M.
      • Booth S.L.
      • Kwoh C.K.
      • et al.
      The association between vitamin K status and knee osteoarthritis features in older adults: the Health, Aging and Body Composition Study.
      . Similarly, vitamin K antagonists such as warfarin, that are frequently prescribed for the prevention of thromboembolic events in patients with atrial fibrillation
      • Holbrook A.
      • Schulman S.
      • Witt D.M.
      • Vandvik P.O.
      • Fish J.
      • Kovacs M.J.
      • et al.
      Evidence-based management of anticoagulant therapy: antithrombotic therapy and prevention of thrombosis, 9th ed: American college of chest physicians evidence-based clinical practice guidelines.
      , have been suggested to predispose to OA
      • Chin K.Y.
      The relationship between vitamin K and osteoarthritis: a review of current evidence.
      . Nonetheless, the direct effect of warfarin on human articular cartilage tissue homeostasis has not been assessed.
      Here we set out to explore MGP gene expression in relation to the OA risk allele rs1800801-T, in a large RNA-sequencing dataset containing both macroscopically preserved and lesioned cartilage
      • Coutinho de Almeida R.
      • Ramos Y.F.M.
      • Mahfouz A.
      • den Hollander W.
      • Lakenberg N.
      • Houtman E.
      • et al.
      RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage.
      and subchondral bone
      • Tuerlings M.
      • van Hoolwerff M.
      • Houtman E.
      • Suchiman E.
      • Lakenberg N.
      • Mei H.
      • et al.
      RNA sequencing reveals interacting key determinants of osteoarthritis acting in subchondral bone and articular cartilage.
      as well as in our recently established full thickness human ex vivo osteochondral explant model
      • Houtman E.
      • van Hoolwerff M.
      • Lakenberg N.
      • Suchiman E.H.D.
      • van der Linden-van der Zwaag E.
      • Nelissen R.
      • et al.
      Human osteochondral explants: reliable biomimetic models to investigate disease mechanisms and develop personalized treatments for osteoarthritis.
      . The latter allowing us to study the effect of the OA risk allele on the dynamic MGP response to different OA related stimuli, such as inflammation (Interleukin one beta (IL-1β)), hypertrophy (Triiodothyronine (T3)) and 65% mechanical stress (65%MS). Moreover, we used the human ex vivo explant model to study the direct effect of vitamin K antagonist, warfarin, on articular cartilage and subchondral bone homeostasis.

      Material and methods

      Sample description

      Human material was obtained from the Research in Articular Osteoarthritis Cartilage (RAAK) biobank as previously described in detail
      • Ramos Y.F.
      • den Hollander W.
      • Bovee J.V.
      • Bomer N.
      • van der Breggen R.
      • Lakenberg N.
      • et al.
      Genes involved in the osteoarthritis process identified through genome wide expression analysis in articular cartilage; the RAAK study.
      . The RAAK study is approved by the medical ethics committee of the Leiden University Medical Center (P08.239/P19.013). In this study, RNA-sequencing data was included of paired macroscopically preserved and lesioned OA cartilage of N = 35 participants
      • Coutinho de Almeida R.
      • Ramos Y.F.M.
      • Mahfouz A.
      • den Hollander W.
      • Lakenberg N.
      • Houtman E.
      • et al.
      RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage.
      and subchondral bone of N = 24 participants
      • Tuerlings M.
      • van Hoolwerff M.
      • Houtman E.
      • Suchiman E.
      • Lakenberg N.
      • Mei H.
      • et al.
      RNA sequencing reveals interacting key determinants of osteoarthritis acting in subchondral bone and articular cartilage.
      for which sample characteristics have previously been described. In total 136 osteochondral explants were harvested from the macroscopically preserved condyle knee joints of N = 18 participants and divided over the different experiment (Supplementary Fig. 1). Multiple osteochondral explants containing both cartilage and bone (diameter of 8 mm) were extracted per participant and washed in sterile PBS before taking into culture. Donor characteristics of osteochondral explants are described in Supplementary Table 1 and study design is described in Supplementary Fig. 1. For additional details on neo-cartilage deposition, RNA and DNA isolations, TaqMan genotyping, RNA sequencing data of cartilage and subchondral bone, unconfined dynamic (cyclic) compression, expression quantitative trait loci (eQTL), Allelic Expression Imbalance (AEI) and data analysis, see the Supplementary Methods.

      Treatment of osteochondral explants

      Explants were cultured in 24 wells plates (Greiner CELLSTAR; Sigma) supplemented with 1.5 ml CDM in a 5% (v/v) CO2 incubator at 37°C. Three days after extraction, explants were treated with either IL-1β (10 ng/ml), triiodothyronine (T3, 10 nM; Sigma) or warfarin (50 μM; Sigma), depicted in Fig. 1(A). Six days after extraction, dynamic unconfined compression was applied to explant tissues using the Mach-1 mechanical testing system (Biomomentum Inc., Laval, QC, Canada) on four subsequent days [Fig. 1(B)]. Mechanical stress was applied at a strain of 65% of cartilage height and at a frequency of 1 Hz, mimicking walking speed [Fig. 1(C)]. Cartilage and bone were separated using a scalpel, snap-frozen in liquid nitrogen, and stored at −80°C for RNA isolation.
      Fig. 1
      Fig. 1Schematic representation of different perturbations applied to osteochondral explants. Osteochondral explants were punched from the still macroscopically preserved looking knee condyle area and taken into culture. [A] Explants were subjected to treatment with IL-1β (10 ng/ml), triiodothyronine (T3; 10 nM) or warfarin (50 μM). [B] Explants received mechanical stresses at a strain of 65% for 10 min per day on four subsequent days. On day 13, cartilage and bone was separated, snap frozen and stored at −80C. [C] Schematic representation of the dynamic (cyclic) compression applied to osteochondral explants. Legend: N=Newton; h=height.

      Reverse transcription and Real-Time PCR

      Real-Time PCR for gene expression was performed with QuantStudio 6 Real-Time PCR system (Applied Biosystems) using Fast Start Sybr Green Master mix (Roche Applied Science). Primer sequences (Table I) used were tested for linear amplification and missing datapoints for genes are summarized in Supplementary Tables 2 and 3. Details on normalization can be found in the Supplementary Methods.
      Table IPrimer sequence used to determine gene expression levels in real-time PCR
      Gene nameForward 5′-3′Reverse 5′-3′
      SDHATGGAGCTGCAGAACCTGATGTGTAGTCTTCCCTGGCATGC
      MGPCGCCCCCAGATTGATAAGTATCTCCTTTGACCCTCACTGC
      SOX9CCCCAACAGATCGCCTACAGCTGGAGTTCTGGTGGTCGGT
      ACANAGAGACTCACACAGTCGAAACAGCCTATGTTACAGTGCTCGCCAGTG
      COL2A1CTACCCCAATCCAGCAAACGTAGGTGATGTTCTGGGAGCCTT
      RUNX2CTGTGGTTACTGTCATGGCGAGGTAGCTACTTGGGGAGGA
      ALPLCAAAGGCTTCTTCTTGCTGGTGCCTGCTTGGCTTTTCCTTCA
      COL1A1GTGCTAAAGGTGCCAATGGTACCAGGTTCACCGCTGTTAC
      COL10A1GGCAACAGCATTATGACCCATGAGATCGATGATGGCACTCC
      MMP3GAGGCATCCACACCCTAGGTTTCAGAAATGGCTGCATCGATT
      MMP13TTGAGCTGGACTCATTGTCGGGAGCCTCTCAGTCATGGAG
      ADAMTS5TGGCTCACGAAATCGGACATGCGCTTATCTTCTGTGGAACC
      COMPACAATGACGGAGTCCCTGACTCTGCATCAAAGTCGTCCTG
      OMDGGACACAACAAATTGAAGCAAGCTGGTGGTAATGTAGTGGGTCA
      BGLAPCCCTCCTGCYYGGACACAAACACACTCCTCGCCCTATTGG
      OGNTGATGAAATGCCCACGTGTCTTTGGTAAGGGTGGTACAGCA
      SPP1GCCAGTTGCAGCCTTCTCAAAAAGCAAATCACTGCAATTCTCA
      TNFRSF11BTTGATGGAAAGCTTACCGGGATCTGGTCACTGGGTTTGCATG
      BMP2TCCATGTGGACGCTCTTTCAAGCAGCAACGCTAGAAGACA
      POSTNTACACTTTGCTGGCACCTGTTTTAAGGAGGCGCTGATCCA

      Statistical analyses

      Differential MGP expression analyses between preserved and lesioned OA cartilage and bone including false discovery rates (FDR) as multiple testing correction for the genome wide analyses were reproduced from Coutinho
      • Coutinho de Almeida R.
      • Ramos Y.F.M.
      • Mahfouz A.
      • den Hollander W.
      • Lakenberg N.
      • Houtman E.
      • et al.
      RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage.
      and Tuerlings
      • Tuerlings M.
      • van Hoolwerff M.
      • Houtman E.
      • Suchiman E.
      • Lakenberg N.
      • Mei H.
      • et al.
      RNA sequencing reveals interacting key determinants of osteoarthritis acting in subchondral bone and articular cartilage.
      , respectively. Description on their study design and sample numbers are in Supplementary Methods. To assess allelic expression imbalance (AEI) we applied our previously published methodology in R
      • den Hollander W.
      • Pulyakhina I.
      • Boer C.
      • Bomer N.
      • van der Breggen R.
      • Arindrarto W.
      • et al.
      Annotating transcriptional effects of genetic variants in disease-relevant tissue: transcriptome-wide allelic imbalance in osteoarthritic cartilage.
      to RNA sequence data of MGP in the larger dataset of cartilage
      • Coutinho de Almeida R.
      • Ramos Y.F.M.
      • Mahfouz A.
      • den Hollander W.
      • Lakenberg N.
      • Houtman E.
      • et al.
      RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage.
      and a dataset of bone
      • Tuerlings M.
      • van Hoolwerff M.
      • Houtman E.
      • Suchiman E.
      • Lakenberg N.
      • Mei H.
      • et al.
      RNA sequencing reveals interacting key determinants of osteoarthritis acting in subchondral bone and articular cartilage.
      which is further outlined in Supplementary Methods.
      To test expression quantitative trait analyses (eQTL) and differential expression of MGP in genotype strata in the current manuscript, we used the variance stabilizing transformation (VST) normalized MGP expression levels of these RNA sequencing datasets and used generalized estimating equations (GEE)
      • Zeger S.L.
      • Liang K.Y.
      Longitudinal data analysis for discrete and continuous outcomes.
      to effectively adjust for dependencies of genotypes among donors by adding a random effect for sample donor. Details of the models applied are outlined in Supplementary Methods.
      MGP expression by RT-qPCR in the in vitro 3D-neo cartilage formation was estimated using a generalized linear mixed model (GLMM) using MGP levels (-ΔCT) as dependent variable and time as repeated measure: MGP level ∼ Time + (1|Donor). In the osteochondral explant models fold changes (FC) of RT-qPCR expression were determined by calculating the log 2 of the -ΔΔCt for each sample (2−ΔΔCT) where FC > 1 is upregulation and FC < 1 is downregulation of treated samples compared to control samples. The reported P-values were determined by applying GEE to -ΔCT values to effectively adjust for dependencies among donors of the explants by adding a random effect for sample donor as we did not have perfect pairs for each analysis. We followed a linear GEE model, with MGP level as dependent variable, treatment as factor and exchangeable working matrix: MGP level ∼ Treatment + (1|Donor)
      • Diggle P.
      • Liang K.-Y.
      • Zeger S.L.
      . Differences in effect sizes between strata was determined by performing unpaired student's t-test on the fold changes corrected for control samples. Warfarin treated osteochondral explants samples were paired hence a paired sample t-test was performed to determine between-group differences and p-values. Except for AEI, Statistical analyses were performed in SPSS statistics 23 (IBM). Outliers were investigated using Grubbs's test and normal distribution was determined using Shapiro–Wilk test and visually inspecting Q–Q plots. The boxplots represent 25th, 50th and 75th percentiles, and whiskers extend to the 95%CI.

      Results

      Expression patterns of MGP in previous established RNA sequencing datasets of preserved and lesioned OA cartilage and subchondral bone

      We used our previously established RNA sequencing dataset of macroscopically preserved and lesioned OA cartilage samples (N = 35 pairs
      • Coutinho de Almeida R.
      • Ramos Y.F.M.
      • Mahfouz A.
      • den Hollander W.
      • Lakenberg N.
      • Houtman E.
      • et al.
      RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage.
      ) and subchondral bone (N = 24 pairs
      • Tuerlings M.
      • van Hoolwerff M.
      • Houtman E.
      • Suchiman E.
      • Lakenberg N.
      • Mei H.
      • et al.
      RNA sequencing reveals interacting key determinants of osteoarthritis acting in subchondral bone and articular cartilage.
      ), to examine differential MGP expression with OA tissue status and with the OA risk SNPs (see Supplementary Methods and Figures). An increased expression of MGP in lesioned compared to preserved OA cartilage was observed (FC = 1.45, 95%CI [1.24; 1.61], P-value = 1.78 × 10−3) and this increase of MGP was genome wide significant (FDR = 0.021). Similarly, MGP was upregulated in lesioned compared to preserved OA subchondral bone (FC = 1.53, 95%CI [1.22; 1.64], P-value = 0.023), but this was not genome wide significant (FDR = 0.12). Together, these results show a robust upregulation of MGP expression with ongoing OA pathophysiology.
      Here we studied whether the MGP differential expression between preserved and lesioned OA tissues was affected by MGP OA risk allele carriership. As shown in Fig. 2(A), the MGP upregulation occurs particularly among risk allele carriers rs1800801-T in lesioned compared to preserved OA cartilage independent of age and sex of donors (OR = 2.70, 95%CI [1.16; 6.29], P-value = 0.021). Notably however, the overall MGP expression remains lower among risk allele carriers rs1800801-T as compared to carriers of the reference allele rs1800801-C. The same effect was observed in subchondral bone, where MGP was found to be upregulated in lesioned compared to preserved tissue only in risk allele carriers rs1800801-T (OR = 3.04, 95%CI [1.24; 7.45], P-value = 0.015) independent of age and sex of donors [Fig. 2(B)].
      Fig. 2
      Fig. 2MGP expression as function of the transcript and OA risk SNP rs1800801. [A] Variance stabilizing transformation (VST) normalized MGP expression levels extracted from the RNA sequencing dataset in preserved and lesioned OA cartilage stratified for rs1800801 genotype CC (nPreserved = 19 vs nLesioned = 16) and CT + TT (nPreserved = 36 vs nLesioned = 27). [B] Variance stabilizing transformation (VST) normalized MGP expression levels extracted from the RNA sequencing dataset in preserved and lesioned OA subchondral bone stratified for rs1800801 genotype CC (nPreserved = 5 vs nLesioned = 5) and CT + TT (nPreserved = 9 vs nLesioned = 9). The boxplots represent the median with the 25th, 50th and 75th percentiles, and whiskers reflect the total range. Independent samples are depicted by black dots in each graph. To adjust for donor variation, P-values were estimated by performing logistic generalized estimation equations, with tissue status as dependent variable and MGP level, age and sex as covariate: Tissue status ∼ MGP level + age + sex + (1|Donor). ∗P ≤ 0.05.
      Next, we attempted to replicate the previously shown AEI of MGP in association with the OA risk SNP rs1800801
      • den Hollander W.
      • Boer C.G.
      • Hart D.J.
      • Yau M.S.
      • Ramos Y.F.M.
      • Metrustry S.
      • et al.
      Genome-wide association and functional studies identify a role for matrix Gla protein in osteoarthritis of the hand.
      in heterozygous individuals in this larger RNA sequencing dataset of preserved and lesioned OA cartilage
      • Coutinho de Almeida R.
      • Ramos Y.F.M.
      • Mahfouz A.
      • den Hollander W.
      • Lakenberg N.
      • Houtman E.
      • et al.
      RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage.
      and a novel dataset of OA subchondral bone
      • Tuerlings M.
      • van Hoolwerff M.
      • Houtman E.
      • Suchiman E.
      • Lakenberg N.
      • Mei H.
      • et al.
      RNA sequencing reveals interacting key determinants of osteoarthritis acting in subchondral bone and articular cartilage.
      . Additionally, we explored whether the effect size in AEI differed in these tissues between preserved and lesioned areas. As shown in Supplementary Fig. 2(A) we confirmed AEI expression of MGP in preserved OA cartilage with the risk-conferring allele rs1800801-T associated to a reduced MGP expression of 10% (95%CI [2.24; 18.64]) relative to the reference allele rs1800801-C. In lesioned OA cartilage the AEI was very comparable with rs1800801-T associated to a reduced MGP expression of 11% (95%CI [2.25; 19.49]) relative to the reference allele rs1800801-C. In subchondral bone, genotype of rs1800801 could not be called thus we used its proxy SNP rs4236 (r2 = 0.93 with rs1800801), which was also investigated previously
      • den Hollander W.
      • Boer C.G.
      • Hart D.J.
      • Yau M.S.
      • Ramos Y.F.M.
      • Metrustry S.
      • et al.
      Genome-wide association and functional studies identify a role for matrix Gla protein in osteoarthritis of the hand.
      . As shown in Supplementary Fig. 2(B), we confirmed AEI of MGP in preserved OA subchondral bone with the risk-conferring allele rs4236-C associated to a reduced MGP expression of 10% (95%CI [5.70; 14.52]) relative to the reference allele rs4236-T. In lesioned OA subchondral bone AEI was very comparable, with rs4236-C associated to a reduced MGP expression of 12% (95%CI [8.80; 14.55]) relative to the reference allele rs4236-T.
      Finally, we analysed MGP expression levels among genotype carriers of one or two of the OA risk alleles rs1800801-T (eQTL) and confirmed that also overall MGP expression in cartilage is reduced in a dose responsive manner with MGP risk alleles, independent of donor, age, sex, and OA status, i.e., preserved or lesioned (OR = 0.73, 95%CI [0.64; 0.84], P-value = 4.00 × 10−6; Supplementary Fig. 3(A)). In subchondral bone we observed a similar pattern, however this was not significant (Supplementary Fig. 3(B)). Together these data confirm that innate lower MGP expression levels confer risk to OA, though its effect seems more pronounced in articular cartilage.

      MGP expression patterns in human in vitro and ex vivo models and as function of OA related cues

      First, we investigated expression of MGP during neo-cartilage formation using a human in vitro 3D pellet culture with primary chondrocytes. As shown in Fig. 3, MGP is expressed in primary chondrocytes (day-0) and increases during cartilage extracellular matrix (ECM) deposition until day-14, suggesting that MGP expression can be considered a marker of neo-cartilage formation.
      Fig. 3
      Fig. 3Gene expression of MGP in an in vitro 3D model of neo-cartilage. Gene expression of MGP in an in vitro 3D chondrocyte pellet model of neo-cartilage formation (N = 3 donors; day 0: n = 2, all other time points: n = 3). Data is depicted as mean expression (-ΔCT) ± standard error of the mean (SEM) and each dot represents a sample of two combined biological duplicates. Statistical analysis was performed by generalized linear mixed model (GLMM) using MGP levels as dependent variable and time as repeated measure: MGP level ∼ Time + (1|Donor).
      Next, we explored dynamic changes in MGP expression in cartilage and subchondral bone in an established human ex vivo osteochondral explant model
      • Houtman E.
      • van Hoolwerff M.
      • Lakenberg N.
      • Suchiman E.H.D.
      • van der Linden-van der Zwaag E.
      • Nelissen R.
      • et al.
      Human osteochondral explants: reliable biomimetic models to investigate disease mechanisms and develop personalized treatments for osteoarthritis.
      as function of OA related stimuli being inflammation (IL-1β), hypertrophy (T3), and 65% mechanical stress (65%MS). As shown in Fig. 4, we observed in cartilage a consistent and significant downregulation of MGP expression after treatment with IL-1β (FC = 0.03, 95%CI [0.02; 0.06], P-value = 4.40 × 10−7), T3 (FC = 0.80, 95%CI [0.56; 0.97], P-value = 0.046), as well as with mechanical stress (FC = 0.65, 95%CI [0.45; 0.85], P = 0.002). Notable, in Fig. 4, is an outlier in the mechanical stress group, however removing this datapoint did not influence our result (FC = 0.67, 95%CI [0.47; 0.87], P-value = 0.046). In subchondral bone we were not able to isolate RNA for IL-1β treated samples and only observed a significant downregulation of MGP expression after treatment with T3 (FC = 0.81, 95%CI [0.52; 1.10], P-value = 0.015).
      Fig. 4
      Fig. 4Gene expression of MGP in response to three different OA relevant cues in cartilage of osteochondral explants. Gene expression of MGP (-ΔCT) in an ex vivo osteochondral explant model in articular cartilage (A, B, C) and subchondral bone (D, E). MGP expression, represented by the housekeeping gene corrected value (-ΔCT), in articular cartilage upon perturbation with [A] IL-1β (nControl = 6 vs ntreated = 6) [B] T3 (nControl = 21 vs ntreated = 21) and [C] posttraumatic OA after 65% MS (nControl = 30 vs ntreated = 23). MGP expression in subchondral bone upon perturbation with [D] T3 (nControl = 16 vs ntreated = 16) and [E] posttraumatic OA after 65% MS (nControl = 23 vs ntreated = 19). The boxplots represent the median with 25th, 50th and 75th percentiles, and whiskers reflect the total range. -ΔCT of each independent sample is depicted by black dots in the graphs. To adjust for donor variation P-values were determined by performing linear generalized estimation equations, with MGP levels as dependent variable and treatment as factor: MGP level ∼ Treatment + (1|Donor). Far out values are represent by the white filled circle (о) which did not affect the result (see main body text) and therefore analysis including this sample is presented in C. ∗P ≤ 0.05, ∗∗P ≤ 0.01, ∗∗∗P ≤ 0.001. Legend: Ctrl=Control; 65% MS=65% mechanical stress; T3=triiodothyronine.

      Changes in MGP expression in the ex vivo OA models as function of the transcript and OA risk SNP rs1800801

      Since general MGP expression was identified to change between preserved and lesioned OA cartilage and subchondral bone, and in a osteochondral explant model to several OA related stimuli (Fig. 4), we next explored whether the OA risk allele rs1800801-T modified these effects. Hereto we investigated the observed dynamic downregulation of MGP, upon inducing hypertrophy (T3 exposure; Fig. 5(A) and (C)) and mechanical stress (65%MS; Fig. 5(B) and (D)) in our ex vivo cartilage explant model stratified by rs1800801 genotypes. For IL-1β treatment, donor numbers were too low to explore the effect of genotype.
      Fig. 5
      Fig. 5MGP expression as function of the transcript and OA risk SNP rs1800801. MGP expression (-ΔCT) in an ex vivo osteochondral explant model stratified by rs1800801 genotype in articular cartilage (A, B) and subchondral bone (C, D). MGP gene expression in articular cartilage upon perturbation with [A] T3 (CC: nControl = 10 vs ntreated = 8; CT: nControl = 11 vs ntreated = 13) and [B] posttraumatic OA upon 65% mechanical stress (CC: nControl = 11 vs ntreated = 7; CT + TT: nControl = 19 vs ntreated = 16). MGP gene expression in subchondral bone upon perturbation with [C] T3 (CC: nControl = 8 vs ntreated = 9; CT: nControl = 8 vs ntreated = 7) and [D] posttraumatic OA upon 65% mechanical stress (CC: nControl = 10 vs ntreated = 7; CT + TT: nControl = 13 vs ntreated = 12). The boxplots represent the median and 25th, 50th and 75th percentiles, and whiskers reflects the total range. Independent samples are depicted by black dots in each graph. Numeric values associated with this Figure are shown in . To adjust for donor variation P-values were determined by performing linear generalized estimation equations, with MGP levels as dependent variable and treatment as factor: MGP level ∼ Treatment + (1|Donor). Far out values are represent by the white filled circle (о) which did not affect the result (see main body text) and therefore analysis including this sample is presented in B. ∗P ≤ 0.05. Legend: 65% MS=Mechanical stress; T3=triiodothyronine.
      In cartilage (Fig. 5(A) and (B)), we observed that downregulation of MGP occurred particularly among carriers of the reference allele rs1800801-C for hypertrophy (T3) (FC = 0.69, 95%CI [0.49; 0.89]) and for mechanical stress (FC = 0.26, 95%CI [0.14; 0.38]) as compared to carriers of the risk allele rs1800801-T for hypertrophy (FC = 0.92, 95%CI [0.57; 1.27]) and for mechanical stress (FC = 0.85, 95%CI [0.61; 1.09]). Also in the data shown in Fig. 5(B) is the previously identified outlier in the mechanical stress group which did not influence our result upon removal (FC = 0.29, 95%CI [0.17–0.41], P-value = 0.045). As shown in Supplementary Table 4, the difference in response (FC) among carriers of the reference allele rs1800801-C relative to carriers of the OA risk allele rs1800801-T is significant for mechanical stress (FC = 0.34, 95% CI [0.28–0.38], P-value = 2.8 × 10−3). Similarly in bone [Fig. 5(C)–(D)], we observed that downregulation of MGP expression in subchondral bone upon hypertrophy induction (T3) was more pronounced among carriers of the reference allele rs1800801-C (FC = 0.50, 95%CI [0.09; 0.91] as compared to the carriers of the risk allele rs1800801-T (FC = 0.82, 95%CI [0.37; 1.27]). This difference, however, did not reach statistical significance (Supplementary Table 4). For mechanical stress no effects were observed in subchondral bone. Together these data suggest that particularly in cartilage the OA risk allele rs1800801-T may have a different response in MGP expression upon OA relevant cues.

      Treatment of osteochondral explants with warfarin

      Since the activation of MGP is dependent on vitamin K and innate lower MGP expression confers risk to OA, we next investigated the direct effect of the vitamin K antagonist warfarin on articular chondrocyte and subchondral bone signalling. Hereto, ex vivo osteochondral explants (n = 15 pairs for cartilage and n = 13 pairs for subchondral bone) were treated with warfarin. The effect of this reduced vitamin K bioavailability on the cartilage homeostasis was determined by measuring chondroprotective genes (SOX9, COL2A1 and ACAN), genes involved in early/late cartilage hypertrophy (RUNX2, ALPL, COL1A1, COL10A1 and MGP) and catabolic genes (MMP3, MMP13 and ADAMTS5). As shown in Fig. 6(A), warfarin exposure to cartilage reduced expression of SOX9 (FC = 0.87, 95%CI [0.77; 0.97], P-value = 0.023) and MMP3 (FC = 0.56, 95%CI [0.43; 0.69], P-value = 1.02 × 10−5), while increasing COL10A1 (FC = 2.26, 95%CI [1.14; 3.38], P-value = 0.045). In addition, RUNX2 (FC = 1.43, 95%CI [0.96; 1.90], P-value = 0.094), a master transcriptional regulator of chondrocyte maturation, and ALPL (FC = 6.21, 95%CI [1.36; 11.06], P-value = 0.059) show a trend towards upregulation in response to warfarin treatment. In subchondral bone, genes involved in matrix formation (COL10A1, RUNX2, ALPL COL1A1, OMD, BGLAP and OGN) and remodelling (MGP, SPP1, TNFRSF11B, BMP2 and POSTN) were measured. As shown in Fig. 6(B), warfarin exposure to subchondral bone significantly reduced expression of the bone formation marker COL1A1 (FC = 0.81, 95%CI [0.59; 1.03], P-value = 0.046) and the remodelling marker POSTN (FC = 0.67, 95%CI [0.42; 0.92], P-value = 0.011). Together these results show that addition of warfarin to aged osteochondral explants resulted in a significant upregulation of hypertrophic signalling among articular chondrocytes and reduced bone formation and altered remodelling signalling.
      Fig. 6
      Fig. 6Gene expression after 10 days of warfarin treatment. mRNA expression of genes depicted as a fold change (FC) following 50 μM warfarin treatment relative to controls in [A] cartilage (nControl = 15 vs ntreated = 15) and [B] subchondral bone (nControl = 13 vs ntreated = 13). Controls are depicted by the dotted line, while each gray dot represents a warfarin treated sample. The number and percentage of missing data points per gene are summarized in . The light gray bars represent the mean ± standard error of the mean (SEM) of the Fold change (FC) in which FC > 1 represent upregulation and FC < 1 represents downregulation of warfarin treated samples relative to its paired control. The x-axis is given in a log2 scale to depict the up and down regulation in the same scale. Differences in gene levels between warfarin exposure and controls were calculated by means of a paired t-test.∗P ≤ 0.05; ∗∗∗P ≤ 0.001.

      Discussion

      In the current paper we explored (dynamic) changes of MGP expression in relation to the OA risk allele rs1800801-T, in preserved and lesioned OA cartilage, as well as in a human ex vivo explant model subjected to OA related stimuli, such as inflammation, hypertrophy and mechanical stress. Furthermore, we studied the direct effect of the frequently used vitamin K antagonist, warfarin, on articular chondrocyte and subchondral bone signalling. In doing so, we confirm that MGP expression, as inhibitor of calcification via high affinity of calcium to its Gla-residues, should be considered a beneficial marker of articular cartilage. Consequently, the significantly upregulated MGP expression with ongoing OA pathophysiology is likely an attempt of chondrocytes to halt the OA associated osteo-induction. Noteworthy is our observation that the OA risk allele may also hamper adequate dynamic change in expression of MGP in response to OA and relevant cues like mechanical stress (65%MS) and this effect was most pronounced in cartilage. Finally, warfarin treatment to the aged human cartilage explants resulted in a significant upregulation of hypertrophic signalling among articular chondrocytes and reduced bone formation while altering remodelling.
      Similar to previous reports
      • den Hollander W.
      • Boer C.G.
      • Hart D.J.
      • Yau M.S.
      • Ramos Y.F.M.
      • Metrustry S.
      • et al.
      Genome-wide association and functional studies identify a role for matrix Gla protein in osteoarthritis of the hand.
      ,
      • Shepherd C.
      • Reese A.E.
      • Reynard L.N.
      • Loughlin J.
      Expression analysis of the osteoarthritis genetic susceptibility mapping to the matrix Gla protein gene MGP.
      , we here confirmed in a large RNA-sequencing dataset, that the OA risk allele rs1800801-T is associated with lower (overall) expression of MGP in articular cartilage (Supplementary Fig. 3(A)). In addition, we confirmed that MGP gene expression is significantly upregulated in both articular cartilage and subchondral bone in OA pathophysiology. Although, our results showed that this effect was mainly driven by carriers of the rs1800801-T OA risk allele in both tissues, the expression of MGP does not reach the level of that in carriers of the reference allele rs1800801-C [Fig. 2(A) and (B)]. We advocate that MGP upregulation with OA pathophysiology in cartilage is an attempt of chondrocytes to compensate for the osteo-inductive effect of low MGP levels and that this is not sufficient among the MGP OA risk allele carriers. On the other hand, the upregulation of MGP in bone may be a marker of active bone resorption, as it was previously found that MGP inhibits mineralization by osteoblasts while increased MGP expression in osteoclasts mark increased osteoclastic commitment
      • Zhang Y.
      • Zhao L.
      • Wang N.
      • Li J.
      • He F.
      • Li X.
      • et al.
      Unexpected role of matrix Gla protein in osteoclasts: inhibiting osteoclast differentiation and bone resorption.
      . Together our data highlights that, similar to vascular calcification and bone loss
      • Vassalle C.
      • Mazzone A.
      Bone loss and vascular calcification: a bi-directional interplay?.
      , also articular cartilage calcification and bone loss in OA could share a common pathogenetic mechanism involving MGP.
      We also explored the dynamic response of MGP in a human ex vivo explant model while applying OA relevant perturbing cues such as inflammation, hypertrophy, and mechanical stress. The strength of our explant model is that it represents physiological relevant aged human articular cartilage prone to OA pathophysiology, hence suitable to study the initial process of OA related cartilage destruction. Moreover, and despite the inherent heterogeneity between donors, we found in cartilage a consistent downregulation of MGP, associated with matrix mineralization, as general response to OA related perturbations (Fig. 4). Additionally, we showed that the rs1800801-T OA risk allele may hamper such innate dynamic change in MGP expression upon stress. A possible mechanism by which the genetic risk variant modifies response to stress lies in the fact that rs1800801 is localized in the transcription factor binding site (POLR2A, CTCF, p300) of the MGP promoter (Supplementary Fig. 4). In addition, the OA risk allele rs1800801-T was shown to reduce expression between 34% and 47% in a luciferase reporter assay and in silico prediction suggested this to be due to a loss of binding site for the transcription factor c-Ets
      • Tunon-Le Poultel D.
      • Cannata-Andia J.B.
      • Roman-Garcia P.
      • Diaz-Lopez J.B.
      • Coto E.
      • Gomez C.
      • et al.
      Association of matrix Gla protein gene functional polymorphisms with loss of bone mineral density and progression of aortic calcification.
      . In the subchondral bone compartment of the human ex vivo explants, the MGP response to the OA like perturbing cues were smaller and less consistent although a similar MGP response appeared for T3 exposure. This is likely the result of (slightly) lower sample sizes but, more importantly, a more complex innate regulation and signalling of MGP in bone as multicellular tissue type. As such, the observed variation in the MGP response in bone remains inconclusive and needs to be repeated in larger sample sizes.
      Upon identifying MGP, encoding an inhibitor of ectopic calcifications, as strong OA risk gene, it was hypothesized that the OA risk was conferred via calcification of cartilage tissue
      • Luo G.
      • Ducy P.
      • McKee M.D.
      • Pinero G.J.
      • Loyer E.
      • Behringer R.R.
      • et al.
      Spontaneous calcification of arteries and cartilage in mice lacking matrix GLA protein.
      ,
      • Yagami K.
      • Suh J.Y.
      • Enomoto-Iwamoto M.
      • Koyama E.
      • Abrams W.R.
      • Shapiro I.M.
      • et al.
      Matrix GLA protein is a developmental regulator of chondrocyte mineralization and, when constitutively expressed, blocks endochondral and intramembranous ossification in the limb.
      . Moreover, as MGP protein is activated by vitamin K dependent carboxylation (c-MGP) this finding underscored the relevance of previously found associations between OA and low vitamin K status
      • Neogi T.
      • Booth S.L.
      • Zhang Y.Q.
      • Jacques P.F.
      • Terkeltaub R.
      • Aliabadi P.
      • et al.
      Low vitamin K status is associated with osteoarthritis in the hand and knee.
      ,
      • Shea M.K.
      • Kritchevsky S.B.
      • Hsu F.C.
      • Nevitt M.
      • Booth S.L.
      • Kwoh C.K.
      • et al.
      The association between vitamin K status and knee osteoarthritis features in older adults: the Health, Aging and Body Composition Study.
      . Here, we showed that exposure of the vitamin K inhibitor warfarin to intact human articular cartilage explants provoked unbeneficial functional chondrocyte signalling towards hypertrophy, as reflected by upregulation of COL10A1 and almost significant upregulation of RUNX2 and ALPL. Moreover, we showed a modest but significant downregulation of SOX9, a transcription factor marking healthy articular cartilage. These observed effects of warfarin on chondrocyte signalling were similar to those previously found during in vitro knockdown of MGP in chondrocyte monolayer cultures
      • Shepherd C.
      • Reese A.E.
      • Reynard L.N.
      • Loughlin J.
      Expression analysis of the osteoarthritis genetic susceptibility mapping to the matrix Gla protein gene MGP.
      . With regard to the seemingly increased MMP13 and reduced MMP3 gene expression, it has been suggested that MMP3 plays a role mainly in healthy cartilage remodeling, while MMP13 more so in pathophysiological processes. This was confirmed by performing a look-up in our RNA-sequencing data set
      • Coutinho de Almeida R.
      • Ramos Y.F.M.
      • Mahfouz A.
      • den Hollander W.
      • Lakenberg N.
      • Houtman E.
      • et al.
      RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage.
      were MMP3 showed a marked downregulation in lesioned compared to preserved OA cartilage. Exposure of warfarin to subchondral bone of osteochondral explants also provoked unbeneficial functional signalling towards reduced bone formation, as reflected by downregulation of COL1A1 and the suggestive downregulation of BGLAP, whereas the upregulation of the osteoclastogenesis inhibitor TNFRSF11B (although not significant) and downregulation of the vitamin K dependent protein POSTN suggests altered bone remodeling likely resulting in bone loss
      • Zhang Y.
      • Zhao L.
      • Wang N.
      • Li J.
      • He F.
      • Li X.
      • et al.
      Unexpected role of matrix Gla protein in osteoclasts: inhibiting osteoclast differentiation and bone resorption.
      ,
      • Bonnet N.
      • Gineyts E.
      • Ammann P.
      • Conway S.J.
      • Garnero P.
      • Ferrari S.
      Periostin deficiency increases bone damage and impairs injury response to fatigue loading in adult mice.
      . Due to the low numbers and heterogeneity of patients, future studies are necessary to investigate if rs1800801 genotype influences response of cells to warfarin. In light of our result we advocate that the frequent prescription of warfarin as vitamin K dependent blood anticoagulant
      • Holbrook A.
      • Schulman S.
      • Witt D.M.
      • Vandvik P.O.
      • Fish J.
      • Kovacs M.J.
      • et al.
      Evidence-based management of anticoagulant therapy: antithrombotic therapy and prevention of thrombosis, 9th ed: American college of chest physicians evidence-based clinical practice guidelines.
      may have clinical consequences in evoking OA comorbidity. As such, the risk of OA comorbidity may be considerably reduced by preferred prescription of non-vitamin K antagonist as anticoagulants
      • Zhu J.
      • Alexander G.C.
      • Nazarian S.
      • Segal J.B.
      • Wu A.W.
      Trends and variation in oral anticoagulant choice in patients with atrial fibrillation, 2010-2017.
      . In addition, vitamin K supplementation should be considered a potential novel OA-modifying treatment option. In this respect, there has been one underpowered clinical trial studying the effect of vitamin K supplementation on OA progression. This ancillary study, originally designed to study vascular calcification, reported no overall beneficial effects of vitamin K supplementation. However, in individuals with insufficient vitamin K levels at baseline a beneficial effect was observed
      • Neogi T.
      • Felson D.T.
      • Sarno R.
      • Booth S.L.
      Vitamin K in hand osteoarthritis: results from a randomised clinical trial.
      .
      Although the human aged macroscopically normal osteochondral explants used in our study may represent physiological relevant human articular cartilage and subchondral bone model, prone to OA pathophysiology, hence suitable to study the initial process of OA related destruction, the model is inherently subject to heterogeneity. Moreover, it does not provide insight into the MGP effect of such environmental perturbations to healthy cartilage and bone. Another limiting factor was the low sample size of T3, IL-1β and warfarin treated explants upon stratifying for rs1800801 genotype, resulting in no or less robust results than upon investigating the response in the larger mechanically stressed group. It should also be noted that the modifying effect of the MGP OA risk allele rs1800801-T as function of the OA status in articular cartilage and subchondral bone was only measured as a static effect i.e., differential expression of MGP between paired preserved and lesioned OA cartilage samples. Finally, the focus of our paper was on exploring gene expression changes of MGP only. Although studying protein levels of MGP as function of the OA risk SNP and the OA relevant cues in joint tissue would be an interesting addition and a preferred next step, such analyses should involve the detection of activated (hence carboxylated) MGP protein.
      Together our data highlight that, similar to the bi-directional interplay of vascular calcification and bone loss in osteoporosis and atherosclerosis
      • Vassalle C.
      • Mazzone A.
      Bone loss and vascular calcification: a bi-directional interplay?.
      , also articular cartilage calcification and bone loss in OA might share a common pathogenetic mechanism likely involving MGP. Moreover, warfarin on human osteochondral explant cultures functionally underscores the previously found association between vitamin K deficiency and OA.

      Contributions

      All authors have made contributions to the completion of this study. Study concept and design: EH, RCA, JM, YFM, IM. Acquisition of material and data: EH, MT, HED, DB, RGHHN. Data analysis: EH, RCA, MT, IM. Preparation of the manuscript: EH, IM. Critical reviewing and approval of the manuscript: All authors.

      Conflicts of interest

      None declared.

      Funding source

      The study was funded by the Dutch Scientific Research council NWO/ZonMW VICI scheme ( nr 91816631/528 ), Dutch Arthritis Society/ReumaNederland ( DAF-15-4-401 and DAA_10_1–402 ) and European Commission Seventh Framework program (TreatOA, 200800 ).

      Patient consent for publication

      Not required.

      Data availability statement

      Data are available on reasonable request.

      Acknowledgements

      We thank all the participants of the RAAK study. The LUMC has and is supporting the RAAK study. We also thank Enrike van der Linden, Robert van der Wal, Peter van Schie, Anika Rabelink-Hoogenstraaten, Shaho Hasan, Maartje Meijer, Daisy Latijnhouwers and Geert Spierenburg for collecting surgical waste material. Data is generated within the scope of the Medical Delta Programs Regenerative Medicine 4D: Generating complex tissues with stem cells and printing technology and Improving Mobility with Technology. We thank the Sequence Analysis Support Core (SASC) of the Leiden University Medical Center for their support.

      Appendix A. Supplementary data

      The following are the Supplementary data to this article:

      References

        • Spector T.D.
        • Cicuttini F.
        • Baker J.
        • Loughlin J.
        • Hart D.
        Genetic influences on osteoarthritis in women: a twin study.
        BMJ. 1996; 312: 940-943
        • Felson D.T.
        Risk factors for osteoarthritis: understanding joint vulnerability.
        Clin Orthop Relat Res. 2004; 427: S16-S21
        • Bos S.D.
        • Bovee J.V.
        • Duijnisveld B.J.
        • Raine E.V.
        • van Dalen W.J.
        • Ramos Y.F.
        • et al.
        Increased type II deiodinase protein in OA-affected cartilage and allelic imbalance of OA risk polymorphism rs225014 at DIO2 in human OA joint tissues.
        Ann Rheum Dis. 2012; 71: 1254-1258
        • Reynard L.N.
        • Bui C.
        • Canty-Laird E.G.
        • Young D.A.
        • Loughlin J.
        Expression of the osteoarthritis-associated gene GDF5 is modulated epigenetically by DNA methylation.
        Hum Mol Genet. 2011; 20: 3450-3460
        • Loughlin J.
        Genetic contribution to osteoarthritis development: current state of evidence.
        Curr Opin Rheumatol. 2015; 27: 284-288
        • Shepherd C.
        • Zhu D.
        • Skelton A.J.
        • Combe J.
        • Threadgold H.
        • Zhu L.
        • et al.
        Functional characterisation of the osteoarthritis genetic risk residing at ALDH1A2 identifies rs12915901 as a key target variant.
        Arthritis Rheum. 2018; 70: 1577-1587
        • den Hollander W.
        • Boer C.G.
        • Hart D.J.
        • Yau M.S.
        • Ramos Y.F.M.
        • Metrustry S.
        • et al.
        Genome-wide association and functional studies identify a role for matrix Gla protein in osteoarthritis of the hand.
        Ann Rheum Dis. 2017; 76: 2046-2053
        • Misra D.
        • Booth S.L.
        • Crosier M.D.
        • Ordovas J.M.
        • Felson D.T.
        • Neogi T.
        Matrix Gla protein polymorphism, but not concentrations, is associated with radiographic hand osteoarthritis.
        J Rheumatol. 2011; 38: 1960-1965
        • Shepherd C.
        • Reese A.E.
        • Reynard L.N.
        • Loughlin J.
        Expression analysis of the osteoarthritis genetic susceptibility mapping to the matrix Gla protein gene MGP.
        Arthritis Res Ther. 2019; 21: 149
        • Luo G.
        • Ducy P.
        • McKee M.D.
        • Pinero G.J.
        • Loyer E.
        • Behringer R.R.
        • et al.
        Spontaneous calcification of arteries and cartilage in mice lacking matrix GLA protein.
        Nature. 1997; 386: 78-81
        • Yagami K.
        • Suh J.Y.
        • Enomoto-Iwamoto M.
        • Koyama E.
        • Abrams W.R.
        • Shapiro I.M.
        • et al.
        Matrix GLA protein is a developmental regulator of chondrocyte mineralization and, when constitutively expressed, blocks endochondral and intramembranous ossification in the limb.
        J Cell Biol. 1999; 147: 1097-1108
        • Tunon-Le Poultel D.
        • Cannata-Andia J.B.
        • Roman-Garcia P.
        • Diaz-Lopez J.B.
        • Coto E.
        • Gomez C.
        • et al.
        Association of matrix Gla protein gene functional polymorphisms with loss of bone mineral density and progression of aortic calcification.
        Osteoporos Int. 2014; 25: 1237-1246
        • Sheng K.
        • Zhang P.
        • Lin W.
        • Cheng J.
        • Li J.
        • Chen J.
        Association of Matrix Gla protein gene (rs1800801, rs1800802, rs4236) polymorphism with vascular calcification and atherosclerotic disease: a meta-analysis.
        Sci Rep. 2017; 7: 8713
        • Neogi T.
        • Booth S.L.
        • Zhang Y.Q.
        • Jacques P.F.
        • Terkeltaub R.
        • Aliabadi P.
        • et al.
        Low vitamin K status is associated with osteoarthritis in the hand and knee.
        Arthritis Rheum. 2006; 54: 1255-1261
        • Shea M.K.
        • Kritchevsky S.B.
        • Hsu F.C.
        • Nevitt M.
        • Booth S.L.
        • Kwoh C.K.
        • et al.
        The association between vitamin K status and knee osteoarthritis features in older adults: the Health, Aging and Body Composition Study.
        Osteoarthritis Cartilage. 2015; 23: 370-378
        • Holbrook A.
        • Schulman S.
        • Witt D.M.
        • Vandvik P.O.
        • Fish J.
        • Kovacs M.J.
        • et al.
        Evidence-based management of anticoagulant therapy: antithrombotic therapy and prevention of thrombosis, 9th ed: American college of chest physicians evidence-based clinical practice guidelines.
        Chest. 2012; 141: e152S-e184S
        • Chin K.Y.
        The relationship between vitamin K and osteoarthritis: a review of current evidence.
        Nutrients. 2020; 12
        • Coutinho de Almeida R.
        • Ramos Y.F.M.
        • Mahfouz A.
        • den Hollander W.
        • Lakenberg N.
        • Houtman E.
        • et al.
        RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage.
        Ann Rheum Dis. 2019; 78: 270-277
        • Tuerlings M.
        • van Hoolwerff M.
        • Houtman E.
        • Suchiman E.
        • Lakenberg N.
        • Mei H.
        • et al.
        RNA sequencing reveals interacting key determinants of osteoarthritis acting in subchondral bone and articular cartilage.
        Arthritis Rheum. 2021; 73: 789-799
        • Houtman E.
        • van Hoolwerff M.
        • Lakenberg N.
        • Suchiman E.H.D.
        • van der Linden-van der Zwaag E.
        • Nelissen R.
        • et al.
        Human osteochondral explants: reliable biomimetic models to investigate disease mechanisms and develop personalized treatments for osteoarthritis.
        Rheumatol Ther. 2021; 20: 499-515
        • Ramos Y.F.
        • den Hollander W.
        • Bovee J.V.
        • Bomer N.
        • van der Breggen R.
        • Lakenberg N.
        • et al.
        Genes involved in the osteoarthritis process identified through genome wide expression analysis in articular cartilage; the RAAK study.
        PLoS One. 2014; 9e103056
        • den Hollander W.
        • Pulyakhina I.
        • Boer C.
        • Bomer N.
        • van der Breggen R.
        • Arindrarto W.
        • et al.
        Annotating transcriptional effects of genetic variants in disease-relevant tissue: transcriptome-wide allelic imbalance in osteoarthritic cartilage.
        Arthritis Rheum. 2019; 71: 561-570
        • Zeger S.L.
        • Liang K.Y.
        Longitudinal data analysis for discrete and continuous outcomes.
        Biometrics. 1986; 42: 121-130
        • Diggle P.
        • Liang K.-Y.
        • Zeger S.L.
        Analysis of Longitudinal Data. xi. Clarendon Press; Oxford University Press, Oxford New York1994: 253
        • Zhang Y.
        • Zhao L.
        • Wang N.
        • Li J.
        • He F.
        • Li X.
        • et al.
        Unexpected role of matrix Gla protein in osteoclasts: inhibiting osteoclast differentiation and bone resorption.
        Mol Cell Biol. 2019; 39: e00012-e00019
        • Vassalle C.
        • Mazzone A.
        Bone loss and vascular calcification: a bi-directional interplay?.
        Vasc Pharmacol. 2016; 86: 77-86
        • Bonnet N.
        • Gineyts E.
        • Ammann P.
        • Conway S.J.
        • Garnero P.
        • Ferrari S.
        Periostin deficiency increases bone damage and impairs injury response to fatigue loading in adult mice.
        PLoS One. 2013; 8e78347
        • Zhu J.
        • Alexander G.C.
        • Nazarian S.
        • Segal J.B.
        • Wu A.W.
        Trends and variation in oral anticoagulant choice in patients with atrial fibrillation, 2010-2017.
        Pharmacotherapy. 2018; 38: 907-920
        • Neogi T.
        • Felson D.T.
        • Sarno R.
        • Booth S.L.
        Vitamin K in hand osteoarthritis: results from a randomised clinical trial.
        Ann Rheum Dis. 2008; 67: 1570-1573