Open Access
Open Peer Review

This article has Open Peer Review reports available.

How does Open Peer Review work?

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

  • Erfan Aref-Eshghi1,
  • Yuhua Zhang1,
  • Ming Liu1,
  • Patricia E. Harper1,
  • Glynn Martin2,
  • Andrew Furey2,
  • Roger Green1,
  • Guang Sun3,
  • Proton Rahman3 and
  • Guangju Zhai1, 4Email author
BMC Musculoskeletal Disorders201516:287

https://doi.org/10.1186/s12891-015-0745-5

Received: 15 June 2015

Accepted: 2 October 2015

Published: 9 October 2015

Abstract

Background

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.

Methods

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.

Results

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.

Conclusions

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.

Keywords

Osteoarthritis Epigenome DNA methylation Skeletal morphogenesis Hip Knee

Background

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.

Methods

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 http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE73626 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 (http://www.illumina.com/products/methylation_450_beadchip_kits.html).

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.

Results

Subjects

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

 

OA-free Hip

Knee OA

P

Hip OA

P

Age

78.2 ± 11.6

65.3 ± 10.6

0.06

64.4 ± 13.8

0.09

BMI

26.0 ± 4.6

34.6 ± 8.3

0.03

32.3 ± 9.5

0.15

Figures are Mean ± standard deviation; the p-vales were obtained using student’s T-test between OA affected and OA-free individuals

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

CpG

Δβ

P-value

Gene symbol

UCSC location group

UCSC island group

Enhancer

cg22669656

−0.22

0.0004

PGS1

Body

 

Yes

cg11905061

0.21

0.0004

AGAP1

Body

S_Shore

 

cg27390206

−0.21

0.0002

BLMH

Body

 

Yes

cg09140531

−0.21

0.0004

   

Yes

cg14223856

−0.21

0.0005

   

Yes

cg13688786

−0.20

0.0003

MYO18A

Body

 

Yes

cg22022821

−0.20

0.0002

    

cg10340048

−0.20

0.0002

   

Yes

cg02464866

−0.20

0.0004

  

N_Shore

Yes

cg05033952

−0.20

0.0001

  

N_Shore

 

cg19629120

−0.19

0.0005

EFCAB6

3'UTR

  

cg04973183

−0.19

0.0004

   

Yes

cg00150785

−0.19

0.0004

  

N_Shore

Yes

cg12027254

−0.19

0.00004

TNRC6C

Body

 

Yes

cg14068309

−0.19

0.0002

EIF2B1

3'UTR

  

cg13556934

−0.18

0.0003

  

N_Shelf

Yes

cg14022778

−0.18

0.0004

FHAD1

Body

  

cg02017450

0.18

0.0004

   

Yes

cg23074762

−0.18

0.0002

CHSY1

Body

 

Yes

cg04228742

−0.18

0.0002

    

cg22203890

−0.18

0.0005

    

cg17611936

0.18

0.0003

PRKAG2

Body

  

cg07107113

−0.17

0.0004

FBLIM1

5'UTR

S_Shore

 

cg12582728

−0.17

0.0004

   

Yes

cg07404223

−0.17

0.0002

   

Yes

cg25002179

−0.16

0.0005

STARD13

5'UTR;Body;TSS1500

 

Yes

cg17025149

−0.16

0.0001

   

Yes

cg26043955

−0.16

0.0004

   

Yes

cg26919145

−0.16

0.0002

LDLRAD3

Body

 

Yes

cg09425279

−0.16

0.0003

  

N_Shelf

 

cg06712559

0.16

0.0004

AGRN

Body

Island

Yes

cg02099390

−0.16

0.0002

OSBPL10

Body

 

Yes

cg11805414

−0.16

0.0005

   

Yes

cg04038680

−0.16

0.0002

SHISA9

Body

 

Yes

cg25341923

−0.15

0.0001

KRTAP4-7

TSS1500

  

cg14728071

−0.15

0.00004

MLLT10

3'UTR

  

cg03667871

−0.15

0.0004

NEK7

TSS1500

N_Shore

 

cg13258453

−0.15

0.0002

   

Yes

cg23010507

−0.15

0.0001

  

S_Shelf

Yes

cg12158488

−0.15

0.0002

   

Yes

aΔβ: difference in methylation value between sample groups (OA cartilage - intact); UCSC: University of California, Santa Cruz; 5’-UTR: 5’-untranslated region; N: North; S: South; TSS 200: within 200 bp of transcription start site

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

Term

Gene count

% of the total genes entered

Genes in each pathway from our result

Fold Enrichment

Bonferroni P-value

GO:0048705 ~ skeletal system morphogenesis

14

9.6

GSC, TBX15, PAX1, GLI3, HOXB3, HOXD9, HOXC8, HOXD8, HOXC9, OSR2, HOXD3, HOXB6, ALX4, RUNX2

14.96

8.39E-09

GO:0048704 ~ embryonic skeletal system morphogenesis

10

6.9

HOXD9, HOXB3, GSC, OSR2, TBX15, HOXC9, HOXD3, HOXB6, ALX4, GLI3

21.00

9.89E-07

GO:0001501 ~ skeletal system development

17

11.7

CYP24A1, GSC, TBX15, GLI3, PAX1, GLI1, HOXD9, HOXB3, HDAC4, HOXC8, HOXD8, HOXC9, OSR2, HOXD3, HOXB6, ALX4, RUNX2

6.38

8.98E-06

GO:0048562 ~ embryonic organ morphogenesis

12

8.3

HOXD9, HOXB3, GSC, OSR2, TBX15, HOXC9, HOXD3, HOXB6, VAX2, ALX4, GLI3, GLI1

10.80

1.43E-05

GO:0048706 ~ embryonic skeletal system development

10

6.8

HOXD9, HOXB3, GSC, OSR2, TBX15, HOXC9, HOXD3, HOXB6, ALX4, GLI3

15.54

1.55E-05

GO:0003002 ~ regionalization

13

8.9

GSC, VAX2, PAX1, GLI3, GLI1, HOXB3, HOXD9, HOXC8, HOXC9, HOXD8, HOXD3, HOXB6, ALX4

7.91

9.34E-05

GO:0043565 ~ sequence-specific DNA bindingb

20

13.8

GSC, ESRRG, VAX2, MEIS1, GLI3, GLI1, HOXD9, HOXB3, HDAC4, HOXC8, HOXD8, MEIS2, HOXC9, HAND2, GATA6, HOXD3, HOXB6, THAP1, ALX4, ETV6

4.03

1.02E-04

GO:0048568 ~ embryonic organ development

12

8.3

HOXD9, HOXB3, GSC, OSR2, TBX15, HOXC9, HOXD3, HOXB6, VAX2, ALX4, GLI3, GLI1

8.35

2.05E-04

Homeoboxb

12

8.3

HOXD9, HOXB3, HOXC8, GSC, HOXD8, MEIS2, HOXC9, HOXD3, HOXB6, VAX2, ALX4, MEIS1

6.96

3.14E-04

GO:0007389 ~ pattern specification process

14

9.6

GSC, VAX2, PAX1, GLI3, GLI1, HOXB3, SEMA5A, HOXD9, HOXC8, HOXD8, HOXC9, HOXD3, HOXB6, ALX4

6.27

3.57E-04

SM00389:HOXb

12

8.3

HOXD9, HOXB3, HOXC8, GSC, HOXD8, MEIS2, HOXC9, HOXD3, HOXB6, VAX2, ALX4, MEIS1

5.65

5.84E-04

IPR017970:Homeobox, conserved siteb

12

8.3

HOXD9, HOXB3, HOXC8, GSC, HOXD8, MEIS2, HOXC9, HOXD3, HOXB6, VAX2, ALX4, MEIS1

6.62

6.07E-04

IPR001356:Homeoboxb

12

8.3

HOXD9, HOXB3, HOXC8, GSC, HOXD8, MEIS2, HOXC9, HOXD3, HOXB6, VAX2, ALX4, MEIS1

6.54

6.87E-04

IPR012287:Homeodomain-relatedb

12

8.3

HOXD9, HOXB3, HOXC8, GSC, HOXD8, MEIS2, HOXC9, HOXD3, HOXB6, VAX2, ALX4, MEIS1

6.46

7.77E-04

GO:0009952 ~ anterior/posterior pattern formationb

10

6.9

HOXD9, HOXB3, HOXC8, HOXD8, HOXC9, HOXD3, HOXB6, ALX4, PAX1, GLI3

8.55

0.002

DNA-binding region:Homeoboxb

10

6.9

HOXD9, HOXB3, HOXC8, GSC, HOXD8, HOXC9, HOXD3, HOXB6, VAX2, ALX4

7.34

0.004

IPR001827:Homeobox protein, antennapedia type, conserved siteb

5

3.4

HOXB3, HOXC8, HOXD8, HOXD3, HOXB6

26.69

0.01

GO:0043009 ~ chordate embryonic developmentb

13

8.9

GSC, TBX15, PAX1, GLI3, HOXB3, HOXD9, HOXC9, OSR2, HAND2, GATA6, HOXD3, HOXB6, ALX4

4.70

0.02

GO:0009792 ~ embryonic development ending in birth or egg hatchingb

13

8.9

GSC, TBX15, PAX1, GLI3, HOXB3, HOXD9, HOXC9, OSR2, HAND2, GATA6, HOXD3, HOXB6, ALX4

4.66

0.02

GO:0003700 ~ transcription factor activityb

21

14.5

TBX15, GSC, ESRRG, VAX2, MEIS1, GLI3, GLI1, HOXD9, HOXB3, HOXC8, HOXD8, MEIS2, HOXC9, HAND2, GATA6, HOXD3, MLLT10, HOXB6, ALX4, ETV6, RUNX2

2.64

0.02

Developmental proteinb

17

11.7

GSC, ZNF521, VAX2, MEIS1, PAX1, GLI1, HOXD9, HOXB3, SEMA5A, HOXC8, HOXD8, HOXC9, HAND2, HOXD3, HOXB6, ROBO2, ALX4

3.06

0.03

aAll of the GO terms above were clustered into one annotation cluster with an overall enrichment score of 3.95. bThe GO terms were only significant before the removal of the genes differentially methylated between hip OA and knee OA

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

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

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

Discussion

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.

Conclusion

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.

Notes

Abbreviations

OA: 

Osteoarthritis

ADAMTS4: 

A disintegrin and metalloproteinase with thrombospondin motifs 4

COL9A1: 

Anabolic molecule collagen, type IX, alpha 1

NOS: 

Nitric oxide synthase

SNP: 

Single nucleotide polymorphism

RUNX1: 

RUNX2, Runt-related transcription factor 1&2

TGFB1: 

Transforming growth factor beta 1

miR-128: 

Micro RNA 128

COL11A2: 

Collagen, type XI, alpha 2

NFOAS: 

Newfoundland Osteoarthritis Study

NL: 

Newfoundland and Labrador

BMI: 

Body mass index

GO: 

Gene ontology

HOX: 

Homeobox

COL4A1: 

Alpha-1 subunit of collagen type IV gene

Declarations

Acknowledgements

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.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), 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 (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Discipline of Genetics, Faculty of Medicine, Memorial University of Newfoundland
(2)
Division of Orthopedics, Faculty of Medicine, Memorial University of Newfoundland
(3)
Disicpline of Medicine, Faculty of Medicine, Memorial University of Newfoundland
(4)
Department of Twin Research & Genetic Epidemiology, King’s College London

References

  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.View ArticlePubMedGoogle Scholar
  2. Kean WF, Kean R, Buchanan WW. Osteoarthritis: symptoms, signs and source of pain. Inflammopharmacology. 2004;12:3–31.View ArticlePubMedGoogle 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.View ArticlePubMedPubMed CentralGoogle Scholar
  4. Felson DT. An update on the pathogenesis and epidemiology of osteoarthritis. Radiol Clin North Am. 2004;42:1–9.View ArticlePubMedGoogle Scholar
  5. Hunter DJ, March L, Sambrook PN. Knee osteoarthritis: the influence of environmental factors. Clin Exp Rheumatol. 2002;20:93–100.PubMedGoogle 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.View ArticlePubMedGoogle Scholar
  7. Zhai G, Aref-Eshghi E. Biomarkers for osteoarthritis: investigation, identification, and prognosis. Curr Biomark Find. 2012;2:19–28.View ArticleGoogle Scholar
  8. Reynard LN, Loughlin J. Genetics and epigenetics of osteoarthritis. Maturitas. 2012;71:200–4.View ArticlePubMedGoogle Scholar
  9. Rodriguez-Paredes M, Esteller M. Cancer epigenetics reaches mainstream oncology. Nat Med. 2011;17:330–9.View ArticlePubMedGoogle 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.View ArticlePubMedPubMed CentralGoogle 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.View ArticlePubMedPubMed CentralGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedPubMed CentralGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedPubMed CentralGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedPubMed CentralGoogle Scholar
  23. Hansen KD, Aryee M, Timp W. minfiData: Example data for the Illumina Methylation 450 k array. R package version 0.7, 1.Google Scholar
  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.View ArticlePubMedPubMed CentralGoogle 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.View ArticlePubMedGoogle Scholar
  26. Laird PW. Principles and challenges of genome-wide DNA methylation analysis. Nat Rev Genet. 2010;11:191–203.View ArticlePubMedGoogle 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.View ArticleGoogle 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.View ArticleGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle 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.PubMedGoogle 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.PubMedGoogle 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.View ArticlePubMedGoogle 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.PubMedGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle Scholar
  40. Wagner EF, Karsenty G. Genetic control of skeletal development. Curr Opin Genet Dev. 2001;11:527–32.View ArticlePubMedGoogle 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.View ArticlePubMedGoogle 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.View ArticlePubMedPubMed CentralGoogle Scholar
  43. Valdes AM, Spector TD. The contribution of genes to osteoarthritis. Rheum Dis Clin North Am. 2008;34:581–603.View ArticlePubMedGoogle Scholar

Copyright

© Aref-Eshghi et al. 2015

Advertisement