Serum snoRNAs as biomarkers for joint ageing and post traumatic osteoarthritis

The development of effective treatments for the age-related disease osteoarthritis and the ability to predict disease progression has been hampered by the lack of biomarkers able to demonstrate the course of the disease. Profiling the expression patterns of small nucleolar RNAs (snoRNAs) in joint ageing and OA may provide diagnostic biomarkers and therapeutic targets. This study determined expression patterns of snoRNAs in joint ageing and OA and examined them as potential biomarkers. Using SnoRNASeq and real-time quantitative PCR (qRT-PCR) we demonstrate snoRNA expression levels in murine ageing and OA joints and serum for the first time. SnoRNASeq identified differential expression (DE) of 6 snoRNAs in young versus old joints and 5 snoRNAs in old sham versus old experimental osteoarthritic joints. In serum we found differential presence of 27 snoRNAs in young versus old serum and 18 snoRNAs in old sham versus old experimental osteoarthritic serum. Confirmatory qRT-PCR analysis demonstrated good correlation with SnoRNASeq findings. Profiling the expression patterns of snoRNAs is the initial step in determining their functional significance in ageing and osteoarthritis, and provides potential diagnostic biomarkers and therapeutic targets. Our results establish snoRNAs as novel markers of musculoskeletal ageing and osteoarthritis.

Osteoarthritis (OA) is an age-related musculoskeletal disease and a common cause of chronic disability worldwide 1 . In addition it is a significant contributor to both individual and socioeconomic burden and the number of disability adapted life years globally 2 . If the deterioration in musculoskeletal health and development of OA can be identified and treated early serious life impairment may be abrogated. Ageing is the time-dependent reduction of functional capacity and stress resistance, associated with an increased risk of morbidity and mortality. The joint and its articular cartilage is particularly affected by ageing 3 . There is evidence that the rate of ageing, that is the 'biological age' , differs significantly between individuals' actual age in years (i.e. the 'chronological age'). Defining markers of joint ageing may enable a prediction of the risk of onset of OA, enabling early intervention. OA is characterised by a non-symptomatic, pre-radiographical phase that if identified would allow earlier diagnosis. However radiographic changes are only evident later in disease progression. Magnetic resonance imaging techniques have been developed for early-stage evaluation of cartilage damage in OA but are expensive and contraindicated in some individuals.
The development of effective treatments for OA and the ability to predict disease progression has been hampered by the lack of substantive biomarkers, able to demonstrate pathological disturbances preceding identifiable tissue alterations. Others have attempted to identify products of tissue turnover in serum and synovial fluid (reviewed 4 ). This has been challenging due to patient and disease heterogeneity and dilution effects either by tissue fluids or with similar products from other joints or diseases. In addition, the variability of antibody assays has been problematic.
SnoRNAs are a class of evolutionary conserved non-coding small guide RNAs of which the majority direct the chemical modification of other RNA substrates, including ribosomal RNAs and spliceosomal RNAs. In addition, some snoRNAs are involved in the regulation of alternative splicing and post-transcriptional modification of OARSI scoring of histological sections of mouse knee joints. For histology, as the total knee joint was used in the study for RNA extraction, joints were collected (into 4% paraformaldehyde) from additional equivalent aged and treated young (n = 8), old (n = 4); sham (n = 5); and DMM (n = 6) mice in order to evaluate the extent of OA. The procedure, surgeon and duration of the studies were identical. Knees were decalcified in 0.5 M ethylenediaminetetraacetic acid (pH 7.4) for 4 weeks at 4 °C and coronally embedded in paraffin. Sectioning, Safranin-O Fast-Green staining and histological scoring (defined as the severity and extent of OA) was undertaken on a scale from 1 to 6 by two blinded independent observers using the OARSI histopathology initiative 18 . All four quadrants of the section (medial tibial plateau, lateral tibial plateau, medial femoral condyle, lateral femoral condyle) were scored individually and added for each histological section. For statistical analyses mean summed score values of joints of 3-5 section per knee joint at 4 depths throughout the joint was determined (thus a maximum score of 24 was possible). Inter-observer variability was calculated using Cohen's Kappa statistics using an online software tool: (http://www.statstodo.com/CohenKappa). RNA isolation, RNA-Seq analysis, cDNA library preparation and sequencing. Total RNA was isolated from equal weights of joints and 500 μ l serum using miRNeasy or RNeasy Serum kits respectively with DNase treatment (all Qiagen, Crawley, UK) to remove residual gDNA. Total RNA integrity (RIN) was confirmed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, USA). Ribosomal RNA was depleted using the Ribo-Zero ™ rRNA Removal Kit (Epicentre, Madison, USA). From 41 samples 100 ng of rRNA-depleted RNA was submitted for library preparation using NEB small RNA library kit (New England Biolabs (NEB), Ipswich, USA). To reduce workflow bias we used tobacco acid pyrophosphatase (Epicentre, Madison, USA) to remove potential 5′ caps found on some snoRNAs. Samples were amplified for 15 cycles, mixed into 3 pools, and size selected. The size-selected material was purified with Ampure beads (Agencourt, Beckman-Coulter, High-Wycombe, UK). SnoRNA sequencing was undertaken on the Illumina HiSeq 2000 platform (Illumina, San Diego, USA) using 100 base paired-end reads.
SnoRNASeq data analysis. Sequence data measured from 5 lanes of an Illumina HiSeq2000 were processed through a number of steps to obtain snoRNA expression values. The processes include basecalling and de-multiplexing of indexed reads using CASAVA version 1.8.2 19 ; adapter and quality trimming using Cutadapt version 1.2.1 20 and Sickle version 1.200 to obtain fastq files of trimmed reads; aligning reads to Ensembl GRCm38.77 mouse genome reference sequences which contains 1,555 annotated snoRNA features using Bowtie2 21 version 2.0.10 with option "-very-sensitive-local"; counting aligned reads against snoRNA features using THSeq-count. The count values were used as snoRNA expression measurements for the DE analysis.
Scientific RepoRts | 7:43558 | DOI: 10.1038/srep43558 DE analysis was performed in R environment using package edgeR 22 . 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; 3-D PCA plots were generated using R function in package plot3D) and correlation analysis; handling library size variation respectively for joint samples and serum samples through data normalisation; formulating data variation using negative binomial distributions; modelling data using a generalised linear model; computing log 2 Fold Change (logFC) values for required contrasts based on model fitting results through contrast fitting approach, assigning P-values to logFC values by LR 23 testing; dealing with the effects of multiple tests using FDR approach to obtain FDR adjusted P-values; and defining significantly DE snoRNAs as those with FDR-adjusted p-value < 5%. Sequence data have been submitted to National Centre for Biotechnology Information Gene Expression Omnibus (NCBI GEO); E-MTAB-4878. RNA isolation, poly(A) cDNA synthesis and snoRNA qRT-PCR. qRT-PCR of snoRNAs was performed 24 . Total RNA was isolated using the mirVana kit. Isolated RNA samples were polyadenylated at 37 °C for 60 minutes in a 50 μ L reaction volume containing 1 μ g RNA and 1.5 U poly(A) polymerase (NEB, Ipswich, USA). 500 μ L lysis binding buffer was added. Then, an equal volume of acid-phenol:chloroform was added, vortexed and samples were centrifuged for 10 minutes and the aqueous phase removed. The poly(A)-tailed total RNA was extracted using the filter cartridge provided by the mirVana kit. To generate poly(A) cDNA, 500 ng poly(A)-tailed RNA and 250 ng RTQ primer (Eurogentec, Seraing, Belgium) ( Table 1) were mixed in a 26 μ L reaction volume, incubated at 65 °C for 10 minutes and annealed at 4 °C for 20 minutes. Reverse transcription was performed with 200 U M-MLV reverse transcriptase, 20 U RNAsin (both Promega, Southampton, UK), 2 μ L dNTP mix (10 mM each; Eurogentec, Seraing, Belgium) and 8 μ L 5x M-MLV buffer (Promega, Sothampton, USA) in a total reaction volume of 40 μ L at 50 °C for 60 minutes. The reverse transcriptase was inactivated at 70 °C for 15 minutes. Finally, 1.5 U of RNAse H (NEB, Ipswich, USA) was added to remove small RNAs. A snoRNA-specific forward primer and a universal reverse primer (RTQ-UNIr, matched to the Tm of each individual snoRNA) were used for the amplification of each snoRNA target (all Eurogentec, Seraing, Belgium) ( Table 1). For each cDNA sample a mix was prepared with Mesagreen qPCR Mastermix Plus for SYBR Green (Eurogentec, Seraing, Belgium) and 300 nM forward and reverse oligonucleotides. An ABI-7300 Detection System was used for amplification using the following protocol: denaturation at 95 °C for 5 minutes, followed by 50 cycles of DNA amplification (15 seconds 95 °C and 45 seconds annealing at 62-68 °C). The annealing temperature was optimized for each snoRNAs target. Serially diluted standard curves were utilized to quantify snoRNA expression and data was normalized to a validated housekeeping snoRNA (joint: U2, young-old serum: SNORD85, old sham-old DMM serum and equine serum: U6).
Validation of SNORD116 as a marker of OA in equine serum. We determined the reproducibility of the expression of SNORD116 in OA serum of another species. There are well-published studies on the application of metacarpophalangeal (MCP) joint OA changes in the horse, a joint with similarities to the human knee joint 25 . We studied equine serum from normal and MCP OA horses collected from eight normal (mean age ± standard deviation 5.3 ± 2.1 years) and four OA (7.5 ± 1.0 years) castrated male thoroughbred horses at post-mortem. Samples were collected under the regulations of the Hong Kong Jockey Club with owner consent and stored at − 80 °C. OA diagnosis was based on histological (modified Mankin) 26 and synovitis scoring 27 of MCP joint tissues. Statistical analysis. For statistical evaluation of histological scoring non-parametric Mann-Whitney-U test was used. Inter-observer agreement of histological scoring systems was calculated using Cohen's kappa coefficient (www.statstodo.com/CohenKappa_Pgm.phpl). qRT-PCR data was log-transformed prior to statistical evaluation with an independent samples t-test. Statistical evaluation was performed between young and old or old sham versus old DMM using GraphPad Prism 5 (San Diego); p-values are indicated.

Results
OARSI scoring of joints. OARSI scoring of joints (mean ± 95%CI) for young and old were 0.5 ± 0.  Fig. 2A. This indicates that major disease responses in the data were contributed from serum samples. A read length distribution graph was generated to highlight the constitutional difference of non-coding RNAs in joint and serum samples at 100 bp or below (Fig. 2C). For this one serum and one joint sample with approximately the same library size were chosen. The plots of the read length distribution for joint and serum Therefore the corresponding reads are generally well annotated miRNAs. Another peak is noted at around 100 bp, which is consistent with non-coding RNAs whose length is a hundred bp and over. Validation of SnoRNASeq data using qRT-PCR. Very little is known about the roles of snoRNAs as biomarker for disease or about the functional roles of snoRNAs in mammalian cell biology. We thus could not define functional criteria to select snoRNAs for validation based on prior knowledge. To stay unbiased we thus selected snoRNAs from the RNAseq data that displayed a moderate fold significant expression difference, higher fold significant expression difference, non-significant expression difference and selected box H/ACA (SNORAs) as well as box C/D (SNORDs) snoRNAs for validation. Levels of candidate snoRNAs for further qRT-qPCR analysis were determined using the original RNA from all donors used to perform the SnoRNASeq experiment. There was good concordance between SnoRNASeq and qRT-PCR platforms ( Table 2 and Fig. 3). SNORD88   was significantly decreased in young versus old joint (Fig. 3A), while SNORA73 was validated to be increased in young versus old joint (Fig. 3A). In agreement with its absence in the DE group (Table 2), SNORA30 was not differentially expressed in young versus old joint (Fig. 3A) and SNORD88 was not differentially expressed in sham versus DMM joint (Fig. 3A). SNORA31, SNORA28, SNORD23 and SNORA73 were confirmed to be significantly increased in young versus old serum (Fig. 3B) and SNORD116, SNORA64 and U3 were significantly increased in sham versus DMM serum (Fig. 3C). SNORD46 was confirmed to be significantly decreased in DMM serum as compared to sham serum (Fig. 3C).

SNORD116 in equine serum in OA.
To confirm increased SNORD116 in OA serum levels in a different species using qRT-PCR, we measured SNORD116 in equine serum samples. Normal equine donors had a Mankin's score 1.25 ± 0.9 (mean ± 95%CI) and OA donors 7.75 ± 7.6. Synovial membrane from normal donors had a synovitis score 27 of 3.25 ± 2.3 and OA donors 3.25 ± 3.2. There was a significant increase in serum expression of SNORD116 in OA compared to normal donors (p = 0.0010) (Fig. 4).

Discussion
SnoRNAs are emerging as important regulators of cell functions, such as alternative splicing 28 , metabolic stress 29 and development of disease; in cancer 6 , Prader-Willi Syndrome and autism 5 . Through expression profiling of snoRNAs using deep-sequencing we reveal novel molecular features relating to joint ageing and OA. Whilst gene expression has been evaluated in animal models of OA, including the rat anterior cruciate transection 30 and meniscal tear models 31 , and mouse DMM model 32 these were primarily interrogating protein-coding genes. Furthermore, apart from the final study these experiments evaluated a single tissue, primarily articular cartilage. In musculoskeletal ageing single tissues have been investigated 13,14 . OA is a process that involves the whole joint as an organ; therefore we undertook our analysis on whole mouse joints, which included cartilage, meniscus, subchondral bone, and joint capsule with synovium.
Mice are considered skeletally mature at around 3-months-old (approximate equivalent of a teenaged human) while a 12-month-old mouse would signify a 40 to 50 year human 33 . Thus, to investigate the effects of joint age we used 8-month-old equivalent to 25 to 28 year human (referred to as young) and 24-month-old mice (referred to as old). To study the development of post-traumatic OA we measured OA severity histologically and analysed snoRNAs expression in joints and serum from 24-month-old mice following DMM. Limitations of this study were that we did not additionally undertake the DMM model in young mice. Additionally the 'old' age in young versus old contrast was not the same age as 'old' age in sham versus DMM. OARSI scores that were observed in the sham group (DMM surgery cohort; 24 months old) were lower than the old group (young versus old cohort; 18 months old). At this moment we cannot explain this difference, but different housing and environments of the two cohorts may have influenced this. The DMM model is a post-injury model in which the histological lesions within the affected joint are similar to those observed in human OA 17 . Mild OA-like pathology was present in the old and sham mice implying that mice at this age are in the early stages of acquiring naturally occurring OA comparable to studies of human knees at the equivalent age of approximately 40 years-old 34 . However the effect of the DMM model exacerbated these changes.
Detection of snoRNAs in serum has previously been demonstrated 12,35,36 . SnoRNAs are serum stable, and although normally resident in the nucleolus it is thought in serum they are present as unidentified protein complexes 12 . It is however not clear whether disease-associated RNAs detected in the circulation result from local tissue disturbances and cell death, or whether they are actively locally secreted via exosomes or microvesicles or are a systemic response upon local tissue damage 8,37,38 . This may even depend on the specific pathology or specific RNA species. Determining the average read-length distribution of representative joint and serum samples (Fig. 2C) revealed an unexpected finding. When looking at the read-length distribution the serum reads specifically contain a large peak of RNAs with a length of approximately 30 nt. The combined realisation that quite a large portion of the reads could not be mapped back to the used genome database, that a large portion of the piRNA class of non-coding RNAs lacks annotation in the used genome database and that piRNAs are typically 30 nt long, makes it is reasonable to think that this explains the presence of the large 30 nt peak specifically in serum. Only recently the presence of piRNAs in human blood has been described 39 and this significant difference in the constitution of small RNAs for joint and serum samples is a phenomenon calling for further investigation.
In the young versus old serum (27 snoRNAs) more snoRNAs were significantly up-regulated compared to the old sham versus old DMM serum (18 snoRNAs) indicating ageing per se had a greater effect on differential snoRNA presence in serum than OA. This could be due to different 'old' ages in the comparisons, but also maybe expected as in ageing serum snoRNAs will be from many tissues whereas in the DMM model we are most likely highlighting primarily OA joint-related snoRNAs. As we were studying the whole joint this may be due to varying expression of snoRNAs in joint tissues, as it is known that there is tissue-specific snoRNAs expression 40 . We are unable to determine the amount each tissue type contributed to overall expression of snoRNAs. Each tissue will vary in its cellularity and hence RNA and snoRNAs content. Thus whilst a novel aspect of this study was that snoRNAs were extracted from the multiple tissues that form the joint, our approach may be less sensitive in detecting snoRNAs that change in a single tissue. However it has the advantage of determining snoRNAs that could be more globally implicated in OA. Despite the potential limitations, we identified a number of potentially interesting snoRNAs for future studies.
SNORA73 was increased in old joint and serum (Table 2, Fig. 3) and represents a potential joint 'biological ageing' marker. A reliable measurement of the state of ageing and a prediction of the risk of the onset of morbidity for chronic age-related diseases such as OA would be beneficial. Such a strategy could serve as a measure of 'biological' age and predict an age-related biological response more accurately than chronological age. SNORA64 was increased and SNORD46 was decreased in DMM serum, but both were not DE in young versus old serum, indicating these snoRNAs as possible OA markers. SNORD18 was increased in serum both in ageing and following DMM (Table 2) signifying that they are affected in ageing and also OA. It is difficult to speculate how much this snoRNA changes in age-related OA versus age from this study. Histology identified mild increase in OARSI score with age and the level of OA changes in the DMM model was mild. It would be beneficial to investigate this snoRNAs further in tissues and serum in more severe OA.
An interesting finding was an increase in SNORD38 in ageing joint. Zhang et al. 12 demonstrated a strong association between serum levels of SNORD38 and severe cartilage damage, in anterior cruciate ligament (ACL) injury, enabling distinction between ACL injury patients from normal donors. They were unable to determine  age effect of SNORD38 in serum from normal donors as it was undetectable. One year post-surgery there was no relationship between donor age and serum SNORD38. Our study demonstrated an increase in SNORD38 in mouse joints (but not serum) with age. The human study did not assess tissue snoRNAs and the serum samples were primarily from middle-aged donors. Our old mice group represents an equivalent older age than that in the human study, which could contribute to this disparity.
The most DE snoRNA in DMM serum was SNORD116. Additionally we demonstrated an increase in serum SNORD116 in horses with MCP OA. We have previously identified SNORD116 as increased in OA compared to normal human cartilage in an array study 15 . A loss of SNORD116 is a significant contribution to the aetiology of the neurodegenerative genetic condition Prader-Willi syndrome (PWS) 41 . This paternally imprinted disorder results in developmental delay and genetic obesity due to hyperphagia 42 . Clinical signs include short stature and low bone mineral density 43 . In a recent mouse transgenic study the loss of PWS critical region (including SNORD116) resulted in reduced bone mineral density (BMD), delayed skeletal development and reduced bone size and osteoblastic suppression 44 . Humans with OA have an increased BMD in affected joints 45 . Thus our findings not only elucidate a potential marker of OA but a snoRNA with a potential role in the pathogenesis of OA.
Previously we have identified DE snoRNAs in ageing cartilage 13 and tendon 14 and together with results in this study we propose that in musculoskeletal tissues snoRNAs potentially modulate the ageing process as previously described 46 . While determining the transcriptomic signature of ageing equine cartilage 13 we found the differential expression of a number of snoRNAs associated with ageing. When comparing the differentially expressed snoRNAs between the equine and this mouse study it is important to realize that in our previous equine study we did not analyse serum or the whole joint, but specifically the articular cartilage. Overlapping differential snoRNA expression between studies was identified for SNORD113, SNORA53, SNORA48 and SNORA5. SNORA53 was decreased in old equine cartilage, but increased in old mouse serum. SNORA48 was decreased in old equine cartilage, and decreased in old DMM mouse serum. SNORA5 was increased in old equine cartilage and increased in old mouse serum. Possibly the best consistency was found for SNORD113, which was decreased in old mouse joint (and old DMM mouse joint) and equine cartilage. The apparent cross-species conservations of differentially expressed snoRNAs in ageing and OA strengthen our belief that snoRNAs could indeed be used as biomarkers. Further analysis of snoRNA expression profiles and detailed genetic studies will give new insights into novel molecular networks in musculoskeletal ageing and common mechanisms in ageing and age-related diseases such as OA.
In mammalians the majority of snoRNAs are encoded within the introns of protein coding or non-coding genes; host-genes 47 . There is evidence that genes which host snoRNAs might contribute to the aetiology of cancer through regulation of cell homeostasis and cancer biology (reviewed 6 ). Potentially both the host-gene and the snoRNAs encoded within them may be important in different situations. For example growth-arrest-specific-5 (GAS-5) (hosts ten C/D box snoRNAs 48 ), a non-coding RNA which accumulates in growth arrested cells, regulates cell death and proliferation by acting as a decoy hormone response element for glucocorticoid receptors thereby inhibiting gene upregulation by activated glucocorticoid receptors 49 . Interesting snoRNA host-genes were identified in this study including transforming growth factor β regulated gene 4, sorting nexin 5 and collagen type 1 (with roles in joint homeostasis). Therefore, an alteration of snoRNA expression may result from changes in transcriptional activity of the host-gene related to joint homeostasis or disease.

Conclusion
Our results implicate specific changes in snoRNA abundance in joint ageing (SNORD88 and SNORD38 were respectively decreased and increased) and OA suggesting the potential use of snoRNAs such as SNORA73 and SNORD23 as a novel biomarker for joint ageing, SNORA64, SNORD46 and SNORD116 for OA, SNORD18 for ageing and OA.