Identification of Co-Expression Network Correlated With Different Periods of Adipogenic and Osteogenic Differentiation of BMSCs by Weighted Gene Co-Expression Network Analysis (WGCNA)

doi:10.21203/rs.3.rs-142525/v1

Download PDF

Research Article

Identification of Co-Expression Network Correlated With Different Periods of Adipogenic and Osteogenic Differentiation of BMSCs by Weighted Gene Co-Expression Network Analysis (WGCNA)

https://doi.org/10.21203/rs.3.rs-142525/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 10 Apr, 2021

Read the published version in BMC Genomics →

You are reading this latest preprint version

Background: The differentiation of bone marrow mesenchymal stem cells (BMSCs) is a complex and dynamic process. The gene expression pattern and mechanism of different periods of adipogenic and osteogenic differentiation remain unclear. Additionally the inaction between these two lineages determination requires further exploration.

Results: Five modules that are most significantly associated with osteogenic or adipogenic differentiation of BMSCs were selected for further investigation. Biological terms, such as ribosome biogenesis, TNF-α signaling pathway, glucose import, fatty acid metabolism along with hub transcript factors, such as PPARG, YY1, and hub miRNAs, such as hsa-mir-26b-5p were enriched in different modules. The expression pattern of 6 hub genes, ADIPOQ, FABP4, SLC7A5, SELPLG, BIRC3, and KLHL30 were validated by RT-qPCR. In the end, cell staining experiments extended the findings of bioinformatics analysis.

Conclusion: This study identified the key genes, biological functions, and regulators of each time point of adipogenic and osteogenic differentiation of BMSCs and provided novel evidence and ideas for further research on the differentiation of BMSCs.

Epigenetics & Genomics

bioinformatics

BMSCs

WGCNA

osteogenic differentiation

adipogenic differentiation

Bone marrow mesenchymal stem cells (BMSCs), a typical multipotent stromal cell, have the capability to differentiate into multiple cell types, including adipocyte, osteoblasts, chondrocytes, and myocytes (Grayson et al., 2015; Polymeri et al., 2016). Dysregulation of its differentiation has been proved to be related to various diseases, such as osteoporosis, which is caused by an imbalance of the differentiation of BMSCs (Adamik et al., 2018; Verma et al., 2002). An increased capacity for adipocyte differentiation but a reduced capacity for osteoblast differentiation raises the susceptibility to brittle fracture for patients who suffer from osteoporosis (Shen et al., 2018).

It is reported that the lineage determination is a delicate balance between adipogenic and osteogenic differentiation of BMSCs (Rauch et al., 2019). The master regulator of adipogenesis, peroxisome proliferator-activated receptor γ (PPARγ), and the hub regulator of osteogenesis, runt-related transcription factor 2 (RUNX2), have been demonstrated to suppress each other(Komori et al., 1997; Tontonoz et al., 1994). However, the differentiation of BMSCs is a dynamic process that contains complex regulation and variation. Little is known about the comprehensive molecular mechanism regarding the whole procedure of adipogenic and osteogenic differentiation, hindering the development of stem cell therapy.

Several studies have individually investigated the potential biomarkers for osteogenic or adipogenic differentiation based on Differentially Expressed Genes (DEGs)(Yang et al., 2019; Zhang et al., 2018). These studies put forward that certain gene, such as LncRNA MALAT1 and MicroRNA-223 could affect the differentiation of BMSCs. Nevertheless, focusing on single gene could not demonstrate the mechanism of whole differentiation process and might neglect some meaningful regulators or pathways throughout the progress. Weighted gene co-expression network analysis (WGCNA) is a bioinformatics algorithm method which is designed to identify highly correlated gene clusters and relate them to biological traits (Langfelder and Horvath, 2008). Rather than concentrating on a single gene or isolated biomarker, WGCNA modularly investigates the co-expressed genes and extracts intramodular hub genes from systems networks, increasing the sensitivity to recognize potential worthwhile targets for biological regulations. WGCNA has been widely used for various genomic applications (Kong et al., 2020; Niemira et al., 2020).

The presented study comprehensively analysed gene expression patterns through WGCNA at each time point of adipogenic and osteogenic differentiation of BMSCs to investigate the relationship between these two lineages determination. We explored the regulation network for every period of differentiation and validated the analysis results via corresponding cell staining and RT-qPCR. Our study provides valuable research implications for the differentiation of BMSCs and prospective therapeutic targets for clinical stem cell therapy.

Data acquisition and pre-processing

The publicly available dataset GSE113253, which contained 170 samples of various mesenchymal stem cells, was downloaded from Gene Expression Omnibus (GEO). Among all the samples, 33 total RNA-sequencing data of bone marrow-derived human mesenchymal stem cells (BMSCs) were selected for further analysis. First of all, the “edgeR” package in R 3.6.3 (R Core Team, 2020, https://www.r-project.org) was applied to filter the low-expression raw data and normalise the sequencing depth difference of all the samples. R is a free software environment for statistical computing and graphics . Furthermore, we conducted hierarchical clustering and principal component analysis (PCA) to eliminate outlier samples accordingly. Finally, 31 samples (14 osteogenic BMSCs, 14 adipogenic BMSCs, and 3 undifferentiated BMSCs) were included in our analysis.

Construction of Weighted Gene Co-expression Network

We constructed the co-expression network through “WGCNA” package under R environment (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/). Firstly, the gene expression file and trait file were transformed into an appropriate format and the soft thresholding power (β value) was filtered based on the calculation of scale-free topological fit index and mean connectivity. The best β value was confirmed with a scale-free fit index bigger than 0.85 as well as the highest mean connectivity by performing a gradient test from 1 to 30. After that, the topological overlap matrix (TOM) was constructed by calculating the topological overlap between pairwise genes, and hierarchical clustering analysis was performed. The co-expression relationships among different modules were analysed and modules with high similarity were merged at the threshold of 0.25.

Module-trait correlation analysis and identification of interesting modules

To excavate interested modules which are highly related to the differentiation of BMSCs, correlation analysis between each module and different time points of adipogenic or osteogenic differentiation were conducted. This relationship is determined by Spearman’s correlation coefficient between module eigengene (ME, the major component of gene expressions in a module) and differentiation traits. Modules which had significant correlations with differentiation traits were selected for further validation.

Subsequently, gene significance (GS) and module membership (MM) were used for intramodular analysis. Gene significance is the relationship between gene expression level and differentiation trait while module membership represents the association between gene expression profile and ME of a given module. The modules which contained genes with a significant correlation between GS and MM were considered meaningful.

Enrichment analysis for biological function and Gene Set Enrichment Analysis

To investigate the biological function of the differentiation-related modules, enrichment analysis such as Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, hallmark gene set, etc., was performed through Metascape (https://metascape.org/gp/), which is a comprehensive gene annotation and analysis resource (Y. Zhou et al., 2019). Terms with a P-value < 0.01, a minimum count of 3, and an enrichment factor > 1.5 were collected and grouped into clusters based on their membership similarities.

For further verification of the function of the modules during the respective differentiation process, we conducted the Gene Set Enrichment Analysis (GSEA) for each time point of differentiation. The gene expression pattern in certain differentiation point was compared to the undifferentiated BMSCs. The criteria for statistical significance were set as p < 0.01 and FDR < 0.25.

Construction of Protein-Protein Interactions network and Hub Gene identification

To study the potential mechanism under the differentiation of BMSCs, we constructed a protein-protein interactions (PPI) network through the Search Tool Retrieval of Interacting Genes / Proteins (STRING) database, which is a trustworthy online database that can predict the physical and functional association between known and predicted protein-protein interactions.

The MCODE package in Cytoscape (version 3.8, https://cytoscape.org/) was applied to identify the hub proteins which are high degree of connectivity in the whole network. Meanwhile, we calculated the intra-modular connectivity and significance for all genes in the interested module. The genes with absolute gene significance > 0.8 and absolute intramodular connectivity > 0.8 were screened to compare with the hub genes from the PPI network. Moreover, the ”edgeR“ package was used to determine the DEGs between differentiated cells and undifferentiated BMSCs at each time point. The overlapped genes between PPI hub protein, intramodular hub genes, and DEGs were deemed to be key candidate genes that may regulate the differentiation of BMSCs.

Transcription factors and target-miRNA interaction analysis for interested modules

To recognize the crucial regulation factors for the differentiation of BMSC, we performed the enrichment analysis of transcription factor binding motifs (TFBMs) through the RcisTarget package under the R environment. The databases, hg19-tss-centered-10kb-7species.mc9nr (species = Homo sapiens, genome = hg19, distance = 10kb around the transcription start site (TSS), number of orthologous species (nOrt) =7, motif collection = Version 9 (mc9nr): 24453 motifs), were utilized to analyse the gene list of each interested module. The significant motifs based on the Normalized Enrichment Score (NES) and the genes with the best enrichment for each motif were outputted. Additionally, the TFs-genes interaction network was constructed by NetworkAnalyst (https://www.networkanalyst.ca/), which is a comprehensive network visual analytics platform for gene expression analysis.

The Tarbase v.8, a decade-long collection of experimentally supported miRNA-gene interactions, was used to seek out the potential miRNAs that can target the mRNAs within interested modules(Karagkouni et al., 2018). For improving the reliability of the predicted miRNAs, only the experimentally validated mRNAs-miRNAs pairs were included. Based on the number of targeted mRNAs, the top 5 miRNAs were extracted and the miRNA-mRNA network was visualized by Cytoscape 3.8.

Cell culture and differentiation

BMSCs extracted from the human femoral head (H. Li et al., 2016) were cultivated in Dulbecco’s Modified Eagle Medium (DMEM; 21885-025; Gibco, Germany) supplemented with 10% fetal bovine serum (FBS; P30-3702; PAN BIOTECH, Germany) and 1% penicillin-streptomycin (15140-148; Gibco, Germany) at 37°C in an incubator with 5% CO₂. The culture medium was changed every 3 days. To induce osteogenic differentiation, the culture medium was replaced by DMEM with 100nM dexamethasone (D2915-100MG; Sigma, Germany), 10nM sodium-β glycerophosphate (G5422-25G; Sigma, Germany), 0.05mM L-ascorbic acid (A8960-5G; Sigma, Germany), 10 % FBS and 1% penicillin-streptomycin. The adipogenic differentiation was induced by DMEM with 1µM dexamethasone, 0.5mM 3‑isobutyl‑1‑methylxanthine (IBMX; I5879; Sigma, Germany), 200µM indomethacin (I7378-5G; Sigma, Germany), 10µM insulin (I9278; Sigma, Germany), 10 % FBS and 1% penicillin-streptomycin. The medium was changed also every 3 days.

ALP staining, Alizarin Red S staining, and Oil Red O staining

After differentiation for a certain time, the BMSCs were fixed with 4% PFA for 30min at room temperature and washed with PBS twice. After that, cells were stained with ALP staining kit (ab242286; Abcam, Germany), Alizarin Red S (C.I. 58005; ROTH, Germany), and Oil red O (O0625-25G; Sigma, Germany) working solution for 15min, respectively. The stained cells were observed and recorded by microscopic analysis .

RNA extraction and Real-Time quantitative PCR (RT-qPCR)

Total RNA of BMSCs and differentiated cells was extracted with RNA-Solv Reagent (R6830-02; Omega, Germany) and was reverse transcribed into complementary DNA by High-Capacity RNA-to-cDNA^TM Kit (4387406; Thermo Fisher Scientific, Germany) according to the manufacture instruction. RT-qPCR was performed on a 7300 Real-Time PCR system (Applied Biosystems, Germany) using SYBR^TM Green Mastermix (4385612; Thermo Fisher Scientific, Germany). The results were normalised to the expression level of GAPDH and the relative expression levels of each gene were calculated by the 2^-ΔΔCtmethod. The sequences of primers are listed in a Supplementary table 1.

Data processing and Weighted Gene Co-expression Network Analysis (WGCNA)

Among 170 samples in GSE113253, a total of 33 samples regarding to RNA-Sequencing data of bone marrow-derived human mesenchymal stem cells were selected and underwent data filtration and normalization (Fig. 1A and B). 2 samples (RNA_14dob_BM_rep3 and RNA_1dAd_BM_rep2) were excluded because of comparatively differing from other subjects after outlier detection (Fig. 1C and D). As for the result, 31 samples were included in WGCNA.

The best soft-thresholding, 22, was chosen to construct an approximately scale-free topological overlap matrix (Fig. 1E). WGCNA identified 18 district co-expression gene modules by the Dynamic Tree Cut algorithm (Fig. 1F). Modules were automatically allocated different colours to distinguish from each other while the genes not clustered were grouped into grey module.

Identification of differentiation process-related modules

Correlation between ME and differentiation time point were evaluated through Spearman’s correlation analysis (Fig. 2A). The heatmap showed that 15 modules were significantly associated with one or more differentiation time points (p < 0.05). 4 modules with the highest correlation coefficient were selected for further research, including red, green, grey60 and tan correlated to adipogenic differentiation 14 days (Ad14d), adipogenic differentiation 4 hours (Ad4h), osteogenic differentiation 4 hours (Ob4h) and osteogenic differentiation 7 days (Ob7d), respectively. The module yellow was also included as it showed both significant relativity to adipogenic and osteogenic differentiation. The results of the intra-modular analysis demonstrated that the genes in each of the 5 modules were distinctly correlated to the corresponding differentiation time point (Fig. 2B-G), which confirmed the crucial role of 5 modules in the network of differentiation of BMSCs.

Biological function annotation of differentiation process-related modules

The enrichment analysis for the biological function of each module was performed through Metascape. The module red, which is related to Ad14d, was involved in the metabolism of lipids, lipid localization, regulation of lipid metabolic process, and PPAR signaling pathway (Fig. 3A and F). This was consistent with the later stage characteristics of adipogenic differentiation. Besides, GSEA revealed that adipogenesis and oxidative phosphorylation were significantly enriched in the 14th day of adipogenic differentiation (Fig. 3K). The most enriched terms of module tan were cell morphogenesis involved in differentiation, muscle system process, divalent inorganic cation homeostasis, and regulation of ion transport (Fig. 3B and G). GSEA indicated the genes in module tan may regulate the mid-term osteogenic differentiation through glucose import and fatty acid metabolism (Fig. 3L).

Ad4h correlated green modules were enriched in ribosome biogenesis and cell part morphogenesis (Fig. 3C and H) which was consistent with the GSEA result (Fig. 3M), indicating the morphology change during the early stage of adipogenic differentiation. Both enrichment analysis (Fig. 3D and I) and GSEA (Fig. 3N) of Ob4h correlated gery60 module showed regulation of cytokine production and TNF-α signaling pathway play important roles in the early stage of osteogenic differentiation. The module yellow, correlated to both adipogenic and osteogenic differentiation were involved in the regulation of neuron differentiation, response to growth factor, and negative regulation of growth (Fig. 3E and J), suggesting that these terms possibly affect the direction of differentiation of BMSCs.

PPI network analysis and hub gene identification

To explore the interaction among genes in each module, PPI network was constructed by STRING database (Supplementary Figure S1-S5). Hub gene clusters that score > 3 within each PPI network were identified using the Cytoscape MCODE plugin (Fig. 4B).

“EdgeR” package was applied to investigate the DEGs between genes at each differentiation time point and undifferentiated BMSCs with the thresholds of P-value < 0.05 and |Fold Change (FC)| > 2.0. In total, 2253 DEGs (1037 upregulated and 1216 downregulated) in Ad14d, 1047 DEGs (472 upregulated and 575 downregulated) in Ob7d, 373 DEGs (208 upregulated and 165 downregulated) in Ob4h and 899 DEGs (416 upregulated and 484 downregulated) in Ad4h were identified, respectively (Fig. 4A).

Intra-module hub genes, which possessed high connectivity in each module, were filtered at the threshold of absolute gene significance > 0.8 and absolute intramodular connectivity > 0.8. The intramodular hub genes in each module were listed in Supplementary Table S2. PPI hub cluster genes and DEGs at each differentiation time point were overlapped with high connectivity genes in the respective module. As showed in Fig. 4C, 15 hub genes were overlapped between module red and Ad14d, include ACACB, GPAM, ADIPO1, FABP4, etc.. 2 hub genes, KLHL30 and MYOZ2, were overlapped between module tan and Ob7d. Module grey60 and OB4h DEGs shared 2 hub genes: BIRC3 and PTGS2. Module green and Ad4h DEGs shared 3 hub genes: MYLIP, SLC7A5, and DLX5. Module yellow overlapped with Ad4h DEGs by BDNF and SELPLG while did not overlap with Ob4h DEGs.

Regulation mechanism analysis of differentiation-related modules

Given the modules were consisted of co-expressed genes, which may be regulated by a common mechanism such as transcription factors and miRNAs. The transcription factor binding motifs (TFMFs) enrichment analysis showed that transcript factor PPARG was the crucial regulator for the red module, which is annotated to the motif transfac_public__M00528. The TFs-genes interaction network analysis for module red indicated HNF4A, YY1, FOXC1 played important roles in the late stage of adipogenic differentiation as well (Fig. 5B). For module tan, ZNF232, which was annotated to the motif taipale__ZNF232_full_RTGTTAAAYGTAGATTAAG_repr, was the master regulator. GATA2, YY1 was also found to be a key component of the TFs network of tan (Fig. 5E), suggesting that they were essential for the regulation of the mid-term of osteogenic differentiation. The transcription factor CREB1 aligned with the motif transfac_pro__M03544 in module green was required for the induction of adipogenic differentiation (Fig. 5A). The transcript factors NFKB1 and RELA which corresponded to motif totransfac_public__M00054 were enriched both in TFMFs and TF-genes interaction analysis in module grey60, revealed their vital role in the induction of osteogenic differentiation (Fig. 5D). For module yellow, which was correlated to both adipogenic and osteogenic differentiation, the transcript factor ELK1 and its binding motif dbcorrdb__ELK1__ENCSR000EFU_1__m2 were significantly enriched (Fig. 5C). All the enriched motifs and their high confidence transcription factors were shown in Supplemental Digital Content (SDC1-SDC5).

MiRNAs were another fundamental regulators through recognition of cognate sequence and participate in transcriptional, translational, or epigenetic processes. Consequently, we constructed a miRNAs-mRNAs interaction network based on experimentally validated miRNA-target pairs in Tarbase v8.0. The top5 miRNAs regulating the most number of genes were highlighted in each module (Fig. 6). Hsa-mir-26b-5p functioned in all 5 modules, indicating its pivotal influence in the differentiation of BMSCs. Other hub miRNAs, including hsa-mir-335-5p, hsa-mir-16-5p, and has-mir-124-3p, also regulated in both adipogenic differentiation-related and osteogenic differentiation-related modules, suggesting the intimate interaction between two directions of the differentiation of BMSCs.

Differentiation-related staining and hub gene validation

As shown in Fig. 7B, the morphology of BMSCs started to gradually transform from arborisation into round from 4 hours of adipogenic differentiation, which is coincident with the biological function enrichment analysis of module green. Lip droplets existed from the 7th day of adipogenic differentiation and continuously accumulated until the 14th day. The result of Oil Red O staining also illustrated the transformation process of adipogenic differentiation (Fig. 7E).

For osteogenic differentiation, BMSCs initiated the morphological change and mineralization deposit after 7 days (Fig. 7A). Alizarin Red S staining demonstrated that mineralization accumulation increased rapidly from the 7th day to the 14th day of osteogenic differentiation (Fig. 7D). However, ALP production of BMSCs raised from 4 hours after osteogenic differentiation and reached the peak at the 7th day (Fig. 7C). After that, the level of ALP began to decrease and disappeared on the 14th day. Such phenomenon revealed that the 7th day was a crucial turning point of osteogenic differentiation, which was corresponding to our WGCNA results.

Finally, we validated the hub gene expression level during the differentiation of BMSCs. 6 hub genes from 5 modules were selected, including ADIPOQ, FABP4, SLC7A5, SELPLG, BIRC3, and KLHL30. All the hub genes showed significant overexpression at their respective differentiation time point. Their expression trend was consistent with the RNA sequencing results in dataset GSE113253 (Fig. 7F and G), which certified that the results of this study are reliable and accurate.

Recent studies demonstrated that the osteogenic and adipogenic lineages could alternate during cell differentiation, indicate the subtle and complicate relevance between them (Jing et al., 2016; C. J. Li et al., 2015). The exact mechanism functioned in the differentiation process remains unclear. In this study, we identified several modules that highly correlated to diverse stages of differentiation and performed enrichment analysis for each module. Hub genes and crucial regulation factors, such as miRNAs and TFs, were recognized from the networks.

The grey60 module was significantly associated with 4 hours after osteogenic differentiation. Enrichment analysis and GSEA for grey60 module revealed that TNF-α via NF-κB signaling pathway played an important role in the early stage of osteogenic differentiation. Interestingly, TNF-α was widely accepted as an inhibitor of osteogenic differentiation and osteogenesis among previous studies(Deng et al., 2018; Zhao et al., 2011). However, some recently contradictory findings suggested the paradoxical effects of TNF-α in the regulation of bone homeostasis (Osta et al., 2014). Daniele et al. demonstrated that the effect of TNF-α on osteogenic differentiation was dose-dependent (Daniele et al., 2017). A low concentration of TNF-α showed enhanced osteogenic differentiation of BMSCs while a high concentration indicated the opposite result. Moreover, Huang et al. pointed out that treatment time also affected the function of TNF-α in osteogenic differentiation(Huang et al., 2011). Besides, the hub gene BIRC3 and the hub transcript factor RELA achieved from our grey60 module could be activated by NF-κB through the TNF-α signaling pathway (Wang et al., 2012). The expression of BIRC3, as showed in our RT-qPCR result, impressively increased at 4 hours of osteogenic differentiation and decreased immediately afterward, which was perfectly consistent with the dual effect of TNF-α in osteogenic differentiation.

Mid-term osteogenic differentiation was correlated to the module tan, within which, transcript factors FOXC1, FOXL1, YY1, GATA2, and miRNAs hsa-mir-335-5p, hsa-mir-92a-3p, hsa-mir-16b-5p, hsa-mir-4459, and hsa-mir-124-3p were master regulators. KLHL30 and MYOZ2 were overlapped between Ob7d DEGs, intramodular hub genes, and hub PPI cluster. Cell morphogenesis involved in differentiation was enriched in this module, which was also observed under a microscope during our differentiation induction. The 7th day of osteogenic differentiation was a watershed for ALP production and mineralization deposit. GSEA results demonstrated that glucose import and fatty acid metabolism were significantly enriched in Ob7d compared to undifferentiated BMSCs. Glucose has been recognized to be an essential nutrient for osteoblasts and many other studies have already focused on the glutamine metabolism in bone homeostasis (Al-Qarakhli et al., 2019; T. Zhou et al., 2019). Fatty acids were also exhibited close implication in osteogenesis (Bao et al., 2020). Our research confirms the importance of nutrient metabolism and proposes a possible time point for them to function during the process of osteogenic differentiation.

According to the results of cell staining and differentiation induction, adipogenesis experienced considerably more obvious and faster changes than osteogenesis at the very early stage of differentiation. The sample cluster analysis and PCA supported that the transcriptome of osteoblasts clustered more closer to BMSCs than that of adipocytes after 4 hours of differentiation. The analysis of module green, which was highly related to Ad4h, showed that ribosome biogenesis was significantly enriched. Besides, ribosome biogenesis was both enriched in Ob4h and Ad4h with GSEA, indicating its decisive role in the fate of cell differentiation. Our conclusion is in agreement with previous studies in other stem cells (Gabut et al., 2020; Sanchez et al., 2016). Three hub genes, SLC7A5, MYLIP, and DLX5 were overlapped between module green and DEGs of Ad4h. Moreover, transcript factor CREB1 annotated to motif transfac_pro__M03544, and transcript factor ATF1, annotated to motif transfac_pro__M07034, were identified from the TFMFs enrichment analysis for module green. Further research on these hub genes and hub transcript factors are required for uncovering the mechanism of early adipogenic differentiation of BMSCs.

Compared to the early stage of adipogenic differentiation, the late period of adipogenesis mainly focuses on dealing with lipids. The top activated biological terms in module red were metabolism of lipids, lipid localization, lipid catabolic process, and regulation of lipid metabolic process. The upset plot showed many genes are overlapped in these terms. A total of 15 hub genes were recognized from module red, among which, PPARG, ADIPOQ, LPL, and FABP4 were already well-known biomarkers of adipocytes (Cristancho and Lazar, 2011; Floresta et al., 2017). Other hub genes in module red, such as LIPE, had intimate interactions with these validated biomarkers, suggesting their unignorable position during the late stage of adipogenesis. LIPE, a kind of lipase encoding gene, has been reported to be involved in various lipid metabolism-related syndrome(Sollier et al., 2021; Wu et al., 2017). However, few study reports the relationship between LIPE and adipogenic differentiation. Our TFMFs enrichment analysis of red module also obtained some motifs that aligned with transcript factor PPARG, including transfac_public__M00528, cisbp__M3785, and cisbp__M6433. These enriched motifs along with the hub miRNAs in the miRNAs-mRNAs network in module red, including hsa-mir-335-5p, hsa-mir-124-3p, hsa-mir-26-5p, hsa-mir-6499-3p, hsa-mir-16-5p, could be the master regulators of the late stage of adipogenic differentiation.

Module yellow was the only module that was significantly correlated to both adipogenic and osteogenic differentiation, which supported the point that these two lineages determination interacted with each other closely. Many common genes were involved in both differentiation and affected the destiny of BMSCs. As a result, the enriched terms in module yellow, for example, regulation of neuron differentiation, response to growth factor, and negative regulation of growth could function in early-stage differentiation of BMSCs. Nevertheless, hub genes, SELPLG and BNDF, were only overlapped between module yellow and DEGs of Ad4h, but not DEGs of Ob4h. This phenomenon may result from the situation we mentioned above that adipogenesis changed faster and greater than osteogenesis.

For the regulation network of differentiation, we found many hub transcription factors and hub miRNAs existed in more than one module, such as YY1, GATA2, E2F1, hsa-mir-26b-5p, hsa-mir-335-5p, and hsa-mir-16-5p, suggesting these master regulators could be active in the whole process of differentiation. Although, constrained by the length of the article and practical situation, some of the findings could not be explained and examined in detail there is no doubt that our study supplied some potential research targets for further investigation in the differentiation of BMSCs despite the limitations.

Taken together, we applied WGCNA for exploring transcriptome data of various time points of differentiation of BMSCs and identified 5 modules that were significantly correlated to different stages of osteogenic or adipogenic differentiation. The pivotal biological terms, hub genes, and master regulators for each time point of differentiation were predicted through bioinformatics analysis Meanwhile, cell experiments were conducted as the verification for some of them. We are looking forward our study could provide novel evidences and ideas for further research on the differentiation of BMSCs.

BMSCs Bone marrow mesenchymal stem cells

WGCNA Weighted gene co-expression network analysis

Ad adipogenic differentiation

ALP Alkaline Phosphatase

DEGs Differentially Expressed Genes

DMEM Dulbecco’s Modified Eagle Medium

FBS fetal bovine serum

GEO Gene Expression Omnibus

GO Gene Ontology

GS Gene significance

GSEA Gene Set Enrichment Analysis

KEGG Kyoto Encyclopedia of Genes and Genomes

MM module membership

NES Normalized Enrichment Score

Ob osteogenic differentiation

PCA Principal Component Analysis

PPI Protein-Protein Interactions

RT-qPCR Real-Time quantitative PCR

TFs transcription factors

TFBMs transcription factor binding motifs

TOM topological overlap matrix

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and materials

The dataset analysed during the current study are available in the GEO repository, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE113253

Competing interests

The authors declare that they have no competing interests.

Funding

Not applicable.

Authors' contributions

Yu Liu contributes to the design of the paper, drafting the manuscript, performing the experiments and data analysis. Markus Tingart contributes to the generation of the idea and reviewing the manuscript. Sophie Lecouturier contributes to performing cell experiment. Jianzhang Li contributes to data collection and analysis. Jörg Eschweiler contributes to the design and generation of the idea and reviewing the manuscript.

Acknowledgements

Not applicable.

Adamik, J., et al. (2018). Osteoblast suppression in multiple myeloma bone disease. Journal of Bone Oncology, 13, 62-70. doi:10.1016/j.jbo.2018.09.001
Al-Qarakhli, A. M. A., et al. (2019). Effects of high glucose conditions on the expansion and differentiation capabilities of mesenchymal stromal cells derived from rat endosteal niche. BMC Mol Cell Biol, 20(1), 51. doi:10.1186/s12860-019-0235-y
Bao, M. Y., et al. (2020). Therapeutic potentials and modulatory mechanisms of fatty acids in bone. Cell Proliferation, 53(2). doi:ARTN e1273510.1111/cpr.12735
Cristancho, A. G., and Lazar, M. A. (2011). Forming functional fat: a growing understanding of adipocyte differentiation. Nature Reviews Molecular Cell Biology, 12(11), 722-734. doi:10.1038/nrm3198
Daniele, S., et al. (2017). Osteogenesis Is Improved by Low Tumor Necrosis Factor Alpha Concentration through the Modulation of Gs-Coupled Receptor Signals. Mol Cell Biol, 37(8). doi:10.1128/MCB.00442-16
Deng, L., et al. (2018). Involvement of microRNA-23b in TNF-alpha-reduced BMSC osteogenic differentiation via targeting runx2. J Bone Miner Metab, 36(6), 648-660. doi:10.1007/s00774-017-0886-8
Floresta, G., et al. (2017). Adipocyte fatty acid binding protein 4 (FABP4) inhibitors. A comprehensive systematic review. European Journal of Medicinal Chemistry, 138, 854-873. doi:10.1016/j.ejmech.2017.07.022
Gabut, M., et al. (2020). Ribosome and Translational Control in Stem Cells. Cells, 9(2). doi:ARTN 49710.3390/cells9020497
Grayson, W. L., et al. (2015). Stromal cells and stem cells in clinical bone regeneration. Nat Rev Endocrinol, 11(3), 140-150. doi:10.1038/nrendo.2014.234
Huang, H., et al. (2011). Dose-specific effects of tumor necrosis factor alpha on osteogenic differentiation of mesenchymal stem cells. Cell Prolif, 44(5), 420-427. doi:10.1111/j.1365-2184.2011.00769.x
Jing, H., et al. (2016). Suppression of EZH2 Prevents the Shift of Osteoporotic MSC Fate to Adipocyte and Enhances Bone Formation During Osteoporosis. Molecular Therapy, 24(2), 217-229. doi:10.1038/mt.2015.152
Karagkouni, D., et al. (2018). DIANA-TarBase v8: a decade-long collection of experimentally supported miRNA-gene interactions. Nucleic Acids Res, 46(D1), D239-D245. doi:10.1093/nar/gkx1141
Komori, T., et al. (1997). Targeted disruption of Cbfa1 results in a complete lack of bone formation owing to maturational arrest of osteoblasts. Cell, 89(5), 755-764. doi:Doi 10.1016/S0092-8674(00)80258-5
Kong, Y., et al. (2020). Identification of Immune-Related Genes Contributing to the Development of Glioblastoma Using Weighted Gene Co-expression Network Analysis. Frontiers in Immunology, 11. doi:ARTN 128110.3389/fimmu.2020.01281
Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. Bmc Bioinformatics, 9. doi:Artn 55910.1186/1471-2105-9-559
Li, C. J., et al. (2015). MicroRNA-188 regulates age-related switch between osteoblast and adipocyte differentiation. Journal of Clinical Investigation, 125(4), 1509-1522. doi:10.1172/Jci77716
Li, H., et al. (2016). Isolation and characterization of primary bone marrow mesenchymal stromal cells. Ann N Y Acad Sci, 1370(1), 109-118. doi:10.1111/nyas.13102
Niemira, M., et al. (2020). Molecular Signature of Subtypes of Non-Small-Cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-Expression Network Analysis (WGCNA). Cancers, 12(1). doi:ARTN 3710.3390/cancers12010037
Osta, B., et al. (2014). Classical and Paradoxical Effects of TNF-alpha on Bone Homeostasis. Frontiers in Immunology, 5, 48. doi:10.3389/fimmu.2014.00048
Polymeri, A., et al. (2016). Bone Marrow Stromal Stem Cells in Tissue Engineering and Regenerative Medicine. Horm Metab Res, 48(11), 700-713. doi:10.1055/s-0042-118458
Rauch, A., et al. (2019). Osteogenesis depends on commissioning of a network of stem cell transcription factors that act as repressors of adipogenesis (vol 51, pg 716, 2019). Nature Genetics, 51(4), 766-766. doi:10.1038/s41588-019-0400-4
Sanchez, C. G., et al. (2016). Regulation of Ribosome Biogenesis and Protein Synthesis Controls Germline Stem Cell Differentiation. Cell Stem Cell, 18(2), 276-290. doi:10.1016/j.stem.2015.11.004
Shen, G. S., et al. (2018). The GDF11-FTO-PPARgamma axis controls the shift of osteoporotic MSC fate to adipocyte and inhibits bone formation during osteoporosis. Biochim Biophys Acta Mol Basis Dis, 1864(12), 3644-3654. doi:10.1016/j.bbadis.2018.09.015
Sollier, C., et al. (2021). LIPE-related lipodystrophic syndrome: clinical features and disease modeling using adipose stem cells. Eur J Endocrinol, 184(1), 155-168. doi:10.1530/EJE-20-1013
Tontonoz, P., et al. (1994). Stimulation of adipogenesis in fibroblasts by PPAR gamma 2, a lipid-activated transcription factor. Cell, 79(7), 1147-1156. doi:10.1016/0092-8674(94)90006-x
Verma, S., et al. (2002). Adipocytic proportion of bone marrow is inversely related to bone formation in osteoporosis. Journal of Clinical Pathology, 55(9), 693-698. doi:DOI 10.1136/jcp.55.9.693
Wang, Y., et al. (2012). Gene network revealed involvements of Birc2, Birc3 and Tnfrsf1a in anti-apoptosis of injured peripheral nerves. PLoS One, 7(9), e43436. doi:10.1371/journal.pone.0043436
Wu, J. W., et al. (2017). Epistatic interaction between the lipase-encoding genes Pnpla2 and Lipe causes liposarcoma in mice. PLoS Genet, 13(5), e1006716. doi:10.1371/journal.pgen.1006716
Yang, X., et al. (2019). LncRNA MALAT1 shuttled by bone marrow-derived mesenchymal stem cells-secreted exosomes alleviates osteoporosis through mediating microRNA-34c/SATB2 axis. Aging (Albany NY), 11(20), 8777-8791. doi:10.18632/aging.102264
Zhang, S., et al. (2018). MicroRNA-223 Suppresses Osteoblast Differentiation by Inhibiting DHRS3. Cell Physiol Biochem, 47(2), 667-679. doi:10.1159/000490021
Zhao, L., et al. (2011). Tumor Necrosis Factor Inhibits Mesenchymal Stem Cell Differentiation into Osteoblasts Via the Ubiquitin E3 Ligase Wwp1. Stem Cells, 29(10), 1601-1610. doi:10.1002/stem.703
Zhou, T., et al. (2019). Glutamine Metabolism Is Essential for Stemness of Bone Marrow Mesenchymal Stem Cells and Bone Homeostasis. Stem Cells International, 2019. doi:Artn 892893410.1155/2019/8928934
Zhou, Y., et al. (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun, 10(1), 1523. doi:10.1038/s41467-019-09234-6

supplementarymaterilas.pdf

Download PDF

Journal Publication

published 10 Apr, 2021

Read the published version in BMC Genomics →

Editorial decision: Major revision
29 Jan, 2021
Reviews received at journal
14 Jan, 2021
Reviewers agreed at journal
08 Jan, 2021
Reviewers invited by journal
08 Jan, 2021
Editor assigned by journal
08 Jan, 2021
Editor invited by journal
08 Jan, 2021
Submission checks completed at journal
08 Jan, 2021
First submitted to journal
07 Jan, 2021

You are reading this latest preprint version

Identification of Co-Expression Network Correlated With Different Periods of Adipogenic and Osteogenic Differentiation of BMSCs by Weighted Gene Co-Expression Network Analysis (WGCNA)

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Materials And Methods

Results

Data processing and Weighted Gene Co-expression Network Analysis (WGCNA)

Identification of differentiation process-related modules

Biological function annotation of differentiation process-related modules

PPI network analysis and hub gene identification

Regulation mechanism analysis of differentiation-related modules

Differentiation-related staining and hub gene validation

Discussion

Conclusion

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1