SnoRNA signatures in cartilage ageing and osteoarthritis

Osteoarthritis presents as a change in the chondrocyte phenotype and an imbalance between anabolic and catabolic processes. Age affects its onset and progression. Small nucleolar RNAs (SnoRNAs) direct chemical modification of RNA substrates to fine-tune spliceosomal and rRNA function, accommodating changing requirements for splicing and protein synthesis during health and disease. Articular cartilage from young, old and OA knees was used in a microarray study to identify alterations in snoRNA expression. Changes in snoRNAs in osteoarthritis-like conditions were studied in chondrocytes using interleukin-1 and osteoarthritic synovial fluid. SNORD26 and SNORD96A knockdown and overexpression were undertaken using antisense oligonucleotides and overexpression plasmids. We identified panels of snoRNAs differentially expressed due to ageing (including SNORD96A, SNORD44) and osteoarthritis (including SNORD26 and SNORD116). In vitro experiments using osteoarthritis-like conditions affected snoRNA expression. Knockdown or overexpression of SNORD26 or SNORD96A resulted in changes in chondrogenic, hypertrophic, rRNA and osteoarthritis related gene expression. We demonstrate that snoRNA expression changes in cartilage ageing, and osteoarthritis and in osteoarthritis-like conditions, and when the expression of these snoRNAs is altered this affects chondrogenic and hypertrophic gene expression. Thus, we propose an additional dimension in the molecular mechanisms underlying cartilage ageing and osteoarthritis through the dysregulation of snoRNAs.

Osteoarthritis (OA) is the most common degenerative joint disease characterised by chondrocyte alterations and an irreversible loss of extracellular matrix (ECM) 1 leading to biomechanical failure. From a biochemical perspective, OA is characterized by uncontrolled synthesis of ECM degrading enzymes resulting in active cartilage breakdown. Chondrocytes are significant secretory cells enabling the synthesis and maintenance of the protein-rich ECM. Changes in the chondrocyte phenotype in OA are considered fundamental pathological mechanisms 2 . During ageing (and joint disease), the chondrocyte's homeostatic balance is disrupted and the rate of collagen and proteoglycan loss from the matrix may exceed deposition rate of newly synthesized molecules (reviewed 3 ) resulting in increased OA risk.
To ensure continuous ECM deposition it is essential chondrocytes control the number and quality of its ribosomes. Ribosomes are cellular nanomachines equipped for conversion of genetic information encoded in mRNAs into proteins 4 . Human 18S, 5.8S and 28S rRNAs contain at least 110 individual 2′O-ribose methylated and 100 pseudouridylated nucleotides 5 . These post-transcriptional modifications (PTMs) fine-tune translational characteristics of the ribosome. Positioning of these modifications is undertaken by small nucleolar RNAs (snoR-NAs) 6 . SnoRNAs are (mainly) intron-derived non-coding RNAs of approximately 50-250 nucleotides long with an important task in the PTMs of rRNA substrates and are classified into box C/D and box H/ACA snoRNAs. Work undertaken in zebrafish suggested decreased snoRNA expression (SNORD26, SNORD44 and SNORD78) reduced the snoRNA-guided methylation of the target nucleotides and that impaired rRNA modification, at a solitary site, led to critical morphological defects and embryonic lethality. The researchers alluded that rRNA modifications play an fundamental role in vertebrate development 7 . This was recently emphasised by a study showing important dynamics in rRNA ribose methylation and snoRNA expression during mouse development 8 . www.nature.com/scientificreports/ Additionally there are many orphan snoRNAs for which no target has been identified, however emerging evidence shows that orphan snoRNAs (and some canonical snoRNAs) fulfill non-canonical functions in alternative splicing 9 , modification of other RNAs 10 , a source for miRNAs 11 , regulation of metabolic stress 12 and gene expression 13 . Furthermore, internal 2′-O-methylation modifications (Nm) are guided by snoRNAs and these Nm sites can regulate mRNA and protein translation 14 .
Atypical snoRNAs expression has been implicated in some disease processes 15 including prostate, lung and breast cancer (reviewed 16 , neurodegenerative disease; Prader-Willi Syndrome 17 and viral infection 18 ). We have identified aberrant expression of groups of snoRNAs in ageing cartilage 19 and diseased tendon 20 . Whilst others have noted age-associated changes in snoRNAs in C. elegans 21 . Recently through snoRNA profiling using deep-sequencing we identified differentially expressed snoRNAs relating to joint ageing and OA 22 . Others have identified SNORD38 and SNORD48 as potential non-age-dependant serum biomarkers for OA progression following cruciate ligament injury 23 .
Expression of snoRNAs in human articular cartilage has not been explored, neither is it known whether snoRNA expression changes in cartilage due to ageing or OA, and whether the expression of individual snoRNAs can significantly influence the chondrocyte phenotype. In this study we therefore tested the hypothesis that with age and OA there is a dysregulation of expression and function of specific snoRNAs. Our study identifies a novel set of potential OA therapeutic targets.

Methods
Study design. The ethical committee approved all methods, which were undertaken following the relevant guidelines and regulations of the University of Maastricht Medical Centre. The studies were completed with prior approval from the University of Maastricht Medical Centre Institutional Review Board. Informed consent was obtained from all subjects. Human cartilage was harvested from male knees at the University of Maastricht Medical Centre. Cartilage was collected at total knee arthroplasty surgery from lateral [protected (P)] and medial [unprotected (U)] femoral condyles (MUMC + IRB; 2017-0183) or from the medial side of the lateral femoral condyle following anterior cruciate ligament repair surgery from young donors (MUMC + IRB; 14-4-038). We use the term 'protected' to denote this is the grossly unaffected cartilages within the knee of the donors used and do not insinuate that the lateral condylar cartilage is not subject to trauma or response to load.
Assessment of disease severity. Old specimens came from patients with a diagnosis of OA on preoperative knee radiographs using Kellgren and Lawrence (KL) Scoring 20 . The Outerbridge Scoring System 21 was applied to U and P samples.
Histological analysis of OA severity. Cartilage biopsies taken adjacent to those used for microarray from the protected and unprotected areas of femoral condyles were collected for histology. Biopsies were fixed in 4% paraformaldehyde and embedded in paraffin wax. 4 μm longitudinal sections were cut on a microtome (Leica Biosystems, UK) and stained with Haematoxylin/Eosin (Leica Biosystems, UK) and Safranin-O/Fast Green. Sections were scored for OA severity on a scale from 1 to 11 by two independent blinded observers using a modified Mankin scoring system 22 . Inter-observer variability was calculated using Cohen's Kappa statistics with online software (https ://www.stats todo.com/Cohen Kappa ).
RnA extraction, microarray. RNA was extracted from cartilage 23 . Twenty-nine Affymetrix miRNA-4.0 microarrays were used. Total RNA samples were quantitated using a Nanodrop spectrophotometer (NanoDrop Technologies). 500 ng of total RNA was labelled using the Affymetrix Flash-Tag Biotin HSR RNA labelling kit according to manufacturer's instructions. Following Flash-Tag labelling the biotin-labelled samples were stored at − 20 °C prior to hybridisation onto Affymetrix GeneChip miRNA 4.0 for 17.5 h at 48 °C 60 rpm in an Affymetrix hybridisation oven 645. Following hybridisation the arrays were washed using Affymetrix Hybridisation wash and stain kit on the GeneChip Fluidics station 450 using fluidics script FS450_0002, and scanned using the Affymetrix GeneChip scanner 3,000 7G. CEL files were generated using the Affymetrix GeneChip Command Console Software, and Expression Console software was used to QC.

RNA isolation, poly (A) cDNA synthesis and snoRNA qRT-PCR. For validation of microarray
findings an independent cohort was used. QRT-PCR of snoRNAs was performed 24 . Total RNA was isolated as described. RNA samples were polyadenylated as previously described 24 . A snoRNA-specific forward primer and a universal reverse primer (RTQ-UNIR, matched to the Tm of each snoRNA) were used for the amplification of each target (all Eurogentec, Seraing, Belgium) (Primer sequences are in Supplementary File 1). SnoRNA qRT-PCR was undertaken and the annealing temperature was optimized for each snoRNA 24 . Standard curves were used to quantify snoRNA expression and data was normalised to a validated housekeeping snoRNA. Steadystate transcript abundance of potential endogenous control genes was measured in the microarray data. Assays for SNORD63, SNORD28, SNORA28, SNORA31, U6, U2, and miR6786 were selected as potential reference genes because their expression was unaltered in the microarrays. Stability of this panel of genes was assessed with a gene stability algorithm 25 using genormPLUS (Biogazelle, Zwijnaarde, Belgium) 26  (Sigma-Aldrich) was used at 10 ng/ml. OA SF (n = 10) was perioperatively collected from the OA patients undergoing total knee arthroplasty and centrifuged to remove cells/debris, aliquoted and stored at − 80 °C. For experiments, SF from 10 patients was pooled and applied on cells in a 20% (v/v) concentration. Details of all donors are in Supplementary file 2. The expression of SNORD116, SNORD96A, SNORD26, SNORD33, SNORD44, SNORD95 and SNORD98 were measured in HAC with qRT-PCR using methods described previously. A panel of chondrocyte phenotype genes was measured in these samples; SOX9, ACAN, COL2A1, RUNX2, COL10A1, MMP13 and COX2. To determine the effect of oxygen tension (5% or 20%) and serum (serum-free or 10% foetal calf serum (FCS)) on snoRNA gene expression non-OA HACs were plated at 30,000 cells/cm 2 in triplicate per donor (passage 1, n = 3 males; mean ± standard deviation age, 25 ± 20), using methods above. RNA was harvested 48 h after plating. The extent of chondrocyte dedifferentiation due to passage on snoRNA expression was assessed using equine articular chondrocytes isolated from grossly normal metacarpophalangeal cartilage (n = 5; mean ± standard deviation age, 4 ± 1). Cells were cultured as described above in 20% oxygen in 10% FCS and passaged at 80% confluence from P0 to P3. Selected snoRNAs DE in the microarray were measured along with chondrogenic, and hypertrophic genes.
topographical changes in snoRnA expression. To assess topographical changes in snoRNA gene expression we utilised cartilage from high and low load areas of the equine metacarpophalangeal joint. Donor RNA (n = 5; age mean ± standard deviation 7.6 ± 0.9 years old) was extracted from cartilage of grossly normal medial and lateral condylar regions of the joint representing high and low load areas 24 . Selected snoRNAs DE in the microarray were measured along with chondrogenic, and hypertrophic genes.
Microarray data analysis. Expression Console software was utilised for array quality control. CEL files were generated with the Affymetrix GeneChip Command Console Software. The snoRNA expression data measured using Affymetrix miRNA 4.0 arrays were preprocessed using Affymetrix Expression Console with optioned method RMA 27 for normalisation. The further statistical analyses were carried out on the 1996 snoRNA and scaRNA probe sets for Homo sapiens extracted from all probes and used to determine the detected and differentially expressed (DE) snoRNAs. In each test the p value of each sample was combined using Fisher's combined p value methods. The expressions were dereplicated to transcript level by averaging replicated probes. The p value associated with the presence of dereplicated expression was assigned by combining replicated probes using Fisher's combined p test. The DE analyses on contrasting sample conditions were performed through linear models using limma package 28 in R. The significance of log fold change (logFC) values for snoRNAs were evaluated using t tests, the p values associated with logFC values were adjusted for multiple testing using the False Discovery Rate (FDR) 29 . Significantly DE was defined as those with an FDR-adjusted p value < 5%.  25 . Expression of SNORD26, SNORD96A and a panel of chondrogenic and hypertrophic markers were measured relative to cyclophilin 26 . Alkaline phosphatase (ALP) was determined following lysis of a parallel well per donor, condition and time point using an ALP activity assay using p-nitrophenyl phosphate (pNPP) as a phosphatase substrate. Data was corrected for DNA content. Prostaglandin E2 (PGE 2 ) secretion was determined in culture medium at 24 and 48 h using an ELISA (Cayman Chemicals, USA) 26 . Statistical analysis. Statistical evaluation of KL and histological scoring, and topographical gene expression results were undertaken using a Mann-Whitney Test. Other data were tested for normality prior to further analysis using Student's t test to compare between two samples, or one-way ANOVA with post hoc Tukey's test to compare between multiple samples. Statistical analysis was performed in GraphPad Prism version 8.03 (Graph-Pad, La Jolla California USA, www.graph pad.com); p values are indicated. Microarray snoRnA analysis overview. RNA was extracted and used in microarray studies from ten donors for P and U and nine for Y. Data quality assessment of the data was good and consistent for all arrays. The outcomes of variation assessment are visualised in Fig. 1A, B. Sample clustering is demonstrated by the correlation coefficient matrix heatmap (Fig. 1A). A Principal Component Analysis (PCA) plot confirmed that young samples were separated from old (P and U) samples. Samples from the P group clustered into two sub-populations; designated P1 and P2 (Fig. 1B) Fig. 2A, B). Figure 2B demonstrates the comparative features within the microarray and qRT-PCR results. Topographical area did not affect selected snoRNA expression (Supplementary File 6). Together these data show that in HAC snoRNAs are DE as a function of ageing and OA.

Selected snoRNAs expression in OA chondrocytes and following IL-1β or synovial fluid treatment of non-OA chondrocytes.
We measured seven snoRNAs DE in the microarray in isolated non-OA and OA HACs and following treatment of a pool of non-OA HACs with IL-1β or OA synovial fluid. There was a significant increase in the age of OA compared to non-OA chondrocytes (p = 0.01). We detected a significant increase in expression of SNORD116 and SNORD26 in OA chondrocytes (increased in P versus U in the micro-  Fig. 2D). There was an increase in expression of SNORD116, SNORD26, SNORD33 and SNORD98 following IL-1β treatment (Fig. 3A). Finally, we assessed the effect of OA SF on a pool of non-OA chondrocytes. There was an increase in expression of SNORD116, SNORD96A, SNORD26, SNORD33, SNORD95 and SNORD98 (Fig. 3B). Upon exposure to IL-1β or OA SF there was reduced expression of chondrogenic genes (COL2A1, ACAN, SOX9) and increase in hypertrophic genes (COL10, RUNX2, MMP13) (Supplementary File 7).

SNORD26 and SNORD96A knockdown.
To determine whether interference with the expression of a single snoRNA induces cellular changes relevant for the chondrocyte phenotype, we reduced expression of SNORD26 or SNORD96A using ASOs in HACs. A significant reduction in SNORD26 (24 h; 60%, 48 h; 61%), and SNORD96A (24 h; 72%, 48 h; 69%) was confirmed following target-specific ASOs transfection. Following both SNORD26 and SNORD96A knockdown at 24 and 48 h there were significant alterations in expression of rRNAs, chondrogenic, hypertrophic and OA-related genes (Fig. 4A-D). To functionally confirm hypertrophyrelated gene expression changes 28 at the protein level following SNORD26 and SNORD96A knockdown, we measured PGE 2 in culture supernatants (as a confirmation of COX2 expression) and ALP activity in cell extracts (as a confirmation of ALPL expression). We found an increase in COX2 expression was functionally accompanied by increased PGE 2 in supernatants and increased ALPL expression accompanied the same alterations in ALP enzymatic activity (Fig. 4E-H).
SNORD26 and SNORD96A overexpression. Reciprocally, we increased expression of SNORD26 and SNORD96A individually into primary chondrocytes by transfection of engineered snoRNA mini-genes and measured gene expression 24 h and 48 h later. We produced a profound overexpression of each snoRNA with similar trends for protein coding gene expression at each time point. However, the most pronounced effect was evident at 48 h. At 24 h there were alterations in the OA chondrocyte-relevant genes COX2 and IL6 (SNORD26) and COX2 and MMP13 (SNORD96A). At 48 h after transfection of the snoRNA mini-genes, for both SNORD26 and SNORD96A we identified significant changes in expression of chondrogenic, hypertrophic and OA-relevant genes and of 18S and 28S rRNAs for SNORD26 (Fig. 5A-D). Together with our knockdown data we demonstrated that alterations in the expression of SNORD26 and SNORD96A in primary HAC profoundly changes the chondrocyte's phenotype and rRNA expression.

Discussion
There is increasing evidence that snoRNAs have important roles in cellular development, homeostasis and disease. Indeed we previously demonstrated novel snoRNAs features in joint ageing and OA in the mouse 22 , and cartilage ageing in the horse 19 . Furthermore, we have reported that the expression of the non-canonical snoRNA RMRP is regulated during chondrogenic differentiation and regulates chondrocyte hypertrophy 29 . The snoRNA field is a terra incognita regarding knowledge on the functional implications of individual snoRNAs in cell biological processes, and this is a novel area of OA research. In OA the balance between cartilage ECM anabolism and catabolism is disrupted and important alterations in the chondrocyte phenotype are evident 2,30 . The maintenance of the ECM demands a sufficient number of functional ribosomes to translate cartilage ECM-related mRNAs into proteins and ribosome functionality may also adapt to OA chondrocyte phenotypic changes. SnoRNAs are essential for ribosome function by PTMs of rRNAs, which critically supports rRNA stability, inter-molecular interactions between rRNA and ribosomal proteins, and are essential for the ribosome's peptidyl transferase activity and mRNA decoding activity 20,28 . Here we identified for the first time that canonical snoRNA expression changes during cartilage ageing and OA and propose an additional dimension in the molecular mechanisms underlying cartilage ageing and OA. These changes in cartilage snoRNA expression appeared to be specific of/ for ageing and OA, since we demonstrated there was little effect on chondrocyte snoRNA expression of hypoxia, serum, passage number and location of cartilage within the joint (Supplementary Files 6,8,9). Our study identified that SNORD26 was increased in OA and SNORD44 and SNORD78 were reduced in ageing. Previously it has been reported that reduction of snord26, snord44 and snord78 leads to morphological abnormalities and lethality in zebra fish development (7). Also the expression of snord78 has been reported to depend on mouse development, with concomitant changes in rRNA target methylation and potential consequences for ribosome specialisation 8 . Additionally, further evidence of the biological relevance of DE snoRNAs Although most of the rRNA work in relation to PTMs and ribosomal fidelity has been undertaken in yeast, rRNA modification patterns are largely maintained throughout evolution, enabling projection from yeast to human 33 . Thus, we hypothesize that the DE of canonical snoRNAs in concert alters ribosome functionality, relevant for maintaining cartilage tissue homeostasis and a healthy chondrocyte phenotype. We were therefore intrigued by the fact that OA-relevant extracellular environments like IL-1β and OA SF altered snoRNA expression and we speculate that age-and OA-associated DE snoRNAs impact ribosomal function and thereby the translational control 34 relevant for cartilage homeostasis.
We measured the expression of selected snoRNAs from the microarray data in an independent cohort and in cultured non-OA and OA chondrocytes. From these selected snoRNAs increased expression in chondrocytes isolated from OA cartilage was evident for SNORD26, which is in concert with increased expression found in OA cartilage. The expression of SNORD96A was heavily reduced due to age in cartilage, and slightly reduced in isolated OA chondrocytes. These findings suggest that differential expression of SNORD26 is predominantly OA-related, whereas expression of SNORD96A is predominately age-related. We investigated the relevance of these two snoRNA in detail. The knockdown of SNORD26 and SNORD96A appeared to have similar gross transcriptional effects, predominately increasing OA, chondrogenic and hypertrophic gene expression and causing an overall upregulation of rRNA expression. SNORD96A knockdown is a relevant condition for the observed age-related reduction in cartilage. The apparent overall deregulation of chondrocyte gene expression following its knockdown indicates an age-related central function of SNORD96A in chondrocyte homeostasis. The overexpression of SNORD26 or SNORD96A led to more distinct snoRNA-specific differences. For both SNORD26 and SNORD96A the effect of overexpression on the mRNAs measured was most pronounced at 48 h after transfection, expected due to the time required for the β-globin gene to be transcribed, splice and 'release' the snoRNA. SNORD26 overexpression is a relevant condition for the here observed OA-related induction of cartilage SNORD26. Indeed its overexpression induced the expression of many OA-related genes, with an apparent uncoupling of expression of chondrogenic genes (increase in SOX9 and ACAN, but a reduction in COL2A1). In contrast to SNORD26, SNORD96A overexpression was less disruptive for chondrocyte transcriptional homeostasis, with mainly ALPL responding and a weak effect on COX2, IL6 and MMP13.
Non-canonical and orphan snoRNAs have no (predicted) function in the PTMs of rRNAs, have a variety of non-canonical functions. SNORD32A, SNORD33 and SNORD35A have non-canonical roles as critical SnoRNA, rRNAs, hypertrophic, chondrogenic and OA genes were determined by RT-qPCR. Data is presented as relative to the control condition (empty vector) and for statistical evaluation an independent samples t test was performed relative to the corresponding control condition using GraphPad Prism. The p values are indicated; *p < 0.05; **p < 0.01; ***p < 0.001. Data represents the mean value and error bars represent standard error mean.
Scientific RepoRtS | (2020) 10:10641 | https://doi.org/10.1038/s41598-020-67446-z www.nature.com/scientificreports/ promoters of metabolic stress. A loss of these snoRNAs causes resistance to lipotoxic and oxidative stress in vitro and prevents prorogation of oxidative stress in vivo 11 . These snoRNAs shuttle to the cytoplasm and trigger cell death in response to oxidative stress. In our study we found an increase in OA and/or ageing of SNORD33 and SNORD35A and SNORD35A. Additionally there was increased expression of SNORD33 following OA SF treatment of HACs. Oxidative stress in cartilage ageing and OA development has been reported 35 . We speculate that the OA and/or ageing-dependent increased expression of these oxidative stress-related snoRNAs is a representation of oxidative stress pathway activity. We identified SNORD116 as another DE non-canonical snoRNA in OA. SNORD116 also responded to IL-1β and OA synovial fluid treatment. In concert with the upregulation of SNORD116 in OA cartilage we previously reported the upregulation of SNORD116 in the mouse DMM model 22 . Microdeletions of the SNORD116 cluster is evident in Prader-Willi Syndrome (PWS), a disease leading to developmental delay and genetic obesity 36 . SNORD116 lacks any significant complementarity with rRNA targets and has been mainly implicated in alternative splicing of specific target genes 13 . The significance of the DE of SNORD116 in OA cartilage remains to be determined, but an inflammation-related regulation of alternative splicing in chondrocytes is expected to be an interesting avenue for further investigation 37 .
In this study we wished to determine snoRNA expression in cartilage ageing and OA. We were unable to source normal cartilage from old joints without any gross signs of OA within the entire joint. Thus, a limitation of this study is that we were unable to collect truly 'normal' old cartilage. Instead we used matched cartilage from the lesser affected lateral condyles of knees removed following TKA, as previously suggested 28 . Indeed, there was an increase in KL, Outerbridge and Modified Mankin' scoring in OA compared to old. Our data also demonstrates that snoRNA expression in chondrocytes residing in minimal OA cartilage have a distinct pattern of gene expression compared to those residing in advanced OA cartilage irrespective of gender or age.
This study comprehensively determined snoRNA signatures in ageing and OA cartilage. Given canonical but also emerging non-canonical functions of snoRNAs in homeostasis and disease, snoRNAs provide a plethora of unexplored molecular routes in chondrocyte biology. This warrants further in-depth studies addressing the cellular functions of specific snoRNAs in chondrocyte pathobiology and explore their potential as novel molecules to target OA treatment.