Widespread epigenomic, transcriptomic and proteomic differences between hip osteophytic and articular chondrocytes in osteoarthritis

Abstract Objectives To identify molecular differences between chondrocytes from osteophytic and articular cartilage tissue from OA patients. Methods We investigated genes and pathways by combining genome-wide DNA methylation, RNA sequencing and quantitative proteomics in isolated primary chondrocytes from the cartilaginous layer of osteophytes and matched areas of low- and high-grade articular cartilage across nine patients with OA undergoing hip replacement surgery. Results Chondrocytes from osteophytic cartilage showed widespread differences to low-grade articular cartilage chondrocytes. These differences were similar to, but more pronounced than, differences between chondrocytes from osteophytic and high-grade articular cartilage, and more pronounced than differences between high- and low-grade articular cartilage. We identified 56 genes with significant differences between osteophytic chondrocytes and low-grade articular cartilage chondrocytes on all three omics levels. Several of these genes have known roles in OA, including ALDH1A2 and cartilage oligomeric matrix protein, which have functional genetic variants associated with OA from genome-wide association studies. An integrative gene ontology enrichment analysis showed that differences between osteophytic and low-grade articular cartilage chondrocytes are associated with extracellular matrix organization, skeletal system development, platelet aggregation and regulation of ERK1 and ERK2 cascade. Conclusion We present a first comprehensive view of the molecular landscape of chondrocytes from osteophytic cartilage as compared with articular cartilage chondrocytes from the same joints in OA. We found robust changes at genes relevant to chondrocyte function, providing insight into biological processes involved in osteophyte development and thus OA progression.


Introduction
OA is a degenerative joint disease characterized clinically by pain and loss of physical function [1]. It is very common, affecting >40% of individuals over the age of 70 years [2]. There is no curative therapy; end-stage disease is treated by joint replacement surgery. This procedure provides an opportunity to directly examine and characterize disease tissue from patients using genomic technologies.
A key feature of OA is cartilage degeneration, and several genomics studies have investigated the molecular characteristics of this process (e.g. reviewed in [3]). However, another important feature that can develop in joints affected by OA is the osteophyte, an area of apparent new tissue formation consisting of a cartilage-topped bony outgrowth. Osteophytes can have a significant clinical impact on both pain and loss of movement and are a typical radiographic feature of OA [4]. Osteophytes arise primarily on the margins of the articular cartilage from cells of the periosteum or synovium by a process of endochondral ossification within newly forming fibrocartilage [5]. While the pathogenesis of osteophytes has been studied in mice [5], there is still much to be learned from the molecular characterization of these important structures, particularly in human joints.
A previous study of osteophytic cartilage in the knee investigated gene expression using a microarray and suggested large differences compared with 'macroscopically intact' articular cartilage [6]. The authors hypothesized that cells in the osteophyte transition between chondrocyte and hypertrophic chondrocyte phenotypes, as opposed to a stable chondrocyte phenotype in articular cartilage.
In this proof-of-concept study, we report the first analysis of hip OA patient tissue using integrated multi-omics across genome-wide DNA methylation, RNA sequencing and quantitative proteomics to obtain a molecular portrait of chondrocytes from the cartilaginous layer of osteophytes. We identify key molecular players linked to this aspect of the disease and its pathogenesis.

Methods
More details, including references, are given in the supplementary Methods, available at Rheumatology online.

Patients and samples
The study complies with the Declaration of Helsinki. Tissue samples were collected under National Research Ethics approval reference 11/EE/0011, Cambridge Biomedical Research Centre Human Research Tissue Bank, Cambridge University Hospitals, UK. Three samples each were collected from nine patients (six women, three men, age 4484 years) undergoing hip joint replacement surgery for OA. All patients provided written informed consent before participation. Cartilage tissue was classified macroscopically for each femoral head as: low-grade, with a smooth surface and no obvious evidence of damage or fibrillation; highgrade, with damaged and fibrillated cartilage; and osteophytic, from the cartilaginous layer of osteophytes located mainly around the margins of the articular surface (sample extraction section of the supplementary methods and supplementary Figs S1 and S2, available at Rheumatology online). In each zone, a cartilage sample was removed, with subsequent extraction of DNA, RNA and protein.
Cartilage from all structural layers was obtained, with care taken to avoid the removal of non-cartilage tissue. Details for histological examination, chondrocyte preparation, as well as extraction of DNA, RNA and protein are described in the supplementary methods, sample extraction section, available at Rheumatology online.

Proteomics
Liquid chromatography -mass spectrometry (LC-MS) analysis was performed on the Dionex Ultimate 3000 UHPLC system coupled with the Orbitrap Fusion Tribrid Mass Spectrometer (Thermo Fisher Scientific, Waltham, MA, USA). Abundance values were normalized by the sum of all protein abundances in a given sample, then log2-transformed and quantile normalized. No protein was detected in only osteophytic or only articular chondrocytes. Hence we restricted the analysis to 4653 proteins that were quantified in all individuals and tissues. Details for sample processing, LC-MS analysis, protein identification and quantification are described in the supplementary methods, proteomics section, available at Rheumatology online.

RNA sequencing
Multiplexed libraries were sequenced on Illumina HiSeq 2000 (Illumina Inc., San Diego, CA, USA; 75 bp paired-end read length) and a cram file was produced for each sample. We obtained transcript-level quantification using salmon 0.7.2 [7]. After quality control, we retained 14 029 genes. Details for sample processing and read quantification are described in the supplementary methods, RNA sequencing section, available at Rheumatology online.

Methylation
Methylation was assayed using the Illumina 450k BeadChip (Illumina Inc., San Diego, CA, USA). The resulting idat files were parsed and QCed using ChAMP [8] in R, yielding 424 705 probes. The probe beta values were normalized using the funnorm method [9] in R, and converted to M-values. Details for Illumina 450k BeadChip processing and quality control are described in the supplementary methods, methylation section, available at Rheumatology online.

Differential analysis
The proteomics, gene expression, and probe methylation differential analyses were carried out using limma [10] in R. We used a within-individual paired sample design, that is, the individual ID as a covariate in the comparison of osteophytic to low-grade, osteophytic to high-grade and high-grade to low-grade tissue. A Benjamini-Hochberg false discovery rate (FDR) was applied to each analysis to correct for multiple testing. Differentially methylated regions (DMRs) were identified using the DMRcate R package [11].

Principal component analysis
Principal component analyses were carried out using the prcomp function in R for each omics level, based on significant differences between osteophytic and low-grade cartilage at 0.1% FDR.

Gene identifier mapping
To map Ensembl gene IDs to gene names and vice versa, we used the assignment in Ensembl (downloaded from Ensembl biomart, GRCh38.p7). We only included instances where a unique Ensembl gene ID corresponded to a unique gene name.

UK Biobank association analysis
We applied MAGMA v1.06 [12] to test the joint association of genetic variants in the 56 genes changed between osteophytic and low-grade articular cartilage on all three molecular levels. We used genetic association data from UK Biobank, including 2396 hospital-diagnosed hip OA cases and 9593 non-OA controls based on ICD 10 and/ or 9 codes; controls were not diagnosed with any musculoskeletal disorders, symptoms or signs. Further details of the dataset, including quality control, are described in the supplementary methods, UK Biobank association analysis section, available at Rheumatology online.
Each gene was assigned the single nucleotide polymorphisms (SNPs) located between the gene's start and stop sites based on NCBI 37.3 gene definitions. For each gene, we used the combined statistic based on the sum of SNP log-P-values and the lowest SNP P-values, as recommended by MAGMA. Linkage disequilibrium (LD) was calculated from a subset of the UK Biobank samples (see supplementary methods, UK Biobank association analysis section, available at Rheumatology online). We then tested whether the 56 genes are more associated with OA than expected by chance, correcting for the potentially confounding effects of sample size, gene size, gene density and the inverse of the mean minor allele count in the gene, as well as the log of these variables, as recommended.
Integrative Gene Ontology gene set analysis Gene Ontology (GO) [13] biological process and molecular function gene annotations were obtained from Ensembl Biomart. We followed the integrative cross-omics analysis from [14], as described in detail in the supplementary methods, Gene Ontology gene-set analysis section, available at Rheumatology online. Briefly, enrichment of each annotation for each of the three omics levels was assessed using a one-sided hypergeometric test, the Pvalues integrated across the three omics levels followed by randomizations. Significance was defined at 5% FDR. We excluded annotations that were enriched in only one of the omics levels, or where fewer than five genes contributed to the enrichment on at least two omics levels.

Widespread differences between osteophytic and articular cartilage
First, we examined genome-wide methylation, gene expression and protein abundance differences between osteophytic and low-grade articular cartilage. At each molecular level, we found widespread differences at 0.1% FDR. In particular, we found significant differences in protein abundance for 942 of 4653 proteins (supplementary Table S1, available at Rheumatology online), in gene expression for 3601 of 14 029 genes (supplementary Table  S2, available at Rheumatology online), and 3161 DMRs that overlapped 3277 genes (supplementary Table S3, available at Rheumatology online).
Second, we also examined molecular differences between chondrocytes from osteophytic and high-grade articular cartilage. Fewer differences were significant at 0.1% FDR: protein abundance differences for 517 of 2653 genes (supplementary Table S1, available at Rheumatology online), gene expression differences for 1512 of 14 029 genes (supplementary Table S2, available at Rheumatology online) and 1113 DMRs overlapping 1113 genes (supplementary Table S4, available at Rheumatology online). However, globally, the gene expression differences between osteophytic and high-grade cartilage were similar to the differences between osteophytic and low-grade articular cartilage ( Fig. 1a and b; for estimates of the log-fold-differences in gene expression and protein abundance, the 95% CIs overlap for 97.9% of genes for gene expression and 97.6% of genes for protein abundance). Moreover, for both gene expression and protein abundance levels, of the genes significant at 0.1% FDR in a given comparison (osteophytic/low-grade or osteophytic/high-grade), 599.5% show the same direction of change in the other comparison and 590% are at least nominally significant. Similarly, of the genes contained in DMRs at 0.1% FDR in one comparison, 593.5% are contained in DMRs at 5% FDR in the other comparison.
This suggests that the differences between osteophytic and low-grade articular cartilage are similar to, but more pronounced than the differences between osteophytic and high-grade articular cartilage. In agreement with this, on each omics level, we found that a principal component analysis based on the significant differences between osteophytic and low-grade articular cartilage separated the osteophytic cartilage from both the lowgrade and the high-grade articular cartilage samples (Fig. 2). Hence we took the comparison of chondrocytes from osteophytic and low-grade articular cartilage as the basis for further analyses.
Differences between osteophytic and low-grade articular cartilage are correlated with differences between high-and low-grade articular cartilage At any given level of significance, we identified far more significant differences in the comparison of osteophytic and low-grade articular cartilage than in the comparison of high-and low-grade articular cartilage. For example, at 0.1% FDR, we found no protein with differential abundance between high-and low-grade articular cartilage (supplementary Table S1, available at Rheumatology online), only one gene with differential RNA expression (supplementary Table S2, available at Rheumatology online), and no differentially methylated regions.
The differences observed between osteophytic and low-grade cartilage are significantly correlated with the differences found between high-and low-grade articular cartilage ( Fig. 1c and d; gene expression: Spearman r = 0.62, protein abundance: r = 0.47; both P < 10 À15 ). https://academic.oup.com/rheumatology FIG. 1 Gene expression and protein abundance differences between osteophytic chondrocytes, low-grade, and highgrade articular chondrocytes (A and B) Differences between osteophytic and low-grade articular chondrocytes are correlated with differences between osteophytic and high-grade articular chondrocytes for protein abundance (A) and gene expression (B). (C and D) Differences between osteophytic and low-grade articular chondrocytes are correlated with differences between high-and low-grade articular chondrocytes for protein abundance (C) and gene expression (D). (E) Differences in gene expression and protein abundance identified between osteophytic and low-grade articular chondrocytes are correlated. Each point represents one gene. Black: genes with significant changes between osteophytic and low-grade articular cartilage at 0.1% FDR. Red: genes with significant changes on both protein and RNA level between osteophytic and low-grade articular chondrocytes at 0.1% FDR.
The direction of difference (increase or decrease) in osteophytic compared with low-grade articular chondrocytes agrees with the direction of difference in high-grade compared with low-grade articular chondrocytes for 74% of genes on mRNA level and for 66% of genes on protein level.
In agreement with this, the vast majority of genes with significant mRNA or protein level differences between osteophytic and low-grade articular cartilage also show the same direction of difference in high-grade compared with low-grade articular cartilage (gene expression: 90.1%, protein abundance: 86.5%; both P < 10 À15 ; supplementary Tables S1 and S2, available at Rheumatology online). However, the magnitude of the log-fold-differences between the high-and low-grade articular chondrocytes is smaller (gene expression: 99.5% of 3601 significant genes and 75.4% of all genes, protein abundance: 99.9% of 942 significant proteins and 73.9% of all proteins).

RNA sequencing and proteomics
There was a significant positive correlation of the differences in gene expression and protein abundance identified between osteophytic and low-grade articular cartilage in the proteomics and the RNA sequencing data ( Fig. 1e; Spearman r = 0.31, P < 10 À15 ; based on 4345 genes present in both). We identified 309 genes with significant differences on both RNA and protein level at 0.1% FDR, 88.7% of these differences were directionally concordant (binomial P < 10 À15 ).

Methylation, RNA sequencing and proteomics
We found 56 genes that showed differences in protein abundance and gene expression levels, and also overlapped a DMR between osteophytic and low-grade articular cartilage ( Table 1). The direction of change on RNA and protein level agreed for 52 of the 56 genes (binomial P < 10 À10 ). For all 56 genes, the direction of difference between osteophytic and high-grade articular chondrocytes is the same as the direction of difference between osteophytic and low-grade articular chondrocytes, and 42 genes also have evidence for difference between osteophytic and high-grade articular chondrocytes across all three molecular levels at 5% FDR or lower (Table 1).
Link to genetic variants associated with OA Of the 56 genes that showed significant differences between osteophytic and low-grade articular cartilage on all three molecular levels, two have been robustly associated with OA in published genome-wide association studies. Both also show directionally concordant, significant RNA and protein level differences between osteophytic and high-grade articular chondrocytes at 0.1% FDR, and are contained in DMRs at 5% FDR.
COMP demonstrates significantly lower gene and protein levels in osteophytic compared with low-grade articular cartilage, and is located in a hyper-methylated region. The c.1141 G > C (p.Asp369His) missense variant in the gene has been found to significantly increase the risk of OA in a study of hip OA patients who underwent joint replacement surgery [15].
ALDH1A2 also displays significantly reduced gene and protein levels in osteophytic compared with low-grade articular cartilage, and is located in a hyper-methylated region. Several genetic variants in and close to ALDH1A2 have been associated with severe hand OA [16]. The most strongly associated variant is rs12907038; the risk allele has been associated with a decrease of ALDH1A2 gene expression, and another associated variant has also been associated with allelic imbalance in ALDH1A2 gene expression [16].
We further tested the joint association of all 56 genes with susceptibility to OA using an unpublished dataset from UK Biobank (2396 hip OA cases, 9593 non-OA FIG. 2 PCA separates osteophytic from low-and high-grade articular cartilage PCA based on proteomics data (A), RNA sequencing data (B) and probe methylation data (C). Each point represents one sample. The PCA was carried out on the proteins or genes with significant differences or within DMRs between osteophytic and low-grade articular cartilage; the plots show that these expression or methylation patterns also separate osteophytic from degraded cartilage. PCA: principal component analysis. controls). There was no significant excess of genetic association in the 56 genes together (P = 0.2116 using MAGMA, see Methods section).

Cross-omics pathway analysis
In an integrative comparison of osteophytic and low-grade articular cartilage using all three molecular levels (see Methods section), we identified 36 GO annotations as significantly associated with the molecular changes at 5% FDR (supplementary Table S5, available at Rheumatology online). The most significant annotations (Table 2; all FDR < 2%) include gene sets with links to OA (reviewed e.g. in [3]), such as extracellular matrix organization and collagen catabolic process; skeletal system development; inflammatory response; positive regulation of the ERK1 and ERK2 cascade; and platelet aggregation. A further annotation with highly significant association (integrative FDR < 2%) and P < 0.05 on each one of the omics levels was endodermal cell differentiation ( Table 2).

Replication of gene expression changes
A previous study [6] examined gene expression differences between osteophytic and low-grade articular cartilage in knee OA patients using a microarray. The authors identified 31 genes with >20-fold change in gene expression and P < 0.005. Of these 31 genes, 18 are present in our RNA sequencing data, and all have directionally concordant effects (supplementary Table S6, available at Rheumatology online; binomial P < 10 À5 ). This includes TF and ALDH1A2, which were identified as significantly different between osteophytic and low-grade articular cartilage on all three omics levels in this study.

Discussion
This study provides a systematic molecular characterization of osteophytic chondrocytes in OA across genomewide methylation, gene and protein expression levels.
We have shown widespread molecular differences between chondrocytes from osteophytic and low-grade articular OA chondrocytes (as previously seen for gene expression in the knee), with similar but smaller differences between osteophytic and high-grade articular OA chondrocytes. By contrast, there were far fewer significant differences between chondrocytes from high-and low-grade articular cartilage for any given FDR. In a direct comparison, we have shown that the differences between osteophytic chondrocytes and those from lowgrade cartilage are positively correlated with the differences between chondrocytes from high-and low-grade articular cartilage. The correlation between these differences is observed despite the morphological and histological dissimilarities between osteophytic and articular tissues. One interpretation would be that, although both tissues are subject to the disease process and the altered internal joint milieu [4], osteophytic chondrocytes are better able to respond in terms of new cartilage production, which may represent attempted joint recovery (e.g. as in the transient phenotype suggested by [6]). Moreover, osteophytes principally develop at the periphery of the articular surface in response to an altered mechanical environment; chondrocyte proliferation is followed by chondrocyte hypertrophy and endochondral ossification within the osteophyte. Several of the genes and pathways identified in this study are implicated in these processes.
Of the 56 genes with differences between osteophytic and low-grade articular cartilage on all three molecular levels, gene expression changes in Transferrin (TF), MMP13 and ALDH1A2 have also been previously identified in the knee [6], indicating prominent involvement of these genes. TF (reduced mRNA and protein levels in osteophytic chondrocytes) is known to be produced in hypertrophic cartilage as a pro-angiogenic molecule [17] and is found in the synovial fluid of OA patients [18]. MMP13 (increased mRNA and protein levels in osteophytic chondrocytes) is a marker of chondrocyte hypertrophy considered important in cartilage degeneration [19,20]. ALDH1A2 is known to have genetic links to OA, discussed in detail below. The 56 genes with significant differences across omics levels also include CILP (reduced mRNA and protein levels in osteophytic chondrocytes), which was previously found to be down-regulated in mechanically induced OA in mice [21]; and IL6 (increased mRNA and protein levels in osteophytic chondrocytes), which has been localized to chondroblasts and preosteoblasts in human osteophytes during endochondral ossification [22]. All of the genes highlighted above also show evidence for differences between chondrocytes from osteophytic and high-grade articular cartilage across all three molecular levels at 5% or lower FDR. Notably, the 56 genes with differences between chondrocytes from osteophytic and low-grade articular cartilage on all three molecular levels also include two genes that harbor genetic variants robustly associated with OA identified from genome wide association studies, ALDH1A2 and COMP. ALDH1A2 (aldehyde dehydrogenase 1 family, member A2 or retinaldehyde) is an enzyme that catalyses the synthesis of retinoic acid, the active derivative of vitamin A (retinol). Vitamin A is involved in postnatal bone health and bone remodelling, with both high and low levels having negative effects [23]. COMP is a constituent of the cartilage matrix, present in the interterritorial matrix and is involved in collagen fibrillogenesis [24]. Notably, the decreased expression of both genes is consistent across all three molecular levels in this study, and consistent with the molecular mechanisms suggested for the associated genetic variants (reduction of gene expression for ALDH1A2 and a missense variant in COMP). We did not find such cross-omics significant changes of ALDH1A2 or COMP between low-and high-grade articular cartilage (with only protein-level changes ALDH1A2 significant at 5% FDR), which could suggest that their action is stronger in osteophytic cartilage, or could be due to the limited power in this discovery study. However, we have replicated gene expression changes of ALDH1A2 in osteophytic cartilage using independent data. Interestingly, the genetic association of ALDH1A2 was identified in severe hand OA [16], which is characterized by node formation.
Using UK Biobank data, we did not find evidence for association of the joint set of the 56 genes with differences between chondrocytes from osteophytic and lowgrade articular cartilage on all three molecular levels. It is possible that some of these genes exhibit molecular changes as a consequence rather than cause of the disease, or are involved in OA progression rather than incidence. The lack of association could also be due to still limited sample size of the genetic data. Larger cohorts will be required to determine the comprehensive set of genetic variants associated with OA, as well as which molecular changes are causal to disease processes.
In the cross-omics GO analysis, several of the gene annotations with highly significant associations have known links to OA: changes in the extracellular matrix, collagen catabolism, inflammation and activity of the ERK cascade are all known interrelated processes taking place in the OA joint [25]. The ERK signalling pathway is important in mesenchymal cell differentiation and can be regulated by mechanical stimuli during joint formation [26,27]. Another annotation with highly significant cross-omics association (integrative FDR < 2%) was endodermal cell differentiation, which is not directly linked to cartilage, but could reflect tissue development factors involved in osteophyte formation.
We did not find cross-omics differences between osteophytic and articular cartilage in previously reported osteophyte development genes such as TGF, PTH and IGF1 [4]. This could be explained by the small sample size of the study, or the fact that the chondrocytes were taken from patients with end-stage hip OA, so may not reflect processes involved early in osteophyte development.
The major strengths of this study are the integration of DNA methylation, gene expression and proteomics data, for a comprehensive overview of the changes across molecular levels, and the precise matching of osteophytic, low-and high-grade articular cartilage samples from the same joint. The latter reduces the possibility of false-positives due to biological differences between individuals. This approach has helped illuminate the molecular basis of OA progression; tissue from healthy individuals and early OA stages would be required to characterize the onset of the disease. The main limitation of this study is the size of the cohort examined here. As such, this study is a proof-of-concept discovery study, and replication in larger independent datasets will be required. The deposition of all data in open repositories also allows the data to be combined with other datasets in the future.
As noted above, the molecular associations identified may be a result of, rather than causal to the disease processes. Nonetheless, they can provide insights into characteristics of osteophytic chondrocytes and into disease progression, and suggest targets for further functional follow-up with translational potential.
In summary, we present the first integrative methylation, gene expression and proteomics study across osteophytic and articular cartilage in hip OA. We have identified multiple genes with significant cross-omics changes between osteophytic and low-grade articular cartilage, including two genes associated with OA through genome-wide association studies.
These findings offer evidence that the study of osteophytic cartilage can provide distinct insight into OA pathogenesis.