Dynamics of DNA methylation during osteogenic differentiation of porcine synovial membrane mesenchymal stem cells from two metabolically distinct breeds

ABSTRACT

Mesenchymal stem cells (MSCs), with the ability to differentiate into osteoblasts, adipocytes, or chondrocytes, show evidence that the donor cell’s metabolic type influences the osteogenic process. Limited knowledge exists on DNA methylation changes during osteogenic differentiation and the impact of diverse donor genetic backgrounds on MSC differentiation. In this study, synovial membrane mesenchymal stem cells (SMSCs) from two pig breeds (Angeln Saddleback, AS; German Landrace, DL) with distinct metabolic phenotypes were isolated, and the methylation pattern of SMSCs during osteogenic induction was investigated. Results showed that most differentially methylated regions (DMRs) were hypomethylated in osteogenic-induced SMSC group. These DMRs were enriched with genes of different osteogenic signalling pathways at different time points including Wnt, ECM, TGFB and BMP signalling pathways. AS pigs consistently exhibited a higher number of hypermethylated DMRs than DL pigs, particularly during the peak of osteogenesis (day 21). Predicting transcription factor motifs in regions of DMRs linked to osteogenic processes and donor breeds revealed influential motifs, including KLF1, NFATC3, ZNF148, ASCL1, FOXI1, and KLF5. These findings contribute to understanding the pattern of methylation changes promoting osteogenic differentiation, emphasizing the substantial role of donor the metabolic type and epigenetic memory of different donors on SMSC differentiation.

KEYWORDS:

• DNA methylation
• Mesenchymal stem cells
• Osteogenic differentiation
• Pig breeds
• Epigenetic pattern

Introduction

Mesenchymal stem cells (MSCs) are multipotent cells with the ability to differentiate into various cell types of the mesodermal lineage, which have been widely used in regenerative medicine. They were first isolated from the bone marrow and then had been found in various other tissues, while the ability of these MSCs to differentiate into different cell types varies according to tissue source. Despite their robust osteogenic ability, MSCs derived from bone marrow are rare and always need invasive procedures. Studies have reported that synovial membrane mesenchymal stem cells (SMSCs) have better osteogenic potential compared with MSCs derived from adipose [1,2]. Therefore, SMSCs have generated great interest in the treatment of bone and joint diseases [3]. We have previously characterized the porcine SMSCs and investigated the differential gene expression patterns during in vitro osteogenic induction [4]. However, the underlying epigenetic patterns that regulate osteogenesis of SMSCs are as yet poorly understood.

Epigenetics refers to no alteration in DNA sequence itself, but to mitotically and meiotically heritable changes in gene expression. DNA methylation on the fifth position of cytosine (5mC) is one of the main epigenetic modifications which predominantly happens as paired symmetric methylation of CpG dinucleotide. DNA methyltransferases (DNMTs), including DNMT1, DNMT3a, and DNMT3b, are responsible for catalysing the methylation of CpG dinucleotides. DNMT1 functions as a maintenance methyltransferase, preserving existing methylation patterns, while DNMT3a and DNMT3b serve as de novo methyltransferases, establishing new methylation marks [5]. Many studies have suggested that DNA methylation is required for osteogenic differentiation [6,7]. It has been observed that the supplement of DNA demethylating agent leads to increased osteogenic performance of MSCs [8,9]. DNA methylation can directly regulate transcription by interfering with the recognition sites of transcription factors (TFs) [10]. Meanwhile, several key TFs, such as RUNX2, PPAR γ2, and SOX9, are thought to regulate osteogenic differentiation of MSCs by DNA methylation [11–13]. Although it is commonly perceived that DNA methylation dynamically regulates diverse biological processes [14], only a few studies have described the global DNA methylation changes in osteogenic lineage commitment [15–17]. More importantly, little attention has been paid to the temporal modulation of genome-scale DNA methylation at the single nucleotide level during osteogenic differentiation. The dynamic characterization of DNA methylation changes can contribute to a better understanding of osteogenic differentiation and bone development.

There is mounting evidence of the complex interplay between intermediary metabolism in bone and whole-body metabolism [18]. There is increasing evidence that the metabolic type of the donor cell plays an important role in the osteogenic process of the mesenchymal stem cells. There has been evidence that the osteogenic capacity of MSCs from donors with obesity and/or metabolic disorders is impaired [19–22]. Pigs are seen as an attractive option for bone and MSCs research [23,24], while recent studies explored the epigenetic differences between MSCs isolated from lean and obese pigs [25–27]. Although these studies showed how extrinsic dietary factors affect the epigenetic landscape at the DNA and histone level, it remains unclear whether the intrinsic epigenetic differences of pig breeds with distinct metabolic patterns can influence the molecular characteristics of MSCs, especially in the process of osteogenic differentiation. Traditional pig breeds, like Angeln Saddleback (AS), tend to have higher fat content compared to modern breeds, such as German Landrace (DL). Understanding the epigenetic differences of MSCs derived from donors with distinct inheritable metabolic traits in osteogenic differentiation will allow to optimization of therapeutic osteogenic applications.

Therefore, in this study, we investigated the genome-wide DNA methylation changes of pig SMSCs during in vitro osteogenic differentiation, taking into account not only different time points but also two pig breeds with different metabolic features. We also performed an integrative analysis on DNA methylation and gene expression data from our previous study [4], as well as transcription factor-binding motifs enrichment analysis, in order to identify latent methylation-driven genes and key transcription factor during in vitro osteogenesis.

Materials and Methods

Isolation, culture and osteogenic differentiation of porcine SMCSs

We used cell samples from our previous study, which were derived from fibrous synovial tissue of the knee joint of piglets of two breeds: German Landrace (DL, n = 3) and Angeln Saddleback (AS, n = 3) [4]. These SMSCs expressed stemness markers CD44, CD29, CD90, and CD105, confirmed by flow cytometry, and could differentiate into osteocytes, adipocytes, and chondrocytes [28]. The isolation methods, culture conditions, morphological aspects, and molecular characterization of synovial mesenchymal stem cells (SMSCs) were detailed in our previous study [4,28]. Our previous study demonstrated that during osteogenic differentiation, mineral deposition was first observed at day 14 and further increased until day 21 as confirmed by the osteogenic markers including ALPL, COL1A1, RUNX2 and NANOG [4]. In brief, after SMSCs at passage 3 became nearly 80% of confluence, cells for osteogenic induction (Ost) were incubated with osteogenic differentiation medium (ODM), while the control undifferentiated (Con) SMSCs remained in complete culture medium (CCM) supplemented with 10% FBS and 1% antibiotic/antimycotic solution. On culturing days (D) 0, 1, 7, 14, and 21, cell samples from each breeds and groups were collected, then pooled across individuals for RNA and DNA extraction. A total of 36 samples were used: differentiating cells (Ost): non-differentiated cells (Con): collected at 5 time points (D 0, 1, 7, 14, 21) × 2 breeds × 2 technical replicates; collected at 4 time points (D 1, 7, 14, 21) × 2 breeds × 2 technical replicates.

Gene expression analysis of DNA methyltransferases (DNMTs) by Quantitative PCR

As previously described, RNAs were extracted by TRI reagent (Sigma-Aldrich) and reverse transcribed using SuperScript II (Invitrogen) and Oligo (dT). The geometric mean of three reference genes (HPRT1, PPIA, and YWHAZ) was used for normalization. The primers sequence was shown in Supplementary Table S1. The statistical analysis was performed using the MIXED procedure in SAS v. 9.4 (SAS Institute Inc.) to analyse the differential gene expression between breeds and the state of SMSCs (states_days) and their interactions with breed. The adjusting for multiple comparisons were done by using the post hoc Tukey–Kramer test.

DNA extraction, library construction, and RRBS

Genomic DNA was isolated using the DNeasy kit (Qiagen). For library preparation, 2 μg of DNA with a 1% spike-in control (un-methylated cl857 Sam7 Lambda DNA, Promega) digested with double enzyme (MspI and TaqI). Fragments were then end-repaired, A-tailed and ligated to c-methylated adapters using a TruSeq Nano DNA Sample Preparation kit (Illumina) according to the manufacturer’s recommendations. The Illumina TruSeq DNA library preparation kit was used for multiplexing samples per sequencing lane. After this, 2% low-range ultra-agarose gels were used to obtain proper sized adapter-ligated DNA fragments (40–240 bp). Bisulphite conversion was performed with the EpiTect Bisulphite kit (Qiagen). The bisulphite-converted libraries were then PCR amplified using a PfuTurbo Cx Hotstart DNA Polymerase kit (Stratagene). The quality of RRBS libraries was assessed using an Agilent DNA 1000 kit (Agilent Technologies). Next-generation sequencing of the RRBS libraries was performed on an Illumina HiSeq2500 at the FBN, Dummerstorf, Germany.

As previously reported, a standard analysis pipeline of RRBS sequencing has been established by our group including base calling and alignment. Briefly, in the raw fastq files, sequence reads with a mean Phred quality (Q-score) > 20, a minimum length of 30 bp without uncertain base-calling of N and adapter sequence contamination were retained. After a 2bp-trimming both 5’ and 3’ fragment ends, the clean reads were mapped to the reference genome (Sscrofa 11.1) using Bowtie2 version 2.2.8. Subsequently, read coverages and methylation percentages of each cytosine were obtained by Bismark version 0.19.0. In this study, only cytosines in CpG context were considered for further analysis, while non-CpG (CHG and CHH where H is A, C, or T) sites were discarded. In addition, to filter out SNPs which can cause bias in methylation call, we also sequenced non-bisulphite-treated version of reduced representation DNA libraries according to previously described procedure [29]. CpG sites overlapping with SNPs were filtered out.

Differentially methylation analysis and annotation

Since CpG methylation is symmetric, read counts from both strands were aggregated to form a single value for further analysis. The CpG sites with less than 2 coverage or more than the 99.9th percentile of coverage in each sample were omitted. The R package DSS (v2.48.0) was used to identify differential methylation loci (DML) and differential methylation regions (DMR). DSS-general implements a statistical method to analyse bisulphite sequencing data from experimental designs, which fits a beta-binomial generalized linear model with an arcsine link function and Wald tests [30]. In this study, the differentiation states (Ost, Con) and time-point of cell collection were considered as a combined factor (factor states_days; 9 levels: ConD0, ConD1, ConD7, ConD14, ConD21, OstD1, OstD7, OstD14, OstD21). A full model, including states_days, breeds (factor breed; 2 levels: AS, DL) and their interaction was fitted using DMLfit.multiFactor function with a smoothing algorithm. Subsequently, the contrast matrixes for comparisons of interest were created and passed to DMLtest.multiFactor function [31]. After obtaining DMLs, the callDMR function was used to get DMRs (parameters: p.threshold = 0.05). The areaStat values provided by DSS were used to rank significant DMRs and determine the direction of methylation (hyper- or hypomethylation). The minimal absolute value of areaStat was 11 in this study. The methylation levels of each DMRs were defined as weighted average of the methylation levels of CpGs located in the DMRs. The weighting factors were the coverage data of CpGs.

The annotation files of NCBI RefSeq (including 77,708 transcripts) and CpG Islands (CGIs) from genome assembly Sscrofa11.1 were downloaded from the UCSC Table Browser. All CpG sites and DMRs annotation to gene features were performed using ChIPseeker package (v3.17) [32]. For genomic features, the promoter regions were defined as ±2 kb from the Transcription Start Site (TSS). For CpG annotation, CGI shore and CGI shelf were defined as 2 kb regions flanking CGIs and 2 kb regions flanking CGI shore, respectively.

The VennDiagram package (v1.7.3) was used to draw the intersection of gene symbols of annotated DMRs from each comparisons. The hyper- and hypomethylated DMRs were plotted against the number of CpGs located in DMRs from each comparisons by EnhancedVolcano package (v1.18.0). We further calculated the mean methylation levels of each DMRs in the same groups from different comparisons. The DMRs with a relative methylation difference ≥ 30% between any of comparisons were considered for visualization by ComplexHeatmap package (v2.16.0) [33].

Functional enrichment analysis

DMRs from each comparisons were uploaded to the Database for Annotation Visualization and Integrated Discovery (DAVID) for functional and pathway enrichment analysis. The biological process (BP) enrichment and KEGG pathway were performed with default setting. Only terms with P-value <0.05 and count > 4 were considered. In addition, terms related to cancer and infection will not be discussed in this study.

Longitudinal pattern recognition

To characterize DNA methylation changing over time, we made use of Short Time-Series Expression Miner (STEM) algorithm to cluster DMRs [34]. The methylation data of each DMRs in SMSCs samples prior to (D0) and after osteogenic induction (D1, 7, 14, and 21) were fed into STEM (v1.3.12) and run by default. The median methylation level of the osteogenic induction from days 1, 7, 14 and 21 were normalized to day 0 data, which was used as the reference. The gene symbols of each significant profile (P-value <0.05) were submitted to IPA for canonical pathways analysis.

Gene expression profiling

The gene expression data of same samples used in methylation data were based on the Affymetrix porcine snowball array (SNOWBALLs520824F). The expression analysis was performed using a subset of data from previously published dataset (GEO accession number: GSE219289) [4]. To explore the relationship between DMRs and its related genes, we used Hmisc package (v5.1–1) to calculate Pearson correlation coefficients (r) and significance levels (p-value). The potential DNA methylation driven genes with negative correlation and p-value < 0.05 were shown.

Transcription factor binding site (TFBS) at DMRs

To explore the association between TFBS and DNA methylation during osteogenic differentiation, we selected promoter DMRs to determine if experimentally defined TFBS are enriched in these DMRs. The TFBS profiles were retrieved from JASPAR2022 CORE non-redundant database [35]. We used the TFBSTools package (v1.38.0) to recognize putative TFBS in these DMRs [36]. The runAme function from memes package (v1.8.0) was performed in the ‘Vs Shuffled’ mode for TFBS motif enrichment [37]. To further narrow down the result, only TFBS included in DMRs with larger than 30% of average methylation differences between any of comparisons were presented.

Results

DNMTs expression profiles and global DNA methylation during osteogenic induction of SMSCs

In the comparisons of Ost group and Con group, besides a reduced expression of DNMT3a on day 14, the expressions of all three DNMTs (DNMT1, DNMT3a, and DNMT3b) were significantly downregulated in osteogenic-induced SMSCs on day 21 (Figure 1a). We also found that SMSCs from AS pigs have higher expression of DNMT1 compared to those from DL pigs, which corresponds to higher DNA methylation in AS (Figure 1b).

Figure 1. (a) Gene expression level of three DNMTs of SMSCs grown in growth medium (Con) and osteogenic medium (Ost) over 21 days. (b) Gene expression level of three DNMTs of SMSCs derived from as and DL pigs. (c) Genic annotation of CpG sites. (d) Genomic annotation of CpG sites to DNA-methylation environments. (e) Principal component analysis (PCA) of methylation data of CpG sites. (f) Venn diagram showing unique and overlapping DMR-associated genes between control SMSCs and osteogenic induced SMSCs at each time point. (g) Venn diagram showing unique and overlapping DMR-associated genes between as and DL prior to and after osteogenic differentiation at each time point.

In total, we obtained 1 billion sequence reads for all samples with an average mapping efficiency of 53.3% (Suplementary Table S2). After quality check and filtering steps, around 2.78 million CpGs were detected in all samples. For genomic annotations, about 45.67% of these CpGs were mapped to the gene body region (introns and exons), while CpGs in the promoter and distal intergenic regions account for 26.85% and 25.07%, respectively (Figure 1c). For CpG annotations, the share of CpGs located in CpG islands was 36.01%, while CpGs in the CpG shores and shelves region account for 17.12% and 6.6%, respectively (Figure 1d). Principal component analysis (PCA) of all CpGs showed that sample segregation was more distinct between pig breeds (AS vs DL) than between SMSCs groups (Ost vs Con) for the first two principal components (Figure 1e). At the global level, median methylation rates were similar between the Ost and Con groups at all-time points, while SMSCs from AS tended to have higher median methylation rates than those from DL. In addition, promoter and CpG islands regions were uniformly hypomethylated across all samples.

Identification of DMRs in response to osteogenic induction

The numbers of DMLs (FDR <0.1) between osteogenic and control SMSCs were 18, 298, 622, and 992 from day 1 to day 21. For the comparisons between AS and DL, we found 888, 937, 1150, 1456, and 1062 DMLs (FDR <0.1) between day 0 and day 21. In view of the fact that methylation at several proximal CpG sites mostly occur in clusters throughout the genome, we will concentrate on the changes of DMRs. The Venn diagram showed that 364, 391, 442, and 517 unique gene from DMRs between Con and Ost were exclusively observed at day 1, 7, 14, and 21, respectively (Figure 1f). In addition, there were 491, 512, 502, 635, and 534 unique genes in DMRs between AS and DL at day 0, 1, 7, 14, and 21 (Figure 1g). The full list of DMRs in different comparisons was shown in Supplementary Table S3.

Differential methylation analysis revealed that the number of DMRs between osteogenic-induced SMSCs and undifferentiated SMSCs consistently increased from 539 on day 1 to 916 on day 21, mainly due to the increase in the number of hypomethylated DMRs in differentiating SMSCs from 292 on day 1 to 643 on day 21 (Figure 2a). Approximately two-thirds of DMRs showed hypomethylation in differentiating SMSCs between day 7 and day 21 relative to undifferentiated SMSCs (Figure 2a). Interestingly, in differentiating cells, the DMRs-related genes GSE1 and SLC2A1 were among the top 5 hypomethylated DMRs at most time points (Figure 2a). Comparing AS and DL, the number of DMRs (around 1400) remains stable during osteogenic differentiation, except for a slight peak at day 14. The SMSCs isolated from AS pigs always had a higher number of hypermethylated DMRs than those from DL pigs, especially at day 1 and day 21 (Figure 2b), which is consistent with the result of a higher global median methylation level in AS. The most common DMRs found at different times were TST, LOC106505785 and LOC106510599 (Figure 2b).

Figure 2. (a) Volcano plots of DMRs between control SMSCs and osteogenic induced SMSCs at each time point. (b) Volcano plots of DMRs between AS and DL prior to and after osteogenic differentiation at each time point.

Functional enrichment of DMRs

The complete results of functional enrichment analysis can be found in Supplementary Table S4. The hierarchical clustering heatmap showed that methylation differences between the state of SMSCs (Figure 3a) and between the different breeds (Figure 3b). During osteogenic induction of SMSCs, both Axon guidance and HIF-1 signalling pathway were enriched at the early stage of osteogenic induction, while extracellular matrix organization was observed from day 1 to day 14 (Figure 4a). Two pathways engaging with extracellular matrix, ECM-receptor interaction and Focal adhesion, were also found at day 14 and day 21, respectively (Figure 4a). To be specific, hypomethylated DMRs associated with extracellular matrix during osteogenic induction were functionally assigned, for instance, FBLN1, ITGB4, HSPG2, and ELN (Figure 3a, Supplementary Table S5). Other pathways that have universal impact on cell growth and differentiation were observed at the late stage of osteogenic induction, such as developmental growth, PI3K-AKt signalling, and transmembrane receptor protein tyrosine kinase signalling pathway (Figure 4a). Furthermore, several pathways that are widely perceived as contributing to osteogenic differentiation were also identified in this study, including Wnt, BMP, and TGFB signalling pathways. Most DMRs in these pathways were also hypomethylated in osteogenic SMSCs compared to undifferentiated SMSCs, covering genes such as, WWTR1, SMAD7, BMPER, and TGFB3 (Figure 3a, Supplementary Table S4).

Figure 3. (a) Methylation heatmap of 209 DMRs with relative methylation difference ≥ 30% in the promoter or gene body regions between control SMSCs and osteogenic induced SMSCs at each time point. (b) Methylation heatmap of 143 DMRs with relative methylation difference ≥ 30% in the promoter region between AS and DL prior to and after osteogenic differentiation at each time point.

Figure 4. (a) BP and KEGG pathway enrichment analysis of DMRs-associated genes between control SMSCs and osteogenic induced SMSCs at each time point. (b) BP and KEGG pathway enrichment analysis of DMRs-associated genes between AS and DL prior to and after osteogenic differentiation at each time point. (c) IPA canonical pathway analysis of DMRs-associated genes from 9 significant time-series profiles detected by STEM.

When comparing osteogenic SMSCs isolated from AS and DL, the Hippo pathway was observed throughout the induction period, except on day 1 (Figure 4b). The lipid and phospholipid pathways, namely lipid and atherosclerosis, phosphatidylinositol-mediated signalling and phosphatidylinositol phosphorylation, were enriched at day 0 and/or day 1. As osteogenesis progressed, the TGFB, Wnt and BMP signalling pathways were enriched at mid and late stages of induction. The hierarchical clustering heatmap showed that the methylation differences between AS and DL in many promoter DMRs persisted over all time points (Figure 3b).

Temporal changes in DNA methylation along osteogenic induction

We observed 9 significant time-varying methylation profiles, with profiles 9 and 41 showing the patterns of progressive decrease and increase of methylation, respectively, over time compared to the control (day 0) (Figure 4c). Profile 7 is enriched for genes associated with most canonical pathways, including CD40, IL-17A in arthritis, BMP, RANK pathways, etc (Figure 4c).

Integration of DNA methylation and gene expression

To identify potential methylation-driven genes associated with osteogenic induction, we used previously collected microarray data from the same sample set. In total, we found 340 DMRs whose associated genes showed a negative correlation between DNA methylation level and transcript abundance, with 156 and 183 genes in the differentiation state (Ost vs. Con) and breed (AS vs. DL) comparisons, respectively (Supplementary Table S5). The heatmap of genes belonging to different biological processes or pathways illustrates the relationships between the expression and methylation patterns (Figure 5a). The top 5 genes in terms of stringency of correlation between methylation and expression, located in promoter or gene body regions, were CRLF3, ELN, AFF1, BCL2L1 and FOXO3, most of which were hypomethylated in osteogenic SMSCs compared to undifferentiated SMSCs at the late stage of induction (Figure 5b, Supplementary Table S6). When comparing AS and DL, the top candidates were ARHGEF6, CARMIL1, NPR2, NTPCR and NUDT5, of which CARMIL1 consistently showed higher methylation than DL in AS at all time points (Figure 5c, Supplementary Table S5).

Figure 5. (a) Genes shown negative correlation between DNA methylation and gene expression and its enriched pathways. (b) Top 5 negative correlations between methylation levels and gene expression in the comparisons of control SMSCs and osteogenic induced SMSCs. (c) Top 5 negative correlations between methylation levels and gene expression in the comparisons of SMSCs from as and DL. (d) Most frequently enriched TF binding motifs in promoter DMRs.

DNA methylation patterns within TFBS

To better understand the regulatory potential of DMRs during osteogenic induction, we searched for consensus sequences representing transcription factor binding sites within the promoter DMRs. Comparison of cell differentiation states revealed 25 promoter DMRs in which TFBSs of 68 non-redundant transcription factors were mapped (Supplementary Table S6). The transcription factors that were enriched more than fivefold were KLF1, NFATC3, ZNF148, ASCL1, ERF:FOXI1 and KLF5 (Figure 5c). In the breed comparisons, 138 promoter DMRs were found in which TFBSs of 105 non-redundant transcription factors were enriched, of which NFATC3, KLF5 and ZNF148 were the most abundant (Figure 5d).

Discussion

We have previously investigated dynamic changes in the transcriptome profile of SMSCs derived from two different genetic backgrounds of porcine breeds during in vitro osteogenesis [4]. In addition to gene expression, DNA methylation, a mechanism of epigenetic regulation, is also associated with cell lineage commitment, including osteogenic differentiation [15,38]. Therefore, in this study, we provide further comprehensive insights into the dynamic changes in DNA methylation patterns related to osteogenic differentiation and the impact of distinct metabolic donor types on porcine SMSCs.

Regarding the global methylation pattern, our study showed that the majority of CpGs had high methylation levels, except those located in promoter regions or CpG islands, which is consistent with many previous studies [5]. Previous studies have shown that DNA methylation dynamics differ in haematopoietic stem cells derived from different sources and cell type differentiation [39,40]. They found changes in DNA methylation associated with commitment to the myeloid and lymphoid lineages [40]. In our study, we found that the origin of donor cells in terms of the breed used had a significant effect on DNA methylation changes, while a lesser change in DNA methylation was linked to a commitment to the osteogenic process.

Most DMRs were hypomethylated in osteogenic SMSCs compared to undifferentiated SMSCs, as found in this study. This reduction in DNA methylation levels coincided with lower expression levels of three key DNA methyltransferases (DNMTs) - DNMT1, DNMT3a, and DNMT3b, particularly in the late stage of osteogenic induction. Since these enzymes are crucial for establishing and maintaining DNA methylation patterns [5], it appears that the low methylation of specific genes contributes to the osteogenic transition of SMSCs. Additionally, numerous studies have reported that DNA methylation inhibitors could enhance the osteogenic capability of MSCs [8,41], and even improve their potential under unfavourable circumstances [42,43].

Similar to our earlier study in gene expression [4], we observed extracellular matrix related pathways throughout the entire induction period. The extracellular matrix is mainly produced by osteoblast before mineralization and could enhance osteogenic differentiation of MSCs [44]. ELN (Elastin) makes up part of the extracellular matrix and regulates MSCs proliferation and differentiation through its biological and mechanical properties [45]. A significant negative correlation between gene expression and methylation level of ELN has been found from day 7 to day 21 of differentiation in this study, suggesting the methylation pattern of ELN may influence its expression during osteogenic induction. Additionally, it has also been assumed that the level of ELN may act as a predictive indicator for the osteogenic ability of MSCs [46].

The commitment of MSCs to the osteogenic lineage requires the involvement of BMP and TGF-β signalling pathways [47,48]. In this study, we noticed that the BMPER, encoding a secreted protein that interacts with BMP functions, was hypermethylated on days 1 and 7. However, due to the biphasic effects of BMPER on BMPs [49,50], its exact role in osteogenic differentiation remains to be investigated. In addition, TGFB3 showed a hypomethylation pattern during the late stage of induction, whereby TGFB3 mRNA levels were upregulated in osteogenic-induced SMSCs at day 21 [4]. It has been reported that the overexpression of TGFB3 contributes to the osteogenic differentiation of MSCs [51]. Interestingly, SMAD7, an antagonist of TGF-β signalling [52], was hypomethylated in our study and showed negatively correlated gene expression. A repressive effect of SMAD7 on the differentiation and mineralization of osteoblasts has been observed [53,54]. On the other hand, SMAD7 knockout mice exhibited impaired ossification [55]. A recent study showed that SMAD7 knockdown in human SMSCs enhanced chondrogenic differentiation and suppressed endochondral ossification [56]. Our study further substantiates the role of SMAD7 during osteogenic and chondrogenic differentiation.

We also identified several upregulated genes, possibly driven by their hypomethylation pattern during the mid and late stages of osteogenic induction, such as CRLF3, AFF1, BCL2L1, and FOXO3. CRLF3 could be associated with neuronal differentiation [57], but its function in osteogenesis is currently unknown. AFF1 was regarded as a negative regulator in osteogenic differentiation [58]. Unexpectedly, our results seem to indicate a positive effect of AFF1 on osteogenic differentiation driven by methylation changes. In addition, the results for two other genes, BCL2L1 and FOXO3, are consistent with the observation that BCLXL and FOXO3 contribute to osteogenic differentiation [59,60]. Further experiments are required to explore the detailed mechanism of DNA methylation in these genes across osteogenesis.

A growing body of evidence showed that metabolic characteristics of donors, for example age and obesity, affect the methylation profile of their MSCs [61,62]. Such differences also occurred in the MSCs of pigs induced by high-fat diet [26,27]. It has also been demonstrated that obesity or metabolic syndrome could impede the osteogenic differentiation of MSCs [19,21,22]. However, it is unclear how the epigenetic changes affect the osteogenic differentiation of MSCs from donors with distinct metabolic phenotypes. DL is a modern pig breed with high lean meat, whereas AS is featured by fatness. Thus, we employed SMSCs from these two breeds to identify DNA methylation changes along with osteogenic induction. Characterization of global methylation shows an increased level of methylation in SMSCs from AS when compared to DL, consistent with that observed in skeletal muscle and similar tissues from obese and lean pig breeds [63,64], while this may be connected with the higher expression level of DNMT1, observed in this study. As a modern breed, DL has been selected for efficiency, lean meat content and commercial purposes, whereas AS has been bred in recent centuries to maintain its native characteristics. Epigenetic modifications have also contributed to breed differentiation, as shown by the comparison of non-induced SMSCs. Selective breeding strategies in the commercial breed (DL) and conservation breeding in the ancient breed (AS) have resulted in the emergence of new alleles and new epialleles. The resulting genetic differentiation affects multiple metabolic and signalling pathways, including developmental patterns. The induction of osteogenesis involves a specific spectrum of genes and molecular pathways, some of which are breed-specific and mediated by epigenetic differences, but there is also a high degree of overlap between breeds.

The osteogenesis-related signalling pathways, including TGF-beta, Wnt, and BMP, were observed exhibiting differentially methylated genes at various time points, suggesting that methylation regulates the osteogenic processes of SMSCs from AS and DL. The most frequently enriched pathway was Hippo signalling in the comparison of AS and DL. However, there was no significant change in the methylation pattern of core components in Hippo signalling [65]. On the other hand, several DMRs enriched in Hippo signalling were also key effectors of WNT signalling, for instance, WNT5A and TCF7L2. These two high-variable DMRs were hypermethylated in AS compared to DL. Experiments have shown that hypermethylation of WNT5A negatively regulated the function of limbal epithelial stem cells from diabetic donors [66]. It has also been demonstrated that the genetic variation of TCF7L2 is highly correlated with obese traits in pigs and type-2 diabetes in humans [67,68]. At the same time, both WNT5A and TCF7L2 were regarded as important regulators of osteoblast formation lately [69,70]. Regarding the identification of possible methylation-driven genes, there was very little evidence about the role of these genes detected in this study with regard to pig breeds or osteogenic differentiation. DNA methylation could directly impede the binding of transcription factors to their respective sites, thereby controlling gene expression [71]. Hence, we examined the TFBS motifs enriched in promoter DMRs. The sequence motif of two members from the krüppel-like factor (KLF) family, KLF1 and KLF5, were detected. KLFs play diverse roles in various cellular processes, including skeletal development and pathologies [72,73]. KLF1 is a vital transcription factor in erythropoiesis, yet there is no literature regarding its function in bone development [72]. KLF5 is expressed in osteoblasts and may mediate endochondral ossification [74], while knockdown of KLF5 inhibited osteogenic differentiation [75]. Moreover, NFATC3, a downstream regulator of NFATC1, was thought to maintain bone homoeostasis through the regulation of RANKL expression [76].

Conclusions

In this study, we found that the most DMRs were hypomethylated in osteogenic SMSCs compared to undifferentiated SMSCs. These findings indicate that the demethylation of specific genes plays an important role in the osteogenic transition process. In addition, we demonstrated the persistence of the breed-specific epigenetic pattern and the acquisition of an osteogenic phenotype in cell culture, while maintaining donor memory. It provides evidence that the differences in osteogenic differentiation depending on the breed origin of the donor tissue are partly mediated by epigenetic differences, which, in addition to genetic differences, are associated with breed-specific differences in metabolic type. Overall, the DNA methylation pattern of SMSCs from AS and DL pigs retains some inherent features during osteogenic induction that may confer a distinct osteogenic capacity.

Authors’ contributions

SL analysed and curated the data and wrote the original draft. PS and SP gave the conception and design. SP supervised and helped with data analysis. FH and NT helped with data analysis. SP, MO, HR and KW helped with writing – revision and editing. All authors have read and approved the published version of the manuscript.

Ethics approval and consent to participate

This study was approved by the Animal Care Committee of the Research Institute for Farm Animal Biology and subject to the approved guidelines for safeguarding good scientific practice according to German Research Foundation (DFG). The measures were taken to minimize pain and discomfort and match with the guidelines laid down by the European Directive 2010/63/EU. The animals underwent no experimental treatment, diagnostic sampling, or any other intervention before slaughter. Animal handling was in line with applicable laws, relevant guidelines, and provisions for ethical regulations.

Acknowledgments

We would like to express special thanks to Joana Bittner, Nicole Gentz, and Annette Jugert for their commitment and excellent technical assistance.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

All RRBS sequencing data have been deposited in the ArrayExpress database at EMBL-EBI (www.ebi.ac.uk/arrayexpress) under the accession number E-MTAB-9913.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592294.2024.2375011

Additional information

Funding

The Research Institute for Farm Animal Biology (FBN) provided own matched funding for experimentation and open access publication.