p53 plays a central role in the development of osteoporosis

Osteoporosis is a metabolic disease affecting 40% of postmenopausal women. It is characterized by decreased bone mass per unit volume and increased risk of fracture. We investigated the molecular mechanism underlying osteoporosis by identifying the genes involved in its development. Osteoporosis-related genes were identified by analyzing RNA microarray data in the GEO database to detect genes differentially expressed in osteoporotic and healthy individuals. Enrichment and protein interaction analyses carried out to identify the hub genes among the deferentially expressed genes revealed TP53, MAPK1, CASP3, CTNNB1, CCND1, NOTCH1, CDK1, IGF1, ERBB2, CYCS to be the top 10 hub genes. In addition, p53 had the highest degree score in the protein-protein interaction network. In vivo and in vitro experiments showed that TP53 gene expression and serum p53 levels were upregulated in osteoporotic patients and a mouse osteoporosis model. The elevated p53 levels were associated with decreases in bone mass, which could be partially reversed by knocking down p53. These findings suggest p53 may play a central role in the development of osteoporosis.


INTRODUCTION
Osteoporosis is a metabolic disease characterized by decreased bone mass per unit volume, despite the bone tissue having normal calcification and a normal ratio of calcium salt and matrix [1]. Osteoporosis can occur in both genders at any age, but it is most common in postmenopausal women. Overall, the disease affects 30% of women and 12% of men at some point in their lifetimes [2], but it affects more than 40% of postmenopausal women [3]. Patients with osteoporosis usually have bone pain and are prone to having a fracture [1]. As one of the most commonly occurring chronic diseases among the elderly, osteoporosis has become a serious problem for public health care systems [4]. The main therapeutic strategies for treating osteoporosis are the use antiresorptive and anabolic drugs. However, most of the drugs in these classes have limitations and adverse side effects [3].
The pathogenesis of osteoporosis has not yet been fully elucidated. Factors that inhibit osteogenesis, promote bone resorption, or cause bone microstructural destruction may play a role in the development of osteoporosis, and a variety of genes may be directly or indirectly involved. For example, mutations of COLI are responsible for osteogenesis imperfecta [5], while AGING allelic variation of TCIRGI are significantly associated with low bone mineral density (BMD) in a certain population [6]. In addition, Xie et al. also suggested that CTNNB1 and TP53 may play crucial roles in primary osteoporosis [7]. Osteoporosis is thus potentially associated with multiple genes and may result from gene-environment interactions [2].
Bioinformatics technology has been used to integrate and analyze big data from public database repositories for several diseases. For instance, bioinformatic methods have demonstrated prevalent alterations in RNA methylation regulators across cancer types. It was concluded that the m 6 A regulators tightly correlate with the activation and inhibition of cancer pathways, and also correlate with prognostically relevant tumor subtypes [8].
In the present study, we applied similar bioinformatic analysis of an osteoporosis microarray dataset retrieved from the Gene Expression Omnibus (GEO) to explore the mechanism underlying osteoporosis. Identification and validation of differentially expressed genes (DEGs) suggest that p53 may play a key role in the development of osteoporosis.

Identification of DEGs
The GSE100609 dataset was obtained from the GEO database. It included gene expression profiles from 4 healthy individuals and 4 osteoporotic patients. Analysis of the dataset using the Morpheus online tool revealed 509 (228 upregulated and 281 downregulated) genes that were differentially expressed between the healthy group and the osteoporotic patients. The top 30 upregulated and downregulated genes are shown in Figure 1.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses
GO term analysis and KEGG pathway enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics tool. The results showed that under the biological processes category, upregulated DEGs in osteoporotic patients were significantly enriched in "Regulation of locomotion," "Regulation of cellular component movement," "Regulation of cell motility," "Anatomical structure formation involved in morphogenesis," and "Movement of cell or subcellular component." On the other hand, the DEGs downregulated in osteoporotic patients were enriched in "Axon extension," "Regulation of cellular component movement," "Neuron projection extension," "Developmental growth involved in morphogenesis," and "Positive regulation of cellular protein metabolic process" (Table 1). The top five KEGG pathways for the DEGs upregulated in osteoporotic patients were "cancer pathway," "small cell lung cancer pathway," "p53 signaling pathway," "Wnt signaling pathway," and "rap1 signaling pathway." The top five KEGG pathways for the DEGs downregulated in osteoporotic patients were "axon guidance pathway," "bacterial invasion of epithelial cells pathway," "African trypanosomiasis pathway," "Alzheimer's disease pathway," and "calcium signaling pathway" ( Table 2 and Figure 2).

Protein-protein interaction (PPI) network construction and module analysis
A PPI network was constructed using STRING ( Figure 3), and the top 10 hub genes with the highest degrees of interaction were determined with Cytoscape. These hub genes were TP53, MAPK1, CASP3, CTNNB1, CCND1, NOTCH1, CDK1, IGF1, ERBB2, and CYCS (Table 3). Gene centrality is represented by indexes of degree, closeness, and betweenness. Degree represented the degree of association of one node with all other nodes in the network. Closeness represented the closeness between a node and other nodes in the network. Betweenness was the number of times a node acted as the shortest bridge between two other nodes. The results for degree, closeness and betweenness are showed in Figure 4. Among the 10 hub genes, TP53 had the highest degree score of 174, indicating it potentially plays an important role during the development of osteoporosis.
The genes involved in the six PPI network modules with MCODE scores ≥5 and >5 nodes (Table 4) were characterized through GO term and KEGG pathway enrichment analyses. In the biological process category, the GO terms "Cellular response to chemical stimulus," "Positive regulation of cellular metabolic process," and "Positive regulation of macromolecule metabolic process" were significantly enriched. In addition, in the molecular function category, "Enzyme binding," "Receptor binding," and "Macromolecular complex binding" were also significantly enriched. Among KEGG pathways, "signaling pathway of cancer pathway," "proteoglycans in cancer pathway," and "p53 signaling pathway" were enriched (Table 5). In Figure  5, the top 10 genes associated with the enriched GO term/KEGG pathway are illustrated using a chord diagram. Expression of the top 50 genes from the 6 modules and their positions on chromosomes are also shown in Figure 6.   Alzheimer's disease pathway 7 5.2E-2 hsa04020 Calcium signaling pathway 7 6.7E-2 *The top 5 terms were selected based on P-values when more than five terms enriched terms were identified in a category.

Downregulating p53 expression may protect against osteoporosis in vitro and in vivo
Because TP53 had the highest degree score among the top 10 hub genes, we hypothesized that p53 plays an important role during osteoporosis development. Consistent with our hypothesis, TP53 gene expression was increased two-fold in the serum samples from osteoporotic patients compared to healthy controls ( Figure 7B). Correspondingly, levels of p53 protein were also significantly increased in osteoporotic patients ( Figure 7A), suggesting a potential role for p53 in osteoporotic development and/or progression.
In cultured human mesenchymal stem cells (hMSCs) transfected with siRNA targeting p53, expression levels of both p53 mRNA and protein were significantly decreased as compared to cells treated with PBS (control) or a negative control siRNA (siRNA-NC) ( Figure 7D, 7E). As expected, expression levels of the osteogenesisrelated genes Collagen Type I Alpha 1 (Col 1a1), Alkaline phosphatase (ALP), Osteocalcin (OCN) and RUNX Family Transcription Factor 2 (RunX2) were all significantly higher in p53 knockdown cells than in control (PBS) cells or cells expressing siRNA-NC ( Figure 7F, 7G).
The effect of p53 on osteoporosis was further investigated in vivo using a murine model. Osteoporosis was induced by oophorectomy the model mice treated with PBS (control), siRNA-p53 or siRNA-NC. Consistent with the in vitro experiments, bone volume       *The top 3 terms were selected based on P-values when more than 3 enriched terms were identified in a category.
(BV), total volume (TV), the BV/TV ratio and BMD were all significantly higher in the siRNA-p53 group than the control or siRNA-NC group ( Figure 8). Thus, downregulating p53 appears to suppress features of osteoporosis both in vitro and in vivo.

DISCUSSION
Osteoporosis is a metabolic disease characterized by low bone mass and microstructural destruction of bone tissue, resulting in increased bone fragility and proneness to fracture [1]. The pathogenesis of osteoporosis has not yet been fully elucidated. Factors that inhibit osteogenesis, promote bone resorption, and/or cause bone microstructural destruction can contribute to the development of osteoporosis.
The normal homeostasis of mature bone is maintained largely through the coordinated actions of regulating hormones and local cytokines [9]. Bone tissue continuously absorbs old bone and generates new bone, thereby maintaining a stable state of bone turnover [10]. With increasing age, however, the rate of bone turnover gradually declines, leading to gradual declines in BMD and bone mineral content (BMC) [11]. In general, the annual rate of BMD loss is about 0.5%, but the rate of BMD loss is usually slower in older men than women, in part because women experience a rapid loss of estrogen [12]. No such rapid hormonal change occurs in men [12]. The BMD loss is accompanied by distortion and destruction of bone microstructure. Some structures such as the trabeculae are unable to maintain normal morphology, resulting in narrowing, thinning, bending, and/or misalignment of trabecular bone, and even micro-injury or microfractures [13]. In some of the worse cases, the bone is completely absorbed and cavities develop. In those areas, the cortical bone becomes thinner, numbers of trabeculae decrease, and fragility increases until spontaneous compression fractures (e.g. in vertebral bodies) or transverse fractures (e.g. in the femoral neck or distal radius) occur [14].

AGING
Osteoporosis is a multifactorial disease [15] and a variety of genes and signaling pathways participate its pathogenesis [16]. In the present study, we analyzed the interactions among involved proteins and ranking the top 10 hub genes. inhibits osteogenesis by affecting the function of MSCs via miRNA signaling pathways [20]. Xie et al. found that CTNNB1 and TP53 are the most upregulated DEGs and may play a crucial role in primary osteoporosis [7]. It was also reported by that signaling in the MAP kinase pathway is differentially activated in MSCs derived from osteoporotic postmenopausal women [17], and that miR-378 overexpression attenuates high glucosemediated suppression of osteogenic differentiation by targeting CASP3 and activating the PI3K/Akt pathway [23]. In addition, Notch appears to suppresses osteoblast differentiation, at least in part by restricting glucose metabolism [18], and inhibiting microRNA-139-5p promotes osteogenic differentiation of BMSCs by activating the Wnt/beta-catenin signaling pathway by AGING targeting NOTCH1 [21]. CDK1, which is essential for osteoblast proliferation, functions as a molecular switch that shifts osteoblast proliferation to maturation [22]. IGF-1 is a key regulator of longitudinal skeletal growth and is an anabolic hormone that influences bone modeling and remodeling throughout life [19]. Finally, sustained activation of ErbB1 and ErbB2 enhances hMSC osteogenesis [24]. All these findings are consistent with the bioinformatics analysis reported in the present study.
P53 had the highest degree score in the PPI network, indicating it may play a central role during the development of osteoporosis. P53 is encoded by the tumor suppressor gene TP53 and suppresses tumor growth by slowing cell growth and division. In the present study, it is found that serum p53 levels are increased in osteoporosis patients, and knocking down p53 partially reversed decreases in BMD in vitro and in vivo. In addition, GO and KEGG enrichment analyses indicate that p53 is involved in "the cancer pathway," "proteoglycans in cancer pathway," and "P53 signaling pathway." We therefore suggest that p53 may contribute to the pathogenesis of osteoporosis via these pathways.
There are a few limitations to this study. First, only one dataset was found for the bioinformatics analysis. Analysis of additional datasets would strengthen the evidence from bioinformatics analyses. Second, we did not perform a subgroup analysis of different osteoporosis subtypes, which may involve different signaling pathways and genes. A future study will focus on investigating how p53 interacts with other osteoporosisrelated pathogenic genes in different osteoporosis subtypes.
In summary, our findings suggest that p53 may play a key role in the development of osteoporosis, and that suppressing some of the activities of p53 may inhibit the development and/or progression of osteoporosis. For that reason, we suggest p53 should be considered a potential therapeutic target for the treatment of osteoporosis.

Search of microarray data and identification of DEGs
The microarray dataset GSE100609 for osteoporosis patients and control subjects were retrieved from the NCBI Gene Expression Omnibus (GEO). The expression data were then uploaded to and analyzed by Morpheus (Morpheus, https://software.broadinstitute.org/ morpheus), a versatile online matrix visualization and analysis software package. The gene expression was visualized using a heatmap. Signal to noise >1 or signal to noise <1 is considered as DEGs.

GO and KEGG pathway enrichment analysis
The DAVID version 6.8 online bioinformatics resources tool was used to perform GO term and KEGG pathway enrichment analyses [25,26].

PPI network analysis
The PPI network was used to further explore the functional interactions between DEGs imported into Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org). Interactions with a combined score >0.5 were identified. The PPI network was then built using Cytoscape software (version 3.7.2). Modules of the PPI networks with a Molecular Complex Detection (MCODE) score ≥5 and >5 nodes were identified using the MCODE plugin in Cytoscape. Function enrichment analysis of DEGs in the top module was performed with DAVID. The DEGs were then ranked based on degree centrality using the Centiscape 2.2 plugin in Cytoscape.

Blood collection
Between May 2016 and June 2018, peripheral blood samples were collected from patients in Shanghai Tongji Hospital (10 healthy volunteers, 10 osteoporosis patients) for gene expression analysis. The patient studies were approved by the Committees of Clinical Ethics in the Tongji Hospital (Tongji University of Medicine, Shanghai, China), and informed consent was obtained from all participants.

Animals and specimen collection
Thirty female C57BL/6 mice (6 weeks old) were obtained from the Department of Animal Science of Tongji Hospital, Tongji University School of Medicine. A mouse osteoporosis model is established using oophorectomy. Following oophorectomy, the 30 mice were randomly divided into three groups: a control group treated with PBS (n=10), a siRNA-NC group treated with negative control siRNA (n=10), and a siRNA-p53 group treated with p53 siRNA (n=10). We administered subcutaneously injected 5 ml of solution per kilogram body weight every 2 days for 10 weeks. The weights of the animals were recorded weekly during the experimental period. All animal experiments adhered to the National Institutes of Health Guide for the Care and Use of Laboratory Animals and were performed according to the protocols approved by the Medical Research Ethics Committee of the Tongji Hospital, Tongji University of Medicine.

Cell culture and transfection
hMSCs were kindly donated by the Tongji Hospital, Tongji University of Medicine, Shanghai, China. Cells were maintained for a maximum of 4 passages in specific medium designed for human mesenchymal stem cells at 37°C in a 5% CO2 incubator. Cells were transfection with 50 nM siRNA-NC or siRNA-p53 (GenePharma, Shanghai, China) using Lipofectamine 3000 (ThermoFisher Scientific, MA, USA, #L3000001) following the manufacturers protocol.