Genomic Comparison of Mechanisms Underlying Zuogui Pill and Estradiol Valerate in Prevention and Treatment of Postmenopausal Osteoporosis

， -， -， -， -

(1. Department of TCM, Shenzhen Longgang Central Hospital， Shenzhen, 518116, China; 2. Suzhou Hospital of Traditional Chinese Medicine Affiliated to Nanjing University of Chinese Medicine, Suzhou, 215003, China)

ABSTRACT: OBJECTIVE To compare the potential genes and pathways involved in the preventive treatment of ovariectomized (OVX)-induced osteoporosis in bone marrow mesenchymal stem cells (BMSCs) obtained from OVX rats treated with Zuogui Pill and Estradiol valerate. METHODS A total of 36 rats were divided into four groups: Model, Sham (Control), Zuogui Pill, and Estradiol valerate. The rats in the last two groups were intragastrically administered with Zuogui Pill or Estradiol valerate, respectively, from 2 weeks post operation to 4, 8, and 12 weeks post operation. The rats in the other two groups were treated with distilled water. The gene expression patterns in BMSCs were investigated, which were subjected to bioinformatics analyses, including identification of differentially expressed genes (DEGs), pathway enrichment and correlation analyses, weighted correlation network analysis (WGCNA), and protein-protein interaction (PPI) network construction. RESULTS At 12 weeks post operation, bone mineral density in OVX model group significantly decreased. Meanwhile, bone mineral densities in Zuogui Pill and Estradiol valerate groups were significantly higher than that in the model group, indicating that the models were constructed successfully, and both traditional Chinese medicine and western treatment could improve the bone density of OVX rats. Numerous DEGs were identified in different groups at different time points. Difference analysis showed that compared with the model group, 405, 403 and 618 DEGs were screened in Zuogui Pill treatment group at 4, 8 and 12 weeks after surgery, respectively, and 78, 326 and 232 DEGs in Estradiol valerate group. Pathway analysis showed that these DEGs were closely related to inflammation- and immunity-associated pathways, such as herpes simplex infection (Il1b and Tnf), the PI3K-Akt signaling pathway (Itga2 and Fgf21), and leukocyte transendothelial migration (Cxcl12). WGCNA identified several key modules that were significantly linked to the PI3K-Akt signaling pathway and inflammatory bowel disease. The PPI network in the Zuogui Pill group included 388 nodes and 133 nodes in the Estradiol valerate group. CONCLUSION The results of this study suggest the importance of inflammation and immunity in the pathogenesis of OVX-induced osteoporosis in bone marrow microenvironment. The identified pathways and genes may serve as therapeutic targets for this disease.

KEY WORDS: Zuogui Pill; Estradiol valerate; inflammation and immunity; postmenopausal osteoporosis; genetics

Osteoporosis is a common and complex skeletal disorder characterized by low bone mass and microarchitectural deterioration of bones because bone resorption occurs more rapidly than bone formation, ultimately leading to bone fragility and fractures[1]. Osteoporotic fractures impose a heavy social and economic burden on patients[2]. Moreover, osteoporosis has become one of the most prevalent bone diseases in postmenopausal women, namely postmenopausal osteoporosis (PMOP)[3]. Researchers have been aware of the link between menopause and osteoporosis since the 1960s[4].

For any disease, prevention is better than cure. Therefore, it is important to develop preventive treatment for PMOP. It is well known that estrogen plays a fundamental role in maintaining skeletal homeostasis, and widely accepted that the development of PMOP in women with ovarian failure is associated with estrogen withdrawal[5]. Thus, estrogen treatment has become the standard therapy for bone loss in postmenopausal women via increasing bone density[4]. Estrogen-containing hormone replacement therapy has been proven to be effective in reversing the impact of menopause on bone density and preventing fractures[6]. However, this treatment is unsatisfied due to its unwanted side effects. Specifically, it increases the risk of venous thrombosis, stroke, breast cancer, and possibly dementia[6-10]. Hence, it is necessary to explore alternative strategies for the prevention and treatment of PMOP.

Traditional Chinese medicine considers that the essence of life is stored in the kidneys, which dominate bones and produce bone marrow. If the kidney essence is insufficient, bone marrow formation will be lacking, leading to fragile bones[11]. Therefore, kidney tonifying and essence supplement may shed some light on the prevention of PMOP. Zuogui Pill is a well-known traditional Chinese prescription used for invigorating the kidneys and nourishing the essence; it has been demonstrated to be efficacious in preventing and treating osteoporosis induced by insufficient kidney essence[12]. However, the underlying molecular mechanisms behind this have not been fully investigated and require further studies.

In this study, we aimed to compare the mechanisms underlying the efficacy of Zuogui Pill and Estradiol valerate in ovariectomized rats in terms of important target cells of osteoporosis, namely, bone marrow mesenchymal stem cells, at the whole-genome expression level. Estradiol valerate and Zuogui Pill administrations were performed once the OVX model was established . Then, we extracted gene expression profiles from their BMSCs. Finally, bioinformatics analyses were employed to evaluate the associated DEGs and pathways.

1 Materials and methods

1.1 Animals

A total of 36 specific-pathogen-free female Sprague-Dawley rats (7 months old), weighing (305±45) g, wereused in this study. Rats were provided by the Center for Comparative Medicine of Nantong University (Jiangsu, China) [SCXK (Su) 2008-0010].

1.2 Drugs

Estradiol valerate, purchased from Bayer Healthcare (Germany), was dissolved in double distilled water (0.01 mg/mL) after grinding. Zuogui Pill was obtained from Nantong University Affiliated Hospital (Jiangsu, China)， which was composed of 24 g ofRadixRehmanniaePreparata(Shu Di), 12 g ofRhizomaDioscoreae(Shan Yao), 12 g ofFructusLycii(Gou Qi Zi), 12 g ofFructusCorni(Shan Yu Rou), 12 g ofSemenCuscutae(Tu Si Zi), 12 g ofCollaPlastiTestutinis(Gui Ban Jiao), 12 g ofCollaCornusCervi(Lu Jiao Jiao), and 9 g ofRadixAchyranthisBidentatae(Chuan Niu Xi). All Chinese herbs were routinely decocted with about 1 000 mL water and finally condensed into 105 mL (1 g crude drug/mL).

1.3 Grouping and modeling

All rats were randomly assigned to one of the following groups: OVX model, sham operation, western medicine (Estradiol valerate), and traditional Chinese medicine (Zuogui Pill) (n=9 for each group). Additionally, each group was further divided into three subgroups: 4, 8, and 12 weeks (drugs adminstrated from 2 weeks post operation to 4, 8, and 12 weeks post operation). The rats in the OVX model, western medicine, and traditional Chinese medicine groups underwent bilateral ovariectomy to induce bone mass loss with estrogen reduced. In contrast, the rats in the sham operation group underwent the removal of the same amount of fat tissue around the ovaries.

1.4 Administration method and sample collection

The rats in the sham and model groups were intragastrically administered with distilled water from 2 weeks post operation to 4, 8, and 12 weeks post operation. In contrast, the rats in the western medicine and traditional Chinese medicine groups were intragastrically administered with 0.1 mg/kg of Estradiol valerate and 5 g/kg of Zuogui Pill, respectively. The rats were killed at 4, 8, and 12 weeks post operation, after which their femurs and tibiae were removed in accordance with previously reported procedures[13]. BMSCs was confirmed by morphology and flow cytometry assays. BMSCs were cultured to the F2generation using the whole bone marrow adherent method. Moreover, flow cytometry assay was used to detect positive and negative markers (CD90, CD45) of BMSCs. Total RNA was extracted from BMSCs using Trizol (Invitrogen, USA). After the quality test (RIN≥8 and 28S/18S>0.7), microarray analysis was carried out according to the protocol from Agilent-014879 Whole Rat Genome Microarray 4×44K G4131F (Biochip National Engineering Research Center of Shanghai, China).

1.5 Bone mineral density measurement

To confirm the success of modeling as well as effects of Estradiol valerate and Zuogui Pill on bone mineral density, lumbar vertebrae 4-6, the left distal femur, and the proximal tibia of the rats in the four groups were sampled at 12 weeks post operation. Bone mineral density was determined using dual-energy X-ray absorptiometry.

1.6 Data preprocessing

Original data were subjected to background correction (normexp method) and normalization (quantile method) using limma[14](http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) in R package. Then, probes were mapped to gene symbols based on the annotation file. When multiple probes corresponded to the same gene symbol, their mean value was served as the expression value of the gene. Finally, a expression matrixes containing 17 071 genes was obtained.

1.7 Differential expression analysis

Genes that were differentially expressed in Estradiol valeratevs. OVX model pair, and Zuogui Pillvs. OVX model pair at three time points were analyzed.P-value was calculated applying the Bayesian approach in the limma package and adjusted using the BH method. An adjustedP-value < 0.05 and |log2FC| ≥ 0.58 were considered as cut-off values.

1.8 Venn diagram of DEGs

Venn diagram analysis was performed to study the differences and similarities of DEGs in different groups by inputting the DEGs and theirlog2FC values into VennPlex[15](http://www.irp.nia.nih.gov/bioinformatics/vennplex.html), which could filtrate important genes quickly.

1.9 Pathway enrichment analysis

The DEGs of the same drug intervention at different time points were combined together. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways tightly associated with the union DEGs were analyzed using clusterProfiler[16]. ClusterProfiler calculated theP-value of gene enrichment using the geometric distribution method in the KEGG.db[17]package.P-value < 0.05 were considered to be significantly enriched.

1.10 Pathway correlation analysis

Pathways that had strong connections to the union DEGs and had at least one gene in common with other pathways were screened out and used to construct a correlation network that reflected the functional correlations among pathways.

1.11 Weighted gene correlation network analysis

WGCNA[18]is a system biology method used to identify clusters and modules of highly correlated genes, and to summarize the intra-modular hub gene as well as calculate module membership measures. WGCNA is also frequently applied to analyze time-series expression profiles, obtaining modules closely related to expression changes over time. In the present study, the DEGs identified from the comparisons of Zuogui Pillvs. OVX model and Zuogui Pillvs. control groups were combined together and subjected to WGCNA. The obtained significant co-expression module (|co-expression coefficient| ≥ 0.6 andP-value < 0.05) was considered to be closely related to Zuogui Pill treatment. Similarly, the Estradiol valerate treatment co-expression module was analyzed.

1.12 Protein-protein interaction network construction

Based on the DEGs involved in the significant co-expression modules, interaction pairs between these DEGs were predicted using the Search Tool for the Retrieval of Interacting Genes (STRING, version 10.0, http://www.string-db.org/)[19]database (score>0.4). The PPI network was then constructed using Cytoscape[20]. The hub genes in the PPI network were detected by calculating the node scores applying a random walk algorithm[21].

2 Results

2.1 Morphological identification and flow cytometry analysis of BMSCs

Morphological identification of BMSCs is shown in Figure 1A-D. Flow cytometry analysis demonstrated that the CD90 positive expression rate of F2 generation cells was higher than 95%, while the CD45 positive expression rate was lower than 5% (Figure 1E), indicating a high level of BMSCs in F2 generation cells.

A. morphology under 50× microscope; B. morphology under 100× microscope; C. mineralized nodules of osteogenic induction; D. adipocytes of lipogenic induction; E. Positive expression rates of CD90 and CD45 in BMSCs F2 cells

Figure1MorphologicalidentificationandflowcytometryanalysisofBMSCs

2.2 Bone mineral density

At 12 weeks post operation, bone mineral densities in lumbar vertebrae 4-6, the left distal femur, and the proximal tibia of rats in the OVX model group was significantly lower than those in the sham operation group. Compared with the sham operation group, no significant difference was found in the Estradiol valerate group, whereas a significant decrease in bone mineral density was detected in the proximal tibia of the rats in the Zuogui Pill group. Meanwhile, bone mineral densities in the Zuogui Pill and Estradiol valerate groups were significantly higher than those in the model group (Table 1).

2.3 Differential expression analysis

After preprocessing, hundreds of DEGs were identified from the comparisons between different groups at different time points. In general, there were more DEGs in the Zuogui Pill group than those in the Estradiol valerate group. Moreover, in the Zuogui Pill group, DEGs at 12 weeks were the most numerous, while in the Estradiol valerate group, the number of DEGs peaked at 8 weeks. The results of this analysis are shown in Table 2. The heatmaps of DEGs in different groups are shown in Figure 2.

Table 1 The bone mineral density of rats in three groups at 12 weeks post-operation

Note：&P<0.05,&&P<0.01 compared with sham operation group;*P<0.05,**P<0.01 compared with model group.

2.4 Venn diagram analysis

The similarities and differences of gene expressions between the Zuogui Pill and Estradiol valerate groups are presented in Venn diagrams (Figure 3). Comparing the Zuogui Pillvs. OVX model pair with the Estradiol valeratevs. model pair, 43, 118, and 61 common genes were identified at 4, 8, and 12 weeks post operation, respectively.

Table 2 The number of differentially expressed genes in different comparison groups at three time points

Note： Up count and down count represent the number of up-regulated or down-regulated genes.

2.5 Pathway enrichment analysis

The pathways strongly linked to the DEGs in the Zuogui Pill and Estradiol valerate groups are shown in Figure 4A and 4B, respectively. As shown in Figure 4A, the DEGs in the Zuogui Pill group were significantly associated with the following pathways: p53 signaling pathway, herpes simplex infection, hematopoietic cell lineage, and neuroactive ligand-receptor interaction. As shown in Figure 4B, the DEGs in the Estradiol valerate group were mainly involved in pathways of leukocyte transendothelial migration, apoptosis, PI3K-Akt signaling, and inflammatory bowel disease (IBD).

2.6 Pathway correlation analysis

The pathway correlation networks of the Zuogui Pill and Estradiol valerate groups were constructed (Figure 5). As shown in Figure 5A, the NOD-like receptor signaling pathway had a high degree of connectivity and interacted with numerous pathways, such as the herpes simplex infection and influenza A pathways. As shown in the network in Figure 5B, pathways in cancer, PI3K-Akt signaling, and hepatitis B had large number of interactions with other pathways.

2.7 WGCNA results

After WGCNA, seven significant modules with high co-expression coefficients were identified from the Zuogui Pill group (Figure 6A), among which six were particularly associated with KEGG pathways, including the Hippo signaling pathway-multiple species, and PI3K-Akt signaling pathway (Figure 6B).

In addition, five significant modules with higher co-expression coefficients were identified from the Estradiol valerate group (Figure 6C), among which four modules were significantly associated with pathways including amoebiasis and primary immunodeficiency (Figure 6D).

2.8 PPI network analysis

Based on the DEGs involved in the significant modules of the Zuogui Pill group, a PPI network with 849 interaction pairs and 388 nodes was constructed (Figure 7A). In this network, RAN-binding protein 2 (Ranbp2) had the highest degree of connectivity (Table 3). Additionally, Ras-related C3 botulinum toxin substrate 2 (Rac2), bone morphogenetic protein 4 (Bmp4), and aldehyde dehydrogenase 1 family member B1 (Aldh1b1) also had high degrees. The PPI network of the Estradiol valerate group comprised 156 interaction pairs and 133 nodes (Figure 7B), with the hub nodes of PBX homeobox 1 (Pbx1), tumor necrosis factor (Tnf), protein phosphatase with EF-hand domain 1 (Ppef1), and others.

3 Discussion

It has been suggested that BMSCs play important roles in the formation, maintenance, and repair of bone tissue[22]. BMSCs can differentiate into osteoblasts, the malfunction of which was implicated in the development of osteoporosis[22]. In the present study, the bone mineral density results revealed that Zuogui Pill induced less effective therapeutic effect than Estradiol valerate on the proximal tibial bone mineral density， which may attribute to their different regulatory genes and pathways. Based on the expression profiles obtained from BMSCs in OVX rats, numerous DEGs were identified in comparisons between the Zuogui Pill group and Estradiol valerate group at different time points. These DEGs were significantly associated with inflammation- and immunity-related pathways, including herpes simplex infection and the PI3K-Akt signaling pathway. Furthermore, WGCNA identified several key modules that significantly connected to the PI3K-Akt signaling pathway and IBD.

Gene expression regulation is a dynamic process. Therefore, it is necessary to identify and characterize changes in gene expression over time[23-24]. Compared with the screening of DEGs at a single time point, gene arrays (time series analysis), which monitor changes in expression patterns over time, are able to capture multidimensional dynamics of complex biological systems, hence, reveal genes associated with a disease or phenotype and explore their roles in time dimension[24]. Therefore, the present study conducted microarray experiments to identify DEGs in time series. The results showed that in the Zuogui Pill group, hundreds of DEGs were identified at three time points, suggesting the rapid regulatory effect of Zuogui Pill on PMOP. Interestingly, the highest number of DEGs was identified at 12 weeks post operation, which suggested that Zuogui Pill exerts efficacy over time. In the Estradiol valerate group, DEGs at 8 and 12 weeks post operation were far more numerous than those at 4 weeks post operation, and the number of DEGs peaked at 8 weeks. It has been documented that the gene regulation of rat BMSCs might be the most complex and diversified at 8 weeks after castration. In this study, pathway analysis showed that in Estradiol valerate group, the DEGs at week 8 were significantly enriched in inflammation-related pathways, such as hepatitis B and PI3K-Akt signaling. This result implied that Estradiol valerate may exert its anti-inflammatory functions at an earlier stage, which further suggested that the therapeutic intervention should not be given after the 8th week.

Table 3 The top 15 gene nodes with high degrees in the PPI network

Recently, an increasing number of studies have indicated that inflammation has a major effect on bone metabolism, resulting in increases in resorption and fracture risk[25]. Numerous pro-inflammatory cytokines have been implicated in the regulation of osteoblasts and osteoclasts, suggesting that inflammation significantly contributes to the etiopathogenesis of osteoporosis[26]. In the present study, pathway enrichment analysis revealed that DEGs identified from the comparison of the Zuogui Pill group with the OVX model group were significantly associated with the pathway of herpes simplex infection, such asIl1bandTnf. Herpes simplex virus is maintained in a latent state in ganglia by CD8+T cells which are necessary for effective immunity[27].Il1bandTnfare cytokines involved in normal inflammatory and immune responses[28]. The fact that the immune system and bone are intimately linked through essential interactions is important in this context[29]. Roodman et al. reported that increased cytokine production in the bone microenvironment expands the osteoclastic pool and improves the activity of mature osteoclasts[30]. Therefore, Zuogui Pill may exert therapeutic effects in the prevention and treatment of OVX-induced PMOP through inflammation- and immunity-associated pathways and genes.

In the current study, seven modules in the Zuogui Pill group, which were significantly involved in the PI3K-Akt signaling pathway, were identified via WGCNA, . The PI3K-Akt signaling pathway is frequently disturbed in human cancers, and this disturbance plays an important role in tumor development[31]. Different stages of tumor development are reported to be closely related to inflammatory responses[32]. The PI3K-Akt signaling pathway has been found to regulate genes associated with inflammation and immunity[33]. Moreover, a variety of drugs have been developed to inhibit inflammation targeting the PI3K-Akt signaling pathway[34]. In our study, integrin subunit alpha 2, a molecule critical for the function of leukocytes in inflammation[35], was found to be enriched in this pathway. Noticeably, the PI3K-AKT signaling pathway has recently been reported to be involved in the regulation of osteoporosis[36]. Combining these findings together, Zuogui Pill may alleviate OVX-induced PMOP through the PI3K-Akt signaling pathway.

In this study, the PI3K-Akt signaling pathway was also strongly associated with the DEGs identified in the comparison of the Estradiol valerate group and the model group (at 8 weeks post operation), such asFgf21. Fgf21 is a newly identified member of the FGF family and implicated in energy balance[37]. Wei et al.[38]reported that Fgf21 overexpression reduced bone mass in mice. Fgf21 was also found to regulate BMSC differentiation. BMSCs are multipotent marrow stromal cells that have the ability to differentiate into various cell types required for tissue regeneration, including chondrocytes and osteoblasts[38]. The malfunction of BMSCs has been reported to be essential for osteoporosis[22,39]. In our study,Fgf21 was downregulated in the Estradiol valeratevs. OVX model group pair, suggesting that it is a target gene for Estradiol valerate induced efficacy in PMOP.

In the comparison of the Eestradiol valerate group and the model group at 12 weeks post operation, DEGs such as C-X-C motif chemokine ligand 12 were significantly linked to inflammation- and immunity-associated pathways including leukocyte transendothelial migration.Cxcl12 encodes a chemotactic cytokine, which belongs to the chemokine family, the members of which can activate leukocytes and are often induced by pro-inflammatory stimuli. Chemokines have been demonstrated to play important roles in bone remodeling and the regulation of osteoclast functions[40]. Yu et al. reported that chemokines contributed to inflammation-induced bone resorption by stimulating the recruitment of osteoclast progenitor cells, the fusion of which forms multinucleated osteoclasts[41]. Thus,Cxcl12 may be a candidate gene in leukocyte transendothelial migration accounted for the preventive treatment of Estradiol valerate in PMOP.

In the Estradiol valerate group, four significant modules were identified by WGCNA, which were associated with the pathways of primary immunodeficiency and amoebiasis. Primary immunodeficiencies are a group of disorders which affect distinct components of the innate and adaptive immune systems, inducing part of the body’s immune system function abnormally[42]. This study found that the pathway of amoebiasis was particularly associated withIl1band Arginase 2 (Arg2).Il1bwas also differently expressed in the Zuogui Pill groupvs. OVX model group pair;Il1bhas been found to play an important role in osteoclast action via a range of pathways[25].Arg2 can compete with inducible nitric oxide synthase for its substrate, and the balance between these two enzymes plays an important role in the regulation of immune responses[43]. AlthoughArg2 has not been reported to be connected to osteoporosis, we speculated thatArg2 andIl1bare candidate genes of Estradiol valerate-produced effects on PMOP given their roles in immune responses in BMSCs. In conclusion, the present study identified numerous inflammation- and immunity-associated pathways and genes in BMSCs in the Zuogui Pill and Estradiol valerate groups, further indicating the importance of inflammation and immunity in the pathogenesis of PMOP at the target cell level. These pathways and genes may serve as preventive and therapeutic targets for this disease.

However,owing to constraints regarding research funding and time, this study focused only on changes at the genetic level. The corresponding proteins, direct participants involved in various in vivo biological processes, also warrant intensive studies. Therefore, in the future, we intend to explore the biological effects of these corresponding proteins with the goal of exploring PMOP and the drug treatment mechanism more comprehensively. It is worth noting that the protein levels and gene expressions may not be completely paralleled, which could complicate our goal of obtaining a comprehensive understanding of the changes at the genetic and protein levels.