Haihong Zhang1Jinghui Feng2Zhiguo Lin1
Shuya Wang1Yan Wang1Siming Dai1Weisi Kong1
Yanli Wang1Zhiyi Zhang1*- 1Department of Rheumatology and Immunology, The First Affiliated Hospital of Harbin Medical University, Harbin, China
- 2Department of Gerontology, The First Affiliated Hospital of Harbin Medical University, Harbin, China
Osteoporosis is a kind of brittle bone disease, which is characterized by a reduction in bone mineral density (BMD). In recent years, a number of genes and pathophysiological mechanisms have been identified for osteoporosis. However, the genes associated with BMD remain to be explored. Toward this end, we integrated multiple osteoporosis microarray datasets to identify and systematically characterize BMD-related genes. By integrating the differentially expressed genes from three osteoporosis microarray datasets, 152 genes show differentially expressed between high and low BMD osteoporosis samples in at least two of the three datasets. Among them, 88 were up-regulated in high BMD samples and 64 were up-regulated in low BMD samples. The expression of ZFP36, JUNB and TMEM8A were increased at high BMD samples in all three datasets. Hub genes were further identified by co-expression network analysis. Functional enrichment analysis showed that the gene up-regulated in high BMD were enriched in immune-related functions, suggesting that the immune system plays an important role in osteoporosis. Our study explored BMD-related genes based on the integration of osteoporosis microarray data, providing guidance to other researchers from a new perspective.
Introduction
Osteoporosis, a common systemic bone disease, can lead to weak bones and increase the risk of fractures (Ensrud and Crandall, 2017; Wilson et al., 2017; Khosla et al., 2018). A third of women and a fifth of men over the age of 50 have broken bones due to osteoporosis (Brown, 2017). Osteoporosis is characterized by low bone mass, microstructure degeneration and reduced bone strength. Patients often have a reduction in bone mineral density (BMD) (Coughlan and Dockery, 2014; Black and Rosen, 2016). Bone homeostasis depends on osteoclast absorption and osteoblast formation. The imbalance of this tightly coupled process can lead to the development of osteoporosis (Chen et al., 2018b). This process involves changes in a variety of signaling pathways, including MAPK signaling pathway, NF-κB pathway, Notch signaling pathway, etc. (Chen et al., 2018a, 2019; Lee and Long, 2018; Wang et al., 2019). Notch signaling pathway regulates the differentiation and function of osteoblasts and osteoclasts and participates in the process of bone reconstruction, activation of notch signaling pathway can inhibit glucose metabolism and osteoblast differentiation of bone marrow mesenchymal progenitor cells (Lee and Long, 2018).
In recent years, some genes and pathophysiological mechanisms have been identified by microarray analysis for patients with osteoporosis (Liu et al., 2005; Lei et al., 2009; Li et al., 2016). The study based on peripheral blood monocyte cells (PBMCs) microarray data of osteoporosis patients revealed the pathophysiological mechanism of osteoporosis, which is characterized by increased recruitment of monocytes into bone and then differentiating into osteoclasts (Liu et al., 2005). Zhou et al. (2018b, 2019) predicted osteoporosis related transcription factors (TF) and long non-coding RNA (lncRNA) via exon arrays. Another study explored osteoporosis-related pathways based on microarray data (Zhou et al., 2018a). The studies mentioned above were based on a limited set of data, and have certain limitations. In order to obtain more robust results, it is critical to integrate multiple datasets for obtaining new insights.
In this study, we collected three microarray datasets of osteoporosis. Differential expression analysis was conducted to identify differentially expressed genes (DEGs) between high and low BMD samples. Then DEGs were integrated to obtain uniformly expressed BMD-related genes. Analysis of gene co-expression networks revealed key genes of different BMD conditions. And enrichment analysis showed that they were enriched in immune-related functions and biological pathways, suggesting a potential role of the immune system in osteoporosis.
Materials and Methods
GEO Datasets
Microarray datasets of high and low BMD osteoporosis samples were downloaded from Gene Expression Omnibus (GEO) database1, including GSE2208, GSE56814 and GSE56815. A total of 172 samples were collected from the three datasets, of which 92 were high BMD samples and 80 were low BMD samples (Table 1).
Table 1. The overview of three GEO datasets.
Identification of Differentially Expressed Genes (DEGs)
The up-regulated genes in high or low BMD samples of each dataset were identified using the R package named limma with a threshold of |log2FoldChange|>0 and P < 0.05. Genes up-regulated in at least two datasets were identified as uniformly expressed BMD-related genes and used for subsequent analysis.
The Construction of Gene Co-expression Network
We calculated the Pearson Correlation Coefficient (PCC) of integrated up-regulated genes in high or low BMD samples (PCC > 0.4 and P < 0.05). For co-expressed gene pairs that appeared in multiple datasets, the averaged PCC was calculated. Then Cytoscape (V.3.8.0) was used for network visualization, and hub genes were identified by cytoHubba plugin (Shannon et al., 2003; Chin et al., 2014).
Enrichment Analysis
We performed GO and KEGG enrichment analysis of genes included in co-expression network using R package “clusterProfiler” (Yu et al., 2012). The significance threshold is 0.05.
Results
Differentially Expressed Genes Between High and Low BMD Osteoporosis Samples
To explore the key genes associated with BMD, we downloaded three osteoporosis microarray datasets from NCBI GEO database. In total, 19718 genes and 172 samples were involved. The number of high and low BMD samples in each dataset was comparable (Table 1).
We performed the analysis of differentially expressed genes between the two types of samples (Figure 1). There were 12.59% (GSE2208), 3.30% (GSE56814), and 18.49% (GSE56815) genes identified as differentially expressed, respectively. In GSE56814, compared to the low BMD samples, there were more genes up-regulated in high BMD samples, while GSE56814 had more genes up-regulated in low BMD samples. In GSE2208, the two types of genes were almost the same in number (Figure 1D).
Figure 1. Differential expression analysis. (A–C) Volcano plot of DEGs in each dataset, red nodes represent upregulation in high BMD samples and blue nodes represent upregulation in low BMD samples. (D) Statistics of two types of DEGs.
Integration of BMD-Related Genes
We integrated the DEGs of three datasets to obtain robust different BMD-related genes. Based on the criterion that the same gene was up-regulated in at least two datasets, we identified 88 and 64 uniformly expressed genes with high and low BMD by UpSetR, respectively (Figures 2A,B and Supplementary Tables 1, 2; Conway et al., 2017). In all datasets, there were three genes, ZFP36, JUNB and TMEM8A shown up-regulated in high BMD (Figure 2C and Table 2). JUNB is a member of c-Jun protein family, which interacts with c-Fos protein family to form transcription factor AP-1. AP-1 can activate osteoclast specific genes (Wagner and Eferl, 2005; Asagiri and Takayanagi, 2007; Hamamura et al., 2015). Only one gene, CACNA2D3, had elevated expression with low BMD in all datasets (Figure 2C and Table 2).
Figure 2. Integration of BMD-related DEGs. (A) UpSetR plot of genes up-regulated in high BMD samples, the number of genes present in at least two datasets is marked (red). (B) UpSetR plot of genes up-regulated in low BMD samples, the number of genes present in at least two datasets is marked (red). (C) Boxplot of genes that expressed consistently across all three datasets.
Table 2. Genes expressed consistently across all three sets of data.
Co-expression Network and Functional Analysis of BMD-Related Genes
Through the similarity of gene expression, the possible interaction between gene products can be analyzed, so as to understand the interaction between genes and find out the core genes. Therefore, we performed a co-expression analysis of integrated up-regulated genes in high and low BMD samples, respectively. We obtained 58 co-expressed gene pairs of high BMD. The top5 hub genes were identified by cytoHubba, a plug-in of Cytoscape, including KDM2A, APH1A, DNPEP, NFKBIB, and TMEM8A. KDM2A was co-expressed with the largest number of genes (Figure 3A). KDM2A can regulate Mesenchymal stem cells (MSCs) osteo/dentinogenic differentiation and cell proliferation (Dong et al., 2013; Gao et al., 2013). Co-expression network of genes that were up-regulated in high BMD samples contained 48 gene nodes. We did functional enrichment analysis on them. In addition to osteoporosis-related functions and pathways, Notch signaling pathway and osteoclast differentiation, enrichment analysis results showed that they were mainly enriched in immune-related functions and pathways (Figure 3B; Boyle et al., 2003; Lee and Long, 2018). Studies have shown that the immune system plays an important role in osteoporosis (Faienza et al., 2013). For example, postmenopausal osteoporosis patients had higher T-cell activity and increased TNFα and RANKL production, which can promote osteoclast differentiation (D’Amelio et al., 2008; Mirza et al., 2010; Kim et al., 2012). In addition, other studies have shown that B lymphocytes were closely related to bone metabolism (Takayanagi et al., 2005; Boyce and Xing, 2008; Breuil et al., 2010).
Figure 3. Co-expression network and function enrichment analysis of genes that are consistently expressed in high BMD samples. (A) The gene co-expression network. Node size represents how many other genes interact with it. Node color indicates in which data sets it is up-regulated in high BMD samples. The darker the edge color, the greater the correlation coefficient. (B) GO and KEGG results of genes in network.
In the co-expression network of genes that were up-regulated in low BMD samples, there were 44 gene pairs, involving 38 genes (Figure 4A). The top5 hub genes were CCT7, DGUOK, MPHOSPH10, RARS, and DMTF1. The results of enrichment analysis showed that low BMD-related genes were mainly enriched in basic biological processes and pathways, including ribosome biogenesis, rRNA processing and translation elongation (Figure 4B).
Figure 4. Co-expression network and function enrichment analysis of genes that are consistently expressed in low BMD samples. (A) The gene co-expression network. Node size represents how many other genes interact with it. Node color indicates in which datasets it is up-regulated in low BMD samples. The darker the edge color, the greater the correlation coefficient. (B) GO and KEGG results of genes in network.
Discussion
In this study, we identified and analyzed BMD-related genes through three osteoporosis microarray datasets. Differential expression analysis identified DEGs between high and low BMD in each dataset. By integration, we screened out 152 uniformly differentially expressed genes. Gene co-expression network analysis further identified key top 5 hub genes. Co-expression genes that had elevated expression in high BMD were enriched in functions of immune systems, suggesting its important potential role in osteoporosis. In addition, they were also enriched in osteoclast differentiation pathway, indicating that high BMD is developing toward low BMD.
In conclusion, we used bioinformatic methods to systematically characterize genes underlying BMD levels based on osteoporosis microarray data. We analyzed osteoporosis from a new perspective and provided new guidance for its diagnosis and treatment.
Data Availability Statement
All datasets presented in this study are included in the article/Supplementary Material.
Author Contributions
ZZ designed the study. HZ, JF, ZL, SW, and YW analyzed the data. SD, WK, and YLW wrote the manuscript. All authors read and approved the manuscript.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.00798/full#supplementary-material
TABLE S1 | The uniformly up-regulated DEGs in high BMD samples.
TABLE S2 | The uniformly up-regulated DEGs in low BMD samples.
Footnotes
References
Asagiri, M., and Takayanagi, H. (2007). The molecular understanding of osteoclast differentiation. Bone 40, 251–264. doi: 10.1016/j.bone.2006.09.023
Black, D. M., and Rosen, C. J. (2016). Clinical practice. postmenopausal osteoporosis. N. Engl. J. Med. 374, 254–262.
Boyce, B. F., and Xing, L. (2008). Bruton and tec: new links in osteoimmunology. Cell Metab 7, 283–285. doi: 10.1016/j.cmet.2008.03.013
Boyle, W. J., Simonet, W. S., and Lacey, D. L. (2003). Osteoclast differentiation and activation. Nature 423, 337–342. doi: 10.1038/nature01658
Breuil, V., Ticchioni, M., Testa, J., Roux, C. H., Ferrari, P., Breittmayer, J. P., et al. (2010). Immune changes in post-menopausal osteoporosis: the Immunos study. Osteoporos. Int. 21, 805–814. doi: 10.1007/s00198-009-1018-7
Chen, E., Liu, G., Zhou, X., Zhang, W., Wang, C., Hu, D., et al. (2018a). Concentration-dependent, dual roles of IL-10 in the osteogenesis of human BMSCs via P38/MAPK and NF-kappaB signaling pathways. FASEB J. 32, 4917–4929. doi: 10.1096/fj.201701256rrr
Chen, X., Wang, Z., Duan, N., Zhu, G., Schwarz, E. M., and Xie, C. (2018b). Osteoblast-osteoclast interactions. Connect. Tissue Res. 59, 99–107. doi: 10.1080/03008207.2017.1290085
Chen, K., Qiu, P., Yuan, Y., Zheng, L., He, J., Wang, C., et al. (2019). Pseurotin a inhibits osteoclastogenesis and prevents ovariectomized-induced bone loss by suppressing reactive oxygen species. Theranostics 9, 1634–1650. doi: 10.7150/thno.30206
Chin, C. H., Chen, S. H., Wu, H. H., Ho, C. W., Ko, M. T., and Lin, C. Y. (2014). cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst. Biol. 8(Suppl. 4):S11. doi: 10.1186/1752-0509-8-S4-S11
Conway, J. R., Lex, A., and Gehlenborg, N. (2017). UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics 33, 2938–2940. doi: 10.1093/bioinformatics/btx364
Coughlan, T., and Dockery, F. (2014). Osteoporosis and fracture risk in older people. Clin. Med. (Lond) 14, 187–191. doi: 10.7861/clinmedicine.14-2-187
D’Amelio, P., Grimaldi, A., Di Bella, S., Brianza, S. Z. M., Cristofaro, M. A., Tamone, C., et al. (2008). Estrogen deficiency increases osteoclastogenesis up-regulating T cells activity: a key mechanism in osteoporosis. Bone 43, 92–100. doi: 10.1016/j.bone.2008.02.017
Dong, R., Yao, R., Du, J., Wang, S., and Fan, Z. (2013). Depletion of histone demethylase KDM2A enhanced the adipogenic and chondrogenic differentiation potentials of stem cells from apical papilla. Exp. Cell Res. 319, 2874–2882. doi: 10.1016/j.yexcr.2013.07.008
Faienza, M. F., Ventura, A., Marzano, F., and Cavallo, L. (2013). Postmenopausal osteoporosis: the role of immune system cells. Clin. Dev. Immunol. 2013:575936.
Gao, R., Dong, R., Du, J., Ma, P., Wang, S., and Fan, Z. (2013). Depletion of histone demethylase KDM2A inhibited cell proliferation of stem cells from apical papilla by de-repression of p15INK4B and p27Kip1. Mol. Cell. Biochem. 379, 115–122. doi: 10.1007/s11010-013-1633-7
Hamamura, K., Chen, A., Tanjung, N., Takigawa, S., Sudo, A., and Yokota, H. (2015). In vitro and in silico analysis of an inhibitory mechanism of osteoclastogenesis by salubrinal and guanabenz. Cell Signal 27, 353–362. doi: 10.1016/j.cellsig.2014.11.020
Khosla, S., Farr, J. N., and Kirkland, J. L. (2018). Inhibiting cellular senescence: a new therapeutic paradigm for age-related osteoporosis. J. Clin. Endocrinol. Metab. 103, 1282–1290. doi: 10.1210/jc.2017-02694
Kim, B. J., Bae, S. J., Lee, S. Y., Lee, Y. S., Baek, J. E., Park, S. Y., et al. (2012). TNF-alpha mediates the stimulation of sclerostin expression in an estrogen-deficient condition. Biochem. Biophys. Res. Commun. 424, 170–175. doi: 10.1016/j.bbrc.2012.06.100
Lee, S. Y., and Long, F. (2018). Notch signaling suppresses glucose metabolism in mesenchymal progenitors to restrict osteoblast differentiation. J. Clin. Invest. 128, 5573–5586. doi: 10.1172/jci96221
Lei, S. F., Wu, S., Li, L. M., Deng, F. Y., Xiao, S. M., Jiang, C., et al. (2009). An in vivo genome wide gene expression study of circulating monocytes suggested GBP1, STAT1 and CXCL10 as novel risk genes for the differentiation of peak bone mass. Bone 44, 1010–1014. doi: 10.1016/j.bone.2008.05.016
Li, J. J., Wang, B. Q., Fei, Q., Yang, Y., and Li, D. (2016). Identification of candidate genes in osteoporosis by integrated microarray analysis. Bone Joint Res. 5, 594–601. doi: 10.1302/2046-3758.512.bjr-2016-0073.r1
Liu, Y. Z., Dvornyk, V., Lu, Y., Shen, H., Lappe, J. M., Recker, R. R., et al. (2005). A novel pathophysiological mechanism for osteoporosis suggested by an in vivo gene expression study of circulating monocytes. J. Biol. Chem. 280, 29011–29016. doi: 10.1074/jbc.m501164200
Mirza, F. S., Padhi, I. D., Raisz, L. G., and Lorenzo, J. A. (2010). Serum sclerostin levels negatively correlate with parathyroid hormone levels and free estrogen index in postmenopausal women. J. Clin. Endocrinol. Metab. 95, 1991–1997. doi: 10.1210/jc.2009-2283
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Takayanagi, H., Sato, K., Takaoka, A., and Taniguchi, T. (2005). Interplay between interferon and other cytokine systems in bone metabolism. Immunol. Rev. 208, 181–193. doi: 10.1111/j.0105-2896.2005.00337.x
Wagner, E. F., and Eferl, R. (2005). Fos/AP-1 proteins in bone and the immune system. Immunol. Rev. 208, 126–140.
Wang, F., Zhang, C., Ge, W., and Zhang, G. (2019). Up-regulated CST5 inhibits bone resorption and activation of osteoclasts in rat models of osteoporosis via suppression of the NF-kappaB pathway. J. Cell Mol. Med. 23, 6744–6754.
Wilson, L. M., Rebholz, C. M., Jirru, E., Liu, M. C., Zhang, A., Gayleard, J., et al. (2017). Benefits and harms of osteoporosis medications in patients with chronic kidney disease: a systematic review and meta-analysis. Ann. Intern. Med. 166, 649–658.
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287.
Zhou, Y., Gao, Y., Xu, C., Shen, H., Tian, Q., and Deng, H. W. (2018a). A novel approach for correction of crosstalk effects in pathway analysis and its application in osteoporosis research. Sci. Rep. 8:668.
Zhou, Y., Zhu, W., Zhang, L., Zeng, Y., Xu, C., Tian, Q., et al. (2018b). Transcriptomic data identified key transcription factors for osteoporosis in caucasian women. Calcif. Tissue Int. 103, 581–588.
Keywords: bone mineral density, osteoporosis, microarray, co-expression, enrichment analysis
Citation: Zhang H, Feng J, Lin Z, Wang S, Wang Y, Dai S, Kong W, Wang Y and Zhang Z (2020) Identification and Analysis of Genes Underlying Bone Mineral Density by Integrating Microarray Data of Osteoporosis. Front. Cell Dev. Biol. 8:798. doi: 10.3389/fcell.2020.00798
Received: 02 July 2020; Accepted: 28 July 2020;
Published: 27 August 2020.
Edited by:
Lei Deng, Central South University, ChinaReviewed by:
Hao Lin, University of Electronic Science and Technology of China, ChinaLiang Yu, Xidian University, China
Copyright © 2020 Zhang, Feng, Lin, Wang, Wang, Dai, Kong, Wang and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhiyi Zhang, zhangzhiyi2014@163.com