- Open Access
Serum microRNA array analysis identifies miR-140-3p, miR-33b-3p and miR-671-3p as potential osteoarthritis biomarkers involved in metabolic processes
Clinical Epigeneticsvolume 9, Article number: 127 (2017)
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.
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 . OA is defined as a heterogeneous disease consisting of a group of distinct joint disorders with common biological, morphological, and clinical outcomes . 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 , 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 . Therefore, the identification of non-invasive 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 . Recent studies have implicated miRNAs 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 .
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 . Several microRNA 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 .
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  have implicated them as attractive molecules for disease prognosis and/or progression . Circulating miRNAs have been demonstrated as important biomarkers in many diseases, as diabetes type II  metabolic syndrome  Alzheimer’s, hypertension, osteoporosis, Parkinson’s disease, and cancer . 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.
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 post-traumatic 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 (miRBase 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 220.127.116.11 (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.
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 .
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 penicillin-streptomycin, and were incubated at 37 °C under a humidified 5% CO2 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) stem-loop RT primer to generate the cDNA as previously described . 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 .
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 . 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 . 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  and DIANA Tools, a collection of web-tools for miRNA functional annotation analysis . Functional annotation was performed with the Webgestalt web-tool  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  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).
Diagnostic value of miRNAs in serum
Receiver Operating Characteristic (ROC) curves were established to evaluate the diagnostic value of deregulated miRNAs in discriminating between OA patients and healthy individuals. ROC analysis revealed that 77 miRNAs had an area under the curve (AUC) > 0.8 at a p < 0.05 significance level, that could thus separate OA from control samples. ROC analysis results are summarized in Additional file 1: Table S1. Among the 77 miRNAs, 7 miRNAs, hsa-miR-33b-3p, hsa-miR-4284, hsa-miR-663a, hsa-miR-150-5p, hsa-miR-1233-3p, hsa-miR-671-3p, and hsa-miR-140-3p, were selected for further validation based on the fact that they were in the top 20 up- or downregulated miRNAs as found in microarray microRNA expression profiling presented in Fig. 1a, b. More specifically, miR-33b-3p and miR-4284 were among the top upregulated, and miR-663a, miR150-5p, miR-1233-3p, miR-671-3p, and miR-140-3p were downregulated.
ROC curve analysis of the seven selected miRNAs is presented in Fig. 2 .
Validation of selected miRNAs expression in serum of OA and healthy individuals by qRT-PCR
We verified by qRT-PCR the expression levels of the seven selected miRNAs (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) in all serum samples. Among these seven miRNAs, 3 miRNAs, hsa-miR-33b-3p, hsa-miR-671-3p, and hsa-miR-140-3p, were found significantly downregulated in OA serum samples compared to controls (p < 0.05) (Fig. 3a ).
Moreover, we found that hsa-miR-33b-3p and hsa-miR-4284 expression levels coincided between microarray analysis and qRT-PCR, exhibiting decreased expression in serum of OA samples compared to controls. Hsa-miR-663a, hsa-miR-150-5p, hsa-miR-1233-3p, hsa-miR-140-3p, and hsa-miR-671-3p, which were marginally over-expressed in microarray analysis, appeared to be downregulated with qRT-PCR. This discrepancy might be due to technical limitations of the microarray methodology (Fig. 3b ).
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, ECM-receptor interaction, and fatty acid biosynthesis. Annotated pathways are summarized in Table 2.
Furthermore, by considering the qRT-PCR data, we focused on predicted target genes, by TargetScan, of the three miRNAs (hsa-miR-140-3p, hsa-miR-33b-3p, and hsa-miR-671-3p) that were found downregulated in the serum of OA patients. It appeared (Fig. 4 ) that these genes are predicted targets of two or more of the three miRNAs. More specifically, 28 genes were co-targeted by at least 2 miRNAs. Among these, INSR and IGF1R were targeted by all three miRNAs, while ADCY7, VANGL1, WNT5A, PPP3R2, TGFBR1, FZD4, BRAF, and CREB5 were co-targeted by hsa-miR-140-3p, and hsa-miR-671-3p, RPTOR, CDKN1A, PRKCB, ADCYA1, PDPK1, MAPK1, ADCY2, TCS1, Kras, CBL, CHRM3, SGL1, CREB1, KCNN3, PRKC1, PPP2R5E, and KCNJ11 were co-targeted by hsa-miR-33b-3p, and hsa-miR-140-3p and KCNN1 were co-targeted by hsa-miR-671-3p and hsa-miR-33b-3p.
Pathway enrichment analysis (Additional file 2: Table S2) revealed that the target genes of hsa-miR-140-3p, hsa-miR-33b-3p,and hsa-miR-671-3p are potentially involved in 48 common signaling pathways, including thyroid hormone synthesis, FoxO signaling pathway, insulin secretion, chemokine signaling pathway, MAPK signaling pathway, PI3K-Akt signaling pathway, estrogen signaling pathway, regulation of lipolysis in adipocytes, glucagon signaling pathway, and calcium signaling pathway, whereas 40 common pathways were predicted common between miR-140-3p and miR-33b-3p, such as mTOR signaling pathway, osteoclast differentiation, Ras signaling pathway, type II diabetes mellitus, adipocytokine signaling pathway, thyroid hormone signaling pathway, insulin signaling pathway, insulin resistance, sphingolipid signaling pathway, sphingolipid metabolism, and ErbB signaling pathway. Furthemore, five pathways were predicted between miR-140-3p and miR-671-3p, including Wnt signaling pathway, endocrine and other factor-regulated, and calcium reabsorption (Fig. 5 ).
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.
miRNA expression in OA and healthy articular cartilage samples
We next evaluated by qRT-PCR the expression levels of the 7 selected miRNAs (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) in 12 OA and 7 healthy articular cartilage samples. The expression levels of hsa-miR-140-3p, hsa-miR-150-5p, and hsa-miR-671-3p were significantly decreased in OA articular cartilage compared to healthy cartilage (p < 0.05) (Fig. 6 ). No significant differences were observed for hsa-miR-33-3p, hsa-miR-4284, hsa-miR-663a, and hsa-miR-1233-3p expression levels between OA and healthy articular cartilage.
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  or adipose tissue depots , 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 . 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 . A more recent study by Beyer et al.  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 miRNAs in OA could be attributed to differences in the study design, as selection of patients, disease status, and differences in ethnicity , 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 .
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 miRNAs 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 miRNAs 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 ).
We, next, validated with qRT-PCR, which offers high accuracy, sensitivity, and dynamic range  the expression levels of 7 selected out of the 77 DE miRNAs, which were among the top 20 up- or downregulated in the microarray screening, namely, 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, in serum and in OA and healthy articular cartilage samples. We found that three miRNAs, hsa-miR-140-3p, hsa-miR-671-3p, and hsa-miR-33b-3p, were significantly downregulated in the serum of OA patients compared to controls, and in addition, hsa-miR-140-3p and hsa-miR-671-3p were also significantly downregulated in OA compared to healthy articular cartilage (Figs. 3a, b and 6 ).
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 and that decreased levels of miR-140 were found in the plasma of patients with morbid obesity influencing TGFbR1 expression levels , 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 . Recently, miR-33b was shown to regulate adipocyte differentiation  and miR-33a different functions of macrophages having a role in development and progress of atherosclerosis . Regarding miR-671, besides its role in cancer cell migration , 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  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 co-targeted by at least two miRNAs. Among common target-genes 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) . Taking into consideration that OA is a metabolic disease , that insulin resistance plays a key role in the metabolic syndrome , and that recently a correlation was demonstrated between radiographic OA severity and insulin resistance , 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, aging-associated 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) .
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 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.
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.
Area under the curve
- CREB5 :
cAMP responsive element binding protein 5
Database for Annotation, Visualization and Integrated Discovery
Fetal bovine serum
False discovery rate
- GLUT4 :
Insulin-sensitive glucose transporter 4
Hierarchical clustering analysis
- IGF1R :
Insulin-like growth factor 1 receptor
- IGFR :
Insulin-like growth factor receptor
- InsR :
High affinity insulin receptor
Myc-associated factor X
- NFkB :
Nuclear Factor kB
Quantitative real-time PCR
- RASGRF1 :
Ras protein-specific guanine nucleotide releasing factor 1
Receiver Operating Characteristic curves
- TGFβR1 :
Transforming growth factor beta receptor 1
- TNFa :
Tumor necrosis factor a
Cross M, Smith E, Hoy D, Nolte S, Ackerman I, Fransen M, et al. The global burden of hip and knee osteoarthritis: estimates from the global burden of disease 2010 study. Ann Rheum Dis. 2014;73(7):1323–30. doi: 10.1136/annrheumdis-2013-204763.
Goldring MB, Goldring SR. Osteoarthritis. J Cell Physiol. 2007;213(3):626–34. doi: 10.1002/jcp.21258.
Loeser RF. Aging and osteoarthritis: the role of chondrocyte senescence and aging changes in the cartilage matrix. Osteoarthr Cartil. 2009;17(8):971–9. doi: 10.1016/j.joca.2009.03.002.
Felson DT. Clinical practice. Osteoarthritis of the knee. N Engl J Med. 2006;354(8):841–8. doi: 10.1056/NEJMcp051726.
Zhuo Q, Yang W, Chen J, Wang Y. Metabolic syndrome meets osteoarthritis. Nat Rev Rheumatol. 2012;8(12):729–37. doi: 10.1038/nrrheum.2012.135.
Kluzek S, Newton JL, Arden NK. Is osteoarthritis a metabolic disorder? Br Med Bull. 2015;115(1):111–21. doi: 10.1093/bmb/ldv028.
Tootsi K, Martson A, Kals J, Paapstel K, Zilmer M. Metabolic factors and oxidative stress in osteoarthritis: a case-control study. Scand J Clin Lab Invest. 2017:1–7. doi: 10.1080/00365513.2017.1354255.
Sharif B, Kopec J, Bansback N, Rahman MM, Flanagan WM, Wong H, et al. Projecting the direct cost burden of osteoarthritis in Canada using a microsimulation model. Osteoarthr Cartil. 2015;23(10):1654–63. doi: 10.1016/j.joca.2015.05.029.
Reichmann WM, Maillefert JF, Hunter DJ, Katz JN, Conaghan PG, Losina E. Responsiveness to change and reliability of measurement of radiographic joint space width in osteoarthritis of the knee: a systematic review. Osteoarthr Cartil. 2011;19(5):550–6. doi: 10.1016/j.joca.2011.01.023.
Gregory RI, Yan KP, Amuthan G, Chendrimada T, Doratotaj B, Cooch N, et al. The microprocessor complex mediates the genesis of microRNAs. Nature. 2004;432(7014):235–40. doi: 10.1038/nature03120.
Chen X, Ba Y, Ma L, Cai X, Yin Y, Wang K, et al. Characterization of microRNAs in serum: a novel class of biomarkers for diagnosis of cancer and other diseases. Cell Res. 2008;18(10):997–1006. doi: 10.1038/cr.2008.282.
Kumar S, Vijayan M, Bhatti JS, Reddy PH. MicroRNAs as peripheral biomarkers in aging and age-related diseases. Prog Mol Biol Transl Sci. 2017;146:47–94. doi: 10.1016/bs.pmbts.2016.12.013.
Miles GD, Seiler M, Rodriguez L, Rajagopal G, Bhanot G. Identifying microRNA/mRNA dysregulations in ovarian cancer. BMC research notes. 2012;5:164. doi: 10.1186/1756-0500-5-164.
Taurino C, Miller WH, McBride MW, McClure JD, Khanin R, Moreno MU, et al. Gene expression profiling in whole blood of patients with coronary artery disease. Clin Sci. 2010;119(8):335–43. doi: 10.1042/CS20100043.
van Meurs JB. Osteoarthritis year in review 2016: genetics, genomics and epigenetics. Osteoarthr Cartil. 2017;25(2):181–9. doi: 10.1016/j.joca.2016.11.011.
Sondag GR, Haqqi TM. The role of MicroRNAs and their targets in osteoarthritis. Curr Rheumatol Rep. 2016;18(8):56. doi: 10.1007/s11926-016-0604-x.
Iliopoulos D, Malizos KN, Oikonomou P, Tsezou A. Integrative microRNA and proteomic approaches identify novel osteoarthritis genes and their collaborative metabolic and inflammatory networks. PLoS One. 2008;3(11):e3740. doi: 10.1371/journal.pone.0003740.
Li YH, Tavallaee G, Tokar T, Nakamura A, Sundararajan K, Weston A, et al. Identification of synovial fluid microRNA signature in knee osteoarthritis: differentiating early- and late-stage knee osteoarthritis. Osteoarthr Cartil. 2016;24(9):1577–86. doi: 10.1016/j.joca.2016.04.019.
Aryal B, Singh AK, Rotllan N, Price N, Fernandez-Hernando C. MicroRNAs and lipid metabolism. Curr Opin Lipidol. 2017;28(3):273–80. doi: 10.1097/MOL.0000000000000420.
Chen X, Liang H, Zhang J, Zen K, Zhang CY. Secreted microRNAs: a new form of intercellular communication. Trends Cell Biol. 2012;22(3):125–32. https://doi.org/10.1016/j.tcb.2011.12.001.
Turchinovich A, Burwinkel B. Distinct AGO1 and AGO2 associated miRNA profiles in human cells and blood plasma. RNA Biol. 2012;9(8):1066–75. https://doi.org/10.4161/rna.21083.
Vickers KC, Palmisano BT, Shoucri BM, Shamburek RD, Remaley AT. MicroRNAs are transported in plasma and delivered to recipient cells by high-density lipoproteins. Nat Cell Biol. 2011;13(4):423–33. doi: 10.1038/ncb2210.
Zhang Y, Liu D, Chen X, Li J, Li L, Bian Z, et al. Secreted monocytic miR-150 enhances targeted endothelial cell migration. Mol Cell. 2010;39(1):133–44. doi: 10.1016/j.molcel.2010.06.010.
Gantier MP, McCoy CE, Rusinova I, Saulep D, Wang D, Xu D, et al. Analysis of microRNA turnover in mammalian cells following Dicer1 ablation. Nucleic Acids Res. 2011;39(13):5692–703. doi: 10.1093/nar/gkr148.
Guay C, Regazzi R. Circulating microRNAs as novel biomarkers for diabetes mellitus. Nat Rev Endocrinol. 2013;9(9):513–21. doi: 10.1038/nrendo.2013.86.
Parrizas M, Novials A. Circulating microRNAs as biomarkers for metabolic disease. Best Pract Res Clin Endocrinol Metab. 2016;30(5):591–601. doi: 10.1016/j.beem.2016.08.001.
Beyer C, Zampetaki A, Lin NY, Kleyer A, Perricone C, Iagnocco A, et al. Signature of circulating microRNAs in osteoarthritis. Ann Rheum Dis. 2015;74(3):e18. doi: 10.1136/annrheumdis-2013-204698.
Borgonio Cuadra VM, Gonzalez-Huerta NC, Romero-Cordoba S, Hidalgo-Miranda A, Miranda-Duarte A. Altered expression of circulating microRNA in plasma of patients with primary osteoarthritis and in silico analysis of their pathways. PLoS One. 2014;9(6):e97690. doi: 10.1371/journal.pone.0097690.
Andersen CL, Jensen JL, Orntoft TF. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004;64(15):5245–50. doi: 10.1158/0008-5472.CAN-04-0496.
Chen C, Ridzon DA, Broomer AJ, Zhou Z, Lee DH, Nguyen JT, et al. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucleic Acids Res. 2005;33(20):e179. doi: 10.1093/nar/gni178.
Zhang D, Zhang M, Wells MT. Multiplicative background correction for spotted microarrays to improve reproducibility. Genet Res. 2006;87(3):195–206. doi: 10.1017/S0016672306008196.
Causton H, Quackenbush J, Brazma A. Microarray gene expression data analysis : a beginner's guide. Malden: Blackwell Publishing; 2003.
Cleveland W. Robust locally weighted regression and smoothing scatterplots. J Am Stat Assoc. 1979;74(Dec.):829–36.
Bolstad BM, Irizarry RA, Astrand M, Speed TP. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics (Oxford, England). 2003;19(2):185–93.
Klipper-Aurbach Y, Wasserman M, Braunspiegel-Weintrob N, Borstein D, Peleg S, Assa S, et al. Mathematical formulae for the prediction of the residual beta cell function during the first two years of disease in children and adolescents with insulin-dependent diabetes mellitus. Med Hypotheses. 1995;45(5):486–90.
Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc Natl Acad Sci U S A. 2003;100(16):9440–5.
Storey JD, Tibshirani R. Statistical methods for identifying differentially expressed genes in DNA microarrays. Methods in molecular biology (Clifton, NJ). 2003;224:149–57.
Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. elife. 2015;4 doi: 10.7554/eLife.05005.
Paraskevopoulou MD, Georgakilas G, Kostoulas N, Vlachos IS, Vergoulis T, Reczko M, et al. DIANA-microT web server v5.0: service integration into miRNA functional analysis workflows. Nucleic Acids Res. 2013;41(Web Server issue):W169–73. doi: 10.1093/nar/gkt393.
Zhang B, Kirov S, Snoddy J. WebGestalt: an integrated system for exploring gene sets in various biological contexts. Nucleic Acids Res. 2005;33(Web Server issue):W741–8. doi: 10.1093/nar/gki475.
Dennis G Jr, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, et al. DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. 2003;4(5):P3.
Laterza OF, Lim L, Garrett-Engele PW, Vlasakova K, Muniappa N, Tanaka WK, et al. Plasma MicroRNAs as sensitive and specific biomarkers of tissue injury. Clin Chem. 2009;55(11):1977–83. doi: 10.1373/clinchem.2009.131797.
Thomou T, Mori MA, Dreyfuss JM, Konishi M, Sakaguchi M, Wolfrum C, et al. Adipose-derived circulating miRNAs regulate gene expression in other tissues. Nature. 2017;542(7642):450–5. doi: 10.1038/nature21365.
Okuhara A, Nakasa T, Shibuya H, Niimoto T, Adachi N, Deie M, et al. Changes in microRNA expression in peripheral mononuclear cells according to the progression of osteoarthritis. Mod Rheumatol. 2012;22(3):446–57. doi: 10.1007/s10165-011-0536-2.
Wang X, Sundquist J, Zoller B, Memon AA, Palmer K, Sundquist K, et al. Determination of 14 circulating microRNAs in Swedes and Iraqis with and without diabetes mellitus type 2. PLoS One. 2014;9(1):e86792. doi: 10.1371/journal.pone.0086792.
Brazma A, Hingamp P, Quackenbush J, Sherlock G, Spellman P, Stoeckert C, et al. Minimum information about a microarray experiment (MIAME)-toward standards for microarray data. Nat Genet. 2001;29(4):365–71. doi: 10.1038/ng1201-365.
Nicolas FE, Pais H, Schwach F, Lindow M, Kauppinen S, Moulton V, et al. Experimental identification of microRNA-140 targets by silencing and overexpressing miR-140. RNA. 2008;14(12):2513–20. doi: 10.1261/rna.1221108.
Miyaki S, Sato T, Inoue A, Otsuki S, Ito Y, Yokoyama S, et al. MicroRNA-140 plays dual roles in both cartilage development and homeostasis. Genes Dev. 2010;24(11):1173–85. doi: 10.1101/gad.1915510.
Karlsen TA, Jakobsen RB, Mikkelsen TS, Brinchmann JE. microRNA-140 targets RALA and regulates chondrogenic differentiation of human mesenchymal stem cells by translational enhancement of SOX9 and ACAN. Stem Cells Dev. 2014;23(3):290–304. doi: 10.1089/scd.2013.0209.
Tardif G, Pelletier JP, Fahmi H, Hum D, Zhang Y, Kapoor M, et al. NFAT3 and TGF-beta/SMAD3 regulate the expression of miR-140 in osteoarthritis. Arthritis research & therapy. 2013;15(6):R197. doi: 10.1186/ar4387.
Gernapudi R, Wolfson B, Zhang Y, Yao Y, Yang P, Asahara H, et al. MicroRNA 140 promotes expression of long noncoding RNA NEAT1 in Adipogenesis. Mol Cell Biol. 2016;36(1):30–8. doi: 10.1128/MCB.00702-15.
Hochberg MC, Tracy JK, Hawkins-Holt M, Flores RH. Comparison of the efficacy of the tumour necrosis factor alpha blocking agents adalimumab, etanercept, and infliximab when added to methotrexate in patients with active rheumatoid arthritis. Ann Rheum Dis. 2003;62(Suppl 2):ii13–6.
Kostopoulou F, Malizos KN, Papathanasiou I, Tsezou A. MicroRNA-33a regulates cholesterol synthesis and cholesterol efflux-related genes in osteoarthritic chondrocytes. Arthritis research & therapy. 2015;17:42. doi: 10.1186/s13075-015-0556-y.
Kuhn BM, Nodzynski T, Errafi S, Bucher R, Gupta S, Aryal B, et al. Flavonol-induced changes in PIN2 polarity and auxin transport in the Arabidopsis thaliana rol1-2 mutant require phosphatase activity. Sci Rep. 2017;7:41906. doi: 10.1038/srep41906.
Dlouha D, Hubacek JA. Regulatory RNAs and cardiovascular disease—with a special focus on circulating microRNAs. Physiol Res. 2017;66(Supplementum 1):S21–38.
Nan A, Chen L, Zhang N, Liu Z, Yang T, Wang Z, et al. A novel regulatory network among LncRpa, CircRar1, MiR-671 and apoptotic genes promotes lead-induced neuronal cell apoptosis. Arch Toxicol. 2017;91(4):1671–84. doi: 10.1007/s00204-016-1837-1.
Tan X, Fu Y, Chen L, Lee W, Lai Y, Rezaei K, et al. miR-671-5p inhibits epithelial-to-mesenchymal transition by downregulating FOXM1 expression in breast cancer. Oncotarget. 2016;7(1):293–307. 10.18632/oncotarget.6344.
Lien GS, Liu JF, Chien MH, Hsu WT, Chang TH, Ku CC, et al. The ability to suppress macrophage-mediated inflammation in orbital fat stem cells is controlled by miR-671-5p. Stem Cell Res Ther. 2014;5(4):97. doi: 10.1186/scrt486.
Li G, Barrett EJ, Wang H, Chai W, Liu Z. Insulin at physiological concentrations selectively activates insulin but not insulin-like growth factor I (IGF-I) or insulin/IGF-I hybrid receptors in endothelial cells. Endocrinology. 2005;146(11):4690–6. doi: 10.1210/en.2005-0505.
Mayr BM, Canettieri G, Montminy MR. Distinct effects of cAMP and mitogenic signals on CREB-binding protein recruitment impart specificity to target gene activation via CREB. Proc Natl Acad Sci U S A. 2001;98(19):10936–41. doi: 10.1073/pnas.191152098.
Mayr B, Montminy M. Transcriptional regulation by the phosphorylation-dependent factor CREB. Nat Rev Mol Cell Biol. 2001;2(8):599–609. doi: 10.1038/35085068.
Qi L, Saberi M, Zmuda E, Wang Y, Altarejos J, Zhang X, et al. Adipocyte CREB promotes insulin resistance in obesity. Cell Metab. 2009;9(3):277–86. doi: 10.1016/j.cmet.2009.01.006.
Authors would like to thank the patients and healthy controls for consenting to give material for this work.
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 firstname.lastname@example.org.
The electronic databases used for pathway/target gene analysis are as follows:
TargetScan http://www.targetscan.org. Accessed 18 June 2017.
Webgestalt web-tool http://www.webgestalt.org/. Accessed 18 June 2017.
DIANA tools http://diana.imis.athena-innovation.gr/DianaTools/. Accessed 18 June 2017.
DAVID https://david.ncifcrf.gov/. Accessed 18 June 2017.
Venn Diagram Generator http://bioinformatics.psb.ugent.be/webtools/Venn/. Accessed 18 June 2017.
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.
Consent for publication
The Authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.