- Open Access
Transcriptome-wide m6A methylome during osteogenic differentiation of human adipose-derived stem cells
Stem Cell Research & Therapy volume 12, Article number: 489 (2021)
Adipose-derived stem cells are frequently used for bone regeneration both in vitro and in vivo. N6-methyladenosine (m6A) is the most abundant post-transcriptional modification on eukaryotic RNAs and plays multifaceted roles in development and diseases. However, the regulatory mechanisms of m6A in osteogenic differentiation of human adipose-derived stem cells (hASCs) remain elusive. The present study aimed to build the transcriptome-wide m6A methylome during the osteogenic differentiation of hASCs.
Materials and methods
hASCs were harvested after being cultured in a basic or osteogenic medium for 7 days, and the osteogenic differentiation was validated by alkaline phosphatase (ALP) and Alizarin Red S staining, ALP activity assay, and qRT-PCR analysis of ALP, RUNX2, BGLAP, SPP1, SP7, and COL1A1 genes. The m6A level was colorimetrically measured, and the expression of m6A regulators was confirmed by qRT-PCR and western blot. Moreover, m6A MeRIP-seq and RNA-seq were performed to build the transcriptome and m6A methylome. Furthermore, bioinformatic analyses including volcano plots, Venn plots, clustering analysis, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, gene sets enrichment analysis, and protein-protein interaction analysis were conducted.
In total, 1145 differentially methylated peaks, 2261 differentially expressed genes, and 671 differentially methylated and expressed genes (DMEGs) were identified. GO and KEGG pathway analyses conducted for these DMEGs revealed extensive and osteogenic biological functions. The “PI3K-Akt signaling pathway”; “MAPK signaling pathway”; “parathyroid hormone synthesis, secretion, and action”; and “p53 signaling pathway” were significantly enriched, and the DMEGs in these pathways were identified as m6A-specific key genes. A protein-protein interaction network based on DMEGs was built, and VEGFA, CD44, MMP2, HGF, and SPARC were speculated as the hub DMEGs.
The total m6A level was reduced with osteogenic differentiation of hASCs. The transcriptome-wide m6A methylome built in the present study indicated quite a few signaling pathways, and hub genes were influenced by m6A modification. Future studies based on these epigenetic clues could promote understanding of the mechanisms of osteogenic differentiation of hASCs.
Large bone defects, caused by trauma, bone loss, and tumors, lead to heavy social and economic burdens. At present, the clinical gold standard for the treatment of skeletal defects is autogenous or allogeneic bone grafts, which suffers disadvantages of a limited amount of harvested bone and donor site morbidity [1, 2]. Bone tissue engineering, built on stem cells or osteoprogenitor cells, osteoconductive biomaterials, and osteoconductive cytokines facilitating cellular proliferation and osteogenic differentiation, offers a promising solution [2, 3]. Human bone marrow-derived mesenchymal stem cells (hBMSCs), naturally resident in the bone marrow, are pluripotent and have been widely applied in bone tissue engineering for years [4, 5]. However, the access to autogenous hBMSCs is invasive and subject to limited cell incidence. What is more, the self-renewal and proliferative ability of hBMSCs weakens with donor aging and diseases such as osteoporosis and arthritis [6,7,8].
Adipose stromal cells, available from lipoaspirate, are a cellular niche containing differentiated cells, committed progenitors, and mesenchymal stem cells . Among the multiple cell types, adipose-derived stem cells (ADSCs) are pluripotent, with adipogenic, chondrogenic, neurogenic, and osteogenic potentials. In the presence of ascorbate, glycerophosphate, and dexamethasone, ADSCs can differentiate into osteoblast-like cells [10,11,12]. Meanwhile, ADSCs have high proliferation rates and immunosuppressive properties and can secrete numerous polypeptides, hormones, and effective growth factors to induce osteogenesis. Moreover, ADSCs can also secrete a mass of suitable chemokines to recruit endogenous stem cells into the bone defect site [12,13,14]. With the increasing incidence of obesity worldwide, subcutaneous adipose tissue is easy to obtain, making human adipose-derived stem cells (hASCs) a first-rate alternative for hBMSCs in bone regenerative medicine .
Eukaryotic RNAs can be modified post-transcriptionally by over 170 modifications, among which, N6-methyladenosine (m6A), mostly resident at the stop codon and the 3′ untranslated region (UTR), is the most abundant [15, 16]. The m6A modification can be dynamically regulated by methyltransferases (also known as “writers”, including METTL3, METTL14, etc.) and demethylases (also known as “erasers,” including FTO, ALKBH5) and can further be recognized by specific RNA-binding proteins (also known as “readers,” including YTHDF1, YTHDF2) [17, 18]. As the first known form of reversible mRNA modification, m6A exerts effects on the stability, alternative splicing, export, cellular localization, and translation of RNAs, therefore influencing gene expression [17,18,19,20]. Till now, this epigenetic mark on RNAs has been extensively reported in virtually all major biological processes, normal development, and various diseases [21,22,23,24,25,26].
In the case of bone biology, m6A modification plays critical roles in bone development as well as bone diseases such as osteoporosis, osteoarthritis, and osteosarcoma [27, 28]. Mechanistically, m6A was found executing dual functions in the osteogenic differentiation of MSCs. On the one hand, m6A methylation enhanced osteogenic differentiation by regulating miR-320, RUNX2, Pth1r, PI3K/Akt, VEGF, or SMAD7/SMURF1 signaling, and inhibited adipogenic differentiation by regulating JAK1, Pth1r, SRSF2/ RUNX1T1, or TRAF4/ PKM2/β-catenin signaling [27,28,29,30]. On the other hand, the m6A methylation was also reported to decrease osteogenic differentiation by regulation of MYD88/NF-κB, miR-7212-5p/FGFR3, or FTO/p-AMPK loop [27, 28, 31]. However, all these findings above are rooted in BMSCs, and little is known about the scenario in ADSCs at the present stage.
In the present study, we aimed to investigate the roles of m6A in the osteogenic differentiation of hASCs. High-throughput sequencing for RNA (RNA-seq) and methylated RNA immunoprecipitation and sequencing (m6A MeRIP-Seq) were performed to comprehensively identify differentially regulated m6A peaks and messenger RNAs between two groups: undifferentiated hASCs (uhASCs) and osteogenically differentiated hASCs (dhASCs). The expressions of osteogenic markers, m6A regulating genes, and hub genes were validated by quantitative real-time polymerase chain reaction (qRT-PCR). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed for the differentially methylated genes (DMGs), differentially expressed genes (DEGs), and differentially methylated and expressed genes (DMEGs) to reveal the biological significance of m6A in osteogenic differentiation of hASCs. Moreover, we also built the protein-protein interaction (PPI) network based on the DMEGs and calculated the hub genes in this process. These findings could provide a new perspective on the molecular mechanisms of osteogenic differentiation in hASCs.
Materials and methods
The primary hASCs from healthy adult donors in this study were obtained from Cyagen Biosciences (Suzhou, China). The cells were cultured with OriCell™ Human ASC Growth Medium (Cyagen) containing 10% fetal bovine serum, 1% penicillin-streptomycin, and 1% glutamine. The medium was renewed every 3 days until the cells were passaged upon approximately 80% confluence. All the cells were cultured in an incubator with 5% CO2 at 37 °C. As our preliminary study indicated, the hASCs expressed a high level of stem cell surface markers CD29, CD44, CD105, and HLA-ABC, while less than 5% expressed CD31, CD34, CD45, and HLA-DR . The third to the seventh generation of hASCs were used for osteogenic differentiation.
The hASCs were seeded at a density of 2.0 × 105 and 1.0 × 105 cells/well into six- and twelve-well plates, respectively. For osteogenic differentiation, hASCs were cultured in OriCell™ hASC Osteogenic Differentiation Medium (Cyagen) containing 10% fetal bovine serum, 1% penicillin-streptomycin, 10 mM β-glycerophosphate, 10 mM glutamine, 50 μM ascorbate, and 0.1 μM dexamethasone. For the control group, hASCs were cultured in the general growth medium as mentioned before. The medium was renewed every 2 days. Both groups were detected by alkaline phosphatase (ALP) and Alizarin Red S (ARS) staining, ALP activity assay, and qRT-PCR analysis of ALP, RUNX2, BGLAP, SPP1, SP7, and COL1A1 genes [33,34,35].
ALP and ARS staining
ALP staining was conducted after 10 days of induction of osteogenic differentiation. Cells were fixed with 4% polyoxymethylene for 15 min and incubated with 0.1 M Tris buffer (pH 9.3) containing 0.25% naphthol AS-BI phosphate (Sigma) and 0.75% Fast Blue BB (Sigma) at 37 °C for 45 min. The Alkaline Phosphatase Assay Kit (Beyotime, China) was adopted according to the manufacturer’s instruction, and the ALP activity was quantified by a spectrophotometer (Thermo Fisher Scientific) at 405 nm.
ARS staining was conducted after 14 days of induction of osteogenic differentiation. Cells were fixed with 4% polyoxymethylene for 15 min and stained with 1% Alizarin Red S (ARS) Stain Solution (Cyagen) at room temperature for 5 min. And then, the solution was removed, and the samples were rinsed 3–4 times with deionized water. Alizarin Red-positive area was analyzed using the ImageJ software (NIH) and shown as a percentage of Alizarin Red-positive area over the total area.
Quantitative real-time polymerase chain reaction analysis
The expressions of osteogenic markers, m6A regulating genes, and hub DMEGs were validated by qRT-PCR. Total RNA from hASCs was extracted using TRIzol Reagent (Invitrogen) according to the manufacturer’s instructions. The RNA purity and concentration were measured with a NanoDrop 2000 (Thermo Fisher Scientific). Reverse transcription was performed by PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara, Japan). qRT-PCR reaction was performed using TB Green™ Premix Ex Taq™ II (Takara) in QuantStudio 3 Real-Time PCR Systems (Thermo Fisher Scientific). Relative gene expression was normalized by GAPDH using a 2−ΔΔCt method. The primers used in this study are listed in Supplementary Table S1.
Cells after 7 days of osteogenic induction were lysed in RIPA buffer (Beyotime, China) on ice. The samples were heated at 95 °C for 5 min in a sample buffer containing 20% SDS-PAGE sample loading buffer, separated on 10% SDS-polyacrylamide gels, and transferred to the PVDF membranes by a wet transfer apparatus (Bio-Rad). The membranes were blocked with 5% BSA for 1 h and then incubated overnight with rabbit METTL3 polyclonal antibody (1:1000, Cell Signaling Technology, Cat No: 96391), METTL14 polyclonal antibody (1:1000, ABclonal, Cat No: A8530), FTO polyclonal antibody (1:1000, Affinity, Cat No: DF8421), ALKBH5 polyclonal antibody (1:2000, Proteintech, Cat No: 16837-1-AP), and α-tubulin rabbit polyclonal antibody (1:2000, Beyotime, Cat No: AF0001), followed by incubation with a goat anti-rabbit IgG secondary antibody HRP conjugated (1:5000, Signaling antibody, Cat No: L3012-2). The antibody-antigen complexes were visualized with Bio-Rad Quantity One (Bio-Rad Laboratories Inc., Hercules, CA, USA).
m6A level quantification
Total RNA from hASCs was extracted using TRIzol Reagent (Invitrogen) according to the manufacturer’s instructions. The RNA purity and concentration were measured with a NanoDrop 2000 (Thermo Fisher Scientific). The m6A RNA Methylation Quantification Kit (P-9005, EpiQuik) was used to measure the m6A content in the total RNAs following the manufacturer’s instruction. The m6A levels were quantified colorimetrically by reading the absorbance of each well at a wavelength of 450 nm, and the m6A level was then calculated based on the standard curve.
m6A MeRIP-Seq and RNA-seq
The m6A MeRIP-Seq and RNA-Seq were performed after 7 days of induced osteogenic differentiation. Total RNA was extracted using TRIzol Reagent (Invitrogen, CA, USA) and was tested for quality and quantity with Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number > 7.0. Then, poly-T oligo-attached magnetic beads (Invitrogen) were used to isolate poly(A) mRNA from the total RNA. Next, the poly(A) mRNA fractions are fragmented into ~ 100-nt-long oligonucleotides and then incubated for 2 h at 4 °C with m6A-specific antibody (No. 202003, Synaptic Systems, Germany) in IP buffer (50 mM Tris-HCl, 750 mM NaCl and 0.5% Igepal CA-630) supplemented with BSA (0.5 μg μl−1). The mixture was then incubated with protein-A beads and eluted with elution buffer (1 × IP buffer and 6.7 mM m6A). Eluted RNA was precipitated by 75% ethanol. Both eluted m6A-containing fragments (IP) and untreated control fragments (input) are converted to the final cDNA library. The average insert size for the paired-end libraries was ~ 100 ± 50 bp. And then the paired-end 2 × 150 bp sequencing was performed on an Illumina Novaseq™ 6000 platform at the LC-BIO Bio-Tech Ltd. (Hangzhou, China) following the vendor’s recommended protocol.
First, we used Cutadapt and perl scripts inhouse to remove the reads containing adaptor contamination, low-quality bases, and undetermined bases. Then, sequence quality was validated by fastp. Next, HISAT2 was used to map reads to the genome of Homo sapiens (version: v96). The R package exomePeak was adopted to process mapped reads of IP and input libraries and identifies m6A peaks with bed or bam format that can be adapted for visualization on the IGV software (http://www.igv.org/). HOMER and MEME were used for de novo and known motif finding followed by localization of the motif with respect to peak summit by perl. Called peaks were annotated by intersection with gene architecture using ChIPseeker. Then, StringTie was used to perform the expression level for all mRNAs from input libraries by calculating FPKM. The differentially expressed mRNAs were selected with the standard of p-value < 0.05 and log2 (fold change) > 1 or < −1 and by the R package edgeR.
Venn analysis was performed to characterize overlapped dysregulated m6A peaks and mRNAs between dhASCs and uhASCs detected by MeRIP-Seq and RNA-Seq. The number of m6A peaks and mRNAs with or without overlapping were shown in different colors in a pie diagram. Up- and downregulated m6A peaks and mRNAs were analyzed separately.
GO and KEGG analyses
GO and KEGG pathway analyses were performed on m6A peaks and mRNAs. GO analysis, including molecular function, biologic process, and cellular component, was performed. p < 0.05 denoted the significance of GO term enrichment. Based on the KEGG database, pathway analysis was performed to analyze the potential functions. GO and pathway enrichment analyses were also performed based on the differentially regulated m6A peaks and mRNAs.
Gene set enrichment analysis
For gene set enrichment analysis (GSEA), gene lists “c2.cp.kegg.v7.3.symbols” were obtained from MSigDB, and the analysis was performed with the GSEA software (http://www.broad.mit.edu/GSEA). In short, we imported our gene lists of interest into the software and examined the significance of the given gene sets. p-values were computed using a bootstrap distribution created by resampling gene sets of the same cardinality.
Construction of the protein-protein interaction network
The protein-protein interaction (PPI) network was constructed to assess the interactions between the differentially methylated and expressed genes (DMEGs) based on the correlation analysis. All the DMEGs were selected to construct the network using Cytoscape 3.8 (Institute of Systems Biology in Seattle). Then, we calculated the hub genes by combining scores from all DMEGs.
Statistical tests were performed using R version 3.6.0, SPSS 24.0 (IBM, NY, USA) and GraphPad Prism 8.0. Numerical data was presented as means ± standard deviation. Differences between the two groups were evaluated for statistical significance by a t-test. The appropriate analysis of variance was used when more than two groups were evaluated, followed by Dunnett’s multiple comparisons test to compare the difference between the groups. p ≤ .05 was considered statistically significant. The R packages VennDiagram, ComplexHeatmap, clusterProfiler, and ggplot2 were used for graphing and data visualization.
Osteogenic differentiation of hASCs was successfully induced.
The ALP activity, validated by ALP staining and ALP activity quantification assay, was significantly higher in the dhASCs group than in the uhASCs group. A similar result regarding mineral nodes forming was observed in ARS staining (Fig. 1a–d). Moreover, compared to the control group, the expressions of ALP, RUNX2, SPP1, SP7, and COL1A1 were significantly increased in the dhASCs group (Fig. 1e). These results confirmed the osteogenic differentiation of dhASCs.
m6A level significantly decreased during osteogenic differentiation of hASCs
To learn about the m6A methylation in the process of osteogenic differentiation of hASCs, we first detected the relative m6A level between the groups. The m6A quantification assay analysis revealed a significantly reduced m6A level during the osteogenic induction of hASCs from day 0 to day 14, with the most significant and steepest reduction at the early stage (Fig. 2a, b). Then, we investigated the expression of m6A regulators. According to the RNA-seq, compared with the uhASCs group, FTO expression was significantly increased in the dhASCs group at day 7. However, the expression of METTL3, METTL14, and ALKBH5 showed no significant differences between the groups (Fig. 2c). Further, the qRT-PCR revealed the expression of METTL3 was significantly elevated at day 7 and day 14, and the expression of FTO was significantly elevated at day 7 and significantly decreased to baseline at day 14 (Fig. 2d). Moreover, the expressions of METTL3, METTL14, FTO, and ALKBH5 between the groups were confirmed by western blotting (Fig. 2e). In a word, these results demonstrated potential effects of m6A in the osteogenic differentiation of hASCs.
m6A peaks were significantly dysregulated during osteogenesis of hASCs
To dissect the underlying roles of m6A in the osteogenesis of hASCs, we performed the m6A MeRIP-Seq to depict the critical picture of m6A methylome. For the whole transcriptome, 39,735 (20431 transcripts) and 44,630 (22789 transcripts) m6A peaks were detected in uhASCs and dhASCs, respectively. For uhASCs, 45.73% of the peaks were detected in the 3′UTR region, 18.39% in the 5′UTR region, 10.6% in the first exon, and 25.28% in the remaining exons (Supplementary Figure S1a-b). For dhASCs, 43.85% of the peaks were detected in the 3′UTR region, 20.48% in the 5′UTR region, 11.06% in the first exon, and 24.61% in the remaining exons (Supplementary Figure S1c-d). According to the MeRIP-Seq, by the standard of p-value < 0.05 and fold change (fc) ≥ 2, 1145 m6A peaks were significantly dysregulated after a 7-day osteogenic induction of hASCs, of which 534 peaks were hypermethylated and 611 hypomethylated (Fig. 3a). In addition, 48.81% of the dysregulated peaks were in the 3′UTR region, 19.46% in the 5′UTR region, 9.77% in the first exon, and 21.96% in the remaining exons (Fig. 3b, c). The top 20 DMGs were listed in Table 1, and the full list of DMGs was shown in Supplementary Table S2. Furthermore, the motif analysis validated m6A RRACH (R was purine, A was m6A, and H was a non-guanine bas) consensus sequence enrichment between uhASCs and dhASCs (Fig. 3d).
The DMGs were associated with extensive and osteogenic biological functions
To figure out the biological significance of m6A methylation in the osteogenesis of hASCs, we performed GO and KEGG pathway analyses for DMGs. GO analysis uncovered that both the hypermethylated and hypomethylated genes in dhASCs were significantly involved with “regulation of transcription, DNA-templated,” “regulation of transcription by RNA polymerase II,” “signal transduction,” “cell differentiation,” and “multicellular organism development” (ontology: biological process); “membrane,” “nucleus,” and “cytoplasm” (ontology: cellular component); and “protein binding,” “metal ion binding,” and “DNA binding” (ontology: molecular function) (Fig. 4a, c). Remarkably, KEGG pathway analysis demonstrated that the hypermethylated genes were significantly associated with signaling pathway “thyroid hormone signaling pathway,” “Rap1 signaling pathway,” and “FoxO signaling pathway” (Fig. 4b). However, the hypomethylated genes were not significantly associated with the common osteogenic signaling pathway (Fig. 4d).
Osteogenic pathways and genes were significantly differentially methylated and regulated
To figure out the alteration of gene expression and signaling pathways in the osteogenesis of hASCs, we performed RNA-seq to figure out the DEGs. By the standard of p-value < 0.05 and fc ≥ 2, compared with uhASCs, 2261 genes were significantly dysregulated in the dhASCs group, of which 1210 were upregulated and 1051 were downregulated (Supplementary Figure S2a). The hierarchical clustering of the top 50 DEGs was shown in Supplementary Figure S2b and the full list of DEGs was shown in Supplementary Table S3. GO analysis and KEGG pathway analysis of DEGs were shown in Supplementary Figure S2b-e. By conjoint analysis of DMGs (by the standard of p-value < 0.05) and DEGs (by the standard of p-value < 0.05 and fold change (fc) ≥ 2), we found 671 differentially methylated and expressed transcripts (598 genes) were mainly divided into four groups, including 140 hypermethylated and upregulated (hyper-up), 208 hypomethylated and downregulated (hypo-down), 105 hypermethylated but downregulated (hyper-down), and 218 hypomethylated but upregulated genes or transcripts (hypo-up) (Fig. 5a, b). The hierarchical clustering of the top 50 DMEGs was shown in Fig. 5c. The top 20 DMEGs were listed in Table 2, and the full list of DMEGs was shown in Supplementary Table S4. GO analysis uncovered that the DMEGs in dhASCs were significantly involved with “signal transduction,” “multicellular organism development,” “regulation of transcription, DNA-templated,” “cell differentiation,” “cell adhesion,” and “regulation of cell population proliferation” (ontology: biological process); “membrane,” “cytoplasm,” and “integral component of membrane” (ontology: cellular component); and “protein binding,” “metal ion binding,” and “DNA binding” (ontology: molecular function) (Fig. 5d). Remarkably, the KEGG pathway analysis demonstrated that the DMEGs were significantly associated with signaling pathway “PI3K-Akt signaling pathway,” “MAPK signaling pathway,” “parathyroid hormone synthesis, secretion and action,” and “p53 signaling pathway” (Fig. 5e). The DMEGs associated with these four signaling pathways were selected and performed a hierarchical cluster analysis (Fig. 5f–i). Furthermore, GSEA analysis of the whole transcriptome revealed these genes were significantly associated with the “Hedgehog signaling pathway” and the “TGF-beta signaling pathway” (Fig. 5j, k).
PPI network and hub genes were discovered among DMEGs
The PPI network of the DMEGs was conducted by STRING database and Cytoscape (Fig. 6a). The network was grouped into four clusters as mentioned before, that is, hyper-up, hypo-down, hyper-down, and hypo-up DMEGs (Fig. 6b–e). GO enrichment analyses were performed for each cluster of DMEGs to elucidate their biological function (Fig. 6b–e). Hub genes in the PPI network were selected using cytoHubba by calculating the protein combining degree (Fig. 6f). And the m6A peak distribution of the top 5 hub genes between the groups was shown in Fig. 6g, and the expressions of the 5 hub DMEGs were shown in Fig. 6h. The PPI network interaction metadata was listed in Supplementary Table S5.
Adipose-derived stem cells (ADSCs), as one form of mesenchymal stem cells with good accessibility and special biological properties, are frequently used for bone regeneration both in vitro and in vivo [11, 12]. According to the current criteria, the isolation of ADSCs produces cell type heterogeneity, and the purity of ADSCs is essential for their engineering use . In the present study, the hASCs were identified with surface markers and successfully induced into osteoblast-like cells [9, 32]. Identifying factors involved in the osteogenic differentiation of ASCs is important for better understanding the mechanism of osteogenic differentiation. Recent studies have shown that m6A modification plays essential roles in stem cell fate decisions in a stage-, state-, and cell type-specific manner [26, 36,37,38,39,40]. In the case of osteogenic differentiation, reduced m6A level in BMSCs inhibited osteogenic differentiation by affecting the PTH/Pth1r signaling axis, the PI3K-AKT signaling pathways, or the alternative splicing of VEGFA [30, 40, 41]. As for preadipocytes, some studies reported FTO promoted adipogenesis by regulating cell cycle-related genes and the JAK2-STAT3-C/EBPb signaling pathway in an m6A-dependent manner in 3T3-L1 cell lines [42,43,44,45]. However, at the present stage, little is known about the roles of m6A modification in the osteogenic differentiation of ADSCs. In this study, using m6A MeRIP-Seq and RNA-seq, we depicted the transcriptome-wide m6A methylome in the process of osteogenic differentiation of hASCs for the first time.
In the present study, we observed that from day 0 to day 14, the level of m6A modification decreased continuously with the most significant and steepest reduction at the early stage. Therefore, we investigated the time point day 7 to explore potential m6A-related signaling pathways. Combined with significant expression changes of m6A regulators METTL3 and FTO, we unveiled the potential effects of m6A modification in osteogenic differentiation of hASCs. The changes of m6A level during osteogenic induction of hASCs appeared stage-specific. This was probably because of the time-dependent expression pattern of m6A regulators. At the early stage, the elevated expression of demethylase FTO predominantly reduced the m6A level. With the osteogenic induction lasting, the expression of FTO decreased while the expression of methyltransferase METTL3 maintained, resulting in a slower reduction of m6A level. Interestingly, this trend was contrary to that of BMSCs where a time-dependent increase of m6A content during osteogenic differentiation was reported . In the cases of BMSCs, from day 0 to day 14, the expression of methyltransferase METTL3 increased while the expression of demethylase FTO and ALKBH5 did not change significantly . The differences between hASCs and BMSCs might reflect the cell type-specific effects of m6A modification in stem cell differentiation [26, 40]. Although both were mesenchymal stem cells, the distinct origins may cause different gene expression pattern and epigenetic modifications, thus producing divergent stemness, proliferation capacity, multilineage differentiation potential, longevity, immunoregulatory properties, and so on [6, 46]. Considering the complex stage- and state-specific roles of m6A modification, the m6A methylome of a prolonged osteogenic induction of hASCs remains to be revealed.
In total, the m6A MeRIP-Seq revealed an abundance of 1.94 and 1.96 peaks per transcript for uhASCs and dhASCs, respectively. Additionally, the m6A peaks of both uhASCs and dhASCs were mostly enriched in the 3′UTR region and near the stop codon; these sites were m6A-specific and consistent with previous studies [15,16,17,18,19, 47]. Furthermore, during the osteogenic differentiation of hASCs, the differentially methylated peaks were significantly enriched in “signal transduction,” “regulation of transcription,” “multicellular organism development,” “cell differentiation,” and so on, indicating the conservative and fundamental roles of m6A in regulating development and cell fate specification [16, 17].
In order to elucidate the mechanism of m6A in affecting osteogenic differentiation of hASCs, we combined the m6A methylome and the transcriptome and found the key signaling pathways that were affected by m6A modification. As with previous studies, several canonical pathways regarding osteogenic differentiation were significantly enriched. The most significant ones were the “PI3K-Akt signaling pathways,” “MAPK signaling pathways,” “parathyroid hormone synthesis, secretion and action,” and “P53 signaling pathways,” in which quite a few key genes were differentially methylated [20, 22, 30, 40, 41, 48]. Moreover, although other common osteogenic pathways were not significantly enriched in the m6A methylome, quite a few key genes in these pathways were significantly methylated. At the same time, we constructed the PPI network of DMEGs and speculated the hub genes based on combining degree, for example, the genes of VEGFA, CD44, MMP2, HGF, and SPARC. All these differentially expressed as well as differentially methylated key genes and hub genes provided the latest and detailed information about how m6A modification affected the osteogenic differentiation of hASCs.
A mass of studies in the past few years have lighted up the biological effects of m6A modification on RNAs [19,20,21,22, 40, 47, 49,50,51,52]. On the one hand, the m6A methylation process is reversable, and this mark on RNA could be written or erased in case of various stimuli and biological factors [17, 47, 49]. On the other hand, m6A can affect RNA processing and metabolism in multifaceted mechanisms, including alternative splicing, alternative polyadenylation, RNA stability, RNA export, RNA degradation, and translation. Therefore, m6A up- or downregulates the gene expression in a complex and context-depended manner [17,18,19]. For this reason, we observed four groups of DMEGs in the present study, those were hyper-up, hypo-down, hyper-down, and hypo-up DMEGs. Our function enrichment analyses suggested that these four clusters of DMEGs were associated with the fundamental and different biological processes. Combined with the methylation and expression pattern of the key DMEGs, future studies based on these molecular clues could give us brand-new knowledge on the osteogenic differentiation of hASCs from an m6A-specific epigenetic perspective.
Total m6A level was reduced with osteogenic differentiation of hASCs, and the transcriptome-wide m6A methylome built in the present study indicated quite a few signaling pathways and hub genes were influenced by m6A modification. In total, 1145 DMGs, 2261 DEGs, and 671 DMEGs were detected. GO and KEGG pathway analysis conducted for these differentially methylated and regulated genes revealed extensive and osteogenic biological functions. The “PI3K-Akt signaling pathway,” “MAPK signaling pathway,” “parathyroid hormone synthesis, secretion and action,” and “p53 signaling pathway” were significantly enriched, and the DMEGs in these pathways were identified as m6A-specific key genes. PPI network based on DMEGs was built, and VEGFA, CD44, MMP2, HGF, and SPARC were speculated as the hub DMEGs. Future studies based on these epigenetic clues could promote the understanding of the mechanisms of osteogenic differentiation of hASCs.
Availability of data and materials
All sequencing data have been deposited into the Sequence Read Archive (SRA) database of NCBI with the identifier PRJNA739515 (https://dataview.ncbi.nlm.nih.gov/object/PRJNA739515?reviewer=jqgj6n7reasr8ei9a20riogqok).
Adipose-derived stem cells
Alizarin Red S
Differentially expressed genes
Osteogenically differentiated hASCs
Differentially methylated and expressed genes
Differentially methylated genes
Gene set enrichment analysis
Human adipose-derived stem cells
Human bone marrow-derived mesenchymal stem cells
Kyoto Encyclopedia of Genes and Genomes
- m6A MeRIP-Seq:
Methylated RNA immunoprecipitation and sequencing
Quantitative real-time polymerase chain reaction
High-throughput sequencing for RNA
Sándor GK, Numminen J, Wolff J, Thesleff T, Miettinen A, Tuovinen VJ, et al. Adipose stem cells used to reconstruct 13 cases with cranio-maxillofacial hard-tissue defects. Stem Cells Transl Med. 2014;3(4):530–40. https://doi.org/10.5966/sctm.2013-0173.
Borrelli MR, Hu MS, Longaker MT, Lorenz HP. Tissue engineering and regenerative medicine in craniofacial reconstruction and facial aesthetics. J Craniofac Surg. 2020;31(1):15–27. https://doi.org/10.1097/SCS.0000000000005840.
Petrovic V, Zivkovic P, Petrovic D, Stefanovic V. Craniofacial bone tissue engineering. Oral Surg Oral Med Oral Pathol Oral Radiol. 2012;114(3):e1–9. https://doi.org/10.1016/j.oooo.2012.02.030.
Jin YZ, Lee JH. Mesenchymal stem cell therapy for bone regeneration. Clin Orthop Surg. 2018;10(3):271–8. https://doi.org/10.4055/cios.2018.10.3.271.
Liao HT, Chen CT. Osteogenic potential: comparison between bone marrow and adipose-derived mesenchymal stem cells. World J Stem Cells. 2014;6(3):288–95. https://doi.org/10.4252/wjsc.v6.i3.288.
Mazini L, Rochette L, Amine M, Malka G. Regenerative capacity of adipose derived stem cells (ADSCs), comparison with mesenchymal stem cells (MSCs). Int J Mol Sci. 2019;20(10).
Ho-Shui-Ling A, Bolander J, Rustom LE, Johnson AW, Luyten FP, Picart C. Bone regeneration strategies: engineered scaffolds, bioactive molecules and stem cells current stage and future perspectives. Biomaterials. 2018;180:143–62. https://doi.org/10.1016/j.biomaterials.2018.07.017.
Szpalski C, Barbaro M, Sagebin F, Warren SM. Bone tissue engineering: current strategies and techniques--part II: cell types. Tissue Eng Part B Rev. 2012;18(4):258–69. https://doi.org/10.1089/ten.teb.2011.0440.
Squillaro T, Peluso G, Galderisi U. Clinical trials with mesenchymal stem cells: an update. Cell Transplant. 2016;25(5):829–48. https://doi.org/10.3727/096368915X689622.
Sheykhhasan M, Wong JKL, Seifalian AM. Human adipose-derived stem cells with great therapeutic potential. Curr Stem Cell Res Ther. 2019;14(7):532–48. https://doi.org/10.2174/1574888X14666190411121528.
Tajima S, Tobita M, Mizuno H. Current status of bone regeneration using adipose-derived stem cells. Histol Histopathol. 2018;33(7):619–27. https://doi.org/10.14670/HH-11-942.
Shafaei H, Kalarestaghi H. Adipose-derived stem cells: an appropriate selection for osteogenic differentiation. J Cell Physiol. 2020;235(11):8371–86. https://doi.org/10.1002/jcp.29681.
Si Z, Wang X, Sun C, Kang Y, Xu J, Wang X, et al. Adipose-derived stem cells: sources, potency, and implications for regenerative therapies. Biomed Pharmacother. 2019;114:108765. https://doi.org/10.1016/j.biopha.2019.108765.
Bacakova L, Zarubova J, Travnickova M, Musilkova J, Pajorova J, Slepicka P, et al. Stem cells: their source, potency and use in regenerative therapies with focus on adipose-derived stem cells - a review. Biotechnol Adv. 2018;36(4):1111–26. https://doi.org/10.1016/j.biotechadv.2018.03.011.
Huang H, Weng H, Chen J. The biogenesis and precise control of RNA m6A methylation. Trends Genet. 2020;36(1):44–52. https://doi.org/10.1016/j.tig.2019.10.011.
Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature. 2012;485(7397):201–6. https://doi.org/10.1038/nature11112.
Zaccara S, Ries RJ, Jaffrey SR. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol. 2019;20(10):608–24. https://doi.org/10.1038/s41580-019-0168-5.
Shi H, Wei J, He C. Where, when, and how: context-dependent functions of RNA methylation writers, readers, and erasers. Mol Cell. 2019;74(4):640–50. https://doi.org/10.1016/j.molcel.2019.04.025.
Jiang X, Liu B, Nie Z, Duan L, Xiong Q, Jin Z, et al. The role of m6A modification in the biological functions and diseases. Signal Transduct Target Ther. 2021;6(1):74. https://doi.org/10.1038/s41392-020-00450-x.
Zhang W, He L, Liu Z, Ren X, Qi L, Wan L, et al. Multifaceted functions and novel insight into the regulatory role of RNA N6-methyladenosine modification in musculoskeletal disorders. Front Cell Dev Biol. 2020;8:870. https://doi.org/10.3389/fcell.2020.00870.
Zhou Z, Lv J, Yu H, Han J, Yang X, Feng D, et al. Mechanism of RNA modification N6-methyladenosine in human cancer. Mol Cancer. 2020;19(1):104. https://doi.org/10.1186/s12943-020-01216-3.
Xu F, Li W, Yang X, Na L, Chen L, Liu G. The roles of epigenetics regulation in bone metabolism and osteoporosis. Front Cell Dev Biol. 2020;8:619301.
Shulman Z, Stern-Ginossar N. The RNA modification N6-methyladenosine as a novel regulator of the immune system. Nat Immunol. 2020;21(5):501–12. https://doi.org/10.1038/s41590-020-0650-4.
McFadden MJ, Horner SM. N6-methyladenosine regulates host responses to viral infection. Trends Biochem Sci. 2020.
Madugalle SU, Meyer K, Wang DO, Bredy TW. RNA N6-methyladenosine and the regulation of RNA localization and function in the brain. Trends Neurosci. 2020;43(12):1011–23. https://doi.org/10.1016/j.tins.2020.09.005.
Vu LP, Pickering B, Cheng Y, Zaccara S, Nguyen D, Minuesa G, et al. The N6-methyladenosine (m6A)-forming enzyme METTL3 controls myeloid differentiation of normal hematopoietic and leukemia cells. Nat Med. 2017;23(11):1369–76. https://doi.org/10.1038/nm.4416.
Chen J, Tian Y, Zhang Q, Ren D, Zhang Q, Yan X, et al. Novel insights into the role of N6-methyladenosine RNA modification in bone pathophysiology. Stem Cells Dev. 2020.
Chen X, Hua W, Huang X, Chen Y, Zhang J, Li G. Regulatory role of RNA N6-methyladenosine modification in bone biology and osteoporosis. Front Endocrinol. 2019;10:911.
Yan G, Yuan Y, He M, Gong R, Lei H, Zhou H, et al. m6A methylation of precursor-miR-320/RUNX2 controls osteogenic potential of bone marrow-derived mesenchymal stem cells. Mol Ther Nucleic Acids. 2020;19:421–36. https://doi.org/10.1016/j.omtn.2019.12.001.
Wu Y, Xie L, Wang M, Xiong Q, Guo Y, Liang Y, et al. Mettl3-mediated m6A RNA methylation regulates the fate of bone marrow mesenchymal stem cells and osteoporosis. Nat Commun. 2018;9(1):4772. https://doi.org/10.1038/s41467-018-06898-4.
Mi B, Xiong Y, Yan C, Chen L, Xue H, Panayi AC, et al. Methyltransferase-like 3-mediated N6-methyladenosine modification of miR-7212-5p drives osteoblast differentiation and fracture healing. J Cell Mol Med. 2020;24(11):6385–96. https://doi.org/10.1111/jcmm.15284.
Luo T, Yang X, Sun Y, Huang X, Zou L, Liu J. Effect of microRNA-20a on osteogenic differentiation of human adipose tissue-derived stem cells. Cells Tissues Organs. 2019;208(3-4):148–57. https://doi.org/10.1159/000506304.
Di Bernardo G, Messina G, Capasso S, Del Gaudio S, Cipollaro M, Peluso G, et al. Sera of overweight people promote in vitro adipocyte differentiation of bone marrow stromal cells. Stem Cell Res Ther. 2014;5(1):4. https://doi.org/10.1186/scrt393.
Xia K, Cen X, Yu L, Huang X, Sun W, Zhao Z, et al. Long noncoding RNA expression profiles during the NEL-like 1 protein-induced osteogenic differentiation. J Cell Physiol. 2020;235(9):6010–22. https://doi.org/10.1002/jcp.29526.
Huang XQ, Cen X, Sun WT, Xia K, Yu LY, Liu J, et al. CircPOMT1 and circMCM3AP inhibit osteogenic differentiation of human adipose-derived stem cells by targeting miR-6881-3p. Am J Transl Res. 2019;11(8):4776–88.
Zhang C, Chen Y, Sun B, Wang L, Yang Y, Ma D, et al. m6A modulates haematopoietic stem and progenitor cell specification. Nature. 2017;549(7671):273–6. https://doi.org/10.1038/nature23883.
Bertero A, Brown S, Madrigal P, Osnato A, Ortmann D, Yiangou L, et al. The SMAD2/3 interactome reveals that TGFβ controls m6A mRNA methylation in pluripotency. Nature. 2018;555(7695):256–9. https://doi.org/10.1038/nature25784.
Yang D, Qiao J, Wang G, Lan Y, Li G, Guo X, et al. N6-methyladenosine modification of lincRNA 1281 is critically required for mESC differentiation potential. Nucleic Acids Res. 2018;46(8):3906–20. https://doi.org/10.1093/nar/gky130.
Lee H, Bao S, Qian Y, Geula S, Leslie J, Zhang C, et al. Stage-specific requirement for Mettl3-dependent m6A mRNA methylation during haematopoietic stem cell differentiation. Nat Cell Biol. 2019;21(6):700–9. https://doi.org/10.1038/s41556-019-0318-1.
Zhang M, Zhai Y, Zhang S, Dai X, Li Z. Roles of N6-methyladenosine m6A in stem cell fate decisions and early embryonic development in mammals. Front Cell Dev Biol. 2020;8:782. https://doi.org/10.3389/fcell.2020.00782.
Tian C, Huang Y, Li Q, Feng Z, Xu Q. Mettl3 regulates osteogenic differentiation and alternative splicing of Vegfa in bone marrow mesenchymal stem cells. Int J Mol Sci. 2019;20(3).
Wu R, Guo G, Bi Z, Liu Y, Zhao Y, Chen N, et al. m6A methylation modulates adipogenesis through JAK2-STAT3-C/EBPβ signaling. Biochim Biophys Acta Gene Regul Mech. 2019;1862(8):796–806. https://doi.org/10.1016/j.bbagrm.2019.06.008.
Liu Q, Zhao Y, Wu R, Jiang Q, Cai M, Bi Z, et al. ZFP217 regulates adipogenesis by controlling mitotic clonal expansion in a METTL3-m6A dependent manner. RNA Biol. 2019;16(12):1785–93. https://doi.org/10.1080/15476286.2019.1658508.
Wu R, Liu Y, Yao Y, Zhao Y, Bi Z, Jiang Q, et al. FTO regulates adipogenesis by controlling cell cycle progression via m6A-YTHDF2 dependent mechanism. Biochim Biophys Acta Mol Cell Biol Lipids. 2018;1863(10):1323–30. https://doi.org/10.1016/j.bbalip.2018.08.008.
Zhao X, Yang Y, Sun BF, Shi Y, Yang X, Xiao W, et al. FTO-dependent demethylation of N6-methyladenosine regulates mRNA splicing and is required for adipogenesis. Cell Res. 2014;24(12):1403–19. https://doi.org/10.1038/cr.2014.151.
Lee RH, Kim B, Choi I, Kim H, Choi HS, Suh K, et al. Characterization and expression analysis of mesenchymal stem cells from human bone marrow and adipose tissue. Cell Physiol Biochem. 2004;14(4-6):311–24. https://doi.org/10.1159/000080341.
He PC, He C. m6A RNA methylation: from mechanisms to therapeutic potential. EMBO J. 2021;40(3):e105977. https://doi.org/10.15252/embj.2020105977.
Zhang Y, Gu X, Li D, Cai L, Xu Q. METTL3 regulates osteoblast differentiation and inflammatory response via Smad signaling and MAPK signaling. Int J Mol Sci. 2019;21(1).
Kim J, Lee G. Metabolic control of m6A RNA modification. Metabolites. 2021;11(2).
Zhang Y, Xu S, Xu G, Gao Y, Li S, Zhang K, et al. Dynamic expression of m6A regulators during multiple human tissue development and cancers. Front Cell Dev Biol. 2020;8:629030.
Zhang L, Hou C, Chen C, Guo Y, Yuan W, Yin D, et al. The role of N6-methyladenosine m6A modification in the regulation of circRNAs. Mol Cancer. 2020;19(1):105. https://doi.org/10.1186/s12943-020-01224-3.
Yi YC, Chen XY, Zhang J, Zhu JS. Novel insights into the interplay between m6A modification and noncoding RNAs in cancer. Mol Cancer. 2020;19(1):121. https://doi.org/10.1186/s12943-020-01233-2.
This study was supported by grants from the National Natural Science Foundation of China (81870743, 81470722, and 81771048) and the Creative Spark Fund of Sichuan University (2018SCUH0007).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
m6A peak calling for uhASCs and dhASCs. a Pie charts showing distribution of m6A peaks of uhASCs in different gene context. b Accumulation of m6A peaks of uhASCs along transcripts. c Pie charts showing distribution of m6A peaks of dhASCs in different gene context. d Accumulation of m6A peaks of dhASCs along transcripts. hASCs, human adipose-derived stem cells; uhASCs, undifferentiated hASCs; dhASCs, osteogenically differentiated hASCs.
Gene expression profile during osteogenesis of hASCs by RNA-seq. a Volcano plots displaying the DEGs (fold change ≥ 2 and p < 0.05). b Hierarchical clustering analysis of the top 50 DEGs (ranking by p value). c GO enrichment analysis for upregulated genes. d KEGG pathway analysis for upregulated genes. e GO enrichment analysis for downregulated genes. f KEGG pathway analysis for downregulated genes. DEGs, differentially expressed genes; hASCs, human adipose-derived stem cells; uhASCs, undifferentiated hASCs; dhASCs, osteogenically differentiated hASCs; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; FC, fold change.
Primers used for qRT-PCR.
Full list of DMGs (uhASCs vs dhASCs).
Full list of DEGs (uhASCs vs dhASCs).
Full list of DMEG (uhASCs vs dhASCs).
The PPI network metadata.
About this article
Cite this article
Sun, W., Song, Y., Xia, K. et al. Transcriptome-wide m6A methylome during osteogenic differentiation of human adipose-derived stem cells. Stem Cell Res Ther 12, 489 (2021). https://doi.org/10.1186/s13287-021-02508-1
- Adipose-derived stem cell
- Osteogenic differentiation