Equine synovial fluid small non-coding RNA signatures in early osteoarthritis

Background Osteoarthritis remains one of the greatest causes of morbidity and mortality in the equine population. The inability to detect pre-clinical changes in osteoarthritis has been a significant impediment to the development of effective therapies against this disease. Synovial fluid represents a potential source of disease-specific small non-coding RNAs (sncRNAs) that could aid in the understanding of the pathogenesis of osteoarthritis. We hypothesised that early stages of osteoarthritis would alter the expression of sncRNAs, facilitating the understanding of the underlying pathogenesis and potentially provide early biomarkers. Methods Small RNA sequencing was performed using synovial fluid from the metacarpophalangeal joints of both control and early osteoarthritic horses. A group of differentially expressed sncRNAs was selected for further validation through qRT-PCR using an independent cohort of synovial fluid samples from control and early osteoarthritic horses. Bioinformatic analysis was performed in order to identify putative targets of the differentially expressed microRNAs and to explore potential associations with specific biological processes. Results Results revealed 22 differentially expressed sncRNAs including 13 microRNAs; miR-10a, miR-223, let7a, miR-99a, miR-23b, miR-378, miR-143 (and six novel microRNAs), four small nuclear RNAs; U2, U5, U11, U12, three small nucleolar RNAs; U13, snoR38, snord96, and one small cajal body-specific RNA; scarna3. Five sncRNAs were validated; miR-223 was significantly reduced in early osteoarthritis and miR-23b, let-7a-2, snord96A and snord13 were significantly upregulated. Significant cellular actions deduced by the differentially expressed microRNAs included apoptosis (P < 0.0003), necrosis (P < 0.0009), autophagy (P < 0.0007) and inflammation (P < 0.00001). A conservatively filtered list of 57 messenger RNA targets was obtained; the top biological processes associated were regulation of cell population proliferation (P < 0.000001), cellular response to chemical stimulus (P < 0.000001) and cell surface receptor signalling pathway (P < 0.000001). Conclusions Synovial fluid sncRNAs may be used as molecular biomarkers for early disease in equine osteoarthritic joints. The biological processes they regulate may play an important role in understanding early osteoarthritis pathogenesis. Characterising these dynamic molecular changes could provide novel insights on the process and mechanism of early osteoarthritis development and is critical for the development of new therapeutic approaches.


Background
Osteoarthritis (OA) remains one of the greatest causes of morbidity and mortality for horses in the UK [1,2]. Additionally, it is the most common disease affecting the joints in humans, and a significant cause of pain and disability worldwide [3]. This degenerative, age-related joint disease is characterised by a progressive degradation of articular cartilage and concomitant structural and functional change of all joint constituents, including the synovial membrane, the subchondral bone and periarticular tissues [4]. Of multifactorial origin, OA is a product of genetic, mechanical and environmental factors such as age, trauma and occupation [4,5]. Despite its high prevalence and significant welfare and economic impact, its pathophysiology remains poorly understood and currently available diagnostic tools can only identify the disease when cartilage has already exceeded its capacity for intrinsic repair, and changes can no longer be reversed [6,7]. As a result, the development of effective treatments is also compromised, and currently recommended therapies are mainly symptomatic.
In the search for molecular biomarkers that could reveal pre-clinical phases of the disease, scientists have focused much of their attention on microRNAs (miRNAs), the best characterised family of small non-coding RNAs. Evolutionarily conserved, these 17-22 nucleotide long molecules regulate gene expression at post-transcriptional level generally by repressing translation or increasing degradation of messenger RNAs (mRNAs). They are involved in different cellular pathways and intercellular communication thus influencing tissue homeostasis [8]. As such, miRNA profiles can be altered as a result of cellular damage and/or tissue injury and altered expression of certain miRNAs is implicated in several diseases, including OA [9][10][11]. miRNAs can promote cell differentiation by modulating expression of catabolic genes; for instance, miR-139 which is increased in OA cartilage inhibits cell proliferation by suppressing expression of insulin-like growth factor 1 receptor (IGF1R) and eukaryotic translation initiation factor 4 gamma 2 (EIF4G2) [12]. Furthermore, miRNA expression can also regulate proinflammatory cytokines, again contributing to altered expression of OA-inducing genes; for example, in lipopolysaccharide (LPS)-treated mouse chondrocytes, inhibition of miR-203 increases apoptosis and further stimulates the production of inflammatory cytokines [13]. Additionally, miRNAs can play cartilage-protective roles; miR-193b-3p inhibits extracellular matrix (ECM) degradation through inhibition of inducible nitric oxide (NO) synthesis [14]; dysregulation of miR-193b-3p can therefore promote cartilage degeneration. Comprehensive reviews on miRNAs involved in osteoblastogenesis and osteoclastogenesis, chondrogenesis and cartilage degradation, synovial inflammation and neurogenesis can be found elsewhere [15][16][17]. miRNAs can be found intracellularly or extracellularly, circulating in virtually any biological fluid in a remarkably stable manner [18][19][20]. Because biological fluids are generally obtainable through minimally invasive techniques, circulating miRNAs are attractive candidates for disease diagnosis, monitoring and prognostication [21,22]. Interest in other classes of small non-coding RNAs such as small nucleolar RNAs (snoRNAs) has recently emerged. Mostly known for their housekeeping functions, snoRNAs have canonical roles in the chemical modification of RNA substrates such as ribosomal RNAs, but can also exhibit miRNA-like activity [23]. Aberrant expression of snoRNAs has also been associated with the development of different diseases and a recent study found alterations in the snoRNA profile of OA joints in mice when compared to healthy controls, highlighting the potential of snoRNAs to be used as novel markers for this disease [24].
Equine miRNAs have been identified in numerous healthy tissues [25,26] and their potential role in different diseases such as osteochondrosis, rhabdomyolysis and insulin resistance has also been investigated [27][28][29]. However, information on miRNA influence on the pathogenesis of equine OA is still lacking. Synovial fluid represents a reliable source of chemical information that can accurately reflect pathological conditions affecting the joint due to its functional proximity within joint tissues [30]. In 2010, Murata et al. investigated the presence and stability of miR-NAs in human synovial fluid for the first time, and found five differentially expressed miRNAs in OA patients compared to healthy controls, supporting the potential use of synovial fluid miRNAs as diagnostic biomarkers [11]. More recently, a screening of 752 miRNAs in synovial fluid from human patients with early-and late-stage OA demonstrated seven upregulated miRNAs in late-stage OA, irrespective of age, gender and body mass index [31]. Intra-articular treatment with hyaluronic acid was shown to modify miRNA expression in OA patients [32]. Although miRNA expression has not yet been investigated in equine OA, a preliminary study has recently described a reproducible method for miRNA isolation from equine synovial fluid and blood plasma [33].
With growing evidence of alterations in small noncoding RNA patterns in the synovial fluid of OA joints, we theorised that early stages of OA would affect these molecules and potentially provide early biomarkers for OA in equine patients. Examining expression of small non-coding RNAs in synovial fluid in early OA may also provide further insights on the pathological changes that occur. Therefore, we investigated the profile of small non-coding RNAs of early equine OA synovial fluid using next generation sequencing.

Macroscopic and histological assessment
The donors used for small RNA sequencing were selected from an elderly population of horses to account for any age-related changes. The ages of the control (mean ± standard deviation; 22 ± 2 years) and early OA (27 ± 7.5 years) groups were not significantly different. Horses included in the control group presented minor macroscopic or histological changes, which are to be expected in healthy older animals. Horses included in the early OA group presented intermediate OA scores and were not obviously lame prior slaughter, suggestive of primary OA.
For samples used for small RNA sequencing there was a significant increase in the macroscopic score between control (1.0 ± 0.5) and early OA (5.4 ± 1.9, P = 0.04). Likewise, there was a significant increase in the histological score between control (2.1 ± 0.7) and early OA (6.1 ± 1.5, P = 0.01) (Additional File 1).
For the independent cohort the ages of the control and early OA groups were not significantly different (Additional File 1). There was a significant increase in the macroscopic score between control (1.75 ± 1.5) and early OA (3.6 ± 0.9, P = 0.04) samples. Similarly, there was a significantly increase in the histological score between control (1.5 ± 1.3) and early OA (5.8 ± 2.5, P = 0.02) (Additional File 1).

Analysis of small RNA sequencing data
Summaries of raw, trimmed reads and mapped reads to the Equus caballus database are in Additional File 2. There were 323 small non-coding RNAs identified. The categories of RNA identified are in Fig. 1a and included small non-coding RNAs; miRNAs, snoRNAs and small nuclear RNAs (snRNAs).

Confirmation of differential gene expression using qRT-PCR
Seven small non-coding RNAs (miR-143, miR-223, miR-99a, miR-23b, let-7a-2, snord96A, snord13) were selected for further validation based on our current work, level of differential expression (P < 0.05 and logFC> 1.2) and following a literature review of differentially expressed genes. An independent cohort of synovial fluid samples was used, comprising of control (n = 6, histological score 1.5 ± 1.3) and early OA (n = 6, histological score 5.8 ± 2.5) synovial fluid samples. In agreement with the sequencing data miR-223 was significantly reduced in early OA and miR-23b, let-7a-2, snord96A and snord13 were significantly increased in early OA (Fig. 3). For two miRNAs, miR-143 and miR-99a-2 quantitative reverse transcription-polymerase chain reaction (qRT-PCR) findings did not validate sequencing findings; despite being decreased in the OA group in our sequencing data, qRT-PCR showed increased expression of both miRNAs in the independent OA group compared to controls, although this was not statistically significant (Fig. 3).

Identification of potential target mRNA genes of the differentially expressed miRNAs
With the goal of exploring potential biological associations with the differentially expressed miRNAs in early OA synovial fluid we undertook an Ingenuity Pathway Analysis (IPA) 'Core Analysis' on these. Interesting features were determined from the gene networks inferred. Significant cellular actions deduced by the differentially expressed miRNAs included apoptosis (P < 0.0003), necrosis (P < 0.0009), autophagy (P < 0.0007) and inflammation (P < 0.00001) (Fig. 4).
Next, we undertook analysis to determine the mRNA targets of the differentially expressed miRNAs. Eight miRNAs were differentially expressed in early OA compared to non-OA controls. Once a conservative filter was applied (only miRNAs with experimentally confirmed or highly conserved predicted targets), miR-let7a-2 and miR-378 were excluded. Six miRNAs remained which collectively putatively target 993 mRNAs. We then additionally added the filters chondrocytes, fibroblast and osteoblasts, removed duplicates and obtained a list of 57 mRNA targets (Additional File 3).
The presumed target mRNAs were input into the gene ontology (GO) tool PANTHER and the biological processes were summarised in REViGO and visualised using Cytoscape (Fig. 5). The top biological processes were regulation of cell population proliferation (false discovery rate (FDR)adjusted P = 6.24E − 13 ), cellular response to chemical stimulus (FDR = 4.54E − 12 ) and cell surface receptor signalling pathway (FDR = 6.39E − 12 ) (Additional File 4).

Discussion
The inability to detect pre-clinical changes in OA has been one of the main impediments to the development of effective therapies against this disease [34]. From a biomarker perspective, profiling synovial fluid circulating locally within the affected joint cavity at an early stage may provide new insights into pathological changes occurring during OA initiation and progression, and ultimately allow for the implementation of new therapeutic approaches. Our study is, to the best of our knowledge, the first to characterise the small non-coding RNA profile of synovial fluid in early OA in horses, providing evidence of a pattern of differential expressed synovial fluid miRNAs and other small non-coding RNAs in early OA synovial fluid when compared to our control samples.
Osteoarthritis is a highly heterogeneous disease and can be broadly divided into primary, naturally occurring OA, which is chronic and associated with age; and posttraumatic OA, usually related to athletic use [35]. Whilst post-traumatic OA is highly prevalent and has a significant economic impact particularly for younger and athletic horses [36], animals over 15 years comprise up to one third of the equine population and represent a growing proportion of referral hospital admissions [37]. Musculoskeletal conditions are a major reason for euthanasia in older horses, suggesting that the social economic burden of age-related osteoarthritis is rising [2]. While trauma might be one of the causative factors of OA in older horses, it is difficult to ascertain whether the body's response to external sources of stress is being affected by ageing or whether age-related changes are actually a predisposing factor for traumatic injuries by increasing their likelihood, for example through altered biomechanics. For this reason, in our study we sought to exclude the cofounding effect of a traumatic injury by selecting older donors originating from the abattoir; although this does not guarantee absence of a traumatic injury, horses can only enter the food chain if they are not obviously lame prior to slaughter, making it more likely that these OA cases arise from age-related molecular mechanisms. Selected donors from both control and early OA groups were age-matched to account for any age-related changes. Donors from the control group presented minor changes in their joints which are to be expected in older horses, hence classifying this group as "control" as opposed to "healthy".
Due to the considerable interest in miRNA-mediated gene regulation in recent years, the list of miRNAs possibly implicated in OA and other joint related pathologies has grown [16]. miRNAs that are differentially expressed in joint tissues of patients with OA are likely to contribute to OA pathophysiology and may be utilised as diagnostic factors [38]. One example is miR-140, which is significantly downregulated in human OA cartilage [10] and is thought to attenuate OA progression by modulating ECM homeostasis [39]; also, dysregulation of miR-140-3p and-5p in synovial fluid has been correlated with OA severity [40].
Among the differentially expressed miRNAs found in our study, miR-23b was significantly increased in the early OA cohort. miR-23b is thought to be involved in OA progression by targeting cartilage-associated protein (CRTAP) and thus influencing cartilage homeostasis [41]. This miRNA has also been shown to positively regulate the chondrogenic differentiation of mesenchymal stem cells by regulating the expression of sexdetermining region Y-Box 9 (SOX9) and protein kinase A (PKA) [42,43].
Likewise, we found let-7a-2 to be upregulated in early OA. In an experiment comparing miRNA expression in synovial fluid from human OA patients undergoing hyaluronic acid treatment, let-7a was significantly upregulated in synovial fluid of OA samples compared to healthy controls; levels of let-7a in affected patients returned to normal after hyaluronan injection [32]. Let-7a is thought to regulate IL-6 receptor (IL6R), and its inhibition can enhance cell proliferation, reduce apoptosis and inhibit inflammatory response in ATDC5 cells in a LPS-induced in vitro model of OA [44]. Members of the let-7 family have often been described in studies involving OA; a large population-based study identified serum let-7e as a promising candidate to predict OA risk, independent of age, sex and body mass index [45]. A recent investigation supported this claim, providing further evidence of decreased expression of let-7e in serum of patients affected with knee OA [46]. The exact roles of miRNAs of the let-7 family remain unclear, but the evidence for their use as biomarkers for OA is growing.
In the above mentioned paper by Xu et al. (2015) [32], miR-223 was also significantly upregulated in synovial fluid of OA patients prior to intra-articular injection of hyaluronan. miR-223 participates in cartilage homeostasis and structure by targeting growth differentiation factor 5 (GDF5) [41]. Early-stage OA patients showed upregulation of miR-223 in peripheral blood mononuclear cells, with its expression decreasing as OA progressed [47]. In our study, we found miR-223 to be downregulated in the synovial fluid of the early OA cohort, which supports the involvement of this miRNA in the early osteoarthritic process. miRNA regulation is complex and differences between our results and previously published literature in human patients may be due to different stages in osteoarthritic process; species variation may also partially justify these disparities. Additionally, an increasing body of evidence demonstrates that long non-coding RNAs (lncRNAs) can act as sponges for microRNAs [48]; a previous study found that the expression of miR-223 was restrained by lncRNA activated by transforming growth factor beta (lncRNA-ATB) [49] which might contribute to variations in miRNA expression.
We have previously shown the involvement of snoR-NAs in cartilage ageing and OA and their potential use as biomarkers for OA [24]. In this study we identified for the first time snord13 and snord96a as highly Fig. 2 Volcano plot of small non-coding RNAs identified represents logFC and -log10 P value. Pink dots represent differentially expressed small non-coding RNAs expressed small non-coding RNAs in early OA. Our previous work in human OA cartilage identified a dysregulation in SNORD96A expression in ageing and OA. In addition, we demonstrated changes in chondrogenic, hypertrophic, ribosomal RNA (rRNA) and OA related gene expression following overexpression and knockdown of SNORD96A in human chondrocytes. Interestingly we also identified an increase in SNORD96A in chondrocytes treated with OA synovial fluid [50]. In another microarray study of young compared to old OA cartilage we identified SNORD13 was increased in OA cartilage [51]. Together these findings indicate that changes in synovial fluid snoRNAs could in part be due to a dysregulation in their expression in cartilage in OA. snoRNAs are emerging with unappreciated functional roles in cell physiology [52] and our results support our earlier work for the potential use of snoRNAs as novel biomarkers in OA [24].
Predicted targets of the miRNAs of interest appear to be involved in processes of inflammation and cellular destruction including necrosis, apoptosis and autophagy, which have been previously shown to contribute to the pathogenesis of OA in human patients through proinflammatory cytokines production [53], synovial inflammation [54] and chondrocyte apoptosis [55]; subchondral bone changes [56] and chondrocyte apoptosis [57] have also been implicated in the pathogenesis of OA in horses.
For example, a disintegrin and metalloproteinase with thrombospondin motif 1 (ADAMTS1) is known to cleave aggrecan, a critical component for cartilage structure [58]; bone morphogenetic protein (BMP) receptor type 1B (BMPRB) is a receptor for BMP, and BMP signalling is essential for chondrocyte proliferation, survival and differentiation [59]; and IL-6R interacts with IL-6, one of the pro-inflammatory cytokines increased in osteoarthritis [60]. Experimental validation of these and other predicted target genes can clarify biological mechanisms behind small non-coding RNAs of interest and elucidate their role in the pathogenesis of OA; this is critical for the success of future interventions, as these molecules can be targeted in a specific manner [61,62].
Profiling circulating, cell-free small non-coding RNAs is generally a challenging task due to the limited amount of RNA present in biofluids, as well as presence of inhibitory compounds which potentially hinder downstream enzymatic processes. However, liquid biopsies for the investigation of non-coding RNA profiles have gained prominence due to their ease of collection and potential use as diagnostic tools. Future studies in this field would benefit from analysing larger cohorts of patients; our study was limited by the availability of joints Fig. 3 Validation of small non-coding RNAs differentially expressed following small RNA sequencing in an independent cohort using qRT-PCR. RNA extracted from the synovial fluid of six healthy control donors and six early OA donors. Histograms of the relative expression calculated using 2^-DCT method using the geometric mean of miR-100 and miR191 as an endogenous control. All qRT-PCR reactions were performed in triplicate. Statistical significance was tested in Graphpad Prism using a Mann Whitney test. Bars represent means with standard error of the mean. P < 0.05; * with early OA, resulting in a small sample size. Notwithstanding, a previous study on the subject of RNA-Seq analysis performance [63] has shown the number of genes called significant increases as the sample number increases; this suggests that for pipelines such as the one used in this study, having a slightly underpowered approach means we are more likely to underestimate rather than overestimate the number of differentially expressed miRNAs. That fact that we were able to validate our findings through qRT-PCR in an independent cohort solidifies our findings, despite the small sample size. Further work is unquestionably needed, yet this experiment enabled us to identify small non-coding RNA changes in the initial and an additional cohort and revealed, for the first time, the potential use of small noncoding RNAs as biomarkers for early OA. These results support the use of synovial fluid small non-coding RNAs as molecular biomarkers for early disease in OA joints. Our future research is currently ascertaining the applicability of these findings in a clinical setting.

Conclusions
This study demonstrates that equine synovial fluid displays a pattern of small non-coding RNA differential expression in early OA when compared to controls, as defined by gross and histological scoring and many of these small non-coding RNAs have previously been demonstrated to have a role in OA. The affected biological cellular processes in response to changing miR-NAs and their target genes might play an important role in early OA pathogenesis. This opens the possibility of a relatively non-invasive method for early detection of OA. Furthermore, characterisation of these dynamic molecular changes could provide novel insights on the process and mechanism of early OA development.

Methods
All reagents were from ThermoFisher Scientific, unless stated.

Sample collection and preparation
Samples were collected from the metacarpophalangeal joints of horses from an abattoir as a by-product of the agricultural industry. Specifically, the Animal (Scientific procedures) Act 1986, Schedule 2, does not define collection from these sources as scientific procedures. Ethical approval was therefore not required.
The joints were aseptically dissected to allow visual inspection of the metacarpus, the proximal phalange and the sesamoids. All joints were photographed and macroscopic changes were scored based on a scoring system as previously described [64]. Synovial fluid was aseptically Fig. 5 The position of differentially expressed miRNAs in the chondrocyte/fibroblast/osteoblast expression network. PANTHER was used to identify gene ontology (GO) biological processes associated with predicted mRNA targets and perform overrepresentation analysis to highlight the GO terms most significantly affected by dysregulated miRNA-mRNA interactions in early OA synovial fluid. GO terms (FDR > 0.05) were summarised and visualised using REViGO and Cytoscape. Allowed similarity setting in REViGO was tiny (0.4). The line width specified the amount of similarity collected directly from the open joint with a 5 ml sterile syringe, immediately placed on sterile microcentrifuge tubes on ice and centrifuged for 10 min at 3000 g and 4°C to remove cells and debris. The supernatant was collected and stored at − 80°C. A cartilage and subchondral bone fragment was collected from the palmar aspect of one of the metacarpal condyles, fixed on paraformaldehyde and sent for histology; histological scoring was performed using the previously described scoring system [65].
Donors were assigned to groups based on the macroscopic and histologic scoring. The control (non-OA) group was comprised of 3 donors with age mean ± standard deviation 22 ± 2 years; while the early OA group was comprised of 3 donors with 22 ± 7.5 years.

RNA isolation, cDNA library preparation and small RNA sequencing
Synovial fluid was treated to reduce viscosity with 1 μg/ ml of hyaluronidase at 37°C for 1 h, centrifuged at 1000 g for 5 min, and supernatant used for total RNA extraction using miRNeasy serum kits (Qiagen, Crawley, UK). The integrity of the RNA was assessed on the Agilent 2100 Bioanalyzer system using an RNA Pico chip. 100 ng samples were submitted for library preparation using NEBNext® Small RNA Library Prep Set for Illumina (New England Biosciences (NEB), Ipswich, USA) but with the addition of a Cap-Clip™ Acid Pyrophosphatase (Cell script, Madison, USA) step to remove any 5′ cap structures [24] and size selected using a range 120-300 bp. This enabled both miRNAs and snoRNAs to be identified in a non-biased approach. The pooled libraries were sequenced on an Illumina HiSeq4000 platform with version 1 chemistry to generate 2 × 150 bp paired-end reads. Data has been submitted to National Centre for Biotechnology Information; accession E-MTAB-8409.

Small RNA sequencing data analysis
Sequence data were processed through a number of steps to obtain non-coding RNA expression values including; basecalling and de-multiplexing of indexed reads using CASAVA version 1.8.2; adapter and quality trimming using Cutadapt version 1.2.1 [66] and Sickle version 1.200 to obtain fastq files of trimmed reads; aligning reads to horse genome reference sequences (release 90) from Ensembl using Tophat version 2.0.10 [67] with option "-g 1"; counting aligned reads using HTSeqcount [68] against the features defined in horse genome GTF file (release 90).
Differential expression analysis was performed in R using package DESeq2 [69]. The processes and technical details of the analysis include; assessing data variation and detecting outlier samples through comparing variations of within and between sample groups using principle component analysis (PCA) and correlation analysis; handling library size variation using DESeq2 default method; formulating data variation using negative binomial distributions; modelling data using a generalised linear model; computing logFC values for control versus early OA based on model fitting results through contrast fitting approach, evaluating the significance of estimated logFC values by Wald test; adjusting the effects of multiple tests using FDR approach [70] to obtain FDR adjusted P-values.
The Ensembl horse genome GTF file release 90 does not have mature miRNA features. We linked the defined miRNA primary transcripts to miRBase horse miRNA GFF3 file by feature's genome coordinates so as to obtain the corresponding mature miRNA.

qRT-PCR validation
Validation of the selected small RNA sequencing results in an independent cohort of equine metacarpophalangeal synovial fluid was undertaken using qRT-PCR. Six control (non-OA), mean ± standard deviation (20.2 ± 2.4 years) and six early OA (20.8 ± 4.1) with macroscopically and histologically graded sample scores similar to those used for sequencing were used. Total RNA was extracted as above. Small non-coding RNAs were chosen based on our current work, level of differential expression (P < 0.05 and logFC> 1.2) and following a literature review of differentially expressed genes. These were miR-143, miR-223, miR-99a, miR-23b, let-7a-2, snord96A and snord13. Primer sequences/assays used can be found in Additional File 5. PolyA cDNA was synthesized using 200 ng RNA and the miScript II RT Kit. A mastermix was prepared using the miScript SYBR Green PCR Kit (Qiagen, Crawley, UK) and the appropriate bespoke designed miScript Primer Assays (Qiagen, Crawley, UK). Real-time PCR was undertaken using a LightCycler® 96 system (Roche). Steady-state transcript abundance of potential endogenous control genes was measured in the small RNA sequencing data. Assays for four genes -miR-181a, miR-100, miR-191a and U6 were selected as potential reference genes because their expression was unaltered in this study. Stability of this panel of genes was assessed by applying a gene stability tool RefFinder [71]. The geometric mean of miR-100 and miR-191a was selected as the stable endogenous control. miR-100 has been previously used as a normaliser in a similar study as it was identified by NormFinder as the most stable [31]. Relative expression levels were normalised to the geometric mean of miR-100 and miR-191 and calculated using the 2^-DCT method [72].

miRNA target prediction and pathway analysis
Potential biological associations of the differentially expressed miRNAs in early OA synovial fluid were identified using IPA (IPA, Qiagen Redwood City, CA, USA) 'Core Analysis'. Canonical pathways, novel networks, and common upstream regulators were then queried. Additionally in order to identify putative miRNA targets, bioinformatic analysis was performed by uploading differentially expressed miRNA data into the MicroRNA Target Filter module within IPA software. This identifies experimentally validated miRNA-mRNA interactions from TarBase, miRecords, and the peer-reviewed biomedical literature, as well as predicted miRNA-mRNA interactions from TargetScan. We used a conservative filter at this point, using only experimentally validated and highly conserved predicted mRNA targets for each miRNA. Targets were then also filtered on the cells chondrocyte, osteoblasts and fibroblasts (the latter two settings were the nearest to bone and synovial cells available for selection), to represent joint cells in contact with synovial fluid.
PANTHER (GO Ontology database 2020-02-21) [73] was used for overrepresentation analysis of the mRNA targets using Fisher's Exact test with FDR correction. This tests whether the input mRNAs associate significantly with specific pathways and generates a list of biological process GO terms. Terms with FDR adjusted P < 0.05 were summarised using REViGO [74] with allowed similarity of 0.4 and visualised using Cytoscape [75].

Statistical analysis
The heatmap and volcano plots were made using Meta-boAnalyst 3.5 (http://www.metaboanalyst.ca) which uses the R package of statistical computing software [76]. For statistical evaluation of gene expression data, following normality testing, Mann-Whitney tests were performed using GraphPad Prism version 8.0 for Windows (Graph-Pad Software, La Jolla California USA, www.graphpad. com); P values are indicated.
Additional file 1. Histograms of age, gross score and Modified Mankin's Score for dependent and independent equine donor cohorts (.tiff). Expressions are means and error bars ± standard error means. Statistical analysis undertaken in GraphPad Prism 8.0 using a Mann Whitney Test. P values *; P < 0.05.
Additional file 2.Summary of raw, trimmed reads and mapped reads (.xlsx). Summary of raw, trimmed reads and mapped reads to Equus caballus database, from analysis of small RNA sequencing data.
Additional file 3. mRNA targets predicted by IPA (.xlsx). List of mRNA targets predicted by bioinformatic analysis with IPA software, using a conservative filter of only experimentally validated and highly conserved predicted mRNA targets for each miRNA. Targets were then also filtered on the cells chondrocyte, osteoblasts and fibroblasts.
Additional file 4. PANTHER GO terms FDR-adjusted P < 0.05 (.xlsx). List of GO terms with FDR-adjusted P < 0.05, obtained with PANTHER overrepresentation analysis of the mRNA targets using Fisher's Exact test.
Additional file 5. Primer sequences/assays used for detection of small non-coding RNAs through qRT-PCR analysis (.xlsx). For miRNAs and snoR-NAs with sequences homologous to human, Qiagen primer assays were used. Remaining miRNA primers were customised using Eurogentec primer design.