In silico analysis of the molecular mechanism of postmenopausal osteoporosis.

Postmenopausal osteoporosis (PO) is a common disease in females >50 years of age worldwide and is becoming an increasing burden to society. The present study aimed to assess the molecular mechanism of PO using bioinformatic methods. The gene expression data from patients with PO and normal controls were downloaded from the ArrayExpress database provided by European Bioinformatics Institute. Following the screening of the differentially expressed genes (DEGs) using the Limma package in R language, Kyoto Encyclopedia of Genes and Genomes pathways enrichment analysis was performed using the Database for Annotation, Visualization and Integrated Discovery online tools. Sequentially, modulators of the DEGs, including transcription factors (TFs) and microRNAs, were predicted by the ChIP Enrichment Analysis databases and WEB‑based GEne SeT AnaLysis Toolkit system, respectively. In addition, the protein‑protein interaction network of DEGs was constructed via the search tool for the retrieval of interacting genes and then the functional modules were further analyzed via the clusterMaker package and The Biological Networks Gene Ontology package within the Cytoscape software. A total of 482 DEGs, including 279 upregulated and 203 downregulated DEGs, were screened out. DEGs were predominantly enriched in the pathways of fatty acid metabolism, cardiac muscle contraction and DNA replication. TFs, including SMAD4, in addition to microRNAs, including the microRNA‑125 (miR‑125) family, miR‑331 and miR‑24, may be the modulators of the DEGs in PO. In addition, the five largest modules were identified with TTN, L1G1, ACADM, UQCRC2 and TRIM63 as the hub proteins, and they were associated with the biological processes of muscle contraction, DNA replication initiation, lipid modification, generation of precursor metabolites and energy, and regulation of acetyl‑CoA biosynthetic process, respectively. SMAD4, CACNG1 and TRIM63 are suggested to be important factors in the molecular mechanisms of PO, and miR‑331 may be novel potential biomarker for PO.

Introduction

Osteoporosis is a common disease, which is characterized by low bone mass and micro-architectural deterioration of bone tissue, leading to increased bone fragility and an increased risk of fracture (1). Postmenopausal osteoporosis (PO) is one of four types of osteoporosis and is suggested to directly result from a lack of endogenous oestrogen in menopausal females (2). This disease affects millions of females >50 years of age worldwide and treatment of PO is placing an increasing economic burden on society.

The etiology of PO is multifactorial. In addition to the effects of estrogen, calcium and other environmental factors on bone structure and fracture, there is a marked genetic effect on osteoporosis risk in postmenopausal women (3). Mullin et al (4) concluded that genetic variation in ARHGEF3 served a role in the determination of bone density in Caucasian females, and proposed that the RhoGTPase-RhoGEF pathway is associated with PO. A previous study identified that chondroadherin is a novel regulator of bone metabolism, which suppresses pre-osteoclast motility and bone resorption, with a potential effect for the treatment of PO (5).

A previous study described the transcriptional alterations in 84 trans-iliac bone biopsies associated with bone mineral density (BMD) variations in postmenopausal women (6). This previous study identified that sclerostin and dickkopf homolog 1, both involved in the Wnt signaling pathway, exhibited a clear correlation and was involved in bone metabolism. Jemtland et al (7) analyzed these data to identify osteoporosis candidate genes and identified that the transcription factor (TF), SOX4, and the bone matrix proteins, MMP13 and MEPE, were all downregulated in osteoporosis. However, the authors focused on the gene expression levels associated with the BMD variation, rather than an intensive analysis of the gene regulation changes and interactions in PO, of which the underlying mechanisms remain to be elucidated.

Therefore, the data mentioned above was obtained and the differentially expressed genes (DEGs) between PO samples and healthy controls were assessed by genome-wide microarray analysis, in addition to performing a more comprehensive analysis, to achieve an improved understanding of the mechanisms of this disease. Various bioinformatic methods were applied to identify the potential modulators, including TFs and microRNAs of the DEGs, and the significant pathways associated with the DEGs, in addition to identification of the functional modules in the interaction network of the DEGs involved in PO.

Materials and methods

Data acquisition and preprocessing

The expression data, numbered E-MEXP-1618 (6), were downloaded from the ArrayExpress database (8) provided by the European Bioinformatics Institute (Saffron Walden, UK). All 66 trans-iliac bone biopsy samples were obtained from postmenopausal females, including 27 osteoporosis patients (mean age, 69.6 years, range, 51.6–86.1 years) and 39 healthy controls (mean age, 61.7 years, range, 49.7–80.9 years).

The primary data was standardized and transformed into expression values using the Robust Multi-array Average algorithm (www.bioconductor.org) (9) in R language (version 2.4.1), based on the microarray platform Affymetrix GeneChip Human Genome U133 plus 2.0 (Affymetrix, Inc., Santa Clara, CA, USA).

DEG screening

The DEGs were screened out by significance analysis using the Empirical Bayes methods within Limma package (10) in R language. The adjusted P-value represents the P-value adjusted using the Benjamini-Hochberg method (11), following Student's t-test, with <0.1 as a cut-off criterion.

Pathway enrichment of DEGs

Pathway enrichment analysis of the DEGs was performed using the Database for Annotation, Visualization and Integrated Discovery (12) online tools version 6.7 based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway databases (13). A false discovery rate <0.05 was set as a cut-off criterion.

Prediction of DEG regulation

To determine the potential transcriptional and post-transcriptional modulators, the ChIP Enrichment Analysis (ChEA) database (http://amp.pharm.mssm.edu/lib/chea.jsp) (14) and the WEB-based GEne SeT AnaLysis Toolkit (WebGestalt) system (http://bioinfo.vanderbilt.edu/webgestalt) (15) were used to predict the TFs and microRNAs of DEGs, respectively. The ChEA database contains interactions describing the regulation of TFs on target genes and P<[0.05/Σ(TFs)] was set as the cut-off criterion. The WebGestalt system includes interactions describing the binding of microRNAs to the 3′ untranslated region of the target genes and an adjusted P-value ≤0.05 was set as the cut-off criterion.

Protein-protein interaction (PPI) network construction and analysis

The PPI pairs of the screened DEGs were analyzed using the Search Tool for the Retrieval of Interacting Genes (STRING) software 9.0 (16). The pairs with combined scores >0.4 were used for the PPI network construction using Cytoscape software 2.8 (17). Furthermore, the modules with close internal communication were screened out with the Markov Cluster (MCL) (18) algorithm in the clusterMaker package (19) within the Cytoscape software. In addition, the biological processes in which the screened modules were enriched were identified by the Biological Networks Gene Ontology (BiNGO) package 2.44 (20) within the Cytoscape software package.

Results

Multiple DEGs are involved in various pathways

A total of 482 DEGs, including 279 upregulated and 203 downregulated DEGs, were screened out in the samples from patients with osteoporosis when compared with the healthy control samples. The DEGs were subjected to KEGG pathway enrichment analysis. As presented in Table I, the upregulated genes were predominantly enriched in the pathway of fatty acid metabolism and the downregulated genes were predominantly enriched in the pathway of DNA replication. Notably, cardiac muscle contraction was also a significant pathway in which the upregulated genes were enriched.

Table I

Key enriched pathways of the differentially expressed genes.

Regulation	Term	Count	P-value	Gene	FDR
Upregulated	hsa00071: fatty acid metabolism	8	4.29−6	ACADSB, ACSL1, ACADM, ADH5, HADH, ACSL3, ACAT1, HADHB	0.000480
	hsa04260: cardiac muscle contraction	8	3.57−4	UQCRC2, ATP1B1, ACTC1, COX7A1, MYH7, MYH6, ATP1A2, CACNG1	0.019784
	hsa00280: valine, leucine and isoleucine degradation	6	8.84−4	ALDH6A1, ACADSB, ACADM, HADH, ACAT1, HADHB	0.032468
	hsa04920: adipocytokine signaling pathway	7	9.69−4	LEP, PRKCQ, PPARA, ACSL1, PRKAA2, ACSL3, CHUK	0.02678
Downregulated	hsa03030: DNA replication	5	5.98−4	MCM7, LIG1, MCM3, RNASEH2C, MCM5	0.047308

FDR, false discovery rate.

Up and downregulated DEGs are modulated by TFs and microRNAs

The transcriptional and post-transcriptional modulators of the DEGs were also predicted. The ChEA analysis included 94 TFs, therefore P<0.0005 was set as the criterion. This analysis demonstrated that the TFs, including HNF4A, SMAD4 and SOX2, were significantly associated with the upregulated DEGs. HNFA4 had the most significant P-value, therefore, may regulate 110 genes. SOX2, which exhibited the second most significant P-value was suggested to regulate 72 genes, followed by SMAD4, which may regulate 59 genes (Fig. 1A). Additionally, FOXP1 and SPI1 were significantly associated with the downregulated DEGs and may regulate 65 and 55 genes, respectively (Fig. 1B).

Figure 1

TF regulatory networks of the DEGs. The TF regulatory network of the (A) upregulated and (B) downregulated DEGs. The results were analyzed using the ChEA database and 94 TFs were analyzed, with the P<[0.05/Σ(TFs)] and videlicet 0.0005 as the criterion. Orange represents upregulation, blue represents downregulation and yellow represents no clear alteration in the expression. TFs are shown as diamonds and DEGs are shown as circles. TF, transcription factor; DEGs, differentially expressed genes.

The WebGestalt analysis indicated that the microRNAs, including the microRNA-125 (miR-125) family, miR-331 and miR-24, potentially modulated the downregulated DEGs. This suggested that these microRNAs may be in an active state in PO. Additionally, the miR-26, miR-15 and miR-200 families were identified to possibly modulate the upregulated DEGs, which suggested their inactive state in PO (Fig. 2).

Figure 2

microRNA regulatory network of the DEGs. The results were analyzed by the WebGestalt system and an adjusted P<0.05 was set as the cut-off criterion. Orange represents upregulation, blue represents downregulation and grey represents no clear alteration in the expression. microRNAs are shown as rounded triangles and DEGs are shown as circles. DEGs, differentially expressed genes.

PPI network construction

Interactions between protein pairs were identified using STRING software. The PPI network of the DEGs in PO was constructed using Cytoscape software while the nodes with no connections were filtered out. A total of 146 DEGs with 271 interactions were detected. TTN, MYOZ2 and LDB3 were identified to possess the highest degrees of connectivity and were observed to be involved in 22, 19 and 17 pairs of interactions, respectively.

Twelve functional modules associated with biological processes

A total of 12 modules, which contained at least three DEGs in the PPI network, were constructed, of which the five largest modules were identified with TTN (Fig. 3A), L1G1 (Fig. 3B), ACADM (Fig. 3C), UQCRC2 (Fig. 3D) and TRIM63 (Fig. 3E) as the hub proteins, respectively. BiNGO analysis demonstrated that these five modules were associated with the biological processes of muscle contraction, DNA-dependent DNA replication initiation, lipid modification, generation of precursor metabolites and energy, and regulation of the acetyl-CoA biosynthetic process from pyruvate, respectively (Table II).

Figure 3

Major functional modules of the protein-protein interaction network of the DEGs. (A) The TTN centered module was associated with muscle contraction; (B) the L1G1 centered module was associated with DNA-dependent DNA replication initiation; (C) the ACADM centered module was associated with lipid modification; (D) the UQCRC2 centered module was associated with the generation of precursor metabolites and energy; (E) the TRIM63 centered module was associated with the regulation of acetyl-CoA biosynthetic process from pyruvate. The modules and their associated biological processes were identified by the MCL algorithm and Biological Networks Gene Ontology package within Cytoscape, respectively. Orange represents upregulation and blue represents downregulation. The sizes of the nodes represent the quantity of the interactions. DEGs, differentially expressed genes.

Table II

Biological processes involved in the five largest functional modules.

Module hub	GO-term (biological process)	P-value	Adjusted P-value	Gene
TNN	Muscle contraction	5.04−17	2.11−14	TRDN, MYOM2, ANKRD2, SMPX, MYH7, ACTN2, MYH6, MYOM1, CACNG1, TTN, MYOT
L1G1	DNA-dependent DNA replication initiation	2.51−7	4.14−5	MCM7, MCM3, MCM5
ACADM	Lipid modification	5.78−9	1.04−6	ACADM, HADH, ACAT1, HADHB
UQCRC2	Generation of precursor metabolites and energy	2.82−8	1.33−6	UQCRC2, NDUFA5, COX7A1, SUCLA2, MDH1
TRIM63	Regulation of acetyl-CoA biosynthetic process from pyruvate	9.66−6	5.77−4	PDK4, DLAT

GO, gene ontology.

Discussion

In order to elucidate the molecular mechanisms of PO, the gene expression profiles were systematically analyzed using bioinformatic approaches. A total of 482 DEGs, including 279 upregulated and 203 downregulated DEGs, were screened out in patients with PO. The biological functions of these DEGs were further assessed based on pathway enrichment data. Further modulator prediction identified the potential TFs and microRNAs, which may regulate the DEGs in PO. In addition, the functional modules in the PPI network of the DEGs were identified, of which certain modules were clearly involved in PO.

TFs prediction in the present study identified that 59 of the upregulated DEGs were targets of SMAD4, which is the only member of common-mediator SMAD (co-SMAD) class. A previous study demonstrated that defects in bone morphogenetic protein (BMP)-SMAD signaling led to bone-associated disorders, including osteoporosis (21). Association of BMPs to BMP receptors on the cell surface leads to the activation of the formation of Smad4 and Smad1/5/8 complexes (22). The complexes subsequently translocate into the nucleus and bind to the consensus DNA sequence to modulate the transcription of BMP target genes (23). In addition, SMAD4 was suggested to function as a transcriptional co-repressor for estrogen receptor α (ERα) by forming a complex when ERα binds to the estrogen-responsive element within the promoters of estrogen target genes (24). Estrogen-associated therapies are widely used for the treatment of PO (25–27) and antiestrogens are able to enhance the endogenous interactions between SMAD4 and ER (24). According to the data from the present study, the expression of SMAD4 itself was upregulated in PO, suggesting that it may be activated by a mechanism of cross-talk between BMP-SMAD signaling and ERα-estrogen interaction to subsequently regulate downstream target genes involved in this disease.

Several microRNAs have been identified to serve important roles in PO, including miR-148a, which promotes osteoclastogenesis (28), miR-133a (29) and miR-422a (30), which are upregulated with low BMD in human circulating monocytes (osteoclast precursors). However, the roles of microRNAs in PO remain to be elucidated. The microRNA prediction in the present study suggested that the miR-24 and miR-125 families may be activated in PO to suppress the expression levels of a series of genes, which is consistent with a previous observation that miR-24 and miR-125b were significantly upregulated in the serum and bone tissue of patients with PO (31). Another microRNA suggested to be activated was miR-331, which may be associated with miR-24 by its interactions with Wolf-Hirschhorn Syndrome Candidate 1 (32). This suggested that miR-331 may be a novel potential biomarker for PO.

In addition, the three hub proteins identified in the PPI network of DEGs were TTN, MYO2 and LDB3, which were demonstrated to be associated in one functional module of muscle contraction. Notably, KEGG pathway enrichment analysis demonstrated that the cardiac muscle contraction pathway was significantly associated with the upregulated DEGs. Voltage-dependent calcium channel γ1 (CACNG1) is a DEG, which was observed to be associated with the muscle contraction module and the cardiac muscle contraction pathway. It has been previously reported that CACNG1 in different cell types may be important in the mechanism of the regulation of Ca2+ channel function (33). A previous study indicated that patients with PO often exhibit hypercalciuria with normal blood Ca2+ levels (34). Therefore, CACNG1 may be important in the underlying mechanisms of PO through Ca2+ regulation in the muscle contraction process.

Another functional module of the PPI network identified in the present study with TRIM63 as the hub, was demonstrated to be associated with the regulation of the acetyl-CoA biosynthetic process from pyruvate, which is an important pathway in the human metabolic process. TRIM63, also termed muscle-specific ring finger protein 1, has been reported as an E3 ubiquitin ligase expressed predominantly in muscular tissue. Azuma et al (35) proposed that the overexpression of TRIM63 increased the expression of an osteoblastic differentiation marker gene, alkaline phosphatase, resulting in reduced proliferation. In addition, TRIM63 was identified to be involved in the two major bone remodeling activities, osteoblastic bone formation and osteoclastic bone resorption (36). According to the present study, the four genes, LMCD1, PDK4, USP13 and GFM1, encoding the proteins which interacted with TRIM63, were all upregulated in PO, indicating a potential synergistic effect of these proteins with TRIM63 in the bone remodeling activities in PO.

In conclusion, the DEGs in PO were screened comparing them with the normal controls and further intensive bioinformatic analysis, including pathway enrichment, modulator prediction of TFs and microRNAs, PPI network analysis and functional module identification was performed on the DEGs. It was suggested that SMAD4, CACNG1 and TRIM63 may have important roles in the molecular mechanism of PO and that miR-331 may be novel potential biomarker for PO. The present study may provide bioinformatic support for further investigations into the mechanisms of PO. However, associated experimental data are necessary to confirm the conclusions of the present study.