Small Non-Coding RNAome of Ageing Chondrocytes

Ageing is a leading risk factor predisposing cartilage to osteoarthritis. However, little research has been conducted on the effect of ageing on the expression of small non-coding RNAs (sncRNAs). RNA from young and old chondrocytes from macroscopically normal equine metacarpophalangeal joints was extracted and subjected to small RNA sequencing (RNA-seq). Differential expression analysis was performed in R using package DESeq2. For transfer RNA (tRNA) fragment analysis, tRNA reads were aligned to horse tRNA sequences using Bowtie2 version 2.2.5. Selected microRNA (miRNAs or miRs) and small nucleolar RNA (snoRNA) findings were validated using real-time quantitative Polymerase Chain Reaction (qRT-PCR) in an extended cohort of equine chondrocytes. tRNA fragments were further investigated in low- and high-grade OA human cartilage tissue. In total, 83 sncRNAs were differentially expressed between young and old equine chondrocytes, including miRNAs, snoRNAs, small nuclear RNAs (snRNAs), and tRNAs. qRT-PCR analysis confirmed findings. tRNA fragment analysis revealed that tRNA halves (tiRNAs), tiRNA-5035-GluCTC and tiRNA-5031-GluCTC-1 were reduced in both high grade OA human cartilage and old equine chondrocytes. For the first time, we have measured the effect of ageing on the expression of sncRNAs in equine chondrocytes. Changes were detected in a number of different sncRNA species. This study supports a role for sncRNAs in ageing cartilage and their potential involvement in age-related cartilage diseases.


Introduction
Articular cartilage is a specialised connective tissue of diarthrodial joints. Its smooth lubricated surface assists joint movement and its mechanical properties facilitate load bearing in the joint. The tissue harbours one cell type; the chondrocyte, and is devoid of blood vessels and nerves, receiving nutrients from synovial fluid and subchondral bone [1]. Articular cartilage is characterised by an extracellular matrix (ECM) consisting of mainly collagen type 2 and proteoglycans, which give the tissue many of its properties. After reaching maturity, cartilage displays a limited repairing capacity as indicated by low chondrocyte proliferation and low collagen turnover [2].
There are a number of factors affecting the homeostatic properties of cartilage such as genetics and obesity [3]. However, ageing is a leading risk factor that predisposes cartilage to pathological changes and disease, such as osteoarthritis (OA), the most common joint disease [3]. These age-related fragments are cleaved at the T loop by Dicer [23] and ANG [24] and tRF-5 fragments are derived from the cleavage of the D-loop by Dicer [23]. The final category of tRFs are i-tRFs, which are internal to the respective tRNA and can straddle the anticodon loop [25]. Limited knowledge of the expression and role of tRNA and tRFs is available in health and disease [26] with even less in musculoskeletal biology [27,28].
In this study, we investigated the expression changes of sncRNAs in chondrocytes isolated from healthy metacarpophalangeal joints of young and old horses. Findings for tRNAs are compared to a human OA data set. Age-related changes may predispose cartilage to disease by altering the complex sncRNA expression profile. This provides a sncRNA-wide insight into age-related targets for future therapeutic approaches.

Materials and Methods
All reagents were from Thermo-Fisher-Scientific, Loughborough, UK unless stated.

Sample Collection and Preparation
Samples were collected 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. Full thickness equine cartilage was removed from the entire surface of macroscopically normal metacarpophalangeal joints of young n = 5 (age mean ± standard deviation; 4 ± 1 years) and old n = 5 (17.4 ± 1.9 years) non-Thoroughbred horses. Scoring of the metacarpophalangeal joint was undertaken using a macroscopic grading system, as previously described [29] and samples with no macroscopic perturbations were selected (combined score of zero). Freshly isolated chondrocytes were removed from all 10 samples of harvested cartilage, as previously described [30], plated to confluence and RNA extracted from two million cells per donor.
In addition to the above samples, RNA from chondrocytes from young n = 2 (age mean ± standard deviation; 0.75 ± 0.3 years) and old n = 6 (age mean ± standard deviation; 19.3 ± 3.6 years) non-Thoroughbred horses were used for validation.

RNA Isolation, cDNA Library Preparation, and Small RNA Sequencing (RNA-seq)
Total RNA including small RNAs was extracted, as previously described [31], and purified using the miRNAeasy kit (Qiagen, Crawley, UK) according to manufacturer's instructions and including an on-column DNAse step to remove residual genomic DNA. The integrity of the RNA was assessed on the Agilent 2100 Bioanalyser system using an RNA Pico chip (Agilent, Stockport, UK). The NEBNext®Small RNA Library Prep Set for Illumina®was used for library preparation (New England Biolabs, Ipswich, MA, USA) but with the addition of a Cap-Clip™ Acid Pyrophosphatase (Cell script, Madison, WI, USA) step to remove potential 5 caps found on some snoRNAs. Samples were amplified for 15 cycles and size was selected. The libraries were sequenced on an Illumina MiSEq platform (Illumina, San Diego, CA, USA) with version 2 chemistry using sequencing by synthesis technology to generate 2 × 150 bp paired-end reads with >12 million clusters per run.

Data Processing
Sequence data were processed through a number of steps to obtain sncRNA expression values including basecalling and de-multiplexing of indexed reads using CASAVA version 1.8.2 [32]; adapter and quality trimming using Cutadapt version 1.2.1 [33] and Sickle version 1.200 to obtain fastq files of trimmed reads; aligning reads to horse genome reference sequences (GCF_002863925.1) using Tophat version 2.0.10 [34] with option "-g 1"; counting aligned reads using HTSeq-count against the annotated features which are combined annotation information from the sources: NCBI Equus caballus 3.0 genome annotation, miRBase horse micro RNA annotation, Rfam snoRNA annotation.
Differential expression analysis was performed in R environment using package DESeq2 [35]. The processes and technical details of the analysis included: 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 the DESeq2 default method; formulating data variation using negative binomial distributions; modelling data using a generalised linear model; computing log2 Fold Change (logFC) values for required contrasts 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 False Discovery Rate (FDR) approach to obtain FDR [36] adjusted P-values.

Pathway Analysis
In order to identify miRNA targets, bioinformatic analysis was performed by uploading differentially expressed miRNA data into the MicroRNA Target Filter module within Ingenuity Pathway Analysis software (IPA) (Qiagen Redwood City, CA, USA) along with previously identified differentially expressed mRNAs from our ageing equine cartilage study following RNA-seq [31]. In IPA we selected miRNA-target genes based on the direction of differential expression (for example, if a miRNA was reduced in expression it was only matched to mRNAs that demonstrated increased expression). We then identified the networks, functions, and canonical pathways of these miRNA-target genes.

qRT-PCR Validation
Validation of the small RNA-seq results was undertaken using real-time quantitative PCR (qRT-PCR) analysis in the samples used for sequencing as well as additional samples. SncRNAs were chosen based on level of differential expression. Total RNA was extracted and quantified as above. cDNA was synthesized using 200 ng RNA and the miScript II RT Kit according to the manufacturer's protocol (Qiagen, Crawley, UK). qRT-PCR mastermix was prepared using the miScript SYBR Green PCR Kit (Qiagen, Crawley, UK) and the appropriate miScript Primer Assay (Qiagen, Crawley, UK) (File S1) using 1 ng/µL cDNA according to manufacturer's guidelines. qRT-PCR was undertaken using a LightCycler®96 system (Roche, Welwyn Garden City, UK). Relative expression levels were normalised to U6 (as this was stable in the small RNA-seq data set) and calculated using the 2-∆Ct method [37].

tRNA Fragment Analysis
Following the alignment of trimmed reads to NCBI horse genome reference sequences (version 3.0) using Tophat version 2.1.0 [38], the candidate tRNA reads were extracted from the BAM files according to whether they overlapped the ranges covered by tRNA features. The read pairs were stitched into RNA fragments using PEAR (version 0.9.10) [39]. The output reads were aligned to horse tRNA sequences (defined in NCBI GCF_002863925.1_EquCab3.0_genomic.gff) using Bowtie2 version 2.2.5. Only the perfectly mapped fragments were extracted and taken as tRNA fragments for further explorations. Finally, statistical analyses were mainly focused on the fragment length and the mapping start location, which generated the length distribution and the mapping start position distribution of observed tRNA fragments, as well as the summary table for observed tRNA fragments and their target tRNAs.

Novel snoRNA Analysis
Putative snoRNAs were detected from the raw pair-end reads using ShortStack 3.8.4 [40] with the setting "-mincov 5", which specifies that the clusters of small RNAs must have at least five alignments. The ShortStack results were subsequently fed into SnoReport 2.0 [41], which uses RNA secondary structure prediction combined with machine learning as the basis to identify the two main classes of snoRNAs; the box H/ACA and the box C/D. Putative snoRNAs were annotated from our experimental small RNA-seq data using ShortStack and SnoReport.

Statistical Analysis
3D PCA score plots and heat maps were carried out using MetaboAnalyst 3.5 [42] (http://www. metaboanalyst.ca) which uses the R package of statistical computing software [43].
For statistical evaluation of qRT-PCR results, following normality testing, Mann-Whitney tests were performed using GraphPad Prism version 7.03 for Windows, (GraphPad Software, La Jolla, CA, USA); p values are indicated.

Human Sample Collection and Preparation
De-identified human OA cartilage approved by Northeast Ohio Medical University (NEOMED) Institutional Review Board and Summa Health Systems, Barberton, Ohio as 'non-human subject study under 45 CFR' was used. Total RNA was extracted from smooth, macroscopically intact human OA cartilage with a Mankin score of 2 or less (n = 1, female, 60 years old) and damaged OA cartilage with a Mankin score of 4 or higher (n = 1, female, 80 years old) using MiRNeasy Kit (Qiagen, Germantown, MD, USA). RNA was quantified using the Nanodrop 1000 Spectrophotometer (Thermo Fisher, Waltham, MA, USA). TapeStation 4200 (Agilent Technologies, Santa Clara, CA, USA) was used to determine RNA integrity using High Sensitivity RNA Screentape analysis kit (Agilent, Santa Clara, CA, USA).

Removal of tRF Modifications from Human Cartilage RNA
The rtStar™ tRF&tiRNA Pretreatment Kit (ArrayStar, Rockville, MD, USA) was used according to the manufacture's description. For cDNA qRT-PCR library construction, tRF modifications were removed from RNA. The kit removes 3 -aminoacyl and 3 -cP for 3 adaptor ligation, phosphorylates 5 -OH for 5 -adaptor ligation, and demethylates m1A, m1G, and m3C for efficient cDNA reverse transcription.

tRF 3 and 5 Adaptor Ligation for Human Cartilage cDNA Synthesis
The rtStar™ First-Strand cDNA Synthesis Kit (ArrayStar, Rockville, MD USA) was used according to the manufacturer's description. The kit sequentially ligates 3 -Adaptor with its 5 -end to the 3 -end of the RNAs, and 5 -Adaptor with its 3 -end to the 5 -end of the RNAs. The non-ligation ends of 3 and 5 Adaptors were blocked by modifications. A universal priming site for reverse transcription was contained within the 3 adaptor. Spike-in RNA was used for monitoring the cDNA synthesis efficiency and as a quantitative reference.

tRF and tRNA Human qPCR Arrays
The nrStar™ Human tRF and tiRNA PCR Array (ArrayStar, Rockville, MD, USA) was used according to the manufacturers description to profile 185 tRF and tiRNA fragments, of which 101 are derived from tRF and tiRNA database [44,45] and the other 84 from recently published papers [46][47][48]. RNA Spike in control and a positive PCR control were used to evaluate PCR efficiency and a genomic DNA control was used to monitor genomic DNA contamination. To profile parent tRNAs, the nrStat Human tRNA PCR Array (ArrayStar, Rockville, MD, USA) was used according to the manufacturer's protocol. This array consisted of 163 PCR-distinguishable nuclear tRNA isodecoders and 22 PCR distinguishable mitochondrial tRNA species covering all anti-codons compiled in GtTNAdb [49,50] and tRNAdb [51] databases. Genomic DNA and positive PCR controls were included to monitor the quality of RNA sample.
qRT-PCR reactions were conducted using Power SYBR Green master mix (Life technologies, Carlsband, CA, USA) on a Step One Plus (Applied Biosystems, Waltham, MA, USA) machine. U6, SNORD43 and SNORD45 were used as endogenous controls for normalisation of tRF, tiRNA and tRNA detection. Target tRNA, tiRNA and tRF levels were determined as fold change differences utilising the ∆∆Ct method [37].

Preliminary Analysis of Small RNA-seq Data
To identify differential expression of sncRNAs in ageing, Illumina MiSeq was utilised. Summaries of raw, trimmed reads, and mapped reads to Equus caballus database can be found in File S2. Reads mapping percentages were between 92 to 94.3%. There were 2128 sncRNAs identified. The categories of non-coding RNAs identified are in Figure 1A and File S3 and included miRNAs, snoRNAs, novel snoRNAs, tRNAs, snRNAs, and long non-coding RNAs (lncRNAs).

Preliminary Analysis of Small RNA-seq Data
To identify differential expression of sncRNAs in ageing, Illumina MiSeq was utilised. Summaries of raw, trimmed reads, and mapped reads to Equus caballus database can be found in File S2. Reads mapping percentages were between 92 to 94.3%. There were 2128 sncRNAs identified. The categories of non-coding RNAs identified are in Figure 1A and File S3 and included miRNAs, snoRNAs, novel snoRNAs, tRNAs, snRNAs, and long non-coding RNAs (lncRNAs). Clustering results shown as a heatmap (distance measure using Euclidean, and clustering algorithm using Ward) for the top 90 molecules.

Age-Related Differential Small Non-Coding RNA Gene Expression
There were no overall differences in the distribution of classes of sncRNAs in ageing. The effect of age on the expression of sncRNAs was weak and the separation of young and old samples was not clear. The 3D PCA plot ( Figure 1B) indicated that few, but very specific changes in the expression of sncRNAs were found, and samples from old donors were more variable than those from young. A heat map of hierarchical clusters of correlations among samples ( Figure 1C) depicts that the sncRNA expression of young and old groups are not very different.

Age-Related Differential Small Non-Coding RNA Gene Expression
There were no overall differences in the distribution of classes of sncRNAs in ageing. The effect of age on the expression of sncRNAs was weak and the separation of young and old samples was not clear. The 3D PCA plot ( Figure 1B) indicated that few, but very specific changes in the expression of sncRNAs were found, and samples from old donors were more variable than those from young. A heat map of hierarchical clusters of correlations among samples ( Figure 1C) depicts that the sncRNA expression of young and old groups are not very different.
There were 83 sncRNAs differentially expressed with age; six snoRNAs, 11 novel snoRNAs, three snRNAs, 31 lncRNAs, 27 tRNAs (p < 0.05) and five miRNAs (FDR-adjusted p < 0.05) ( Table 1). Data is deposited on NCBI GEO, accession; E-MTAB-8112. Table 1. Differentially expressed sncRNAs in ageing chondrocytes. Log 2 fold change values were derived with young as the reference group. A positive log 2 fold change equates to higher expression in old, whereas a negative log 2 fold change equates to lower expression in old. All are significant at p < 0.05, except for miRNAs that are significant with a false-discovery rate (FDR-)adjusted p value of p < 0.05.

Age Specific miRNA Interactome
To generate an age-specific miRNA interactome of the most likely miRNA-mRNA target pairs, analysis was performed to identify miRNA targets of the five differentially expressed miRNAs from this study. In IPA, the five miRNAs (miR-143, miR-145, miR-181b, miR-122, miR-148a) were paired with 351 protein coding genes differentially expressed in our previous RNA-seq study on ageing equine cartilage [31]. In File S4, the miRNA-mRNA pairings in the correct direction (miR increase and mRNA decrease or vice versa) are shown, including the target predictions and /or experimental validation in the respective database. Four of the five miRNAs interacted in IPA with 31 different differentially expressed target genes reflecting that miRNAs target many mRNAs. These mRNAs targeted by the miRNAs were used in IPA as network-eligible molecules and overlaid onto molecular networks based on information from the IPA Database. Networks for the four miRNAs (miR-148a, miR-122, miR-143, and miR-181b) were generated based on connectivity (Figure 2A-D). Interesting age-related features were determined from the gene networks inferred. Among the top canonical pathways were hepatic fibrosis (p = 1.51 × 10 −3 ), glycoprotein 6 (GP6) signalling (p = 5.67 × 10 −4 ), and osteoarthritis pathway (p = 2.93 × 10 −3 ). The top diseases and cellular functions associated with this network are shown in Table 2 and Figure 3. All IPA results are in File S5.

tRNA Fragment Changes in Equine Ageing
As we identified differentially expressed tRNAs in ageing equine chondrocytes, we undertook additional analysis to identify tiRNAs and tRFs, as increasing experimental evidence suggests their

tRNA Fragment Changes in Equine Ageing
As we identified differentially expressed tRNAs in ageing equine chondrocytes, we undertook additional analysis to identify tiRNAs and tRFs, as increasing experimental evidence suggests their functional roles in OA [27]. Figure 5 shows cumulative density of tRNA fragment length, alignment length, gene counts, and map start position. On the assessment of data variation using PCA, samples in group "old" did not scatter very closely. Differential expression analysis of all young versus all old did not identify any differentially expressed tRFs (FDR < 0.05). Therefore, we reanalysed the data with the old group divided into two subgroups based on PCA findings: 'old 1' (samples 6,7,8) and 'old 2' (samples 9 and 10) for further differential expression analysis. There were 81 differentially expressed tiRNAs/tRFs; 44 higher in 'old 2' and 37 lower in 'old 2' compared to young (File S6). functional roles in OA [27]. Figure 5 shows cumulative density of tRNA fragment length, alignment length, gene counts, and map start position. On the assessment of data variation using PCA, samples in group "old" did not scatter very closely. Differential expression analysis of all young versus all old did not identify any differentially expressed tRFs (FDR < 0.05). Therefore, we reanalysed the data with the old group divided into two subgroups based on PCA findings: 'old 1' (samples 6,7,8) and 'old 2' (samples 9 and 10) for further differential expression analysis. There were 81 differentially expressed tiRNAs/tRFs; 44 higher in 'old 2' and 37 lower in 'old 2' compared to young (File S6).

Human tRNA and tRF Profiles Compared to Equine tRNA and tRF Profiles
In both human and equine samples, 26 parent tRNAs were detected. Of these 26 tRNAs, 13 tRNAs were induced and 13 were reduced in cartilage from high grade OA compared to cartilage from low grade OA. In equine chondrocytes, 6 parent tRNAs were expressed higher and 20 parent tRNAs were expressed lower in old samples compared to young. Four parent tRNAs were induced in both high grade OA cartilage compared to low grade OA cartilage and in old versus young equine samples. Eleven parent tRNAs were reduced in high grade OA cartilage compared to low grade OA cartilage and in the old versus young equine samples ( Figure 6A).
In both human and equine samples, the tRF-5 fragment known as tRF293/294 and 10 tiRNAs were detected ( Figure 6B). In high grade OA compared to low grade OA cartilage, seven tiRNAs and tRF293/294 were induced and three tiRNAs were reduced. In old 1 versus young equine, five tiRNAs and tRF293/294 were induced and five tiRNAs were reduced. In old 2 versus young equine, two tiRNAs were induced and eight tiRNAs and tRF293/294 were reduced. Of these tRNA fragments, three tiRNAs and tRF293/294 were induced and tiRNA-5029-GlyGCC-3 was reduced in high versus low grade OA cartilage and in old 1 versus young equine. Two tiRNAs were induced and three

Human tRNA and tRF Profiles Compared to Equine tRNA and tRF Profiles
In both human and equine samples, 26 parent tRNAs were detected. Of these 26 tRNAs, 13 tRNAs were induced and 13 were reduced in cartilage from high grade OA compared to cartilage from low grade OA. In equine chondrocytes, 6 parent tRNAs were expressed higher and 20 parent tRNAs were expressed lower in old samples compared to young. Four parent tRNAs were induced in both high grade OA cartilage compared to low grade OA cartilage and in old versus young equine samples. Eleven parent tRNAs were reduced in high grade OA cartilage compared to low grade OA cartilage and in the old versus young equine samples ( Figure 6A). tiRNAs were reduced in high versus low grade OA cartilage and in old 2 versus young equine. In high versus low grade OA cartilage and old 1 versus young equine tiRNA-5029-GlyGCC was reduced.  In both human and equine samples, the tRF-5 fragment known as tRF293/294 and 10 tiRNAs were detected ( Figure 6B). In high grade OA compared to low grade OA cartilage, seven tiRNAs and tRF293/294 were induced and three tiRNAs were reduced. In old 1 versus young equine, five tiRNAs and tRF293/294 were induced and five tiRNAs were reduced. In old 2 versus young equine, two tiRNAs were induced and eight tiRNAs and tRF293/294 were reduced. Of these tRNA fragments, three tiRNAs and tRF293/294 were induced and tiRNA-5029-GlyGCC-3 was reduced in high versus low grade OA cartilage and in old 1 versus young equine. Two tiRNAs were induced and three tiRNAs were reduced in high versus low grade OA cartilage and in old 2 versus young equine. In high versus low grade OA cartilage and old 1 versus young equine tiRNA-5029-GlyGCC was reduced.

Discussion
Our study investigated the changing sncRNA landscape in ageing chondrocytes. Several risk factors exist that influence cartilage health and chondrocyte homeostasis. Among them, ageing is one of the leading risk factors contributing to cartilage-related diseases, such as OA [52]. Many studies have shown that ageing can affect cartilage in different ways, both at cellular and molecular level. Increased chondrocyte death, apoptosis, and a shift towards a catabolic profile have been observed in aged chondrocytes [53]. Additional age-related changes in articular cartilage include increased chondrocyte senescence [54], oxidative stress [55], and changes in the composition and structure of ECM [56]. Although the underlying molecular causes of these changes are not completely understood, it is hypothesised that aged chondrocytes respond differently to various stimuli, such as growth factors, [53,57] and demonstrate altered molecular signatures [6].
SncRNAs are a subset of epigenetic modifiers and their role in cartilage ageing has been studied increasingly in the last decade [6,31,58]. In the current study, we have used small RNA-seq to identify alterations in sncRNAs between young and old equine chondrocytes. The horse is a good model to study musculoskeletal ageing and disease as we could assess the whole joint for pathological perturbations during tissue collection. It is very challenging to source aged human cartilage that has no OA changes, whereas this is easily undertaken in equine samples. Moreover, the horse has been used as a model of OA and there has been significant research on equine joint anatomy and pathophysiology [59,60].
Within our set of differentially expressed sncRNAs, miRNAs are the best studied in musculoskeletal ageing and cartilage. In old chondrocytes, we identified two miRNAs with higher expression; miR-122 and miR-148a, and three miRNAs with lower expression; miR-143, miR-145, and miR-181b. Of these five miRNAs, all except miR-148a were validated in an extended cohort of young and old equine chondrocytes with qRT-PCR. MiR-122 has been researched extensively in liver [61,62], but its role in musculoskeletal ageing is less clear. MiR-122 was decreased in the serum and plasma of patients with osteoporosis, the most common age-related bone disease [63], but was significantly upregulated in senescent human fibroblasts [64] and was shown to upregulate p53 which is induced in senescence [65]. MiR-143 was downregulated in muscle satellite cells from old mice and primary myoblasts from old humans and mice [66]. In addition, circulating miR-143 was upregulated in young individuals following resistant exercise, but was downregulated in older individuals after resistant exercise [67]. MiR-145 was downregulated in old OA patients [68] as well as in experimental OA rat chondrocytes treated with tumour necrosis factor (TNF) [69] and, finally, miR-181b was downregulated in skeletal muscle of old rhesus monkeys [70].
To further investigate the potential role of the differentially expressed miRNAs identified in this study, we used IPA to combine them with the differentially expressed mRNAs from our previous equine cartilage study [31]. IPA miRNA 'Target Filter and Expression Pairing' identified 31 potential target genes. IPA core analysis of these genes revealed canonical pathways associated with cartilage physiology, such as role of chondrocytes in rheumatoid arthritis, OA-related pathways and bone morphogenic protein (BMP) signalling, all of which have been reported to change with ageing [56,71,72]. Moreover, top diseases and disorders linked to these genes, as identified by IPA, included skeletal and muscular disorders and connective tissue disorders. Of note, follistatin (FST) which was upregulated in old equine cartilage and was predicted by IPA as a target of the downregulated miR-143, was overexpressed in human OA chondrocytes [73] and canine OA cartilage [74] and was induced by telomere shortening [75], a process associated with ageing. Moreover, TNF ligand superfamily member 11 (TNFSF11), also known as receptor activator of nuclear factor kappa-B ligand (RANKL) was upregulated in old equine cartilage and was identified as a predicted target of the downregulated miR-181b. Higher expression of TNFSF11/RANKL, which correlated with bone loss, was reported in old C57BL/6 mice [76], rabbits with chronic antigen-induced arthritis [77], and in human high grade OA cartilage [78]. These results demonstrate the adverse effect of ageing on miRNA levels and their potential use as biomarkers or therapeutic targets for age-related musculoskeletal diseases.
Six snoRNAs were identified as differentially expressed due to ageing in chondrocytes. This conserved class of non-coding RNAs are principally characterised as guiding site-specific post-transcriptional modifications in ribosomal RNA [79] (canonical snoRNAs), but can also modify additional classes of RNAs including other snoRNAs, tRNAs and mRNAs; so called non-canonical snoRNAs [80,81]. Examples of age-related snoRNAs in equine ageing chondrocytes with canonical functions include snora71 and snord29. Novel non-canonical functions reported for snoRNAs including the modulation of alternative splicing [82], an essential involvement in stress response pathways [83] and the modulation of mRNA 3 end processing [84]. Similar to miRNAs, snoRNAs are emerging as important regulators of cellular function and disease development [15], related to their ability to fine-tune ribosomes accommodating changing requirements for protein production during development, normal function, and disease [85]. We have previously identified a role for snoRNAs in cartilage ageing and OA [16] and their potential use as biomarkers for OA [9]. Interestingly there was a reduction in snord113/114 in ageing chondrocytes which agrees with our previous findings in equine cartilage ageing [31]. We have also previously demonstrated that SNORD113 was reduced in ageing human knee cartilage but increased in OA [16] and increased in human anterior cruciate ligament [17]. SNORD113/114 are located in imprinted loci and may play a role in the evolution and/or mechanism of epigenetic imprinting. Although belonging to the C/D box class of snoRNAs which direct site-specific 2 -O-ribose methylation of substrate RNAs, they differ from other C/D box snoRNAs in their tissue specific expression profiles (including fibroblast, chondrocytes and osteoblasts) and the lack of known substrate RNA complementarity. This currently classifies them as orphan snoRNAs as they are not predicted to guide the 2 -O-ribose methylation but have novel, unknown roles [86]. Additionally, SNORD113 functions as a tumour suppressor in hepatic cell carcinoma by reducing cell growth, and it inactivates the phosphorylation of extracellular signal-regulated kinase (ERK) 1/2 and mothers against decapentaplegic homolog (SMAD) 2/3 in mitogen-activated protein kinase (MAPK)/ERK and transforming growth factor beta (TGF-β) pathways [87]. Together, our snoRNA findings indicate that age-related changes in chondrocyte snoRNAs could have important implications through both canonical and non-canonical snoRNA routes. This is the first study to detect tRNA and tRNA fragments in equine chondrocytes and to compare these findings with tRNA and tRNA fragments detected in human OA cartilage. For our tRNA data, old donors clustered into two groups and further analysis was undertaken with these subgroups. There were no apparent differences in these subgroups with regards to age or sex and scores were all zero. The parent tRNA Cys-GCA was found to be increased in both aged equine chondrocytes and high grade OA human cartilage samples. tRNA Cys-GCA levels have previously been reported to be increased in human chondrocytes induced with the cytokine interleukin (IL) 1 beta resulting in the production of the tRNA fragment tRF-3003a, a type 3-tRF produced by the cleavage of Cys-GCA. tRF-3003a has been shown to post-transcriptionally regulate the Janus Kinase 3 (JAK3) expression through sequence complementarity via the Argonaute (AGO) / RNA-induced silencing complex (RISC) in human chondrocytes [27].
Many of these cytokines are known to play important roles in synovial inflammation during OA pathogenesis [88,89].
The tRNA fragments detected in equine samples that matched with human samples consisted of 5 tiRNA and tRF-5 fragments. Many of the equine fragments detected did not fall into the classical tRF-3, tRF-5, or tiRNA size ranges and instead may likely be i-tRFs, which are internal to the respective tRNA and can straddle the anticodon loop. 5 tiRNA, 3 tiRNA, tRF-3, and tRF-5 fragments were detected in the low-grade OA cartilage and in the high grade OA cartilage. In our human studies, base modifications found on tRNAs and tRF/tiRNA fragments that would normally block reverse transcription were removed and this may account for some of the differences found between the equine and human tRNA/tRF profiles.
The importance of modulation of tRNA levels and tRNA fragments in articular cartilage homeostasis remains an unexplored area. This is the first evidence that aged equine samples have changes in the expression of specific tRNAs and tRFs when compared to young equine samples. We report for the first time several 5 tiRNAs, such as tiRNA Glu-TTC and tiRNA His-GTG, were induced in old compared to young equine chondrocytes and in high grade compared to low grade human OA cartilage. Previous reports have shown that 5 tiRNAs can be produced by cell stress in mammalian cells and these 5 tiRNA half fragments may have a role in inhibiting cell translation and could be involved in stress granule formation [90]. Further studies are required to find the mechanism by which these fragments are produced and whether the changes in the profile of fragments found in old compared to young equine chondrocytes or high compared to low grade OA cartilage potentially contribute to the development of OA.
We are aware our study has a number of limitations. The effect of ageing between young and old equine chondrocytes was small on the differential expression of sncRNAs. However, it is likely that we are therefore interrogating highly specific changes that are age dependent. Furthermore, we cannot rule out changes related to the use of chondrocytes instead of cartilage tissue. Even though chondrocytes of low passage were used, collagenase digestion and plating of cells could have affected their phenotype and gene expression. The choice of chondrocytes over cartilage tissue was made based on RNA from cartilage tissue is of low quality (in our hands) and heavily contaminated with proteoglycans [91], usually making it challenging for sequencing. In our previous study interrogating snoRNAs in human knee cartilage ageing and OA [16] we utilized cartilage tissue as opposed to chondrocytes. The old non-OA cartilage was derived from the grossly normal condyle of OA joints. In the study we found a number of age and OA specific snoRNAs, but we could not define the cartilage as truly normal as it was derived from an OA joint. Nevertheless, we identified two age related, no species specific snoRNAs that were DE in both our studies, snord113 and snord29. The use of the Illumina MiSeq platform would have contributed to the low number of differentially expressed sncRNAs as it offers less depth coverage compared to other platforms, such as the HiSeq platform. Finally, we have used a relatively small number of samples per group. Given the use of primary cells and the degree of variability observed, especially for the old group, the inclusion of five samples pre group may have contributed to the small number of differentially expressed sncRNAs in ageing equine chondrocytes on our study.

Conclusions
For the first time we have described, using unbiased methods, the effect of ageing on the expression of sncRNAs in equine chondrocytes. We detected variable classes of sncRNAs (the small non-coding RNAome) in young and old chondrocytes, which are differentially abundant, indicating that there are multiple levels of epigenetic control in cartilage and chondrocyte ageing. Among them, there were miRNAs, which are predicted to play a role in the development of the musculoskeletal system and in skeletal disorders. In addition, the current study is one of the few studies that have investigated tRNAs and tRNA fragments, in an attempt to uncover novel molecular signatures in aged and diseased chondrocytes/cartilage that could be useful in the future as therapeutic targets. Further research is needed to elucidate the role and function of these molecules and their potential link to disease.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/21/16/ 5675/s1. File S1. .xls, Small RNA assay details; File S2. .xls, MiSeq read mapping details; File S3. .xls, Types of small non-coding RNAs identified in MiSeq data; File S4. .xls, miRNA-mRNA pairings in the correct direction; File S5. .xls, IPA analysis summary; File S6. .xls, differentially expressed tRNA halves/tRFs Author Contributions: M.J.P., P.B., T.J.M.W., J.A.G. and T.M.H. conceived the study and designed the work; P.B., P.D., J.A.G., M.J.P., Y.A.K. acquired and analysed the data, Y.F., X.L., interpreted the data; all authors have drafted the work. All authors have approved the submitted version. All authors have agreed both to be personally accountable for the author's own contributions and to ensure that questions related to the accuracy or integrity of any part of the work, even ones in which the author was not personally involved, are appropriately investigated, resolved, and the resolution documented in the literature. All authors have read and agreed to the published version of the manuscript.