Introduction

Osteoporosis is a common disorder caused by the imbalance between osteoblastic bone formation and osteoclastic bone resorption [1]. Age is proportional to the risk of development of osteoporosis [2]. Due to the absence of estrogen, women after menopause were at the highest risk of osteoporosis according to a previous study [3]. Certainly, the absence of androgenic hormones could also cause osteoporosis in men [4]. Additionally, other factors also lead to osteoporosis, such as metabolic diseases, anorexia nervosa, thyroid and renal dysfunctions, and dietary as well as lifestyle habits such as low calcium intake or immobilization [5].

Osteoporosis can lead to a reduction in bone mass, deterioration in bone microarchitecture, susceptibility to skeletal fragility, and increased risk of fracture [6, 7]. Patients suffering osteoporotic fractures, particularly spinal fracture, are characterized by high morbidity and mortality, and significantly reduced quality of life. It has been reported that up to one-third of patients will sustain a new fracture within 5 years after the initial fracture [8]. Although anti-osteoporosis drugs reduce the risk of osteoporotic fractures by 20–70% in clinical trials depending on the drug and fracture type, persistence with osteoporosis therapy is poor, with one-year persistence ranging from 18% to 78% in real-world studies [914]. Postmenopausal women over age 55 are susceptible to osteoporosis-related spinal fracture [15]. Therefore, understanding the key mechanisms underlying osteoporotic fractures is crucial for treating spinal fracture.

The wingless-related integration site (Wnt)/β-catenin signaling pathway is the key pathway of bone metabolism to regulate bone mass. Defective Wnt signaling causes several monogenic skeletal disorders such as osteoporosis-pseudoglioma syndrome, van Buchem disease, and sclerosteosis [1618]. For example, WNT7B enhanced the ability of bone formation by increasing osteoblast activity to increase bone mass [19]. Glucocorticoids depressed bone formation by inhibiting the Wnt/β-catenin signaling pathway [20]. Laine et al. found that the mutation of Wnt1 could decrease the activity of the Wnt/β-catenin signaling pathway in bone, leading to a decreased number of bone cells, disrupted bone formation, low bone mass, and skeletal fragility [21]. If the skeletal fragility occurred in vertebrae, the spinal fracture might be caused by simple movements such as coughing or sneezing. Mäkitie et al. also reported that impaired WNT/β-catenin signaling progressively changed the spinal structures, which increased the risk of compression fractures, especially after the age of 50 [22].

In this study, the expression profiles of mRNA and miRNA after manipulation of Wnt signaling were obtained from the GEO database. Then, bioinformatic analysis including GO enrichment, KEGG enrichment, Reactome pathway, and PPI network analysis was performed to analyze the key pathway, miRNAs, and genes after mutation of the Wnt signaling pathway in osteoporosis and osteoporotic fractures. In the long term, our study should contribute to the treatment of osteoporotic fractures, especially spinal fracture.

Material and methods

Clinical samples

A total of 42 patients diagnosed with osteoporosis and spinal fracture between September 2018 and March 2020 and 45 age-matched healthy donors with unintentional spinal fractures (serving as normal controls) were enrolled in this study. Small needle bone biopsies from the spine were obtained from the subjects. All participants signed written informed consent forms before biopsy collection. This study has been approved by the Ethics Committee of the Affiliated Huai’an Hospital of Xuzhou Medical University and The Second People’s Hospital of Huai’an.

Data collection and array data analysis

Two expression profile data sets, GSE34747 and GSE103473, relating to osteoporosis and peripheral and spinal fracture were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/gds/). The differentially expressed genes (DEGs) of GSE34747 between Wnt activation samples (n = 3) and normal samples (n = 3) were identified using the impute R package. The differentially expressed miRNAs of GSE103473 between Wnt1 mutation samples (n = 12) and normal samples (n = 12) were identified using the impute R package. The DEGs and differentially expressed miRNAs were selected with the log |fold change| value ≥ 1 and P-value < 0.05. The shared genes between DEGs of GSE34747 and the target genes of miRNAs of GSE103473 were identified using Venny 2.1.0.

Construction and analysis of protein-protein interaction (PPI) network

To construct the PPI network, DEGs were uploaded to STRING (https://string-db.org/) and Metascape (http://metascape.org/gp/index.html#/main/step1). STRING is an online tool to predict and visualize the PPI, which includes direct and indirect associations. Metascape is an online gene annotation and analysis tool, which can analyze and visualize the PPI network. The analysis of PPI by the Metascape algorithm depends on BioGrid, InWeb_IM, and OmniPath databases. Molecular Complex Detection (MCODE) was applied to identify connected network components.

Gene ontology (GO) enrichment and pathway enrichment analysis

GO enrichment of DEGs including biological process, molecular function, and cellular component was analyzed using STRING and Metascape. The Reactome pathway database is a relational database of signaling and metabolic molecules. STRING PPI network construction was performed to analyze the Reactome pathways of DEGs. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways containing the information of the network of genes or molecules were also analyzed by STRING and Metascape.

Quantitative real-time PCR (RT-qPCR)

Total RNA was extracted from bone tissues from spines using TRIzol reagent (Invitrogen, USA) and quantified using NanoDrop 2000 (Thermo Fisher Scientific, USA). Then 2 µg of RNA was subjected to reverse transcription PCR to generate cDNA through the use of PrimeScriptVR RT reagent Kit (Takara, Japan). Then qRT-PCR was conducted to detect the expression of target genes using SYBR Premix Ex Taq (Takara, Japan). GAPDH was applied as the internal control, and the gene expression was calculated using the 2–ΔΔCT method. The measurement data were shown as mean ± standard deviation (SD), and the difference between two groups was analyzed by Student’s t-test using GraphPad Prism 8.0 (GraphPad Software, USA). Statistical significance was defined as p < 0.05.

Results

GO enrichment and Reactome pathway enrichment of 562 DEGs using STRING

Wnt signaling was associated with osteoporosis and could be activated by lithium. GSE34747 including the LiCl-stimulated samples (Wnt activation samples) and normal samples was analyzed by R software. Hierarchical clustering analysis showed that the datasets were well clustered: most genes in Wnt activation samples and normal samples tended to be grouped in two clusters, with minimal overlap between groups (Figure 1 A). For functional enrichment analysis by STRING, 562 DEGs were finally selected with the fold change value ≥ 1 and P-value < 0.05, and the complicated PPI network was displayed (Figure 1 B). GO enrichment showed that the biological process of 562 DEGs was associated with the chemical stimulus, the molecular function of 562 DEGs was associated with receptor activity, and cellular component was associated with membrane (Figure 1 C). Meanwhile, Reactome pathway analysis revealed that G-protein-coupled receptor (GPCR) was a key signaling pathway which has been proved to be related to osteoporosis. These results identified 562 DEGs that may be associated with membrane protein-related signaling pathways.

Figure 1

Functional enrichment analysis of 562 DEGs in GSE34747. A – Heat map of DEGs in GSE34747. Red color indicates upregulated genes while green indicates downregulated genes. B – PPI network for DEGs constructed using STRING. C – GO enrichment and Reactome pathway enrichment of DEGs analyzed by STRING

DEGs – differentially expressed genes, PPI – protein-protein interactions, GO – Gene Ontology, FDR – false discovery rate.

https://www.archivesofmedicalscience.com/f/fulltexts/125210/AMS-22-2-125210-g001_min.jpg

Analysis of process enrichment, pathway enrichment, and PPI network of DEGs using Metascape

To further identify the function of 562 DEGs, another algorithm, Metascape, was used to analyze and visualize the key processes and pathways. As shown in Figure 2 A, the calcium signaling pathway, consistent with the GO analysis results from STRING, was the key pathway. In addition, the PPI network constructed by Metascape displayed 5 MCODEs (Figure 2 B). The top 3 MCODEs were calcium signaling pathway, cAMP signaling pathway, and anterograde trans-synaptic signaling. The results from Metascape identified that the calcium signaling pathway and cAMP signaling pathway were the key pathways.

Figure 2

Analysis of GO enrichment, KEGG enrichment, and PPI network of 562 DEGs using Metascape. A – The top 20 enriched pathways and biological processes are shown in different colors. B – Construction of PPI network using Metascape. The top 3 MCODEs are displayed. DEGs, differentially expressed genes

PPI – protein-protein interactions, GO – Gene Ontology, KEGG – Kyoto Encyclopedia of Genes and Genomes, MCODE – Molecular Complex Detection.

https://www.archivesofmedicalscience.com/f/fulltexts/125210/AMS-22-2-125210-g002_min.jpg

Identification of overlapping genes between target genes of miRNAs and DEGs

GSE103473 was the miRNA profile of spinal fracture involving the Wnt1 mutation samples and normal samples. Hierarchical clustering analysis showed that most miRNAs in Wnt mutation samples and normal samples tended to be grouped in two clusters, while there was some degree of overlapping (Figure 3 A). miR-34a-5p, miR-22-3p, miR-143-5p, miR-18a-3p, miR-31-5p, and miR-223-3p were the top 6 differentially expressed miRNAs of GSE103473. Venny 2.1.0 was then used to select the overlapping genes between the target genes of the top 6 differentially expressed miRNAs and the formerly identified 562 DEGs of GSE34747. Due to the highest number of overlapping genes in miR-18a-3p, miR-18a-3p and 75 overlapping genes were identified as being associated with osteoporosis and spinal fracture (Figures 3 B–G).

Figure 3

Target genes of top 6 differentially expressed miRNAs were screened out using R software and Venny 2.1.0. A – Heat map of differentially expressed miRNAs in GSE103473 using R software. B–G – Venny 2.1.0 was applied to identify the target genes of miRNAs which also belonged to DEGs in GSE34747

DEGs – differentially expressed genes. -, low expression. +, high expression.

https://www.archivesofmedicalscience.com/f/fulltexts/125210/AMS-22-2-125210-g003_min.jpg

Identification of key genes using STRING and Metascape

To explore the biological functions of 75 genes involved in osteoporosis and peripheral and spinal fracture, STRING was first used to construct the PPI network of the 75 genes. The PPI network analysis showed that ADCY5, ADCY2, MLLT4, and GRIA1 were the genes associated with the cAMP signaling pathway (Figure 4 A). Similar to the result of STRING, process and pathway analysis by Metascape also revealed that the cAMP signaling pathway was the key pathway (Figure 4 B). Comparing the results of STRING with Metascape, ADCY2, ADCY5, and GRIA1 were identified as significant genes associated with the cAMP signaling pathway and possibly participate in osteoporosis and spinal fracture following Wnt signaling manipulation (Figure 4 C).

Figure 4

ADCY2, ADCY5, and GRIA1 were identified as key genes by STRING, Metascape and Venny 2.1.0. A – The PPI network of 72 target genes of miR-18a-3p was constructed using STRING. Genes involved in the cAMP signaling pathway are shown. B – cAMP signaling pathway was identified as the key pathway by Metascape analysis. Different colors represent different processes and pathways. C – Common genes (ADCY2, ADCY5, and GRIA1) from Metascape and STRING involved in cAMP signaling pathway identified using Venny 2.1.0. PPI, protein-protein interactions

https://www.archivesofmedicalscience.com/f/fulltexts/125210/AMS-22-2-125210-g004_min.jpg

Expression of ADCY2, ADCY5, and GRIA1 in osteoporotic patients with spinal fracture

To further investigate the association of ADCY2, ADCY5, and GRIA1 with spinal fracture, we collected bone tissues from osteoporotic patients with spinal fracture (n = 42) and age-matched healthy donors (n = 45), and examined the expression of ADCY2, ADCY5, and GRIA1 mRNA by RT-qPCR. The results showed approximately 50% downregulation of ADCY2 mRNA in the bone tissues of osteoporotic patients compared with the healthy donors (normal group) (Figure 5 A). Meanwhile, the expression of ADCY5 and GRIA1 mRNA exhibited 60% downregulation in the bone tissues of osteoporotic patients (Figures 5 B–C). These results indicated that ADCY2, ADCY5, and GRIA1 may play a regulatory role in the occurrence of osteoporosis with spinal fracture.

Figure 5

ADCY2, ADCY5, and GRIA1 mRNA expression in osteoporotic patients with spinal fracture. A–C – Relative expression of ADCY2 (A), ADCY5 (B), and GRIA1 (C) mRNA in bone tissues of osteoporotic patients with spinal fracture and healthy controls was detected by RT-qPCR. Data are shown as mean ± SD. Statistics analysis was conducted using the Mann-Whitney test. Normal: healthy control group

https://www.archivesofmedicalscience.com/f/fulltexts/125210/AMS-22-2-125210-g005_min.jpg

Discussion

Spinal fracture induced by osteoporosis has become a health burden of the aging population, especially in postmenopausal women over 55. Activation of the Wnt/β-catenin signaling pathway has been proved to prevent osteoblast and osteocyte apoptosis, so it plays a negative role in osteoporosis [20]. In this study, the GO analysis and Reactome pathway analysis revealed that 562 DEGs may be associated with membrane protein-related signaling pathways, especially the calcium signaling and cAMP signaling pathways. By analyzing the miRNA microarray, the 75 target genes of miR-18a-3p were screened out for further PPI network construction and GO term enrichments. Both STRING and Metascape enrichments identified that the cAMP signaling pathway was a crucial pathway. By comparing the results of STRING and Metascape, ADCY2, ADCY5, and GRIA1 were identified as key genes participating in osteoporosis and spinal fracture after the manipulation of Wnt signaling. Lastly, the aberrantly downregulated ADCY2, ADCY5, and GRIA1 in osteoporotic patients with spinal fracture further suggested the potential role of the three genes in the pathogenesis of osteoporotic spinal fracture.

The Wnt/β-catenin signaling pathway promoted the regeneration of osseous tissue by stimulating proliferation and differentiation of osteoblasts [23, 24]. The Wnt proteins are secreted glycoproteins, which can stimulate the signaling pathway by binding to LRP5/6 and the co-receptor Fizzled [25]. Then, receptors including Dsh, Axin, and APC can inhibit the activity of glycogen synthase kinase 3 (GSK3) to prevent the phosphorylation of β-catenin [26]. The phosphorylation of β-catenin results in the degradation of β-catenin so that the Wnt/β-catenin is inactivated. Kim et al. found that β-catenin expression in bone tissues from patients suffering from osteoporotic fractures was reduced, indicating that the decrease of β-catenin could cause osteoporotic fractures [24]. Additionally, the Wnt/β-catenin signaling pathway has also been reported to participate in the regulation of chondrocyte proliferation and apoptosis in osteoarthritis [27]. In our study, impaired Wnt signaling led to significant miR-18a-3p upregulation that possibly participated in osteoporosis and spinal fracture. Seventy-five genes that were both target genes of miR-18a-3p and DEGs caused by proactive Wnt activation underwent STRING and Metascape interrogation, which demonstrated that the cAMP signaling pathway was the key pathway that may be associated with Wnt activation or mutation. We concluded that activation or inactivation of the Wnt signaling pathway could affect cAMP signaling, consequently affecting osteoporosis and spinal fracture processes.

Cyclic 3′,5′-adenosine monophosphate (cAMP) is an important second messenger in bone homeostasis, which plays as a prominent role in determining the fate of cells. The intracellular cAMP level may be elevated by activating GPCR, a major mediator of bone remodeling, by inhibiting osteoblast apoptosis and enhancing osteoblast differentiation [28]. In the present study, the Reactome pathway analysis revealed that signaling by GPCR was closely associated with Wnt activation, which was consistent with the enrichment analysis results for 75 overlapping genes. cAMP has been proved to inhibit osteoblast proliferation by suppressing the MAP kinase pathway [29]. In addition, a network pharmacological study demonstrated that the cAMP signaling pathway was one of the critical pathways closely related to bone formation and resorption [30]. Weivoda et al. investigated the relationship between cAMP and Wnt pathways and found that Wnt3a suppressed osteoclast differentiation by activating the cAMP/PKA pathway [31]. Together with the previous studies and our bioinformatic analysis, we inferred that the cAMP signaling pathway is closely associated with Wnt signaling, and thus may participate in osteoporotic fractures.

Adenylate cyclase, a family of membrane bound enzymes that catalyze the formation of cyclic AMP from ATP under the stimulation of G-protein signaling, is a direct regulator of cAMP signaling pathway [32]. The adenylate cyclase family consists of nine members, named ADCY1–ADCY9 [33]. It has been shown that the blockage of cAMP/PKA/CREB signaling through inhibiting the activity of adenylate cyclase can repress icariin-induced osteogenesis [34]. Another report stated that the stimulation of adenylyl cyclase could activate cAMP-mediated MAPK signaling and induce the expression of Runx2 in osteoblasts to accelerate bone regeneration [35]. Given the potential regulation of cAMP signaling in osteoporosis, the key catalytic enzyme of cAMP signaling, adenylate cyclase, can be speculated to play a regulatory role in the development and progression of osteoporotic spinal fracture. Particularly, ADCY6 was demonstrated to promote the proliferation and differentiation of osteoblasts in osteoporotic rats through activating the Rap1/MAPK signaling pathway [36]. It was also once reported that ADCY3 had a positive effect on bone formation [37]. Interestingly, a study on ADCY5 knock-out mouse models suggested that ADCY5 protected the mice from bone density reduction and susceptibility to fractures of aging [38]. ADCY2 was also identified to be a potential regulator in osteoporotic spinal fracture by our study; however, there is no study supporting this. Based on the close relationship between ADCY2 and cAMP signaling, we believe that it is worth studying the effects of ADCY2 in osteoporotic spinal fracture. On the other hand, only glutamate ionotropic receptor AMPA type subunit 1 (GRIA1), belonging to a family of AMPA receptors, has been proved to be a tumor suppressor gene in human osteosarcoma [39]. Therefore, the identified genes ADCY2, ADCY5, and GRIA1 are of high significance with respect to osteoporosis and potentially spinal fracture, and thus need to be further studied, which may provide new therapeutic strategies for osteoporotic fractures.

During our search on gene expression profiling GEO data series in spinal fracture, we found only the GSE34747 dataset. The limited number of samples is indeed a drawback; however, no other profiling data are available, and we currently do not have sufficient funds to conduct our own gene profiling experiment. We will certainly consider conducting gene expression profiling of our own data in the future. In addition, we have provided the PCR verification results, which to some extent makes the results more reliable.

In conclusion, we identified that the cAMP signaling pathway was associated with the activation or inactivation of the Wnt signaling pathway in osteoporotic fractures. Meanwhile, ADCY2, ADCY5, and GRIA1 were associated with osteoporotic fractures involving the Wnt pathway and cAMP pathway, although studies on these genes in osteoporotic fractures remain limited. Our findings may provide novel therapeutic strategies for osteoporotic fractures.