Association of circulating microRNAs with prevalent and incident knee osteoarthritis in women: the OFELY study

Objectives In the context of the scarcity of biomarkers for knee osteoarthritis (OA), we examined the associations of prevalent and incident OA with the expression levels of serum miRNAs in subjects with and without OA. Methods With a next-generation sequencing approach, we compared the miRome expression of 10 women with knee OA and 10 age-matched healthy subjects. By real-time qPCR, we analyzed the expression levels of 19 miRNAs at baseline selecting 43 women with prevalent knee OA (Kellgren Lawrence score of 2/3), 23 women with incident knee OA over a 4-year follow-up and 67 healthy subjects without prevalent or incident OA matched for age and body mass index. Results Serum miR-146a-5p was significantly increased in the group of prevalent knee OA compared with controls (relative quantification (RQ); median [Interquartile range] 1.12 [0.73; 1.46] vs 0.85 [0.62; 1.03], p = 0.015). The likelihood of prevalent knee OA was significantly increased (odds ratio [95% confidence interval (CI)] 1.83 [1.21–2.77], p = 0.004) for each quartile increase in serum miR-146a-5p. The women with miR-146a-5p levels above the median (0.851) had a higher risk of prevalent knee OA compared to those below the median [95% CI] 4.62 [1.85–11.5], p = 0.001. Moreover, we found a significant association between the baseline level of serum miR-186-5p and the risk of incident knee OA (Q4 vs Q1–3; odds ratio [95% CI] 6.13 [1.14–32.9], p = 0.034). Conclusion We showed for the first time that miR-146a-5p and miR-186-5p are significantly associated with prevalent and incident knee OA, respectively.


Introduction
Osteoarthritis (OA) is the most frequent chronic musculoskeletal disease affecting approximately 40% of adults aged 70 years and over [1]. It is considered as a slowly progressive disease degrading all tissues of the affected joint, appearing as a destruction of the cartilage, the hallmark of OA, but also a mild-to-moderate synovial inflammation and an alteration of the subchondral bone structure. So far, the gold standard to estimate the extent of the disease has been plain radiography, but the poor sensitivity and the relatively large precision error lead to disease detection only when significant cartilage degradation has already occurred. Measuring circulating biomarkers appears as a potential non-invasive approach to diagnose the disease and prognosticate its evolution. Research has focused on the structural components of the cartilage matrix, especially type II collagen degradation markers [2]. Because of the large overlap between biomarker concentrations in controls and OA patients, the diagnosis of OA in an individual patient cannot be made by currently available biomarkers. Given the limitations of the tools that are currently available for OA assessment, there is considerable interest in the identification of specific biological markers that reflect quantitative and dynamic variations in joint tissue remodeling. For this purpose, microRNAs (miRNAs) are new potential targets. These small non-coding RNAs of 22-28 nucleotides in length can silence gene expression by binding to target messenger RNA (mRNA) repressing the translation. This regulation process is redundant because 1 miRNA can inhibit a large number of mRNA, and 1 gene can be targeted by multiple miRNAs. The miRNAs are important players in maintaining the health of all joint tissues including the cartilage, ligament, tendon, and muscle. Some of them have a potential role in regulating homeostasis of joint tissues in the context of OA [3]. Several in vitro and in vivo studies have reported the miRNA involvement in the OA onset and progression, by targeting cartilage-associated genes [4]. They regulate the expression of genes involved in pathways altered in OA chondrocytes such as apoptosis [5], expression levels of MMPs and ADAMTS [6][7][8], and chondrocyte signaling [9,10]. To date, it has been reported that about 80 miRNAs are involved in the pathology of OA and that their circulating level in patients might reflect the underlying disease state [4]. Indeed, circulating miR-NAs are considered to be attractive molecules to serve as biomarkers because they are easily accessible, stable in biofluids, and quantifiable with a high degree of sensitivity and specificity. Circulating miRNAs are protected from degradation because they circulate bound to proteins such as Argonaute 2 and nucleoplasmin or inside microvesicles such as exosomes [11][12][13]. Finally, although the expression of miRNAs has been well documented in OA cartilage, it remains difficult to understand to which extent circulating miRNAs reflect the tissue alteration. Moreover, the aberrant expression levels previously reported in circulating miRNAs differ considerably between clinical studies due to the different characteristics of these settings [14][15][16][17][18][19][20][21][22].
In this study, a first screening phase was conducted by comparing the serum miRome in subjects with (n = 10) and without (n = 10) knee OA in the OFELY cohort by next-generation sequencing (NGS) analysis. Then, on the basis of our NGS data and on the previously reported function of miRNAs in the OA process, both in basic and/or clinical research, we have selected 19 candidate miRNAs for quantitative real-time polymerase chain reaction (RT-qPCR) analysis. By conducting the validation phase of these candidate miRNAs in the OFELY cohort, we report that miR-146a-5p and miR-186-5p are significantly associated with prevalent and incident knee OA, respectively.

Patients selected from the OFELY cohort
The study group for both steps, NGS and RT-qPCR, included French women belonging to a population-based cohort. These women were participants in a prospective investigation of the determinants of bone loss, the Os des FEmmes de LYon (OFELY) study. This cohort has previously been described in detail elsewhere [23,24]. In the present analysis, we studied the miRNA expression in women with knee OA at baseline for this analysis, 8 years after the recruitment of the cohort (ninth followup visit). The primary or post-traumatic etiology of knee OA was not known, but the selected women had no other skeletal diseases. The evaluation of OA disease was performed by radiography for knee OA and spine disc degeneration, by clinical examination for hand OA, and by a questionnaire for the hip OA (for details on the OA evaluation, see [8]). After baseline OA assessment, women were followed prospectively for 4 years.

Assessment of knee OA
Radiographs of both knees were obtained in all women. Radiologic evaluation of the knees consisted of bilateral posteroanterior weight-bearing knee radiographs with fixed flexion using the SynaFlex X-ray Positioning Frame (Synarc, San Francisco, CA), as previously described [25]. Radiographs were obtained in a single radiography unit by the same staff of two technicians using a previously described standardized protocol [26]. The severity of OA was performed and graded according to Kellgren and Lawrence (KL) classification at baseline and 4 years later [27]. Prevalent knee OA was defined by a KL score higher or equal to 2 at baseline. Incident OA was defined by a KL score higher or equal to 2 at year 4 and a KL score < 2 at baseline. All knee radiographs were scored by a single trained rheumatologist (ES-R). Measurements were made paired but not blinded to sequence, which has been shown not to modify the sensitivity to change [28].
The miRNA analysis by next-generation sequencing Selection of patients The expression levels of miRNAs in the serum were measured in ten women with knee OA (KL score of 2 and 3) and in ten healthy women without OA at any site (for details, see [24]). Both groups were matched for age (healthy 61.9 ± 3.03 years and OA 63.9 ± 3.4 years, p = 0.17). The flow chart for the design of the experimental study is reported in Fig. 1.

NGS analysis
Total RNA was extracted from 400 μl of serum with the miRCURY Biofluids extraction kit (Exiqon®, Denmark) and analyzed by small RNA sequencing on a NextSeq500 sequencing instrument (Illumina platform). After RNA conversion into miRNA NGS libraries using the NEBNEXT library generation kit, cDNA was pre-amplified prior to library purification and quantification (for detailed procedures of extraction, library preparation and purification, normalization and quality controls, see Additional file 1: Data S1). Measurements were expressed as tags per million (TPM) in which the number of reads for a particular miRNA is divided by the total number of mapped reads in a sample and multiplied by 10 6 . The differential expression of miRNA between the OA and control groups was expressed by the TPM ratio converted as the Log 2 (fold change) (Log 2 FC).
The miRNA analysis by real-time quantitative polymerase chain reaction Selection of patients All women were postmenopausal and belonged to the OFELY cohort. The expression levels of selected miR-NAs were measured, for the prevalent OA, in 43 women with prevalent knee OA (KL score of 2 and 3; early and intermediate knee OA, age 68.3 ± 6.6 years, body mass index (BMI) 26.6 ± 4.4 kg/m 2 ) and in 42 healthy women without knee OA matched for age and BMI, and for the incident OA, in 23 women with incident knee OA over the next 4 years (age 68.4 ± 8 years, BMI 25.2 ± 4 kg/m 2 ) and in 25 healthy subjects without incident knee OA matched for age and BMI ( Table 1).

Preselection of miRNAs
The first selection of candidate miRNAs was based on the data obtained from the initial screening phase. We have selected miR-139-5p, miR-200a-3p, and miR-1299 in regard to their significant differential expression between the groups after Benjamini-Hochberg false discovery rate (FDR) correction. We have also added miR-375 for its high level of expression and a tendency to be differentially expressed, despite the absence of significant difference after FDR correction. Then, a second miRNA selection was performed by a critical assessment of literature data until December 2016 to highlight miRNAs that are differentially expressed in the biofluids of patients with knee OA vs controls. We have considered the criteria for patient selection, the size of the groups studied, the technical approaches used for screening and validation of miRNAs in biofluids, and the involvement of these differentially expressed miRNAs in the regulation of metabolic process contributing to OA pathology (see Table 2). Thus, we have further selected 15 candidate miRNAs: let-7e-5p, 16-5p, 29a-3p, 29b-3p, 29c-3p, 93-5p, 126-3p, 132-3p, 146a-5p, 184, 186-5p, 195-5p, 199a-3p, 345-5p, and 885-5p.

RNA isolation and miRNA RT-qPCR analysis
Total RNA was extracted from 200 μl of serum (miR-CURY RNA isolation kit for biofluids, EXIQON, Danemark) by adding a lysis solution containing cel-miR-39-3p as exogenous control, followed by protein precipitation and miRNA purification on a silica column. miRNAs quantification was performed by the TaqMan® Advanced miRNA technology (Applied Biosystems, Thermo Fisher Scientific) in which miRNA was reverse-transcribed and cDNA pre-amplified with the TaqMan Advanced miRNA-cDNA synthesis kit, followed by analysis on pre-designed miRNA TaqMan arrays (TLDA) according to the manufacturer's protocols (Applied Biosystems, CA, USA). Amplified cDNAs (15 μl) were mixed with 75 μl of TaqMan Fast Advanced Mastermix and water (60 μl), and 100 μl of the sample was added to each array tank of the card. The TLDA wells were pre-spotted with TaqMan®Advanced miRNA assays that allow the quantification in duplicate of 19 miRNAs (Table 3) by RT-qPCR reaction on a QuantStudio7 flex (Applied Biosystems) according to the manufacturer's protocol. The Ct values recorded by the software Expression Suite were normalized by the Ct mean of 3 endogenous miRNAs for the relative quantification (RQ) of miRNA levels as fold change (FC) = Log 2 (2 −ΔΔCT ) (for protocol details, see Additional file 1: Data S1).

Statistical analysis
A miRNA with a p value < 0.05 and a false discovery rate of 5% [Benjamini-Hochberg false discovery rate (FDR) correction for NGS approach] were considered as differentially expressed. Wilcoxon tests were used to compare the miRNA levels between women with and without OA because of the skewness of the data. We have examined the likelihood of OA, expressed as odds ratios (ORs) and 95% confidence intervals (CIs), per quartile increase in miRNA levels in a logistic regression model. All statistical analyses were performed using Stata 12 (Stata Corp LP, College Station, TX, USA).
Validation: differential expression of candidate miRNAs in serum of patients with knee OA and controls Prevalent knee OA and miRNA-146a-5p When considered as a continuous variable, miR-146a-5p was significantly increased in the group of prevalent knee OA compared with controls (relative quantification ( Incident knee OA and miRNA-186-5p When considered as a continuous variable, serum miR-186-5p levels were not significantly associated with an increased risk of incident knee OA (p = 0.09). However, in an analysis across quartiles, we found a significant association between baseline miR-186-5p levels and the likelihood of incident knee OA for each quartile increase (odds ratio [95% CI] 1.71 [1.00-2.95], p = 0.049).   Women in the upper quartile were 6.13 times more likely to develop incident knee OA over the next 4 years compared to those in the three others quartiles (Q4 vs Q1-3; odds ratio [95% CI] 6.13 [1.14-32.9], p = 0.034) ( Table 5 and Fig. 2b).

Discussion
In the present work, we have studied the serum miRNA profiling of women with knee OA to identify potential non-invasive biomarkers of the disease. We first screened the serum miRNAs with significant differences in the expression between women with knee OA and healthy controls by a NGS approach. Then, we have conducted a validation study in the OFELY cohort of women having an incident or prevalent knee OA with a long-time follow-up. In this, we analyzed by RT-qPCR whether the 3 differentially expressed miRNAs (139-5p, 200a-3p, 1299) in patients with knee OA in the screening step were also aberrantly expressed in the women with incident and prevalent knee OA of the validation cohort. In addition, miR-375 has been chosen because of its high expression level and a tendency to be differentially expressed in our NGS approach despite the absence of significant difference after FDR correction. Moreover, we have measured the expression levels of 15 additional miRNAs previously highlighted as having a dysregulated expression in patients with knee OA. The miRs-29 family (a, b, c) and miR-199a-3p were chosen for their activity as key regulators of chondrocyte gene expression with aberrant expression in OA cartilage [29][30][31]. The miRNAs let-7e-5p, 16-5p, 93-5p, 126-3p, 132-3p, 146a-5p, 184, 186-5p, 195-5p, 345-5p, and 885-5p were selected because previous studies analyzing circulating miRNAs in the serum and plasma [14,17,18] have reported their dysregulated expression in patients with knee OA (see Table 4).
The NGS approach revealed that miR-139-5p, miR-1299, and miR-200a-3p had levels of expression significantly different between OA patients and controls. However, we did not validate these 3 miRNAs in the second step of the analysis when we measured their level in the largest validation cohort. This lack of validation might  [21]) and/or differences in the specificity and sensibility of the platforms used (Illumina for this study; Borgonio-Cuadra [17] and Beyer [18]: Applied Biosystems; Kong [21]: Affimetrix). It appears that the previous studies performed with biological fluids are difficult to compare each other because of the heterogeneity of the clinical situations (see Table 2). Nevertheless, our study can be compared to that of Borgonio-Cuadra et al. [17], particularly the patient selection is based on the same criteria, radiologic KL scores of 2 and 3 and a BMI smaller than 27. However, among the ten miRNAs in common, only miR-146a-5p (miR-146) showed significant overexpression, emphasizing the challenges faced with patient heterogeneity. The reasons for this discrepancy could be the ethnic origin, Mexican vs European, enhancing the differences in exposure to environmental factors; the sex of the participants, men and women vs women only; the biological fluid used-plasma vs serum; or the mean age of participants, 55 vs 68 years.
Nevertheless, we found that miR-146a-5p expression is significantly increased in the OA group of postmenopausal women compared to healthy controls independently of age and BMI. Moreover, we reported for the first time that the risk of prevalent knee OA was significantly increased for each quartile increase in the serum level of miR-146a-5p. However, the interpretation of this increase remains challenging in regard to the dual role of miR-146a-5p in the pathophysiology of osteoarthritis. MiR-146 is expressed in chondrocytes which begin to undergo degenerative changes [32]; its expression is stimulated by pro-inflammatory cytokines in an NF-κBdependent manner [33]. It may contribute to OA by impairing the TGF-β signaling pathway targeting Smad4 and increasing apoptosis [5,34]. Recently, Zhang et al. showed that miR-146 facilitates OA progression by targeting Camk2d and Ppp3r2 that are required to maintain the phenotype of mature chondrocytes [33].
In contrast, miR-146 could play a protective role in OA. It is induced during hypoxia by HIF-1 and promotes When we compared the OA vs healthy groups, 22 miRNAs showed differential expression (p < 0.05) between controls and OA patients. After Benjamini-Hochberg false discovery rate (FDR) correction, hsa-miR-139-5p, hsa-miR-1299, and hsa-miR-200a-3p remained significantly different between OA patients and controls (p < 0.05, FDR at 5%) autophagy by decreasing the Bcl-2 expression, an autophagy inhibitor. Autophagy has a protective function for cartilage during OA [35]. Moreover, miR-146 exerts negative control on inflammatory responses by suppressing IRAK1 and TRAF6, two molecular targets of miR-146 impairing NF-kβ activity and suppressing NF-kβ target genes such as IL-6 and TNF-α [36,37]. Recently, Zhong et al. [38] showed that miR-146 could increase the proliferation and inhibit apoptosis of OA chondrocytes by inhibiting TRAF6 expression and suppressing the activation of the NF-κB signaling pathway. Finally, Guan et al. [39] showed that miR-146 is activated by cyclic loads in the physiological range but suppressed by mechanical overload in human articular chondrocytes. They proposed that miR-146 protects joint cartilage from degeneration through inhibition of the feedback loop of Notch-1/IL-6 during aging or trauma. It must be noted that circulating miR-146 has been highlighted in two other studies analyzing the peripheral blood mononuclear cells. Okuhara et al. [16] found it overexpressed in patients suffering from knee OA compared to healthy controls and highly expressed in OA peripheral blood mononuclear cells from patients with early OA. Soyocak et al. [20] observed that, in patients classified according to the KL score, miR146a expression increased in the progressive stages (grades 3 and 4).
Finally, miR-146 is one of the four miRNAs with miR-223, miR-16, and miR-30b that were differently expressed in the cartilage and in the circulation with a role in the cartilage homeostasis [4]. Moreover, miR-146 is also involved in the inhibition of osteoclastogenesis [40] suggesting that its measurement in the serum could reflect the deregulations affecting the cartilage and bone tissues. Although the activity of this miRNA seems dependent on the disease phases, it appears as a key element in the OA pathophysiology, and our results suggest that miR-146a-5p could be a new biological marker for knee OA.
We report for the first time a significant association, independent of age and BMI, of serum miR-186-5p and incident knee OA risk in a population of postmenopausal women followed prospectively for 4 years. This result is in accordance with Borgonio-Cuadra et al. [17] showing that miR-186-5p was overexpressed in the plasma of early-stage OA patients compared to controls. In silico analysis has suggested that potential pathways regulated by miR-186 could be signaling by PDGF, developmental biology, membrane trafficking, and collagen formation. However, Li et al. [19] showed that miR-186 originated from the synovial membrane and belongs to a panel of seven miRNAs that were significantly upregulated in late-stage OA synovial fluid compared to early- The miRNA levels were analyzed in the serum from individuals of the validation cohort, and differences in the serum miRNA levels are reported according to the prevalent or incident OA status. The relative miRNA level was used for the statistical comparisons. IQR = interquartile range [(+) = women with prevalent or incident knee OA, (−) = controls] stage irrespective of age, gender, and BMI. The reasons for this discrepancy are not clear but could be the geographic origin of patients, North America vs Europe, or the biological fluids analyzed, synovial fluid vs serum/ plasma. Moreover, patients from Li's study were selected from a group with degenerative meniscal tears undergoing arthroscopic surgery (early OA) or total knee replacement (late-stage OA) while Borgonio-Cuadra et al. [17] and us have selected patients on the radiographic score only, without any consideration regarding the OA etiology.
Together with the serum association of miR-146a-5p and prevalent OA, these data reinforce the notion that the miRNAs involved in OA pathophysiology vary according to the stage of the disease [19]. Finally, miR-186 has been essentially studied in cancer, Alzheimer's disease, neuropathic disorder, and acute myocardial infarction [41]. It can be noticed that the last 3 diseases have been shown to be linked to OA. Recently, Weber et al. [42] demonstrated an association between OA and the risk of dementia. The miR-186 belonged to a set of 27 differentially expressed miRNAs between patients suffering from Alzheimer's disease and controls [43]. The pathways potentially downregulated are linked with neuronal synaptic functions. Moreover, neuropathic pain is at least partially responsible for pain in OA [44]. Interestingly, miR-186 negatively regulated CXCL13, also known as B lymphocyte chemoattractant, and its downregulation in spinal neurons after spinal nerve ligation causes upregulation of CLCX13 to drive neuropathic pain [45]. Finally, several cross-sectional studies have shown a strong association between OA and cardiovascular diseases [46], and the overexpression of miR-186 in plasma have the potential to be used for the early detection of myocardial infarction [47]. However, for the present analysis, no woman suffered from dementia. Moreover, we had incomplete information regarding cardiovascular events, and we did not have data concerning pain in the OFELY study. So, the potential links between miR-186-5p, cardiovascular diseases, dementia, and OA could not be evaluated.
Our study has strengths and limitations. We investigated the association between the serum miRNAs expression and prevalence and for the first time incidence of knee OA in a well-characterized population of women followed prospectively over a long period of time of 4 years. We have chosen to carry out the validation step of our study by analysis of duplicate samples on TaqMan microRNA arrays, in order to reduce erratic manipulations and improve the reproducibility between the samples. Moreover, our study combines the validation of our NGS results and the first attempt to replicate previously published results, a process that Fig. 2 Risk of prevalent and incident OA according to the median of miR-146a-5p and the quartiles of miR-186-5p, respectively. Multiple logistic regression analyses to determine the ability of serum miR-146a-5p and miR-186-5p levels to predict prevalent or incident OA, respectively. a Prevalent OA: percentage of women with osteoarthritis under and upper the median of miR-146a-5p (number of women with prevalent OA − total number of women). b Incident OA: percentage of women with osteoarthritis in the first three quartiles vs the upper quartile (number of women with incident OA − total number of women). CI, confidence interval; OD, odds ratio had been lacking in the miRNA research field. Finally, we have indicated that women selected for the analyses of prevalent and incident knee OA had also OA at other sites (Table 1). This is a key point because we showed in the OFELY cohort that at this mean age, 68 ± 8 years, the prevalence of OA was high with 75% and 88% of women having osteophytosis at the lumbar spine and thoracic spine, respectively [48,49] in agreement with previous prevalence data [50,51]. Consequently, in both groups, knee OA and healthy women, we have verified that the presence of OA at other sites did not cause any bias between the groups (see Table 1) in order to decrease as much as possible the effects of this unavoidable situation of generalized OA on the detection of miRNAs associated with knee OA. However, we studied knee OA in women only, and the number of participants is limited. It has been demonstrated that gender may influence the level of the expression of miRNAs [52]. Consequently, our results need to be confirmed by additional larger studies and particularly in males. It can be noticed that miR-146 and miR-186 have been quantified in serum patients because the miR-3p forms of both miRNAs were almost undetectable in the 20 tested sera in the screening phase (data not shown). This lack of detection in the serum is probably related to thermodynamic instability of the 3p arms and to their subsequent degradation in the cellular p-bodies. Finally, the molecular mechanisms associating serum miR-146 and miR-186 to OA remain to be further clarified, for example, the role of the singlenucleotide polymorphisms in the promoter region of miR-146 gene [53].
Collectively, our results show that miR-146a-5p is increased in women with mild to moderate knee OA compared to healthy women. Importantly, miR-186-5p is also increased in those women who will develop radiographic knee OA over the next 4 years; therefore, this miRNA has the potential to detect preclinical knee OA.