This article is mentioned in:

Introduction

Osteoporosis is the most common and serious skeletal disorder among the elderly. Symptomatic osteoporosis occurs due to a decreased bone mineral density (BMD) leading to reduced bone strength and an increased risk of fractures (1). Low bone mass in the elderly is highly dependent on their peak bone mass (PBM) as young adults (2). Therefore, it is necessary to understand and identify the risk factors for impaired PBM in young and middle-aged adults.

Osteopenia may result from an imbalance between increased bone resorption and decreased bone formation (3,4). Bone resorption involves the dissolution of bone mineral and degradation of the organic bone matrix. These two functions are performed by osteoclasts. Osteoclasts are members of the monocyte/macrophage lineage and are formed by multiple instances of cellular fusion of their mononuclear precursors (5). Monocytes differentiate into osteoclasts in the presence of various molecular signals (6). RANKL, one of the most frequently studied, is a ligand for the receptor activator of nuclear factor-κB (NF-κB; RANK) on osteoclast precursor cells (7). RANKL/RANK signaling activates four pathways that mediate osteoclast formation; NF-κB, c-fos and calcineurin/NFATc1 and three pathways that mediate osteoclast activation; Src and MKK6/p38/MITF and survival; Src and extracellular signal-regulated kinase (8). Osteoblasts produce and secrete osteoprotegerin, a decoy receptor that binds to RANKL and blocks RANKL/RANK interactions and hence suppresses the ability of RANK to increase bone resorption (9). Previous studies have shown that blood monocytes also produce a wide variety of inflammatory factors and transcription factors involved in bone metabolism, including interleukin-1 (10), tumor necrosis factor-α (TNF-α) (11), interleukin-6 (12), platelet-derived growth factor (13), transforming growth factor-β (14), resolvinE1 (15), runt-related transcription factor 2 (Runx2; 16), guanylate binding protein 1 (GBP1), signal transducer and activator of transcription 1 (STAT1), CXC chemokine ligand 10 (CXCL10) (17), chemokine receptor 3, histidine decarboxylase and glucocorticoid receptor genes (18).

However, it is unknown whether other mechanisms regulating these factors are significant in the ability of monocytes to affect bone metabolism. Since biological processes are mediated by multiple, co-regulated genes working in synchrony, certain unknown genes may be assigned potential biological functions when studied in gene sets with known genes and ontology groups (19). Thus, the objective of this study was to screen the differential gene expression in monocytes using a high-throughput microarray platform and to explore gene expression profiles comprehensively by grouping individual differentially expressed genes (DEGs) into gene sets and Gene Ontology (GO) terms. The DEGs between high and low PBM samples were grouped into nine gene sets using the graph-clustering approach. GO term enrichment analysis was applied to identify the relevant molecular functions in response to an impaired PBM. The current study revealed that the DEGs, as precursors of osteoclasts, are functionally involved in the immune response. The stimulus response may contribute to differential osteoclastogenesis, leading to differential PBM levels.

Materials and methods

Affymetrix microarray data

Circulating monocyte affymetrix microarray datasets were accessible from the National Center for Biotechnology Information Gene Expression Omnibus (GEO) data repository (http://www.ncbi.nlm.nih.gov/geo/) using the series accession number GSE7158. Fourteen subjects with extremely high PBM levels and 12 subjects with extremely low PBM levels were selected for DNA microarray experiments. All the recruited volunteers signed an informed consent form prior to entering this project.

Statistical analysis

The limma method (20) was used to identify DEGs. The raw expression datasets from all conditions were normalized using the Robust Multiarray Average (RMA) method with the default settings implemented in Bioconductor and then the linear model was constructed. DEGs with a fold change >1.5 and P<0.05 were selected.

The Pearson correlation coefficient (r) was used to compare the potential correlations between DEGs. Statistical significance was set at r>0.95 and P<0.05. All statistical tests were performed using R language (21).

Network analyses and graph clustering

To identify co-expressed groups, DPClus, a graph clustering algorithm that extracts densely connected nodes as a cluster, was used (22). DPClus is based on the density and periphery tracking of clusters and is freely available from http://kanaya.naist.jp/DPClus/. In the current study, the overlapping mode with the DPClus settings were used. The parameter settings of cluster properties were set; density values were set to 0.5 (23) and minimum cluster size was set to 2.

GO term enrichment analysis

The GO (24) project is a major bioinformatics initiative with the aim of standardizing the representation of genes and gene product attributes across species and databases. The project provides a controlled vocabulary of terms for describing gene product characteristics and gene product annotation data from GO Consortium members, as well as tools to access and process this data.

The DAVID tool (25) was used to identify overrepresented GO terms in biological process. P<0.05 and counts of >2 were set as the threshold for the analysis using the hypergeometric distribution.

Results

Differential gene expression profiling and co-expression network construction

GSE7158 microarray datasets were publicly available from the GEO database. Following microarray analysis, a total of 49 genes were selected as DEGs with a fold change >1.5 and P<0.05. The expression profiling of these 49 DEGs is presented in Fig. 1.

Figure 1 The expression profiling of 49 DEGs. Each row represents the samples and each line represents the expression values of the DEGs. Black indicates low expression, gray indicates medium expression and white indicates high expression. DEGs, differentially expressed genes.

To form the correlations between DEGs, r>0.7 and P<0.05 were selected as the cut-off points. A correlation network was constructed with a total of 159 correlations among 49 DEGs (Fig. 2).

Figure 2 Co-expression network of osteoporosis. A total of 159 correlations with r>0.7 and P<0.05 are exhibited. The nodes indicate the DEGs and the links indicate the high correlation among the DEGs. DEGs, differentially expressed genes.

Graph clustering identifies modules significantly enriched for DEGs contained in GO term pathways

At r>0.7, DPClus (22) identified 9 clusters in the correlation network for osteoporosis, ranging in size from 3–14 genes. Clusters 1, 2, 7, 8 and 9 were connected as they shared the same genes. For example, one gene (epithelial stromal interaction 1, EPSTI1) was shared between clusters 1 and 8; three genes (suppressor of cytokine signaling, SOCS3; superoxide dismutase, SOD2 and activating transcription factor 3, ATF3) were shared between cluster 2 and 7 and one gene (adrenomedullin, ADM) was shared between clusters 2 and 9. The higher the number of genes shared, the more connectivity among them (corresponding to the thicker lines; Fig. 3).

Figure 3 Graph clustering of correlated molecules in osteoporosis (threshold r≥0.7). Using the DPClus algorithm, 9 clusters were extracted for osteoporosis. The internal nodes of the clusters are connected by gray edges; neighboring clusters are connected by black edges.

To assess the significance of the obtained clusters, the overrepresented GO terms were used. Enrichment analysis was performed using the hypergeometrical distribution to find the significant GO term enrichment pathways. In accordance with the graph clustering results, the genes in clusters 1 and 8 were enriched in similar pathways, including immune responses and circulatory system processes. The genes in clusters 2, 7 and 9 were enriched in similar pathways regulating apoptosis and responding to various stimuli, including insulin, hypoxia, nutrients, drugs, radiation and hormones (Table I). Clusters 2 and 7 had the most similar GO term enrichment pathways. These GO biological processes may be relevant to the differentiation of monocytes into osteoclasts.

Table I List of enriched GO terms in clusters 1, 2, 7, 8 and 9 detected by DPClus.

Table I

List of enriched GO terms in clusters 1, 2, 7, 8 and 9 detected by DPClus.

Category	Term	Description	Count	P-value	FDR
Cluster 1	GO:0009615	Response to virus	4	5.86e−5	0.006833
	GO:0006955	Immune response	5	0.00110057	0.062388
	GO:0006952	Defense response	4	0.00882795	0.292359
Cluster 2	GO:0032868	Response to insulin stimulus	3	0.00147149	0.378515
	GO:0007568	Aging	3	0.00177687	0.249653
	GO:0001666	Response to hypoxia	3	0.00262243	0.246268
	GO:0070482	Response to oxygen levels	3	0.00289863	0.208958
	GO:0043434	Response to peptide hormone stimulus	3	0.00344654	0.19991
	GO:0006915	Apoptosis	4	0.00415134	0.200641
	GO:0012501	Programmed cell death	4	0.00432956	0.181443
	GO:0031667	Response to nutrient levels	3	0.00557639	0.202104
	GO:0006916	Anti-apoptosis	3	0.00608264	0.196651
	GO:0008219	Cell death	4	0.00684652	0.199007
	GO:0009991	Response to extracellular stimulus	3	0.00691087	0.184238
	GO:0016265	Death	4	0.00698058	0.171843
	GO:0042981	Regulation of apoptosis	4	0.00934738	0.208115
	GO:0043067	Regulation of programmed cell death	4	0.00960751	0.19967
	GO:0010941	Regulation of cell death	4	0.00970618	0.189438
	GO:0010332	Response to gamma radiation	2	0.01352424	0.240339
	GO:0031100	Organ regeneration	2	0.01527641	0.253599
	GO:0048666	Neuron development	3	0.01586432	0.249457
	GO:0043066	Negative regulation of apoptosis	3	0.01722442	0.255741
	GO:0043069	Negative regulation of programmed cell death	3	0.01768878	0.250411
	GO:0060548	Negative regulation of cell death	3	0.01778231	0.241164
	GO:0009725	Response to hormone stimulus	3	0.01844308	0.23914
	GO:0009628	Response to abiotic stimulus	3	0.01853835	0.231094
	GO:0006873	Cellular ion homeostasis	3	0.01911445	0.228747
	GO:0055082	Cellular chemical homeostasis	3	0.01969829	0.226664
	GO:0009719	Response to endogenous stimulus	3	0.0222131	0.24351
	GO:0050801	Ion homeostasis	3	0.02262762	0.239518
	GO:0030182	Neuron differentiation	3	0.02573112	0.259709
	GO:0019725	Cellular homeostasis	3	0.02888801	0.27855
	GO:0048878	Chemical homeostasis	3	0.03440423	0.314043
	GO:0010212	Response to ionizing radiation	2	0.03494492	0.309692
	GO:0031099	Regeneration	2	0.0400934	0.338357
	GO:0032496	Response to lipopolysaccharide	2	0.04464963	0.360509
	GO:0009266	Response to temperature stimulus	2	0.04805437	0.373652
	GO:0002237	Response to molecule of bacterial origin	2	0.04975276	0.375599
Cluster 7	GO:0070482	Response to oxygen levels	3	3.21e−4	0.079885
	GO:0009314	Response to radiation	3	6.46e−4	0.080292
	GO:0042493	Response to drug	3	7.53e−4	0.062989
	GO:0009991	Response to extracellular stimulus	3	7.81e−4	0.049354
	GO:0055093	Response to hyperoxia	2	0.00177318	0.087833
	GO:0043066	Negative regulation of apoptosis	3	0.00201309	0.08331
	GO:0043069	Negative regulation of programmed cell death	3	0.00206992	0.073801
	GO:0060548	Negative regulation of cell death	3	0.00208138	0.06523
	GO:0009628	Response to abiotic stimulus	3	0.00217417	0.060715
	GO:0010332	Response to gamma radiation	2	0.00509224	0.123857
	GO:0031100	Organ regeneration	2	0.00575517	0.12707
	GO:0048145	Regulation of fibroblast proliferation	2	0.00774219	0.154437
	GO:0042127	Regulation of cell proliferation	3	0.0097487	0.177311
	GO:0042981	Regulation of apoptosis	3	0.01016581	0.172238
	GO:0043067	Regulation of programmed cell death	3	0.01036499	0.164649
	GO:0010941	Regulation of cell death	3	0.01044016	0.156241
	GO:0010212	Response to ionizing radiation	2	0.01324778	0.183871
	GO:0031099	Regeneration	2	0.0152248	0.198085
	GO:0007568	Aging	2	0.02419781	0.283882
	GO:0014070	Response to organic cyclic substance	2	0.02659589	0.294663
	GO:0001666	Response to hypoxia	2	0.02942492	0.308128
	GO:0048545	Response to steroid hormone stimulus	2	0.04197997	0.396429
	GO:0031667	Response to nutrient levels	2	0.0430572	0.390801
	GO:0010035	Response to inorganic substance	2	0.0447791	0.39006
	GO:0006916	Anti-apoptosis	2	0.04499419	0.379328
Cluster 8	GO:0008015	Blood circulation	2	0.01374926	0.839184
	GO:0003013	Circulatory system process	2	0.01374926	0.839184
Cluster 9	GO:0051384	Response to glucocorticoid stimulus	2	0.01149882	0.709892
	GO:0031960	Response to corticosteroid stimulus	2	0.01252751	0.490566
	GO:0048545	Response to steroid hormone stimulus	2	0.02818517	0.6393

[i] GO, gene ontology; FDR, false discovery rate.

Discussion

In the current study, differential expression profiling was systematically investigated and its possible role in the differentiation of osteoclasts was explored. A total of 49 DEGs were identified and correlated to produce 159 network connections. These DEGs were assigned into nine clusters using the graph clustering method in response to different PBM levels. A total of 14 genes were included in cluster 1 [GBP1; interferon-induced protein with tetratricopeptide repeats 2, IFIT2; eukaryotic translation initiation factor 2-α kinase 2, EIF2AK2; interferon-induced protein 44, IFI44; IFI44L; DEAD (Asp-Glu-Ala-Asp) box polypeptide 60, DDX60; HECT and RLD domain containing E3 ubiquitin protein ligase 5, HERC5; radical S-adenosyl methionine domain containing 2, RSAD2; sialic acid binding Ig-like lectin 1, sialoadhesin, SIGLEC1; cytidine monophosphate kinase 2, CMPK2; EPSTI1; interferon, α-inducible protein 6, IFI6; CXCL10; and myxovirus resistance 1, interferon-inducible protein p78, MX1] and 3 genes were involved in cluster 8 (STAT1; EPSTI1; and serpin peptidase inhibitor, clade G, SERPING1). Notably, cluster 8 was connected with all the genes of cluster 1 by STAT1 and EPSTI1 in order to be involved in immune responses and circulatory system processes, as demonstrated in previous studies.

The immune system has been correlated with bone resorption through a complex interaction involving T and B lymphocytes, dendritic cells (DCs), cytokines and cell-cell interactions (26). There is strong evidence that STAT1 is significant in bone metabolism as STAT1 has been reported to be upregulated in the femur tissue of osteoporotic mice (27) and humans (18). STAT1 may serve as a primary mediator of interferon (IFN) signaling pathways involving osteoclast differentiation. Through the p38 MAPK pathway, RANKL stimulates the serine phosphorylation of STAT1, resulting in the migration and adhesion of osteoclast precursors (28). STAT1 interacts with Runx2, an essential transcription factor for osteoblast differentiation, in its latent form in the cytoplasm, thereby inhibiting the nuclear localization of Runx2. This function of STAT1 does not require the Tyr 701 that is phosphorylated when STAT1 becomes a transcriptional activator (29).

The GBP1 gene is also predicted to be involved in bone metabolism or osteoclast differentiation (30) in a STAT1-dependent manner (31). The sumoylation-defective STAT1 mutant exhibits increased induction of GBP1 and transporters associated with antigen presentation 1 (TAP1) transcription (32). The mutation in the STAT1 gene dramatically reduces the inducibility of the GBP1 and TAP1 genes by IFN (33). In this study, STAT1 and GBP1 directly interacted with each other (Figs. 2 and 3).

Chemokines have a potential role in the regulation of osteoclast functions. For example, IFN-γ-inducible protein-10 (CXCL10) is expressed in human osteoclasts with changing expression levels during osteoclast differentiation (34). CXCL10 has been suggested to contribute to osteoclastogenesis by increasing RANKL expression in CD4+ T cells in an animal model of rheumatoid arthritis (35). Notably, previous studies have shown that osteoblasts secrete IFN-β in response to viral infections and that endogenous IFN-β induces CXCL10 and IFI44L production via an IFN-α/β receptor-STAT1 pathway (36,37).

EIF2AK2 is also reported to interact with STAT1 and increase its degradation. Reduction of EIF2AK2 activity also reduces RUNX2 activity and murine osteoblast differentiation (38,39). Therefore, it appears illogical that EIF2AK2 is upregulated in human osteoblasts following IFN-β treatment which results in an inhibition of mineralization (40).

Ten genes were included in cluster 2 [ADM; early growth response 2, EGR2; BCL2-related protein A1, BCL2A1; chromosome 9 open reading frame 72, C9orf72; pleckstrin homology-like domain, family A, member 2, PHLDA2; ATF3; SOCS3; SOD2; 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3, PFKFB3; and pleckstrin homology domain containing, family G (with RhoGef domain) member 2, PLEKHG2], five genes were included in cluster 7 (SOCS3; SOD2; ATF3; cyclin-dependent kinase inhibitor 1A, CDKN1A; and ankyrin repeat domain 28, ANKRD28) and three genes were included in cluster 9 (ADM; interleukin 1 receptor antagonist, IL1RN; and ring finger protein 19B, RNF19B). Cluster 7 was connected with all the genes of cluster 2 and 9 by SOCS3, SOD2 and ATF3. Cluster 9 was connected with all the genes of clusters 2 and 7 by ADM and IL1RN. These findings indicate that SOCS3, SOD2, ATF3 and ADM are significant genes for responding to various stimuli, including insulin, hypoxia, nutrients, drugs, radiation and hormones regulating apoptosis.

The SOCS3 family are cytoplasmic adaptor proteins that negatively regulate various cytokine responses in leukocytes. SOCS3 overexpression augments TGF-β, TNF-α and RANKL-induced osteoclast formation, priming precursors to the osteoclast lineage by suppressing specific anti-osteoclastic JAK/STAT signals (41). Zhang et al demonstrated that a higher SOCS3 expression level is associated with RANKL-mediated alveolar bone loss and enhances CD11c+ DC-derived osteoclastogenesis in vivo and in vitro. The reduced expression of functional SOCS3 in CD11c+ DCs results in significantly lower osteoclastogenesis and dendritic cell-derived osteoclasts development during immune interactions with T cells, based on TRAP expression and bone resorptive activity (42). In SOCS3-deficient bone marrow-derived monocytes, the expression levels of TNF-receptor-associated factor-6 and IκB are drastically reduced. The receptor activation of NF-κB ligand-induced IκB phosphorylation is severely impaired, indicating that SOCS3 regulates osteoclastogenesis by blocking the inhibitory effect of inflammatory cytokines on receptor activation of the NF-κB ligand-mediated osteoclast differentiation signals (43).

ADM is a 52-amino acid peptide first described in a human phaeochromocytoma but has since been identified in numerous tissues, including the bone (44). Systemic administration of ADM stimulates the proliferation of osteoblasts and promotes bone growth (45). Treatment with ADM significantly blunts the apoptosis of serum-deprived osteoblastic cells, evaluated by caspase-3 activity, DNA fragmentation quantification and Annexin V-FITC labeling. This effect is eliminated by calcitonin-related polypeptide α (CGRP1) and insulin-like growth factor-I (46). The selective inhibitor of MAPK kinase (MEK), PD98059, also eliminates the protective effect of ADM on apoptosis and prevents ADM activation of ERK1/2. These data show that ADM acts as a survival factor in osteoblastic cells via a CGRP1 receptor-MEK-ERK pathway, which provides further understanding on the physiological function of ADM in osteoblasts (47).

The SOD2 gene encodes a free radical-scavenging enzyme that removes superoxidate and catalyzes the production of hydrogen peroxide. Oxidative stress is significant in the pathogenesis of osteoporosis (48). Previous studies have revealed that SOD2 is significantly upregulated in circulating monocytes at the mRNA and protein level in vivo in Chinese patients with low versus high hip BMD levels (49). Women with postmenopausal osteoporosis have significantly higher plasma SOD enzyme activity levels than those in controls (50). This indicates that SOD2 is significant in the pathogenesis of osteoporosis, promoting osteoclast differentiation, formation and activity (51).

EGR2 is a highly conserved transcription factor involved in bone remodeling. The upregulation of EGR2 is involved in the biological affinity of titanium for osteogenic cells and in the promotion of osteoblast differentiation (52). Macrophage colony-stimulating factor activates MEK/ERK and induces the MEK-dependent expression of the immediate early gene EGR2. Inhibition of either MEK1/2 or EGR2 increases osteoclast apoptosis (53). Previous studies have revealed a novel role for EGR2 in postnatal skeletal metabolism. EGR2+/− mice reveal a low bone mass phenotype. EGR2 silencing in pre-osteoclasts increases the expression of cFms and the response to macrophage colony-stimulating factor, leading to a cell-autonomous stimulation of cell-cycle progression. Thus, the anti-mitogenic role of EGR2 in pre-osteoclasts is the predominant mechanism underlying the low bone mass phenotype of EGR2-deficient mice (54).

The osteoporotic state increases ATF3 expression in dorsal root ganglia neurons innervating L3 vertebrae (55). BCL2A1, an anti-apoptotic activated macrophage protein, is also heavily overexpressed in osteolysis patients, providing a possible mechanism for the persistence of the particle-laden cells expressing macrophage phenotype activation markers (56).

In conclusion, the present findings shed new light on the biology of osteoporosis and have implications for future research. The changes in the immune system (GBP1, STAT1, CXCL10 and EIF2AK2) and stimulus response (SOCS3, SOD2, ATF3, ADM EGR2 and BCL2A1) may be associated with osteoclast differentiation. This study provides a number of candidate genes that warrant further investigation, including DDX60, HERC5, RSAD2, SIGLEC1, CMPK2, MX1, SERPING1, EPSTI1, C9orf72, PHLDA2, PFKFB3, PLEKHG2, ANKRD28, IL1RN and RNF19B.

References