Age is an important factor in the development of osteoarthritis. Microarray studies provide insight into cartilage aging but do not reveal the full transcriptomic phenotype of chondrocytes such as small noncoding RNAs, pseudogenes, and microRNAs. RNA-Seq is a powerful technique for the interrogation of large numbers of transcripts including nonprotein coding RNAs. The aim of the study was to characterise molecular mechanisms associated with age-related changes in gene signatures.
Methods
RNA for gene expression analysis using RNA-Seq and real-time PCR analysis was isolated from macroscopically normal cartilage of the metacarpophalangeal joints of eight horses; four young donors (4 years old) and four old donors (>15 years old). RNA sequence libraries were prepared following ribosomal RNA depletion and sequencing was undertaken using the Illumina HiSeq 2000 platform. Differentially expressed genes were defined using Benjamini-Hochberg false discovery rate correction with a generalised linear model likelihood ratio test (P < 0.05, expression ratios ± 1.4 log2 fold-change). Ingenuity pathway analysis enabled networks, functional analyses and canonical pathways from differentially expressed genes to be determined.
Results
In total, the expression of 396 transcribed elements including mRNAs, small noncoding RNAs, pseudogenes, and a single microRNA was significantly different in old compared with young cartilage (± 1.4 log2 fold-change, P < 0.05). Of these, 93 were at higher levels in the older cartilage and 303 were at lower levels in the older cartilage. There was an over-representation of genes with reduced expression relating to extracellular matrix, degradative proteases, matrix synthetic enzymes, cytokines and growth factors in cartilage derived from older donors compared with young donors. In addition, there was a reduction in Wnt signalling in ageing cartilage.
Conclusion
There was an age-related dysregulation of matrix, anabolic and catabolic cartilage factors. This study has increased our knowledge of transcriptional networks in cartilage ageing by providing a global view of the transcriptome.
Introduction
Ageing presents huge challenges for society because whilst the lifespan increases, the quality of life faced by individuals in old age is often poor [1]. The musculoskeletal system in particular is severely affected by the ageing process, with many tissues undergoing changes that lead to loss of function and frailty. Articular cartilage is susceptible to age-related diseases, such as osteoarthritis (OA), although it is not an inevitable result of ageing but rather a consequence of a complex inter-relationship between age and further predisposing factors such as obesity [2], injury [3], genetics [4] and anatomical configuration [5].
A number of studies have interrogated ageing cartilage in order to elucidate the underlying mechanisms that contribute to OA. An age-related reduction in response to insulin-like growth factor in rats resulted in a decline in synthetic activity [6]. Furthermore, using whole mouse joints, Loeser and colleagues demonstrated that there was a reduction in extracellular matrix (ECM) gene expression in older sham-operated mice following surgical destabilisation of the medial meniscus [7]. A characteristic of ageing articular cartilage is the reduction in the number of chondrocytes within the tissue [8, 9] and there is evidence of chondrocyte senescence [10]. Chondrocyte senescence is believed to be one cause of a decline in the ability of chondrocytes to respond to growth factors; resulting in the anabolic/catabolic imbalance evident in OA [11]. One of the consequences of cell senescence is an alteration in cell phenotype [12] characterised by increased production of cytokines and growth factors. The increase in ageing chondrocytes expressing this phenotype has been proposed to contribute to cartilage ageing and, given the rise in cytokine production in OA, could directly connect ageing to OA development [13]. Furthermore, there is evidence for the role of oxidative damage in cartilage ageing from reactive oxygen species [14, 15], which can result in damage to cartilage DNA [16], whilst a link between reactive oxygen species and development of OA has also been established [17]. Hence, the outcome of ageing on chondrocyte function is an inability to maintain homeostasis when stressed.
There is a need to examine and understand the processes and mechanisms involved specifically in cartilage ageing. Whilst some insights into cartilage ageing have been learnt from transcriptome profiling studies in ageing joints using microarrays [7], these data did not identify a specific chondrocyte phenotype associated with ageing alone. Limitations in coverage and sensitivity mean that a significant part of the chondrocyte ageing transcriptomic phenotype is as yet poorly defined. Advances in high-throughput sequencing methodologies are allowing a new approach to studying transcriptomes: massively parallel sequencing of short reads derived from mRNAs known as RNA-Seq [18]. Compared with microarray technologies, RNA-Seq is demonstrated to enable more accurate quantification of gene expression levels [19]. Furthermore, RNA-Seq is an effective approach for gene expression profiling in ageing tissues with a greater dynamic range and the ability to detect noncoding RNAs [20].
Here we examine the effect of ageing on gene expression in cartilage. Using RNA-Seq analysis of RNA extracted from whole cartilage of young and old equine donors, we elucidate the differential transcriptional signatures associated with ageing and identify some of the molecular mechanisms associated with these changes.
Methods
Sample collection and preparation
Samples were collected as a byproduct 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 for this study. Full-thickness equine cartilage from the entire surface of macroscopically normal metacarpophalangeal joints of eight horses was collected from an abattoir. Horses selected were non-Thoroughbred leisure horses. No exercise history was available for the donors. Macroscopic scoring of the metacarpophalangeal joint was measured using a macroscopic grading system as described previously [21] and samples with no macroscopic perturbations were selected (combined score of zero). Subsequent RNA-Seq experiments were undertaken on normal cartilage from four young horses (4 years old) and four old horses (>15 years old).
RNA extraction
Cartilage from both articular condyles was removed from the underlying subchondral bone with a scalpel blade under sterile conditions into RNAlater (Sigma-Aldrich, Dorset, UK) according to the manufacturer's instructions. Cartilage was pulverised into a powder with a dismembranator (Mikro-S, Sartorius, Melsungen, Germany) following freezing in liquid nitrogen prior to addition of Tri Reagent (Ambion, Warrington, UK). RNA was extracted using the guanidium-thiocyanate-phenol-chloroform technique, as described previously [22]. Briefly, 20 volumes of Tri Reagent were added to the powdered cartilage tissue and incubated at room temperature for 30 minutes. Following centrifugation at 12,000×g for 10 minutes at 4°C, 200 μl chloroform was added to the supernatant, mixed and incubated at room temperature for 10 minutes. The aqueous phase was then precipitated following centrifugation at 12,000×g for 10 minutes at 4°C using 70% ethanol. RNA was purified using RNeasy spin columns (Qiagen, Crawley, UK) with on-column DNase treatment (Ambion) to remove residual gDNA according to the manufacturer's instructions. RNA was quantified using a Nanodrop ND-100 spectrophotometer (Labtech, Uckfield, UK) and assessed for purity by ultraviolet absorbance measurements at 260 nm and 280 nm.
RNA-Seq analysis: cDNA library preparation and sequencing
Eight libraries were prepared representing four animals from two groups, young (n = 4) and old (n = 4). Total RNA was analysed by the Centre for Genomic Research, University of Liverpool, for RNA-Seq library preparation and sequencing using the Illumina HiSeq 2000 platform (Illumina Inc., San Diego, CA, USA). Total RNA integrity was confirmed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Ribosomal RNA was depleted from eight total RNA samples using the Ribo-Zero™ rRNA Removal Kit (Human/Mouse/Rat; EpiCentre, Madison, WI, USA) following the manufacturer's instructions. cDNA libraries were prepared with the ScriptSeq v2 RNA-Seq library preparation kit (Epicentre) using 50 ng ribosomal-depleted RNA as the starting material and following the manufacturer's protocols. Briefly, ribosomal RNA-depleted sample was fragmented using an RNA fragmentation solution prior to cDNA synthesis. Fragment size of the final libraries and pooled libraries was confirmed using the Agilent 2100 Bioanalyzer software in the smear analysis function.
Fragmented RNA was reverse transcribed using random-sequence primers containing a tagging sequence at their 5' ends. The 3' tagging was accomplished using the Terminal-Tagging Oligo, which features a random nucleotide sequence at its 3' end, a tagging sequence at its 5' end and a 3'-blocking group on the 3'-terminal nucleotide. Terminal-Tagging Oligo randomly annealed to the cDNA, including to the 3' end of the cDNA. Purification of the di-tagged cDNA was undertaken with AMPure™ XP (Agencourt, Beckmann-Coulter, Beverly, MA, USA). The di-tagged cDNA underwent 15 cycles of amplification using polymerase chain reaction (PCR) primer pairs that annealed to the tagging sequences of the di-tagged cDNA. Excess nucleotides and PCR primers were removed from the library using AMPure™ XP (Agencourt, Beckmann-Coulter). The final pooled library was diluted to 8 pmol before hybridisation. The dilute library (120 μl) was hybridised on each of three HiSeq lanes.
Data processing
The 100-base-pair paired-end reads obtained by RNA-Seq were compiled using manufacturer-provided pipeline software (CASAVA 1.8.2; Illumina Inc., San Diego, CA, USA). Reads were then aligned onto the equine chromosomes with TOPHAT 1.3.2 (John Hopkins University, Baltimore, MD, USA) using default settings. Only uniquely mapped reads retained with less than two mismatches were used for analysis. Quality control of the reads in each lane was undertaken with FASTQC [23].
The R (version 2.15.1) Bioconductor package edgeR (version 2.13.0) [24] was used to identify differentially expressed genes. edgeR models data as a negative binomial distribution to account for biological and technical variation using a generalisation of the Poisson distribution model. Prior to assessing differential expression, data were normalised across libraries using the trimmed mean of M values normalisation method [25]. Genes were deemed differentially expressed with Benjamini-Hochberg false discovery rate-corrected P < 0.05 and fold-change ≥1.4 log2 [26] using a generalised linear model likelihood ratio test. This represents a 50% linear fold-change; that is, log21.4 = 0.5 or 50%. Statistical analysis on mapped reads was undertaken with a custom Perl script. All sequence data produced in this study have been submitted to the National Centre for Biotechnology Information GEO under Array Express [GEO:E-MTAB-1386].
Gene ontology and ingenuity pathway analysis
Owing to the minimal annotation for the equine genome, equine genes were converted to their human Ensembl orthologs prior to bioinformatics analysis. Functional analysis of age-related differentially expressed genes was undertaken to evaluate the differences in gene expression due to age. The functional analysis and clustering tool from the Database for Annotation, Visualisation, and Integrated Discovery (DAVID bioinformatics resources 6.7) was used [27].
Networks, functional analyses, and canonical pathways were generated through the use of ingenuity pathway analysis (IPA; Ingenuity Systems, Redwood City, CA, USA) on the list of differentially expressed genes with value-adjusted P < 0.05 and ± 1.4 log2 fold regulation. Gene symbols were used as identifiers and the Ingenuity Knowledge Base gene was used as a reference for pathway analysis. For network generation, a dataset containing gene identifiers and corresponding expression values was uploaded into the application. Default settings were used to identify molecules whose expression was significantly differentially regulated. These molecules were overlaid onto a global molecular network contained in the Ingenuity Knowledge Base. Networks of network-eligible molecules were then algorithmically generated based on their connectivity. The functional analysis identified the biological functions and diseases that were most significant to the dataset. A right-tailed Fisher's exact test was used to calculate P values. Canonical pathways analysis identified the pathways from the IPA library of canonical pathways that were most significant to the dataset.
Real-time polymerase chain reaction
Samples of RNA from the same pools used for the RNA-Seq analysis were used for real-time (RT)-PCR. M-MLV reverse transcriptase and random hexamer oligonucleotides were used to synthesise cDNA from 1 μg RNA (both from Promega, Southampton, UK) in a 25 μl reaction. PCR was performed on 1 μl of 10× diluted cDNA, employing a final concentration of 300 nM each primer in 20 μl reaction volumes on an ABI 7700 Sequence Detector using a SYBR Green PCR mastermix (Applied Biosystems, Paisley, Scotland, UK). Exon-spanning primer sequences were used that had been validated in previous publications [28, 29] or were designed for this study using Primer-Blast; National Centre for Biotechnology Information BLAST searches were performed for all sequences to confirm gene specificity. Oligonucleotide primers were supplied by Eurogentec (Seraing, Belgium). Steady-state transcript abundance of potential endogenous control genes was measured in the RNAseq data. Assays for four genes - glyceraldehyde-3-phosphate dehydrogenase (GAPDH), TATA box binding protein, beta-actin, and 18 ribosomal RNS - were selected as potential reference genes because their expression was unaltered. Stability of this panel of genes was assessed by applying a gene stability algorithm [30] using genormPLUS (Biogazelle, Zwijnaarde, Belgium) [31]. GAPDH was selected as the most stable endogenous control gene. Relative expression levels were normalised to GAPDH and calculated using the 2-ΔCt method [32]. Standard curves were generated from fivefold serial dilutions for each assay to confirm that all efficiencies were acceptable; within 5% of GAPDH and R2 > 0.98. Primers pairs used in this study are presented in Table 1. RT-PCR analysis data were log10 transformed to ensure normal distribution and then analysed using Student's t test.
Table 1
Gene primer sequences used in RNA-Seq validation
Gene
Accession code
Primer sequence
GAPDHa
AF157626
F: GCATCGTGGAGGGACTCA
R: GCCACATCTTCCCAGAGG
TBPa
XM_001502211
F: TGCTGCTGTAATCATGAGGGTAA
R: TCCCGTGCACACCATTTTC
ACTBa
AF035774
F: CCAGCACGATGAAGATCAAG
R: GTGGACAATGAGGCCAGAAT
18Sa
AJ311673
F: GGCGTCCCCCAACTTCTTA
R: GGGCATCACAGACCTGTTATTG
RUNX2
XM_001502519
F: TCCCTGAACTCTGCACCAAG
R: GCCAGGTAGGAGGGGTAAGA
IL7R
NM_001081942
F: GGCTATGCACAGAATGGAGACT
R: CAACTGGCTGTAGCACGAGA
SRPX
XM_001489643
F: CTGAGAACAAGGGCG-TTGC
R: CCGGAGCGTTGAGTTTGC
ACSL5
XM_001915998
F: CCTGGGCTCCTATCTCTTGC
R: CGGAGATGATCCACTCTGGC
DKK
NM_001267802
F: TAGAACCCTGGGACCTCTGG
R: GTGTCACTTTGCAAGCCTGG
ADAMTS4b
NM_001111299
F: CAGCCTGGCTCCTTCAAAAA
R: CCGCAGGAATAGTGACCACAT
COL1A1a
O46388
F: GACTGGCAACCTCAAGAAGG
R: CAATATCCAAGGGAGCCACA
COL2A1a
NM_001081764
F: TCAAGTCCCTCAACAACCAGAT
R: GTCAATCCAGTAGTCTCCGCTCTT
COL10A1a
XM_001504101
F: TGCCCAGTGGACAGGTTTCT
R: GTCTTTTCGTTTCTAGTCAGATTTTGAA
MMP1b
NM_001081847
F: GGTGAAGGAAGGTCAAGTTCTGAT
R: AGTCTTCTACTTTGGAAAAGAGCTTCTCT
MMP13b
NM_001081804
F: CTGGAGCTGGGCACCTACTG
R: ATTTGCCTGAGTCATTATGAACAAGAT
IL1βb
NM_001082526
F: GAGCCCAATCTTCAACATCTATGG
R: CAGGCTTGGTAAAAGGACTTGGTAT
TNFαb
NM_001081819
F: GCTCCAGACGGTGCTTGTG
R: GCCGATCACCCCAAAGTG
TGFβb
NM_001081849
F: CCCTGCCCCTACATTTGGA
R: CGGGTTGTGCTGGTTGTACA
ACSL5, acyl-CoA synthetase long-chain family member 5; ACTB, beta-actin; ADAMTS4, a disintegrin and metalloproteinase with thrombospondin motifs 4; COL1A1, collagen type I, alpha 1; COL2A1, collagen type II, alpha 1; COL10A1, collagen type X, alpha 1; DKK1, dickkopf homolog 1; F, forward; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; IL-1β, interleukin 1 beta; IL7R, interleukin 7 receptor; MMP1, matrix metalloproteinase 1; MMP-13, matrix metalloproteinase 13; R, receptor; RUNX2, Runt-related transcription factor 2; 18S, 18 ribosomal RNS; SRPX, Sushi repeat-containing protein; TBP, TATA box binding protein; TGFβ, transforming growth factor beta; TNFα, tumour necrosis factor alpha;. aPrimer pairs previously published in [29]. bPrimer pairs previously published in [28].
Statistical analysis
The analyses were undertaken using the software edgeR [24], S-Plus (Tibco Software Inc., Palo Alto, CA, USA), SPSS (IBM, Portsmouth, Hampshire, UK) and Excel (Microsoft, Redmond, WA, USA).
Results
Preliminary analysis of RNA-Seq data
Approximately 116 million to 235 million reads were obtained per sample. Low-quality reads were eliminated, resulting in 7 million to 58 million mapped reads (equal to 6.5 to 35% of the total reads). In total, 3 million to 49 million uniquely mapped read pairs were obtained per sample and aligned to the reference sequence of the equine genome (Equus caballus; EquCab2.56.pep [33].
Identical reads mapped to the same genomic position were retained as duplicates because these were potentially real reads. The number of genes per read were normalised to reads per kilobase of exon model per million mappable reads; the values were therefore considered the final expression level for each gene [34]. Using the E. caballus database, analysis demonstrated that in total 16,635 genes (from a total of 25,180 genes) were expressed in cartilage, which represented 66% of the equine genome. These data were used for subsequent analysis and are comparable with other recent RNA-Seq studies [35].
Age-related differential gene expression in cartilage
A multidimensional scaling plot (Figure 1A) revealed that data were clustered tightly in two groups: one for older donors, and one for younger donors.
Figure 1
Principal component analysis and volcano plot of differentially abundant transcripts. (A) Principal component analysis revealed the greatest variability in RNA-Seq data was due to the age of the donor. (B) A set of differentially expressed genes between young and old cartilage was discovered. Using the common dispersion in edgeR [24], 396 differentially expressed genes were identified with P < 0.05 (red). To enable expression of all genes to be visualised simultaneously, a smear plot was produced. The smear at the left-most edge allows visualisation of genes with zero counts in one of the groups. This was undertaken as if the total counts in one group are zero, the log fold-change is technically infinite, and the log concentration is negative infinity.
Alterations in gene expression between young and old cartilage demonstrated significant age-related changes. There were 396 genes differentially expressed with the criteria P < 0.05 and ± 1.4 log2 fold-change (Figure 1B); 93 were at higher levels in the older cartilage and 303 were at lower levels in the older cartilage. Table 2 represents the top 10 genes most differentially expressed up and down in the young horses compared with the older horses.
Table 2
Genes with the highest and lowest log2 fold-change when comparing RNA from young and old cartilage
Gene symbol
Gene name
Log2 fold-change
Qvalue
Genes with increased expression in young cartilage
CPZ
Carboxypeptidase Z
32.09
1.60 × 10-8
C18orf8
Chromosome 8 open reading frame 4
31.09
1.40 × 10-4
SRPX
Sushi repeat-containing protein SRPX
30.69
1.94 × 10-4
CYP1A1
Cytochrome P450, family 1, subfamily A, polypeptide 1
Solute carrier family 22 (extraneuronal monoamine transporter), member 3
-3.73
4.68 × 10-4
TOX3
TOX high mobility group box family member 3
-3.86
4.42 × 10-3
RELN
Reelin
-4.57
1.87 × 10-5
COCH
Coagulation factor C homolog, cochlin (Limulus polyphemus)
-4.57
1.49 × 10-4
DKK1
Dickkopf homolog 1
-4.92
6.23 × 10-4
LINGO1
Leucine-rich repeat and immunoglobulin domain containing 1
-5.09
2.01 × 10-2
SKA1
Spindle and KT associated 1
-5.55
1.55 × 10-2
RORA
RAR-related orphan receptor B
-5.98
3.26 × 10-11
Log2 fold-change and Q value (adjusted P value) were determined in edgeR [24]. A logarithm to the base 2 of 30 is approximately a linear fold-change of 4.9. Shown are the 10 genes with highest and lowest expression in old compared with young cartilage samples.
The top 25 differentially expressed genes are represented in Figure 2. The National Centre for Biotechnology Information [GEO:E-MTAB-1386] contains a complete list of all genes mapped. The subset of 93 genes that were significantly higher in older donors contained six small nuclear (SNORA)/nucleolar (SNORD) RNAs, 12 pseudogenes, 11 genes that were not identified and a single microRNA (miRNA), miR-21. Thus, 60 known protein coding genes were differentially expressed as higher in the older cartilage. Within the group where gene expression was lower in old compared with young cartilage, nine genes were SNORAs/SNORDs, one was a pseudogene and three were not identified, giving 292 known protein coding genes that were reduced in abundance in older cartilage. Table 3 presents SNORA and SNORDs that displayed age-related differential expression. Thus, 352 genes were used in downstream DAVID and IPA analysis.
Figure 2
Top 25 differentially expressed genes in cartilage ageing. The heat map illustrates the 25 most highly upregulated and downregulated genes in cartilage. The counts represent raw counts for each donor. Significance was set at P < 0.05 and ± 1.4 log2 fold-change in gene expression based on mapped reads following normalisation and statistical testing in edgeR [24]. Orange, less counts; white, greater number of counts.
Table 3
SNORDs and SNORAs identified as being differentially expressed in ageing cartilage.
Complex of snRNA and protein subunits that removes introns from a transcribed pre-mRNA
2.3
Young
U2 splicesomal RNA
Splicesome
Complex of snRNA and protein subunits that removes introns from a transcribed pre-mRNA
2.3
Young
SNORA40
H/ACA box class
H/ACA family of pseudouridylation guide snoRNAs
28SrRNA
-1.5
Old
SNORA5
H/ACA box class
H/ACA family of pseudouridylation guide snoRNAs
Not identified
-1.9
Old
SNORA30/SNORA37
H/ACA box class
H/ACA family of pseudouridylation guide snoRNAs
28SrRNA
-2.7
Old
Small nucleolar RNA U89
H/ACA box class
H/ACA family of pseudouridylation guide snoRNAs
Not identified
-1.4
Old
U4 splicesomal RNA
Splicesome
Complex of snRNA and protein subunits that removes introns from a transcribed pre-mRNA
-2.2
Old
U6 splicesomal RNA
Splicesome
Complex of snRNA and protein subunits that removes introns from a transcribed pre-mRNA
-2.1
Old
The class of action and target of these RNAs are shown with higher differential gene expression (DGE) in young or old cartilage. rRNA, ribosomal RNA; snoRNA, small nucleolar RNA.
Age-related changes in important cartilage genes
There was a reduction in the expression of 42 genes relating to the ECM, degradative proteases, matrix synthetic enzymes, cytokines and growth factors in cartilage derived from older donors compared with young donors. In comparison, there was an increase in only three ECM genes (COL10A1, COL25A1 and lubricin) together with a single growth factor (fibroblast growth factor 9) in older donors (Table 4).
Table 4
Older cartilage demonstrated reduced expression of many important cartilage genes compared with young cartilage.
Gene class
Gene name
Gene symbol
Log2 fold-change
Qvalue
DGE higher in old
ECM
Collagen, type X, alpha 1
COL10A1
-1.40
3.09 × 10-2
Collagen, type XXV, alpha 1
COL25A1
-2.78
3.62 × 10-4
Lubricin
CSPG4
-2.25
2.53 × 10-5
Growth factor/cytokine
Fibroblast growth factor 9
FGF9
-3.33
4.16 × 10-4
DGE lower in old
Protease
A disintegrin and metallopeptidase with thrombospondin 12
ADAMTS12
2.20
2.31 × 10-2
A disintegrin and metallopeptidase with thrombospondin 2
ADAMTS2
5.41
1.43 × 10-11
A disintegrin and metallopeptidase with thrombospondin 4
ADAMTS4
1.94
8.47 × 10-4
Matrix metallopeptidase 1
MMP1
1.67
6.00 × 10-3
Matrix metallopeptidase 13
MMP13
2.01
3.07 × 10-4
Plasminogen activator inhibitor-1
SERPINE1
4.24
4.04 × 10-16
Plasminogen activator, tissue
PLAT
1.86
4.19 × 10-3
Matrix enzyme
Chondroitin sulfate synthase 3
CHSY3
1.54
2.10 × 10-2
Hyaluronan synthase 3
HAS3
1.64
1.65 × 10-2
Procollagen C-endopeptidase enhancer
PCOLCE
1.90
1.89 × 10-3
ECM
Asporin
ASPN
1.55
2.72 × 10-2
Biglycan
BGN
1.47
2.25 × 10-2
Cartilage intermediate layer protein 2
CILP2
4.92
5.97 × 10-19
Chondroadherin
CHAD
2.55
7.89 × 10-6
Collagen alpha 1(V) chain
COL5A1
3.32
9.05 × 10-11
Collagen, type I, alpha 1
COL1A1
6.55
1.49 × 10-28
Collagen, type I, alpha 2
COL1A2
5.57
4.45 × 10-25
Collagen, type II, alpha 1
COL2A1
6.53
1.29 × 10-32
Collagen, type III, alpha 1
COL3A1
5.11
1.04 × 10-22
Collagen, type IV, alpha 1
COL4A1
30.11
6.16 × 10-3
Collagen, type IV, alpha 5
COL4A5
4.07
2.76 × 10-6
Collagen, type IX, alpha 1
COL9A1
8.02
9.70 × 10-32
Collagen, type IX, alpha 2
COL9A2
4.04
2.47 × 10-15
Collagen, type VIII, alpha 1
COL8A1
3.66
3.79 × 10-10
Collagen, type XI, alpha 1
COL11A1
2.01
2.79 × 10-4
Collagen, type XII, alpha 1
COL12A1
1.68
5.11 × 10-3
Collagen, type XIII, alpha 1
COL13A1
4.14
1.23 × 10-4
Collagen, type XIV, alpha 1
COL14A1
3.66
9.07 × 10-13
Collagen, type XV, alpha 1
COL15A1
2.21
7.11 × 10-5
Collagen, type XVI, alpha 1
COL16A1
2.33
1.28 × 10-5
Fibulin 1
FBLN1
4.83
4.44 × 10-9
Fibulin-7
FBLN7
1.66
1.69 × 10-2
Matrilin 2
MATN2
4.56
3.95 × 10-17
Matrilin 4
MATN4
3.69
6.65 × 10-8
Procollagen V, alpha 2
COL5A2
4.03
2.69 × 10-15
Thrombospondin 2
THBS22
2.41
1.26 × 10-4
Thrombospondin 3
THBS3
2.10
1.54 × 10-4
Growth factor/cytokine
Fibroblast growth factor 12
FGF12
2.34
1.55 × 10-2
Interleukin-11
IL11
1.48
2.46 × 10-2
Interleukin-8
IL8
1.65
1.41 × 10-2
Interleukin-1b
IL1B
6.26
8.91 × 10-10
Tumour necrosis factor, alpha-induced protein 3
TNFAIP3
1.73
1.49 × 10-2
The table illustrates significant differential gene expression (DGE) in young and old cartilage of important cartilage, extracellular matrix (ECM), cytokines and growth factors, proteases (causing cartilage degradation) and matrix enzymes (involved in matrix synthesis). Significance was set at P < 0.05 and ± 1.4 log2 fold-change in gene expression based on mapped reads following normalisation and statistical testing in edgeR [24].
Gene ontology analysis of differentially expressed genes to characterise transcriptomic signatures in cartilage ageing
DAVID analysis of all differentially expressed genes included annotations for cell adhesion and the ECM (see Additional file 1). The genes most differentially expressed, with reduced expression in cartilage from older donors, included two involved in Wnt signalling: carboxypeptidase Z and chromosome 8 open reading frame 4. Furthermore, the abundance of three other genes involved in Wnt signalling (secreted frizzled-related protein 2, Wnt11 and Wnt inhibitory factor-1) were also reduced in old cartilage. Interestingly, of the genes expressed in higher levels in older cartilage, one of the most highly regulated was the negative regulator of Wnt signalling, dickkopf homolog 1 (DKK1). DAVID analysis of this group revealed annotations for skeletal and cartilage development, and immune response.
Differential expressed genes and network analysis
Both sets of differentially expressed genes associated with ageing were analysed together in IPA with the following criteria; P < 0.05 and ± 1.4 log2 fold-change. Network-eligible molecules were overlaid onto molecular networks based on information from the ingenuity pathway knowledge database. Networks were then generated based on connectivity. (Additional file 2 contains all identified networks and their respective molecules.) Interesting age-related features were determined from the gene networks inferred. According to the top-scoring network, the differentially expressed genes were from connective tissue disorders, such as collagens COL12A1, COL16A1, COL1A1, and COL25A1 plus leucine-rich repeat and immunoglobulin domain containing 1 (LINGO), transforming growth factor beta (TGFβ)-induced 68 kDa and coclin (COCH) (Figure 3A).
Figure 3
Top-scoring networks derived from the 352 genes differentially expressed in ageing. (A) Ingenuity pathway analysis (IPA) identified connective tissue disorders as the principle associated network functions with scores of 43. (B) The second top-scoring network was a further connective tissue disorder with scores of 35. (C) IPA identified ageing significantly affected the connective tissue development and function network in ageing cartilage. Figures are graphical representations between molecules identified in our data in their respective networks. Green nodes, upregulated gene expression in older cartilage; red nodes, downregulated gene expression in older cartilage. Intensity of colour is related to higher fold-change. Key to the main features in the networks is shown.
Other networks significantly enriched also related to a further network in connective tissue disorders that contained genes including collagens COL10A1, COL11A1 and COL2A1 plus a disintegrin and metalloproteinase with thrombospondin motifs-2 (ADAMTS-2) and fibulin-1 (FBLN1) (Figure 3B). Additionally, a connective tissue development network was also significantly affected. The genes most affected in this network included acyl-synthetase long chain family member 5 (ACSL5), phosphate-regulating neutral endopeptidase (PHEX) and DKK1 (Figure 3C).
Significant IPA canonical pathways are demonstrated in Table 5 and the associated molecules of the top canonical pathways identified are in Additional file 3. These include atherosclerosis signalling, prothrombin activation and rheumatoid arthritis.
Table 5
IPA canonical pathways were significantly affected in ageing cartilage
Name of canonical pathway
Pvalue
Ratio
Atherosclerosis signalling
3.80 × 10-9
15/136 (0.11)
Role of osteoblasts, osteoclasts and chondrocytes in rheumatoid arthritis
3.41 × 10-6
16/238 (0.067)
Intrinsic prothrombin activation
9.82 × 10-6
6/35 (0.171)
Hepatic fibrosis and stellate cell activation
9.92 × 10-6
12/146 (0.082)
Role of macrophages, fibroblasts and endothelial cells in rheumatoid arthritis
1.73 × 10-4
16/333 (0.048)
The significance of the association between the dataset and the canonical pathway was measured using a ratio of the number of molecules from the dataset that mapped to the pathway divided by the total number of molecules that map to the canonical pathway is displayed. Fisher's exact test was used to calculate P values.
Confirmation of differential gene expression using real-time PCR measurements of selected genes
To validate the RNA-Seq technology, 14 genes were selected to measure using reverse transcription and RT-PCR based on differences noted in the arrays and/or their potential importance in the OA process. This was performed on the original RNA from all donors used to perform the RNA-Seq experiment (Table 6). Genes were selected based on differences noted in the RNA-Seq results. All genes were found to have comparable results with RNA-Seq data; for instance, genes identified as having an increase in expression in older samples in the RNA-Seq experiment also gave increased expression relative to GAPDH following RT-PCR. Statistical significance was tested using Student's t test. Two genes whose expressions were not significantly altered in RNA-Seq results - tumour necrosis factor alpha and transforming growth factor β (TGFβ) - were also unaltered when assessed with RT-PCR.
Table 6
Real-time polymerase chain reaction analysis of 14 selected genes reveals good correlation with RNA-Seq results
Gene
RNA-Seq results
RT-PCR results
Differential expression
Significant log2 fold-change
Qvalue
Age
2-ΔCt log2 fold-change
Pvalue
Mean Ct value, young
Mean Ct value, old
Young
Old
DKK1
Higher in old
-4.9
0.0006
0.0044 ± 0.006
0.0338 ± 0.20
-3.0
0.024
30.9
28.8
COL10
Higher in old
-1.4
0.03
0.0006 ± 0.0001
0.0013 ± 0.0009
-1.0
0.26
32.4
32.3
RUNX2
Lower in old
2.2
0.038
0.0093 ± 0.009
0.0034 ± 0.002
1.5
0.002
29
30.8
SRPX
Lower in old
1.6
0.009
0.0046 ± 0.005
0.0006 ± 0.0005
3.0
0.007
30
33.8
ACSL5
Lower in old
30.1
0.02
0.0155 ± 0.0002
0.0062 ± 0.005
1.3
0.09
29.5
29.6
IL7R
Lower in old
30.1
0.06
0.0005 ± 0.0002
0.0001 ± 0.0001
2.3
0.06
33.9
37
COL2A1
Lower in old
6.5
2.29 × 10-32
63.3246 ± 54.7
1.3165 ± 1.12
5.7
0.04
18.3
22.4
COL1A1
Lower in old
6.5
1.48 × 10-28
3.4815 ± 1.57
0.0278 ± 0.02
7.0
0.15
24.3
27.5
MMP1
Lower in old
1.6
0.006
0.7093 ± 0.21
0.4027 ± 0.14
0.8
0.05
23.1
24
MMP13
Lower in old
2
0.0003
0.2323 ± 0.15
0.044 ± 0.04
2.4
0.1
26.1
26.8
ADAMTS-4
Lower in old
1.9
0.0008
0.5121 ± 0.35
0.1345 ± 0.05
1.9
0.07
24
24.9
IL1β
Lower in old
6.2
8.91 × 10-10
0.0057 ± 0.005
0.0004 ± 0.0003
3.9
0.05
30.1
33.4
TNFα
No change
Not significant
1
0.0041 ± 0.004
0.001 ± 0.0003
2.1
0.28
31.1
32.1
TGFβ
No change
Not significant
1
1.2865 ± 0.23
2.3124 ± 1.24
-0.8
0.15
21.9
21.4
Values for real-time polymerase chain reaction (RT-PCR) are the mean ± standard deviation of relative expression levels normalised to expression of glyceraldehyde-3-phosphate dehydrogenase. Statistical significance was tested using Student's t test. Log2 fold-change of 2-ΔCt values is shown for comparison. Average threshold cycle (Ct) values for young and old donors are demonstrated to indicate robustness of expression. DKK1, dickkopf homolog 1; COL10, collagen type X; RUNX2, Runt-related transcription factor 2; SRPX, Sushi repeat-containing protein; ACSL5, acyl-CoA synthetase long-chain family member 5; IL7R, interleukin 7 receptor; COL2A1, collagen type II, alpha 1; COL1A1, collagen type I, alpha 1; MMP1, matrix metalloproteinase 1; MMP-13, matrix metalloproteinase 13; ADAMTS4, a disintegrin and metalloproteinase with thrombospondin motifs 4; IL-1β, interleukin 1 beta; TNFα, tumour necrosis factor alpha; TGFβ, transforming growth factor beta.
In addition, quantitative RT-PCR was undertaken for the 14 genes on a different set of donors to those used in the RNASeq study in order to validate our findings young (4 years old, n = 4) and old (>15 years old, n = 4) (Table 7). All genes were found to have comparable results.
Table 7
Real-time PCR analysis of 14 selected genes using a different set of donors reveals similar correlation with RNA-Seq results
Gene name
Differential expression
Age
Pvalue
Young
Old
DKK1
Higher in old
0.0004 ± 0.0006
0.008 ± 0.0006
0.04
COL10
Higher in old
7.66 × 10-5 ± 3.05 × 10-5
0.000148 ± 7.13 × 10-5
0.1
RUNX2
Lower in old
0.002 ± 0.0004
0.0005 ± 0.0004
0.02
SRPX
Lower in old
0.0025 ± 0.0002
0.001 ± 0.002
0.05
ACSL5
Lower in old
0.004 ± 0.002
0.002 ± 0.02
0.13
IL7R
Lower in old
0.001 ± 0.0002
0.0007 ± 0.0004
0.07
COL2A1
Lower in old
32.6 ± 24.5
0.2 ± 0.13
0.04
COL1A1
Lower in old
0.71 ± 0.11
0.02 ± 0.02
0.0003
MMP1
Lower in old
0.28 ± 0.17
0.07 ± 0.03
0.03
MMP13
Lower in old
0.08 ± 0.12
0.02 ± 0.02
0.28
ADAMTS-4
Lower in old
0.07 ± 0.05
0.03 ± 0.02
0.07
IL1β
Lower in old
0.0001 ± 0.0002
2.66 × 10-5 ± 1.14 × 10-5
0.03
TNFα
No change
0.0001 ± 7.7 × 10-5
9.8 × 10-5 ± 4.5 × 10-5
0.4
TBFβ
No change
0.58 ± 0.14
1.06 ± 0.5
0.12
Values for real-time polymerase chain reaction (RT-PCR) are the mean ± standard deviation of relative expression levels normalised to expression of glyceraldehyde-3-phosphate dehydrogenase. Statistical significance was tested using Student's t test. DKK1, dickkopf homolog 1; COL10, collagen type X; RUNX2, Runt-related transcription factor 2; SRPX, Sushi repeat-containing protein; ACSL5, acyl-CoA synthetase long-chain family member 5; IL7R, interleukin 7 receptor; COL2A1, collagen type II, alpha 1; COL1A1, collagen type I, alpha 1; MMP1, matrix metalloproteinase 1; MMP-13, matrix metalloproteinase 13; ADAMTS4, a disintegrin and metalloproteinase with thrombospondin motifs 4; IL-1β, interleukin 1 beta; TNFα, tumour necrosis factor alpha; TGFβ, transforming growth factor beta.
Discussion
Ageing has an important role in the development of OA by making the joint more susceptible to OA risk factors. To provide interventions to prevent age-related changes and reduce the risk of developing OA, the underlying mechanisms involved in age-related changes of cartilage require elucidation. Characterisation of both young and old cartilage at the molecular level is essential for identifying the important signalling pathways in OA development. In the present study, we used the RNA-Seq technique to undertake deep transcriptome profiling of young and old cartilage. This is the first time that, to our knowledge, this technique has been used to interrogate transcriptional changes in cartilage ageing and, importantly, validation studies using RT-PCR demonstrated high correlation between methodologies and demonstrated reproducibility using a different donor set.
This study built on previous findings that identified a reduction in matrix gene expression with joint ageing [7]. We took a single tissue, articular cartilage, and undertook RNA-Seq in order to interrogate a greater range of genes for differential expression. Not surprisingly, our experiments identified that the age of the donor accounted for the principal variability in the data. The major findings of this study were as follows: the age-related gene expression changes identified were most notably involving reduced differential gene expression in older cartilage; there was an over-representation of genes with reduced expression relating to the ECM, degradative proteases, matrix synthetic enzymes, cytokines and growth factors in cartilage derived from older donors compared with young donors; cartilage ageing caused a decrease in many important Wnt signalling genes; IPA revealed that the top-scoring network for differentially expressed genes was from connective tissue disorders and connective tissue development; IPA also demonstrated significant canonical pathways for atherosclerosis signalling, prothrombin activation and rheumatoid arthritis; and there was differential expression of pseudogenes and small noncoding RNAs in cartilage ageing with increased expression of 12 pseudogenes and six noncoding RNAs in older cartilage, and of one pseudogene and nine small noncoding RNAs in younger cartilage.
Equine tissue was readily obtained, enabling collection of cartilage samples from macroscopically normal, skeletally mature young and aged horses. Importantly, the horse suffers clinical joint diseases similar to man (reviewed [36]), and as such has been used as a model for naturally occurring OA [37] due to extensive knowledge of its pathogenesis and clinical experience of the disease [38]. Indeed, the incidence of equine metacarpophalangeal OA in young racehorses [39] in training is similar to the incidence of post-traumatic OA in man [40]. Additionally, the articular cartilage thickness is also comparable between species [41].
For young horses one year is equivalent to about 3.5 years of a human [42, 43]. The rate of equine ageing relative to equivalent human age is greatest within the first two years of life and decreases after the horse reaches maturity at 4 years of age [44]. Hence, horses >15 years old, as used in this study, are likely to equate to humans older than 52 years. The average lifespan of a horse is 25 to 30 years and so it is possible that the obvious differences in lifespan may yield significant differences in the effect of ageing amongst animal species due to cumulative lifetime load. However, whilst the work in this study may not be directly applied to humans, it does enable an insight into human cartilage ageing by studying a population at skeletal maturity to one beyond the middle age equivalent in man.
This study utilised the entire articular surface of distal metacarpal III bone. High and low load-bearing cartilage was thus used. An assessment of macroscopic changes revealed no abnormalities in our samples. Previous studies indicated a high correlation between gross scoring and Mankin's grading in equine cartilage from the distal metacarpal III bone [28, 45]. To validate that the RNA extracted from the harvested tissue was articular cartilage, the expression level of several genes typically expressed and those of bone were measured. There was a high expression of articular cartilage genes only (data not shown).
Previous studies have identified a number of age-related changes in chondrocyte metabolism (reviewed in [46]). Most of these studies demonstrate changes at the protein level, such as an age-related decline in matrix production when equine chondrocytes were stimulated with TGFβ1 [47]. Others have provided evidence for a chondrocyte senescence secretory phenotype in ageing, demonstrated by an increase in cytokines [48, 49] along with matrix metalloproteinase (MMP) production and a reduction in growth factors [50, 51]. These studies did not interrogate transcript changes and of course simple deduction of protein from mRNA expression is insufficient because post-translational regulation, small noncoding RNAs, decay differences in mRNA and proteins, and locations or molecular associations of proteins affect overall protein levels [52]. However, a recent whole mouse-joint study demonstrated a reduction in matrix genes with age [7] in agreement with our findings. Furthermore, a study of equine articular cartilage concluded that although there was no change in the age-related expression of MMP-13 there was a reduction in MMP-3 and interleukin (IL)-1β gene expression in cartilage from older donors [53]. Annotations of genes at reduced levels in older samples included many relating to the ECM, degradative proteases, matrix synthetic enzymes, cytokines and growth factors. In contrast, within these annotations those at higher levels in older cartilage were very small: COLX, COLXXV, lubricin and fibroblast growth factor 9.
There appears to be an age-related failure of matrix, anabolic and catabolic cartilage factors. This is interesting because a recent study on postnatal and skeletally mature equine cartilage identified a reduction in collagens, matrix modelling and noncollagenous matrix transcripts with age [54]. ADAMTS-4 expression was reduced in the older cartilage in this study, which is in agreement with findings in ageing rat cartilage [55]. In contrast, previous studies have demonstrated an increase in IL-7 in ageing chondrocytes and in response to fibronectin fragments or IL-1 [49]. Although our experiment did not identify IL-7, interestingly one of the most downregulated genes identified in this study was the IL-7 receptor. A reduction in IL-7 receptor signalling in ageing β-progenitor cells has been demonstrated previously to result in ageing-like gene expression profiles [56]. Also, whereas other studies have demonstrated an increase in IL-1 [48] (where an increase in IL-1 protein was evident in older cultured human chondrocytes) and MMP-13 [48, 57] in ageing human cartilage, this study identified an age-related decline in their transcript abundance. However, one MMP-13 study looked at catabolic responsiveness with age whilst another used immunolocalisation of MMP-13 to identify protein. These two factors are not always related [58]. Whilst differences could also be attributed to our age classification of young and old and species distinctions, increased matrix enzymes (MMP-1, MMP-13) and cytokines such as IL-1, IL-8 and IL-11 identified in younger cartilage could be due to increased turnover. Interestingly a recent study identified that low innate capacity to produce IL-1β and IL-6 was associated with the absence of OA in old age [59]. The reduction in IL-1β evident in older cartilage may represent a protective mechanism against OA.
We noted in cartilage derived from old donors that there was primarily a reduction in the expression of some key Wnt signalling genes plus an increase in the Wnt antagonist DKK1 and a reduction in RUNX2, a downstream target of Wnt. Wnt signalling is active in adult cartilage, with deregulation being detrimental, resulting in age-associated joint pathologies due to excessive remodelling and degradation [60]. This signalling pathway has also been found to both regulate matrix synthesis in chondrocyte cell lines [61] and stimulate catabolic genes such as MMP-13 and ADAMTS-4 in chondrocytes [62]. A recent study demonstrated a potential protective function of Wnt in ageing. The activation of the Wnt pathway inhibited IL-1-mediated MMP-13 expression in human chondrocytes through the direct interaction between nuclear factor-κB and β-catenin [63]. One study has linked Wnt signalling with chondrocyte hypertrophy through RUNX2 activation [64], whilst elsewhere it was shown that DKK1 is a major player in the cessation of hypertrophic differentiation that can contribute to OA [65]. Interestingly, COL10A1, a marker of chondrocyte hypertrophy, was increased in old cartilage. However, COL10A1 has also been identified in the transitional zone of cartilage and may have a role in the modification of collagen fibril arrangement [66]. A recent study in mesenchymal stems cells derived from OA patients found that COL10A1 downregulation played a role in the establishment of a defective cartilage matrix in OA [67]. It would seem that this increased expression with ageing is not through the Wnt signalling interaction with subsequent RunX2 activation as described previously [64]. Further credence is given to this hypothesis by our findings that alkaline phosphatase expression, also regulated through RunX2, was downregulated in old cartilage. Overall Wnt signalling is involved in maintenance of cartilage, and the dysregulation event here in ageing may be an important episode. Interfering with the pathway may contribute to improvements in cartilage regeneration.
Using IPA, this study identified age-related changes in pathways and processes including connective tissue disorders and development in which a significant number of genes, regulated both strongly and subtly, were enriched. This is not remarkable given the number of matrix genes differentially identified in the study. Care should also be taken in overinterpretation of this finding because some of the genes in this network are minor components of cartilage, such as COL12A, COL16A, COL25A, LINGO and COCH. Canonical pathways identified as significantly affected by ageing, such as the role of osteoblasts and osteoclasts in rheumatoid arthritis, were not surprising. Interestingly, age-affected atherosclerosis signalling pathways follow the differential expression of a mixture of proteases and lipoproteins. In ageing cartilage, further studies to investigate the significance of this are clearly required.
One advantage for the use of RNA-Seq to undertake differential gene expression studies is that other sets of RNA molecules from the transcriptome can be identified, such as nonprotein coding RNAs (for example, miRNA and small nucleolar RNA (snoRNA)) that constitute a significant part of the transcriptome [68] as well as pseudogenes.
Pseudogenes provide a novel tier of gene regulation through the generation of endogenous silencing RNA or miRNA binding sites, which act as decoys for miRNAs [69]. Indeed some miRNAs have been demonstrated to target the genes [70]. It is hypothesised that pseudogenes act as post-transcriptional regulators of the corresponding parental gene [71]. Whilst possessing very similar sequences to their counterpart coding genes, they are unable to be transcribed due to mutation/deletion or insertion of nucleotides. Transcription of pseudogenes has tissue specificity and can be activated or reduced in disease, indicating a possible functional role in cells [72]. Interestingly, pseudogenes have been identified as increasing with age, such as pseudogene cyclin D2 in the ovary [73]. Whilst this study identified the differential expression of pseudogenes in cartilage of different ages, it is not known whether these are functional or have relevance to cartilage ageing. Recent work by the Encyclopaedia of DNA Elements (ENCODE) Consortium identified that 8% of the pseudogenes in the human genome are functional [74], and so with the publication of GENCODE, a reference human genome annotation for The ENCODE Project [75], more light may be shed relating to the role of pseudogenes in cartilage ageing in the near future. Pseudogenes thus present an interesting area for future research in cartilage ageing and disease.
The methodology used here does not enrich for miRNAs. To increase the identifications of small miRNAs using RNA-Seq, specific techniques are used for their enrichment in conjunction with additional miRNA abundance quantification algorithms. A single miRNA, miR-21, was however identified as increased in ageing cartilage. miRNAs are short noncoding RNAs that regulate the translation [76] and/or degradation of target message [77]. miR-21 has been implicated in inflammation [78], cancers including osteosarcomas [79], and hypomethylation [80]. The role of miR-21 in cartilage is not fully elucidated, although a study in rats found that it promoted increased proliferation and matrix synthesis in chondrocytes embedded in atelocollagen gel [81]. However, our finding is interesting because epigenetic changes such as hypomethylation occur with ageing, a risk factor contributing to several age-related pathologies [82].
A further set of small noncoding RNAs, snoRNAs - a class of small guide RNAs found in the nucleolus - were also identified in the study. The snoRNAs direct chemical modification of other RNAs, and like miRNAs are emerging as important regulators of cellular function and disease development. There are two principle classes: the C/D box snoRNAs (SNORDs) and H/ACA box snoRNAs (SNORAs), which are associated with methylation and pseudouridylation of ribosomal and other RNAs. In addition, RNase MRP and RNaseP are the only members of a further special class of snoRNAs [83]. Both were significantly reduced in older cartilage in this study. Interestingly, mutations in RNase MRP cause cartilage hair hypoplasia in which patients display dwarfism [84]. In recent work, RNase MRP was identified as a regulator of chondrocyte hypertrophy, demonstrating functional cross talk with chondrogenic pathways [85]. snoRNAs fine-tune the ribosome to accommodate changing requirements for protein production during development, normal function and disease [86]. Indeed, control of snoRNA expression may play a pivotal role in the regulation of high protein-producing cells such as chondrocytes, as demonstrated by the phenotypes of ribosomopathies [87]. Whilst there are very few studies into the significance of snoRNAS in cartilage ageing or disease, a recent study proposed the use of serum snoRNA U38 and U48 as biomarkers of early cartilage damage. These snoRNAs was detected in serum following anterior cruciate ligament injury, but were not associated with normal ageing [88]. The snoRNA transcriptome signatures in ageing cartilage provide an interesting set of genes for further studies to determine their role in ageing.
Conclusions
A major strength of this study is that it represents the first application of RNA-Seq technology for transcriptomic studies in cartilage ageing. The study has increased our knowledge of transcriptional networks by providing a global view of the transcriptome. The molecular signatures described in this paper reflect a combination of degenerative processes and transcriptional responses to the process of ageing. This analysis further supports the use of next-generation sequencing as an ideal quantitative framework to study pathways and networks as an integrated system in order to understand the complex processes of cartilage ageing.
Abbreviations
ADAM-TS:
a disintegrin and metalloproteinase with thrombospondin
COCH:
Coclin
DAVID:
Database for Annotation, Visualisation, and Integrated Discovery
DKK1:
Dickkopf homolog
ECM:
extracellular matrix
GAPDH:
glyceraldehyde-3-phosphate dehydrogenase
IL:
interleukin
IPA:
ingenuity pathway analysis
LINGO:
leucine-rich repeat and immunoglobulin domain containing
miRNA:
microRNA: MMP: matrix metalloproteinase
OA:
osteoarthritis
PCR:
polymerase chain reaction
RNA-Seq:
RNA-sequencing
RT:
real time
snoRNA:
small nucleolar RNA
TGFβ:
transforming growth factor beta
Wnt:
Wingless and Int signalling pathway.
Declarations
Acknowledgements
This study was supported by a Wellcome Trust Veterinary Integrated Research Fellowship (WT088557MA). The authors would like to thank Tim Welting, Maastricht University for his guidance with the snoRNA data.
Authors’ Affiliations
(1)
Comparative Musculoskeletal Biology, Institute of Ageing and Chronic Disease, University of Liverpool
(2)
Centre for Genomic Research, Institute of Integrative Biology, University of Liverpool
References
Beard JR, Biggs S, Bloom BD, Fried LP, Hogan L, Kalache A, Olshansky SJ: Global Population Ageing: Peril or Promise. 2011, Geneva: Forum WE (World Economic Forum)Google Scholar
Cooper C, Snow S, McAlindon TE, Kellingray S, Stuart B, Coggon D, Dieppe PA: Risk factors for the incidence and progression of radiographic knee osteoarthritis. Arthritis Rheum. 2000, 43: 995-1000. 10.1002/1529-0131(200005)43:5<995::AID-ANR6>3.0.CO;2-1.View ArticlePubMedGoogle Scholar
Samilson RL, Prieto V: Dislocation arthropathy of the shoulder. J Bone Joint Surg Am. 1983, 65: 456-460.PubMedGoogle Scholar
Ma WJ, Guo X, Liu JT, Liu RY, Hu JW, Sun AG, Yu YX, Lammi MJ: Proteomic changes in articular cartilage of human endemic osteoarthritis in China. Proteomics. 2011, 11: 2881-2890. 10.1002/pmic.201000636.View ArticlePubMedGoogle Scholar
Martin JA, Ellerbroek SM, Buckwalter JA: Age-related decline in chondrocyte response to insulin-like growth factor-I: the role of growth factor binding proteins. J Orthop Res. 1997, 15: 491-498. 10.1002/jor.1100150403.View ArticlePubMedGoogle Scholar
Loeser RF, Olex AL, McNulty MA, Carlson CS, Callahan MF, Ferguson CM, Chou J, Leng X, Fetrow JS: Microarray analysis reveals age-related differences in gene expression during the development of osteoarthritis in mice. Arthritis Rheum. 2012, 64: 705-717. 10.1002/art.33388.PubMed CentralView ArticlePubMedGoogle Scholar
Aigner T, Kim HA, Roach HI: Apoptosis in osteoarthritis. Rheum Dis Clin North Am. 2004, 30: 639-653. 10.1016/j.rdc.2004.04.002. xiView ArticlePubMedGoogle Scholar
Adams CS, Horton WE: Chondrocyte apoptosis increases with age in the articular cartilage of adult animals. Anat Rec. 1998, 250: 418-425. 10.1002/(SICI)1097-0185(199804)250:4<418::AID-AR4>3.0.CO;2-T.View ArticlePubMedGoogle Scholar
Martin JA, Buckwalter JA: Telomere erosion and senescence in human articular cartilage chondrocytes. J Gerontol A Biol Sci Med Sci. 2001, 56: B172-B179. 10.1093/gerona/56.4.B172.View ArticlePubMedGoogle Scholar
Mueller MB, Tuan RS: Anabolic/catabolic balance in pathogenesis of osteoarthritis: identifying molecular targets. PM R. 2011, 3 (6 Suppl 1): S3-S11.View ArticlePubMedGoogle Scholar
Campisi J: Senescent cells, tumor suppression, and organismal aging: good citizens, bad neighbors. Cell. 2005, 120: 513-522. 10.1016/j.cell.2005.02.003.View ArticlePubMedGoogle Scholar
Loeser RF: Age-related changes in the musculoskeletal system and the development of osteoarthritis. Clin Geriatr Med. 2010, 26: 371-386. 10.1016/j.cger.2010.03.002.PubMed CentralView ArticlePubMedGoogle Scholar
Jallali N, Ridha H, Thrasivoulou C, Underwood C, Butler PE, Cowen T: Vulnerability to ROS-induced cell death in ageing articular cartilage: the role of antioxidant enzyme activity. Osteoarthritis Cartilage. 2005, 13: 614-622. 10.1016/j.joca.2005.02.011.View ArticlePubMedGoogle Scholar
Loeser RF, Carlson CS, Del Carlo M, Cole A: Detection of nitrotyrosine in aging and osteoarthritic cartilage: correlation of oxidative damage with the presence of interleukin-1beta and with chondrocyte resistance to insulin-like growth factor 1. Arthritis Rheum. 2002, 46: 2349-2357. 10.1002/art.10496.View ArticlePubMedGoogle Scholar
Chen AF, Davies CM, De Lin M, Fermor B: Oxidative DNA damage in osteoarthritic porcine articular cartilage. J Cell Physiol. 2008, 217: 828-833. 10.1002/jcp.21562.PubMed CentralView ArticlePubMedGoogle Scholar
Kurz B, Jost B, Schunke M: Dietary vitamins and selenium diminish the development of mechanically induced osteoarthritis and increase the expression of antioxidative enzymes in the knee joint of STR/1N mice. Osteoarthritis Cartilage. 2002, 10: 119-126. 10.1053/joca.2001.0489.View ArticlePubMedGoogle Scholar
Matkovich SJ, Zhang Y, Van Booven DJ, Dorn GW: Deep mRNA sequencing for in vivo functional analysis of cardiac transcriptional regulators: application to Galphaq. Circ Res. 2010, 106: 1459-1467. 10.1161/CIRCRESAHA.110.217513.PubMed CentralView ArticlePubMedGoogle Scholar
de Magalhaes JP, Finch CE, Janssens G: Next-generation sequencing in aging research: emerging applications, problems, pitfalls and possible solutions. Ageing Res Rev. 2010, 9: 315-323. 10.1016/j.arr.2009.10.006.PubMed CentralView ArticlePubMedGoogle Scholar
Kawcak CE, Frisbie DD, Werpy NM, Park RD, McIlwraith CW: Effects of exercise vs experimental osteoarthritis on imaging outcomes. Osteoarthritis Cartilage. 2008, 16: 1519-1525. 10.1016/j.joca.2008.04.015.View ArticlePubMedGoogle Scholar
Robinson MD, McCarthy DJ, Smyth GK: edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010, 26: 139-140. 10.1093/bioinformatics/btp616.PubMed CentralView ArticlePubMedGoogle Scholar
Fu YR, Yi ZJ, Guan SZ, Zhang SY, Li M: Proteomic analysis of sputum in patients with active pulmonary tuberculosis. Clin Microbiol Infect. 2012, 18: 1241-1247.View ArticlePubMedGoogle Scholar
Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. Methodology. 1995, 57: 289-300.Google Scholar
Huang da W, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009, 4: 44-57.View ArticlePubMedGoogle Scholar
Barr ED: The association of bone and cartilage in matrix proteolysis of articular cartilage, and its role in palmar/plantar osteochondral dosease in the racing thoroughbred. PhD thesis. 2010, University of LiverpoolGoogle Scholar
Taylor SE, Vaughan-Thomas A, Clements DN, Pinchbeck G, Macrory LC, Smith RK, Clegg PD: Gene expression markers of tendon fibroblasts in normal and diseased tissue compared to monolayer and three dimensional culture systems. BMC Musculoskelet Disord. 2009, 10: 27-10.1186/1471-2474-10-27.PubMed CentralView ArticlePubMedGoogle Scholar
Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F: Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002, 3: RESEARCH0034-PubMed CentralView ArticlePubMedGoogle Scholar
Hellemans J, Mortier G, De Paepe A, Speleman F, Vandesompele J: qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 2007, 8: R19-10.1186/gb-2007-8-2-r19.PubMed CentralView ArticlePubMedGoogle Scholar
Livak KJ, Schmittgen TD: Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001, 25: 402-408. 10.1006/meth.2001.1262.View ArticlePubMedGoogle Scholar
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008, 5: 621-628. 10.1038/nmeth.1226.View ArticlePubMedGoogle Scholar
O'Loughlin A, Lynn DJ, McGee M, Doyle S, McCabe M, Earley B: Transcriptomic analysis of the stress response to weaning at housing in bovine leukocytes using RNA-seq technology. BMC Genomics. 2012, 13: 250-10.1186/1471-2164-13-250.PubMed CentralView ArticlePubMedGoogle Scholar
Innes JF, Clegg P: Comparative rheumatology: what can be learnt from naturally occurring musculoskeletal disorders in domestic animals?. Rheumatology (Oxford). 2010, 49: 1030-1039. 10.1093/rheumatology/kep465.View ArticleGoogle Scholar
McIlwraith CW, Frisbie DD, Kawcak CE: The horse as a model of naturally occurring osteoarthritis. Bone Joint Res. 2012, 1: 297-309. 10.1302/2046-3758.111.2000132.PubMed CentralView ArticlePubMedGoogle Scholar
Goodrich LR, Nixon AJ: Medical treatment of osteoarthritis in the horse - a review. Vet J. 2006, 171: 51-69. 10.1016/j.tvjl.2004.07.008.View ArticlePubMedGoogle Scholar
Neundorf RH, Lowerison MB, Cruz AM, Thomason JJ, McEwen BJ, Hurtig MB: Determination of the prevalence and severity of metacarpophalangeal joint osteoarthritis in Thoroughbred racehorses via quantitative macroscopic evaluation. Am J Vet Res. 2010, 71: 1284-1293. 10.2460/ajvr.71.11.1284.View ArticlePubMedGoogle Scholar
Dirschl DR, Marsh JL, Buckwalter JA, Gelberman R, Olson SA, Brown TD, Llinias A: Articular fractures. J Am Acad Orthop Surg. 2004, 12: 416-423.PubMedGoogle Scholar
Brommer H, Laasanen MS, Brama PA, van Weeren PR, Helminen HJ, Jurvelin JS: Functional consequences of cartilage degeneration in the equine metacarpophalangeal joint: quantitative assessment of cartilage stiffness. Equine Vet J. 2005, 37: 462-467.View ArticlePubMedGoogle Scholar
Mauderly JL, Hahn FF: The effects of age on lung function and structure of adult animals. Adv Vet Sci Comp Med. 1982, 26: 35-77.PubMedGoogle Scholar
Paradis MR: Demographics of health and disease in the geriatric horse. Vet Clin North Am Equine Pract. 2002, 18: 391-401. 10.1016/S0749-0739(02)00021-4.View ArticlePubMedGoogle Scholar
Peffers MJ: Protein and transcriptomic signatures of cartilage ageing and disease. PhD thesis. 2013, University of LiverpoolGoogle Scholar
Loeser RF: Aging and osteoarthritis: the role of chondrocyte senescence and aging changes in the cartilage matrix. Osteoarthritis Cartilage. 2009, 17: 971-979. 10.1016/j.joca.2009.03.002.PubMed CentralView ArticlePubMedGoogle Scholar
Iqbal J, Dudhia J, Bird JL, Bayliss MT: Age-related effects of TGF-beta on proteoglycan synthesis in equine articular cartilage. Biochem Biophys Res Commun. 2000, 274: 467-471. 10.1006/bbrc.2000.3167.View ArticlePubMedGoogle Scholar
Forsyth CB, Cole A, Murphy G, Bienias JL, Im HJ, Loeser RF: Increased matrix metalloproteinase-13 production with aging by human articular chondrocytes in response to catabolic stimuli. J Gerontol A Biol Sci Med Sci. 2005, 60: 1118-1124. 10.1093/gerona/60.9.1118.PubMed CentralView ArticlePubMedGoogle Scholar
Long D, Blake S, Song XY, Lark M, Loeser RF: Human articular chondrocytes produce IL-7 and respond to IL-7 with increased production of matrix metalloproteinase-13. Arthritis Res Ther. 2008, 10: R23-10.1186/ar2376.PubMed CentralView ArticlePubMedGoogle Scholar
Blaney Davidson EN, Scharstuhl A, Vitters EL, van der Kraan PM, van den Berg WB: Reduced transforming growth factor-beta signaling in cartilage of old mice: role in impaired repair capacity. Arthritis Res Ther. 2005, 7: R1338-R1347. 10.1186/ar1833.PubMed CentralView ArticlePubMedGoogle Scholar
Baldi P, Long AD: A Bayesian framework for the analysis of microarray expression data: regularized t-test and statistical inferences of gene changes. Bioinformatics. 2001, 17: 509-519. 10.1093/bioinformatics/17.6.509.View ArticlePubMedGoogle Scholar
Trumble TN, Trotter GW, Oxford JR, McIlwraith CW, Cammarata S, Goodnight JL, Billinghurst RC, Frisbie DD: Synovial fluid gelatinase concentrations and matrix metalloproteinase and cytokine expression in naturally occurring joint disease in horses. Am J Vet Res. 2001, 62: 1467-1477. 10.2460/ajvr.2001.62.1467.View ArticlePubMedGoogle Scholar
Mitani H, Takahashi I, Onodera K, Bae JW, Sato T, Takahashi N, Sasano Y, Igarashi K: Comparison of age-dependent expression of aggrecan and ADAMTSs in mandibular condylar cartilage, tibial growth plate, and articular cartilage in rats. Histochem Cell Biol. 2006, 126: 371-380. 10.1007/s00418-006-0171-8.View ArticlePubMedGoogle Scholar
Curtis H, Marusyk A, Zaberezhnyy V, Casas M, Adane B, Merz A, Serkova N, DeGregori J: Aging-associated alterations in IL-7 receptor signaling and inflammation promote declining B-lymphopoiesis and increased leukemogenesis. J Immunol. 2012, 188: 109-Google Scholar
Billinghurst RC, Dahlberg L, Ionescu M, Reiner A, Bourne R, Rorabeck C, Mitchell P, Hambor J, Diekmann O, Tschesche H, Van Wart H, Poole A: Enhanced cleavage of type II collagen by collagenases in osteoarthritic articular cartilage. J Clin Invest. 1997, 99: 1534-1545. 10.1172/JCI119316.PubMed CentralView ArticlePubMedGoogle Scholar
Greenbaum D, Colangelo C, Williams K, Gerstein M: Comparing protein abundance and mRNA expression levels on a genomic scale. Genome Biol. 2003, 4: 117-10.1186/gb-2003-4-9-117.PubMed CentralView ArticlePubMedGoogle Scholar
Goekoop RJ, Kloppenburg M, Kroon HM, Frolich M, Huizinga TW, Westendorp RG, Gussekloo J: Low innate production of interleukin-1beta and interleukin-6 is associated with the absence of osteoarthritis in old age. Osteoarthritis Cartilage. 2010, 18: 942-947. 10.1016/j.joca.2010.03.016.View ArticlePubMedGoogle Scholar
Yates KE, Shortkroff S, Reish RG: Wnt influence on chondrocyte differentiation and cartilage function. DNA Cell Biol. 2005, 24: 446-457. 10.1089/dna.2005.24.446.View ArticlePubMedGoogle Scholar
Zhu M, Chen M, Zuscik M, Wu Q, Wang YJ, Rosier RN, O'Keefe RJ, Chen D: Inhibition of beta-catenin signaling in articular chondrocytes results in articular cartilage destruction. Arthritis Rheum. 2008, 58: 2053-2064. 10.1002/art.23614.PubMed CentralView ArticlePubMedGoogle Scholar
Yuasa T, Otani T, Koike T, Iwamoto M, Enomoto-Iwamoto M: Wnt/beta-catenin signaling stimulates matrix catabolic genes and activity in articular chondrocytes: its possible role in joint degeneration. Lab Invest. 2008, 88: 264-274. 10.1038/labinvest.3700747.View ArticlePubMedGoogle Scholar
Ma B, van Blitterswijk CA, Karperien M: A Wnt/beta-catenin negative feedback loop inhibits interleukin-1-induced matrix metalloproteinase expression in human articular chondrocytes. Arthritis Rheum. 2012, 64: 2589-2600. 10.1002/art.34425.View ArticlePubMedGoogle Scholar
Dong YF, Soung do Y, Schwarz EM, O'Keefe RJ, Drissi H: Wnt induction of chondrocyte hypertrophy through the Runx2 transcription factor. J Cell Physiol. 2006, 208: 77-86. 10.1002/jcp.20656.View ArticlePubMedGoogle Scholar
Buckland J: Osteoarthritis: control of human cartilage hypertrophic differentiation. Nat Rev Rheumatol. 2012, 8: 368-View ArticlePubMedGoogle Scholar
Lammi PE, Lammi MJ, Hyttinen MM, Panula H, Kiviranta I, Helminen HJ: Site-specific immunostaining for type X collagen in noncalcified articular cartilage of canine stifle knee joint. Bone. 2002, 31: 690-696. 10.1016/S8756-3282(02)00904-3.View ArticlePubMedGoogle Scholar
Lamas JR, Rodriguez-Rodriguez L, Vigo AG, Alvarez-Lafuente R, Lopez-Romero P, Marco F, Camafeita E, Dopazo A, Callejas S, Villafuertes E, Hoyas J, Tornero-Esteban M, Urcelay E, Fernandez-Gutierreq B: Large-scale gene expression in bone marrow mesenchymal stem cells: a putative role for COL10A1 in osteoarthritis. Ann Rheum Dis. 2010, 69: 1880-1885. 10.1136/ard.2009.122564.View ArticlePubMedGoogle Scholar
Kapranov P, Cheng J, Dike S, Nix DA, Duttagupta R, Willingham AT, Stadler PF, Hertel J, Hackermuller J, Hofacker IL, Bell I, Cheung E, Drenkow J, Dumais E, Patel S, Helt G, Ganes M, Ghosh S, Piccolboni A, Sementchenko V, Tammana H, Gingeras T: RNA maps reveal new RNA classes and a possible function for pervasive transcription. Science. 2007, 316: 1484-1488. 10.1126/science.1138341.View ArticlePubMedGoogle Scholar
Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP: A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language?. Cell. 2011, 146: 353-358. 10.1016/j.cell.2011.07.014.PubMed CentralView ArticlePubMedGoogle Scholar
Poliseno L, Salmena L, Zhang J, Carver B, Haveman WJ, Pandolfi PP: A coding-independent function of gene and pseudogene mRNAs regulates tumour biology. Nature. 2010, 465: 1033-1038. 10.1038/nature09144.PubMed CentralView ArticlePubMedGoogle Scholar
Muro EM, Mah N, Andrade-Navarro MA: Functional evidence of post-transcriptional regulation by pseudogenes. Biochimie. 2011, 93: 1916-1921. 10.1016/j.biochi.2011.07.024.View ArticlePubMedGoogle Scholar
Zheng D, Frankish A, Baertsch R, Kapranov P, Reymond A, Choo SW, Lu Y, Denoeud F, Antonarakis SE, Snyder M, Ruan Y, Wei CL, Gingeras TR, Guigo R, Harrow J, Gerstein MB: Pseudogenes in the ENCODE regions: consensus annotation, analysis of transcription, and evolution. Genome Res. 2007, 17: 839-851. 10.1101/gr.5586307.PubMed CentralView ArticlePubMedGoogle Scholar
Choi D, Yoon S, Lee E, Hwang S, Yoon B, Lee J: The expression of pseudogene cyclin D2 mRNA in the human ovary may be a novel marker for decreased ovarian function associated with the aging process. J Assist Reprod Genet. 2001, 18: 110-113. 10.1023/A:1026538826392.PubMed CentralView ArticlePubMedGoogle Scholar
Pei B, Sisu C, Frankish A, Howald C, Habegger L, Mu XJ, Harte R, Balasubramanian S, Tanzer A, Diekhans M, Reymond A, T J, Harrow J, Gerstein MB: The GENCODE pseudogene resource. Genome Biol. 2012, 1: R51-View ArticleGoogle Scholar
Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, Aken BL, Barrell D, Zadissa A, Searle S, Barnes I, Bignell A, Boychenko V, Hunt T, Kay M, Mukherjee G, Rajan J, Despacio-Reyes G, Saunders G, Steward C, Harte R, Lin M, Howald C, Tanzer A, Derrien T, Chrast J, Walters N, Balasubramanian S, Pei B, Makitie O, et al: GENCODE: the reference human genome annotation for The ENCODE Project. Genome Res. 2012, 22: 1760-1774. 10.1101/gr.135350.111.PubMed CentralView ArticlePubMedGoogle Scholar
Berninger P, Gaidatzis D, van Nimwegen E, Zavolan M: Computational analysis of small RNA cloning data. Methods. 2008, 44: 13-21. 10.1016/j.ymeth.2007.10.002.View ArticlePubMedGoogle Scholar
Ziyan W, Shuhua Y, Xiufang W, Xiaoyun L: MicroRNA-21 is involved in osteosarcoma cell invasion and migration. Med Oncol. 2011, 28: 1469-1474. 10.1007/s12032-010-9563-7.View ArticlePubMedGoogle Scholar
Pan W, Zhu S, Yuan M, Cui H, Wang L, Luo X, Li J, Zhou H, Tang Y, Shen N: MicroRNA-21 and microRNA-148a contribute to DNA hypomethylation in lupus CD4+ T cells by directly and indirectly targeting DNA methyltransferase 1. J Immunol. 2010, 184: 6773-6781. 10.4049/jimmunol.0904060.View ArticlePubMedGoogle Scholar
Kongcharoensombat W, Nakasa T, Ishikawa M, Nakamae A, Deie M, Adachi N, Mohamed A, Ochi M: The effect of microRNA-21 on proliferation and matrix synthesis of chondrocytes embedded in atelocollagen gel. Knee Surg Sports Traumatol Arthrosc. 2010, 18: 1679-1684. 10.1007/s00167-010-1111-7.View ArticlePubMedGoogle Scholar
Fraga MF, Esteller M: Epigenetics and aging: the targets and the marks. Trends Genet. 2007, 23: 413-418. 10.1016/j.tig.2007.05.008.View ArticlePubMedGoogle Scholar
Ridanpaa M, van Eenennaam H, Pelin K, Chadwick R, Johnson C, Yuan B, vanVenrooij W, Pruijn G, Salmela R, Rockas S, Makitie O, Kaitila I, de la Chapelle A: Mutations in the RNA component of RNase MRP cause a pleiotropic human disease, cartilage-hair hypoplasia. Cell. 2001, 104: 195-203. 10.1016/S0092-8674(01)00205-7.View ArticlePubMedGoogle Scholar
Caron MM, Steinbusch M, Reicherter K, Mattijssen S, Surtel DA, van Rhijn LW, Pruijn GJ, Lausch E, Zabel B, Welting TJ: RNase MRP is a novel regulator of endochondral ossification. Osteoarthritis Cartilage. 2013, 21 (Suppl):Google Scholar
Zhang L, Yang M, Marks P, White LM, Hurtig M, Mi QS, Divine G, Gibson G: Serum non-coding RNAs as biomarkers for osteoarthritis progression after ACL injury. Osteoarthritis Cartilage. 2012, 20: 1631-1637. 10.1016/j.joca.2012.08.016.PubMed CentralView ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
