Front. Cell Dev. Biol. Frontiers in Cell and Developmental Biology Front. Cell Dev. Biol. 2296-634X Frontiers Media S.A. 10.3389/fcell.2020.00798 Cell and Developmental Biology Original Research Identification and Analysis of Genes Underlying Bone Mineral Density by Integrating Microarray Data of Osteoporosis Zhang Haihong 1 Feng Jinghui 2 Lin Zhiguo 1 Wang Shuya 1 Wang Yan 1 Dai Siming 1 Kong Weisi 1 Wang Yanli 1 Zhang Zhiyi 1 * 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

Edited by: Lei Deng, Central South University, China

Reviewed by: Hao Lin, University of Electronic Science and Technology of China, China; Liang Yu, Xidian University, China

 *Correspondence: Zhiyi Zhang, zhangzhiyi2014@163.com

This article was submitted to Molecular Medicine, a section of the journal Frontiers in Cell and Developmental Biology

 27 08 2020 2020 8 798 02 07 2020 28 07 2020 Copyright © 2020 Zhang, Feng, Lin, Wang, Wang, Dai, Kong, Wang and Zhang. 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.

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.

 bone mineral density osteoporosis microarray co-expression enrichment analysis Health and Family Planning Commission of Heilongjiang Province10.13039/501100008978 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).

The overview of three GEO datasets.

Dataset Sample number
 Gene number
High BMD Low BMD
GSE2208 10 9 6353
GSE56814 42 31 17321
GSE56815 40 40 13515
Total 92 80 19718
The total number of gene does not include genes that appear multiple times. 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).

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).

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.

Genes expressed consistently across all three sets of data.

Gene GSE2208 GSE56814 GSE56815
log2FC P log2FC P log2FC P
ZFP36 0.355 0.024 0.216 0.042 0.274 0.004
JUNB 0.512 0.009 0.249 0.049 0.305 0.006
TMEM8A 0.279 0.024 0.134 0.021 0.211 0.013
CACNA2D3 −0.256 0.015 −0.111 0.015 −0.268 0.001
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).

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).

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

The uniformly up-regulated DEGs in high BMD samples.

The uniformly up-regulated DEGs in low BMD samples.

 References Asagiri M. Takayanagi H. (2007). The molecular understanding of osteoclast differentiation. Bone 40 251–264. 10.1016/j.bone.2006.09.023 17098490 Black D. M. Rosen C. J. (2016). Clinical practice. postmenopausal osteoporosis. N. Engl. J. Med. 374 254–262. Boyce B. F. Xing L. (2008). Bruton and tec: new links in osteoimmunology. Cell Metab 7 283–285. 10.1016/j.cmet.2008.03.013 18396132 Boyle W. J. Simonet W. S. Lacey D. L. (2003). Osteoclast differentiation and activation. Nature 423 337–342. 10.1038/nature01658 12748652 Breuil V. Ticchioni M. Testa J. Roux C. H. Ferrari P. Breittmayer J. P. (2010). Immune changes in post-menopausal osteoporosis: the Immunos study. Osteoporos. Int. 21 805–814. 10.1007/s00198-009-1018-7 19876583 Brown C. (2017). Osteoporosis: staying strong. Nature 550 S15–S17. Chen E. Liu G. Zhou X. Zhang W. Wang C. Hu D. (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. 10.1096/fj.201701256rrr 29630408 Chen X. Wang Z. Duan N. Zhu G. Schwarz E. M. Xie C. (2018b). Osteoblast-osteoclast interactions. Connect. Tissue Res. 59 99–107. 10.1080/03008207.2017.1290085 28324674 Chen K. Qiu P. Yuan Y. Zheng L. He J. Wang C. (2019). Pseurotin a inhibits osteoclastogenesis and prevents ovariectomized-induced bone loss by suppressing reactive oxygen species. Theranostics 9 1634–1650. 10.7150/thno.30206 31037128 Chin C. H. Chen S. H. Wu H. H. Ho C. W. Ko M. T. Lin C. Y. (2014). cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst. Biol. 8(Suppl. 4):S11. 10.1186/1752-0509-8-S4-S11 25521941 Conway J. R. Lex A. Gehlenborg N. (2017). UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics 33 2938–2940. 10.1093/bioinformatics/btx364 28645171 Coughlan T. Dockery F. (2014). Osteoporosis and fracture risk in older people. Clin. Med. (Lond) 14 187–191. 10.7861/clinmedicine.14-2-187 24715132 D’Amelio P. Grimaldi A. Di Bella S. Brianza S. Z. M. Cristofaro M. A. Tamone C. (2008). Estrogen deficiency increases osteoclastogenesis up-regulating T cells activity: a key mechanism in osteoporosis. Bone 43 92–100. 10.1016/j.bone.2008.02.017 18407820 Dong R. Yao R. Du J. Wang S. 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. 10.1016/j.yexcr.2013.07.008 23872478 Ensrud K. E. Crandall C. J. (2017). Osteoporosis. Ann. Intern. Med. 167 ITC17–ITC32. Faienza M. F. Ventura A. Marzano F. 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. 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. 10.1007/s11010-013-1633-7 23559091 Hamamura K. Chen A. Tanjung N. Takigawa S. Sudo A. Yokota H. (2015). In vitro and in silico analysis of an inhibitory mechanism of osteoclastogenesis by salubrinal and guanabenz. Cell Signal 27 353–362. 10.1016/j.cellsig.2014.11.020 25435425 Khosla S. Farr J. N. Kirkland J. L. (2018). Inhibiting cellular senescence: a new therapeutic paradigm for age-related osteoporosis. J. Clin. Endocrinol. Metab. 103 1282–1290. 10.1210/jc.2017-02694 29425296 Kim B. J. Bae S. J. Lee S. Y. Lee Y. S. Baek J. E. Park S. Y. (2012). TNF-alpha mediates the stimulation of sclerostin expression in an estrogen-deficient condition. Biochem. Biophys. Res. Commun. 424 170–175. 10.1016/j.bbrc.2012.06.100 22735261 Lee S. Y. Long F. (2018). Notch signaling suppresses glucose metabolism in mesenchymal progenitors to restrict osteoblast differentiation. J. Clin. Invest. 128 5573–5586. 10.1172/jci96221 30284985 Lei S. F. Wu S. Li L. M. Deng F. Y. Xiao S. M. Jiang C. (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. 10.1016/j.bone.2008.05.016 19223260 Li J. J. Wang B. Q. Fei Q. Yang Y. Li D. (2016). Identification of candidate genes in osteoporosis by integrated microarray analysis. Bone Joint Res. 5 594–601. 10.1302/2046-3758.512.bjr-2016-0073.r1 27908864 Liu Y. Z. Dvornyk V. Lu Y. Shen H. Lappe J. M. Recker R. R. (2005). A novel pathophysiological mechanism for osteoporosis suggested by an in vivo gene expression study of circulating monocytes. J. Biol. Chem. 280 29011–29016. 10.1074/jbc.m501164200 15965235 Mirza F. S. Padhi I. D. Raisz L. G. 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. 10.1210/jc.2009-2283 20156921 Shannon P. Markiel A. Ozier O. Baliga N. S. Wang J. T. Ramage D. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13 2498–2504. 10.1101/gr.1239303 14597658 Takayanagi H. Sato K. Takaoka A. Taniguchi T. (2005). Interplay between interferon and other cytokine systems in bone metabolism. Immunol. Rev. 208 181–193. 10.1111/j.0105-2896.2005.00337.x 16313349 Wagner E. F. Eferl R. (2005). Fos/AP-1 proteins in bone and the immune system. Immunol. Rev. 208 126–140. Wang F. Zhang C. Ge W. 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. (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. 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. 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. (2018b). Transcriptomic data identified key transcription factors for osteoporosis in caucasian women. Calcif. Tissue Int. 103 581–588. Zhou Y. Xu C. Zhu W. He H. Zhang L. Tang B. (2019). Long noncoding RNA analyses for osteoporosis risk in caucasian women. Calcif. Tissue Int. 105 183–192.

https://www.ncbi.nlm.nih.gov/geo/