Decoding the Regulatory Landscape of Ageing in Musculoskeletal Engineered Tissues Using Genome-Wide DNA Methylation and RNASeq

Mesenchymal stem cells (MSC) are capable of multipotent differentiation into connective tissues and as such are an attractive source for autologous cell-based regenerative medicine and tissue engineering. Epigenetic mechanisms, like DNA methylation, contribute to the changes in gene expression in ageing. However there was a lack of sufficient knowledge of the role that differential methylation plays during chondrogenic, osteogenic and tenogenic differentiation from ageing MSCs. This study undertook genome level determination of the effects of DNA methylation on expression in engineered tissues from chronologically aged MSCs. We compiled unique DNA methylation signatures from chondrogenic, osteogenic, and tenogenic engineered tissues derived from young; n = 4 (21.8 years ± 2.4 SD) and old; n = 4 (65.5 years±8.3SD) human MSCs donors using the Illumina HumanMethylation 450 Beadchip arrays and compared these to gene expression by RNA sequencing. Unique and common signatures of global DNA methylation were identified. There were 201, 67 and 32 chondrogenic, osteogenic and tenogenic age-related DE protein-coding genes respectively. Findings inferred the nature of the transcript networks was predominantly for ‘cell death and survival’, ‘cell morphology’, and ‘cell growth and proliferation’. Further studies are required to validate if this gene expression effect translates to cell events. Alternative splicing (AS) was dysregulated in ageing with 119, 21 and 9 differential splicing events identified in chondrogenic, osteogenic and tenogenic respectively, and enrichment in genes associated principally with metabolic processes. Gene ontology analysis of differentially methylated loci indicated age-related enrichment for all engineered tissue types in ‘skeletal system morphogenesis’, ‘regulation of cell proliferation’ and ‘regulation of transcription’ suggesting that dynamic epigenetic modifications may occur in genes associated with shared and distinct pathways dependent upon engineered tissue type. An altered phenotype in engineered tissues was observed with ageing at numerous levels. These changes represent novel insights into the ageing process, with implications for stem cell therapies in older patients. In addition we have identified a number of tissue-dependant pathways, which warrant further studies.


Introduction
The limited ability of articular cartilage, bone and tendon to regenerate has prompted the development of cell-based tissue engineering techniques. One cell therapy option is mesenchymal stem cells (MSC); a heterogeneous population of multi-potent cells with the ability to differentiation into tissues including cartilage, bone and tendon, thus accommodating tissue repair and homeostasis. The principles of tissue engineering involve a multifarious interaction of factors, and knowledge of the extent MSC phenotype and differentiation capacity alter with ageing is limited. Subsequently, any loss in functionality with age would have profound consequences for the maintenance of tissue viability and the quality of tissues. MSCs have been utilised in clinical trials of cell therapies for cartilage repair and osteoarthritis (reviewed [1]), bone fracture treatment [2] and in a limited number of tendon therapies [3]. However, the therapeutic efficiency of MSCs for clinical applications remains limited, possibly due to the attenuation of their regenerative potential in aged patients with chronic diseases.
Advancing age is a prominent risk factor that is closely linked with the onset and progression of diseases such as osteoarthritis, osteoporosis and tendinopathy. Understanding the influence that ageing has on chondrogenic, osteogenic and tenogenic progenitor cells such as MSCs is important in determining how these processes affect their capacity to differentiate into functional chondrocytes, osteoblasts and tenocytes for use in therapeutic applications. A model using MSCs derived from young and old donors to musculoskeletal engineered tissues could aid in our understanding of musculoskeletal ageing.
To understand the underlying mechanisms that are responsible for age-related changes in musculoskeletal engineered tissues, a number of studies have been undertaken on ageing MSCs (reviewed [4]), as well as the differentiation potential of tissue engineered cartilage [5] and bone [6], though no studies have addressed these questions in tendon.
There are a few studies investigating the effect of chronological age of donor MSC on engineered tissues, some with conflicting findings. One study found a reduction in glycosaminoglycans in chondrogenesis with age [7] whereas another experiment using a wider donor agerange found no change [8]. Contrasting results of the chondrogenic differentiation potential of adult MSCs has been described with one study reporting age independence [9], whilst another demonstrated a negative correlation with advancing age in male but not female donors [10]. A study demonstrated that foetal and adult MSCs are differentially regulated by transforming growth factor-β stimuli to activate the onset of chondrogenesis, suggesting that discrete agerelated mechanisms direct chondrogenic regulation following development and postnatal maturation [11].
Age-related changes in osteogenesis have been more widely studied. Osteogenic progenitors in MSCs derived from rat bone marrow demonstrated an age-related decline. In addition MSCs from young rats had a significantly greater bone formation capability in vivo compared with aged rats [12]. The osteogenic potential of MSCs is independent of advancing age in adult human donors [13]. However, a decline in the osteogenic precursor population, due to accelerated senescence and lower rate of population doublings in MSCs isolated from older donors suggests a reduction in osteoblast formation. This may contribute to the age-related reduction in bone formation in the elderly [14]. In age-related studies of osteogenic differentiation one group identified an increase in alkaline phosphatase with age [15] whilst another demonstrated a decrease [16]. It is thought these discrepancies could be due to the heterogeneous population which is propagated within and amongst donor populations.
Few studies have investigated the effects of ageing MSCs on tendon tissue engineering. However, one study on human tendon stem cells from aged tendons described reduced proliferation capacity and premature entry into senescence [17]. A recent study in rat tendonderived stem cells from older donors demonstrated earlier entry into senescence which was postulated to be due to a reduction in the levels of miRNA-135a, a ROCK-1 targeting micro-RNA (miR) that blocks entry into senescence pathways. This may be due to a reduction in miR-135a, which binds to ROCK-1 and inhibits entry into senescence in young tendon. Thus a decrease in miR-135a in older tendon may be the cause of reduced stem cell proliferation, selfrenewal and tenogenic differentiation [18].
The advent of global DNA methylation arrays and RNASeq studies have made it possible to explore gene methylation and/or expression during cell development [19], tissue differentiation [20], disease [21] and ageing [22,23]. In addition, the global relationship between gene methylation and expression can now be investigated in ageing [24]. Whilst global methylation and RNASeq are powerful tools to study methylation variation and transcription changes, no joint analysis with these two types of data have been reported in tissue engineering. Tissue engineering aims to develop biomimetic tissues that recapitulate biological, structural and functional characteristics of native tissue. Thus age-related changes have potential implications for the tissue engineering strategies used for enhancing musculoskeletal repair. In this study we evaluated and compared the methylome and transcriptome of chondrogenic, osteogenic and tenogenic engineered tissues derived from young and old human bone marrow derived MSCs in order to determine similar and distinct changes with ageing. In doing so we have identified areas for future research to improve functionality of ageing MSC derived engineered tissues.

Materials and Methods
All chemicals are supplied by Sigma unless stated otherwise.

Cell Culture and Differentiation
Human MSCs from young; n = 4 (21.8years±2.4SD) and old; n = 4 (65.5years±8.3SD) donors (Stem Cell Technologies, Grenoble, France and Promocell, Heidelberg, Germany), grown to passage 4 and each donor each differentiated into chondrogenic [25], osteogenic [26] and tenogenic [27] tissues as previously described and used in all subsequent experiments [28]. All tissue culture was undertaken in 5% oxygen and tissues harvested at 21 days (osteogenic) and 28 days (chondrogenic and tenogenic). All cells were purchased and thus ethical approval was not required.

Validation of differentiation
Differentiation of chondrogenic and osteogenic engineered tissues was assessed by comparing to MSCs treated identically except with maintenance media (complete Dulbecco's Modified Eagles Media (Gibco)) using histology, and quantitative real-time PCR (qRT-PCR). Calcium depositions were determined in osteogenic tissues using Alizarin red staining as previously described [29]. Chondrogenic pellets were paraffin embedded and 4 μm sections taken and further stained with Alcian Blue/Nuclear Fast Red. Tendon engineered tissues were fixed in 4% paraformaldehyde, longitudinally embedded in paraffin and 4μm sectioned on polylysine slides. Staining was undertaken with Masson's Trichrome [30].
Transmission electron microscopy (TEM) of tendon tissues were performed by fixation in 2.5% glutaraldehyde in 0.1M sodium cacodylate buffer followed for 8 hours, followed by buffer washing procedure and second fixation and contrast stain with 0.1% osmium tetroxide for 90 minutes. Samples were stained with 8% uranyl acetate in 0.69% maleic acid for 90 minutes, dehydrated in ascending ethanol concentrations and embedded in epoxy resin. 60-90 nm sections were cut with a Reichert-Jung Ultracut on an ultramicrotome using a diamond knife, mounted on 200 mesh copper grids and stained with 'Reynold's Lead citrate' stain for 4 minutes. Images were viewed in Philips EM208S Transmission Electron Microscope at 80k.
RNA was extracted from all assays and converted to cDNA to analyse lineage-specific gene expression markers using qRT-PCR relative to GAPDH [22]. All primer sequences are in S1 file.

RNA Data processing
The RNASeq data was processed as previously described [22,24]. Concise details are in S2 file. Data was assessed using pairwise comparisons, correlation heatmaps and PCA plots and outliers removed accordingly. For RNASeq and smallSeq differentially expressed genes (DEGs) and transcripts were extracted by applying the threshold of false discovery rate (FDR) adjusted pvalues < 0.05, generated using the Benjamini  Genomic DNA isolation, bisulphite treatment and methylation profiling Genomic DNA was extracted using the SureSelect gDNA Extraction Kit (Agilent, Santa Clara, USA) according to manufacturer's instructions, 500 ng of genomic DNA was then bisulphite converted using the EZ-96 DNA Methylation Kit (Zymo Research, Irvine, USA). DNA methylation profiling of the samples was carried out by Cambridge Genomic Services (Cambridge, UK), using the Illumina Infinium HumanMethylation450 Beadchip array (Illumina, Inc., San Diego, USA).

Methylation data processing
GenomeStudio (Illumina Inc., San Diego,USA) was used to extract the raw data. GenomeStudio provides the methylation data as β values: β = M/(M + U) (M represents the fluorescent signal of the methylation probe; U represents the signal of the un-methylated probe). β values range from 0 (no methylation) to 1 (100% methylation). The raw methylation data was processed using R (version 3.0.1) and the Watermelon package (version 2.12) as has been previously described [24,33]. Probes with a detection P value > 0.01 were removed. Age-related differential methylation was defined as Benjamini-Hochberg corrected P value [32] < 0.01 or <0.05 (differentially methylated loci (DML) and gene/CpG island/promoter respectively) and a mean methylation difference (Δ β score) 0.15 (15%), as previously reported [34].

Functional analysis of transcriptomic and methylation data
To determine gene ontology, functional analyses, networks, canonical pathways and proteinprotein interactions of age-related differentially expressed genes and methylated genes we performed analyses using Panther Classification System [35] and the functional analysis and clustering tool from the Database for Annotation, Visualisation, and Integrated Discovery (DAVID bioinformatics resources 6.7) [36] (using expressed genes as a reference), Ingenuity Pathway Analysis (IPA) [37]. Targetscan v6.2 was used to identify potential miR [38].
Relative gene expression using real-time polymerase chain reaction (qRT-PCR) qRT-PCR was undertaken on engineered tissues from similarly sourced MSCs at P4 from an independent cohort from those used for the RNA-Seq analysis young; n = 3 (22.2years±2.3SD) and old; n = 3 (64.8years±6.6SD). Primers were either validated in previous publications [23,39] and supplied by Eurogentec (Seraing, Belgium) or designed and validated commercially (Primer Design, Southampton, UK). Steady-state transcript abundance of potential endogenous control genes was measured in the RNA-Seq data. Assays for four genes-Glutaraldehyde dehydrogenase (GAPDH), ribosomal protein 13 (RPS8), ribosomal protein 13 (RPS13), and ribosomal protein 16 (RPS16) were selected as potential reference genes as their expression was unaltered. Stability of this panel of genes was assessed by applying a gene stability algorithm [40]. RPS8 was selected as the most stable endogenous control gene. For miRs cDNA synthesis was performed using 100 ng RNA and miRscript RT kit II (Qiagen, Crawley, UK) according to the manufacturer's protocol. qRT-PCR analysis was performed using miRScript SybrGreen Mastermix (Qiagen, Crawley, UK) using Rnu-6 as the endogenous control. Relative expression levels were calculated by using the 2−ΔCt method [40]. Data were analysed statistically using GraphPad Prism 6 (GraphPad Software, San Diego, CA, USA) following normality testing using a Mann-Whitney test at a 0.05 level of significance.

Characterisation of engineered tissues
To evaluate chondrogenesis markers of chondrocytes were assessed; Alcian Blue staining for glycosaminoglycans and aggrecan gene expression. There was an increase in Alcian Blue staining and aggrecan (ACAN) expression (Fig 1A and 1B). Osteogenesis was evaluated with Alizarin Red staining. There was a significant increase in staining with Alizarin Red both visually and using quantitative analysis (Fig 1C and 1D) demonstrating osteogenesis. Tenogenic differentiation was evaluated histologically using Masson's Trichrome staining indicating areas of organised and disorganised collagen fibril formation within the tissues which were confirmed with TEM and gene expression of COL1A1 (Fig 1E and 1F), Serpin peptidase inhibitor F (SER-PINF1) ( Fig 1G) and thrombospondin 4 (THBS4) (Fig 1H) [41]. There were no age-related differences in the differentiation markers measured (data not shown).

Overview of RNASeq and smallSeq data
For the RNASeq data an average of 68.53 million pairs of 100-bp paired-end reads per sample was generated that aligned to the human reference sequence. Based on data variation assessment one young and one old sample from chondrogenic were classed as outliers and removed from subsequent data analysis. Mapping results are summarised in the S3 file. Of the 63,152 human genes, between 39.9% and 47% had at least one read aligned [42]. This is similar to the output of other RNA-Seq studies [43] In the smallSeq data an average of 12.2 million 50bp single-end reads was generated. This represented an average of 44% of reads mapped. Many of the 4206 human small RNAs were mapped with at least one read; 21.5-38.5% within all samples. Mapping results are summarised in S4 file and are similar to other small RNASeq studies [44,45]. Reads were used to estimate . Gene expression of aggrecan following chondrogenic differentiation, young and old donors combined. Data are represented as 2^-DCT compared with GAPDH. Box and whisker plots represent the median and 25-75 percentiles. Statistical evaluation was undertaken using Mann Whitney-U test (n = 6). c; Osteogenic differentiation from MSCs was confirmed with Alizarin Red S staining at day 21 to visualize mineralized bone matrix following extraction of the calcified mineral from the stained monolayer at low pH (young donor shown). d; Box and whisker plot showing quantitative results of Alizarin red staining of all donors, statistical significance Mann-Whitney-U test p<0.001 (n = 12). e; Histology images of a tendon engineered tissue (young donor shown) made from MSCs stained with Masson's Trichrome to identify collagenous matrix. Image was captured at x4 magnification and x10 magnification inset (upper image); scale bar is 100μm. Example of more organised areas of collagen is marked on the inset image. Lower image depicts ultrastructural analysis using scanning transmission electron microscopy. The presence of aligned extracellular collagen fibrils (A) and less organised collagen (B) are inset in red; scale bar is 1μm. Tenogenic differentiation was also evaluated with using gene expression of f; COL1A1, g; SERPINF1 and h; THBS4. Data from all donors are represented as 2^-DCT compared with GAPDH. Statistical evaluation was undertaken using Mann Whitney-U test (n = 8). small RNA transcript expression of all samples using FPKM in order to identify the most abundant miRs and small nucleolar RNAs (snoRNAs). S5 file demonstrates the expression levels of the entire data set and highlights the top 10 highly expressed small RNAs genes within each class.
Identification of differentially expressed genes and differentially spliced genes using RNASeq For RNASeq a principal component analysis (PCA) plot (Fig 2A) of log 2 gene expression data identified age-related biological variation within all engineered tissue groups. Hierarchical clustering using a sample to sample distance matrix identified clustering principally by engineered tissue type (Fig 2B).
Sets of age-related differentially expressed (DE) genes were identified including proteincoding RNA, long non-coding RNA (lnc), small nucleolar RNA (snoRNA), small nuclear RNA (snRNA), pseudogenes and miRs ( In total 94190±3005 chondrogenic, 116105 ±3008 osteogenic, and 113075±5346 tenogenic isoforms were detected (mean±standard deviation)). No isoforms were differentially expressed. However, using Cuffdiff to calculate changes in the relative splice abundances by quantifying the square root of the Jensen-Shannon divergence on primary transcripts with at least two isoforms, identified 119, 21 and 9 differential splicing events in chondrogenic, osteogenic and tenogenic tissues respectively (alternative splicing (AS)) (S9 file). These included small nucleolar RNAs, long non-coding RNAs and pseudogenes.
For the smallSeq PCA of log 2 gene expression data indicated the age-effect was weak ( Fig  2D and 2E). The greatest variability was due to engineered tissue type. There were no agerelated DE small RNAs in chondrogenic tissues. In osteogenic tissues, there were 36 DE miRs (all reduced were derived from old MSCs) and three DE snoRNAs and in tenogenic engineered tissues three miRs were DE ( Fig 2F, Table 3). The donor age-associated DE of several miRs in the osteogenic and tenogenic tissues was validated using qPCR (Table 4). Validated miRs were chosen based on our own and published data with regards to the relevance to the osteogenicand tenogenic-related processes.
Reproducibility of RNASeq results with an independent platform is high [22,23]. Nevertheless we selected genes (mRNA and miRNA) DE and assessed their expression levels with qRT-PCR analyses for each engineered tissue type. There was good correlation between the deep sequencing analyses and qRT-PCR results (Table 5) reflecting the accuracy and reliability of deep sequencing analyses.

Gene ontology (GO) and IPA analysis of DEGs and AS genes
For each engineered tissue type age-related DEGs (adjusted P<0.05 and 1.4 log 2 fold change) were analysed in DAVID. Significant annotations included shared terms 'glycoprotein' and 'extracellular matrix' (ECM) for chondrogenic and osteogenic. In addition the terms 'growth factor' and 'secreted' were identified for chondrogenic and osteogenic respectively. For tenogenic 'developmental protein' and 'homeobox' were significantly enriched. The DEGs were next input into IPA. This inferred the nature of the engineered tissue protein-coding transcript networks was predominantly for the functions 'cell death and survival', 'cell morphology', and 'cell growth and proliferation'. The canonical pathways significantly identified from the gene sets are shown in Table 6. The top networks identified for each engineered tissue type were 'developmental disorders, hereditary disorders and metabolic disease' for chondrogenic ( Fig  3A); 'cellular growth and proliferation, cell development' and 'morphology' for osteogenic ( Fig  3B); and 'embryonic and organismal development' for tenogenic ( Fig 3C).
GO analyses using PANTHER indicated enrichment in genes associated principally with metabolic processes in all engineered tissue type genes undergoing AS ( Fig 4A). The chondrogenic and tenogenic AS gene sets was then analysed with IPA. For chondrogenic tissues the top

Gene pairing analysis of DE miRs and DE RNAs
The expression patterns of DE miRs and mRNA were further analysed using IPA by investigating opposite fold-change direction (up/down or down/up), following the canonical miR-mRNA target expression paradigm with moderate to high confidence. Potentially relevant miR-mRNA signatures involved in the age-related changes were identified; 16 for osteogenic and one for tenogenic (S10 file). Using PANTHER the mRNA in which related miRs were identified in osteogenic tissues were enriched in the cellular components ECM (52% of genes) and enriched for the functions 'binding (44%).

Comparison of the DNA methylome in ageing MSCs
Unsupervised hierarchical clustering revealed that young and old samples are distinguished by their DNA methylome in all engineered tissue types ( Fig 5). Technical triplicate replicates were included for a single donor for chondrogenic and osteogenic young donors and correlation was excellent. Significant age-related differentially methylated loci (DML), both tissue specific and common CpGs, were identified in all engineered tissue types (Table 7). S11 file contains all DML for each engineered tissue type and S12 file identifies these at the site, promoter, gene and CpG island level. In all engineered tissue groups and regions hypomethylation in old samples was dominant.
Gene ontology analysis of genes containing DML indicated age-related enrichment for all engineered tissue types in skeletal system morphogenesis, regulation of cell proliferation and regulation of transcription.  To identify the canonical pathways, biological function, and networks that were affected by the differentially methylated genes, we used IPA analysis. Results suggest that dynamic epigenetic modifications may occur in genes associated with a number of shared and distinct pathways dependent upon engineered tissue type. The top 10 genes with increased and decreased methylation levels based on Beta values (methylation difference) are listed in Table 8.
Canonical pathways were analysed based on the ratio of input genes to the total number of reference genes in the corresponding pathways in the IPA knowledge bases. Fisher's exact test was employed to calculate the P values to determine significant associations between the DM genes and the canonical pathways. The top five canonical pathways for each engineered tissue type are in Table 9. Then we used IPA comparison analysis to visualise downstream effects analysis results across each engineered tissue type simultaneously. This identified diseases and biological functions predicted to increase or decrease related to age-affected DML through functional scores (Fig 6A). Interestingly in all engineered tissue types the function 'differentiation of cells' was activated (Fig 6B) whereas the 'cell survival' network was only affected in chondrogenic and osteogenic tissues. The network 'congenital anomalies of the musculoskeletal system' was activated in tenogenic but inhibited in chondrogenic and osteogenic. The most significant network for each engineered tissue was 'skeletal and muscular development and function'.
Next we wanted to determine to what degree the age-related gene expression differences among engineered tissue types are affected by epigenetic changes. The methylation of gene promoters and/or enhancers is known to correlate with decreased gene expression, contrastingly methylation within non-enhancer regions of the gene body correlates with increased gene expression [46]. Therefore for chondrogenic, osteogenic and tenogenic tissues we compared DEGs from RNASeq with location of DMLs and found 10, 13 and 4 genes identified in both data sets for comparison (Table 10).

Discussion
Adult MSCs are an appealing source for cell-based treatment for musculoskeletal diseases and injury [47]. Our previous work in bone-marrow-derived MSCs using a systems biology approach demonstrated an altered phenotype in MSC ageing at a number of levels, implicating roles for inflammageing and mitochondrial ageing [48]. The changes identified represented novel insights into the ageing process, with implications for stem cell therapies in older patients. Tissue engineering aims to develop biomimetic tissues that recapitulate biological, structural and functional characteristics of native tissue. However there is sparse global information available on the effect of donor age on engineered musculoskeletal tissues at the transcriptome and methylome level. Age-related changes have potential implications for the tissue-engineering strategies used for enhancing musculoskeletal repair. Furthermore the study of musculoskeletal ageing in bone, cartilage and tendon are usually undertaken in seclusion and it is frequently difficult to attain aged matched tissue samples in humans. Therefore we propose our approach as a potential model of musculoskeletal ageing that could be probed further to identify factors that may aid in recapitulation of a younger tissue phenotype.
It is known that the site of MSC extraction can affect cell behaviour [49] we therefore used MSCs derived from alveolar bone. Additionally, as low oxygen tension improves MSC vitality and metabolic state in culture [50] all MSCs and then subsequently tissues were cultured in 5% oxygen tension. Standard methods of engineered tissue characterisation were undertaken following chondrogenic and osteogenic differentiation.
Transcriptome profiling is a key method for functional characterization of cells and tissues. However one challenge of this study was the integration of the different types of data. We used gene ontology, network analysis and annotated of the pathways identified. However, due to limited information on the donors we were unable to integrate key environmental factors such as nutrition and other lifestyle features which can alter molecular measurements. Other potential limitations of the study include small sample size and significance threshold filtering which can affect the subsequent pathway/network analysis.
Epigenetic processes have been implicated in age-related musculoskeletal diseases such as osteoarthritis [21] and osteoporosis (reviewed [51]). This study identified a number of epigenetic molecular classes including small non-coding RNAs; small nucleolar RNAs, small cajal body RNA (scaRNAs), miRs and lncRNAs.
Our study found weak age-related effects on expression at the miR level with no DE small RNAs in chondrogenic engineered tissues. The low mapping in chondrogenic samples implies that RNA populations or fragments other than the targeted small RNA were the input material into the library prep workflows. Further investigation did indeed demonstrate that the low percentage of alignments to the small RNA reference dataset corresponds to a high mapping percentage to rRNA in chondrogenic samples. However the percentage of mapping to rRNA for old and young samples was similar resulting in a lower sequencing depth which may reduce the statistical power in differential expression analysis. This effect was roughly similar for the two sample groups. Therefore, the count values are usable values, though the statistical power may be weaker due to lower sequencing depth. Compared to tenogenic and osteogenic, no DE miRNA detected from chondrogenic tissue is best explained by either no existence of DE miRNA, or existing DE miRNA was not detected due to lack enough statistical power. Among the miR expression of which was differentially expressed in the osteogenic tissues from adult and old donors, miRs with known function in bone biology were validated: let-7 [52], miR-21 [53], miR-30 [54], miR-96 [55], miR-27 [56], and miR-140 [57]. Interestingly, among predicted targets of these miRs are genes and proteins regulated in MSC from adult and old donors as shown. For example, MMP16, predicted target of miRs miR-27, miR-30 and miR-140, is an important protein regulating bone homeostasis through regulating osteocyte differentiation [58]. Other genes of interest, predicted to be target of more than one of the validated miRs, include members of the ADAMTS family, key to cartilage and bone homeostasis [59], interleukin 18 with important role in bone metabolism [60]. Several genes with a not yet established function in MSC or bone biology, however reported to be expressed in these cells or tissues have also been characterised as DE in this study and are predicted target genes of miRs here validated, for example desmoglein 2 [61].   Interestingly, the expression of all miRs in the osteogenic tissues was downregulated in tissues from older donors. This may be due to defective miRNA biogenesis machinery [62] or decreased ability of the MSCs from older donors to undergo the osteogenic differentiation pathway [63]. Interestingly, the lower levels of expression of the enzyme associated with miR production, Dicer, in MSCs have been associated with their decreased differentiation potential [62].
We have also validated DE of miRs: miR-500 and miR-548j in the tenogenic tissues from young and older donors. It has been shown that miRs may play a role in tendon homeostasis [64,65], however little is still understood about the role of miRs in tenogenic differentiation or tendon homeostasis. Based on target prediction databases), miRs: mir-500 and miR-548j may be regulating processes associated with matrix remodelling which are important in both tendon  Table 9. The five significant canonical pathways related to changes in the methylation patterns for each tissues type. The log (p-value) of each pathway was determined using Fisher's exact test. The ratios were calculated as the number of input molecules mapped to a specific pathway divided by the total number of molecules in the given pathway. formation and maintenance, as well as healing. Interestingly, miR-548j is predicted to target peroxisome proliferator-activated receptor gamma (PPARG), a gene differentially expressed in tenogenic tissues from young and older MSCs donors. PPARG has been shown to be involved in tendon healing [66] further indicating the potential involvement of miR-548j in tendon repair.
To summarise, we have validated DE of miRs and their predicted target genes in the osteogenic and tenogenic tissues from young and older donors that may be associated with the decreased function of MSC from older donors and of relevance to MSC-based therapies.
LncRNAs play important roles in age-related diseases. Evidence is emerging that lncRNAs affect the molecular processes that underlie age-associated phenotypes. LncRNAs modulate gene expression patterns at the transcriptional, post-transcriptional and post-translational Diseases and biological functions identified from the sets of DM loci input into IPA for chondrogenic, osteogenic and tenogenic engineered tissues. A. Heatmap of the top 20 diseases and biological functions identified using IPA comparison analysis with significant activation z scores (infers the activation state of regulation). Scale relates to activation Z scores were green is a positive activation z-score (activated) and red is a negative score (inhibited). B. A cell differentiation network was identified in all engineered tissue types. The network shown includes DM genes identified in tenogenic tissues. C. The network 'congenital anomaly of the musculoskeletal system' was activated in tenogenic tissues. In networks red genes relates to those hypomethylated and green hypermethylated in tissues derived from older MSCs. level. They affect many cellular processes relevant to ageing biology such as proliferation, differentiation and senescence (reviewed [67]). We identified a catalogue of lncRNAs for further work to define their roles in musculoskeletal ageing as although studies suggest the majority are functional only a few established biological relevance [68]. In tenogenic tissues we identified XIST as having a reduced expression in older tissues. XIST, responsible for imprinting controls epigenetic changes through DNA methylation and declines in senescence; though its function in this is unknown [69].
In transcriptomic studies we used gene ontology and network analysis tools to study pathways affected by MSC donor age. However, there are a few interesting findings for the some individual genes. Chondrogenic tissues were the most affected engineered tissue type with age demonstrated by the number of DEGs, whilst tenogenic were the least age-affected engineered tissue. Whilst principally large 'omics' datasets are analysed using network analysis to understand the overall effects of expression changes in the tissue, there were a number of interesting findings at the individual gene level that warrant discussion. In chondrogenic tissues the most DE gene was neurotrophin-3 (NTF3); highly expressed in young. This was not reflected in DML across the gene. NTF3 is an important gene in arthritic processes within the joint [70], produced by inflammatory cells of the nervous system as well as connective tissue [71], with survival-promoting and trophic effects on chondrocytes [72]. There is also a down-regulation of NTF expression in chondrocytes in arthritis [73]. Age-related changes in mouse brain have also been reported [74]. In osteogenic and tenogenic tissues ALX homeobox-1 had the most reduced expression in old similar to ageing MSCs [48]. It is important in skeletal development and we previously demonstrated an increased expression in old tendon [23]. Along with Homeobox (Hox) B7 (lower in old) and Mab-21-like 2 (higher in old) it was DE in all engineered tissue types and ageing MSCs [48]. HOX genes have been implicated in ageing of tissues [75] including tendon [23]. Furthermore HOX genes are required for tissue appropriate regeneration [76] and may be involved in the timing of ageing [49]. HOX-mediated transcriptional memory may reduce stem cell-mediated tissue regeneration [77]. Therefore this has special relevance to tissue engineering and musculoskeletal repair in ageing marking them as an interesting gene for further work in tissue engineering using MSCs. Pathway analysis identified similar age-related changes at the molecular and cellular function level from input DE genes for the functions 'cell death and survival', 'cell morphology', and 'cell growth and proliferation'. This suggests that although the DEG may be different between engineered tissue types the age-related pathways involved at this level are similar. Interestingly in ageing MSCs we demonstrated age-related changes in gene profiles included differences in cell proliferation, signalling, function and maintenance suggesting an age-related loss in MSCs ability to respond to biological cues [48]. Thus these changes seem to have impacted on all classes of engineered tissues.
The top canonical pathways in chondrogenic and tenogenic tissues were related to oxidative stress similar to that previously identified in tenocytes exposed to extracellular glucose [78] and ageing chondrocytes [79]. Mitochondrial dysfunction and oxidative phosphorylation were the most significant in chondrogenic tissues similar to ageing MSCs [48]. One hallmark of ageing is mitochondrial dysfunction and alterations in redox balance could account for the observed reduction in cellular function associated with age [80]. The mitochondria represent an important source of cellular ROS and recent evidence suggests that age-related oxidative stress can disrupt normal physiological cell signalling pathways. Human chondrocytes isolated from older adult chondrocytes exhibited increased peroxiredoxin (PRX) hyperoxidation, particularly for mitochondrial PRXs, when compared to younger chondrocytes in a recent study. PRXs represent a key antioxidant system, and oxidative stress mediated hyperoxidation leads to inhibition of PRX function, and a reduction in antioxidant capacity. PRX hyperoxidation was associated with inhibition of downstream pro-survival signalling and up-regulation of prodeath signalling which led to chondrocyte cell death [79]. Additionally, oxidative stress mediated inhibition of pro-survival cell signalling has also been demonstrated in another recent study supporting the theory that high levels of oxidative stress, such that are observed in ageing tissues, could lead to oxidative stress mediated alterations in cell signalling and contribute to the ageing phenotype and the development of age-related disease such as osteoarthritis [81]. The current study also identified significant alterations in cell signalling pathways in chondrogenic tissues. Thus the use of older MSC donors for cartilage tissue engineering could result in increased oxidative stress within the engineered tissue which could alter normal physiological signal transduction which could have significant consequences for the synthesis and degradation of cartilage matrix components [82]. Peroxisome proliferator-activated receptor signalling (PPAR) are key regulators in various age-related processes related to oxidative stress and energy metabolism. PPAR signalling was the dominant pathway in tenogenic tissues. As PPAR signalling has roles in cell proliferation, differentiation and tissue remodelling [83], and these pathways were also identified through a network of DE genes, this could have detrimental effects on engineered tendon from older MSCs. Furthermore PPAR signalling affects the impairment in mitochondrial biogenesis demonstrated in OA chondrocytes [84].
For osteogenic tissues age-related changes in genes involved in VDR/RXR (vitamin D receptor (VDR)-9-cis-retinoic acid receptor (RXR)) were identified. The classical actions of vitamin D3 are through this signalling pathway facilitating transcription of genes important in bone for the expression of several proteins including osteopontin [85] and in osteoblasts transcription of nuclear factor-kappaB ligand (RANK-L); important for the activation and differentiation of the osteoclasts [86]. A change in VDR expression with ageing has been reported in rat bone [87] and a reduction in ageing mice osteoblasts [88]. A significant effect on cell differentiation and survival is evident following a reduction in VDR activity in bone. Furthermore a decrease in VDR may be partially responsible for increased levels of apoptosis in ageing osteoblasts [88], together with reduction in bone mineralization proteins; osteopontin and osteocalcin [89]. In osteogenesis from ageing MSCs there are alterations in osteocalcin expression with negative effects on proliferation and differentiation capacity of BMSCs in culture [6]. Another study demonstrated osteogenic potential of ageing MSCs was reduced as measured by Alizarin Red staining (which stains calcium deposits) [90]. Together these findings suggest that the reduced osteogenic potential of ageing MSCs could in part be due to a dysregulation of VDR/RXR signalling. The most significant canonical pathways related to changes in the methylation patterns for osteogenic tissues was active AMP-activated protein kinase (AMPK) signalling. AMPK is highly affected by age and may be a crucial cell signalling pathway that regulates cell function. Age related decline in AMPK plays a key role in the observed loss of function, disordered bioenergetics and metabolism observed in ageing cells and likely contributes to age related disease [91]. Indeed OA chondrocytes are deficient in the metabolic biosensors active AMPK [84].
In our previous cartilage ageing equine transcriptomics study [22] we identified an over representation of DEGs (principally reduced in ageing) involved in pathways for extracellular matrix, degradative proteases and matrix synthetic enzymes and cytokines/growth factors and Wnt signalling. In our aged tissue engineered cartilage DEGs principle pathways involved mitochondrial dysfunction and oxidative phosphorylation. In both ageing cartilage and ageing tissue engineered cartilage the pathway hepatic fibrosis was identified with the genes COL2, COL9 and COL24 DE in both data sets. However for cartilage this was evident as reduced expression in ageing whereas in tissue-engineered cartilage they were increase in ageing. Thus the age-related changes in cartilage are not reflected using tissue-engineered cartilage. Our Achilles tendon ageing RNASeq study [23] determined the top networks for DEGs were from cellular function, cellular growth, and cellular cycling pathways. In tissue-engineered tendon the main pathways from DEGs were signalling pathways with large proportion of DEGs being transcription factors. There were no overlaps in the pathways with ageing between tendon and tissue-engineered tendon. It seems that from these results that age-related changes in cartilage and tendon are not reflected at the transcriptome level to those in MSC-derived tissue engineered cartilage and tendon. This could be due to the method/conditions of tissue generation, a more embryonic phenotype in engineered tissues [92](which may be altered with the addition of loading or growth factors [93]) or disparate ageing mechanism's. Thus whilst tissueengineering may provide too artificial a systems to study age-related changes in MSC biology per-se, our findings are important in relation to therapeutic cell source decisions.
In tissues derived from ageing MSCs we identified changes in expression level and differences in the relative balance of splice products. AS; the production of multiple mRNA isoforms from a single gene due to alternative choice of exons or splice sites during pre-mRNA splicing is a major source of protein diversity for higher organisms, and is frequently regulated in a tissue-specific manner. It is estimated that up to 90 percent of human genes have multiple isoforms [94]. Splice variants from the same gene can produce proteins with discrete properties and diverse (including antagonistic) functions. Furthermore, a number of genetic mutations involved in human disease have been mapped to changes in splicing signals or sequences that regulate splicing. Thus, an understanding of changes in splicing patterns is critical to a comprehensive understanding of biological regulation and disease mechanisms. There is a growing interest in the role of AS in normal tissues [95], development [96] and disease (reviewed [97]), but little is known on its role in MSC ageing and tissue engineering. Changes in AS may have a major impact on cell survival [98] and post-translational modifications [99].
Our study demonstrated that donor MSC age has an effect on splicing events in all engineered tissue types, similar to an ageing study in peripheral blood leukocytes [100], suggesting that modification of mRNA processing may be a feature of human ageing. GO ontology demonstrated enrichment in genes associated principally with metabolic processes in genes undergoing AS in all engineered tissue types. AS in metabolic processes is a frequent phenomenon in diseases such as cancer [101] but also in ageing brain [102] and MSCs [48]. Further analysis of the tenogenic AS genes identified carbohydrate and lipid metabolism as significant metabolic pathways affected in ageing. In tenogenic tissuesthe principal network identified in IPA was cell signalling, interaction function and maintenance. This suggests that similar to some cancer cells [98] in tendon tissues derived from ageing MSCs there may be expression of isoforms that alter cell survival. We have previously observed an age-associated disruption to the balance of alternatively expressed isoforms for selected genes in tendon ageing [23]. In AS affected genes for chondrogenic tissues with age were related to cell death and survival, and growth and proliferation of connective tissue. Interestingly a recent study found pyruvate dehydrogenase kinase isoform 2-mediated alternative splicing switches hypoxia-inducible death protein from cell death to cell survival in cancer cells [103].
Though there was little overlap between DEG and DML those that were displayed a good correlation of DNA methylation with differentially expressed genes (promoter increased methylation; reduced gene expression, gene body increased methylation; increased gene expression). There are a number of possibilities as to why correlation between gene expression and DML was poor. DNA methylation is stable so the methylation changes evident may be associated with ancestral gene expression differences. For instance in the study we differentiated MSCs for 21-28 days. There may have been a gene expression changes between day 0 and 7 accompanied by a methylation change. The gene expression could then return back to basal day 0 levels at day 28 when gene expression was measured, but the methylation change remains. Thus a DNA methylation change may not be accompanied by gene expression change as although gene expression did alter, it is not different at the time points measured. Another possibility is that gene enhancers can be located within the gene body of a different gene. Thus a DML within the gene body of gene A may actually be within the enhancer of gene B and so the DML may be associated with a change in gene B but not gene A in which we assessed the effect of the DML on gene A. Finally gene body methylation can be associated with alternative splicing and transcription from alternative/cryptic transcription start site. These may not be have been detected in our RNASeq analysis.
DNA methylation provides a stable and heritable gene regulatory mechanism enabling cells to alter the expression of many genes. Alterations in DML have previously been seen in ageing tissues [104] as a global hypomethylation, as well as aged MSCs [48,49]. Significant age-related DMLs, both tissue specific and common were identified in all tissue types. DNA methylation also has a role in genomic imprinting (the epigenetic occurrence by which certain genes are expressed in a parent-of-origin-given manner) by regulating the differential expression of maternal and paternal imprinted genes. It is also important in DNA damage/repair and genomic instability [105]. Furthermore a number of disease have been associated with aberrant DNA methylation (reviewed [106]). Thus alterations in methylation in engineered tissues may have further consequences than gene expression changes alone. For instance altered methylation may affect the DNA damage response resulting in senescence and apoptosis [107]. Further work is required to decipher the impact of the DNA methylation changes evident in this study on DNA damage and genomic instability in musculoskeletal engineered tissues with ageing.
Methylation has dual and opposing roles in the gene expression regulation. In promoter regions DNA methylation is correlated with transcriptional repression, while in gene bodies it is generally associated with high expression levels [49]. This paradox emphasizes the diverse involvement of methylation in specific genomic and cellular contexts. We described common age-related pathways from DML of skeletal system morphogenesis, regulation of transcription regulation (both principally due to DML in HOX genes) and cell proliferation (also identified in RNASeq data).
The network 'skeletal and muscular development and function' significant in all tissues due to DMLs was principally due to the enrichment of homeotic or HOX genes, similar to studies in ageing MSCs [24,49]. HOX gene expression is tightly regulated temporally during vertebrate development. An association between HOX genes and longevity has been previously proposed [108]. HOXB7 (chondrogenic), HOXB6 and HOXB7 (osteogenic) and HOXA3 and HOXB6 (tenogenic) were also DEGs in the RNASeq study. However, for HOXB7 (osteogenic) and HOXA3 (tenogenic) DEG did not correlate with DML position. Despite this, the dysregulation of HOX genes at the mRNA and epigenetic level consolidate their role in both MSC ageing and in engineered tissues from ageing MSCs.
An interesting signature of age-related DML was that between 40% (chondrogenic and osteogenic) and 50% (tenogenic) of the top 20 most DML for each engineered tissue type were transcription factors. DML not only transcriptionally regulates gene expression but in ageing we identified significant regulation of the transcription factors. In our study of ageing tendon we identified an overrepresentation of DE transcription factors in ageing tendon [23] suggesting they may have a pivotal role in tendon ageing and tissue engineered tendon. Conversely in our cartilage ageing study there were few age-related DE transcription factors. Furthermore many age-related DML identified in all engineered tissue types resulted in 'differentiation of cells' being highlighted. Site specific CpG site methylation changes can affect MSC chondrogenic differentiation [109] whilst alterations in MSC potential have been previously noted [110]. However 'cell survival' networks relating to DML were affected in chondrogenic and osteogenic tissues only suggesting distinct age-related methylation patterns between tissue types with potential distinct consequences to engineered tissues.
There was an age-related activation in the function differentiation in all tissues due to DMLs. The DM at some of these loci could be responsible for changes in differentiation potential previously described in chondrogenesis and osteogenesis from ageing MSCs. Studies have described a more rapid decline in differentiation potential for osteoblastic and chondrogenic lineages relative to adipogenic differentiation [90,111]. No studies have investigated tenogenic differentiation potential with age. Our results would suggest based on the amount of gene expression changes, AS and DML alterations that tenogenic differentiation is likely to be more maintained with ageing compared to chondrogenic and osteogenic lineages. However this hypothesis needs validation. Additional parameters such as histone modification would benefit the analysis but this was beyond the financial scope of the study. Dynamic histone methylation can contribute to the ageing process through influencing DNA repair and transcriptional regulation of ageing processes (reviewed [112]).

Conclusion
Context-dependent DNA methylation plays a critical role in regulating gene transcription, thereby serving as an important epigenetic marker or regulator in many biological activities. Taken together our RNASeq and DNA methylation results in engineered tissues suggest an age-related loss in the differentiated cells ability to respond to biological cues. Age-related cellular dysfunctions have been hypothesized to results in musculoskeletal age-related diseases such as OA, osteoporosis and tendinopathy. These results have significant implications for therapeutic cell source decisions (autologous or allogeneic) revealing the necessity of approaches to improve functionality of ageing MSCs.
Supporting Information S1 File.

Acknowledgments
We thank Luke Tregilgas for transcript stabilisation analysis, Sophie Sutton for assistance with the gene expression validation and analysis and Aaron Dale and Louise Reynard for critical evaluation of the manuscript. We also thank Danae Zamboulis for her assistance with image production.