Weighted gene correlation network analysis reveals novel biomarkers associated with mesenchymal stromal cell differentiation in early phase

Introduction

Osteoporosis is a systemic metabolic bone disease characterized by reduction in bone mass and degeneration of bone microstructure, which makes the bone brittle and prone to fracture. Osteoporosis has become an important public health problem as the global population ages. Therefore, revealing the molecular mechanisms of osteoporosis and developing effective and preventive treatments would be crucial to the human well-being. Bone remodeling occurs at discrete sites called bone remodeling units, wherein mineralized bone resorption occurs by osteoclasts that break down the tissue in bones and release the minerals to the blood. Osteoblasts are then recruited to the site after osteoclasts undergo apoptosis, and a new bone is eventually formed (Compston, McClung & Leslie, 2019). The imbalance between bone resorption and bone formation can lead to various diseases, such as osteopenia, osteopetrosis and osteoporosis. While osteoclasts are derived from hematopoietic stem cell precursors, osteoblasts and adipocytes are derived from bone marrow mesenchymal stem cells (MSCs) (Chen et al., 2016b). MSCs are mesoderm-derived multipotent cells that adhere to culture dish surface and can proliferate and differentiate into adipocytes, osteoblasts and chondrocytes in vitro. Aging is associated with reduced osteogenic (OS) potential and increased fat formation by MSCs, which can lead to senile osteoporosis (Bergman et al., 1996).

With the advent of the era of big data, a variety of biological public databases have emerged, which provide a large amount of accessible genomic and clinical data for basic and translational studies. Today, the Gene Expression Omnibus (GEO) database is considered the most comprehensive public repository of large-scale genomics data. These databases allow for mining meaningful genomic changes and for discovering biological mechanisms that are involved in disease development and progression (Gruebner et al., 2017). In addition, mining public databases can provide insights into the functional networks of unknown genes, and thus, confirming gene expression trends that may be overlooked in a single experiment. Besides, large-scale analysis can validate experimental findings, provide supportive evidence, help design a research plan, and test an already set hypothesis (Lachmann et al., 2018).

Weighted Gene Co-Expression Network Analysis (WGCNA) is a widely used technique for transcriptomic data analysis (Langfelder & Horvath, 2008). It is a new systems biology approach that is frequently used to study the time course of cell differentiation and biological development (Li et al., 2017; Liu et al., 2017; Tao et al., 2018). WGCNA can be used to construct gene expression networks by clustering highly correlated genes into modules. These biologically relevant modules usually contain functionally related genes and may include key driver genes that can serve as potential diagnostic and prognostic biomarkers, or as therapeutic targets (Chen et al., 2016a; Chou et al., 2014).

In this study, we constructed a WGCNA network using expression data of MSCs during OS and adipogenic (AD) differentiation. Biological enrichment analysis on modules of interest and their corresponding hub genes identified potential key driver genes that could be used as biomarkers and therapeutic targets of osteoporosis.

Materials and Methods

Results

Co-expression network construction

A total of 66 samples and 19,610 genes from the GSE80614 data set were used for network construction. One of the most critical parameters in WGCNA network construction is the power value, which affects the independence and average connectivity of the co-expression modules. The degree of independence reached 0.95 and the average connectivity was highest using a power equal to 7 (Figs. 1 and 2). The weighted network was then constructed based on the scale-free topology criteria. Eventually, 22 modules were detected using the dynamic tree cutting approach (Fig. 3).

Gene co-expression modules were associated with clinical traits

After correlating the modules with clinical traits, a high correlation was observed with Phase I, Phase II and AD differentiation stages of mesenchymal cells (Fig. 4). The darkolivegreen module (r = 0.86, P = 2e⁻²⁰) had the highest correlation with Phase I, followed by the brown module (r = 0.78, p = 2e⁻¹⁴). The darkturquoise module (r = 0.81, p = 3e⁻¹⁶) had the highest correlation with Phase II of differentiation. On the other hand, the orangered4 module had the highest correlation (r = 0.84, p = 2e⁻¹⁸) with the AD phase. The inter- and intra-modular gene correlation of the four modules was then plotted based on GS and MM (Figs. 5A–5D). GO and KEGG enrichment analysis was then performed, which showed that 893 genes from the darkolivegreen module were mainly enriched for the p130Cas-ER-alpha-cSrc-kinase-PI3-kinase p85-subunit complex and tissue morphogenesis (Fig. 6A, P < 0.01). The 2,879 genes in the brown module were primarily enriched for Cell Cycle and Cell Cycle Checkpoints (Fig. 6B, P < 0.001). On the other hand, 437 genes in the darkturquoise module were associated with the UDP-glucuronate metabolic process and with the term “Metallothioneins bind metals” (Fig. 6C, P < 0.01). Finally, 47 genes of the orangered4 module were principally enriched for leukocyte apoptotic processes and the regulation of ion transmembrane transporter activity (Fig. 6D, P < 0.01).

Identification of module-related hub genes and lncRNAs

Numerous studies have shown that genes with higher MM and GS could be potential candidates for further research (Fuller et al., 2007; Oldham, Horvath & Geschwind, 2006; Wang et al., 2017). Hence, based on MM > 0.6 and GS > 0.6, 174 genes were identified in the darkolivegreen module, 608 genes and 15 lncRNAs (DLEU1, EPHA5-AS1, H19, LINC00284, LINC00839, LINC00921, LINC01119, LINC01213, LINC01291, LINC01616, PVT1, SCARNA9, SMAD5-AS1, STX18-AS1, TTTY15) in the brown module, 42 genes in the darkturquoise module, and 12 genes in the orangered4 module. The Metascape website was used to perform enrichment analysis on the brown module hub genes. The results showed that the brown module is mainly related to cell cycle regulation, which is in line with the original findings published using this data set (Van de Peppel et al., 2017). Besides, it has been shown that cell cycle regulation is associated with the differentiation of bone marrow MSCs (Batsali et al., 2017; Boucher et al., 2016). In the enrichment results, the column named “Description” represents detailed pathways and biological processes, while the corresponding pathway-related genes can be found in the column named “Hits” (S1). The results showed that many pathways and biological processes are related to mesenchymal cell differentiation, stem cell differentiation, and bone development. Cytoscape was then used to identify hub genes in the brown module using the MCODE plug-in, wherein the highest ranked module was selected, which had 24 nodes, 276 edges, and an average MCODE score of 23 (Fig. 7).

Gene expression analysis of Phase I group

The Phase I group included 30 bone marrow MSCs samples during OS and AD differentiation. The bone marrow MSCs were induced to differentiate for 0.5, 1, 2 and 3 h, respectively. Each time point was then compared with the 0 h samples to obtain DEGs. The volcano plots of OS and AD differentiation were created using the R software (Fig. 8).

Hierarchical clustering of DEGs

A total of 197 DEGs were screened, including 117 upregulated and 31 downregulated DEGs during AD differentiation, and 98 upregulated and 21 downregulated DEGs during OS differentiation (S2). Among these, in the 0.5 h vs 0 h comparison, eight genes were upregulated after 0.5 h of AD differentiation and 13 genes were upregulated after 0.5 h of OS differentiation. In the 1 h vs 0 h group, 37 upregulated and two downregulated genes were associated with AD differentiation, while 32 differentially upregulated genes were associated with OS differentiation. In the 2 h vs 0 h comparison, 72 and 13 genes were upregulated and downregulated, respectively, after 2 h of AD differentiation, while 66 upregulated and 9 downregulated genes were differentially expressed after 2 h of OS differentiation. In the 3 h vs 0 h comparison, 94 upregulated and 32 downregulated genes were associated with AD differentiation, while 74 upregulated and 19 downregulated genes were associated with OS differentiation. Hierarchical clustering of the DEGs during OS differentiation and AD differentiation was performed using the R software (Figs. 9 and 10).

Identification of genes involved in the differentiation of bone marrow MSCs

Enrichment analysis showed that the brown module was associated with the differentiation of bone marrow MSCs. Interestingly, a set of 13 candidate hub genes were common between the brown module genes and the DEGs from Phase I (Fig. 11).

Transcription factors and miRNAs screening

The transcription factor regulatory network was constructed and visualized using Cytoscape software. The network included 2 transcription factors, (ZNF740, FOS), and the above-mentioned 13 candidate hub genes (Fig. 12). ZNF740 and FOS were connected with 11 and 13 hub genes, respectively. Similarly, the miRNA regulatory network was also constructed for the hub genes using TargetScan, miRTarBase and miRDB databases (Fig. 13; Table 1).

Table 1:

The miRNA regulatory network was constructed using TargetScan, miRTarBase and miRDB databases.

A total of 110 miRNAs were identified using public miRNA databases.

Gene	miRNA
BHLHE40	hsa-miR-374a-5p	hsa-miR-8485	hsa-miR-4469
	hsa-miR-454-3p
FOXQ1	hsa-miR-378j	hsa-miR-6507-3p	hsa-miR-6839-5p
	hsa-miR-422a	hsa-miR-6780a-3p	hsa-miR-520h
	hsa-miR-1271-5p	hsa-miR-520g-3p	hsa-miR-506-3p
HBEGF	hsa-miR-4659a-3p	hsa-miR-6868-3p	hsa-miR-4659b-3p
	hsa-miR-30a-3p	hsa-miR-30e-3p	hsa-miR-4254
	hsa-miR-194-5p	hsa-miR-4778-3p	hsa-miR-6081
	hsa-miR-30d-3p	hsa-miR-132-3p
IL6	hsa-miR-149-5p	hsa-let-7f-5p	hsa-let-7a-5p
	hsa-let-7c-5p	hsa-miR-98-5p
LIF	hsa-miR-3180-5p	hsa-miR-624-5p	hsa-miR-4286
	hsa-miR-660-3p	hsa-miR-6778-3p	hsa-miR-5193
	hsa-miR-3922-5p	hsa-miR-6873-3p
PTGS2	hsa-miR-146a-5p	hsa-miR-26b-5p	hsa-miR-26a-5p
	hsa-miR-132-3p
SERPINE1	hsa-miR-30c-5p	hsa-miR-3145-5p	hsa-miR-148a-3p
	hsa-miR-301a-3p	hsa-miR-519d-5p	hsa-miR-143-3p
	hsa-miR-30b-5p	hsa-miR-486-5p	hsa-miR-192-5p
	hsa-miR-145-5p
SGK1	hsa-miR-576-3p	hsa-miR-2115-5p	hsa-miR-5682
	hsa-miR-5692a	hsa-miR-4641	hsa-miR-125a-3p
	hsa-miR-3662	hsa-miR-548aw	hsa-miR-548x-5p
	hsa-miR-1468-3p	hsa-miR-133b	hsa-miR-3202
	hsa-miR-19a-3p	hsa-miR-4753-3p	hsa-miR-29c-3p
	hsa-miR-19b-3p	hsa-miR-6079	hsa-miR-548aj-5p
	hsa-miR-6871-3p	hsa-miR-4639-5p	hsa-miR-29b-3p
	hsa-miR-520a-5p	hsa-miR-4533	hsa-miR-548f-5p
	hsa-miR-525-5p	hsa-miR-5571-5p	hsa-miR-29a-3p
	hsa-miR-365b-3p	hsa-miR-548g-5p	hsa-miR-552-3p
	hsa-miR-365a-3p
SPRY2	hsa-miR-23a-3p	hsa-miR-122-5p	hsa-miR-21-5p
	hsa-miR-27a-3p
ZNF281	hsa-miR-518a-3p	hsa-miR-6499-3p	hsa-miR-23a-3p
	hsa-miR-3140-3p	hsa-miR-6868-3p	hsa-miR-1248
	hsa-miR-4430	hsa-miR-3910	hsa-miR-6856-3p
	hsa-miR-4761-3p	hsa-miR-6857-5p	hsa-miR-5681a
	hsa-miR-518b	hsa-miR-495-5p	hsa-miR-23c
	hsa-miR-33a-5p	hsa-miR-23b-3p	hsa-miR-5197-5p
	hsa-miR-518f-3p	hsa-miR-451b	hsa-miR-7844-5p
	hsa-miR-4517	hsa-miR-488-3p	hsa-miR-518d-3p
	hsa-miR-518c-3p	hsa-miR-3652	hsa-miR-33b-5P

DOI: 10.7717/peerj.8907/table-1

Note:

DEGs, differentially expressed genes; miRNA, microRNAs.

GSEA analysis on the candidate genes

GSEA was applied to reveal the potential role of candidate genes (FOXQ1, SGK1) during MSC differentiation. Our results showed that the JAK-STAT signaling pathway was positively correlated with FOXQ1 expression (NES > 1.5, NOM P-val (nominal P value) < 0.05 and FDR < 0.25) (Fig. 14A). On the other hand, three gene sets were positively correlated with SGK1 expression, including the “RAS”, the “JAK-STAT signaling pathway” and the “Phosphatidylinositol signaling system” (Figs. 14B–14D).

Discussion

Osteoporosis can increase fracture risk, which is a growing concern as the population ages. MSCs are mesoderm-derived multipotent cells that can differentiate into adipocytes, osteoblasts and chondrocytes. Transformation of MSCs into adipocytes, rather than osteoblasts, has been associated with the development of osteoporosis. Indeed, the microenvironment of the bone marrow cavity changes with age, which can lead to increase in fat cells formation and inhibition of bone formation. During this process, adipose tissue gradually replaces bone tissue, leading to osteoporosis. Therefore, understanding the molecular mechanisms of osteogenesis and MSC differentiation is particularly imperative. In this study, we used WGCNA to extract gene modules that are associated with OS and AD differentiation of bone marrow MSCs, and to identify potential biomarkers and therapeutic targets of osteoporosis.

First, we constructed a co-expression network containing 22 modules, which were then compared with DEGs during Phase I (0–3 h), Phase II (6–24 h), Phase III (48–96 h), AD and OS differentiation of MSCs. The brown module was highly enriched in genes involved in cell cycle-related pathways and other processes related to MSC differentiation, such as mitotic nuclear division (Dudakovic et al., 2014), DNA repair (Dudakovic et al., 2014), and chromosome maintenance (Oktar et al., 2011). Interestingly, the enrichment results showed that numerous genes were related to osteoblast differentiation, bone mineralization, bone remodeling, bone resorption, and skeletal system development. Then, we analyzed 608 genes that were selected based on the parameters of MM and GS to obtain hub genes with the highest scores (MM > 0.85, GS > 0.85). Among these genes, the NRG1 gene has been shown to mediate cell–cell signaling and play a critical role in the growth and development of multiple organ systems. Besides, NRG1 was shown to be upregulated by Wnt3a during Wnt3a-induced osteoblast differentiation in primary human MSCs (Jullien et al., 2012). On the other hand, ITGA2 was shown to induce the expression of WNT5A to promote OS differentiation of human MSCs (Olivares-Navarrete et al., 2011). SERPINB2 is a TGF-β-responsive gene that plays a negative regulatory role in the differentiation of human bone stromal cells (Elsafadi et al., 2017). PRKCA, another hub gene, was demonstrated to regulate bone architecture and osteoblast activity (Galea et al., 2014). FAM20C encodes a family member of secreted protein kinases, which can bind calcium and phosphorylate proteins involved in bone mineralization (Liu et al., 2018). RBPJ, a primary nuclear mediator of Notch signaling, has also been associated with osteogenesis (S3) (Wang et al., 2010).

In order to screen for candidate genes that can serve as potential biomarkers or therapeutic targets, the Phase I samples, associated with the brown module, were compared with the control samples, which led to the identification of 119 DEGs related to OS differentiation and 148 DEGs that were associated with AD differentiation. The identified DEGs were then compared with the brown module, which led to the identification of 13 common genes, including BHLHE40, FOSB (Schneider et al., 2017), CSRNP1, HBEGF (Li et al., 2019a), CLCF1 (Nahle et al., 2019), SPRY2 (Schneider et al., 2017), IL6 (Munir et al., 2017), PTGS2 (Feigenson et al., 2017), SGK1, FOXQ1, LIF (Huat et al., 2014), SERPINE1 (Karbiener et al., 2011) and ZNF281. Of interest, the previously identified OS differentiation-related genes FOSB, SPRY2 and PTGS2 were upregulated, while HBEGF, CLCF1 and LIF were downregulated. In addition, IL6 and SERPINE1, involved in AD differentiation of bone marrow MSCs, were downregulated. Among the identified candidate genes, BHLHE40, CSRNP1, SGK1, FOXQ1 and ZNF281 could be novel potential targets related to OS and AD MSC differentiation. Moreover, 110 miRNAs and 2 transcription factors (ZNF740, FOS) were identified using public databases as candidate upstream regulators of the network’s hub genes. We found that both ZNF740 and FOS genes were differentially expressed during Phase I of differentiation, while only FOS was upregulated during OS and AD differentiation. Besides, FOS was also a member of the brown module and was differentially expressed in the validation data set. The transcription factor FOS belongs to the FOS gene family, which has been described as a regulator of cell proliferation, differentiation and transformation. On the other hand, another data set of MSC differentiation was used to validate four candidate genes with significant differential expression (IL6, HBEGF, FOXQ1 and SGK1). Finally, GSEA was performed to identify candidate pathways associated with FOXQ1 and SGK1 expression, which showed that FOXQ1 is associated with JAK-STAT signaling pathway, while SGK1 is mainly associated with the RAS. Studies have shown that inhibition of the JAK-STAT signaling pathway promoted OS differentiation (Levy et al., 2010). The RAS regulated ossification and angiogenesis of bones, as well as OS differentiation (Durik, Sevá Pessôa & Roks, 2012). On the other hand, the phosphatidylinositol signaling system could regulate the OS and AD differentiation of bone marrow MSCs (Song et al., 2017).

SGK1 is a member of the PI3-kinase pathway, which regulates cell survival, proliferation and differentiation (Lang et al., 2006). In addition, we observed that genes in the darkolivegreen module were associated with cell cycle arrest, cell proliferation and apoptotic processes, but also were related to PI3-kinase and MAPK signaling pathway that could also participate in osteoblast differentiation. In fact, the PI3-kinase signaling pathway is downstream to the JAK-STAT and RAS signaling pathways. Therefore, we may infer that the brown and the darkolivegreen modules could be complementary, or common, with different co-expression patterns. Among the genes associated with the PI3K pathway in the darkolivegreen module, ESR1 encodes the estrogen receptor ERα that regulates the classic estrogen signaling pathway. Studies have shown that ERα deficiency inhibited OS differentiation and promoted adipocytic differentiation of bone MSCs (Li et al., 2019b). In breast cancer, ERα expression was associated with mTORC1-mediated phosphorylation of SGK1-Ser422 (Hall et al., 2012). These results indicate a possible interaction between ESR1 and SGK1 during MSC differentiation. In addition, SRC, a non-receptor tyrosine kinase, is a key signaling molecule in bone metabolism that enhances OS differentiation through phosphorylation of Osterix (Choi et al., 2015). Studies have shown that in some human diseases, SGK1 was the key mediator of SRC-induced transformation (Ma et al., 2019). So, we may assume that SRC can regulate the expression of SGK1 during MSC differentiation. In conclusion, our study suggests that FOXQ1 and SGK1 may regulate the OS and AD differentiation of bone marrow MSCs through these signaling pathways, which requires further experimental validation in the future.

Conclusion

In summary, WGCNA and DEGs analysis were used in this study, which showed that correlated genes in the brown module could be critical for the differentiation of bone marrow MSCs. A total of two transcription factors and two hub genes, which are potential targets of 110 miRNAs, were associated with OS and AD differentiation of MSCs. These key drivers may serve as potential biomarkers for the diagnosis and prognosis of osteoporosis, which need to be validated by more studies in the future.

Supplemental Information