Bioinformatics Analysis Reveals the Altered Gene Expression of Patients with Postmenopausal Osteoporosis Using Liuweidihuang Pills Treatment

Postmenopausal osteoporosis (PMOP), as well as its associated increased risk for fragility fracture, is one of the most disabling consequences of aging in women. This present study aimed to identify candidate genes that involve pathogenesis of PMOP and the therapeutic mechanism of Liuweidihuang (LWDH) pills on PMOP. We integrated microarray datasets of PMOP derived from the Gene Expression Omnibus (GEO) to screen differentially expressed genes (DEGs) between PMOP and normal controls as well as patients with PMOP and patients after treatment of LWDH pills. GO and KEGG enrichment analysis for DEGs were performed. The shared DEGs, associated with both the pathogenesis of PMOP and the therapeutic mechanism of LWDH, were further analyzed by protein-protein interaction (PPI) network. Quantitative real-time polymerase chain reaction (qRT-PCR) was performed to verify the DEGs obtained by our integrated analysis. Compared with normal controls, 1732 DEGs in PMOP were obtained with p<0.05. According to the qRT-PCR results, expression of ATF2, FBXW7, RDX, and RBBP4 was consistent with that in our integrated analysis, generally. GO and KEGG enrichment analysis showed that those DEGs were significantly enriched in regulation of transcription, DNA-dependent, cytoplasm, protein binding, and MAPK signaling pathway. A total of 58 shared DEGs in PMOP versus normal control and in patients with PMOP versus patients after LWDH treatment were identified, which had opposite expression trend in these two comparisons. In the PPI network, CSNK2A1, ATF2, and FBXW7 were three hub proteins. Three genes including ATF2, FBXW7, and RDX were speculated to be therapeutic targets of LWDH for PMOP based on BATMAN-TCM database. We speculated that three genes of ATF2, FBXW7, and RDX may play crucial roles in both pathogenesis of PMOP and therapeutic mechanism of LWDH on PMOP. Our results may provide clues for the molecular pathogenesis of PMOP and offer new possibilities for treatment of PMOP.


Introduction
Postmenopausal osteoporosis (PMOP) often occurs due to the simultaneous interaction of independent predisposing factors, including aging and continuous calcium loss [1]. In the developed country, 9-38% of females are affected by PMOP. It is estimated that osteoporosis-associated fractures affect one-third of women in their lifetime, which is a public health concern [2]. Recently, both in animal models and in humans that were affected by osteoporosis newly, evidences of the relationship between immune system and bone have been accumulated. Estrogen is a well-known regulator of the immune system and T-cell functions. Bone loss induced by estrogen deficiency in menopause is a complex effect of a multitude of pathways and cytokines working in a cooperative fashion to regulate osteoclastogenesis and osteoblastogenesis [3][4][5][6]. Some hypotheses have been done by previous literature data. Recent studies have shown that obesity and osteoporosis share common genetic and environmental factors. Wnt family members and Wnt inhibitor molecules were involved not only in the regulation of differentiation of multipotent mesenchymal stem cells into osteoblasts and adipocytes but also in the coordination of the switch toward osteo-or adipogenesis fate within bone marrow [7].
Microarray studies have been used to explore the pathogenesis of PMOP. Expression profiling with larger sample size could be obtained by integrated microarrays analysis, which contributes to identifying more accurately differentially expressed genes (DEGs) of PMOP than an individual microarray. In this study, we performed an integrated analysis of five gene expression datasets to identify the DEGs between PMOP and normal controls and explore the molecular mechanism of PMOP. Functional annotation and PPI network construction were performed to explore the biological functions of these DEGs.
Various medications, including alendronate, etidronate, risedronate, and strontium ranelate, have been employed to prevent osteoporotic fragility fractures in patients with PMOP [8]. MF Faienza et al. demonstrated that denosumab and romosozumab have promise in the treatment of osteoporosis and their different mechanisms of action compared to existing antiosteoporotic drugs may permit alternative strategies for osteoporosis treatment down the line [9]. In recent years, research of traditional Chinese medicine has confirmed that incidence of osteoporosis was closely associated with kidney. Tonics for kidney, such as Liuweidihuang (LWDH) pills, were demonstrated to exert an effect on diabetes mellitus [10], breast carcinoma [11], and protection of testes in mice [12]. However, few studies were focused on the effect of LWDH on human osteoporosis. In the paper of Manyu Li [13], it is evaluated the effects of morroniside and loganin isolated from LWDH pills on the proliferation, differentiation, and apoptosis of MC3T3-E1 (mouse embryonic osteoblast precursor cells). Morroniside and loganin were found to promote the differentiation of osteoblast precursor cells and inhibit their apoptosis, thereby reducing bone resorption. Ji-rong Ge et al. found that CLCF1 is an important gene associated with Shen-yin deficiency of PMOP and the therapeutic effects of LWDH may be mediated by upregulated expression of CLCF1 and activation of the JAK/STAT signaling pathway [14]. In light of the potent curative effect of LWDH, exploring the candidate genes associated with both pathogenesis of PMOP and therapeutic mechanism of LWDH on PMOP could provide clues for the mechanism of PMOP and development of novel therapeutic strategies for PMOP.

Materials and Methods
. . Eligible Gene Expression Profiles of PMOP. We selected gene expression datasets of PMOP on the Gene Expression Omnibus database (GEO, http://www.ncbi.nlm.nih.gov/geo) which is the largest database of high-throughput gene expression data [15]. Search keywords in PMOP versus normal control were "postmenopausal osteoporosis" [MeSH Terms] OR "postmenopausal osteoporosis" [All Fields]) AND "gse" [Filter]. Datasets that meet the following criteria will be included in our study. (1) The selected dataset must be genomewide mRNA transcriptome data. (2) These data are derived from postmenopausal osteoporosis patients or blood samples (without drug stimulation or transfection) of patients with osteoporosis before and after LWDH pills treatment was considered. (3) Five sets of high-throughput mRNAs were screened by the standardization or original data set.
. . Identification of DEGs Associated with PMOP or LWDH Treatment. Firstly, we analyzed all the selected databases individually. Then, to minimize the heterogeneity among different datasets enrolled in integrated analysis, normalization and log2 transformation were performed for the raw data. Finally, an R package, metaMA package, was used to combine data from multiple microarray datasets. Individual p-values were calculated and multiple comparison correction FDR was obtained by using Benjamini and Hochberg method. Genes with FDR < 0.05 were screened out and regarded as DEGs [16].
. . qRT-PCR Confirmation for PMOP. Based on the results of high-throughput transcriptome data integration analysis, four differentially expressed genes were screened as candidate genes (PMOP versus normal control); blood samples (6 cases) from PMOP patients and healthy individuals were collected. Six blood samples were frozen at 80 ∘ C within 2 hours after blood withdrawal. After thawing the frozen samples at room temperature, 1 ml of the sample was used to perform RNA isolation with Trizol reagent (Invitrogen, China) according to the manufacturer's instructions. By using SuperScript5 III Reverse Transcriptase (Invitrogen, China), we generated cDNA from 1 g extracted RNA. With Power SYBR5 Green PCR Master Mix (Applied Biosystems, USA), we performed quantitative PCR in an ABI 7500 real-time PCR system. Relative gene expression was analyzed using 2 −ΔΔCt method. The human GAPDH gene was used as endogenous control for mRNA expression in analysis. The sequence of the primers used in RT PCR experiment was shown in Supplementary Table 1.

. . Functional Annotation of DEGs in PMOP and Common
DEGs. To identify the function and the potential pathway of DEGs in PMOP and 58 common DEGs, Gene Ontology (GO) classification (molecular functions, biological processes, and cellular component) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment were performed by using the online software GeneCodis (http://genecodis.cnb.csic.es/analysis) [17]. FDR < 0.05 was defined as the criteria of statistical significance. Fifty-eight genes were subjected to GO enrichment and KEGG enrichment analysis using the R language (GSEABase package), Pvalue < 0.01.

. . Protein-Protein Interaction (PPI) Network Construction.
To further research the biological functions of 58 common DEGs, we constructed the PPI network by using Biological General Repository for Interaction Datasets (BioGRID) (http://thebiogrid.org/) and Cytoscape [18]. Based on the existing protein interaction with opposite effects on PMOP versus normal control and LWDH data of BioGRID database, Cytoscape was used to search for common and regulatory genes treatment versus PMOP. After removing the two nondifferentially expressed genes, the protein network interaction map was drawn. The network consists of nodes and edges. The nodes in the network represented the proteins and the lines represented the interactions between protein [19].

. . Target Prediction of LWDH in BATMAN-TCM Database.
To define the key targets of LWDH, we retrieved BATMAN-TCM database, which is the first online Bioinformatics Analysis Tool for molecular mechanism of Traditional Chinese Medicine (http://bionet.ncpsb.org/batman-tcm/) [20]. LWDH was composed of five herbs, such as ZE XIE, SHAN YAO, MU DAN PI, SHU DI HUANG, JIU YU ROU, and FU LING. We input the herb list denoted by Pinyin name of 'ZE XIE, SHAN YAO, MU DAN PI, SHU DI HUANG, JIU YU ROU, and FU LING' with the following default parameters: predicted candidate targets (including known targets) with scores not smaller than score cutoff = 20 for each ingredient are presented and used for further bioinformatics analyses. The targets of Liuweidihuang were then obtained by a combination of the targets of ZE XIE, SHAN YAO, MU DAN PI, SHU DI HUANG, JIU YU ROU, and FU LING.

. . Differential Expression Analysis of DEGs in PMOP.
Four gene expression microarray datasets (GSE100609, GSE56815, GSE13850, and GSE7429) and a microarray dataset associated with therapeutic mechanism of LWDH on PMOP (GSE57273) were enrolled in our study. The clinical features of the enrolled subjects were shown in Supplementary Table 2. In total, the ratio of PMOP versus normal control was 54:60. Compared with the normal controls, 1732 DEGs in PMOP were obtained with p<0.05; among which, 918 genes were upregulated and 814 genes were downregulated. The hierarchical clustering heatmaps of DEGs in PMOP versus normal control were shown in Figure 1 . . qRT-PCR Confirmation for PMOP. To manifest the expression of integrated analysis for PMOP, we compared the results of qRT-PCR and integrated analysis of four genes including ATF2, FBXW7, RBBP4, and RDX. In Figure 4, the expressions of ATF2, FBXW7, RDX, and RBBP4 were all consistent with our integrated analysis. These 4 selected DEGs for qRT-PCR confirmation were derived from the top 10 degrees. Based on previous studies, ATF2, FBXW7, and RDX were osteoporosis-related DEGs. Saul D [21] found that RBBP4 (involved in age-related memory loss) was increased in tibia callus after ovariectomy.  . . PPI Network of Common DEGs. To identify potential interactions between DEGs, a PPI network was constructed. The results identified 264 nodes (genes) and 347 edges of DEGs in Figure 5. All points were differentially expressed in both groups, with blue ovals representing the common and opposite genes, diamonds representing DEGs in PMOP versus normal control, and inverted triangles representing DEGs in LWDH treatment versus PMOP, of which red represents upregulated mRNA and green represents downregulated mRNA. Among them, the higher degrees of genes were CSNK2A1 (degree=66), ATF2 (degree=33), FBXW7 (degree=26), RBBP4 (degree=23), MAP3K1 (degree=20), SRRM1 (degree =19), NCOA1 (degree=15), RDX (degree=11), and APP (degree=10), among which, the hub proteins were CSNK2A1, ATF2, and FBXW7.

Discussion
Postmenopausal osteoporosis (PMOP) is a multifactorial chronic disease. With the increase of the elderly population, the incidence of osteoporosis is on the rise, which is a health issue of concern in the world [6]. Since the effective treatment is deficient, exploring the mechanism and finding novel therapeutic strategies for PMOP have long been an essential goal. In this study, we obtained 1732 DEGs between PMOP and normal controls by integrated microarray analysis. Based on the function annotation analysis, we found that MAPK was a significantly enriched pathway for DEGs in PMOP. And previous study [22] has reported that MAPK signaling pathways are crucial in regulating osteogenesis, a genetic disorder affecting the bones. Lizhi Xing et al. [22] observed the upregulation of MAPK in the ovariectomyinduced bone loss. According to the functional annotation and PMOP-specific protein-protein interaction (PPI) network, four key DEGs (ATF2, FBXW7, RDX, and IGF2R) were speculated to be closely associated with PMOP.
ATF2 is a well-known osteoclastogenic transcription factor successfully identified by motif discovery analysis [23].
Osteoclast (also known as bone-resorbing cells) is a kind of bone tissue composition which is responsible for bone resorption. During bone resorption, the active osteoclasts secrete acidic substances and enzymes to break down the mineralized bone matrix and reduce the bone mass in this area [24]. In menopausal women, regulated by hormone, the activity of osteoclasts was accelerated so that the bone wall was thin. In our study, ATF2 is upregulated in PMOP group compared with normal controls, which supported previous studies. In addition, ATF2 is a hub protein of PMOP-specific PPI network which emphasized its importance in PMOP.
Radixin (RDX) is a cytoskeletal protein that may be important in linking actin to the plasma membrane. Among the ERM family members, only radixin is shown to be directly bound to the barbed ends of actin filaments. Previous studies have clarified the important roles of ERMs in various cellular processes such as epithelial morphogenesis, cell migration, Figure 5: Protein-protein network of 58 common DEGs. All points were differentially expressed in both groups; the blue ovals represented the common and oppositely expressed genes; the diamonds represented the differentially expressed genes in PMOP versus normal control; the inverted triangle represents the differentially expressed genes in LWDH treatment versus PMOP, of which red represented upregulated mRNA and green represented downregulated mRNA. and cell signal transduction [25]. Phosphorylated ERM protein plays an important role in cell growth, metabolism, and antiapoptosis and can promote cell migration and filamentous actin recombination and promote cell growth [26]. Upregulated RDX was found to be upregulated in PMOP compared with normal controls which were a potential regulator of PMOP.
FBXW7, an F-box protein component of the cullin-RING ubiquitin E3 ligase superfamilies, can recruit specific substrates to the RING E3 core through their variable carboxyterminal protein-interaction domains. FBXW7 mediates the polyubiquitination and proteasomal degradation of active Notch. The study of Mi Yang [27] demonstrated that endogenous levels of FBXW7 protein were regulated by overexpression or inhibition of miR-497B195. In addition, overexpression of miR-497B195 was found to promote osteogenesis in endothelial cells of mice. In our study, FBXW7 is upregulated in PMOP group compared with normal controls which was speculated to involve with PMOP regulated by miR-497B195. Further research was needed to confirm this finding.
And previous studies evaluated the relationship between serum insulin-like growth factor-binding protein profiles and bone mineral density measurements in PMOP [28]. Insulinlike growth factor-binding protein, also known as IGFBP, serves as a carrier protein for insulin-like growth factor, which binds to insulin-like growth factor 2 receptor (IGF2R) [29]. The serum insulin-like growth factor-binding protein-2 ratio in postmenopausal women with osteoporosis was significantly higher (p < 0.02) than that in normal healthy postmenopausal women, but the serum insulin-like growth factor binding protein-3 ratio in women with osteoporosis was significantly lower (p < 0.01). Wuster et al. suggested that there is a positive correlation between bone mineral density and IGFBP-3 in postmenopausal women with osteoporosis [30].
Furthermore, three upregulated genes (ATF, RDX, and FBXW7) were found to be downregulated in patients with PMOP after treatment of LWDH in this present study. Based on the BATMAN-TCM database, the potential targets of four herbs of "LWDH (ZE XIE, SHAN YAO, MU DAN PI, SHU DI HUANG and FU LING)" were predicted. All the above three genes (ATF, RDB, and FBXW7) were predicted targets of LWDH which provided evidence for our finding.

Conclusion
Taken together, MAPK was found to be associated with PMOP. Four genes (ATF, RDB, FBXW7, and IGF2R) might be potential regulators and biomarkers of PMOP. Moreover, three genes (ATF, RDB, and FBXW7) were speculated to be closely associated with both the pathogenesis of PMOP and the therapeutic mechanism of LWDH for PMOP. Our finding contributed to providing clues for exploring mechanism and novel strategy of drug design for PMOP. However, small sample size of RT-PCR was a limitation for our study, and further research with large sample size was needed to confirm this conclusion.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.