Study on Key Genes and Regulatory Networks Associated with Osteoporosis by Microarray Technology

Abstract

Objective: DNA microarray data of patients with osteoporosis were compared with that of healthy people to identify key genes and thus disclose the underlying regulatory network. Methods: Microarray dataset GSE35958 was downloaded from the Gene Expression Omnibus database, including five gene chips from patients with primary osteoporosis and four from age-matching nonosteoporosis controls. Raw data were preprocessed and differentially expressed genes (DEGs) were identified by the t-test. Then, function and pathway annotations were given by gene ontology (GO) and KEGG. The regulatory network for the DEGs was established from the aspects of transcription factors and microRNAs (miRNAs). The regulators of the miRNAs were also predicted by the MATCH algorithm. Results: A total of 274 DEGs were obtained with 47 significantly over-represented GO terms and 2 KEGG pathways. Transcriptional and post-transcriptional regulatory networks were established for the DEGs. Moreover, upstream regulators of the miRNAs were also obtained. Conclusion: A range of genes, which might be implicated in the development of osteoporosis were obtained in the present study. Our findings are of possible benefit for the understanding of the unsolved regulatory mechanisms, and future clinical diagnosis as well as treatment.

Introduction

Osteoporosis is a systemic skeletal disease characterized by loss of bone mass and it leads to increased fractures. Sex steroid deficiency, aging, and glucocorticoid excess are thought to be the major risk factors (Manolagas, 2000; Sun et al., 2013). The morbidity in women is higher than that in men because women experience a rapid phase of loss after menopause due to the loss of estrogen (Cummings et al., 1985).

Several therapies have been on the market and more are in development for osteoporosis (Diab and Watts, 2012; Silva-Fernandez et al., 2012; Walsh, 2012). However, it still lacks premonitory symptoms, and is usually diagnosed when fractures happen. The measurement of bone mineral density is the most important diagnostic method currently available (Kanis, 2002; Jergas, 2011). In addition to reducing bone mineral density, microarchitectural deterioration and gene expression alterations also occur (Lin and Lane, 2004). Therefore, it is urgent to elucidate the regulatory mechanism from the molecular level, which would benefit the early diagnosis and prevention of the disease.

In the present study, microarray technology was utilized to identify differentially expressed genes (DEGs) in primary osteoporosis compared with normal controls. Function and pathway annotation were given for these DEGs. To further elucidate their roles in the pathogenesis, the relationships between two important regulators, transcription factors (TFs) and microRNAs (miRNAs), and their target genes were retrieved, and the regulatory networks were established, which demonstrate the molecular pathways involved in osteoporosis and thus provide a guide for future drug development.

Materials and Methods

Microarray data

Microarray dataset GSE35958 was downloaded from the GEO database (www.ncbi.nlm.nih.gov/geo/), including five genechips from patients with primary osteoporosis (minor fracture of the femoral neck) and four from age-matching nonosteoporosis controls. The age is between 79 and 94 years. RNA samples were extracted from bone marrow. Raw data were collected in the following platform: GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array.

TF, miRNA, and their target genes

The TRANSFAC database collects information about transcription in eukaryotes, and a total of 895 TFs and 1847 TF-gene regulatory relationships are included. Information about miRNA comes from miRBase database, and target genes are predicted by 7 different algorithms: PicTar DIANA-microT, miRanda the miRBase RNAhybrid RNA22, TargetScan. To reduce the false positive, relationships predicted by at least three algorithms were regarded as confident.

Screening of DEGs

Microarray data were processed by mapping probes to genes. The average was determined as the expression level when more than one probe was mapped to the same gene, and such probes were removed, while more than one gene was mapped. Then, mRNA expression data underwent log2 conversion. DEGs were screened out by the t-test and multiple test correction. The false discovery rate (FDR)<0.05 and fold change>2 were set as the cutoffs.

Function and pathway annotations

To describe the biological roles of the DEGs, function and pathway annotations were performed by gene ontology (GO) and KEGG, respectively. Package Gostats of R was used for these tasks.

Establishment of regulatory network

TFs and miRNAs associated with the DEGs were retrieved to disclose the pathogenesis of osteoporosis from the aspects of transcription and post-transcription. Relationships between TFs and genes, extracted from the TRANSFACC database, were used to establish a regulatory network. Information about miRNAs and target genes was collected and the corresponding network was set up. Moreover, upstream regulators of these miRNAs were identified from the UCSC database using the MATCH algorithm. The following parameters were adopted in the calculation: vertebrates and minimize false positive matches. All the information was integrated to comprehensively reflect the pathogenic process of osteoporosis.

Results

Identification of DEGs

Differential expression analysis was conducted between patients with osteoporosis and controls. A total of 274 DEGs were obtained according to the criteria: fold change>2 and FDR<0.05.

To reveal the biological functions of the DEGs, GO and KEGG annotation are given (Table 1). A total of 47 GO terms and 2 KEGG pathways were significantly over-represented (p<0.05).

Table 1.

Functional Annotations for the Regulatory Network

ID	Term	p-Value
GO:0001568	Blood vessel development	0.006522
GO:0001944	Vasculature development	0.007596
GO:0002237	Response to molecule of bacterial origin	0.030101
GO:0006638	Neutral lipid metabolic process	5.79E-04
GO:0006639	Acylglycerol metabolic process	0.004421
GO:0006641	Triglyceride metabolic process	0.002743
GO:0006916	Antiapoptosis	0.02305
GO:0006937	Regulation of muscle contraction	0.016866
GO:0007090	Regulation of S phase of mitotic cell cycle	0.019651
GO:0007167	Enzyme-linked receptor protein signaling pathway	0.007332
GO:0007420	Brain development	0.044958
GO:0007584	Response to nutrient	0.003087
GO:0009267	Cellular response to starvation	0.009301
GO:0010565	Regulation of cellular ketone metabolic process	0.007576
GO:0010604	Positive regulation of macromolecule metabolic process	0.045125
GO:0010605	Negative regulation of macromolecule metabolic process	0.043418
GO:0019217	Regulation of fatty acid metabolic process	0.004421
GO:0019752	Carboxylic acid metabolic process	0.041381
GO:0030003	Cellular cation homeostasis	0.023417
GO:0030100	Regulation of endocytosis	0.009598
GO:0030879	Mammary gland development	0.010728
GO:0031324	Negative regulation of cellular metabolic process	0.037481
GO:0031325	Positive regulation of cellular metabolic process	0.017014
GO:0031326	Regulation of cellular biosynthetic process	0.049099
GO:0031399	Regulation of protein modification process	0.02007
GO:0031400	Negative regulation of protein modification process	0.005807
GO:0031401	Positive regulation of protein modification process	0.043354
GO:0031669	Cellular response to nutrient levels	3.17E-04
GO:0032268	Regulation of cellular protein metabolic process	0.013198
GO:0032269	Negative regulation of cellular protein metabolic process	0.011766
GO:0032270	Positive regulation of cellular protein metabolic process	0.014727
GO:0032496	Response to lipopolysaccharide	0.02106
GO:0033261	Regulation of S phase	0.030027
GO:0042594	Response to starvation	0.019497
GO:0042981	Regulation of apoptosis	0.026403
GO:0043067	Regulation of programmed cell death	0.028766
GO:0043627	Response to estrogen stimulus	0.014501
GO:0046486	Glycerolipid metabolic process	0.001713
GO:0048514	Blood vessel morphogenesis	0.025885
GO:0050801	Ion homeostasis	0.024785
GO:0050880	Regulation of blood vessel size	0.035708
GO:0051246	Regulation of protein metabolic process	0.001299
GO:0051247	Positive regulation of protein metabolic process	0.018513
GO:0051248	Negative regulation of protein metabolic process	0.004155
GO:0055082	Cellular chemical homeostasis	0.015278
GO:0060627	Regulation of vesicle-mediated transport	0.042467
GO:0060696	Regulation of phospholipid catabolic process	0.027149
Kegg:04510	Focal adhesion	8.25E-04
Kegg:05120	Epithelial cell signaling in Helicobacter pylori infection	0.032778

GO, gene ontology.

Regulatory relationships

The TFs regulating DEGs were retrieved from the TRANSFAC database. Interactions between 49 TFs and 17 DEGs were obtained. As shown in Figure 1A, most of the TFs had only one target DEG.

FIG. 1.

(A) The distribution for the number of target genes of transcription factors (TFs). (B) The distribution for the number of target genes of microRNA (miRNAs).

miRNAs were acquired from a miRBase database. Target genes were predicted by seven algorithms (PicTar, DIANA-microT, miRanda, miRBase, RNAhybrid, RNA22, and TargetScan), and 52828 miRNA-target gene relationships were presented in at least three algorithms and were considered as of high confidence. Then, 1014 relationships involving DEGs were extracted, containing 355 miRNAs and 137 target DEGs. The distribution for the number of target genes corresponding to one miRNA was shown in Figure 1B.

Since the genomic annotation for miRNAs is incomplete, the primary transcripts of most miRNAs are still unknown. In the present study, sequence information of the precursor miRNAs was acquired from the UCSC database and the MATCH algorithm was adopted to identify upstream regulators. Finally, 40 regulatory relationships, including 14 TFs and 8 miRNAs were revealed.

Regulatory network

The transcriptional regulatory networks (Fig. 2) and post-transcriptional regulatory networks were established based upon the relevant information. Three modules were extracted from the network using the MCODE algorithm to better describe their functions (Fig. 3). The most significantly over-represented biological process for each module was biosynthesis, an enzyme-linked receptor protein signaling pathway and apoptosis regulation. Besides, regulators for the miRNAs in the modules were uncovered (Fig. 4). According to the prediction results, three miRNAs in module 1 were regulated by POU2F1 (Fig. 4).

FIG. 2.

The TF-mediated transcriptional regulation networks. The round node represents the differentially expressed gene (DEG), and the triangular node represents the TF.

FIG. 3.

The function modules from miRNA-involved regulatory network. The rhombic node represents miRNA, the round node represents the differentially expressed gene.

FIG. 4.

Predicted regulators for the miRNAs in the three modules. The rhombic node represents miRNA, the round node represents the differentially expressed gene, and the triangular node represents the TF.

Discussion

In the present study, microarray data were analyzed to globally identify DEGs associated with osteoporosis. GO analysis showed that 47 terms were significantly over-represented. Then, a TF-mediated transcriptional regulatory network and a miRNA-mediated post-transcriptional regulatory network were established to demonstrate the molecular mechanisms. Moreover, regulators of the miRNAs were also predicted and the network was generated. Our study comprehensively described the regulatory network from different aspects.

The TF-target gene regulatory network provides much information. Cytochrome P450, family 1, subfamily A, polypeptide 1 (CYP1A1) is a member of the cytochrome P450 superfamily of enzymes, which catalyzes many reactions involved in drug metabolism and synthesis of cholesterol, steroids, and other lipids. Its polymorphisms have effects on estrogen metabolism and thus are linked with bone density (Napoli et al., 2005; Quan et al., 2009). Its TFs include the aryl hydrocarbon receptor (AHR), nuclear transcription factor Y-alpha, and Ku80 (XRCC5). A previous study indicates that AHR signals have a critical effect on the differentiation of osteoblasts and influence bone modeling in vivo (Naruse et al., 2002). Rayalam et al. (2011) consider that the AHR plays an important role in diverse cellular functions like detoxification, adipogenesis, and adipocyte and osteoblast differentiation. It could be speculated that CYP1A1 is a key mediator. Ku80 is also linked to osteoporosis as Ku80 mutant mice exhibit several characteristics of aging, including osteoporosis (Hasty et al., 2003; Watanabe and Hishiya, 2005). Nitric oxide synthase 3 is involved in the synthesis of nitric oxide, which acts as a biologic mediator in several processes, such as bone loss through osteoblasts as well as osteomyelitis (Asensi et al., 2007).

The miRNA-target gene regulatory network is also beneficial in promoting our understanding about the underlying mechanisms. In module 1, both the pleckstrin homology domain interacting protein (PHIP) and the structure-specific recognition protein 1 (SSRP1) were regulated by three miRNAs: miR-195, miR-497, and miR-15a. A number of studies have indicated the linkage between insulin and osteoporosis (Thrailkill et al., 2005; Abdulameer et al., 2012; Xia et al., 2012). Thus, we speculated that PHIP might be implicated in osteoporosis as it could bind the pleckstrin homology domain of the insulin receptor substrate-1 (Farhang-Fallah et al., 2000) and modulate insulin signaling. Maire et al. (2011) report that hsa-miR-195 is upregulated in osteosarcoma compared with normal osteoblasts. In addition, the myeloblastosis oncogene (MYB) was predicted to be the regulator of miR-195. The study by Willing et al. (1998) shows that the MYB has an influence on the estrogen receptor 1 (ESR1), which is closely associated with bone density and fractures, as well as the pathogenesis of osteoporosis (Styrkarsdottir et al., 2009). BMP and activin membrane-bound inhibitor homolog (BAMBI) was included in module 2 and regulated by miR-19a and miR-19b. BAMBI is a transmembrane glycoprotein related to the type I receptors of the transforming growth factor-beta family and thus takes a part in some processes. It may bind to BMP and thus is associated with bone formation (Li, 2008). Kitazawa et al. (2005) further suggest that BAMBI promoter hypermethylation may be involved in metaplastic bone formation. Ubiquitin-mediated protein degradation has been confirmed to be linked with bone homeostasis (Yamashita et al., 2005; Gunther Sillero et al., 2008) and CUL2 is an important component of the E3 ubiquitin-protein ligase complex. Thus, altered expression in primary osteoporosis is really worth further investigation.

Although much information about the involvement of DEGs and their regulators was provided in the study, more microarray data are needed to completely present all the molecular mechanisms. Anyway, our study is beneficial for understanding the pathogenesis of osteoporosis and also provides clues for future biomarker discovery and drug development.

Footnotes

Acknowledgment

This study was supported by the National Natural Science Foundation of China (81000779).

Author Disclosure Statement

No competing financial interests exist.

References

Abdulameer

, Sulaiman

, Hassali

et al. 2012. Osteoporosis and type 2 diabetes mellitus: what do we know, and what we can do? Patient Prefer Adherence, 6:435-448.

Asensi

, Montes

, Valle

et al. 2007. The NOS3 (27-bp repeat, intron 4) polymorphism is associated with susceptibility to osteomyelitis. Nitric Oxide, 16:44-53.

Cummings

, Kelsey

, Nevitt

et al. 1985. Epidemiology of osteoporosis and osteoporotic fractures. Epidemiol Rev, 7:178-208.

Diab

, Watts

. 2012. Bisphosphonates in the treatment of osteoporosis. Endocrinol Metab Clin North Am, 41:487-506.

Farhang-Fallah

, Yin

, Trentin

et al. 2000. Cloning and characterization of PHIP, a novel insulin receptor substrate-1 pleckstrin homology domain interacting protein. J Biol Chem, 275:40492-40497.

Gunther Sillero

, de Diego

, Perez-Zuniga

et al. 2008. Synthesis of bisphosphonate derivatives of ATP by T4 DNA ligase, ubiquitin activating enzyme (E1) and other ligases. Biochem Pharmacol, 75:1959-1965.

Hasty

, Campisi

, Hoeijmakers

et al. 2003. Aging and genome maintenance: lessons from the mouse? Science, 299:1355-1359.

Jergas

. 2011. Radiological diagnosis of osteoporosis—assessment of fracture risk. MMW Fortschr Med, 153:33-36.

Kanis

. 2002. Diagnosis of osteoporosis and assessment of fracture risk. Lancet, 359:1929-1936.

10.

Kitazawa

, Kitazawa

, Obayashi

et al. 2005. Desmoid tumor with ossification in chest wall: possible involvement of BAMBI promoter hypermethylation in metaplastic bone formation. J Bone Miner Res, 20:1472-1477.

11.

. 2008. Bone morphogenetic protein-Smad pathway as drug targets for osteoporosis and cancer therapy. Endocr Metab Immune Disord Drug Targets, 8:208-219.

12.

Lin

, Lane

. 2004. Osteoporosis: a review. Clin Orthop Relat Res, 425:126-134.

13.

Maire

, Martin

, Yoshimoto

et al. 2011. Analysis of miRNA-gene expression-genomic profiles reveals complex mechanisms of microRNA deregulation in osteosarcoma. Cancer Genet, 204:138-146.

14.

Manolagas

. 2000. Birth and death of bone cells: basic regulatory mechanisms and implications for the pathogenesis and treatment of osteoporosis. Endocr Rev, 21:115-137.

15.

Napoli

, Villareal

, Mumm

et al. 2005. Effect of CYP1A1 gene polymorphisms on estrogen metabolism and bone density. J Bone Miner Res, 20:232-239.

16.

Naruse

, Ishihara

, Miyagawa-Tomita

et al. 2002. 3-Methylcholanthrene, which binds to the arylhydrocarbon receptor, inhibits proliferation and differentiation of osteoblasts in vitro and ossification in vivo. Endocrinology, 143:3575-3581.

17.

Quan

, Yahata

, Tamura

et al. 2009. Relationship between single nucleotide polymorphisms in CYP1A1 and CYP1B1 genes and the bone mineral density and serum lipid profiles in postmenopausal Japanese women taking hormone therapy. Menopause, 16:171-176.

18.

Rayalam

, Yang

, Della-Fera

et al. 2011. Novel molecular targets for prevention of obesity and osteoporosis. J Nutr Biochem, 22:1099-1104.

19.

Silva-Fernandez

, Rosario

, Martinez-Lopez

et al. 2012. Denosumab for the treatment of osteoporosis: a systematic literature review. Reumatol Clin, 9:42-52.

20.

Styrkarsdottir

, Halldorsson

, Gretarsdottir

et al. 2009. New sequence variants associated with bone mineral density. Nat Genet, 41:15-17.

21.

Sun

, Guo

, Chen

et al. 2013. Significant pathways detection in Osteoporosis based on the bibliometric network. Eur Rev Med Pharmacol Sci, 17:1-7.

22.

Thrailkill

, Lumpkin

Jr. , Bunn

et al. 2005. Is insulin an anabolic agent in bone? Dissecting the diabetic bone for clues. Am J Physiol Endocrinol Metab, 289:E735-E745.

23.

Walsh

. 2012. Treatment options for postmenopausal osteoporosis. Ther Clin Risk Manag, 8:367-368.

24.

Watanabe

, Hishiya

. 2005. Mouse models of senile osteoporosis. Mol Aspects Med, 26:221-231.

25.

Willing

, Sowers

, Aron

et al. 1998. Bone mineral density and its change in white women: estrogen and vitamin D receptor genotypes and their interaction. J Bone Miner Res, 13:695-705.

26.

Xia

, Zhong

, Huang

et al. 2012. The relationship between insulin resistance and osteoporosis in elderly male type 2 diabetes mellitus and diabetic nephropathy. Ann Endocrinol (Paris), 73:546-551.

27.

Yamashita

, Ying

, Zhang

et al. 2005. Ubiquitin ligase Smurf1 controls osteoblast activity and bone homeostasis by targeting MEKK2 for degradation. Cell, 121:101-113.