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.

1. 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-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.

2. Materials and Methods

2.1. 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 genome-wide 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.

2.2. 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].

2.3. 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 SuperScript® III Reverse Transcriptase (Invitrogen, China), we generated cDNA from 1μg extracted RNA. With Power SYBR® 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.

2.4. 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), P-value < 0.01.

2.5. 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].

2.6. 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.

3. Results

3.1. 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. Likewise, there were 1350 DEGs (322 upregulated and 1027 downregulated) in patients with PMOP versus patients after LWDH treatment. The hierarchical clustering heatmaps of DEGs in patients with PMOP versus patients after LWDH treatment were shown in Figure 2.

Figure 1

The heatmap of DEGs in PMOP versus normal control. Row and column were used to represent DEGs and datasets. The color scale represented the expression levels of DEGs.

Figure 2

The heatmap of DEGs in LWDH treatment versus PMOP. Row and column were used to represent DEGs and datasets. The color scale represented the expression levels of DEGs.

3.2. Functional Annotation of DEGs in PMOP

In Figure 3, GO and KEGG enrichment analysis showed that DEGs were significantly enriched in the following terms. Biological processes: regulation of transcription, DNA-dependent (p=4.31E-22), and signal transduction (p=4.44E-16). Cell component: cytoplasm (p=4.14E-108) and nucleus (p=1.06E-95). Molecular function: protein binding (p=4.19E-115) and metal ion binding (p=9.37E-27). KEGG pathway enrichment results showed that DEGs were enriched in the Osteoclast differentiation (p=0.00931908) and MAPK signaling pathway (p=2.58E-09).

Figure 3

Go and KEGG functional enrichment of DEGs in PMOP. (a) Biological progress; (b) cellular component; (c) molecular function; and (d) KEGG functional enrichment. The y-axis represented GO and KEGG terms and the x-axis shows represented counts of DEmRNAs enriched in GO and KEGG terms.

3.3. 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.

Figure 4

qRT-PCR verification results for PMOP. The X-axis indicated eight selected DEGs and the Y-axis indicated log 2 (fold change) for qRT-PCR results.

3.4. Coanalysis of PMOP versus Normal Control and LWDH Treatment versus PMOP

Combining the results of PMOP versus normal control and LWDH treatment versus PMOP, we obtained 58 common DEGs; among which, 43 genes were upregulated in PMOP versus normal control and downregulated in LWDH treatment versus PMOP; 15 genes were upregulated in LWDH treatment versus PMOP and downregulated in PMOP versus normal control. The details of these 58 common DEGs were shown in Table 1.

Table 1

Coexpressed and opposite trend of DEGs in the two groups.

Id	Symbol	case vs normal		before vs after
		Raw P-Value	Up-down	P-Value	Up-down
2181	ACSL3	0.044463305	up	5.81E-05	down
23519	ANP32D	0.031671352	up	0.00020335	down
1386	ATF2	0.015288615	up	0.000105124	down
892	CCNC	0.037220624	up	4.23E-05	down
1370	CPN2	0.009788448	up	1.57E-05	down
1457	CSNK2A1	0.029560249	up	0.00035395	down
79848	CSPP1	0.011683224	up	1.11E-05	down
6374	CXCL5	0.027209429	up	6.23E-05	down
80777	CYB5B	0.000178896	up	0.000622171	down
23741	EID1	0.043185211	up	0.000455124	down
9669	EIF5B	0.014971999	up	0.00042694	down
23085	ERC1	0.032754729	up	0.000160588	down
30001	ERO1L	0.0154858	up	0.000625767	down
55294	FBXW7	0.037818674	up	0.000138841	down
2914	GRM4	0.023673494	up	0.000502142	down
2935	GSPT1	0.004988472	up	2.63E-05	down
3482	IGF2R	0.043016644	up	0.000342073	down
7850	IL1R2	0.025284329	up	3.93E-06	down
9929	JOSD1	0.013589301	up	0.000203472	down
51621	KLF13	0.020717689	up	7.15E-05	down
4214	MAP3K1	0.022243844	up	8.40E-05	down
79649	MAP7D3	0.003259237	up	0.000418951	down
158747	MOSPD2	0.001561148	up	0.000121954	down
7703	PCGF2	0.013528324	up	6.33E-05	down
5294	PIK3CG	0.001627058	up	0.00030047	down
59338	PLEKHA1	0.049102137	up	0.000452568	down
59339	PLEKHA2	0.007489418	up	1.26E-05	down
133619	PRRC1	0.022743976	up	0.00052731	down
5928	RBBP4	0.00010116	up	4.70E-05	down
5962	RDX	0.049953292	up	0.000486787	down
55680	RUFY2	0.001194843	up	0.000266002	down
5552	SRGN	0.009208736	up	0.000209717	down
6717	SRI	0.037181628	up	0.000508307	down
10250	SRRM1	0.000952409	up	6.51E-05	down
9262	STK17B	0.016461631	up	3.91E-05	down
26136	TES	0.000594943	up	0.000259901	down
10959	TMED2	0.046585899	up	0.000567562	down
54765	TRIM44	0.000191712	up	0.000565748	down
389898	UBE2NL	0.013171719	up	0.00048327	down
54904	WHSC1L1	0.020186972	up	1.24E-05	down
81555	YIPF5	0.018233854	up	0.000196894	down
10472	ZNF238	0.002670727	up	7.55E-05	down
7639	ZNF85	0.028282272	up	0.000212913	down
55811	ADCY10	0.012001458	down	0.000502024	up
55009	C19orf24	0.006292312	down	0.000610074	up
51372	CCDC72	0.000667113	down	0.00060282	up
2323	FLT3LG	0.017012385	down	0.000637474	up
2794	GNL1	0.012495179	down	0.000573231	up
55716	LMBR1L	0.001157605	down	0.000600187	up
8648	NCOA1	0.016190136	down	0.000192176	up
51070	NOSIP	0.008546637	down	0.000285444	up
54623	PAF1	0.034489911	down	0.000236501	up
51177	PLEKHO1	0.036853037	down	0.000331093	up
5439	POLR2J	0.026864334	down	0.000213443	up
6165	RPL35A	0.041712982	down	0.00055368	up
6827	SUPT4H1	0.016431549	down	0.000265733	up
54958	TMEM160	0.022504594	down	0.000255735	up
389136	VGLL3	0.049059686	down	0.000604032	up

3.5. Functional Annotation of Common DEGs

In Supplementary Figure 1, several GO categories of up- and downregulated common DEGs were enriched in GO terms associated with regulation of biological process (p=9.84E-05) and biological regulation (p=0.000669575) in biological process (BP), intracellular part (p=0.003479728) and intracellular part (p=0.0023611) in cellular component (CC), and protein binding (p=0.00011324), enzyme binding (p=0.00052695), and catalytic activity (p=0.003357531) in molecular functions (MF). The most significantly different biological process (BP) of upregulated and downregulated DEGs was the regulation of biological process. And the most significantly different molecular functions (MF) were the protein binding.

3.6. 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.

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.

3.7. Predicted Gene Targets of LWDH

Using the BATMAN-TCM database, the potential mRNA targets of the six herbs of LWDH (ZE XIE, SHAN YAO, MU DAN PI, SHU DI HUANG, JIU YU ROU, and FU LING) were predicted. There were 177, 478, 73, 42, 0, and 503 targets for ZE XIE, SHAN YAO, MU DAN PI, SHU DI HUANG, JIU YU ROU, and FU LING, respectively. Among them, there were 10 common mRNA targets (ESR1, FGFR2, MED1, PGR, PRKCB, PTGS1, PTGS2, TRIM24, VDR, and WNT4) for the six herbs. Compared with the DEGs after treatment of LWDH, ATF2, FBXW7, and RDX were significantly differentially expressed in patients after treatment, and they are also the gene targets of LWDH. And, in the PPI network, ATF2, FBXW7, and RDX were among the genes which have high degree.

4. 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 ovariectomy-induced 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, 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 carboxy-terminal 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]. Insulin-like 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.

5. 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.