Transcriptional Associations of Osteoarthritis ... - Wiley Online Library

4 downloads 241 Views 467KB Size Report
Apr 14, 2015 - Transcriptional Associations of Osteoarthritis-Mediated. Loss of Epigenetic ... articular cartilage and toward subsequent development of treatments targeting .... carried out in R software version 3.0.2. Analyses were cor-.
ARTHRITIS & RHEUMATOLOGY Vol. 67, No. 8, August 2015, pp 2108–2116 DOI 10.1002/art.39162 C 2015, American College of Rheumatology V

Transcriptional Associations of Osteoarthritis-Mediated Loss of Epigenetic Control in Articular Cartilage Wouter den Hollander,1 Yolande F. M. Ramos,1 Nils Bomer,1 Stefan Elzinga,1 Ruud van der Breggen,1 Nico Lakenberg,1 Wesley J. de Dijcker,1 H. Eka D. Suchiman,1 Bouke J. Duijnisveld,1 Jeanine J. Houwing-Duistermaat,1 P. Eline Slagboom,2 Steffan D. Bos,2 Rob G. H. H. Nelissen,1 and Ingrid Meulenbelt2 Objective. To identify osteoarthritis (OA) progression– modulating pathways in articular cartilage and their respective regulatory epigenetic and genetic determinants in end-stage disease. Methods. Transcriptional activity of CpG was assessed using gene expression data and DNA methylation data for preserved and lesional articular cartilage samples. Disease-responsive transcriptionally active CpG were identified by means of differential methylation between preserved and lesional cartilage. Transcriptionally relevant genetic determinants were addressed by means of single-nucleotide polymorphisms (SNPs) proximal to the OA-responsive transcriptionally active CpG. Statistical analyses were corrected for age, sex, joint, and technical

covariates. A random effect was included to correct for possible correlations between paired samples. Results. Of 9,838 transcribed genes in articular cartilage, 2,324 correlated with the methylation status of 3,748 transcriptionally active CpG; both negative (n 5 1,741) and positive (n 5 2,007) correlations were observed. Hypomethylation and hypermethylation (false discovery rate of 0.05) were observed for 62 and 25 transcriptionally active CpG, respectively, covering 70 unique genes. Enrichment for developmental and extracellular matrix maintenance pathways indicated possible reactivation of endochondral ossification. Finally, we observed 31 and 26 genes for which methylation and expression, respectively, were additionally affected by genetic variation. Conclusion. We identified tissue-specific genes involved in OA disease progression, reflected by genetic and pathologic epigenetic regulation of transcription, primarily at genes involved in development. Therefore, transcriptionally active SNPs near these genes may serve as putative susceptibility alleles. Our results constitute an important step toward understanding the reported widespread epigenetic changes occurring in OA articular cartilage and toward subsequent development of treatments targeting disease-driving pathways.

Supported by the European Union Seventh Framework Programme (Project TreatOA and Project IDEAL; FP7/2007-2011 under grant agreements 200800 and 259679, respectively). The Research Arthritis and Articular Cartilage cohort studies were supported by Leiden University Medical Center and the Dutch Arthritis Association (DAA 101-402 and Reumafonds LRR), with additional support from the Centre for Medical Systems Biology and The Netherlands Consortium for Healthy Aging within the framework of The Netherlands Genomics Initiative. 1 Wouter den Hollander, MSc, Yolande F. M. Ramos, PhD, Nils Bomer, MSc, Stefan Elzinga, BSc, Ruud van der Breggen, BSc, Nico Lakenberg, BSc, Wesley J. de Dijcker, BSc, H. Eka D. Suchiman, MSc, Bouke J. Duijnisveld, MD, MSc, Jeanine J. HouwingDuistermaat, PhD, Rob G. H. H. Nelissen, MD, PhD: Leiden University Medical Center, Leiden, The Netherlands; 2P. Eline Slagboom, PhD, Steffan D. Bos, PhD, Ingrid Meulenbelt, PhD: Leiden University Medical Center, Leiden, The Netherlands, and The Netherlands Genomics Initiative–sponsored Netherlands Consortium for Healthy Aging, Rotterdam, The Netherlands. Dr. Nelissen receives royalties for contributions to a Dutch educational book (Leerboek Orthopedie). Address correspondence to Wouter den Hollander, MSc, Department of Medical Statistics and Bioinformatics, Section of Molecular Epidemiology, Leiden University Medical Center, PO Box 9600, 2300 RC Leiden, The Netherlands. E-mail: [email protected]. Submitted for publication November 27, 2014; accepted in revised form April 14, 2015.

Osteoarthritis (OA) is the most prevalent arthritis among the elderly (1) and is currently recognized as a disease of the whole joint (2). A well-described hallmark of OA is articular cartilage degradation (3). The single cell type present in articular cartilage is the articular chondrocyte, which is a highly specialized, maturation-arrested, nonproliferating cell. To ensure articular cartilage integrity throughout life, the articular chondrocyte needs to adapt its behavior in response to external signals, such as mechanical stress, aging, or microtraumas (4). To facilitate 2108

LOSS OF EPIGENETIC CONTROL IN OA ARTICULAR CARTILAGE

these adaptations, the chondrocyte requires phenotypic plasticity with proper dynamic control of gene expression to shift between active metabolic and maturation-arrested states. In this respect, chondrocytes in OA cartilage were shown to have lost their maturation-arrested state, to have regained their growth plate morphology, and to have started to proliferate, while degrading and calcifying the articular cartilage matrix (5,6). The chondrocyte phenotype is likely maintained through epigenetic control of gene expression such as DNA methylation, a biochemical process used by cells to adapt to environmental challenges such as age or disease by dynamic control of gene expression (7,8). In this respect, methylome-wide studies of articular cartilage in OA have revealed numerous loci differentially methylated between healthy and diseased tissue, while only a small minority of these loci were subsequently studied in terms of gene expression differences (9–12). Therefore, up until now it has remained unclear to what extent the large number of differentially methylated CpG in OA confer relevant gene expression changes in articular cartilage. Moreover, a growing body of literature describes how aberrant gene expression is influenced (in addition to being influenced by DNA methylation) by risk alleles in complex genetic diseases, a mechanism outlined previously in OA (13–17). These reports underscore the need to combine multiple levels of genome-wide data to gain a more robust understanding of the transcriptional processes that occur with complex genetic diseases, such as OA. In a previous report we described functional DNA methylation differences between knee and hip articular cartilage, independent of OA pathophysiology (10). Although the entire epigenomic profile of knee and hip articular cartilage is primarily defined by differentially methylated regions between the 2 joints, reports in the literature suggest that there are highly gene-specific DNA methylation changes in association with OA onset and progression (9,11,12,14–17). Therefore, in the current study we set out to identify gene-specific DNA methylation differences, independent of the joint, between preserved and lesional cartilage in patients undergoing total joint replacement surgery due to primary end-stage OA. Moreover, we combined DNA methylation changes with a previously assessed gene expression data set of overlapping samples (18) to assess OA-related changes in the epigenetically regulated transcriptome. Finally, by integrating the results with genomewide single-nucleotide polymorphism (SNP) data, we aimed to identify OA-relevant, tissue-specific genetic variants that influence gene expression in articular cartilage. The applied consecutive stepwise approach would identify novel OA susceptibility genes as well as their respective

2109

transcriptional determinants. To our knowledge, this is the first study to combine genetic, epigenomic, and transcriptomic data comprehensively to gain a functional understanding of joint-independent DNA methylation changes in relation to OA pathophysiology. MATERIALS AND METHODS The Research Arthritis and Articular Cartilage (RAAK) cohort. Ethics approval was obtained from the medical ethics committee of Leiden University Medical Center (medical ethics approval no. P08.239), and informed consent was obtained from all participants. Participant details are listed in Supplementary Table 1, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.39162/abstract. (For sampling details, see refs. 10, 14, and 18.) Methylation data. Methylation data were obtained and processed as described previously (10). Briefly, DNA was isolated using a Promega Wizard Genomic DNA Purification kit according to the manufacturer’s protocol. Next, DNA was treated with bisulfite using a Zymo Research EZ DNA Methylation kit. DNA methylation was assessed using Illumina Infinium HumanMethylation450 BeadChips. Samples were randomly dispersed, while we made sure that sample pairs were on the same chip. Using the minfi and lumi R packages (www.r-project. org), the methylation data set was filtered for probes that contained SNPs or mapped ambiguously to the genome, and color channels were separately quantile-normalized. Validation and replication using the EpiTyper platform were performed as previously reported (14). Primer sequences are listed in Supplementary Table 2, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.39162/abstract. Expression data. Normalized expression data from the RAAK study were processed and normalized as described previously (GEO accession no. GSE57218) (18). Quantitative reverse transcription–polymerase chain reaction validation primers are listed in Supplementary Table 2, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley. com/doi/10.1002/art.39162/abstract. Genotype data. Using Illumina HumanOmniExpressExome chips, genome-wide genotyping data were constructed for 216 samples from the RAAK study. SNPs with a call rate of ,95%, Hardy-Weinberg equilibrium of ,1024, minor allele frequency of ,0.01, or located on the sex chromosomes were removed prior to imputation together with Leiden Longevity Study data against the 1000 Genomes version 3 March 2012 reference panel (http://csg.sph.umich.edu/abecasis/MACH/ download/1000G.2012-03-14.html) (19). Next, SNPs that were homozygous in all the overlapping samples with the methylation (n 5 23) and expression (n 5 24) data sets were removed prior to analyses, as were SNPs that did not meet an imputation quality of 0.4 (20,21). Statistical analysis. All statistical procedures were carried out in R software version 3.0.2. Analyses were corrected for technical covariates as well as sex, joint, and age. To correct for putative correlations between preserved and OA articular cartilage from the same joint, a random effect for patient identifier (ID) was included using the lme4 package (1|ID) (22) (Figure 1). Correction for multiple testing was done for each gene separately (i), depending on the number of CpG (j) and/or SNPs (k) mapped to it (Figure 1). Multiple

2110

DEN HOLLANDER ET AL

Figure 1. Overview of the applied analysis strategy. The random effect for patient ID included in each analysis step is indicated by (1|ID). Integration of the different data levels is done by running linear mixed models for each CpG j and/or each single-nucleotide polymorphism (SNP) k within 10 kb of each gene i. t-CpG 5 transcriptionally active CpG; OA 5 osteoarthritis.

testing in the differential methylation analysis was corrected genome-wide using the Benjamini-Holm method. Methylation measurements are reported as beta values (10,23). CpG were considered differentially methylated when the mean paired difference was at least |Db| . 0.05, as smaller differences would be hard to address statistically and/or interpret biologically. Pathway enrichment was performed using the online annotation tools DAVID (http://david.abcc.ncifcrf.gov/) and STRING-DB (http://string-db. org/). A full analysis summary scheme is shown in Figure 1.

RESULTS Transcriptionally active CpG dinucleotides in articular cartilage. We recently assessed the late-stage transcriptomic profile of articular cartilage from patients who underwent total joint replacement surgery due to primary OA (GEO accession no. GSE57218) (18) (see Supplementary Table 1, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/ doi/10.1002/art.39162/abstract). For the 13,227 probes (covering 9,838 unique genes) (Figure 1) that were expressed to a detectable extent in articular cartilage, we set out to explore whether they associated with DNA

methylation of proximal CpG. To identify CpG relevant to articular cartilage in terms of transcriptional association, we correlated DNA methylation data for CpG within 10 kb of annotated genes with respective gene expression data for 13 pairs of samples (from 4 knees and 9 hips) of preserved and lesional articular cartilage. After multiple testing correction for the number of CpG for each individual gene, we observed 3,748 CpG that correlated significantly with proximal gene expression, covering a total of 2,324 unique genes (24%) (Figure 1) (see Supplementary Table 3, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art. 39162/abstract), hereafter termed transcriptionally active CpG. Notably, negative correlations (e.g., SPINT2, CILP, BFSP1, TMEM140) (Figures 2A–D) and positive correlations (e.g., COL1A2, THBS2, MSX1, RUNX3) (Figures 2E–H) were observed for 1,741 and 2,007 transcriptionally active CpG, respectively. OA-associated DNA methylation changes at transcriptionally active CpG. Next, we determined in a total of 31 pairs of samples (from 17 knees and 14 hips) which of the detected 3,748 transcriptionally active CpG

LOSS OF EPIGENETIC CONTROL IN OA ARTICULAR CARTILAGE

2111

Figure 2. Examples of cartilage-expressed genes epigenetically regulated by DNA methylation. Shown are preserved cartilage samples (blue), lesional cartilage samples (red), knee cartilage samples (circles), and hip cartilage samples (triangles). A–D, Examples of genes down-regulated upon increased DNA methylation. E–H, Examples of genes whose expression is positively correlated with increased DNA methylation. Values are the model estimate (dashed line) and 95% confidence interval. cg 5 Illumina CpG identifier.

were sensitive to the ongoing OA disease process, as reflected by differential methylation between paired samples of preserved and lesional cartilage. We observed a total of 5,282 differentially methylated CpG (false discovery rate of ,0.05, |Db|. 0.05) (Figure 3A) (see Supplementary Table 4, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/ doi/10.1002/art.39162/abstract), of which 2,188 were hypermethylated and 3,094 were hypomethylated. Among these OA-associated CpG, 151 overlapped with the observed transcriptionally active CpG, covering a total of 117 unique genes (Figure 1). Hypermethylation was observed in 59 OA-responsive transcriptionally active CpG covering 46 genes, while hypomethylation was seen in 92 OA-responsive transcriptionally active CpG covering 75 genes. Among these were genes known to be involved in OA pathophysiology (e.g., FOXA2, RUNX1, COL6A3, CD44) (Figures 3B–E) as well as multiple genes not previously reported (e.g., UACA, DLX5, DYSF, IGFBP7) (Figures 3F–I). Next, to focus solely on genes whose expression is involved in OA progression, we selected transcriptionally active CpG–regulated genes whose expression was also significantly different between preserved and lesional tissue. As such, we continued with 25 hypermethylated and 62 hypomethylated transcriptionally active CpG covering 70 unique genes (Figure 1). Subsequent gene enrichment analysis revealed significant enrichment among the 70 genes for pathways previously re-

ported to be implicated in OA pathophysiology, such as extracellular matrix (ECM) maintenance and developmental processes (see Supplementary Table 5 and Supplementary Figure 1, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/ doi/10.1002/art.39162/abstract). Technical validation and biologic replication of OA-responsive transcriptionally active CpG dinucleotides. Using the EpiTyper platform, a commonly used technique for measuring DNA methylation (14), we set out to technically validate 8 CpG in 17 pairs of preserved and lesional samples. We found a high degree of similarity between the 2 techniques (EpiTyper and HumanMethylation450 BeadChips), reflected by large Pearson’s correlation coefficients (mean r . 0.85) (see Supplementary Figure 2, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.39162/ abstract). Next, we addressed the previously observed relationship between the replicated CpG and their respective gene expression. Except for IGFBP7, we were able to validate the transcriptional involvement of all selected transcriptionally active CpG and/or the disease-associated dysregulation of their respective genes (see Supplementary Figure 3, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art. 39162/abstract). For biologic validation, DNA methylation of the selected CpG was measured in an additional 31 pairs of preserved and lesional cartilage samples. All CpG showed highly similar, significant DNA methylation

2112

DEN HOLLANDER ET AL

Figure 3. A, Volcano plot showing cutoffs (vertical and horizontal dashed lines) used to identify all CpG differentially methylated between preserved and lesional cartilage. Significant (false discovery rate of ,0.05) differentially methylated (|Db|. 0.05) CpG are depicted as green dots. B–E, Significant differential methylation between preserved and lesional cartilage in known osteoarthritis (OA)–associated genes. F–I, Significant differential methylation between preserved and lesional cartilage in genes not previously implicated in OA. Each preserved sample is set at 0 and depicted in blue, while its paired lesional sample is depicted in red. Knee joint samples are shown as circles; hip joint samples are shown as triangles. cg 5 Illumina CpG identifier.

changes as were seen in both the discovery and validation samples (see Supplementary Figure 3). Influence of genetic factors on transcriptionally active CpG methylation and expression in articular cartilage. Finally, we investigated the stable regulatory genetic environment, as reflected by transcriptionally active SNPs in proximity to the 70 genes. The presence of such SNPs may causally affect cartilage homeostasis of epigenetically controlled genes and potentially confer susceptibility to OA. The genotypes of all SNPs (dbSNP build 138) 10 kb upstream and downstream of the 70 genes were assessed in 23 sample pairs of the methylation data set. Using multivariate analysis with methylation as the dependent variable, we identified 36 OA-responsive transcriptionally active CpG that were significantly affected by at least 1 SNP (Figure 1) (see Supplementary Table 6, available on the Arthritis & Rheumatology web site at http:// onlinelibrary.wiley.com/doi/10.1002/art.39162/abstract), covering 31 unique genes. In parallel, using multivariate analysis with expression as the dependent variable, we explored whether epigenetic regulation of the 70 genes was additionally affected by the alleles of proximal SNPs. As such, we observed 26 genes whose expression was modulated by the local genetic background in conjunction with 28 transcriptionally active CpG (Figure 1) (see Supplementary Table 7, available on the Arthritis & Rheumatology web site at http://onlinelibrary.wiley.com/doi/10.1002/art.39162/ abstract). For the genes ESR, NAV2, and WLS, we observed transcriptionally active CpG and SNPs that modulated gene expression jointly.

As an example, 3 notable genes that we have observed to have transcriptomic, epigenetic, and genetic involvement in OA progression are VIT (Figures 4A–D), ROR2 (Figures 4E–H), and WLS (Figures 4I–M). All genes were differentially expressed between preserved and lesional cartilage (Figures 4A, E, and I), which was modulated by differential DNA methylation (Figures 4B, F, and J), also reflected by significant differential methylation at the respective transcriptionally active CpG (Figures 4C, G, and K). Moreover, rs11884419 and rs13292198 influenced gene expression and transcriptionally active CpG methylation of VIT and ROR2, respectively (Figures 4D and H). Additionally, rs12028757 jointly affected transcriptionally active CpG methylation and WLS expression (Figures 4L and M). DISCUSSION The current study encompasses the first comprehensive multilevel integration of genome-wide data to gain a more accurate understanding of OA-associated changes in the epigenetically regulated transcriptome of articular chondrocytes. By the stepwise integration of transcriptomic and epigenetic data in relation to cartilage OA severity, we identified 70 unique genes with OA-responsive transcriptionally active CpG likely affecting expression in articular cartilage. Subsequent pathway analyses showed significant enrichment for genes that act within skeletal development. Moreover, we have shown that for 31 and 26 OA cartilage–relevant genes, methylation and

LOSS OF EPIGENETIC CONTROL IN OA ARTICULAR CARTILAGE

2113

Figure 4. Examples of cumulative evidence for putative causal involvement of VIT, ROR2, and WLS in osteoarthritis pathophysiology. Shown are preserved cartilage samples (blue), lesional cartilage samples (red), knee cartilage samples (circles), and hip cartilage samples (triangles). In paired samples, each preserved sample is set at 0. A, E, and I, Significant differential methylation was observed between preserved and lesional cartilage for all 3 genes. B, F, and J, A significant direct relationship between expression and respective CpG methylation was observed. C, G, and K, As expected, a significant difference in methylation was observed between preserved and lesional tissue as well. D, H, L, and M, In conjunction with DNA methylation, expression was regulated by proximal single-nucleotide polymorphisms. In B, D, F, H, J, L, and M, values are the model estimate (dashed line) and 95% confidence interval. FC 5 fold change; cg 5 Illumina CpG identifier.

2114

expression, respectively, are additionally affected by genetic variation proximal to these genes. Although the observed enrichment of OA-responsive transcriptionally active CpG among genes within developmental pathways marks either disease progression or an adaptation of the preserved cartilage to the adjacent lesional tissue, our data show that changes in epigenetically regulated control of developmental genes and OA progression are markedly linked. This could indicate specific, dynamic regulation of expression of the genes in these pathways by challenged articular chondrocytes in an attempt to cope with end-stage OA. Alternatively, it could indicate that chondrocytes in end-stage disease have lost their ability to epigenetically control expression of genes essential to skeletal development, consequently regaining their growth plate morphology and starting to proliferate, with resulting debilitation of cartilage, a well-described hallmark of OA. The latter hypothesis is supported by the observed difference of epigenetic control of genes associated with skeletal development in OA, such as VIT, ROR2, and WLS. Also important in this respect, comprehensive genome-wide searches for genetic variants conferring risk of OA have resulted in robust genome-wide significant signals at genes implicated in these developmental pathways (24). In the current study, we present genetic and epigenetic loci that are functionally relevant for OAresponsive transcriptionally active CpG and cartilageexpressed genes. Ideally, these should now be followed up as candidate genes in large genome-wide association data sets to investigate whether these variants indeed confer a relatively large number of small effects that are responsible for the missing heritability observed in OA. We observed a relatively small number of differentially expressed genes that are regulated by transcriptionally active CpG and/or SNPs. Although this could be due to statistical power and/or small effect sizes, it unquestionably highlights the importance of combining epigenomic data (or gene-specific epigenetic data for that matter) with other types of molecular data to gain a robust understanding of and to biologically interpret the observed differences. Furthermore, we did not observe established genetic OA susceptibility or OA-related epigenetic loci, which implies that our transcriptional, tissue-relevant approach offers additional, compelling knowledge about genes involved in mature articular cartilage homeostasis and late OA compared to traditional genome-wide association approaches. Furthermore, the epigenetic and transcriptional effects of OA susceptibility genes such as GDF5 (15) and DIO2 (14) are relatively subtle, whereas we aimed to select genes with larger effects in late OA. Also, we may have missed

DEN HOLLANDER ET AL

long-distance or trans-acting transcriptionally active CpG or SNPs as a result of our applied 10-kb cutoff. While we do not disregard the possible impact of long-distance, transcriptionally relevant loci, the measure of effect will likely be correlated inversely with the genomic distance. In order to study these effects accurately, sample sizes larger than those in our study would be required. Unfortunately, we were unable to address the earlier observed gene-specific, transcriptionally relevant and OA-associated differentially methylated CpG, as the applied methylation array lacks the density and consequently does not measure these. VIT, ROR2, and WLS are notable examples for which we have presented functional epigenetic, genetic, and transcriptional differences (Figure 4) depending on the late pathophysiologic state of articular cartilage (25,26). VIT is a relatively understudied gene in both cartilage biology and OA. Nonetheless, proteomic analysis of mouse hip cartilage revealed its involvement in cartilage development (27). More specifically, vitrin, the protein product of VIT, contains a von Willebrand factor A domain and is involved in ECM integrin signaling (28). Expression of ROR2 drives chondrocyte expansion (29) and is known to be involved in regulating the tumor necrosis factor receptor superfamily member 11B:tumor necrosis factor superfamily member 11 protein ratio (commonly referred to as the osteoprotegerin [OPG]:RANKL ratio) in articular chondrocytes (30), a well-described disrupted pathway in OA (31). Down-regulation of ROR2 inhibits the chondrocyte’s regenerative capacities, while disruption of the OPG:RANKL ratio has been shown to induce calcification and bone formation (29,31,32). Another major player in joint development and cartilage biology is the Wnt pathway (33–36), in which ROR2 (37,38) and WLS (39,40) as well as a number of OA susceptibility genes (34,41,42) are situated. While the role of Wnt signaling is evident in cartilage development and OA, WLS is specifically involved in the endochondral ossification process (39). While OA-related differences in methylation in articular cartilage have also been reported by others (9,11,12), our results imply that changes in epigenetic control lead to expression differences only at a limited number of genes, more specifically, at genes involved either in maintaining the chondrocyte phenotype or in adversely pursuing the endochondral ossification lineage. Moreover, the detected local SNPs that affected either methylation or expression of epigenetically controlled genes in articular cartilage may inherently affect proper cartilage homeostasis and potentially affect OA susceptibility. In this regard, we observed SNPs that influenced DNA methylation at transcriptionally active

LOSS OF EPIGENETIC CONTROL IN OA ARTICULAR CARTILAGE

CpG, while no direct relationship was observed between the respective genotypes and gene expression. A large number of factors likely obscure the direct regulatory mechanisms among the local genetic background, transcriptionally active CpG methylation, and gene expression. Of note, these mechanisms likely arbitrate differential expression among the genes in which we did observe differential transcriptionally active CpG methylation but no difference in expression. SNPs and transcriptionally active CpG that appear to solely affect transcriptionally active CpG methylation or expression, respectively, are still relevant; however, the transcriptional effects of these variants should be addressed in larger consortia. Further mechanistic studies, such as performing longitudinal measurements in animal experiments or actively perturbing the relevant genes in cell systems, are needed to accurately address the hypothesis. In conclusion, we have shown that OA-related epigenetic differences need to be integrated with other sources (e.g., genomic and transcriptomic) of molecular data to enhance our understanding of the pathophysiologic processes of OA. Furthermore, by integrating multiple layers of genome-wide data, we have identified genes such as VIT, ROR2, and WLS that are likely modulating OA pathophysiology and possibly reflect the loss of the chondrocyte’s maturation-arrested state. Although targeting DNA methylation seems to be an unlikely basis for developing treatments, it serves to deepen our understanding of the complex transcriptomic changes in OA articular cartilage.

2115

4. 5. 6.

7. 8.

9.

10.

11.

12.

13.

14.

ACKNOWLEDGMENTS We thank all participants in the RAAK study.

15.

AUTHOR CONTRIBUTIONS All authors were involved in drafting the article or revising it critically for important intellectual content, and all authors approved the final version to be published. Mr. den Hollander had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis. Study conception and design. Den Hollander, Houwing-Duistermaat, Slagboom, Bos, Meulenbelt. Acquisition of data. Den Hollander, Ramos, Bomer, Elzinga, van der Breggen, Lakenberg, de Dijcker, Suchiman, Duijnisveld, Bos, Nelissen. Analysis and interpretation of data. Den Hollander, Ramos, Bomer, Houwing-Duistermaat, Meulenbelt.

REFERENCES 1. Bijlsma JW, Berenbaum F, Lafeber FP. Osteoarthritis: an update with relevance for clinical practice. Lancet 2011;377:2115–26. 2. Loeser RF, Goldring SR, Scanzello CR, Goldring MB. Osteoarthritis: a disease of the joint as an organ [review]. Arthritis Rheum 2012;64:1697–707. 3. Chen YA, Lemire M, Choufani S, Butcher DT, Grafodatskaya D, Zanke BW, et al. Discovery of cross-reactive probes and

16.

17.

18.

19.

20.

polymorphic CpGs in the Illumina Infinium HumanMethylation450 microarray. Epigenetics 2013;8:203–9. Dvir-Ginzberg M, Reich E. Chopping off the chondrocyte proteome. Biomarkers 2014;1–7. E-pub ahead of print. Sun MM, Beier F. Chondrocyte hypertrophy in skeletal development, growth, and disease. Birth Defects Res C Embryo Today 2014;102:74–82. Hosaka Y, Saito T, Sugita S, Hikata T, Kobayashi H, Fukai A, et al. Notch signaling in chondrocytes modulates endochondral ossification and osteoarthritis development. Proc Natl Acad Sci U S A 2013;110:1875–80. Laurent L, Wong E, Li G, Huynh T, Tsirigos A, Ong CT, et al. Dynamic changes in the human methylome during differentiation. Genome Res 2010;20:320–31. Slieker RC, Bos SD, Goeman JJ, Bovee JV, Talens RP, van der Breggen R, et al. Identification and systematic annotation of tissue-specific differentially methylated regions using the Illumina 450k array. Epigenetics Chromatin 2013;6:26. 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. 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. 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. Fernandez-Tajes J, Soto-Hermida A, Vazquez-Mosquera ME, Cortes-Pereira E, Mosquera A, Fernandez-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. Bell JT, Pai AA, Pickrell JK, Gaffney DJ, Pique-Regi R, Degner JF, et al. DNA methylation patterns associate with genetic and gene expression variation in HapMap cell lines [published erratum appears in Genome Biol 2011;12:405]. Genome Biol 2011;12:R10. Bomer N, den Hollander W, Ramos YF, Bos SD, van der Breggen R, Lakenberg N, et al. Underlying molecular mechanisms of DIO2 susceptibility in symptomatic osteoarthritis. Ann Rheum Dis 2014. E-pub ahead of print. Reynard LN, Bui C, Syddall CM, Loughlin J. CpG methylation regulates allelic expression of GDF5 by modulating binding of SP1 and SP3 repressor proteins to the osteoarthritis susceptibility SNP rs143383. Hum Genet 2014;133:1059–73. Hashimoto K, Oreffo RO, Gibson MB, Goldring MB, Roach HI. DNA demethylation at specific CpG sites in the IL1B promoter in response to inflammatory cytokines in human articular chondrocytes. Arthritis Rheum 2009;60:3303–13. Hashimoto K, Otero M, Imagawa K, de Andres MC, Coico JM, Roach HI, et al. Regulated transcription of human matrix metalloproteinase 13 (MMP13) and interleukin-1b (IL1B) genes in chondrocytes depends on methylation of specific proximal promoter CpG sites. J Biol Chem 2013;288:10061–72. Ramos YF, den Hollander W, Bovee JV, 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;9:e103056. Schoenmaker M, de Craen AJ, de Meijer PH, Beekman M, Blauw GJ, Slagboom PE, et al. Evidence of genetic enrichment for exceptional survival using a family approach: the Leiden Longevity Study. Eur J Hum Genet 2006;14:79–84. Ramos YF, Metrustry S, Arden N, Bay-Jensen AC, Beekman M, de Craen AJ, et al. Meta-analysis identifies loci affecting levels of the potential osteoarthritis biomarkers sCOMP and uCTX-II with genome wide significance. J Med Genet 2014;51: 596–604.

2116

21. Wood AR, Esko T, Yang J, Vedantam S, Pers TH, Gustafsson S, et al. Defining the role of common variation in the genomic and biological architecture of adult human height. Nat Genet 2014;46:1173–86. 22. Douglas Bates D, Martin Maechler M, Ben Bolker B. lme4: linear mixed-effects models using S4 classes. R package version 0.9999l/r1692. 2012. URL: http://cran.r-project.org/web/packages/lme4/ index.html. 23. arcOGEN Consortium, arcOGEN Collaborators. Identification of new susceptibility loci for osteoarthritis (arcOGEN): a genome-wide association study. Lancet 2012;380:815–23. 24. Bos SD, Slagboom PE, Meulenbelt I. New insights into osteoarthritis: early developmental features of an ageing-related disease. Curr Opin Rheumatol 2008;20:553–9. 25. Pullig O, Weseloh G, Swoboda B. Expression of type VI collagen in normal and osteoarthritic human cartilage. Osteoarthritis Cartilage 1999;7:191–202. 26. Soder S, Hambach L, Lissner R, Kirchner T, Aigner T. Ultrastructural localization of type VI collagen in normal adult and osteoarthritic human articular cartilage. Osteoarthritis Cartilage 2002;10:464–70. 27. Wilson R, Norris EL, Brachvogel B, Angelucci C, Zivkovic S, Gordon L, et al. Changes in the chondrocyte and extracellular matrix proteome during post-natal mouse cartilage development. Mol Cell Proteomics 2012;11:M111.014159. 28. Manabe R, Tsutsui K, Yamada T, Kimura M, Nakano I, Shimono C, et al. Transcriptome-based systematic identification of extracellular matrix proteins. Proc Natl Acad Sci U S A 2008; 105:12849–54. 29. DeChiara TM, Kimble RB, Poueymirou WT, Rojas J, Masiakowski P, Valenzuela DM, et al. Ror2, encoding a receptor-like tyrosine kinase, is required for cartilage and growth plate development. Nat Genet 2000;24:271–4. 30. Ermakov S, Trofimov S, Malkin I, Livshits G. A significant association exists between receptor tyrosine kinase-like orphan receptor 2 gene variants and the OPG/RANKL ratio in human plasma. Osteoporos Int 2012;23:1899–907. 31. Ramos YF, Bos SD, van der Breggen R, Kloppenburg M, Ye K, Lameijer EW, et al. A gain of function mutation in TNFRSF11B encoding osteoprotegerin causes osteoarthritis with chondrocalcinosis. Ann Rheum Dis 2014. E-pub ahead of print.

DEN HOLLANDER ET AL

32. Liu Y, Bhat RA, Seestaller-Wehr LM, Fukayama S, Mangine A, Moran RA, et al. The orphan receptor tyrosine kinase Ror2 promotes osteoblast differentiation and enhances ex vivo bone formation. Mol Endocrinol 2007;21:376–87. 33. Kruger C, Kappen C. Expression of cartilage developmental genes in Hoxc8- and Hoxd4-transgenic mice. PLoS One 2010;5: e8978. 34. Castano Betancourt MC, Cailotto F, Kerkhof HJ, Cornelis FM, Doherty SA, Hart DJ, et al. Genome-wide association and functional studies identify the DOT1L gene to be involved in cartilage thickness and hip osteoarthritis. Proc Natl Acad Sci U S A 2012;109:8218–23. 35. Leijten JC, Emons J, Sticht C, van Gool S, Decker E, Uitterlinden A, et al. Gremlin 1, frizzled-related protein, and Dkk-1 are key regulators of human articular cartilage homeostasis. Arthritis Rheum 2012;64:3302–12. 36. Lories RJ, Corr M, Lane NE. To Wnt or not to Wnt: the bone and joint health dilemma. Nat Rev Rheumatol 2013;9:328–39. 37. Ford CE, Qian Ma SS, Quadir A, Ward RL. The dual role of the novel Wnt receptor tyrosine kinase, ROR2, in human carcinogenesis. Int J Cancer 2013;133:779–87. 38. Fukuyo S, Yamaoka K, Sonomoto K, Oshita K, Okada Y, Saito K, et al. IL-6-accelerated calcification by induction of ROR2 in human adipose tissue-derived mesenchymal stem cells is STAT3 dependent. Rheumatology (Oxford) 2014;53: 1282–90. 39. Lu C, Wan Y, Cao J, Zhu X, Yu J, Zhou R, et al. Wnt-mediated reciprocal regulation between cartilage and bone development during endochondral ossification. Bone 2013;53:566–74. 40. Wang LT, Wang SJ, Hsu SH. Functional characterization of mammalian Wntless homolog in mammalian system. Kaohsiung J Med Sci 2012;28:355–61. 41. Evangelou E, Chapman K, Meulenbelt I, Karassa FB, Loughlin J, Carr A, et al. Large-scale analysis of association between GDF5 and FRZB variants and osteoarthritis of the hip, knee, and hand. Arthritis Rheum 2009;60:1710–21. 42. Loughlin J, Dowling B, Chapman K, Marcelline L, Mustafa Z, Southam L, et al. Functional variants within the secreted frizzled-related protein 3 gene are associated with hip osteoarthritis in females. Proc Natl Acad Sci U S A 2004;101: 9757–62.