Introduction

Differentiation of mesenchymal stem cells (MSCs) into chondrocytes known as chondrogenesis is an important process during development, as it is a pre-requisite for skeletogenesis ¹. Chondrocytes perform the vital roles of generating and maintaining cartilage during embryogenesis, but are also found in many adult tissues, including articular cartilage ². The master transcription factor responsible for triggering chondrogenesis is SOX9 and during this process chondrocytes secrete anabolic proteins such as Collagen Type II and aggrecan, which are respectively encoded by their protein sub-units COL2A1 and ACAN, and these fibres constitute a significant functional portion of cartilaginous extracellular matrix (ECM) ^{2, 3}. Over time and influenced by a range of factors such as age, epigenetics and mechanical stress, chondrocytes express increasing levels of catabolic proteins such as MMP13 and ADAMTS5 which respectively degrade COL2A1 and ACAN ³. During skeletogenesis, endochondral ossification is triggered, and chondrocytes will naturally become hypertrophic and differentiate into osteoblasts for healthy bone growth. However, when chondrocytes within articular cartilage become hypertrophic, this leads to detrimental effects, as the cartilage maintenance is greatly diminished, leading to stiffer, less resilient, and joints less able to bear force. Ultimately the chronic loss of cartilage contributes to the extremely debilitating condition of osteoarthritis (OA), which has few treatment options, such as pain relief drugs or surgery, and no cure ^{4, 5}. Research into the regulation of chondrogenesis could identify regulatory processes and potential drug targets for OA for pain relief or to delay the progress of hypertrophy. Numerous studies have identified significantly altered microRNAs (miRNAs) during chondrogenesis, which includes investigation into the two important chondrocyte selective miRNAs, miR-140 and miR-455 ⁶. miRNAs are short non-coding RNAs - roughly 19-22nt long, that negatively regulate target mRNAs ^{7, 8}. 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 ^{9, 10}. miRNA-mRNA interactions are complex because a single mRNA can be targeted by many miRNAs, and a single miRNA can target many mRNAs. The complexity of miRNA-mRNA interactions led us to perform a combined bioinformatic, experimental and systems biology approach to better understand the relationship between microRNAs 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 of fill experimental gaps and predict novel interactions such as the indirect influence miR-199a/b-5p have over miR-140-5p.

Results

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

A 14-day chondrogenesis time series dataset we generated previously was re-analysed with the TimiRGeN R/ Bioconductor package - a novel tool we developed to analyse longitudinal miRNA-mRNA expression datasets ^{6, 11}. Prior to using TimiRGeN, the transcriptomic data underwent timepoint based differential expression analysis, using the 0 time-point 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 (Table 1). Eight signalling pathways were found to 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 ^12–14. 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.

Table 1.

hsa-miR-140-5p was the most positively changing miRNA in the following seven pathways: Adipogenesis, Clear cell renal cell carcinoma pathways, Epidermal growth factor/Epidermal growth factor receptor, Endochondral ossification, Endochondral ossification with Skeletal Dysplasia, Gastrin signalling pathway and Vascular endothelial growth factor-A/Vascular endothelial growth factor receptor 2 signalling pathway, while hsa-miR-199b-5p was the second most positively changing gene in the following six pathways: Clear cell renal cell carcinoma pathways, Epidermal growth factor/Epidermal growth factor receptor, Endochondral ossification, Endochondral ossification with Skeletal Dysplasia, Gastrin signalling pathway and Vascular endothelial growth factor-A/Vascular endothelial growth factor receptor 2 signalling pathway. Other genes such as hsa-miR-455-5p and BMP2 were also of interest, but we focused on hsa-miR-199b-5p since it was the comparatively under-researched miRNA within the context of chondrogenesis. Furthermore, our previous work identified miR-199b-5p and its paralogue miR-199a-5p as significantly upregulated during chondrogenesis, therefore we chose to focus on both paralogues of miR-199-5p 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-199a-5p overexpression on chondrogenic biomarkers ACAN, COL2A1 and SOX9 by qPCR experiments and GAG (glycosaminoglycan) levels biochemically. Significant upregulation was seen in ACAN, COL2A1 and GAG levels (Figure 2A). In contrast, when miR-199b-5p expression was inhibited ACAN, COL2A1 and SOX9 were significantly downregulated at multiple time points, especially at the earlier times (days 0-3). GAG levels were also significantly decreased, by ∼60%, at day 7 (Figure 2B). Our results indicated both miR-199a-5p and miR-199b-5p contributed towards positive regulatory roles during chondrogenesis.

Figure 2.

FZD6, ITGA3 and CAV1 were the most significantly upregulated miR-199a/b-5p targets in MSCs during chondrogenesis.

Multiple mRNA targets of miR-199a/b-5p are likely regulated during chondrogenesis. To elucidate which genes are most affected by miR-199 knock-down 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). From figure 2, we knew the chondrogenic markers were most negatively affected in the early stages of chondrogenesis, so 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 ¹⁵. 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 differentially expressed genes (adjusted P values < 0.05). Similar comparisons at day 1 of chondrogenesis revealed inhibition of miR-199a-5p or miR-199b-5p respectively resulted in 674 and 817 differentially expressed genes. 25 and 341 genes intersected between day 0 and day 1 chondrogenesis inhibition experiments, indicating that both microRNAs may share functional repertoire of targets. Contrasting the day 1 undifferentiated control samples over the day 0 undifferentiated control samples lead to 4391 upregulated and 4168 downregulated significantly differentially expressed genes (< 0.05 adjusted p value). Positively changing genes included well studied pro-chondrogenic markers such as COL2A1 (log2fc =11.6), ACAN (log2fc = 9.18) and SOX9 (log2fc = 3.37). The hundred most significantly differentially expressed genes (Supplementary Table 3) were analysed with enrichr and revealed the most likely cell type to be chondrocytes according to the Panglao database. Gene Ontology (GO) term analysis confirmed that the cells were differentiating towards chondrocytes with terms such as extracellular matrix organization and collagen-containing extracellular matrix significantly (adjusted P values < 0.05) enriched. The miR-199a-5p and miR-199b-5p knockdown experiments showed that well known chondrogenesis genes such as COL2A1, ACAN and SOX9 were significantly downregulated in day one samples, when contrasted to day 0 samples, which strongly indicates the progression to chondrogenesis has been hindered by miR-199a-5p and miR-199b-5p knockdown.

Gene Ontology (GO) term analysis was conducted to identify biological processes linked to miR-199a/b-5p inhibition (Figure 3A). Interestingly several gene sets 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. The positively enriched and significantly differentially expressed genes which intersected from the knockdown experiments were further analysed using the miRNAtap R package to predict miR-199a-5p and miR-199b-5p targets¹⁶. miRNAtap scored each predicted miRNA-mRNA interaction recorded for miR-199a-5p and miR-199b-5p from five popular miRNA prediction tools (e.g. TargetScans) to rank the most likely targets. We found twenty-one likely mRNA targets were predicted to be targeted by both miRNAs: ABHD17C, ATP13A2, CAV1, CTSL, DDR1, FZD6, GIT1, HIF1A, HK2, HSPA5, ITGA3, M6PR, MYH9, NECTIN2, NINL, PDE4D, PLXND1, PXN, SLC35A3, SLC35D1 and VPS26A. Of the 21 targets, FZD6 and NINL were the only predicted target genes found to be significantly differentially expressed in the MSCs pre-differentiation (day 0) following miR-199a-5p or miR-199b-5p inhibition.

Figure 3.

To examine the expression patterns of the 21 genes found through our bioinformatic analysis (Figure 3C), the genes were also explored in our microarray study (Supplementary Table 2) ⁶. In addition, we used SkeletalVis to contrast the behaviour of the 21 genes in other MSC derived chondrogenesis studies (Supplementary Table 3) ^17–19. 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. Interestingly, FZD6 and NINL were the only genes significantly upregulated during miR-199a/b-5p inhibition at day 0, but NINL was not further explored as they were not significantly upregulated at day 1 following miR-199a/b-5p inhibition. Size of the points reflect the mean read count of the gene, the colour represents the adjusted p.value and the shape represent if the point is from day or day 1.

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

To test the effects miR-199a-5p and miR-199b-5p had on the predicted targets that were 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 wanted to evaluate if knockdown of one of the miRNAs would affect expression of the other paralogue. To check this, we tested if knockdown of one miRNA would lead to similar or different effects to knockdown 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 knockdowns by measuring the expression of the two most significantly enriched predicted miR-199a/b-5p targets from our bioinformatic analysis – FZD6 and ITGA3 (Figure 4A1). We saw significant enrichment of FZD6 and ITGA3 after miR-199a-5p knockdown and miR-199a/b-5p knockdown, however we generally saw no statistical difference between knockdown of miR-199a-5p knockdown and knockdown of miR-199a/b-5p. 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 enriched post miR-199b-5p knockdown (Figure 4A2). Unlike FZD6 and ITGA3, CAV1 has previously been established to be a direct binding partner of miR-199a-5p ²⁰. To identify if FZD6 and ITGA3 were also direct miRNA-mRNA interactions we cloned the 3’ UTR regions of FZD6 and ITGA3 directly downstream of the firefly luciferase gene in the vector pGL3-basic which also contains a control renilla luciferase gene for transfection efficiency normalisation (Figure 4B1). These plasmids were transiently transfected into the chondrocytic cell line SW1353 along with control or miR-199a-5p mimics. Both FZD6 and ITGA3 3’UTRs reduced the firefly/renilla ratio levels when compared to the empty control vector, suggesting these contain potentially repressive elements. Indeed, the co-introduction of an excess of miR-199a-5p into the cells further repressed the 3’UTR containing vectors firefly/renilla ratio but not that of the control vector, suggesting both genes are direct miR199a-5p targets. Furthermore, pmiRGLO experiments showed after miR-199a-5p overexpression, FZD6 and ITGA3 levels increased significantly. Mutations added to the 3’-UTR lead to lower negative effects on the genes, and even had the reverse effect on mut-ITGA3.

Figure 4.

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 knockdown experiments. Finally, we used events to simulate knockdowns of miR-199a-5p or miR-199b-5p (Figure 5B). The experiments shown in figures 2A, 2B and 4A were used to parameterise the model, and since the experiments were performed in a staggered manner, we can use to model to make predictions to fill experimental gaps. For example, using this model we predicted the dynamics of the chondrogenesis biomarkers and CAV1 after miR-199a-5p knockdowns and the dynamics of FZD6 and ITGA3 after miR-199b-5p knockdown. 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). Most species, have a MSE of lower than 3 during the different conditions (∼84%) (3/14) indicating many of the experimental dynamics were captured by this model.

Figure 5.

We added further detail to enhance our model based on data from literature. This included the addition of Tumour 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 major 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, the enhanced chondrogenesis model enabled us to address some weaknesses in our initial model. Firstly, by simulating the miR-199a-5p and miR-199b-5p knockdown 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 knockdown were flat. Furthermore, our CAV1 dynamics were also flat during miR-199b-5p knockdown. Figure 2B shows, from miR-199b-5p knockdown, 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 knockdown experiments, we simulated the knockdown 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 may have been alluded to in our RNAseq experiment, but not further explored such as MYH9, 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 knockdown, 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. This meant the model could also predict the dynamics of different omics of SOX9 – but we lacked the data to provide any evidence for our predictions. Finally, our initial model used simplistic functions, with few parameters and we found using functions with additional parameters improved our model dynamics (Supplementary Materials) (Figure 5C, Figure 5D).

The enhanced model improved dynamics for several model objects e.g. GAG, CAV1, SOX9, ACAN, COL2A1. Generally, our MSE values were better in our enhanced model for the chondrogenesis biomarkers and GAG levels, although overall the MSE values were not improved for the miR-199a/b-5p targets. We found that now 80% (3/15) of the MSE values were below 3, indicating our enhanced model provided a good fit to most of the experimental data. 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 and important genes such as the SMAD2/3 signalling pathways and 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. To assist with future hypothesis generation, these pathways which were beyond the scope of this project are shown in the GRN (Supplementary Figure 2).

Discussion

In this paper we performed experimental and computational analyses to better understand 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 ^{21, 22}. 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 generation ⁶. 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 the TGFB triggering phosphorylation of CAV1 via SRC-kinase, and the TGFB pathway which is of great importance in chondrogenesis ²⁰. In contrast to CAV1, FZD6 and ITGA3 have been less well studied as miR-199a/b-5p targets, and though, miR-199a-5p-FZD6 has been predicted previously, our luciferase assays show for the first time that FZD6 has been validated as a target of miR-199a-5p-FZD ^{22, 27}. 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 5D) 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) ^29–31. 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 ^{6, 32–34}. Precise mechanisms of how FZD6 and ITGA3 regulated chondrogenesis are unclear, with potentially implicated pathways were 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 ^{35, 36}. 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 ^{37, 38}. 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, and also 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 ^{30, 40–42}. It is likely that other mRNA targets of miR-199a/b-5p also contributed 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 2), which showcases how the models we created could be enhanced with such data.

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 and provide 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 ^{43, 44}. 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 time point specific gene filtration was set as <0.05 and adjusted P values were used for filtration. Time point 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 entrezgene IDs using BiomaRt ^{47, 48}. 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.

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

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 knockdown (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 ound 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.

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 2).

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 knockdown 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, 1/n. ∑ⁿ_i=1 (yi − yi)². Model formulas and parameters are in the Supplementary materials. Data from COPASI was imported into R for plotting.

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 will be in excel files which will be found stored with 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].

Acknowledgements

Contributions

All authors provided substantial contributions to conception or design of the study: J.S., D.A.Y. Substantial contributions to drafting the manuscript: K.P, M.J.B., D.S., and D.A.Y. Substantial contributions to data acquisition: M.J.B. All authors contributed to data analysis or interpretation, but this was predominately driven by K.P., J.S., and M.J.B. All authors contributed to revising the manuscript critically for important intellectual content and approved the final manuscript.

Conflicts of interest

None.