Purpose: We aimed to explore the molecule mechanism underlying postmenopausal osteoporosis (PMOP) and smoking-related PMOP.
Methods: The microarray data of GSE13850 was downloaded from Gene Expression Omnibus (GEO) database. The differential expression genes (DEGs) between high and low bone mineral density (BMD) samples (DEG1) and postmenopausal female smokers with high and low BMD (DEG2) were analyzed. The GO (Gene ontology) and immune enrichment analysis were performed for investigating the functions of DEGs. Then the protein-protein interaction (PPI) network was constructed by Cytoscape software.
Results: A total of 51 genes in DEG1 and 86 genes in DEG2 were selected. GO term significantly enriched by DEG1 was immune response. GO terms enriched by DEG2 were cell differentiation and MAPKK (Mitogen-activated protein kinase kinase) activity. The PPI network was established with 737 nodes and 774 edges. FLNA (filamin A, alpha), TGFB1 (transforming growth factor, beta 1) and ATM (ataxia telangiectasia mutated) were hub nodes in the network.
Conclusions: The immune system plays a key role in PMOP. The genes of FLNA and TGFB1 may be candidate genes in the pathogenesis of PMOP. MAPKK and ATM may be key factors in smoking-related PMOP.
Key words: postmenopausal osteoporosis, smoking-related postmenopausal osteoporosis, differential expression genes, immune system
Postmenopausal osteoporosis (PMOP) is a prevalent skeletal disease caused by estrogen deficiency, characterized by reduction in bone mineral density (BMD) and microarchitectural deterioration of bone tissue, leading to increased bone fragility and susceptibility to fractures [1-3]. It most common occurs in women after menopause and fractures are the most dangerous aspect of it. Chronic pain and debilitating acute in the elderly is often attributed to fractures from PMOP and they can lead to further disability and even mortality . Patients sustaining a hip fracture have been up to a 33% mortality rate . For most elderly populations, many factors including Vitamin D deficiency , malnutrition  and highdietary protein  will lead to PMOP. Besides, it has been reported that tobacco smoking inhibited the activity of osteoblasts (OBs), and was an independent risk factor for osteoporosis [9,10]. But the mechanisms underlying smoking-associated PMOP remain poorly understood.
BMD is one of important standards to diagnose PMOP. Osteoporosis is defined by theWorld Health Organization(WHO) as a BMD of 2.5standard deviationsor more below the mean peak bone mass (average of young, healthy adults) . Previous studies indicate that the underlying mechanism in all cases of PMOP is an imbalance between bone formation andbone resorption [11-13]. Net loss of bone mass will result from either reduced OB bone formation or excessive osteoclast (OC) bone resorption . The activation of OCs is regulated by various molecular signals, such as receptor activator fornuclear factor κBligand (RANKL) , and Osteoprotegerin . Besides, many studies reported that cytokines (e.g.,interleukins and TNF) and immune cells (e.g., T cells and B cells) may play a key role in regulating bone absorption in PMOP [11,17-19]. However, PMOP is a complex process and the pathogenic mechanisms in immune were still not clearly demonstrated.
In this study, we analyzed the differentially expressed genes (DEGs) using B cells between high and low BMD samples to gain better insights into PMOP. Moreover, B cells of postmenopausal female smokers with high and low BMD were also used to search for DEGs to explain tobacco smoking effects on PMOP. Gene Ontology (GO) and immune enrichment analysis of DEGs were performed. Additionally, the protein-protein interaction (PPI) network was constructed. We expect the study can give a systematic perspective to understanding the mechanisms of PMOP.
Data and methods
Affymetrix microarray data
The gene expression profile data of GSE13850 was obtained from GEO (Gene Expression Omnibus) database (http://www.ncbi.nlm.nih.gov/geo/) which is based on the Affymetrix Human Genome U133A Array . B cells were isolated from whole blood from 20 postmenopausal female smokers, 10 with high BMD and 10 with low BMD. Likewise, B cells were isolated from whole blood from each of 20 women, 10 with high BMD and 10 with low BMD.
Data preprocessing and differential expression analysis
The probe-level data in CEL files were converted into expression measures and performed background correction and quartile data normalization by the robust multiarray average (RMA)  algorithm and subjected to logarithmic transformation with defaulted parameters in R affy package . The limma package  in R was used to identify differentially expressed genes (DEGs). DEG1 was the DEGs from high and low BMD samples; DEG2 was the DEGs from postmenopausal smokers with high and low BMD. The Benjamin and Hochberg (BH)  method was used to adjust the raw p-values. The DEGs only with the |logfoldchange (logFC) | larger than 0.5 and an adjusted p-value less than 0.1 were selected.
DAVID (Database for Annotation, Visualization and Integrated Discovery) is an analytic tool that provides an integrated functional annotation for large list of protein or genes. We used the DAVID to identify GO categories in biological process (BP) based on the hypergeometric distribution with the p-value < 0.005.
Immune enrichment analysis
In order to ascertain if DEGs is enriched in the immune system, we used the information from Immunome database  with hypergeometric distribution method. P < 0.05 was used as the cutoff criterion for the immune enrichment analysis.
The PPI network was constructed among the DEGs and hub genes using information from human protein reference database (HPRD, http://www.hprd.org/) [29,30]. DEG1 and DEG2 were combined for constructing PPI network using the Cytoscape .
Based on the dataset GSE13850, a total of 51 genes in DEG1 and 86 genes in DEG2 with |logFC| > 0.5 and an adjusted p-value < 0.1 were identified.
Gene Ontology enrichment analysis
In order to evaluate the biological functions in POMP, GO enrichment analysis in BP was performed for all DEGs. Total 61 terms in DEG1 and 5 terms in DEG2 were significantly over-represented among these genes. The top five functions of DEG1 and DEG2 were shown in Table 1. In DEG1, most of GO terms were related to immune response, such as “immune effector process”, “cell activation involved in immune response” and “lymphocyte activation involved in immune response”. In DEG2, only one GO term was related to immune response and others were related with cell differentiation (e.g., “regulation of erythrocyte differentiation” and “regulation of megakaryocyte differentiation”), “activation of MAPKK (Mitogen-activated protein kinase kinase) activity” and “histone mRNA catabolic process”.
Immune enrichment analysis
Total three immune functions in DEG1 and two in DEG2 were significant enriched (Table 2). Both in DEG1 and DEG2, the functions were related to “CD molecules” and “Chemokines and receptor”. In these functions, IL4R(interleukin 4 receptor) was the highest frequency genes.
Construction of Protein-Protein Interaction (PPI) Network
The PPI network was established with 737 nodes and 774 edges. In this network, the proteins GNAI2 (Guanine nucleotide-binding protein G), FLNA (Filamin A, alpha), TGFB1 (transforming growth factor, beta 1), JUP (junction plakoglobin) and ATM (Ataxia Telangiectasia Mutated) with high degrees formed a local network (Figure 1). FLNA, TFGB1 and JUP were the co-DEG between DEG1 and DEG2; GNAI2 belonged to DEG1 and ATM belonged to DEG2.
PMOP is a bone disease that is characterized by a decrease in bone density and mass which can lead to fracture . The risk of fracture can reduce living quality for elderly people and increase the risk of mortality . Thus, there is an urgent need to understand the pathogenesis underlying this disease. Besides, tobacco smoking is related with osteoporosis and the mechanisms underlying smoking-associated PMOP is also needed to further study. In this study, we used the gene expression profile downloaded from GEO to identify DEGs using bioinformatics analysis. We identified 51 DEGs between high and low BMD samples, and 86 DEGs between postmenopausal smoker samples with high and low BMD. GO and immune enrichment analysis results showed that immune system and cytokines (e.g., IL4R) played an important role in PMOP.
Bone resorption and remodeling occur throughout life in all bony animals and several key factors that influence the balance between bone forming and bone resorbing . RANKL played an important role in this balance. RANKL is a member of thetumor necrosis factor(TNF)cytokinefamily. In the immune system, RANKL derived from multiple cells, including osteoblasts, stromal cells, and B-and T-cells . T cell activation was reported to induce expression of RANKL and lead to bone loss and an increase of osteoclastogenesis [36,37].Previous studies showed that the concentration of RANKL was increased in postmenopausal women compared with premenopausal women . RANKL is a key factor for OC differentiation, maturation, and activation [39,40]. The increasing OC formation will lead to prolonged bone loss and deeper resorption cavities increasing the fragility of bone . Besides, other cytokines including TNF-α and IL-4 were shown to exert most of their osteoclastogenic activity by inducing RANKLexpression on OBs . TNF-α inhibited the differentiation of OBs . In IL-4-overexpressing mice, the results exhibited a decrease in bone formation and differentiated OBs on the bone surface .
From the result of PPI network, we could find that many proteins, such as FLNA, TGFB1 and ATM were hub nodes in the network.
FLNA is an actin-binding protein that crosslinks actin filaments and regulates signal transduction during cell migration . FLNA is required for osteoclastogenesis by regulating monocyte migration via Rho GTPases . Thus, decreased FLNA can inhibit OCs formation.TGFB1 is a multifunctional set peptides that controlsproliferation,differentiation, and other functions in many cell types . It is characterized as both an immunosuppressive and osteotropic cytokine. TGFB1 down-modulates RANKL expression in OBs, but it may increaseOC formation via action on OC precursors . Therefore, FLNA and TGFB1 were potential target genes in the pathogenesis of PMOP.
Previous studies have been reported that tobacco smoking inhibits the activity of OBs, and is an independent risk factor for osteoporosis [9,10]. Therefore, GO enrichment analysis was performed for the smoking-related DEGs. We found that most of DEGs in GO terms were associated with MAPKK, which suggesting that PMOP was influenced by smoking through activation of MAPKK activity. MAPKK is akinaseenzymewhichphosphorylatesmitogen-activated protein kinase. There are seven genes (MAP2K1, MAP2K3 (aka MKK3), MAP2K6 (aka MKK6) and so on) belong to MAPKKs. Previous study was consistent with our analysis. Koch et al. reported that cigarette smoking enhances MAPK activation to increase lipopolysaccharide-induced IL-8 production in alveolar macrophages.  David et al. showed that MKK3 likely mediated direct effects on OCs through RANKL signaling, MKK6 likely contributed to the production of pro-inflammatory cytokines that promoted OCs formation . Besides, the result of PPI network showed that ATM belonged to DEG2. ATM is a member of the phosphoinositide 3 kinase protein family . The expression of ATM was associated with tobacco smoke exposure in esophageal cancertissues and benzo[a]pyrene diol epoxide in cell lines . In ATM knockout (AtmKO) mice, a stem cell defect due to lack of IGF (Insulin-Like Growth Factor) signaling caused the osteopenic phenotype . Taking all these factors into consideration, we suggested that MAPKK and ATM may be key factors in smoking-related PMOP.
In conclusion, our study showed that immune system played an important role in PMOP. GO enrichment analysis and PPI networks revealed several hub genes that may be involved in the pathogenesis of PMOP and smoking-related PMOP. These findings can enhance our understanding on the molecular mechanism of PMOP. However, further experiments with larger sample size are still needed to confirm our result.