Changes of Gene Expression Profile during Differentiation of Equine Bone Marrow-Derived MSCs towards Adipocytes and Chondrocytes*

Bone marrow (BM) was obtained from slaughtered horses located in Lesser Poland Voivodeship, according to the protocols described by Opiela et al. (2013, 2016). The BM samples were obtained from 9 horses (6 females and 3 males) representing Polish coldblood breeds. The horses were between 5 and 19 months of age (average 10.5 months), with an average weight of 364 kg. BM biopsy was performed with a dedicated needle (Medical Device Technologies, Inc.) from the ilium and ischium of the horse. The injection site was disinfected with disinfectant before biopsy. For BM biopsy, the injection needle was inserted through the muscle sheath into the ilium to a depth of approximately 1 cm. A syringe containing 1 mL of heparin was then attached to the injection needle, and bone marrow aspiration was initiated. After aspiration, the BM sample was transferred from the syringe to a 50 mL Falcon tube with 3 mL of heparin.

For the isolation of mesenchymal stem cells (MSCs), bone marrow samples were suspended in phosphate-buffered saline (PBS) (Biomed, Lublin) at a 1:1 ratio and layered over Ficoll-Paque (Stem Cells Technologies, USA) using a density gradient centrifugation (30 minutes at 500 × g). The mononuclear cell fraction was collected and washed twice in PBS. Finally, the cells were centrifuged for an additional 5 minutes at 500 × g. The pelleted cells were suspended in culture medium comprising low-glucose Dulbecco's modified Eagle's medium (DMEM) (Merck, Germany) supplemented with 10% foetal bovine serum (FBS; Merck, Germany), 1% antibiotic antimycotic solution (AAS; Merck, Germany) and 1% GlutaMAX (100X) (Gibco, USA). The cells were cultivated in 17 ml of culture medium in T-75 cm² tissue culture flasks. The culture medium was changed every three days. The adherent cells were allowed to form colonies until they reached 90% confluence (Opiela et al., 2013, 2016). The MSCs were maintained under standard conditions at 37°C in an atmosphere enriched with 5% CO₂ and 100% humidity under normoxia (21% O₂) in a HERAcell^® CO₂ incubator (Thermo Scientific) (Opiela et al., 2013, 2016).

After reaching approximately 90% confluence, the cells were rinsed with PBS and then digested with a 0.25% trypsin/EDTA solution (Merck, Germany). Then, after the trypsin was inactivated with FBS solution, the cells were collected and centrifuged for 5 min at 500 × g. The cell pellets were suspended in 0.5 ml of 10% dimethyl sulfoxide (DMSO, Merck, Germany) and transferred to cryovials, which were placed in Mr Frosty™ Freezing Container (Nalgene, Thermo Scientific) filled with 100% isopropyl alcohol. Next, the container was placed in a low-temperature freezer (−80°C) for 24 hours. The mentioned container model guarantees repeatable freezing of cells at a rate of −1°C/minute. Centrifugation was carried out at room temperature at every stage of the procedure. Finally, the cells were stored in liquid nitrogen (−196°C) (Fletcher et al., 2008).

The evaluation of specific gene expression markers was conducted using RT² Profiler PCR Arrays (Qiagen, Hilden, Germany). Briefly, RNA was isolated using the miRNeasy^® Mini Kit (Qiagen) following the manufacturer's instructions. Five hundred nanograms of total RNA was reverse transcribed using the RT² First Strand Kit (Qiagen). Next, the cDNA was combined with an appropriate amount of RT² SYBR Green Master Mix (Qiagen) and distributed into the wells of the RT² Profiler PCR Array (Qiagen). The relative expression was determined using the ΔΔCT method and five reference genes (ACTB, GAPDH, B2M, HPRT1 and RPL13A), which were included in the arrays by the manufacturer, on the dedicated web portal www.SABiosciences.com/pcrarraydataanalysis.php. The results are presented in Supplementary Data S1.

Flow cytometry-based analysis was performed to assess the phenotype of the MSCs. The mesenchymal stem cells were examined for specific positive (CD44 (Bio-Rad cat. MCA1082PE), CD29 (BioLegend cat. 303015), CD90 (Santa Cruz Biotechnology cat. SC6071)) and negative (CD34 (BD Biosciences cat. 555621)) and CD13 (Agrobio cat. ARG23408)) protein markers were analysed via flow cytometry.

MSC samples were washed with CaCl₂− and MgCl₂− free Dulbecco's modified phosphate-buffered saline (DPBS; Sigma-Aldrich, Merck, Poznań, Poland) and separated from the monolayer using 0.25% trypsin-ethylenediaminetetraacetic acid (trypsin/EDTA) (Sigma-Aldrich, Merck, Poznań, Poland). Dissociated cells were centrifuged and washed with DPBS supplemented with 2% BSA and incubated for 30 minutes at 37°C in the dark with antigen-specific antibodies (antibody concentration: 1 µg/~1 million cells for CD29; 2 µg/~1 million cells for other markers) conjugated with the appropriate dyes: CD 44-PE, CD29-AlexaFluor 488, CD90-PE-Cy5.5, CD34-FITC and CD13-PE (dosage according to the manufacturer's instructions). The cells were then washed three times in DPBS and resuspended in 500 mL of DPBS. The fluorescence of the antibody-labelled MSC samples was measured using a CytoFlex flow cytometer (2015 Beckman Coulter, CA, USA). The dyes were excited by blue light (488 nm) induced by a 15 mW argon laser. The emission wavelengths in the green fluorescence spectrum ranging from 480 to 530 nm were quantified in the FITC channel. Data collection and analysis were performed using CytExpert 2.4 software. To exclude contaminants, cells were acquired using forwards scatter (FSC) and side scatter (SSC), and an MSC cell population gate was defined. Positive and negative populations were determined relative to the unstained sample. Images of the stained and unstained samples were overlaid on a graph (FITC green fluorescence/SSC side scatter), and then the portion comprising the unstained sample was subtracted to obtain the remaining positive population. The results are presented in Supplementary Data S1.

After stable MSCs were established, the cells were subjected to differentiation towards adipocytes in culture medium comprised of low-glucose DMEM supplemented with 15% FBS, 1% GlutaMAX (100X), 1% AAS and the appropriate supplements.

A total of 33 µM biotin, 0.1 µM dexamethasone, 1 µM insulin, 200 µM indomethacin, 17 µM D-panthotenic acid, 10 µg/ml transferrin human, 0.2 nM triiodothyronine, sodium salt, and 250 µM 3-isobutyl-1-methylxantin were used as supplements for the culture medium during the 3-day adipogenic differentiation of MSCs. For the next 18 days, the cells were cultured in adipogenic differentiation maintenance medium consisting of 33 µM biotin, 0.1 µM dexamethasone, 1 µM insulin, 200 µM indomethacin, 17 µM D-panthotenic acid, 10 µg/ml transferrin human, 0.2 nM triiodothyronine, and sodium salt according to the protocols of Pittenger et al. (1999) and Quevedo et al. (2009). All reagents used in this study were procured from Merck, Germany. The ADM primary cells were stored in liquid nitrogen.

After stable MSCs were established, the cells were subjected to differentiation towards chondrocytes in culture medium comprised of high-glucose DMEM supplemented with 15% FBS, 1% GlutaMAX (100X), 1% AAS and the appropriate supplements according to Branly et al. (2017). The chondrogenic differentiation medium supplements used were as follows: 100 nM dexamethasone (Merck, Germany), 50 μg/ml L-ascorbic acid-2-phosphate (Merck, Germany), 40 μg/ml L-proline (Merck, Germany), 1 nM sodium pyruvate (Merck, Germany), 1% insulin transferrin selenium (Gibco, USA) and 10 ng/ml transforming growth factor β1 (TGF-β1; Miltenyi Biotec) (Branly et al., 2017). The cells were cultured for 21 days.

The culture of MSCs as a negative control was performed in low-glucose DMEM supplemented with 10% FBS, 1% GlutaMAX (100X) and 1% AAS (100×) (Marion and Mao, 2006). The CDM primary cells were stored in liquid nitrogen.

The correct differentiation of MSCs was confirmed via the use of specific dyes. Before staining, differentiated MSCs and negative control MSCs were fixed in 4% paraformaldehyde (Merck, Germany) for 1 hour and then washed in PBS. Adipocytes were revealed by the intracellular accumulation of lipid droplets using Oil red O (Merck, Germany) staining and identified by their bright red colour. Briefly, the cells were fixed, washed and stained with 0.3% Oil red O in isopropanol and incubated for 10 min at room temperature (Branly et al., 2017). Chondrocytes were incubated with a 1% Alcian blue solution in 3% acetic acid, pH 2.5 (Merck, Germany), for 30 min. The sulphated proteoglycans were stained blue. Protocols were performed according to Rentsch et al. (2010). All of the abovementioned staining methods were applied to MSCs, which constituted a negative control for differentiated cells, and no characteristic staining was observed. The results are presented in Supplementary Data S1.

RNA was isolated from differentiated lines (n=3 from each ADM and CDM line) and the control MSC line (n=3) via a PureLink™ RNA Mini Kit (Ambion) with DNAse digestion between washes. Its quality and quantity were measured on a NanoDrop2000 and a TapeStation 2200 system (RNA tapes, Agilent). Acceptable samples had an RNA integrity number (RIN) over 9.5. In accordance with the manufacturer's protocol, cDNA library cells were prepared from 300 ng aliquots of purified RNA samples via the VAHTS Universal V8 RNA-seq Library Prep Kit for Illumina. The final cDNA library concentration was 10 nM. The quality and quantity of the cDNA libraries were assessed via D1000 tapes on a TapeStation 2200 and a Quibit. The cDNA was sequenced over 150 PE cycles on a HiSeq 3000 System (Illumina, San Diego, CA, USA) at Admera Health Biopharma Services.

The raw data were initially subjected to quality assessment using FastQC software (Andrews, 2010) and subsequently underwent adapter removal and filtering for low-quality reads (those with a Phred quality score below 20 and a read length of less than 36) via Flexbar software (Dodt et al., 2012). Following this preprocessing step, the filtered reads were aligned to the Equus caballus (EquCab3) genome using STAR software (Dobin et al., 2013). After alignment, the mapped reads were annotated and quantified against specific gene thresholds utilizing the Ensembl gtf file version 100, facilitated by HTSeq-count software (Anders et al., 2015). Differential expression analysis was then conducted utilizing Deseq2 software (Love et al., 2014).

The functional annotation clustering of the transcriptomic results was performed using the functional annotation tool DAVID, and protein–protein interactions (PPI) were assessed using the String database.

The qPCR method and Pearson correlation were applied to validate the RNA-seq results. For validation, 9 DEGs – BMP7, BMPR1B, CLDN1, IHH, SOX9, PHOSPHO1, LEP, FABP4, and GERP1 – were used to confirm the RNA-seq results (Supplementary Table S1). The B2M and RPL13 genes were used as endogenous controls. cDNA was prepared from each RNA sample (250 ng) using a high-capacity RNA-to-cDNA™ Kit (Applied Bio-systems, Thermo Fisher Scientific) according to the manufacturer's protocol. The real-time PCRs for each sample were performed in triplicate with the use of the sensitive CiTi Mix EvaGreen^® (A&A Biotechnology) according to the protocol, with an annealing temperature of 58°C. The Pearson correlation coefficients and P value (<0.05) were generated using the SAS Enterprise 7.1 tool.

The preparation of RNA-Seq libraries was conducted for a total of 9 samples, which were stratified into three study groups: MSC, CDM, and ADM. Overall, the sequencing procedure generated raw read counts ranging from 33.5 million to 73.2 million per sample, with an average of 57.4 million. Notably, approximately 71.8% of the filtered reads exhibited unique alignment to the reference genome, and of these uniquely aligned reads, 72.4% were successfully assigned to gene thresholds sourced from the Ensembl annotation file (version 100). Detailed information regarding sample PCA clustering (Supplementary Data S2) and sequencing read and mapping statistics can be found in Supplementary Table S2.

The differential expression analysis identified 2325 DEGs between the MSC and ADM groups, whereas 803 DEGs were detected between the MCS and CDM primary cells. According to the comparison of both DEG sets, the expression of 320 genes was significantly modified despite the direction of differentiation (Figure 1).

Figure 1.

Venn diagram of DEGs identified between the analysed primary MSC, ADM and CDM cells

The top 30 DEGs were selected to show the molecular network most modified during differentiation towards equine adipocyte primary cells. The top 30 genes had log₂-fold changes between 13.09 and 7.49, which corresponded to differences in the expression profiles between the groups from 8723.9 to 180.2 FC. The selected genes are involved in the regulation of the cellular macromolecule biosynthetic process, fat cell differentiation, regulation of the lipid biosynthetic process and positive regulation of the developmental process (Figure 2). The identified top DEGs represent genes related to adipogenesis: MEOX1 (−9.54) and CYYR1 (−8.06) were downregulated and CYP24A1 (13.09); LEP (9.37); TBX5 (11.66); FABP4 (11.44); and GPER1 (8.17) were upregulated. Moreover, GPER1, which is responsible for increasing the number of adipocytes and lipid accumulation, was found to be involved in the greatest number of GO terms (Figure 3). In turn, two of the top DEGs, FABP4 and MMP1, belong to the PPAR signalling pathway and are also considered critical triggers of adipogenesis.

Figure 3.

The involvement of the identified top DEGs (MSC vs. ADM cells) with significant GO terms (FC values according to the RNA-seq results) (R package)

Interestingly, 21 of the 30 top DEGs were expressed in only one group (ADM or MSC) (Figure 4). The complete turning on and turning off of gene expression may indicate a great impact on cell gene expression changes during differentiation processes. The gene expression pattern confirmed that the observed changes occurred via both up- and downregulation mechanisms.

Figure 4.

Heatmap of the DEGs between MSC and ADM primary cells identified as the top DEGs and showing their expression in only one group. A unit variance scale was used with no logarithmic transformation. (R software, ClustVis package)

The top 30 DEGs were the genes with the most altered expression during differentiation towards equine chondrocyte primary cells (log₂fold change between 10.04 and 5.53, which corresponds to differences in the expression profile between groups from 1359.79 to 46.29 FC). The highest log₂FC values between the MCS and CDM groups were observed for the upregulated genes CLDN16 (10.40), HPSE2 (8.94) and BMP7 (8.28) and the downregulated genes CYP4F119 (7.50) and NTRK2 (6.72).

The results revealed 27 DEGs whose expression was detected in only one group (MSC or CDM). During differentiation towards CDMs, the expression of 12 genes was activated, and the expression of 15 genes was completely silenced (Figure 6). Moreover, 19 of these DEGs were detected as the most differentiated (CLDN16, HPSE2, BMP7, CYP4F119, RAPGEF4, APOB, NTRK2, IHH, CDKL1, GPR87, PROK2, FAM110D, BMPR1B, SCNN1B, CCR7, CCDC182, CYRIA, TMEM74B and CCDC69), among which 4 DEGs are involved in the regulation of cell differentiation (IHH, APOB, BMPR1B and NTRK2) and specifically in chondrocyte differentiation (IHH and BMPR1B).

Figure 6.

Heatmap of the DEGs between MSCs and CDM primary cells identified as the top DEGs and showing their expression in only one group. A unit variance scale was used with no logarithmic transformation. (R software, ClustVis package)

The comparison of the sets of DEGs detected for the two types of cell differentiation allowed for us to identify 320 genes that were differentially expressed between the two treatments (Supplementary Data S3). The GO analysis revealed several terms enriched despite the type of treatment, related to the regulation of signalling (GO:0023051; 61 DEGs), the regulation of cell communication (GO:0010646; 60 DEGs), the regulation of signal transduction (GO:0009966; 56 DEGs), the regulation of cell population proliferation (GO:0042127; 33 DEGs), the response to growth factors (GO:0070848; 22 DEGs) and the cell surface receptor signalling pathway (GO:0007166; 48). Moreover, the macromolecule localization (GO:0033036; 51) and protein phosphorylation (GO:0006468; 31) GO terms were identified. The results confirmed that some universal changes occurred in equine cells during differentiation and were related to cell growth and development, communication, and intra-cellular organization.

Compared with the relative quantity from real-time PCR and normalized counts from the NGS data, the qPCR results were strongly correlated with the RNA-seq data. High correlation coefficients were obtained for all the analysed DEGs, but they were not significant for two of them – LEP and GERP1 (Table 3). For the 7 validated DEGs, the correlation coefficients ranged from 0.95 to 0.99, which confirmed the accuracy of the whole-transcriptome approach.

The risks and challenges resulting from MSC differentiation and subsequent clinical therapy require extensive studies. The implementation of the application of MSC-based therapies and the safety of primary cells derived from MSCs indicate the need for molecular research, such as whole-transcriptome studies. However, investigations of the gene expression profiles of equine cells derived from MSCs under in vitro conditions are still limited. The present investigation was conducted with the aim of expanding the knowledge in the field of chondrocytes and adipocyte cells obtained from equine bone marrow-derived MSCs.

Comparative transcriptome analysis revealed greater gene expression changes during differentiation towards adipocytes (2325 DEGs) than during differentiation into chondrocytes (803 DEGs). We detected 320 genes whose expression was modified during MSC differentiation. Overrepresentation of the following basic cell functions was identified: the regulation of signalling, cell communication, signalling, and proliferation. The identified GOs involved a panel of genes responsible for pluripotency, proliferation, and differentiation – SOX9, PAX6, BMP4, GPC3, LIMS2, ROBO1, and VASH1. Previous research in a mouse model demonstrated the critical function of SOX9 during chondrocyte differentiation (Murakami et al., 2000) via alterations in gene expression in chondrocyte cells (Leung et al., 2011). Moreover, the BMP4-SOX9 regulatory axis has been proposed as the main pathway involved in ASC differentiation into chondrocytes (Jang et al., 2016). In rats, SOX9 regulates the survival of cells and the adipogenic differentiation of multipotent mesenchymal stem cells (Stöckl et al., 2013). In turn, the latest report showed that overexpression of the PAX6 gene is a molecular factor enabling the differentiation of bone marrow MSCs into limbal epithelial stem cells (Gao et al., 2023). Our results confirmed that SOX9 together with PAX6, BMP4, and an identified set of DEGs also play essential roles in equine MSC differentiation, probably on a similar basis as in other species.

The identification of gene expression changes that occur following the differentiation of equine MSCs towards adipocytes confirmed the accuracy of the cell specialization process. MSC-derived adipocytes exhibit unique genetic patterns, which affirms that we generated fat cells. After the MSC differentiation treatments, the gene encoding CCAAT enhancer binding protein alpha (CEBPA) was significantly upregulated in the adipocyte primary cells. The role of the CEBPA gene has been associated with adipocyte differentiation in numerous studies. CEBPA induces adipogenesis via the activation of several adipocyte-specific genes (Li et al., 2024), including PPARγ (Rosen et al., 2002; Madsen et al., 2014), and, as a result, promotes the differentiation of preadipocytes into adipocytes. The transcriptome analysis results revealed increased expression, not only of CEBPA but also of the PPARG gene, the critical factor in adipogenesis regulation (Zhang et al., 2023). Importantly, a study performed in humans revealed CEBPA as one of the adipogenesis markers that determines adipose tissue remodelling (Matulewicz et al., 2017). The authors proposed several genes (CEBPA, ADIPOQ, IRS1, IRS2, SLC2A4, and MMP9) as markers of adipogenesis but not the inflammation process. Our results revealed significant upregulation of three genes, CEBPA, ADIPOQ and SLC2A4, in ADMs, confirming the specificity of differentiation. The deregulation of MMP9 was not detected, but other genes from the MMP family – MMP1 (identified as the top differentially expressed gene log₂FC 9.43), MMP12, MMP16, and MMP21 – were found to be significantly modified during differentiation into equine ADMs.

Our results revealed a significant increase in LEP and FABP4 expression in equine adipocytes derived from MSCs. Notably, both genes were expressed only in the adipocyte primary cells, which means that the use of targeted differentiation procedures induced their expression in ADM cells de novo. Lin et al. (2023) presented genes essential for adipogenesis modifications that occur in humans, which are partially consistent with our results. Similar to the findings of Lin et al. (2023), we detected significant modifications of the PPARG, FABP4, DGAT2, ADIPOQ, LEP, CEBPA and MEOX1 genes in the obtained adipocyte primary cells. Additionally, two other upregulated genes identified in adipocytes, GPER1 and TBX5, have been investigated previously in terms of the adipocyte differentiation process. In other species, G protein-coupled oestrogen receptor 1 (GPER1) promotes early adipogenesis and cell proliferation and increases the number of adipocytes (Chappell et al., 2018; Zhao et al., 2022). Similarly, TBX5 is involved in the initial stages of adipocyte differentiation, and its methylation profile in humans determines fat deposition (Bradford et al., 2019), underscoring the importance of the epigenetic nature of biological processes relevant to differentiation. On the basis of our results and the Equus caballus GO reference background, we confirmed that both genes, GPER1 and TBX5, are involved in fat cell differentiation and the regulation of lipid biosynthetic processes in horses.

Interestingly, the most differentially expressed gene in equine adipocytes was CYP24A1, which encodes cytochrome P450 family 24 subfamily A member 1. The identified log₂FC for this gene was 13.09, and the difference in expression level between MSC and ADM primary cells was 1273.46 vs. 0. The limited literature data indicate that CYP24A1 controls and improves adipogenesis through vitamin D signalling (Nimitphong et al., 2020). Significantly increased CYP24A1 expression was demonstrated in human preadipocytes and newly differentiated adipocytes in response to D3 (Nimitphong et al., 2020). This greater CYP24A1 upregulation observed in the equine adipocyte line without the administration of vitamin D may suggest that this gene is critical for adipocyte differentiation and proliferation in this species, but it may act through different molecular mechanisms, which should be investigated in the future.

The investigation of gene expression changes that occurred in equine adipocytes during differentiation revealed several universal mechanisms related to adipogenesis, confirmed the accuracy of differentiation treatment, and expanded the knowledge necessary for the development of an in vitro cell differentiation strategy.

The presented reports focused on two directions of differentiation that allow for adipocyte and chondrocyte primary cells to be obtained. A whole-transcriptome comparison between the initial MSCs and the differentiated chondrocytes revealed the significant deregulation of 803 genes associated with the regulation of cell differentiation and apoptosis, cell adhesion, communication and developmental processes. One of the identified overrepresented biological processes was ossification, represented by 35 DEGs, including the most upregulated genes, IHH, PHOSPHO1 and SRGN, and the downregulated genes BMPR1B, IL6 and MN1. BMP signalling is crucial in the initial phases of chondrogenesis (Yoon et al., 2006). A study in an in vitro model revealed that BMPR1B together with BMPR1A are functionally redundant during early chondrogenesis and that silencing their expression is necessary for in vivo chondrocyte proliferation, viability, and maturation chondrogenesis (Yoon et al., 2005). Our results strongly support these findings through the detection of decreased BMPR1B expression in equine chondrocyte cells. In turn, we revealed that the upregulation of BMP2, BMP4 and BMP7, which are considered key factors (BMP-2,4,7), triggers the differentiation of MSCs into chondrocytes and is vital for the conversion of prehypertrophic chondrocytes to hypertrophic chondrocytes (Chen et al., 2021; Wu et al., 2024). In our study, the transcript level of BMP7 was increased 312-fold in chondrocytes originating from MSCs. These findings confirmed that the BMP family, especially BMP7, promotes chondrogenesis in horses, similar to humans (Zhou et al., 2011). Furthermore, a recent report revealed that BMP7 suppressed hypertrophic and fibrotic chondrocyte phenotypes and demonstrated advantageous clinical effects (Ripmeester et al., 2021). Together with our findings, these results indicate that BMP7 can be considered a positive genetic marker of the process of chondrocyte differentiation.

Moreover, Chen et al. (2021) reported the significant involvement of SOX9 in controlling the differentiation of bone marrow MSCs into chondrocytes, which is in accordance with our results. The downregulation of BMPR1B during the maintenance of bone mass and the transduction of BMP signalling was also confirmed by Shi et al. (2016). Additionally, the authors proposed that the ACVR1 gene is an important factor during chondrocyte differentiation. In addition to the modification of BMPR1B, the obtained data allowed for us to detect significant downregulation of ACVR2, the second receptor of activin (activin A receptor type 2A). To date, the ACVR1 gene has been associated mainly with chondrocyte differentiation, whereas ACVR2 has been proposed to be a negative regulator of bone mass via the control of osteoblast proliferation (Goh et al., 2017). The literature concerning ACVR2 function is still limited. Our findings strongly indicate the possible role of the ACVR2 gene in equine chondrocyte growth and differentiation together with BMP signalling pathways, similar to other species.

The equine chondrocytes originating from MSCs presented increased expression of the IHH gene. The Indian hedgehog protein is secreted by prehypertrophic chondrocytes and may activate chondrocyte proliferation via the C2H2-type zinc finger protein subclass of the Gli family or inhibit terminal differentiation by the IHH/PTHrP loop (Hilton et al., 2005; Jiang et al., 2020; Fan et al., 2022). Our results revealed the upregulation of the IHH gene and the downregulation of GLIS3 (a nuclear protein belonging to the GLI-like zinc finger), which may indicate that it still delayed the proliferation mechanism of chondrocytes. The obtained data indicated that the IHH signalling pathway might play a key role in equine cell remodelling towards chondrocyte differentiation.

The second most upregulated gene in equine chondrocytes from MSCs that is involved in ossification GO was PHOSPHO1, a bone-specific phosphatase (44-fold greater expression in CDMs than in MSCs). Although PHOSPHO1 is essential for bone matrix mineralization (Suchacki et al., 2020), its elevated expression was observed in growing chondrocytes. Another study indicated that PHOSPHO1 plays an important role in chondrocyte-derived matrix vesicle biogenesis (Roberts et al., 2007). Interestingly, the claudin 16 gene, CLDN16, was increased 1359-fold in chondrocyte cells. In mice, claudin 12 promotes chondrocyte differentiation and cartilage growth and development (Xing et al., 2022). In turn, claudin 11 is overexpressed in chondrocytes, which results in the activation of articular cartilage markers (Lindsey et al., 2019). In humans, according to the String dataset and Hu et al. (2016), there is an interaction between CLDN16 and PHOSPHO2. The present results indicated the possibility of a molecular interaction between CLDN16 and PHOSPHO1 in equine chondrocytes, which, owing to the level of the observed changes, may be critical for controlling chondrogenesis in horses.

While this study provides valuable insights into the transcriptomic changes accompanying the differentiation of equine bone marrow-derived MSCs into adipocytes and chondrocytes, several limitations should be acknowledged: the use of in vitro conditions that may not fully replicate the in vivo environment; analyses limited to a single time point; potential MSC population heterogeneity; and the lack of accompanying epigenetic and proteomic data to complement the transcriptomic findings.