Skip to main content
  • Research article
  • Open access
  • Published:

Genome-wide DNA methylation study of hip and knee cartilage reveals embryonic organ and skeletal system morphogenesis as major pathways involved in osteoarthritis



Evidence suggests that epigenetics plays a role in osteoarthrits (OA). The aim of the study was to describethe genome wide DNA methylation changes in hip and knee OA and identify novel genes and pathwaysinvolved in OA by comparing the DNA methylome of the hip and knee osteoarthritic cartilage tissues withthose of OA-free individuals.


Cartilage samples were collected from hip or knee joint replacement patients either due to primary OA or hip fractures as controls. DNA was extracted from the collected cartilage and assayed by Illumina Infinium HumanMethylation450 BeadChip array, which allows for the analysis of >480,000 CpG sites. Student T-test was conducted for each CpG site and those sites with at least 10 % methylation difference and a p value <0.0005 were defined as differentially methylated regions (DMRs) for OA. A sub-analysis was also done for hip and knee OA separately. DAVID v6.7 was used for the functional annotation clustering of the DMR genes. Clustering analysis was done using multiple dimensional scaling and hierarchical clustering methods.


The study included 5 patients with hip OA, 6 patients with knee OA and 7 hip cartilage samples from OA-free individuals. The comparisons of hip, knee and combined hip/knee OA patients with controls resulted in 26, 72, and 103 DMRs, respectively. The comparison between hip and knee OA revealed 67 DMRs. The overall number of the sites after considering the overlaps was 239, among which 151 sites were annotated to 145 genes. One-fifth of these genes were reported in previous studies. The functional annotation clustering of the identified genes revealed clusters significantly enriched in skeletal system morphogenesis and development. The analysis revealed significant difference among OA and OA-free cartilage, but less different between hip OA and knee OA.


We found that a number of CpG sites and genes across the genome were differentially methylated in OA patients, a remarkable portion of which seem to be involved in potential etiologic mechanisms of OA. Genes involved in skeletal developmental pathways and embryonic organ morphogenesis may be a potential area for further OA studies.

Peer Review reports


Osteoarthritis (OA), affecting 250 million people worldwide, is the most common form of arthritis [1]. It is characterized by gradual loss of articular cartilage and subchondral bone changes, presents with joint pain, stiffness, joint deformity and disability [2], and imposes a great socio-economic burden on societies, mainly as a result of hip and/or knee involvement [3].

OA is a multifactorial condition resulting from the combination and interaction of natural and environmental factors [4]. Age, gender, obesity, previous joint injury, mal-alignments, and genetics are known to be of the major risk factors for OA; yet, the etiology of OA remains incompletely elucidated [5]. In pathology, an imbalance between catabolism and anabolism of the molecules in the cartilage extracellular matrix is a major finding in OA [6]. Since these changes are suggested to result from an altered gene expression related to epigenetic modifications of the OA candidate genes [7], it is hypothesized that epigenetic changes in chondrocytes could be a key factor of OA pathogenesis [8]. DNA methylation is by far the most extensively studied epigenetic regulator in complex diseases, and it has long been thought that its changes plays a key role in the onset and progression of complex diseases by linking the genetic and environmental risk factors [9].

To date, only few studies have been undertaken to examine the role of DNA methylation in OA. Candidate gene studies have shown the up-regulation of catabolic factors of Matrix metallopeptidase 9 (MMP9), Matrix metallopeptidase 13 (MMP13), Leptin receptor (LEPR), and A disintegrin and metalloproteinase with thrombospondin motifs 4 (ADAMTS4) [7], as well as the down-regulation of anabolic molecule collagen, type IX, alpha 1 (COL9A1) [10], caused by the promoter hypo- and hyper- methylations of the corresponding genes. Demethylation of an enhancer element within the nitric oxide synthase (NOS) gene is shown to increase transcription through elevated binding of the transcription factor NF-κB, which leads to suppression of the synthesis of cartilage matrix [11]. DNA methylation can also modulate the effect of OA genetic susceptibility loci; for instance, the effect of single nucleotide polymorphism (SNP) rs143383 in GDF5 -- the most replicated genetic association locus in OA - is thought to be caused by the methylation level variability of the CpG dinucleotide created at the location of the SNP, leading to altered expression of the gene [12]. The handful of genome wide methylation studies performed to date have also identified several potential candidate genes including runt-related transcription factor 1&2 (RUNX1, RUNX2), transforming growth factor beta 1 (TGFB1), micro RNA 128 (miR-128) and collagen, type XI, alpha 2 COL11A2 [13], suggesting the involvement of inflammation and immunity in OA pathogenesis [14]. Despite the invaluable information obtained about the pathogenesis of complex diseases from epigenetic studies, the area still remains as one of the least investigated fields in OA research.

In the present study, we conducted a genome wide DNA methylation analysis in OA-free and OA-affected cartilage from human hips and knees using the Illumina Infinium HumanMethylation450 BeadChip in the hope of providing novel insights into the pathogenesis and treatment of OA.


Samples and patient’s information

The study was part of the ongoing Newfoundland Osteoarthritis Study (NFOAS) that was initiated in 2011, aiming at identifying novel genetic, epigenetic, and biochemical markers for OA [1517]. OA patients were recruited from those who underwent total knee or hip joint replacement due to primary OA between November 2011 and December 2013 in St. Clare’s Mercy Hospital and Health Science Centre General Hospital in St. John’s, the capital city of Newfoundland and Labrador (NL), Canada. OA-free controls were recruited in the same hospitals from those who underwent hemiarthroplasty of the hip due to hip fracture with no evidence of OA. OA diagnosis was made based on the American College of Rheumatology criteria [18, 19] and the judgement of the attending orthopaedic surgeons. Cartilage samples were collected from the articular surfaces of the tibial plateau or femoral head where the OA lesion occurred. The pathology report of the cartilage following the surgery was reviewed for all subjects to ensure the consistency of the diagnosis and the status of cartilage degeneration among the control subjects.

Demographic information was obtained by a self-administered questionnaire with the help of the research staff if necessary. Anthropometric data including height and weight was retrieved from their hospital admission and medical records and body mass index (BMI) was calculated by dividing weight in kilograms by squared height in meters. Age was calculated at the time of the surgery.

DNA extraction

Four pieces (~200 mg each) of cartilage tissues were retained from either tibial plateau or femoral heads during the surgery. The samples were then flash frozen and stored in liquid nitrogen until the experiment. Up to 200 mg frozen cartilage tissue was transferred to the homogenizing cylinder together with 1 ml TRIzole lysis reagent and 200 μl guanidine thiocyanate and homogenized using a cryogenic mill (Spex Freezer Mill, model 6770, Metuchen, New Jersey, USA) with the following parameters: two cycles of 2 min grinding at maximum frequency with 10 min cooling down between grinding cycles. The homogenate was then transferred to a new 2 ml RNase free tube and incubated for 5 min at room temperature. Then, 200 μl chloroform was added and the mix was vortexed vigorously, before being incubated for 2–3 min, followed by a centrifugation at 12,000xg at −4 °C for 15 min. Following centrifugation, the sample separated into 3 phases: the aqueous phase containing RNA, the interphase, and the organic phase containing DNA. The DNA was extracted using Phenol-Chloroform method from the interphase and organic phase.

DNA methylation profiling

DNA methylation assay was conducted using the HumanMethylation 450 Bead-Chip microarray (Illumina, San Diego, California, USA), which analyzes the methylation status of 485,000 methylation sites throughout the genome, covering 99 % of RefSeq genes at an average of 17 CpG sites per gene across the 5-UTR, gene promoter regions, first exon, gene body, and 3-UTR, and covering 96 % of University of California, Santa Cruz-defined CpG islands and their flanking regions. Briefly, DNA is first bisulfite converted, which results in unmethylated cytosines being converted to uracils, whereas methylated cytosines are not converted. The bisulfite-converted DNA is amplified, fragmented and hybridized to the arrays. For each CpG site, methylation levels are measured by probes attached to beads, one each for unmethylated and methylated DNA, followed by allele-specific base extension that includes a fluorescent label. Different labels are used for the T (unmethylated) or C (methylated) alleles. The array is fluorescently stained and scanned, and the intensities of the unmethylated and methylated bead types are measured. DNA methylation values, described as “beta values (β)”, are recorded for each locus in each sample and represent the ratio of the intensity of the methylated bead type to the combined locus intensity. The array has several features that make it a powerful option for genome-wide DNA methylation profiling: (i) multi sample format allows for interrogation of 12 samples on a single BeadChip i.e. it is high-throughout and cost-effective; (ii) low sample input (500 ng genomic DNA); (iii) reproducibility (>0.98 between technical replicates) [20, 21]. The genome-wide DNA methylation data is available at with accession number GSE73626.

Statistical analysis

R packages minfi (version 1.6.0) [22] and minfiData (version 0.3.4) [23] was used in R version 3.1.2 to convert signal intensity data into methylation data. β values were calculated as M/(M + U), where M represents the fluorescent signal of the methylation probe and U represents the methylation signal of the unmethylated probe. The β values range from 0 (no methylation) to 1 (100 % methylation) [24].

For the purpose of quality control, CpG probes with detection p-value above 0.01, those located in sex chromosomes and at SNPs, and those with deviation from bimodal distribution were removed from further analysis. In order to eliminate the difference caused by two type design probes (type I and type II), beta-mixture quantile normalization (BMIQ) method [25] was used to normalize the raw methylation level data.

To identify DMRs, the average beta values were compared between the groups of interest using a student T-test assuming equal variances [26]. Given the small sample size, we were not able to correct for multiple testing as none of our tests reached the strict threshold of genome wide significance. Instead, to minimize false positives, we reported the loci with at least 10 % methylation difference and p-value ≤ 0.0005. The identified loci were examined using a linear regression model to determine if they were associated with age. Genomic annotation of DMRs was carried out using the Infinium HumanMethylation450 BeadChip annotation file (

Gene ontology analysis

DAVID bioinformatics database functional tool [27] was used to identify the enriched gene ontology (GO) terms. Gene symbols were used as input for the analysis. Medium classification stringency was used. Enriched GO terms with a Bonferroni corrected p-value of less than 0.05 were reported.

Phenotype clustering

Hierarchical clustering and multiple dimensional scaling were performed on all individuals after genome wide methylation distance reduction. The distance was calculated between samples taking into account the genome wide methylation levels of each sample. Then a matrix containing the pairwise distances between the samples was created. Using singular value decomposition, the matrix was transformed into three matrixes, two of which were orthogonal (U and V) and one was diagonal (D). A new matrix was generated using the multiplication of V and D, containing 12 rows representing each sample and 24 columns representing each dimension. Every two dimensions were plotted against each other. Clustering were further illustrated using heat maps of the top 800 loci with the greatest variations in methylation levels across the entire study population after the same methodology of dimension reduction. To obtain these 800 sites, the cross population variance for each CpG site was calculated, the clustering was performed for the top 200, 400, 600, 800, and above sites with the highest variance. The top 800 samples resulted in the best visual grouping of the three phenotypes while the numbers beyond this figure did not change the pattern.

Ethics statement

The study protocol was approved by the Health Research Ethics Authority (HREA) of Newfoundland and Labrador and a written informed consent was obtained from all the participants.



The study included cartilage samples from 5 patients with hip OA, 6 patients with knee OA, and 7 OA-free hip controls. All subjects were females. The OA-free population were on average about 10 years older than the affected group and had a lower BMI. Table 1 shows the characteristics of the study population.

Table 1 Characteristics of study population

Differentially methylated loci

A total of 384,266 CpG sites were included in the final analysis after quality control. A total of 72, 26, and 103 CpG sites was identified from the comparison of knee OA, hip OA, and combined knee/hip OA versus hip controls, respectively. The comparison of hip and knee OA resulted in 67 CpG sites. After removing the overlaps between these analyses, a total of 239 CpG sites showed more than 10 % difference in β values among the comparison groups with all p <0.0005. Methylation levels of these sites were not associated with age. Almost half of the sites (53 %) showed hypo-methylation and the remainders represented hyper-methylation in OA compared to controls. This was, however, reversed among the sites with the highest methylation difference, since most of them were hypo-methylated in OA as seen in Table 2. Additional file 1: Table S1 presents the β value differences in each comparison. Among the reported sites, 151 sites were annotated to 145 genes; 119 of the sites are located in CpG islands; 79 sites are located in enhancers, 46 in regions with regulatory features, and 28 in DNAse Hypersensitivity sites. From the sites annotated to genes, the majority (46 %) were located in gene bodies, 11 % were in 5’UTR, 11 % were in 3’UTR, 5 % were in the first exon, and 27 % were located within 1500 bp upstream the transcription start site. Table 2 shows the DMRs with β value differences above 15 % between knee/hip OA and OA-free cartilage. The complete list of the DMRs is presented in the Additional file 1: Table S1.

Table 2 Top CpG sites differentially methylated in knee/hip OA compared to OA-free cartilagea

Functional annotation clustering of the differentially methylated loci

The annotation clustering was conducted for the 145 genes. Table 3 shows those GO terms with a Bonferroni corrected p-value < 0.05. The most significant terms are related to skeletal and embryonic organ system development and homeobox (HOX). By increasing the classification stringency, the GO terms were classified into 31 clusters, among which only two yielded significant Bonferroni corrected p-values including Embryoinic and skeletal system development and HOX genes (enrichment scores 6.48 and 5.52, respectively). The analysis was repeated after the removal of 33 genes which were only identified in the comparison of knee OA and hip OA. Similar to previous clustering, the GO terms included skeletal and embryonic system development, but no HOX genes.

Table 3 Enrichment clustering of the differentially methylated genesa

Clustering of hip OA, knee OA and OA-free cartilage

We used multiple dimensional scaling and hierarchical clustering to classify the three phenotypes in the study, i.e. hip OA, knee OA, and OA-free hip cartilage. Due to the small sample size, the classification was not perfect; however, some trends were observed which are worthy of consideration. As it is seen in the plots (Fig. 1), samples from each phenotype tend to cluster together, although a few outliers exist. Overall, OA-free hip cartilage samples tend to be different from hip OA and knee OA samples. Although hip OA and knee OA samples are very close together, the similarity to OA-free hip cartilage is more seen in hip OA rather than knee OA samples. Scaling beyond the 2nd dimension was not informative (not shown). Similar patterns are observed from the heat map and dendogram as shown in Fig. 2.

Fig. 1
figure 1

Multiple dimensions scaling of hip OA, knee OA, and OA-free hip cartilage. *Similarities between hip OA, knee OA, and OA-free cartilage, drawn from log-spectral decompositions for each subject as represented in the two-dimensional space by multiple dimensional scaling (MDS). Each dot represents one sample. Colors represent the type of involvement and the site samples obtained. X- and Y- axes represent the first and the second dimension reductions

Fig. 2
figure 2

Hierarchical clustering and heat map of hip OA, knee OA, and OA-free controls*. *Top: Cluster dendogram was created using the genome wide information; Bottom: Heat map shows the top 800 CpG sites with the most variation across hip OA, knee OA, and OA-free hip cartilage samples. Rows represent CpG sites. Columns represent samples. Dark blue indicates hypermethylation and light blue/white indicates hypomethylation


Our study is one of the few reports on the status of genome wide methylation of DNA from OA-free and OA affected human cartilage. We found a number of CpG sites differentially methylated in hip and knee OA, identified the pathways enriched in the sites, and attempted to classify hip and knee OA and OA-free cartilage according to their genome wide DNA methylation profiling.

The majority of the CpG sites we identified were novel and only about one fifth of them were reported by the previous epigenetic studies of OA (Additional file 1: Table S2) [13, 14, 28, 29]. Similarly, most of the genes differentially methylated were not known to play a role in OA, although several of them were previously reported as candidate genes to OA or other bone metabolic conditions. They include cg07902192 in RUNX2, involved in the regulation of matrix metallopeptidase 13 in OA cartilage [30], cg03050981 in LEPR, associated with knee OA [31], bone marrow density and bone hemostasis [32], cg05516020 in CLCN7, associated with bone marrow density [33], cg17279365 in ESRRG, associated with multiple bone disease phenotypes [34], cg14340103 in IL21, being upregulated in synovial biopsies of rheumatoid arthritis patients [35], cg18637380 in MTHFD1, associated with response to osteosarcoma chemotherapy [36], cg10629004 in PAX1, associated with congenital scoliosis [37], and cg10908116 in the alpha-1 subunit of collagen type IV gene (COL4A1).

Although most of the DMR genes in our study were of unknown significance in OA, the functional analysis revealed their enrichment in relevant pathways; i.e. skeletal and embryonic organ system development and homeobox. The latter was only found from the genes differentially methylated in the comparison of knee OA and hip OA methylation. Since homeobox genes are responsible for the body segmentation procedure and specification of lower limbs from upper limbs, it is likely that the DMRs from knee OA and hip OA comparison represent the processes required for body segment specifications rather than the underlying genetic difference between hip OA and knee OA. The same gene ontology term has also been reported by Den Hollander et al. who made a comparison between genome wide methylation of articular cartilage DNA from hip OA and knee OA [28]. The small number of genome wide methylation studies of OA have reported inflammation and immunity [14, 29], transforming growth factor beta signalling [13], and KEGG and developmental pathways [38]. Our study, however, points out the involvement of skeletal system development in OA , which is in accordance with findings of Delgado-Calle J et al. who reported the enrichment of genes associated with the development of the appendicular skeleton and limb morphogenesis in a genome wide methylation study of femoral bone [39]. In this process the anatomical and physical structures of the skeleton are generated and organized. Skeletal shape, which is tightly regulated by genetics [40], is suggested as a possible mechanism for the influence of genetics in OA incidence. An abnormal center-edge angle and acetabular dysplasia are shown to be associated with an increased risk of hip OA [41], and a significant difference in the shape of the intercondylar notch between the OA and non-OA individuals is reported [42]. Wnt signalling and bone morphogenetic proteins are among the pathways involved in OA which also control skeletal development in animal models, and it is suspected that their mechanism of action in OA could be due to their effect on skeletal shape [43].

Our clustering analysis of the three phenotypes in the study clearly shows that OA-affected cartilage has a distinct methylation profiling compared to OA-free cartilage. The hip OA and knee OA clustering, however, is not perfect and is suggesting that although hip OA and knee OA could have different epigenomic landscape, they are very similar to each other. The only minor differences observed in hip OA and knee OA might only be joint specific differences that could have been observed if OA-free knee and hip cartilage were studied. The only comparison of the hip OA and knee OA was done by Den Hollander et al. [37] who successfully grouped hip OA and knee OA into separate clusters. The main conclusion in that study was based on the OA-affected cartilage, and similar to our study the major pathways they identified enriched in the DMRs between hip OA and knee OA, was homeobox, which strengthens the hypothesis that the observed differences in hip and knee OA might be due to differences in the joints rather than the disease status.

Our study is limited by several factors. Due to technical issues we studied a small sample size, and as the result, none of the DMRs reported reached Bonferroni corrected significance. The cases and controls were different in age, which might influence the status of DNA methylation. The controls were selected from the population of patients with possible osteoporosis who may not necessarily represent healthy cartilage epigenome. In addition, we did not have information on the pathological scoring of the OA joints, which could partially explain why the clustering was not perfect. Due to financial issues we did not validate the methylation levels using alternative methods such as pyrosequencing and did not perform functional experiments to add more to the mechanism of involvement. These limitations will likely tackled by future studies attempting to replicate and further studying our findings. Despite these limitations, the trends observed in the study were informative and add to the current knowledge on the pathogenesis of OA.


Through a genome wide methylation study of OA-free and OA-affected human cartilage, we were able to identify a number of CpG sites with methylation changes in OA. We also reported that genes involved in skeletal system morphogenesis are differentially methylated in OA and might be candidate genes for further OA studies. We found a small difference between the overall landscapes of hip OA and knee OA; however, OA-free hip cartilage samples were significantly different from the OA-affected ones. Our findings shed light on the pathophysiology of OA and can pave the road for further research in the field.





A disintegrin and metalloproteinase with thrombospondin motifs 4


Anabolic molecule collagen, type IX, alpha 1


Nitric oxide synthase


Single nucleotide polymorphism


RUNX2, Runt-related transcription factor 1&2


Transforming growth factor beta 1


Micro RNA 128


Collagen, type XI, alpha 2


Newfoundland Osteoarthritis Study


Newfoundland and Labrador


Body mass index


Gene ontology




Alpha-1 subunit of collagen type IV gene


  1. Vos T, Flaxman AD, Naghavi M, Lozano R, Michaud C, Ezzati M, et al. Years lived with disability (YLDs) for 1160 sequelae of 289 diseases and injuries 1990–2010: a systematic analysis for the Global Burden of Disease Study 2010. Lancet. 2012;380:2163–96.

    Article  PubMed  Google Scholar 

  2. Kean WF, Kean R, Buchanan WW. Osteoarthritis: symptoms, signs and source of pain. Inflammopharmacology. 2004;12:3–31.

    Article  CAS  PubMed  Google Scholar 

  3. Guccione AA, Felson DT, Anderson JJ, Anthony JM, Zhang Y, Wilson PW, et al. The effects of specific medical conditions on the functional limitations of elders in the Framingham Study. Am J Public Health. 1994;84:351–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Felson DT. An update on the pathogenesis and epidemiology of osteoarthritis. Radiol Clin North Am. 2004;42:1–9.

    Article  PubMed  Google Scholar 

  5. Hunter DJ, March L, Sambrook PN. Knee osteoarthritis: the influence of environmental factors. Clin Exp Rheumatol. 2002;20:93–100.

    CAS  PubMed  Google Scholar 

  6. Martel-Pelletier J, Boileau C, Pelletier JP, Roughley PJ. Cartilage in normal and osteoarthritis conditions. Best Pract Res Clin Rheumatol. 2008;22:351–84.

    Article  CAS  PubMed  Google Scholar 

  7. Zhai G, Aref-Eshghi E. Biomarkers for osteoarthritis: investigation, identification, and prognosis. Curr Biomark Find. 2012;2:19–28.

    Article  Google Scholar 

  8. Reynard LN, Loughlin J. Genetics and epigenetics of osteoarthritis. Maturitas. 2012;71:200–4.

    Article  CAS  PubMed  Google Scholar 

  9. Rodriguez-Paredes M, Esteller M. Cancer epigenetics reaches mainstream oncology. Nat Med. 2011;17:330–9.

    Article  CAS  PubMed  Google Scholar 

  10. Imagawa K, de Andrés MC, Hashimoto K, Itoi E, Otero M, Roach HI, et al. Association of reduced type IX collagen gene expression in human osteoarthritic chondrocytes with epigenetic silencing by DNA hypermethylation. Arthritis Rheumatol. 2014;66:3040–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. De Andres MC, Imagawa K, Hashimoto K, Gonzalez A, Roach HI, Goldring MB, et al. Loss of methylation in CpG sites in the NF-κB enhancer elements of inducible nitric oxide synthase is responsible for gene induction in human articular chondrocytes. Arthritis Rheum. 2013;65:732–42.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Reynard LN, Bui C, Canty-Laird EG, Young DA, Loughlin J. Expression of the osteoarthritis-associated gene GDF5 is modulated epigenetically by DNA methylation. Hum Mol Genet. 2011;20:3450–60.

    Article  CAS  PubMed  Google Scholar 

  13. Jeffries MA, Donica M, Baker LW, Stevenson ME, Annan AC, Humphrey MB, et al. Genome-wide DNA methylation study identifies significant epigenomic changes in osteoarthritic cartilage. Arthritis Rheumatol. 2014;66:2804–15.

    Article  CAS  PubMed  Google Scholar 

  14. Rushton MD, Reynard LN, Barter MJ, Refaie R, Rankin KS, Young DA, et al. Characterization of the cartilage DNA methylome in knee and hip osteoarthritis. Arthritis Rheumatol. 2014;66:2450–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Aref-Eshghi E, Rahman P, Zhang H, Martin G, Furey A, Green R, et al. Attempt to replicate the published osteoarthritis-associated genetic variants in the Newfoundland& Labrador Population. J Orthopedics Rheumatol. 2014;1:5.

    Google Scholar 

  16. Zhang W, Likhodii S, Aref-Eshghi E, Zhang Y, Harper PE, Randell E, et al. Relationship between blood plasma and synovial fluid metabolite concentrations in patients with osteoarthritis. J Rheumatol. 2015;42:859–65.

    Article  CAS  PubMed  Google Scholar 

  17. Zhang W, Likhodii S, Zhang Y, Aref-Eshghi E, Harper PE, Randell E, et al. Classification of osteoarthritis phenotypes by metabolomics analysis. BMJ Open. 2014;4:e006286. doi:10.1371/journal.pone.0097786.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Altman R, Alarcón G, Appelrouth D, Bloch D, Borenstein D, Brandt K, et al. The American College of Rheumatology criteria for the classification and reporting of osteoarthritis of the hip. Arthritis Rheum. 1991;34:505–14.

    Article  CAS  PubMed  Google Scholar 

  19. Altman R, Asch E, Bloch D, Bole G, Borenstein D, Brandt K, et al. Development of the criteria for the classification and reporting of osteoarthritis, classification of osteoarthritis of the knee. Arthritis Rheum. 1986;29:1039–49.

    Article  CAS  PubMed  Google Scholar 

  20. Dedeurwaerder S, Defrance M, Calonne E, Denis H, Sotiriou C, Fuks F. Evaluation of the Infinium Methylation 450 K technology. Epigenomics. 2011;3:771–84.

    Article  CAS  PubMed  Google Scholar 

  21. Bibikova M, Barnes B, Tsan C, Ho V, Klotzle B, Le JM, et al. High density DNA methylation array with single CpG site resolution. Genomics. 2011;98(4):288–95.

    Article  CAS  PubMed  Google Scholar 

  22. Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP, Hansen KD, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics. 2014;30(10):1363–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Hansen KD, Aryee M, Timp W. minfiData: Example data for the Illumina Methylation 450 k array. R package version 0.7, 1.

  24. Du P, Zhang X, Huang CC, Jafari N, Kibbe WA, Hou L, et al. Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC bioinformatics. 2010;11(1):587.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Teschendorff AE, Marabita F, Lechner M, Bartlett T, Tegner J, Gomez-Cabrero D, et al. A beta-mixture quantile normalization method for correcting probe design bias in Illumina Infinium 450 k DNA methylation data. Bioinformatics. 2013;29:189–96.

    Article  CAS  PubMed  Google Scholar 

  26. Laird PW. Principles and challenges of genome-wide DNA methylation analysis. Nat Rev Genet. 2010;11:191–203.

    Article  CAS  PubMed  Google Scholar 

  27. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID Bioinformatics Resources. Nature Protoc. 2009;4:44–57.

    Article  CAS  Google Scholar 

  28. den Hollander W, Ramos YF, Bos SD, Bomer N, van der Breggen R, Lakenberg N, et al. Knee and hip articular cartilage have distinct epigenomic landscapes: implications for future cartilage regeneration approaches. Ann Rheum Dis. 2014;73:2208–12.

    Article  Google Scholar 

  29. Fernández-Tajes J, Soto-Hermida A, Vázquez-Mosquera ME, Cortés-Pereira E, Mosquera A, Fernández-Moreno M, et al. Genome-wide DNA methylation analysis of articular chondrocytes reveals a cluster of osteoarthritic patients. Ann Rheum Dis. 2014;73:668–77.

    Article  PubMed  Google Scholar 

  30. Wang X, Manner PA, Horner A, Shum L, Tuan RS, Nuckolls GH. Regulation of MMP-13 expression by RUNX2 and FGF2 in osteoarthritic cartilage. Osteoarthritis Cartilage. 2004;12:963–73.

    Article  PubMed  Google Scholar 

  31. Ma XJ, Guo HH, Hao SW, Sun SX, Yang XC, Yu B, et al. Association of single nucleotide polymorphisms (SNPs) in leptin receptor gene with knee osteoarthritis in the Ningxia Hui population. Yi Chuan. 2013;35:359–64.

    Article  CAS  PubMed  Google Scholar 

  32. Ye XL, Lu CF. Association of polymorphisms in the leptin and leptin receptor genes with inflammatory mediators in patients with osteoporosis. Endocrine. 2013;44:481–8.

    Article  CAS  PubMed  Google Scholar 

  33. Pettersson U, Albagha OM, Mirolo M, Taranta A, Frattini A, McGuigan FE, et al. Polymorphisms of the CLCN7 gene are associated with BMD in women. J Bone Miner Res. 2005;20:1960–7.

    Article  CAS  PubMed  Google Scholar 

  34. Elfassihi L, Giroux S, Bureau A, Laflamme N, Cole DE, Rousseau F. Association with replication between estrogen-related receptor gamma (ESRRgamma) polymorphisms and bone phenotypes in women of European ancestry. J Bone Miner Res. 2010;25:901–11.

    CAS  PubMed  Google Scholar 

  35. Ettinger R, Kuchen S, Lipsky PE. Interleukin 21 as a target of intervention in autoimmune disease. Ann Rheum Dis. 2008;67 Suppl 3:iii83–6.

    CAS  PubMed  Google Scholar 

  36. Windsor RE, Strauss SJ, Kallis C, Wood NE, Whelan JS. Germline genetic polymorphisms may influence chemotherapy response and disease outcome in osteosarcoma: a pilot study. Cancer. 2012;118:1856–67.

    Article  CAS  PubMed  Google Scholar 

  37. Fei Q, Wu ZH, Yuan SM, Wang H, Zhou X, Liu Z, et al. Association of PAX1 gene polymorphisms with susceptibility to congenital scoliosis in Chinese Han population. Zhonghua Yi Xue Za Zhi. 2008;88:2597–602.

    CAS  PubMed  Google Scholar 

  38. Moazedi-Fuerst FC, Hofner M, Gruber G, Weinhaeusel A, Stradner MH, Angerer H, et al. Epigenetic differences in human cartilage between mild and severe OA. J Orthop Res. 2014;32:1636–45.

    Article  CAS  PubMed  Google Scholar 

  39. Delgado-Calle J, Fernández AF, Sainz J, Zarrabeitia MT, Sañudo C, García-Renedo R, et al. Genome-wide profiling of bone reveals differentially methylated regions in osteoporosis and osteoarthritis. Arthritis Rheum. 2013;65:197–205.

    Article  CAS  PubMed  Google Scholar 

  40. Wagner EF, Karsenty G. Genetic control of skeletal development. Curr Opin Genet Dev. 2001;11:527–32.

    Article  CAS  PubMed  Google Scholar 

  41. Lane NE, Lin P, Christiansen L, Gore LR, Williams EN, Hochberg MC, et al. Association of mild acetabular dysplasia with an increased risk of incident hip osteoarthritis in elderly white women: the Study of Osteoporotic Fractures. Arthritis Rheum. 2000;43:400–4.

    Article  CAS  PubMed  Google Scholar 

  42. Shepstone L, Rogers J, Kirwan JR, Silverman BW. Shape of the intercondylar notch of the human femur: a comparison of osteoarthritic and non-osteoarthritic bones from a skeletal sample. Ann Rheum Dis. 2001;60:968–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Valdes AM, Spector TD. The contribution of genes to osteoarthritis. Rheum Dis Clin North Am. 2008;34:581–603.

    Article  PubMed  Google Scholar 

Download references


We thank all the study participants who made this study possible, and all the staff in the hospitals who helped us in the collection of samples. Epigenotyping was done by McGill University and Génome Québec Innovation Centre. The study was supported by Canadian Institutes of Health Research, Research Development Corporation of Newfoundland & Labrador, and Memorial University of Newfoundland.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Guangju Zhai.

Additional information

Competing interest

The authors do not have any competing interest to declare.

Authors’ contribution

Study Design: EA, GZ; Data collection: EA, ML, PH, GM, AF, GZ; Statistical analysis: EA, YZ, GZ; Interpretation of the results: EA, GM, AF, RG, GS, PR, GZ; Manuscript writing: EA, GZ; Critical comments on the manuscript: GM, AF, RG, GS, PR. All authors read and approved the final manuscript.

Authors’ information

Roger Green died August 9th, 2015.

Roger Green deceased.

Additional file

Additional file 1: Table S1.

CpG sites differentially methylated in knee/hip OA compared to healthy cartilage*. Table 2. The genes and CpG sites commonly reported between our study and previous epigenome-wide studies of OA. (DOCX 36 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Aref-Eshghi, E., Zhang, Y., Liu, M. et al. Genome-wide DNA methylation study of hip and knee cartilage reveals embryonic organ and skeletal system morphogenesis as major pathways involved in osteoarthritis. BMC Musculoskelet Disord 16, 287 (2015).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: