p53 plays a central 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.

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 osteogenesis-related 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 (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.

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://www.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% CO₂ 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.

qRT-PCR analysis

qRT-PCR analysis was performed as described previously [27]. Briefly, Trizol reagent (ThermoFisher Scientific, MA, USA, #L15596026) was applied to extract total RNA. cDNA was then synthesized using a qPCR RT Master Mix kit (Toyobo, Osaka, Japan). Relative expression levels of mRNA were calculated using the 2-ΔΔCt method (Ct of GAPDH minus the Ct of the target genes) and normalized to the level of U6 mRNA determined as the 2-ΔΔCt. The primer sequences used were as follows: for p53, CAGCACATGACGGAGGTTGT (forward) and TCATCCAAATACTCCACACGC (reverse); for Col 1a1, TGGCAAAGATGGACTCAACG (forward) and TCACGGTCACGAACCACATT (reverse); for ALP, GCTCTGGAAAGTCCTTCAAAGC (forward) and TCTTCTTCCCTGGACACTGCC (reverse); for OCN, TCACACTCCTCGCCCTATTG (forward) and CTCCTGAAAGCCGATGTGGT (reverse); for Runx2, CTACTATGGCACTTCGTCAGGAT (forward) and ATCAGCGTCAACACCATCATT (reverse); and for GADPH, GGAAGCTTGTCATCAATGGAAATC (forward) and TGATGACCCTTTTGGCTCCC (reverse).

Western blotting analysis

Cell lysates were prepared using NETN buffer (20 mM Tris HCl [pH 8.0], 100 mM NaCl, 1 mM EDTA and 0.5% Nonidet P-40) and were resolved by SDS-PAGE. Proteins were transferred to PVDF membranes, which were blocked in 5% skim milk overnight at 4°C. The membranes were probed with primary antibodies and labeled with HRP-conjugated secondary antibodies (Aspen, Johannesburg, South Africa, #AS1058). Chemiluminescence detection systems (Canon, Tokyo, Japan, #LiDE110) were used to visualize the protein bands. The antibodies used were as follows: anti-collagen I (1:500, Abcam, MA, USA, #ab34710), anti-ALP (1:1,000, Abcam, MA, USA, #ab95462), anti-Osteocalcin (1:500, Abcam, MA, USA, #ab93876), anti-RunX2 (1:500, Abcam, MA, USA, #ab23981), anti-p53 (1:5,00, #ab131442, Abcam, MA, USA), and anti-GAPDH (1:10,000, Abcam, MA, USA, #ab37168).

Micro-CT analysis

All samples were scanned using a BRUKER SkyScan 1276 scanner microCT system (BRUKER, Karlsruhe, Germany) to provide images at 2400 views, 5 frames/view, 37 kV, and 121mA. The bone trabeculae selected were approximately 0.2 mm below the proximal humerus growth plate and measured 0.5 mm in length. After scanning, calluses were preserved at -80°C for miRNA extraction, qRT-PCR, and western blot analysis. Measurement parameters were as follows: bone volume (BV), total volume (TV), BV/TV, and BMD.

Statistical analysis

GraphPad Prism 8.0 (GraphPad Software, Inc.) was used for statistical analyses. The data are presented as the mean ± standard deviation (SD). Comparisons between two groups were made using Student’s t-test. Comparisons among three groups were made using one-way ANOVA with post-hoc Tukey’s test.