Serum microRNA array analysis identifies miR-140-3p, miR-33b-3p and miR-671-3p as potential osteoarthritis biomarkers involved in metabolic processes

MicroRNAs (miRNAs) in circulation have emerged as promising biomarkers. In this study, we aimed to identify a circulating miRNA signature for osteoarthritis (OA) patients and in combination with bioinformatics analysis to evaluate the utility of selected differentially expressed miRNAs in the serum as potential OA biomarkers. Serum samples were collected from 12 primary OA patients, and 12 healthy individuals were screened using the Agilent Human miRNA Microarray platform interrogating 2549 miRNAs. Receiver Operating Characteristic (ROC) curves were constructed to evaluate the diagnostic performance of the deregulated miRNAs. Expression levels of selected miRNAs were validated by quantitative real-time PCR (qRT-PCR) in all serum and in articular cartilage samples from OA patients (n = 12) and healthy individuals (n = 7). Bioinformatics analysis was used to investigate the involved pathways and target genes for the above miRNAs. We identified 279 differentially expressed miRNAs in the serum of OA patients compared to controls. Two hundred and five miRNAs (73.5%) were upregulated and 74 (26.5%) downregulated. ROC analysis revealed that 77 miRNAs had area under the curve (AUC) > 0.8 and p < 0.05. Bioinformatics analysis in the 77 miRNAs revealed that their target genes were involved in multiple signaling pathways associated with OA, among which FoxO, mTOR, Wnt, pI3K/akt, TGF-β signaling pathways, ECM-receptor interaction, and fatty acid biosynthesis. qRT-PCR validation in seven selected out of the 77 miRNAs revealed 3 significantly downregulated miRNAs (hsa-miR-33b-3p, hsa-miR-671-3p, and hsa-miR-140-3p) in the serum of OA patients, which were in silico predicted to be enriched in pathways involved in metabolic processes. Target-gene analysis of hsa-miR-140-3p, hsa-miR-33b-3p, and hsa-miR-671-3p revealed that InsR and IGFR1 were common targets of all three miRNAs, highlighting their involvement in regulation of metabolic processes that contribute to OA pathology. Hsa-miR-140-3p and hsa-miR-671-3p expression levels were consistently downregulated in articular cartilage of OA patients compared to healthy individuals. A serum miRNA signature was established for the first time using high density resolution miR-arrays in OA patients. We identified a three-miRNA signature, hsa-miR-140-3p, hsa-miR-671-3p, and hsa-miR-33b-3p, in the serum of OA patients, predicted to regulate metabolic processes, which could serve as a potential biomarker for the evaluation of OA risk and progression.


Background
Osteoarthritis (OA) is the most common chronic degenerative joint disease and a leading cause of pain and disability. Its prevalence is steadily increasing and is expected to be the greatest cause of disability by 2030 [1]. OA is defined as a heterogeneous disease consisting of a group of distinct joint disorders with common biological, morphological, and clinical outcomes [2]. To date, several risk factors, including aging, obesity, sex, joint injury, and heredity have been implicated in disease development and progression [3,4]. Recently, an association has been suggested between OA and metabolic syndrome or metabolic risk factors, as hypertension, overweight, dyslipidemia, and insulin resistance, and a new OA subtype, metabolic OA, has been introduced [5][6][7]. As the above metabolic factors can be modified, it becomes evident that understanding the molecular mechanisms which are involved in OA pathogenesis is essential for the discovery of OA personalized therapies. However, until now, despite the significant burden that OA imposes in patients and in healthcare systems [8], there is absence of effective non-arthroplasty treatment options for disease management. One major challenge for the development of effective therapy for OA is the detection of the disease at an early stage. To date, the standard method for disease diagnosis and severity is the radiograph, which appears to be limited to detection in late disease stages and has weakness in monitoring disease progression [9]. Therefore, the identification of noninvasive and sensitive serum biomarkers will assist OA diagnosis, prognosis, and early treatment before the onset of radiographic findings.
MicroRNAs (miRNAs) are a novel group of universally present small non-coding RNAs, approximately 22-28 nucleotides in length, which regulate gene expression at the post-transcriptional level through mRNA degradation or translational repression. Each miRNA can regulate a large number of genes by binding mainly to the 3′UTR of mRNA targets [10]. Recent studies have implicated miR-NAs as promising clinical biomarkers in various diseases including cancer, diabetes, and cardiovascular disease [11][12][13][14]. MiRNAs have also been associated with processes involved in OA pathogenesis, such as cartilage homeostasis, extracellular matrix regulation, endochondral ossification, bone metabolism, apoptosis, autophagy, inflammation, and lipid metabolism [15,16]. Our group was among the pioneers to suggest the involvement of miRNAs in OA pathogenesis by demonstrating the differential expression of 16 miRNAs in OA compared to healthy articular cartilage [17].
Recently, the use of microarrays for the global characterization of miRNA expression profiling is becoming a strong research tool that provides useful information on the pathogenesis of many diseases [12]. Several micro-RNA profiling studies have been performed in different OA tissues and have led to the identification of numerous differentially expressed miRNAs between OA joint tissues and controls [17,18]. However, in the majority of these studies, the number of miRNAs tested was limited, ranging from 157 to 723 miRNAs, compared to 2588 mature miRs recently released from the Sanger Institute and University of Manchester miRbase (release June 2014) allowing for many more miRNAs to be identified [19].
Moreover, miRNAs have been discovered in circulation, however, the origin and the biological functions of the circulating miRNAs in a disease remain poorly understood. Evidence suggests that circulating miRNAs could be transferred in the circulation through microvesicles, exosomes, Ago protein complexes, or HDL and that they are either byproducts of the cellular content or are mediating the inter-cellular signaling reflecting the pathophysiological state of the cell [20][21][22][23]. Alterations in their expression levels in abnormal conditions along with their stability in circulation and surprisingly long half-life of approximately 5 days in plasma [24] have implicated them as attractive molecules for disease prognosis and/or progression [11]. Circulating miRNAs have been demonstrated as important biomarkers in many diseases, as diabetes type II [25] metabolic syndrome [26] Alzheimer's, hypertension, osteoporosis, Parkinson's disease, and cancer [12]. Recently, a limited number of studies highlighted the role of circulating miRNAs as potential biomarkers for assessment of prognosis and/or progression in osteoarthritis [23,27,28]. However, a systematic analysis of serum miRNA profiling in osteoarthritis has not yet been performed.
In the present study, we performed miRNA profiling using high-density miRNA-arrays in serum of OA patients and in combination with bioinformatics analysis evaluated the utility of selected differentially expressed miRNAs as potential OA biomarkers.

Patient samples
Serum samples were collected from 12 patients with primary OA (9 females, 3 males, ages 69.83 ± 4.83 years) undergoing total knee replacement surgery at the Orthopaedics Department of the University Hospital of Larissa. Radiographs were obtained before surgery and were graded using the Kellgren-Lawrence system according to the following criteria: grade 1 (doubtful narrowing of joint space and possible osteophytes), grade 2 (definite osteophytes and possible narrowing of joint space), grade 3 (moderate multiple osteophytes, definite narrowing of joint space and some sclerosis, and possible deformity of bone ends), and grade 4 (large osteophytes, marked narrowing of joint space, severe sclerosis, and definite deformity of bone ends). All OA patients had a Kellgren-Lawrence grade ≥ 3. The assessment of the radiographs by two independent expert observers was blinded. Patients with rheumatoid arthritis and other autoimmune disease as well as chondrodysplasias, infection-induced OA, and posttraumatic OA were excluded from the study. As controls, serum samples were obtained from 12 healthy individuals (6 females, 6 males, ages 64,25 ± 5,04 years) undergoing knee fracture repair surgery, with no history of joint disease, who did not show clinical manifestations compatible with OA when specifically explored by radiography. The clinical characteristics of OA patients and the control cohort are shown in Table 1. All individuals who participated in the study were of Greek origin. The study protocol conformed to the ethical guidelines of the 1975 Declaration of Helsinki as reflected in a priori approval by the local ethical committee of the University Hospital of Larissa.
Articular cartilage samples (n = 12) were collected from femoral condyles and tibial plateaus of the same OA patients undergoing total knee replacement surgery. Normal articular cartilage was obtained from small free cartilaginous fragments of 7 out of the 12 healthy individuals undergoing knee fracture repair surgery. Signed informed consent was obtained from all OA and healthy individuals prior to surgery and/or blood drawing.

RNA extraction and microRNA expression profiling by microarray in serum of OA patients and healthy individuals
Total RNA was extracted from serum samples using the RiboEXTMLS kit (GeneAll, Seoul, Korea). RNA quantity and quality were evaluated using a spectrophotometer (NanoDrop® ND-1000 UV-Vis, Nanogen Inc.). Total RNA samples were spiked using the microRNA Spike-In Kit (Agilent Technologies Inc., Santa Clara, CA, USA) to assess the labeling and hybridization efficiencies. Labeling and hybridization were performed using the miRNA complete Labeling and Hybridization Kit (Agilent Technologies Inc., Santa Clara, CA, USA) according to manufacturer's instructions. Samples were hybridized to the SurePrint G3 Human miRNA, 8X60K platform (miR-Base release 21.0, Agilent Technologies Inc., Santa Clara, CA, USA) containing probes for the detection of 2549 human miRNAs. Images were scanned using Agilent Microarray Scanner (G2565CA) controlled by Agilent Scan Control 7.0 software. The signal after background subtraction was exported directly into Agilent Feature Extraction Software version 4.0.1.21 (Agilent Technologies Inc., Santa Clara, CA, USA). Normalization was performed using quantile algorithm. MicroRNAs were considered as differentially expressed (DE) if they obtained a p value < 0.05 and a false discovery rate (FDR) ≤ 0.05.

ROC analysis
Receiver Operating Characteristic curves (ROC) analysis was performed using the MATLAB® simulation environment and was used to calculate the area under the curve (AUC) value along with standard error and 95% confidence intervals. ROC curves were considered significant with AUC value > 0.8 and a p value < 0.05.
Quantitative real-time polymerase chain reaction (qRT-PCR) validation of selected miRNAs in serum of OA and healthy individuals The expression levels of 7 selected miRNAs screened with miRNA microarrays, as mature hsa-miR-33b-3p, hsa-miR-4284, hsa-miR-663a, hsa-miR-150-5p, hsa-miR-1233-3p, hsa-miR-140-3p, and hsa-miR-671-3p, were evaluated in serum samples from 12 OA patients and 12 healthy controls. More specifically, reverse transcription was carried out with 100 ng of total RNA using MiScript II Reverse Transcription Kit (QIAGEN Inc., Valencia, CA, USA). The relative quantification of selected differentially expressed miRNAs was performed by qRT-PCR reaction with the miScript SYBR® Green PCR Kit (QIAGEN Inc., Valencia, CA, USA) and MiScript Primer Assays (QIAGEN Inc., Valencia, CA, USA), using an ABI 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Reactions were performed in duplicate. The relative expression levels were calculated using the 2 −ΔΔCT method. Normalized gene expression values for each gene were generated based on cycle threshold (CT) values for each of the genes. Hsa-miR-25-1 was used as endogenous control, as it exhibited minimum variance between the controls and patients in miRNA microarray analysis. Also, its stability was confirmed using the NormFinder software [29].

Primary cultures of OA and normal and articular chondrocytes
Articular cartilage was dissected and subjected to digestion with 1 mg/ml pronase (Roche Applied Science, Mannheim, Germany) for 30 min and then the sample was centrifuged and the pellet was subjected to digestion with 1 mg/ml collagenase P (Roche Applied Science, Mannheim, Germany) for 3 h at 37°C. Chondrocytes were counted and checked for viability using trypan blue staining. More than 95% of the cells were viable after isolation. Chondrocytes were cultured with Dulbecco's modified Eagles medium/Ham's F-12 (DMEM/F-12) (GIBCO, BRL, UK) plus 5% fetal bovine serum (FBS, GIBCO, BRL, UK) and 100 U/ml penicillinstreptomycin, and were incubated at 37°C under a humidified 5% CO 2 atmosphere until reaching confluence. Chondrocytes were kept in culture for one passage, while types I and II collagen ratios were evaluated in all samples to exclude dedifferentiation events.

RNA extraction and quantification of miRNA expression from OA and healthy articular cartilage
Total cellular RNA was extracted from cultured chondrocytes using Trizol reagent (Invitrogen, Life Technologies, Paisley, UK). RNA quantity and quality was evaluated using a spectrophotometer (NanoDrop® ND-1000 UV-Vis, Nanogen Inc.). Reverse transcription was conducted using the SuperScript III Reverse Transcriptase kit (Invitrogen/ Thermo Fisher Scientific Inc., USA), using primers specific for hsa-miR-33b-3p, has-miR-4284, hsa-miR-663a, hsa-miR-150-5p, hsa-miR-1233-3p, hsa-miR-140-3p, hsa-miR-671-3p, and U6 small nuclear RNA (RNU6B) stemloop RT primer to generate the cDNA as previously described [30]. MiRNA expression was verified by qRT-PCR. Reactions were done in duplicate. The reactions were performed using Power SYBR Green PCR Master Mix (Applied Biosystems) and specific primers (forward and reverse) for each miRNA. For the quantification of the relative expression of each miRNA, the threshold cycle (Ct) values were normalized against the endogenous reference U6. The 2 −ΔΔCT method was used [30].

Statistical analyses and data analyses
Multiparameter analyses were performed with the MATLAB ® simulation environment (The Mathworks Inc., Natick, MA, USA). In particular, data pre-processing was performed in Microsoft Excel ® and data processing was performed within the Matlab® v.7.6.0 (The MathWorks Inc. Natick, MA, USA) computing environment, using the Bioinformatics Toolbox. For background correction, the well-performing multiplicative background correction (MBC) approach was followed [31]. In terms of filtering, as low signal intensity measurements are less reliable in terms of the impact of noise on them, than high gene expression measurements, an intensity-dependent filtering [32,33] with an absolute threshold value of 10 was used, in order to exclude low quality features. A threshold of 1 was set as a cut-off value, meaning that spot intensity should be at least the same as that of the background. The quantile normalization method was used for further processing [34]. In order to identify potentially differentially expressed (DE) genes between samples and among genes of the same experiment a two-tailed Student's t test was used to test the mean differences between the two groups. Genes that received a p < 0.05 were considered as DE. Further on, the false discovery rate (FDR has been calculated as described previously [35][36][37]. There was a FDR of 0.007% for p < 0.05. Continuous variables are expressed as median ± standard deviation unless indicated differently. k-means and hierarchical clustering (HCL) analyses of microRNA expression were performed using all differentially expressed (DE) miRNAs with Euclidian distance.
RNA target prediction for selected miRNAs was performed with TargetScan v. 7.0 [38] and DIANA Tools, a collection of web-tools for miRNA functional annotation analysis [39]. Functional annotation was performed with the Webgestalt web-tool [40] and DIANA tools. Furthermore, KEGG pathway enrichment analysis was conducted for the target genes using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tools [41] with the cut-off criterion of p < 0.05.
Venn diagrams were generated with the Venn Diagram Calculator from the Bioinformatics, Evolutionary Genomics Group of the University of Gent.

MicroRNA microarray profiling in serum of OA patients and healthy individuals
Among the 2549 miRNAs tested using the Agilent 8 × 60 K miRNA-array platform, a total of 279 miRNAs were found to be differentially expressed (205 upregulated and 74 downregulated) (FDR < 0.008, p < 0.05) in the serum of OA patients compared to healthy individuals (Fig. 1a, b).
Classification methods were used in order to detect patterns in gene expression. Unsupervised two-way hierarchical clustering (HCL) with Euclidian distance and k-means analyses could not discriminate accurately between OA and normal serum samples (data not shown).
ROC curve analysis of the seven selected miRNAs is presented in Fig. 2.

Pathway analysis and target gene prediction
Pathway enrichment analysis of the 77 miRNAs that were found significant in ROC analysis revealed that these DE miRNAs are potentially involved in pathways associated with OA pathogenesis, as among others, thyroid hormone, FoxO, mTOR, sphingolipoid, MAPK, PI3K-Akt, Wnt, ErbB, estrogen and TGF-beta signaling pathways, ECMreceptor interaction, and fatty acid biosynthesis. Annotated pathways are summarized in Table 2.

Disease annotation
Based on the previous list of the revealed miRNAs, we further investigated their association with known diseases. It appeared that the miRNAs we identified participate in arthritis, joint diseases, and rheumatic disease. Interestingly, our identified miRNAs included hsa-miR-140-3p, which participated in arthritis, joint diseases, and rheumatoid disease. On average, hsa-miR-140-3p was found to be downregulated in all OA samples. The results of the Disease Annotation are summarized in Table 3.

Discussion
To our knowledge, this is the first study to establish a potential global serum miRNA signature in OA patients using a high-resolution microarray technology interrogating 2549 miRNAs. We identified a three-miRNA signature, including hsa-miR-140-3p, hsa-miR-671-3p, and hsa-miR-33b-3p, which was downregulated in the serum of OA patients compared to controls and in silico predicted to be involved in regulating metabolic factors. Recent evidence points that circulating miRNAs could derive from tissue injury [42] or adipose tissue depots [43], it could thus be speculated that serum miRNAs in OA could potentially have derived from cartilage degeneration, apoptosis, or inflammation, and miRNA expression profiling in OA serum could thus constitute a disease fingerprint.
To date, limited studies have been conducted, mainly using low-density arrays and qRT-PCR for the identification of circulating miRNA expression profile in serum, plasma or blood of OA patients, with inconsistent results among different studies. More specifically, Okuhara et al., investigated the expression of hsa-miR-146a, hsa-miR-155, hsa-miR-181a, and hsa-miR-223 by qRT-PCR in peripheral blood mononuclear cells of OA patients and healthy controls and correlated their expression to disease progression [44]. In a subsequent study, the expression of 380 miRNAs using TaqMan low density arrays was evaluated in the plasma of patients with primary osteoarthritis and controls and 12 miRNAs (hsa-miR-16, hsa-miR-20b, hsa-miR-29c, hsa-miR-30b, hsa-miR-93, hsa-miR-126, hsa-miR-146a, hsa-miR-184, hsa-miR-186, hsa-miR-195, hsa-miR-345, and hsa-miR-885) were identified being over-expressed in OA patients [28]. A more recent study by Beyer et al. [27] identified let-7e in the serum of a large population-based cohort of OA patients necessitating arthroplasty as a potential predictor for severe knee or hip osteoarthritis; however, its expression was not confirmed in OA articular cartilage samples. The inability among the published studies to reveal consistent differentially expressed circulating miR-NAs in OA could be attributed to differences in the study design, as selection of patients, disease status, and differences in ethnicity [45], and to different platforms and bioinformatics tools used in each study.
In the current setting among the 2549 miRNAs screened using a high-density microarray platform, we identified 279 miRNAs of which 205 upregulated and 74 downregulated, that were differentially expressed in serum between OA patients and healthy individuals. Hierarchical clustering analysis (HCA) using the miRNA signatures failed to demonstrate clear differences between OA and control serum samples. This inability of HCA might be due to the small sample size or to other technical issues, as the low dynamic range of miRNA microarrays compared to qRT-PCR to accurately identify fold-changes for miRNAs present in both high and low abundance [46].
In order to evaluate the diagnostic role of the circulating miRNAs in OA serum, we established ROC curves using a strict filtering approach. We determined 77 miR-NAs that showed significant sensitivity and specificity with area AUC values > 0.8, suggesting those miRNAs as potential OA biomarkers. Many of these miRNAs have yet unknown roles that remain to be revealed. To evaluate the biological information provided by the 77 miR-NAs through ROC analysis, we performed pathway enrichment analysis, which revealed that their target genes were commonly enriched in OA related pathways, such as ECM-receptor interaction, mTOR, PI3K/Akt, Wnt, TGF-β and adipocytokine signaling pathways, insulin resistance, FoxO, autophagy, and fatty acid metabolism ( Table 2).
Hsa-miR-140 plays an important role in cartilage homeostasis and chondrogenesis [47][48][49] and has been previously reported with differential expression in OA cartilage [44,50] but not in serum. Recently, it was  shown that miR-140 is required for adipogenesis [51]and that decreased levels of miR-140 were found in the plasma of patients with morbid obesity influencing TGFbR1 expression levels [51], suggesting a new role of miR-140 in metabolic processes and obesity, main contributors in OA development. Disease Association Annotation in our study revealed that among the differentially expressed miRNAs, hsa-miR-140-3p was involved in arthritic diseases. Regarding hsa-miR-671-3p and hsa-miR-33b, there are no reports on their involvement in OA pathogenesis. Human miR-33 has two isoforms, has-miR-33a and has-miR-33b, which share the same seed sequence, have the same predicted targets, and are both regulating lipid and cholesterol metabolism [19,52]. Our group has highlighted the role of hsa-miR-33a in regulating cholesterol transport genes [53]. Recently, miR-33b was shown to regulate adipocyte differentiation [54] and miR-33a different functions of macrophages having a role in development and progress of atherosclerosis [55]. Regarding miR-671, besides its role in cancer cell migration [56], it was recently indicated that it regulates apoptotic genes such as caspase 8, p38, Myc-associated factor X (MAX), and Ras protein-specific guanine nucleotide releasing factor 1 (RASGRF1) in neurons of mice [56] and that it suppresses macrophage-mediated inflammation in orbital fat-derived stem cells by upregulating IL-IRA and TNFRII expressions [57,58] highlighting its potential role in apoptosis and inflammation, both processes involved in OA. Target gene analysis of hsa-miR-140-3p, hsa-miR-33b-3p, and hsa-miR-671-3p revealed that 28 genes were cotargeted by at least two miRNAs. Among common targetgenes were insulin-like growth factor 1 receptor (IGF1R), insulin receptor (InsR), transforming growth factor beta receptor 1 (TGFβR1), cAMP responsive element binding protein 5 (CREB5), nuclear factor kB (NFkB), and tumor necrosis factor a (TNFa), genes enriched in inflammatory and metabolic processes as revealed by pathway analysis (Figs. 4 and 5). An interesting finding in the in silico prediction analysis was that InsR and IGFR1 were common targets of all three hsa-miR-140-3p, hsa-miR-33b-3p, and hsa-miR-671. Insulin is mediating its effects by binding to the specific high affinity insulin receptor (InsR), or with lower affinity to the structural-and functionally related insulin-like growth factor receptor (IGFR) [59]. Taking into consideration that OA is a metabolic disease [6], that insulin resistance plays a key role in the metabolic syndrome [3], and that recently a correlation was demonstrated between radiographic OA severity and insulin resistance [7], our in silico prediction analysis adds to the above evidence, suggesting a possible role of insulin resistance in OA. Along the same thought process, another interesting finding in our study was that cAMP responsive element binding protein 5 (CREB5), co-targeted by hsa-miR-140-3p and hsa-miR-671-3p, was significantly enriched in insulin secretion and TNF signaling pathways. CREB5 is a primary regulator of adipogenesis, that is also involved in regulating innate immune system, agingassociated inflammation and glucose homeostasis [60,61] and under obese conditions promotes insulin resistance by activating the transcriptional repressor ATF3 and by downregulating the expression of adiponectin as well as insulin-sensitive glucose transporter 4 (GLUT4) [62].
All above suggest that the in silico predictions in the present study are highlighting the possible implication of circulating hsa-miR-140-3p, hsa-miR-33b-3p, and hsa-miR-671-3p as novel serum-based biomarkers for Fig. 6 Diagram of relative miRNA expression levels by qRT-PCR in normal (n = 7) and osteoarthritic (OA) chondrocytes (n = 12). U6 was used for normalization of the real-time PCR data. All miRNAs appeared to be downregulated with respect to control samples. Fold change has been calculated as the log2 transformed ratio of the 2 −ΔΔCt of OA samples over control samples (asterisks depict a significant difference between the 2 −ΔΔCt of OA samples and 2 −ΔΔCt of control samples in cartilage tissues at the p < 0.05 level). Hsa-miR-140-3p, hsa-miR-671-3p, and hsa-miR-150p were significantly downregulated in serum samples of OA patients compared to healthy individuals osteoarthritis prognosis and underlining their involvement in the regulation of metabolic processes that contribute to OA pathology. However, functional studies need to be conducted to verify our bioinformatics analysis and to shed light on the metabolic pathways through which the above miRNA signature leads to OA onset and development. In addition, standardized protocols for circulating miRNA isolation and analysis in larger cohorts of different ethnicities must be established before their use in clinical practice.

Conclusions
A serum miRNA signature was established, for the first time, using high-density resolution miR-arrays in OA patients, and their deregulation was verified in articular cartilage samples. We identified a three-miRNA signature, hsa-miR-140-3p, hsa-miR-671-3p, and hsa-miR-33b-3p, in serum predicted to regulate metabolic processes, which could serve as a potential biomarker for the evaluation of OA risk and progression. HCA: Hierarchical clustering analysis; HCL: Hierarchical clustering; HDL: Highdensity lipoprotein; IGF1R: Insulin-like growth factor 1 receptor; IGFR: Insulinlike growth factor receptor; InsR: High affinity insulin receptor; MAX: Mycassociated factor X; miRNAs: MicroRNAs; NFkB: Nuclear Factor kB; OA: Osteoarthritis; qRT-PCR: Quantitative real-time PCR; RASGRF1: Ras proteinspecific guanine nucleotide releasing factor 1; ROC: Receiver Operating Characteristic curves; TGFβR1: Transforming growth factor beta receptor 1; TNFa: Tumor necrosis factor a Acknowledgements Authors would like to thank the patients and healthy controls for consenting to give material for this work.

Funding
The project was funded by a fellowship of Excellence IKY/SIEMENS given to E. Ntoumou.

Availability of data and materials
The majority of data generated or analyzed during this study are included in this manuscript (and its supplementary information files). Additional datasets (micro-RNA array raw data) used and/or analyzed during the current study are available upon reasonable request from the corresponding author atsezou@med.uth.gr. The electronic databases used for pathway/target gene analysis are as follows: TargetScan http://www.targetscan.org. Accessed 18 June 2017. LG was involved in data acquisition, bioinformatics, and statistical analysis. BM was involved in microarray screening, NE in experimental and bioinformatics analysis and drafting the manuscript, PM in acquisition and interpretation of data, TM in article conception and results interpretation and TA was involved in article conception, experimental design, interpretation, and revisiting critically the article for important intellectual content and approved the final version of the manuscript. All authors read and approved the final manuscript.
Ethics approval and consent to participate Consent was obtained from each participant. The study protocol conformed to the ethical guidelines of the 1975 Declaration of Helsinki as reflected in a priori approval by the local ethical committee of the University Hospital of Larissa.

Competing interests
The Authors declare that they have no competing interests.