DNA methylation and noncoding RNA in OA: Recent findings and methodological advances.

Introduction: Osteoarthritis (OA) is a chronic musculoskeletal disease characterized by progressive loss of joint function. Historically, it has been characterized as a disease caused by mechanical trauma, so-called 'wear and tear'. Over the past two decades, it has come to be understood as a complex systemic disorder involving gene-environmental interactions. Epigenetic changes have been increasingly implicated. Recent improvements in microarray and next-generation sequencing (NGS) technologies have allowed for ever more complex evaluations of epigenetic aberrations associated with the development and progression of OA.
Methods: A systematic review was conducted in the Pubmed database. We curated studies that presented the results of DNA methylation and noncoding RNA research in human OA and OA animal models since 1985.
Results: Herein, we discuss recent findings and methodological advancements in OA epigenetics, including a discussion of DNA methylation, including microarray and NGS studies, and noncoding RNAs. Beyond cartilage, we also highlight studies in subchondral bone and peripheral blood mononuclear cells, which highlight widespread and potentially clinically important alterations in epigenetic patterns seen in OA patients. Finally, we discuss epigenetic editing approaches in the context of OA. Conclusions: Although a substantial body of literature has already been published in OA, much is still unknown. Future OA epigenetics studies will no doubt continue to broaden our understanding of underlying pathophysiology and perhaps offer novel diagnostics and/or treatments for human OA.

Introduction

Osteoarthritis (OA) is a debilitating musculoskeletal disease that is characterized by a progressive loss of joint function. Patients who suffer from OA often exhibit chronic, significant pain, mobility loss, and functional impairments involving diarthrodial joints (hips, knees, and hands). As average life expectancy rises, OA has emerged as the most rapidly growing major health condition and the leading cause of chronic disability in the United States [1,2]. Our understanding of OA in recent years has changed drastically. Once thought as simple ‘wear and tear’, it is now well understood that a variety of processes, including inappropriate joint loading [3], dysfunctional healing responses to cartilage injury and low-level chronic inflammation [4], both in the synovium [5] and systemically [6] play significant roles. The contributions of genetic variation to OA susceptibility are varied depending on the joint involved; large twin studies have estimated the heritability of severe hip OA to average 73%, whereas the heritability of severe knee OA is far less at 45% [7].

A variety of non-genetic factors have been demonstrated as key contributors to the development of OA, including age, mechanical trauma, and local inflammatory processes [3]. Among many cellular mechanisms responsible for the physiological integration of these environmental signals, epigenetics has emerged as one of the most significant, and has been linked to both development and progression of OA. Epigenetics refers to the study of any heritable changes in gene expression caused by mechanisms other than genomic DNA mutations. Through epigenetics, gene transcription patterns are altered as a response to various internal and environmental signals. The canonical epigenetic control mechanisms include cytosine genomic DNA methylation, noncoding RNA, and histone post-translational modifications [8].

DNA methylation is a process during which DNA methyltransferase (DNMT, see Table 1 for list of abbreviations) enzymes transfer a methyl group from S-adenosyl methionine (SAM, the ‘methyl donor’) onto the C5 position of cytosine (5 mC) occurring most often in CpG dinucleotides. In vertebrates to date there have been three DNMTs identified (DNMT1, DNMT3a, DNMT3b) [9]. DNMT1 is the maintenance methyltransferase, functioning to maintain existing DNA methylation patterns during DNA replication, whereas DNMT3a and DNMT3b play a role in de novo creation of DNA methylation patterns. Hypermethylation in gene regulatory regions generally suppresses gene transcription whereas hypomethylation generally enhances it. DNA methylation changes are a fluid process with changes occurring throughout an organism's lifespan. However, several months or years after birth, there is a tremendous rise in DNA methylation levels [10] where the majority of methylation patterns are ‘set’ in differentiated tissues. Active ‘demethylation’ of cytosine ‘occurs when ten-eleven translocases (TET) convert 5 mC into 5-hydroxymethylcytosine (5 hmC), 5-formylcytosine (5 fC), or 5-carboxylcytosine (5caC) intermediates which are rapidly converted to unmethylated cytosine [11-13], thereby facilitating gene expression.

Table 1

Abbreviations used in this review and commonly featured in OA epigenetics research.

Abbreviation	Definition
DMP/DML	Differentially methylated position or differentially methylated loci, refers to a statistically significant difference in DNA methylation within a target tissue or cell subset of an individual cytosine within a CpG oligodeoxynucleotide.
DMR	Differentially methylated region, refers to a statistically significant difference n DNA methylation within a target tissue or cell subset of a CpG-dense region. Various definitions of DMR sizes exist; a commonly-used criteria would define a ‘bandwidth’ of 1000 nucleotides with a standard deviation of 500 nucleotides, including a CpG density of at least 2 CpG sites (CpG ‘bandwidth’) within 10 nucleotides of each other.
NGS	Next generation sequencing, refers to a suite of technologies allowing sequencing of an entire genome at relatively low cost. The most commonly used of these high-throughput DNA sequencing methods is Illumina sequencing by synthesis, operating at an accuracy of ~99.9% with reads on the order of 50-500bp.
SNP	Single nucleotide polymorphism, refers to a germline substitution of a single nucleotide at a given position within the genome.
CpG	Cytosine-guanosine dinucleotide, the location within the genome wherein DNA methylation occurs (5′ carbon position of cytosine) within fully differentiated tissues.
DMM	Disruption of the medial meniscus, a surgical technique commonly used to induce OA-like changes within a mouse stifle (knee) joint
5-hmC	5-hydroxymethylcytosine, a ‘transition state’ between methylated cytosine and unmethylated cytosine. The conversion of 5-mC to 5-hmC is catalyzed by the TET1 enzyme.

Among the noncoding RNAs, the most commonly studied are microRNAs (miRNAs), long noncoding RNAs (lncRNA), Piwi-interacting RNAs (piRNA), small nucleolar RNAs (snoRNAs) [14,15], and circular RNAs (circRNAs) [16] MicroRNAs are the most widely studied non-coding RNAs that have the ability to bind to specific messenger RNAs inducing their cleavage or degradation that ultimately alters protein expression. These short 22 nucleotide-long sequences can alter gene expression by binding to the 3’ tail of mRNAs thus targeting them for degradation. Furthermore, miRNAs can also bind and destabilize mRNAs without inducing degradation [17]. In contrast to miRNAs, lncRNAs, as their name suggests, are long RNAs encompassing a more heterogeneous group of noncoding RNAs that are more than 200 nucleotides in length. lncRNAs perform numerous activities in the body [18] such as using RNA-protein interactions, recruitment of regulatory complexes, and acting as direct local regulators [19] of gene expression.

In this review article, we seek to systematically review the recent literature on altered epigenetic patterns associated with OA in humans and animal models, with a specific focus on DNA methylation and noncoding RNAs. As opposed to reviewing studies chronologically, we will instead group publications based on the methodology used, including array-based and next generation sequencing-based approaches to quantification of DNA methylation as well as array-based and sequencing-based noncoding RNA studies.

Brief methods: DNA methylation quantitation by sodium bisulfite treatment

Only a handful of emerging sequencing technologies have the ability to read DNA methylation patterns. To the present day, the gold standard remains sodium bisulfite conversion of genomic DNA. Treatment of genomic DNA with sodium bisulfite results in deamination of unmethylated cytosines into uracils, whereas methylated and hydroxymethylated (5 mC and 5 hmC) cytosines are left unconverted. Following subsequent PCR amplification, uracil is converted to thymine, while methylated and hydroxymethylated cytosines are amplified as cytosines. After this process, the treated DNA can be sequenced either by PCR pyrosequencing methods, microarrays, or by large parallel sequencing methods which fall into the category of ‘next generation sequencing’, discussed later in this article. Each of these technologies ultimately differentiates methylated from unmethylated cytosine based on the ratio of cytosine to thymine in the final bisulfite-converted and sequenced product, allowing quantification of DNA methylation patterns at a single nucleotide resolution [20, 21].

DNA methylation analyses in OA

Historical context

Alterations of DNA methylation patterns as a feature of OA pathogenesis has been suggested since the mid-1980s, when it was found that type I and type II collagen genes were epigenetically regulated in chick embryos [22]. The first evaluation of DNA methylation in OA specifically was performed in 2005, when several key catabolic genes were shown to be demethylated in OA, including ADAMTS4, MMP3, MMP9, and MMP13 [23]. The adipokine leptin was later added to the list of differentially methylated genes in OA [24]. More recent single-gene epigenetic studies have included the OA susceptibility gene GDF5 [25], the Wnt and BMP-Smad signaling gene sclerostin (SOST) [26], and many other key transcription factors, catabolic, and anabolic factors including SOX9, MMP13, and COL9A1 [27-29]. Candidate analysis of inflammation-related genes have also demonstrated differential methylation within cartilage, including interleukin 8 (IL8) [30], interleukin 1 (IL-1β) [31,32].

Starting in the 1990s, DNA microarray technology emerged as both a time and cost-saving solution for large-scale genetic analyses, with subsequent adaptation of this technology to epigenetic analysis. The field has been dominated by Illumina's Infinium BeadArray and BeadChip (Illumina, San Diego, CA, USA) technology since the release of their Infinium HumanMethylation27 BeadArray in 2009 [33]. This first-generation chip provided quantification of roughly 27,000 CpG sites throughout the genome (covering nearly half of known coding regions). As research in epigenetics expanded, it was discovered that the DNA methylations occurring outside of the CpG islands were also vital factors contributing to the alteration of gene expression [34]. To encompass these regions, Illumina created an expanded Infinium HumanMethylation450 BeadChip that allowed surveillance of around 470,000 CpG sites throughout the genome in the early 2010s. However, that still only covered about 2% of CpG loci in the entire human genome. The latest release from Illumina is the Infinium HumanMethylationEPIC BeadChip that covers nearly 850, 000 CpG sites of the genome along with, non-CpG differentially methylated sites, enhancers, open chromatin, transcription factor binding sites, and miRNA promoter regions. The Illumina microarrays consist of two types of probes/beads. The M bead type is designed to measure the methylated loci, while the U bead type measures unmethylated loci. Total methylation levels of all the covered loci are calculated by measuring the ratio of the fluorescent signals emitted from methylated over unmethylated sites generally expressed by a beta-value (quantitative methylation ratio that ranges from 0 to 1).

Although not as widely utilized, Agilent Technologies also produces microarrays for genetic and epigenetic analyses. Agilent Human Promoter Microarrays (Agilent Technologies, Santa Clara, CA, USA) currently come in three sizes covering different numbers of CpG sites each. All of their arrays cover roughly 21,000 human genes as represented by RefSeq with the only difference being the resolution at which the genes are analyzed. The largest in size is their SurePrint G3 Human Promoter Kit 1 × 1M that covers 966,092 distinct biological features, followed by SurePrint G3 Human promoter 2 × 400 K Kit (covers 414,043 distinct biological features), with the smallest being Human Promoter Microarray Kit 1 × 244 K (covers 243,504 distinct biological features).

Infinium HumanMethylation27 BeadChip

The development of the 27 k microarray allowed, for the first time, a relatively inexpensive and highly reproducible way to quantify genome-wide DNA methylation levels. In 2014, Fernández-Tajes et al. [35] were the first to apply this technology to determine methylome changes in OA cartilage (Overview of differentially methylated pathways from studies included in this article in Fig. 1, timeline of microarray publications Fig. 2). The study identified a total of 91 differentially methylated probes between patients with OA and healthy controls. Among those, the runt-related transcription factor-1 (RUNX1) was the most hypomethylated, while msh homeobox-1 (MSX1) was the most hypermethylated in OA patients relative to controls. Runx1 was recently noted as the gene that plays a vital role in articular cartilage maintenance by enhancing matrix production [36]. Interestingly, in supervised clustering analysis, this study also identified a subset of OA patients that had a variety of epigenetic changes of genes related to inflammatory regulations, activation motifs, and cytokine production.

Fig. 1.

Venn diagram demonstrating gene pathways and families differentially methylated among various knee and hip OA groups.

Fig. 2.

Timeline of major advances in methodologies for studying DNA methylation using array-based techniques, with article references.

Infinium HumanMethylation450 BeadChip

Following this initial report, three other papers were published in 2014 and included a wider coverage of the DNA methylation sites leveraging the 450 k microarray (Fig. 2). With over 450,000 CpG sites now covered, there was a clear demonstration that there are epigenetic differences between the hip and knee OA samples along with differences between the OA samples and healthy cadaveric controls [37-39]. One important finding that den Hollander and colleagues discovered was that four HOX clusters were differentially methylated between knee and hip articular cartilage; they went on to confirm differential gene expression among these genes within human cartilage samples as well. Rushton et al. [37] conducted a similar study using hip and knee OA human cartilage samples and found that hip OA samples are clearly distinguishable from neck-of-femur fracture control patients, and identified 5322 differentially methylated position (DMP). Those DMPs were associated with catabolic genes involved in breakdown of cartilage extracellular matrix (ECM), including ADAMTS2, ADAMTS4, ADAMTS5, ADAMTS10, ADAMTS17, MMP13, MMP16. Additionally, genes involving various cartilage homeostatic pathways were also reported as epigenetically altered including members of TGFβ signaling pathway. Importantly, they confirmed a subset of their patient samples, both hips and knees, exhibited differential methylation of immunologically related genes, similar to the Fernandez-Tajes et al. report [35].

Our group published a study using the same 450 k arrays with hip OA patient samples using a unique internally controlled design, where we compared eroded to intact regions within the same joint. We identified 550 DM sites, the majority (378) hypomethylated [38]. Of importance was that DM sites identified were very concentrated in enhancer regions (87%). Analysis of DMPs and their corresponding genes revealed FURIN (proprotein convertase [40]) as the most hypomethylated, followed by RUNX1 (previously identified as hypomethylated by Fernandes-Tajes et al.), DLX5, and HTRA1. Hypermethylated DMPs included COL11A2 (which plays a role in maintenance of the diameter and cohesive strength of cartilage matrix [41]), along with the growth factor receptor FGFR2. Furthermore, DM genes were enriched with upstream regulators including miR-128, miR-27a (MMP regulator), and miR-9 (chondroblast survival regulator).

We then expanded beyond cartilage and applied the same methodology to a DNA methylation analysis of subchondral bone, comparing eroded and intact sections within the same joint as well as comparing subchondral bone to overlying cartilage [42]. Somewhat surprisingly, we found roughly an order of magnitude more differential methylation in subchondral bone compared to cartilage, although many of the gene pathways were similar to our (and others') previous analyses, including EIF2C2 (RISC - associated), TGFB3, NFATC1, HDAC4, CD6, HOXA7 HOXA2. We again identified a differential methylation within a number of previously reported OA susceptibility genes from the Human Genome Epidemiology Network (HuGENet) [43] including FTO, LRP5, IGFBP7, COL13A1, and NCOR2.

These findings were corroborated and expanded upon by Zhang and colleagues [44]. In their examination of subchondral bone methylation patterns, there were 72 significant DMPs found between tibial plateau regions that had intermediate and no cartilage erosion. In comparison of eroded to intact regions, 397 significant DMPs were identified, while intermediate vs. intact regions had 257 DMPs. Similar to our findings, subchondral bone methylation analysis showed enrichment in the enhancer regions. GO analysis elucidated novel differentially methylated genes (DMGs) associated with OA progression such as WNT11, vital component in osteogenesis [45]. Among the DMGs, the homeobox gene family was found to be altered in all 3 comparison groups which are interestingly also altered in intermediate-stage to late-stage OA cartilage comparison [46]. An additional 11 genes showed both methylation and expression changes, however, no consistent methylation-expression pattern was identified. Matching the subchondral and cartilage methylation patterns revealed cAMP-mediated signaling, G-protein coupled receptor signaling and tRNA splicing as the top canonical pathways. Furthermore, comparison of intermediate-stage to late-stage bone and cartilage groups revealed 111 shared DMPs and 41 DMGs signifying common pathological alterations at the late stage of OA. Among the 41 DMGs, enrichments were found in morphogenesis, skeletal system development, and homeobox transcription factors centered on TGFB node all of which play various roles in metabolism of cartilage and subchondral bone.

A number of additional studies have been subsequently carried out over the past several years leveraging the Infinium HumanMethylation450 BeadChip in knee and hip OA. A recurrent theme among these studies was the finding of differential methylation within enhancer regions that play a vital role in the OA phenotype development [47]. Among these studies, Alvarez-Garcoa and collegues used the data gained from the array and validated with in vitro experiments the link between the DNA methylation and gene expression thus showing a subset of 6 transcription factors to be significantly hypermethylated and down-regulated in OA cartilage (ATOH8, MAFF, NCOR2, TBX4, ZBTB16, ZHX2) [48]; following treatment with the DNA methylation inhibitor decitabine, expression of all 6 transcription factors increased in vitro.

Over the past several years, OA epigenetics researchers have begun adding context to the growing body of OA epigenetics data by performing multi-omics analyses of relevant OA tissues. A good example of this was a study performed by Steinberg and colleagues [49] in 2017, which combined epigenetic analysis via 450 k array with gene expression via RNA-Seq and proteomics by mass spectrometry of eroded and intact cartilage from knee and hip OA samples (29 knee samples and 9 hip samples). They identified 209 differentially abundant proteins, 349 differentially expressed genes and 9896 DMPs. Among 49 genes identified as differentially regulated in OA in at least 2 -omics domains, 16 of those were never previously implicated in OA progression. Pathways associated with OA included extracellular matrix degradation, collagen catabolism and angiogenesis. In the final analysis, three genes stood out as differentially regulated across all three -omics levels: AQP1, COL1A1 and CLEC3B, all three previously associated with OA either in humans or animal models.

OA is not a strictly human disease, although curiously it does not affect all mammals. To date, there have been two studies utilizing the Illumina 450 k arrays for OA studies in non-human primates, both conducted on baboons, which develop OA during aging at a similar rate to human patients. The first, published by Housman and colleagues in 2019, studied knee OA from five adult female baboons compared to animals without macroscopic OA. Approximately 44% of 450 k probes aligned to the baboon genome. Out of those matched probes, filtering was necessary to produce reliable results, after which 6 significant DMPs were identified between OA cartilage samples and healthy cartilage samples. Curiously, the most differentially methylated genes included RFXAP and RUNX1 genes which have been previously linked to OA pathophysiology. Additionally, probes linked to KLHL26, RFXAP, MIR497, MIR195, ELF1, ACSL1, and CMIP genes were found in baboons but were never previously linked to human OA, therefore meriting further inspection. Housman and colleagues also published a newer expanded Baboon study using the Human Methylation EPIC (850 k) array from Illumina, discussed in the next section.

With so many studies of knee and hip OA using 450 k array, it has become increasingly difficult to make sense of the voluminous (and sometimes disparate) data regarding DNA methylation aberrations and OA pathogenesis. Many variables, particularly during data analysis, contribute to these differences, including differential filtering of methylation probe data, various methods of statistically significant P value thresholds and corrections and various methods of preprocessing data. An expertly conceived review article by Reynard and colleagues [50] attempted to rectify various 450 k OA data. In it, the authors identified 34.1% similarly DM sites in knee and hip OA comparing the Rushton et al., 2014 [37] and den Hollander et al., 2014 [39] studies. Additionally, the study found ~71% similarity of DM sites between intact and eroded OA hip cartilage when compared to DM sites of intact and eroded hips and knees combined [38,51]. Further, they identified overlaps in a third of DM sites of healthy and eroded knee cartilages compared to similar hip cartilage [37,48] specimens. An overview of the pathway overlaps identified in this review is presented in Fig. 1. These findings both strengthened the case for an underlying epigenetic contribution to OA pathology, as substantial overlap was noted in both individual DM sites as well as DM-associated gene pathways, as well as highlighting the need for consensus processing procedures to reduce inter-laboratory variability in future large-scale epigenetic analyses.

Infinium HumanMethylation850 (EPIC) BeadChip

Only a couple of studies have so far been published using the newer Illumina EPIC array, which expands coverage to >850,000 CpG sites throughout the genome (Fig. 1). The first by Housman and colleagues examined genome-wide DNA methylation patterns in right distal femoral cartilage and subchondral bone of adult baboons and compared them to healthy controls [52]. Overall, they found that cartilage sample comparison showed more DMPs and DMRs than bone tissue, a curious finding given previous data regarding expanded DNA methylation aberrancies in human subchondral bone compared to cartilage [42,44]. There were a total of 39 DMPs that separated healthy bone from OA bone, while 4298 DMPs were noted in comparison of healthy to OA cartilage. After gene ontological (GO) analysis of those cartilage DMPs, several vital OA-related biological processes were enriched. Among those were skeletal system, connective tissue, and cartilage development. Although many of the CpGs associated with OA in this study diverged from previously published human OA epigenetic associations, the study managed to find overlapping loci that were conserved between the species and play a role in organ system development, cell differentiation, inflammatory responses, and Wnt signaling.

The second was work by our group, a pilot study using peripheral blood mononuclear cell DNA methylation patterns to develop models predictive of OA progression [53]. In it, we identified 58 individuals from the large US-based osteoarthritis initiative (OAI) cohort who had evidence for radiographic progression in the 24 months following their entry into the cohort and matched with 58 non-progressive controls. Using a combination of 450 k and EPIC arrays, global DNA methylation patterns in peripheral blood mononuclear cells from baseline blood draws were quantified. Using a machine learning-based approach, we then developed and tested methylation-based models to predict future progression; models based on DNA methylation M-values (a −log10--transformed beta value) exhibited significant accuracy in discriminating future progressors from nonprogressors (mean accuracy 73%, area under receiver operator curve ROC-AUC = 0.81) compared to previous biochemical biomarker models (e.g. the Foundation for the NIH Osteoarthritis Biomarker Consortium, ROC-AUC = 0.631 [54]).

Agilent Human promoter microarray 244 K

Although the majority of DNA methylation microarray studies have been performed using Illumina products, there has been one study published utilizing the Agilent promoter microarray, a whole-genome DNA methylation microarray designed to analyze DNA derived from methylated DNA immunoprecipitation (MDIP) approaches, which covers ~240,000 CpG sites across 28,500 CpG islands. One study that utilized this approach examined the knee cartilage methylation difference between the severe and mild cases of human knee OA [55]. When low vs. high scoring OARSI knees were compared, 1214 DM genetic targets were discovered with 1070 of those hypermethylated and 144 hypomethylated. Among the genetic targets of importance were IL32, TGF-beta2, RUNX3, BMP7, WNT11, IL1RN. Interestingly, BMPs are important contributors to cartilage development [56], and BMP2, 4, and 7, in particular, regulate chondrocyte proliferation and matrix synthesis with the assistance of MMP13 [57,58]. Meanwhile, TGF-beta plays a significant role in SOX9 mRNA expression [59] which is one of the key modulators of OA pathogenesis [60]. Reassuringly, many of these genes have been previously found as differentially methylated using other epigenetic approaches, including the Illumina arrays.

Noncoding RNA analyses in OA

Agilent arrays

We will move now to a discussion of another canonical epigenetic control mechanism, noncoding RNA (Table 2). Originally, microRNA (miRNA) levels were quantified using qRT-PCR and individual primers. However, arrays by Ilumina and Agilent were developed and subsequently entered widespread use, they allow screening of hundreds of miRNA sequences simultaneously. The first large-scale array related to OA was done by Iliopoulos et al. [61]. Using the TaqMan array, they defined the expression pattern of 365 miRNAs from 33 hip and knee OA cartilage samples. They found 16 differentially expressed miRNAs in knee OA, with 9 of them up- and 7 downregulated. Interestingly, this study correlated Body Mass Index (BMI) with miRNA expression changes in OA patients; among the 16 differentially expressed miRNAs, 5 were correlated to BMI. Among these 5, miRNA-22 and miRNA-103 were positively correlated, while miRNA-25, miRNA-337, and miRNA-29a were negatively correlated. This study also performed protein analysis of various miRNA-targeted genes, identifying PPARA, BMP7, IL1B, ITGA5, and SREBP1 proteins as correlated to BMI. A key finding of this study was regulation of BMP7 and PPARA by miRNA-22; furthermore, when miRNA-22 was inhibited, PPARA and BMP7 were upregulated which blocked inflammatory signaling by IL1B, while also inhibiting catabolic gene expression including MMP13. Ultimately, by matching the differentially expressed miRNAs with their putative protein targets, 17 miRNAs were identified as likely pathogenic, including miR-22:PPARA/BMP7, miR-140:ADAMTS5, miR-16:TPM2, miR223-GDF5, miR-509:SOX9.

Table 2

Summary of key OA-related noncoding RNAs.

ncRNA	Expression level of miRNA in OA	Tissue of interest	Putative miR function and effects on OA
miRNA-140	Decreased	Cartilage	miR140−/− mice demonstrate an age-related OA phenotype [64]. miR-140 intraarticular injection into trauma-induced early-stage rat knee OA reduces progression [115]
miR-9	Decreased	Cartilage	Regulates MMP13. Injection of miR-9 agomir intraarticularly into OA rats reduces MMP13 expression and reduces OA pathology [77]. Exosomal miR-9-5p injected intraarticularly into a rat OA model reduces inflammation and OA histopathology [76].
miR-34a	Increased	Plasma, cartilage, synovium	Enhances chondrocyte apoptosis and senescence [116]. Increased in obese OA patients' and high-fat diet-induced obese mouse serum. Injection of miR-34a-5p mimic intraarticularly induces an OA phenotype, whereas miR-134a-5p antisense oligonucleotide intraarticular injection protects mice from developing OA following DMM surgery [70].
miR-146a-5p	Decreased in end-stage OA cartilage, increased in early-stage OA [117]. Increased in OA patient serum [66]	Cartilage, serum	Potential serum biomarker of OA [66]. May regulate cartilage homeostasis via Camk2d and Ppp3r2 [118].
miR-335–3pmiR-199a-5pmiR-671–3pmiR-1260 bmiR-191–3pmiR-335–5pmiR-543	Both increased and decreased.	Plasma	miRNAs associated in large assays with early OA or OA progression, potential biomarker(s) for OA diagnosis and/or prediction of progression [68].

Another study that used large-scale analysis of miRNA involvement in OA was done by Diaz-Prado et al. [62]. Using the larger Agilent Human miRNA Microarray chips, they analyzed 723 miRNAs in cartilage micropellets cultured from human knee OA patient samples. Seven miRNAs showed statistically significant differential expression with one of them being upregulated in OA (hsa-miR-483) and 6 downregulated (hsa-miR-149, hsa-miR-582–3p, hsa-miR-1227, hsa-miR-634, hsa-miR-576–5p, hsa-miR-641). Further gene ontology analysis showed that these were associated with TGF-beta, Wnt signaling, MAPK expression, mTOR signaling, focal adhesion, and regulation of actin cytoskeleton, all pathways previously associated with OA.

A seminal OA microarray-based miRNA study was published by Ntoumou et al. [63] in 2017. This study analyzed miRNA patterns in the serum of OA patients as potential diagnostic biomarkers. They scanned 2549 miRNAs using Agilent's SurePrint G3 Human miRNA, 8 × 60K microarray (Agilent Technologies Inc., Santa Clara, CA, USA). Among those, 279 miRNAs were differentially expressed in serum (205 up- and 74 downregulated) of OA patients compared to healthy controls. The study failed to produce a significant hierarchical clustering and was unable to show a clear difference between OA and control serum samples. However, in the past, it has been reported that hsa-miR-140 plays a vital role in homeostasis and chondrogenesis [64] and has been differentially expressed in OA cartilage [65], but not in the serum of OA patients. This study found hsa-miR-140–3p differentially expressed in serum which opens the door to the potential use of miRNAs as an OA diagnostic tool.

Several other microRNA surveys have since been conducted. In 2020, Rosseau and colleagues evaluated the circulating miRome of 10 women with knee OA and 10 age-matched healthy controls. They identified serum miR-146-a-5p as significantly increased in women with OA compared to controls; further, they noted the likelihood of OA development was linked to the level of miR-146a-5p; that is, for each quartile increase in serum miR-146a-5p levels, the odds ratio for OA development incrased by 1.866. A 2019 study by de Almeida and colleagues evaluated miRNAs within human cartilage, comparing macroscopically eroded and intact cartilage from both OA knees and hips, using the more powerful technique of RNA sequencing [67]. They identified 142 miRNAs as differentially expressed; following computational prioritization, 62 miRNAs were localized in a regulatory network. This network was then characterized by pathway enrichment analysis as related to ‘nervous system development’ with a particularly enriched target of NTF3, a neuronal gene related to nerve growth factor, itself the target of OA therapeutics in clinical trials. A study published by Ali and colleagues in 2020 also used a next generation sequencing approach to the identification of differentially expressed miRNAs within plasma from a much larger cohort (n = 91) of patients with early or late knee OA, defined radiographically [68]. They identified 97 miRNAs which were highly associated with OA progression (e.g. 85% or more of plasma samples demonstrated a change in the median expression in the early OA vs. late OA groups); among these, seven miRNAs were present in at least 95% of samples, including has-miR-335–3p, has-miR-199a-5p, has-miR-671–3p, has-miR-1260 b, has-miR-191–3p, hisa-miR-335–5p, and has-miR-543.

A growing body of literature has also demonstrated the pathogenic and/or disease-modifying potential of specific miRNA. For example, knockdown via antisense oligonucleotide of miR-181a-5p attenuates murine OA [69], whereas miRNA-34a-5p promotes joint destruction in murine models [70] Two homologous miRNAs, miR-204 and miR-211, work in concert to maintain joint homeostasis and prevent OA development, principally through controlling Runx2 expression [71]. miRNA-128a worsens murine OA by suppressing chondrocyte autophagy via controlling Atg12 expression [72] and in human chondrocytes via interactions with WISP1 suppressing apoptosis and inflammation [73]. miR-29a reduces synovitis in synovium of human knee OA patients and has been shown in vitro to target VEGF [74]. Several lines of evidence have demonstrated that miR-9 is associated with OA development in human patients [75], and in a rat OA model, intraarticular injection of miR-9-5p reduces OA pathology by downregulating SDC1 [76] and suppresses the key matrix metalloproteinase MMP13. Many other miRNAs have been linked, either in vitro or in vivo with OA development (recent miRNA-specific review articles [78,79]).

Although miRNA are the most commonly studied noncoding RNA species, many others have been described. For example, Fu and colleagues [80] in 2015 evaluated long noncoding RNA (lncRNA) expression patterns in eroded cartilage from knee OA patients compared to disease-free cartilage from control subjects using the Human lncRNA Array v2.0 (8 × 60 k) Arraystar microarray chip, which can detect 33,045 lncRNAs and 30,215 mRNAs. They identified 15,099 lncRNAs as differentially expressed in OA cartilage with 3007 of them at least two-fold upregulated and 1707 downregulated in comparison to healthy controls. Using these data, they then constructed a coding-non-coding gene coexpression network and identified 48 lncRNAs that were predicted to have more than five cis-regulated targeted genes. Gene ontology analysis revealed upregulated transcripts enriched for anatomical structure development, cell periphery, and lipid binding, while downregulated transcripts were enriched for metabolic processes, intracellular parts, and RNA binding.

Circular RNAs have also been evaluated in OA, and function principally as ‘sponges’, binding to and removing various miRNAs. CDR1as binds miR-641 and promotes OA progression in human patients [81]. CircRNA-CER is a sponge for miR-136, thereby increasing the expression of the matrix metalloproteinase and key OA effector MMP13 [82]. Hsa_circ_0045,714 binds miR-139 b and is downregulated in OA, thereby decreasing expression of collagen and aggrecan [83]; similarly, CircSERPINE2 is downregulated in OA, typically binds to miR-1271 which targets ERG, resulting in increased expression of matrix metalloproteinases MMP3, MMP13, and ADAMTS-4 [84]. Numerous other associations have also been described, as outlined by a recent review article [85].

Next generation sequencing (NGS)-based approaches to epigenetics analysis

Next-generation sequencing (NGS) has emerged as a useful tool in the study of epigenetics over the past decade (Fig. 3). Compared to array-based methods, NGS approaches offer several advantages. First, the design of methylation arrays requires prior knowledge of the region to be targeted and includes only a small subset of genomic CpG sites, generally in regions thought to be important for gene regulation. NGS approaches, particularly whole-genome bisulfite sequencing, overcomes these limitations by providing coverage of nearly the entire genome. Disadvantages of NGS include a requirement for increased sequencing depth, as bisulfite treatment reduces the genomic code to a 3-base (ATG) rather than a 4-base (CATG) system, difficulty in sequencing through certain genomic locations, and increased cost. There are several variations of NGS-based epigenetic analyses, which we discuss in more detail below.

Fig. 3.

Timeline of major advances in methodologies for studying DNA methylation using next generation sequencing-based techniques, with article references.

Reduced representation bisulfite sequencing (RRBS)

RRBS was one of the first NGS-based epigenetic methodologies to be applied to OA. This approach relies on restriction enzyme digestion, most commonly with MspI (targeting CCGG motifs), followed by adapter ligation and sequencing [86]. This approach substantially reduces the total number of nucleotides to sequence by approximately 99% while still capturing most promoter regions and repeating sequences [87]. Although cost-effective and practical, (requiring less than 300 ng of DNA input) this approach to whole-genome sequencing only recognizes 40–220 base pair regions and leaves out a plethora of cleaved regions that are lengthier. Ultimately, this approach still provides a better depth and coverage than bisulfite sequencing and array approaches [88,89] thus making it a valuable tool in epigenetic research.

The first study to utilize RRBS in OA cartilage examination was done by Bonin et al., in 2016 and identified more than 1000 DM regions (39,322 DM sites) in eroded knee OA cartilage compared to intact cartilage from the same joint, with 442 of those being hypermethylated and 560 hypomethylated (Fig. 3) [90]. Using a sliding window technique, 19 DMRs were found to contain 5 or more DM sites in first exons/introns and gene promoters. The group went on to confirm gene expression effects of differential DNA methylation via both qPCR of primary cartilage tissue, and in vitro in chondrocyte cell line culture following treatment with the DNA methylation inhibitor 5-Azacytidine. Five of 6 hypermethylated genes were noted to have reduced expression, although only 1 of 6 hypomethylated genes had corresponding increases in gene expression level.

Analyses of hydroxymethylation in OA

Until very recently, epigenetic analyses of human disease have focused on only one modification of genomic DNA: methylation. However, other covalent modifications exist as intermediaries between ‘fully’ methylated and unmethylated cytosine states during active DNA demethylation and include 5′-hydroxymethylation (5-hmC), -formylation, and -carbosylation. Although 5′-formylcytosine and 5′-carboxylcytosine are rapidly converted to cytosine by the action of TET and TDG enzymes, 5-hmC is more stable and is thought to exist as a sort of transition state between epigenetic activation and inactivation [91], although recent studies suggest it may exist as a stable and heritable epigenetic mark [92]. The presence of 5-hmC varies widely between tissues, as frequent as ~0.7% in brain tissue and as low as ~0.05% in heart tissue [93].

Quantitation of hydroxymethylation involves methods such as Oxidative Bisulfite sequencing (oxBS-seq), Tet-assisted Bisulfite Sequencing (TAB-seq), reduced representation hydroxymethylation profiling (RRHP), and 5-hmC pull-down sequencing. OxBS-Seq relies on chemical conversion of hydroxymethylated cytosine to formylcytosine, which is then converted via bisulfite treatment to cytosine. ‘Traditional’ bisulfite sequencing (see introduction earlier in this article) is performed in parallel, and hydroxymethylated sites are computed as the difference between oxBS sequences and bisulfite-treated sequences. Conversely, in TAB-Seq, genomic 5-hmC is first ‘protected’ with glucosylation, then all methylated cytosines are converted using Tet1 enzyme; remaining cytosines are interpreted as hydroxymethylated in the original sample. Both approaches have substantial disadvantages such as loss of DNA due to degradation in the additional oxidative step (oxBS-seq), and reliance on the reduced efficiency of TET enzyme compared to ‘traditional’ bisulfite conversion (leaving less than 10% of methylated residues unconverted compared to <1% unconverted) [94,95]. Analogous to the RRBS discussed previously, researchers have also developed reduced representation hydroxymethylation profiling (RRHP), which combines protection via glucosylation of 5-hmC sites followed by enzymatic digestion and sequencing using RRBS protocols. RRHP allows mapping of 5 hmC sites at a single-base resolution that avoids harsh chemical conversion processes with increased sequencing quality [96].

Only a handful of publications have evaluated sequence-specific hydroxymethylation patterns in the context of OA. Following on a 2014 publication where they demonstrated global increases in 5-hmC among OA patient cartilage (Fig. 3) [97], in 2015, Taylor et al. [98] analyzed 5-hmC patterns in monolayer chondrocyte cultures from patients with a history of knee OA and disease-free articular chondrocytes. They demonstrated widespread increases in 5-hmC associated with OA, particularly in gene bodies and intergenic regions; in total, 70,591 regions demonstrated differential hydroxymethylation. The investigators coupled their hydroxymethylation analysis with gene expression via microarray and demonstrated that increases in 5-hmC in gene bodies was associated with transcriptional upregulation of a number of key OA-associated genes including MMP3, LRP5, GDF5, and COL11A1. In 2020, the same group expanded their analysis of 5-hmC to an OA animal model (disruption of the medial meniscus or DMM). Similar to their human studies, they identified a gain of ~40,000 differentially hydroxymethylated cytosine sites throughout the mouse genome following OA induction [99]. Importantly, these epigenomic changes occurred at 6 weeks following surgery, prior to the onset of histopathologic OA. Also similar to human analyses, the bulk of these sites were located within gene bodies with only a small fraction occurring within promoter regions. Furthermore, when TET1 was pharmacologically knocked down in vitro in human OA cartilage explants, a 25% reduction of 5-hmC was produced.

Pitfalls in DNA methylation and noncoding RNA OA research

Although we now have a large number of studies reporting both differential DNA methylation and noncoding RNA results throughout the OA field, there are a number of hurdles that must be overcome and pitfalls that must be avoided both in analyzing these epigenetic data within the context of a single study and in comparing and synthesizing data from multiple studies. In human studies, the tissue studied has a significant impact on both DNA methylation and noncoding RNAs. Roughly half of human OA cartilage studies, for example, include comparisons between eroded and intact sections of the same joint, whereas the remainder compare eroded or intact sections from OA samples to disease-free specimens. As it is both quite invasive and likely OA-inducing to obtain cartilage specimens from living healthy donors, these disease-free controls are routinely sourced either from cadavers (raising concerns about epigenetic shifts due to anoxia in the post-mortem period) or from patients undergoing joint replacement for a non-OA condition, e.g. neck-of-femur fracture patients (raising concerns about epigenetic shifts due to other underlying medical conditions, such as osteoporosis).

From a data analysis standpoint, a key concern is setting the threshhold of what one characterizes as ‘statistically significantly different’ within epigenetics data. A particular problem within array-based DNA methylation approaches and next-generation bisulfite sequencing, the question of where the false discovery rate (FDR) correction threshhold for statistical significance is placed is a contentious one, mirroring debate within the larger scientific community regarding type 1 error correction [100]. Generally speaking, rigorous Bonferroni or Sidak adjustment are thought overly conservative and screen out too many true positive results. The Benjamini-Hochberg FDR method is the most commonly reported among OA literature, but it does not account for any correlation among DMPs within underlying DNA methylation datasets [101]. Some have gone so far as to recommend no FDR correction at all for array-based DNA methylation analyses; rather, simply setting a significant threshhold of e.g. P < 9E-8 [102] using a Student t-test. A more sophisticated view, and one that is likely to become more prevalent in the coming years, is to analyze regional differences in DNA methylation, thereby taking into account local linkage in methylation when computing statistical significance. This method, commonly referred to as bump hunting, combines both appropriate FDR control and high power for detection of differential methylation across the genome [103,104], although questions about bump stability and comparability among studies remain.

Another contentious debate in the OA epigenetics field is the threshhold for defining meaningful change in mean DNA methylation values (either among DMLs or DMRs), typically denoted as ‘delta-beta’ when describing DMLs. This term describes how ‘different’ DNA methylation values must be in order to have a biologically relevant effect, for which we unfortunately do not have empirical evidence. Generally speaking, researchers define this threshhold at between 10 and 15% difference in mean DNA methylation (beta) values between groups of interest; higher or more strict threshholds will result in a reduced number of reportable DMLs or DMRs, whereas lower values will result in increases. Cellular composition can also affect these threshholds significantly; for example, analyses of human cartilage, which consists almost exclusively of one cell type (chondrocytes) tend to be more homogeneous in their DNA methylation patterns and therefore more amenable to higher delta-beta threshholds, whereas methylation differences within a poorly-represented cell mixture (e.g. regulatory T cells in a mixed peripheral blood sample) may be diluted within the larger sample and necessitate a lower threshhold. Along these same lines, careful consideration must be paid in analyses of mixed cellularity to ensure that DNA methylation (or noncoding RNA) group differences represent true changes in the underlying epigenetic state rather than variation in cellular composition within the examined groups.

A final potential pitfall will become ever more important in the coming years with an increasing number of research projects utilizing deep sequencing of either DNA methylation or noncoding RNA research. Success in these large-scale analyses is predicated both on depth of sequencing; that is, how many nucleotide strands within each sample align to a particular location of interest, and breadth of sequencing; that is, how well the entire genome (or defined portions of the genome) were covered in a given sequencing run. DNA methylation analysis via deep sequencing using the bisulfite method presents a unique challenge in that it converts the genome to a 3-base code rather than 4-base, as unmethylated cytosines and cytosines not within CpG dinucleotides are converted to uracil. Therefore, sequence read alignment becomes more difficult and requires longer overlapping regions, thereby necessitating increased read length and sequencing depth to achieve adequate quality. Whole-genome sequencing of certain regions of the genome can also prove difficult, particularly in regions that are adenine/thymine or guanosine/cytosine-rich [105]. Amplification-free library preparation is one of the best ways to address this shortcoming, assusming that adequate DNA is available from the source tissue to obviate the need for amplification [106]. Researchers and readers alike should be cognizant of these potential sources for bias within studies, and should take time to carefully examine the specific methodologies utilized in these large-scale epigenetic analyses of OA samples, particularly when comparing studies of dissimilar methodological approach.

DNA methylation editing in OA

Despite the numerous individual CpG sites and gene pathways identified as differentially methylated in association with OA, studies elucidating the specific, causal relationships between these epigenetic changes and OA pathogenesis are still rare. These studies on OA epigenetic editing can be divided into two broad categories: analyses of global DNA methylation alteration via either pharmacologic inhibition of various key epigenetic effectors (i.e. the DNA methyltransferases) or via genetic knockouts of these effectors, or specific DNA methylation modifications, spurred by recent advances in DNA localization and targeting. From a global DNA methylation modification perspective, in 2020 Smeriglio and colleagues published a pivotal article evaluating OA in Tet1 knockout mice. This work followed up on their previous description of increases in 5-hydroxymethyl cytosine (5-hmC, the ‘intermediate’ between methylated and unmethylated cytosine states) in OA [97,98]. Remarkably, Tet1−/− mice were prevented from 5-hmC modification at 98% of the loci associated with DMM-induced OA. Furthermore, administration of 2-hydroxyclutarate, a small molecule inhibitor of Tet1, reduced OA pathology in mice and modified the expression of OA effectors in human chondrocyte culture.

Recent technological advances in epigenetic editing have, over the past several years, allowed investigators to evaluate the downstream effects of altering specific CpG loci in the context of OA. These advances hinged on the identification and further development of highly specific DNA localization/binding moieties which allow a tethered epigenetic effector to have a local effect. These include zinc finger proteins (ZFPs) [107], endogenous DNA binding proteins found in most eukaryotes that consist of ~30 amino acids, generally constructed as six-mers targeting ~18 base pairs of DNA [108]. This system is limited by construct complexity and lack of modularity, and has not yet been applied to OA studies. A second DNA targeting technology are transcription activator-like effectors (TALEs), proteins isolated from Xanthomas bacteria following plant infection [109]. Their natural function is to activate expression of specific plant genes, thereby aiding in infection. They consist of a central 33–34 amino acid repeat segment responsible for DNA binding and sequence specificity [110]. Much like ZFPs, construction of TALE arrays is technically challening, non-modular, and expensive, and thus have not been applied to OA studies. A third, much more widely applied DNA binding technology is based on the Clustered Regularly Interspersed Short Palindromic Repeats (CRISPR)/Cas9 system. Originally identified as a bacterial defense system, CRISPR/Cas9 consists of an endonuclease protein targeted to a genomic location via the binding of a short segment of RNA (guide RNA, gRNA) which provides sequence specificity [111]. Although used widely for genetic studies, particularly by producing site-specific DNA double strand breaks and point mutations/insertions, nuclease-defective Cas9 (dCas9) has emerged as an invaluable tool in epigenetic editing experiments, owing to its modularity, ease of targeting (by synthesizing various gRNAs), and high degree of sequence specificity [112].

So far, OA efforts utilizing dCas9 as an epigenetic editor have been restricted to in vitro studies. First, Parker et al. published last year a genome-wide methylation-expression quantitative loci (meQTL) analysis [113]. mQTLs are genomic loci whereby differential methylation of CpG sites act in concert with genomic polymorphisms to affect gene expression. Although several meQTLs were identified, the group went on to demonstrate that forced DNA demethylation of a specific CpG site (cg20220242) within the gene RWDD2B, via a dCas9-TET1 fusion protein-expressing human chondrocyte cell line, not only increased RWDD2B expression but could modulate the impact of the OA-associated SNP rs6516886. Earlier this year, the same laboratory published a similar study demonstrating the relationship between DNA methylation at three CpG sites within the COLGALT2 gene enhancer and the OA-associated SNP rs11583641 [114]. These two studies demonstrate the power of epigenetic editing techniques to confirm and expand on previous OA DNA methylation hypotheses, and are no doubt harbingers of many more such investigations in the future.

Conclusion

Through recent advances in sequencing technology, new epigenetic analysis methodologies, and reagent and sequencing costs, epigenetic research in OA has been exponentially increasing. Although OA pathogenesis is still not fully elucidated, increasing evidence points to a combination of genetic and environmental factors, with epigenetic modifications playing a potentially vital role in integrating these two factors. Over the past several years, our understanding of the potential epigenetic underpinnings of OA have expanded both through increasing depth of DNA methylation and noncoding RNA sequencing resolution as well as increased sample sizes of OA patient studies and greater integration of new study data with previously existing datasets. Particularly, integrative analyses have begun to shed light on the ways in which the underlying genome, DNA methylation, noncoding RNA, and chromatin accessibility interact with each other and conspire to increase OA risk both in human patients and animal models. Future studies in OA epigenetics will doubtless reveal additional insights into disease pathogenesis and potential biomarkers. Additionally, next-generation epigenetic editing approaches will no doubt expand beyond in vitro studies to in vivo work and help to further solidify our understanding DNA methylation modifications in the development and progression of OA. These insights will help both researchers in the laboratory and clinicians treating OA patients to identify patients sooner and with higher fidelity, advancing the fight against this devastating disease.