Introduction

Differentiation of mesenchymal stem cells (MSCs) into chondrocytes occurs by the process of chondrogenesis¹. This is an important process during development, as it is a pre-requisite for skeletogenesis. Chondrocytes perform the vital role of generating cartilage during embryogenesis, but also maintain cartilage throughout life, including at the ends of long bones in articulating joints². The master transcription factor responsible for chondrogenesis is SOX9 and during this process chondrocytes secrete anabolic proteins such as type II collagen and aggrecan, encoded by COL2A1 and ACAN, which constitute a significant functional portion of the cartilaginous extracellular matrix (ECM) ^{2, 3}. With increasing age and influenced by a mixture of factors such as (epi)genetics, obesity and mechanical injury/stress, chondrocytes will increasingly express catabolic proteins such as matrix metalloproteinases which degrade the cartilage ECM³. Ultimately the chronic loss of cartilage contributes to the extremely debilitating condition osteoarthritis (OA), which remains incurable with treatment options limited to pain relief medication and end-stage joint replacement surgery ^{4, 5}.

miRNAs are short non-coding RNAs - roughly 19-22nt long, that negatively regulate target mRNAs ^{6, 7}. Mammalian miRNA-mRNA interactions occur via complementary sequence specificity between a 7-8nt stretch of the miRNA known as the seed sequence, and positions along the 3’UTR of the target mRNA ^{8, 9}. miRNA-mRNA interactions are complex because a single mRNA can be targeted by many miRNAs, and a single miRNA can target many mRNAs. Our previous work identified significantly altered miRNAs during chondrogenesis and demonstrated the role of miR-140 in regulating chondrocyte gene expression ¹⁰. miR-140 has emerged as a vital regulator of chondrogenesis, cartilage and OA, and have been hypothesised as effective drug targets due to their pro-chondrogenic regulation ^11–14. Additional miRNAs, such as miR-455, have also been shown to be important regulators in maintaining healthy cartilage and protecting against OA development ^{15, 16}.

To identify further important miRNAs from this dataset, and to overcome the complexity of miRNA-mRNA interactions, we performed a combined bioinformatic, experimental and systems biology approach to better understand the relationship between miRNAs which may be important to chondrogenesis, and their predicted targets. We identified miR-199a-5p and miR-199b-5p as pro-chondrogenic miRNAs. Just as with miR-140-5p and miR-455, we anticipated miR-199a-5p and miR-199b-5p (miR-199a/b-5p) to target multiple mRNAs during chondrogenesis. To this end we used experimental, informatic and literary data to build kinetic models to explain how miR-199a/b-5p regulated chondrogenesis. The models included three targets identified through RNAseq analysis (FZD6, ITGA3 and CAV1), and were used to make predictions to fill experimental gaps and predict novel interactions between miR-199a/b-5p and the chondrogenesis machinery..

Results

Analysis with TimiRGeN identified miR-199b-5p to be upregulated during chondrogenesis

Our previously generated 14-day chondrogenesis time series dataset was analysed with the TimiRGeN R/ Bioconductor package - a novel tool we developed to analyse longitudinal miRNA-mRNA expression datasets ^{10, 17}. Prior to using TimiRGeN, the transcriptomic data underwent timepoint based differential expression analysis, using the day 0 time-point (MSCs) as the control group for chondrogenesis samples measured at days: 1, 3, 6, 10, 14, and thus created differential expression data over five time-points. TimiRGeN identified signalling pathways of interest over the 14-day time series data (Supplementary Table 1). Eight signalling pathways were found enriched in at least three of the five time points (Figure 1). The TimiRGeN pipeline was then used to predict miRNA-mRNA interactions that may regulate each of the eight pathways. miRNA-mRNA interactions were kept if the miRNA and mRNA involved in the interaction has a Pearson correlation <-0.75 across the time-series and if the interaction was predicted in at least two of the following three databases: TargetScan, miRDB or miRTarBase ^18–20. To identify which miRNAs involved in the miRNA-mRNA interactions were positively changing in each of the eight signalling pathways we scaled the Log2FC values from limma ²¹. By using scaled Log2FC values, we could highlight magnitude of change, rather than total change.

Figure 1.

miR-140-5p was the most positively changing miRNA in the following seven pathways: Adipogenesis, Clear cell renal cell carcinoma pathways (CCRCCP), Epidermal growth factor/Epidermal growth factor receptor (EGF/EGFR), Endochondral ossification, Endochondral ossification with Skeletal Dysplasia (EOSD), Gastrin signalling pathway and Vascular endothelial growth factor-A/Vascular endothelial growth factor receptor 2 signalling pathway (VEGFA-VEGFR2). miR-199b-5p was the second most positively changing microRNA in the following six pathways: CCRCCP, EGF/EGFR, Endochondral ossification, EOSD, Gastrin signalling pathway and VEGFA-VEGFR2. Other miRNA/genes such as hsa-miR-455-5p and BMP2 were also of interest, but we focused on hsa-miR-199b-5p and its paralogue miR-199a-5p since these are comparatively under-researched miRNAs within the context of chondrogenesis ¹⁰.

Expression of chondrogenic biomarkers and glycosaminoglycan levels change over time after altering miR-199a-5p or miR-199b-5p expression

To identify if miR-199a-5p and miR-199b-5p regulate chondrogenesis we performed MSC chondrogenic differentiation and measured the consequences of miR-199 overexpression and inhibition (Supplementary Figure 1) on chondrogenic biomarkers ACAN, COL2A1 and SOX9 and glycosaminoglycan (GAG) levels. Significant upregulation was seen in ACAN, COL2A1 and GAG levels following miR-199a-5p overexpression (Figure 2A). In contrast, when miR-199a-5p expression was inhibited ACAN, COL2A1 and SOX9 were significantly downregulated at multiple time points (Figure 2B). GAG levels were also significantly decreased, by ∼40%, at day 7. When miR-199b-5p was inhibited, similar significant downregulation of ACAN, COL2A1 and GAG levels occurred (Figure 2C). Inhibition of miR-199a-5p and -199b-5p together caused a more consistent reduction of ACAN, COL2A1 and SOX9 expression.

Figure 2.

FZD6, ITGA3 and CAV1 are the most significantly upregulated miR-199a/b-5p targets during MSC chondrogenesis

Multiple mRNA targets of miR-199a/b-5p are likely regulated during chondrogenesis. To elucidate which genes are most affected by miR-199 inhibition we performed an RNAseq experiment to identify candidate targets during the early phase of MSC chondrogenesis comparing MSC samples which were undifferentiated (day 0) and MSC samples which were in the early stages of chondrogenesis (day 1). Chondrogenic markers were downregulated by miR-199 inhibition in the early stages of chondrogenesis (Figure 1), so we reasoned that identifying mRNAs which are positively enriched during the first few days of chondrogenesis may identify the most important mRNA targets of miR-199a/b-5p. We chose to inhibit miR-199a-5p or miR-199b-5p since supraphysiological overexpression of miRNA mimics can lead to spurious findings ²². Initial comparison between day 1 differentiated and day 0 undifferentiated control samples identified 4391 upregulated and 4168 downregulated significantly differentially expressed genes (DEG) (< 0.05 adjusted p value). Positively changing genes included COL2A1 (log2fc =11.6), ACAN (log2fc = 9.18) and SOX9 (log2fc = 3.37). Gene Ontology (GO) term analysis on the upregulated genes confirmed that the cells were differentiating towards chondrocytes with terms such as skeletal system development, extracellular matrix organization and regulation of cartilage development significantly (adjusted P values < 0.05) enriched. To identify miR-199-regulated genes, undifferentiated control MSC samples were contrasted against undifferentiated MSCs with either miR-199a-5p or miR-199b-5p inhibition, which respectively resulted in 87 and 46 significantly DEGs (adjusted P values < 0.05) (Figure 3A). Similar comparisons at day 1 of chondrogenesis revealed inhibition of miR-199a-5p or miR-199b-5p respectively resulted in 674 and 817 DEGs. 25 and 341 genes intersected between day 0 and day 1 chondrogenesis inhibition experiments, indicating that both microRNAs may share functional repertoire of targets (Chi-squared observed vs expected day 0 = 0.003477, day 1 = 0.000624). COL2A1, ACAN and SOX9 were significantly lower in expression in day 1 samples after miR-199 inhibition, validating the negative impact on chondrogenesis. GO term analysis was conducted to identify biological processes linked to miR-199a/b-5p inhibition (Figure 3B). Interestingly several terms associated with chondrogenesis were suppressed/less activated during day 1 analyses, including: extracellular matrix constituents, skeletal system morphogenesis, extracellular matrix and collagen containing extracellular matrix.

Figure 3.

The significantly differentially expressed genes were also analysed using the miRNAtap R package to predict miR-199a-5p and miR-199b-5p targets (Figure 3A). 21 potential mRNA targets were predicted to be targeted by both miRNAs, and in alphabetical order they were: ABHD17C, ATP13A2, CAV1, CTSL, DDR1, FZD6, GIT1, HIF1A, HK2, HSPA5, ITGA3, M6PR, MYH9, NECTIN2, NINL, PDE4D, PLXND1, PXN, SLC35A3, SLC35D1, VPS26A. FZD6, ITGA3 and CAV1, MYH9, NECTIN2.The expression patterns of these 21 miR-199 target genes were also explored in our microarray study (Supplementary Table 3) ¹⁰. SkeletalVis was used to contrast the behavior of the 21 genes in other MSC derived chondrogenesis studies (Supplementary Table 4) ^23–25. The tables showed many of genes found through our bioinformatic analysis to be significantly downregulated during chondrogenesis in multiple studies, and this included CAV1, FZD6, ITGA3, and MYH9. We decided to further explore if FZD6, ITGA3 and CAV1 were true targets of miR-199a/b-5p because in terms of adjusted p values, these three showed the most significant change (Figure 3B)

miR-199a/b-5p targets were upregulated by miR-199a/b-5p inhibition

To confirm the effects miR-199a-5p and miR-199b-5p had on the targets identified through RNAseq, we performed a series of inhibition experiments. Since the miRNAs shared the same seed site sequence (nucleotides 2-8; 5’-CCAGUGUU-3’), we tested if inhibition of one miRNA would lead to similar or different effects to inhibition of both miRNAs. For this, we picked to supress expression of miR-199a-5p as it was the more highly expressed of the two paralogues and contrasted it to inhibition of miR-199a/b-5p. We tested the effects of the inhibition by measuring the expression of the three most significantly enriched predicted miR-199a/b-5p targets from our bioinformatic analysis – FZD6, ITGA3 and CAV1 (Figure 4A1). We saw significant upregulation of FZD6, ITGA3 and CAV1 after miR-199a-5p inhibition and miR-199a/b-5p inhibition. There was no consistent improvement of the combination of miR-199a-5p and miR-199b-5p over miR-199a-5p alone. We also tested the third most significant predicted target of miR-199a/b-5p - CAV1, and how its’ expression changed after miR-199b-5p inhibition over the same time course, and we saw CAV1 levels were significantly increased post miR-199b-5p inhibition ²⁶(Figure 4A2).

Figure 4.

Unlike FZD6 and ITGA3, CAV1 has previously been established to be a direct target of miR-199a-5p ²⁶. To identify if FZD6 and ITGA3 were also direct miRNA targets we cloned the 3’UTR regions of FZD6 and ITGA3 directly downstream of the firefly luciferase gene and demonstrated reduced expression compared to empty control vector, suggesting these contain potentially repressive elements (Figure 4B1). Introduction of an excess of miR-199a-5p into the cells further repressed expression suggesting both FZD6 and ITGA3 3’UTRs are direct miR199a-5p targets. This was confirmed by mutation of the predicted miR-199-5p seed sequence within the 3’-UTR of each gene which reduced the extent of inhibition of luciferase levels caused by the miRNA (Figure 4B2).

Kinetic modelling creates an in silico demonstration of how miR-199a-5p/miR-199b-5p regulates chondrogenesis

We attempted to capture the complexity of how miR-199a/b-5p regulated chondrogensis using an in silico model (Figure 5A). Only using the experimental data presented in this manuscript, and our previous microarray work, we developed a model to demonstrate the relationships between chondrogenic biomarkers and the targets of miR-199a/b-5p we identified through RNAseq and subsequent inhibition experiments. Finally, we used events to simulate inhibition of miR-199a-5p or miR-199b-5p (Figure 5B). The experiments shown in figures 2A, 2C and 4A were used to parameterise the model, and since the experiments were performed in a staggered manner, we can use the model to make predictions to fill experimental gaps. Using this model, we predicted the dynamics of the chondrogenesis biomarkers and CAV1 after miR-199a-5p inhibition and the dynamics of FZD6 and ITGA3 after miR-199b-5p inhibition. Most objects within the model were based on experimental data, and the differences between the experimental data and simulated data are calculated by Mean Squared Error (MSE). In the initial model, 15/18 modelled objects with experimental data had an MSE of lower than 3, indicating that most of the experimental data was captured by the model, and the average MSE for the model was 15.96.

Figure 5.

To enhance our initial model, we added further detail to increase biological relevance. This included the addition of Transforming Growth Factor Beta (TGFB) as a trigger for chondrogenesis initiation – just as we had during our chondrogenesis MSC experiments. We found strong mechanistic links between TGFB signalling and CAV1 expression via SRC kinase. Based on our previous work we also included miR-140-5p as it has been proven as a vital regulator of chondrogenesis, and with our kinetic model we predicted how miR-199a-5p or miR-199b-5p inhibition indirectly affected miR-140-5p. Further to this, we identified several flaws in our initial model, which we attempted to rectify using the enhanced chondrogenesis model. Firstly, by simulating the miR-199a-5p and miR-199b-5p inhibition to last until day 7, we saw that SOX9 mRNA, COL2A1 mRNA and ACAN mRNA dynamics between days 1 and 7 of miR-199b-5p inhibition were flat. Furthermore, our CAV1 dynamics were also flat during miR-199b-5p inhibition. Figure 2C shows, from miR-199b-5p inhibition, the chondrogenesis biomarkers had a sharp decrease followed by a steady rise until day 7, and to match the effect of the miR-199b-5p inhibition experiments, we simulated the inhibition to last until day 4.5 instead of day 7. Doing so increased the similarities between our experimental data and simulations. Secondly, we wanted to include other miR-199a/b-5p targets which were alluded to in our RNAseq experiment, but not further explored, such as MYH9 and PDE4D. We added a blackbox named ‘OtherTargets’ to represent other targets which miR-199a/b-5p regulate during chondrogenesis. Also, to delay the decrease in ACAN, COL2A1 and GAG levels after miR-199a/b-5p inhibition, we slowed down the interactions between the miR-199a/b-5p targets and SOX9 by including SOX9 protein and SOX9 phospho-protein as objects in the model (Figure 5C, Figure 5D). In the enhanced model, 14/18 (77.7%) of the objects with experimental data had an MSE of over three, indicating the model – even with additional data from the literature, still captured much of the experimental data. Also, the average MSE was 12.08, indicating an improvement over the initial model.

From the enhanced model we saw improved dynamics for several model objects e.g., GAG, CAV1, SOX9, ACAN, COL2A1. Our MSE values which were used to quantify how similar our models’ simulations were in contrast to our experimental data. Generally, our MSE values were better in our enhanced model for the chondrogenesis biomarkers and GAG levels, though our MSE values were better in our initial model for the miR-199a/b-5p targets. Overall, we saw the dynamics improved in the enhanced model, and in addition, we could make predictions on how miR-140-5p levels would indirectly be influenced from knock-down of miR-199a-5p or miR-199b-5p and such a multi-miRNA model has not yet been created before. Further pathways such as the Wnt signalling pathway may also play an important role in this system but could not be explored via modelling with our current level of data. We developed a hypothesis rich GRN to display further pathways which were not in the scope of this project (Supplementary Figure 3).

Discussion

Our initial work in this area used a combined experimental and bioinformatics approach to identify and study the roles of miR-140 and miR-455 which were highly important to chondrogenesis ^{10, 13}. This current work extends this, to identify other miRNAs of regulatory importance during chondrogenesis utilising a recently developed tool (TimiRGeN R/Bioconductor package)¹⁷. Here, we combined experimental, bioinformatic and systems biology approaches to identify and study the role of miR-199a/b-5p during chondrogenesis. The MIR199 family were identified as pro-chondrogenic in mouse mesenchymal C3H10T1/2 cells, reportedly by targeting SMAD1 (miR-199a-3p) and JAG1 (miR-199a/b-5p), respectively ^{27, 28}. Having previously shown the substantial upregulation of miR-199a/b-5p during chondrogenesis, we now show through both loss- and gain of function experiments that miR-199a/b-5p also promotes SOX9, COL2A1 and ACAN expression in human MSC chondrogenesis thus enhancing cartilage formation ¹⁰. RNAseq analysis identified key targets of miR-199a/b-5p including FZD6, ITGA3 and CAV1, which were experimentally validated and subsequently incorporated into an in silico kinetic model.

miR-199a and miR-199b are vertebrate-specific miRNAs which exhibit an expression pattern associated with mesenchymal tissues and the skeleton ²⁹. Antisense transcript Dmn3os encodes miR199a and miR214 and its deletion in mice causes skeletal defects including short stature and cranial deformity ³⁰. Similar phenotypic consequences caused by the loss of DNM3OS and therefore MIR199A and MIR214 have also been reported in humans ³¹. miR-214 has since been reported to negatively impact on chondrocyte differentiation, through targeting ATF4 ³². Thus, our demonstration of the pro-chondrogenic nature of miR-199a-5p and miR-199b-5p in human MSCs, in addition to the previous demonstration in mouse C3H10T1/2 cells, further supports the role of these miRNAs in chondrogenesis and formation of the skeleton.

We show for the first-time functional evidence for the miR-199b-5p-CAV1 interaction to occur in human MSCs. Several mechanisms reported in the literature may also support how miR-199a/b-5p-CAV1 could regulate chondrogenesis, such as TGFB triggering phosphorylation of CAV1 via SRC-kinase ²⁶. In contrast to CAV1, FZD6 and ITGA3 have been less well studied as miR-199a/b-5p targets though, miR-199a-5p-FZD6 has been predicted previously^{28, 33}. Our luciferase assays validate for the first time that FZD6 is a target of miR-199a-5p. Previous work by our group has validated FZD6 as a target of miR-140-5p, so it is likely a highly important gene in chondrogenesis. ITGA3 meanwhile has been found to be a target of miR-199a-5p in neck and head cancer cells, and our results confirm this interaction in humans MSCs ³⁴.

The presented model topology (Figure 5A, 5C) was based on our experimental work. Our initial model (Figure 5A) was re-worked to better match the experimental patterns seen and to include additional genes which we can make prediction from. The enhanced chondrogenesis model (Figure 5C) was initiated with TGFB3, which was used to initiate chondrogenesis in our experiments and was therefore used as the proxy for chondrogenic initiation. As such, TGFB3 promoted miR-199a/b-5p levels to match our microarray data, though based on our previous data TGFB may not directly regulate miR-199a/b-5p levels ¹⁰. TGFB3 also induced CAV1 via SRC kinase and induced SOX9 via the SMAD2-SMAD3 pathway, which then stimulated SOX9 protein production (SMAD2/3 were not included in the models) ^35–37. However, our microarray data showed CAV1 gene expression decreased during early chondrogenesis, but then increased again – perhaps indicating CAV1 has a smaller negative regulative effect on chondrogenesis, or only effects early chondrogenesis. Gene expression of FZD6, ITGA3, CAV1 and the OtherTargets blackbox all had inverse relationships with miR-199a-5p and miR-199b-5p, and FZD6 was also negatively regulated by miR-140-5p. SOX9 mRNA, SOX9 and P-SOX9 were all treated as separate objects in this model. P-SOX9 promoted COL2A1, ACAN and miR-140-5p expression. ACAN contributed directly, and to the greatest extent, to GAG levels (since most cartilage GAGs are post-translationally added to Aggrecan, the protein product of ACAN) which was used as the phenotypic level output for the model and a proxy for chondrogenesis progression ^{10, 38–40}. Precise mechanisms of how FZD6 and ITGA3 regulated chondrogenesis are unclear, with potentially implicated pathways included in the larger GRN (Supplementary Figure 1). FZD6 is a transmembrane protein which functions as a receptor for WNT signalling proteins to, most commonly, activate the non-canonical planar cell polarity (PCP) pathway. However, elements of the Wnt signalling pathway have been implicated to act antagonistically to chondrogenesis ^{41, 42}. ITGA3 is a cell surface integrin which forms a heterodimer with ITGB1 to form α3β1 heterodimers through which chondrocyte-fibronectin ECM connections can be created ^{43, 44}. Conditional deletion of Itgb1 in cartilage impacts profoundly on skeletogenesis in mice ⁴⁵. α3 integrins have been found to increase in OA cartilage, though a direct mechanism between ITGA3 regulating SOX9 is not clear. Regulation of SOX9 via CAV1 has been more studied, which has shown that the relationship between SOX9 and CAV1 is complex and requires further testing. CAV1 may be affecting SOX9 levels, via its activation of RHOA/ROCK1 signalling, which leads to phosphorylation of SOX9. RHOA/ROCK1 inhibition has been shown to both increase levels of SOX9 and chondrogenesis biomarkers in mouse ATDC5 cells, but also decrease levels of SOX9 and chondrogenesis biomarkers in 3D-chondrocytes ^{36, 46–48}. It is likely that other mRNA targets of miR-199a/b-5p also contributed toward chondrogenesis regulation, such as the other genes identified as miR-199a/b-5p targets from the RNAseq analysis e.g. MYH9, NECTIN2. Based on this limitation, a “blackbox” called OtherTargets was added to the kinetic model to represent all other anti-chondrogenic miR-199a/b-5p targets during chondrogenesis which were not explored in this study. A major limitation of the kinetic models is that we were not able to provide any multi-omic data – as no protein level, or phosphor-protein data were available. We generated a broader GRN (Supplementary Figure 3), which showcases how the models we created could be enhanced with such data.

While the in silico model can serve as a resource for researchers interested in this system, there were certain genes such as CAV1 and FZD6 which proved difficult to model. At the time of building the models we assumed, miR-199a-5p inhibition would lead to a greater effect than miR-199b-5p inhibition, due to miR-199a-5p being more abundant. However, we clearly see now this was a misconception and miR-199b-5p inhibition leads to a greater decrease in GAG levels and a greater increase in FZD6, ITGA3 and CAV1 levels. This could be because miR-199b-5p increases by a greater magnitude than miR-199a-5p (Figure 1B), therefore the inhibition of miR-199b-5p has a bigger effect on chondrogenesis. Such timeseries analysis would have been unavailable by only using differential expression analysis, and the use of the TimiRGeN R package was highly useful in finding this microRNA during our re-analysis.

Our results validate miR-199a/b-5p interacting with FZD6, ITGA3 and CAV1 and for miR-199a/b-5p to provide vital pro-chondrogenic regulatory effects, as observed previously for miR-140 and miR-455. Deletion of miR-140 in both humans and mice affects skeletal development ^49,50. Further, from in vivo mouse models miR-140 and miR-455 were additionally shown to be pivotal in protecting from OA. Recently, intra-articular injection of a miR-199a-5p mimic has been shown to reduce cartilage damage in a rat post-traumatic OA model 51. Mouse genetic studies examining the loss of both -199a/b-5p, specifically in cartilage, are required to better understand the function of these miRNAs in skeletogenesis and chondrocyte development. Our results demonstrate early interest and provide a detailed kinetic model to aid researchers interested in this important topic.

Conclusion

Our combined bioinformatic, laboratory, and systems biology approach was a multi-faceted approach to explore miR-199a-5p and miR-199b-5p as pro-chondrogenic regulators. Based on our bioinformatic analysis, the three most significantly positively changing miR-199a/b-5p were predicted targets FZD6, ITGA3 and CAV. Laboratory experiments validated these as direct miR-199a/b-5p targets and confirmed that miR-199a/b-5p positively regulate chondrogenesis. However, the complex nature of miRNA function means there are likely multiple mRNA targets of miR-199a/b-5p which may work synergistically to modulate chondrogenesis. The GRN and kinetic models were created to capture the behaviours of the system which act as a useful resource for further experimental design.

Methods

Data Processing and differential expression analysis

mRNA and miRNA data were produced by Illumina and Exiqon microarray technologies respectively ¹⁰. mRNA data was processed using the lumi R package and the miRNA data was processed using the affy R package ^52,53. limma was then used to perform pairwise differential expression (DE) analysis ²¹. Here the zero-time point was used as the common denominator for all DE analyses. The time point based DE analyses were: D1/D0, D3/D0, D6/D0, D10/D0 and D14/D0. This type of time point based DE analysis was ideal for pair-wise differential expression analysis aproach, as explain by Spies et al (2019) ⁵⁴. Genes which were significantly differentially expressed (BH adjusted P values <0.05) in at least one of the DE analyses have their adjusted P values and log2FC values were kept for analysis by the TimiRGeN R package. All data wrangling and processing took place in R.

Analysis with the TimiRGeN R package

Dataframes containing mRNA and miRNA DE results are analysed using the combined mode of TimiRGeN ¹⁷. The threshold for timepoint specific gene filtration was set as <0.05 and adjusted P values were used for filtration. Timepoint specific pathway enrichment used microarray probe IDs as the background dataset. Our data was from microarrays so for more accurate Overrepresentation analysis we required the microarray specific gene lists to use as the background set of genes. GPL10558 (mRNA) and GPL11434 (miRNA) were downloaded from GEO and the probes were annotated with entrez gene IDs using BiomaRt ^55,56. Eight signalling pathways were enriched for at least three of the five time point based DE analyses performed. miRNA-mRNA interactions were retained if they had a negative Pearson correlation of < −0.75 and if they were identified two of the three target databases used by TimiRGeN.

Human bone marrow MSC culture

Human bone marrow MSCs (from seven donors, 18-25 years of age) were isolated from human bone marrow mononuclear cells (Lonza Biosciences) and cultured and phenotype-tested as described previously ¹⁰. Experiments were performed using cells between passage 2 and 10, and all experiments were repeated with cells from a minimum of three donors.

Chondrogenic differentiation

MSCs were re-suspended in chondrogenic culture medium consisting of high-glucose DMEM containing 100 µg/ml sodium pyruvate, 10 ng/ml TGFβ3, 100 nM dexamethasone, 1× ITS-1 premix, 40 µg/ml proline and 25 µg/ml ascorbate-2-phosphate. 5×10⁵ MSCs in 500 µl chondrogenic medium were pipetted into 15ml falcon tubes or 5×10⁴ MSCs in 150 µl chondrogenic medium were pipetted into a UV-sterilised V-bottom 96-well plate, and then centrifuged at 500 g for 5 min. Media were replaced every 2 or 3 days for up to 7 days.

Dimethylmethylene blue assay

Chondrogenic pellets and transwell discs were digested with papain (10 U/ml) at 60°C ⁵⁷. The sulphated glycosaminoglycan (GAG) content was measured by 1,9-dimethylmethylene blue (DMB) binding (Sigma) using chondroitin sulphate (Sigma) as standard ⁵⁸.

RNA and miRNA Extraction and Real-Time Reverse Transcription PCR

MSC chondrogenic pellets were disrupted in Ambion Cells-to-cDNA II Cell Lysis buffer (for real-time RT-PCR) or mirVana miRNA Isolation Kit Phenol (for RNA-seq) (both Life Technologies) using a small disposable plastic pestle and an aliquot of Molecular Grinding Resin (G-Biosciences/Genotech). Total RNA was then extracted and converted to cDNA using MMLV reverse transcriptase (Invitrogen) and TaqMan real-time RT-PCR was performed, and gene expression levels were calculated as described previously ⁵⁹. Primer sequences and assay details can be found in Supporting Information Materials. For single miRNA-specific analysis, RNA was reverse-transcribed with Applied Biosystems TaqMan MicroRNA Reverse Transcription Kit (Life Technologies) and real-time RT-PCR performed with TaqMan MicroRNA assays (Life Technologies). All values are presented as the mean ± SEM of replicates in pooled experiments. For experiments with multiple MSC donors statistical testing was performed using a matched paired two-tailed Student’s t test on log-transformed values to account for non-normal distribution. Primer details are within Supplementary Table 5.

miRNA Mimic/Inhibitor Transfection in MSC

For modulation of miR-199 levels in MSC, Dharmacon miRIDIAN mimics (C-300607) or miRIDIAN hairpin inhibitors (IH-300607) were transfected into 40%–50% confluent MSC using Dharmafect 1 lipid reagent (Horizon Discovery) at 100 nM final concentration. Analysis was performed in comparison with Dharmacon miRIDIAN miRNA mimic nontargeting Control #2 (CN-002000-01) or Dharmacon miRIDIAN miRNA hairpin inhibitor nontargeting Control #2 (IN-002005-01). For all experiments, cells were subject to a single transfection prior to induction of MSC differentiation.

Cloning and Plasmid Transfection in SW1353 Cells

SW1353 cells were purchased from ATCC Full length miRNA target 3’UTRs were amplified from human genomic DNA using PCR primers (Supplementary Table 4) to enable Clontech In-Fusion HD cloning (Takara Bio Europe, Saint-Germain-en-Laye, France) into the pmirGLO Dual-Luciferase miRNA Target Expression Vector (Promega, Southampton, U.K.) following the manufacturer’s instructions. Mutation of the miRNA seed-binding sites was performed using the QuikChange II Site-Directed Mutagenesis Kit (Agilent Technologies) (Supplementary Figure 2). All vectors were sequence verified. SW1353 chondrosarcoma cells were plated overnight in 96-well plates at 50% confluence (1.5 x 10⁴ cells/cm²). Cells were first transfected with 3′UTR luciferase constructs (10 ng) using FuGENE HD transfection reagent (Promega) for 4 hours then transfected with Dharmacon miR-199a-5p mimic (50 nM) or miRNA mimic nontargeting Control #2 using Dharmafect 1. After 24 hours of transfection, the SW1353 cells were washed and lysed using Reporter Lysis Buffer (Promega) and firefly and renilla luciferase levels determined using the Promega Dual-Luciferase Reporter Assay System and a GloMax 96 Microplate Luminometer (Promega).

RNAseq

RNA isolated as described above was quality assessed with the Agilent Technology 4200 TapeStation system using an RNA screentape assay (Agilent). cDNA libraries were generated using the Illumina TruSeq Stranded mRNA protocol, combinatorial dual index adapters were used to multiplex/pool libraries. Single-read sequencing, 76 cycles (75+1 cycle for index sequences), on an Illumina NextSeq500 instrument using a high-output 75 cycle kit. Kallisto was used for alignment free RNAseq processing ⁶⁰. Tximport was used to import the RNAseq data into R for further processing ⁶¹.

RNAseq DE and miRNA target identification

Time-matched MSC miRNA inhibition (miR-199a-5p and miR-199b-5p) and MSC controls were contrasted by DE using DESeq2 ⁶². miRNAtap was used to score and identify all potential mRNA targets of both miR-199a-5p and miR-199b-5p ⁶³. Using several target databases, potential mRNA targets were scored and ranked. Low-level scoring miRNA targets (50 or below) and negatively changing genes from the time series dataset were removed, leaving 21 genes that were significantly (adjusted P values < 0.05) upregulated during day 0 or day 1 following miR-199a-5p and miR-199b-5p inhibition.

Chi-squared tests

To determine if the number of overlapping genes found from differential expression analysis of 199a/b inhibition at days 0 and 1 were significant chi-squared tests were performed. Observed numbers were used to determine the estimated numbers. Expected number of differentially expressed overlapping/ non-overlapping genes = (Row Total * Column Total/ Grand Total). Chi-squared p-value calculation is performed with the following formula (Observed – Expected)^2 / Expected.

GRN development

The initial system schematic was made in BioRender (Figure 5A). For the enhanced model, we used CellDesigner, SMBL style GRNs are created to represent the biological processed of interest ⁶⁴. A literature search provided information of GRN topology. Once the whole GRN was created (Figure 5C), a more advanced GRN was created to hypothesise further pathways and important players in the system which could not be modelled as we lacked the experimental data to do so (Supplementary Figure 3).

Kinetic modelling

Base data (microarray experiments), and the validation data (Figure 2, 4) from qPCR data was converted to numbers compatible to the calibration data from the microarrays using the following formula KD/C ∗ M, KD being the mean miR-199a-5p or miR-199b-5p value, C being the mean control value and M being the mean microarray value. The initial conditions (zero time point) for the calibration and validation datasets were assumed to be the same so the model can have a single initial condition for each validated species. Thus, KD/C =1 was fixed for each zero-time point, for each species where validation data was available. For GAG levels, the control level at day 7 was treated as 100%, and the change in GAG expression during miR-199b-5p inhibition was measured in contrast to the control to calculate the percentage change.

Species selected from the modelled GRN were modelled using COPASI ⁶⁵. Calibration was performed using parameter estimation via the Particle swarm algorithm. Inhibition of miR-199a-5p and miR-199b-5p were simulated using events, where the miRNA in question was reduced by 90-95%, until day 7 in the initial model and day 4.5 in the enhanced model. Parameters were altered using sliders to make the model perform miR-199a-5p or miR-199b-5p inhibition behaviour. Mean squared error (MSE) was calculated between actual and simulated data where possible, . Model formulas and parameters are in the Supplementary materials. Data from COPASI was imported into R for plotting. All model parameters and equations have been recorded in the Supplementary Materials and the model has been uploaded onto the Biomodels public repository ⁶⁶.

Data availability

All scripts will be found in the following github directory MIR199ab5p-Chondrogenesis-Modelling-Paper/ at main · Krutik6/MIR199ab5p-Chondrogenesis-Modelling-Paper · GitHub. The enhanced kinetic model can be found in the following biomodels repository MODEL2305010001. RNAseq data can be made available on request. Experimental data in excel files can be found in the supplementary materials.

Funding

K.P., I.M.C., M.J.B. and D.Y. were supported by the Dunhill Medical Trust [R476/0516]. D.P.S. was supported by Novo Nordisk Fonden Challenge Programme: Harnessing the Power of Big Data to Address the Societal Challenge of Aging [NNF17OC0027812]. C.P. and D.Y. were supported by the MRC and Versus Arthritis as part of the Medical Research Council Versus Arthritis Centre for Integrated Research into Musculoskeletal Ageing (CIMA) [MR/P020941/1 and MR/R502182/1]. J.S. received funding from Versus Arthritis [22043].

Article and author information

Author details

Krutik Patel

Campus for Ageing and Vitality, Biosciences Institute, Newcastle University, Newcastle upon-Tyne, NE4 5PL, UK

Contribution: Bioinformatics analysis, Systems Biology analysis, Modelling, Methodology Visualization, Writing – original draft, Writing – review and editing.

Contributed equally with: Matt J Barter.

Completing interested: None to declare.

Matt J Barter

Regenerative Medicine, Stem Cells, Transplantation, Biosciences Institute, Newcastle University, Newcastle, upon, UK Tyne, NE1 4EP

Contribution: Experimental design, Experimental work, Data Processing, Writing – original draft, Writing – review and editing.

Contributed equally with: Krutik Patel.

Completing interested: None to declare.

Jamie Soul

Regenerative Medicine, Stem Cells, Transplantation, Biosciences Institute, Newcastle University, Newcastle, upon, UK Tyne, NE1 4EP

Contribution: Bioinformatic Support, Writing – original draft, Writing – review and editing.

Completing interested: None to declare.

Peter Clark

Campus for Ageing and Vitality, Biosciences Institute, Newcastle University, Newcastle upon-Tyne, NE4 5PL, UK

Contribution: Modelling Support, Writing – original draft, Writing – review and editing.

Completing interested: None to declare.

Carole Proctor

Campus for Ageing and Vitality, Biosciences Institute, Newcastle University, Newcastle upon-Tyne, NE4 5PL, UK

Contribution: Supervision, Writing – original draft, Writing – review and editing.

Completing interested: None to declare.

Ian M Clark

School of Biological Sciences, University of East Anglia, Norwich, NR4 7TJ, UK.

Contribution: Supervision, Funding, Writing – original draft, Writing – review and editing.

Completing interested: None to declare.

David A Young

Regenerative Medicine, Stem Cells, Transplantation, Biosciences Institute, Newcastle University, Newcastle, upon, UK Tyne, NE1 4EP

Contribution: Project Leadership, Project Co-ordination, Supervision, Funding, Writing – original draft, Writing – review and editing.

Completing interested: None to declare.

Daryl P Shanley

Campus for Ageing and Vitality, Biosciences Institute, Newcastle University, Newcastle upon-Tyne, NE4 5PL, UK

Contribution: Project Leadership, Project Co-ordination, Supervision, Funding, Writing – original draft, Writing – review and editing.

Completing interested: None to declare.

Conflicts of interest

None.