Research Paper Volume 12, Issue 10 pp 9500—9514
MiR-16-5p regulates postmenopausal osteoporosis by directly targeting VEGFA
- 1 Department of Orthopedic Surgery, Tongji Hospital, Tongji University School of Medicine, Shanghai 200065, China
- 2 Department of Orthopedic Surgery, Brigham and Women’s Hospital, Harvard Medical School, Boston, MA 02115, USA
- 3 Department of Orthopedics, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan 430022, China
Received: February 8, 2020 Accepted: March 31, 2020 Published: May 19, 2020https://doi.org/10.18632/aging.103223
How to Cite
Copyright © 2020 Yu et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
In this study, we used bioinformatics tools, and experiments with patient tissues and human mesenchymal stem cells (hMSCs) to identify differentially regulated genes (DEGs) and microRNAs (miRNAs) that promote postmenopausal osteoporosis. By analyzing the GSE56815 dataset from the NCBI GEO database, we identified 638 DEGs, including 371 upregulated and 267 downregulated genes, in postmenopausal women with low bone density. Enrichment and protein-protein interaction network analyses showed that TP53, RPS27A, and VEGFA were the top three hub genes with the highest degree of betweenness and closeness centrality. TargetScanHuman and DIANA software analyses and dual luciferase reporter assays confirmed that miR-16a-5p directly targets the 3’UTR of VEGFA. Postmenopausal patients with osteoporosis showed higher miR-16-5p and lower VEGFA levels than those without osteoporosis (n=10 each). VEGFA levels were higher in miR-16-5p knockdown hMSCs and were reduced in miR-16-5p-overexpressing hMSCs. mRNA expression of osteogenic markers, ALP, OCN, and RUNX2, as well as calcium deposition based on Alizarin red staining, correlated inversely with miR-16-5p levels and correlated positively with VEGFA levels. These findings suggest that miR-16-5p suppresses osteogenesis by inhibiting VEGFA expression and is a promising target for postmenopausal osteoporosis therapy.
Osteoporosis is characterized by reduced bone mass and weakened bone micro-architecture, resulting in an increased risk of fractures . Osteoporosis is a common age-related disease in post-menopausal women and the elderly . It affects the quality of life of nearly 200 million people worldwide and is a significant burden on the public healthcare systems . Nearly 40% of women suffer from osteoporosis and sustain fractures of the hip, spine, or the forearm during their lifetime . While there are several medications to treat osteoporosis based on their symptoms and severity, there is no gold standard treatment established, as yet . Therefore, it is important to understand the genetic mechanisms underlying osteoporosis so that new therapeutic targets can be identified . Recent studies have identified SQRDL and PPWD1 genes as risk factors that are associated with osteoporosis [6, 7].
miRNAs are small noncoding RNAs between 21–25 nucleotides in length that regulate protein expression through post-transcriptional gene silencing [8–10]. Previous studies have identified differentially expressed miRNAs in osteoporosis patients, including miR-21, miR-133a, miR-152-3p, miR-30e-5p, miR-140-5p, miR-324-3p, miR-19b-3p, miR-335-5p, miR-19a-3p, miR-550a-3p, and miR-422a [11–13]. Moreover, several signaling pathways that are involved in osteoporosis are regulated by microRNAs such as miR-133a, miR-218, miR-618, miR-27a, miR-214-5p, and miR-203 [14–18]. However, the mechanistic details of how these miRNAs influence osteoporosis are not known.
Bioinformatics is an interdisciplinary field that combines computer science and biology to research, analyze and interpret large sets of biological data. Computational tools are routinely used to model biological processes, predict disease mechanisms and generate experimental hypotheses . Bioinformatics tools are also widely used to screen potential genes related to several human diseases, which are then validated by comprehensive follow-up experimental studies [20, 21].
In this study, we analyzed microarray data from the GSE56815 dataset to identify DEGs that are related to osteoporosis. Furthermore, we performed enrichment analysis and constructed an interaction network of the DEGs to identify hub genes. Furthermore, we used the TargetScanHuman and DIANA software to identify microRNAs that potentially target the hub genes and validated these findings in patient tissues and in vitro experiments using human mesenchymal stem cells.
DEGs in postmenopausal osteoporosis
We analyzed the microarray data from the GSE56815 dataset in the NCBI GEO database, which includes gene expression data on the Affymetrix Human Genome U133A Array (GPL96) platform of 20 postmenopausal women each with high or low bone mineral density (BMD). We identified 638 DEGs, including 371 upregulated and 267 downregulated genes in patients with low BMD compared to those with high BMD (Figure 1B). The heat map showing hierarchical clustering of all DEGs and the volcano plot of the DEGs are shown in Figure 1.
Figure 1. Differentially expressed genes in postmenopausal osteoporosis patient samples. (A) The heat map shows hierarchical clustering of differentially expressed gene expression in postmenopausal patients with or without osteoporosis (n=20 each) using the GSE56815 dataset. (B) Volcano plot shows differentially expressed genes postmenopausal patients with or without osteoporosis (n=20 each) using the GSE56815 dataset.
Construction of the protein-protein interaction (PPI) network of the DEGs
Figure 2 shows the top 100 DEGs from the 3 modules and their chromosomal positions. The top 10 hub genes were TP53, RPS27A, VEGFA, MAPK8, CDC42, CREBBP, SIRT1, RPL35A, RPL30, and SNRPG (Table 1). Figure 3 shows the degree, betweenness, and closeness centrality of the top 10 hub genes. The top three hub genes with the highest degree, betweenness and closeness centrality are TP53, RPS27A, and VEGFA.
Table 1. Degree of the top 10 genes in the protein-protein interaction network.
|Gene ID||Gene name||Degree|
|TP53||Tumor protein p53||91|
|RPS27A||Ribosomal Protein S27a||60|
|VEGFA||Vascular Endothelial Growth Factor A||58|
|MAPK8||Mitogen-Activated Protein Kinase 8||42|
|CDC42||Cell Division Cycle 42||40|
|CREBBP||CREB binding protein||38|
|RPL35A||Ribosomal Protein L35a||31|
|RPL30||Ribosomal Protein L30||30|
|SNRPG||Small Nuclear Ribonucleoprotein Polypeptide G||30|
Figure 2. Circular visualization of chromosomal positions and connectivity of the top 100 genes. The circular map represents all the chromosomes and the lines from each gene point to their specific chromosomal locations. The names of the DEGs are shown in the outer circle and the three hub genes are shown in red. The red and green colors in the outer heatmap represent DEGs with high and low degree, respectively. The pink and blue colors in the middle heatmap represent DEGs with high and low betweenness, respectively. The brown and black colors in the inner most heatmap represent DEGs with high and low closeness, respectively.
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis
Table 2 and Figure 4 shows the results of the GO and KEGG pathway enrichment analysis of the hub genes and DEGs using Database for Annotation, Visualization and Integrated Discovery version (DAVID). The top 3 GO terms related to biological processes were cellular response to hypoxia, cellular response to decreased oxygen levels, and cellular response to oxygen levels. The top 3 GO terms related to molecular functions were P53 binding, heterocyclic compound binding, and organic cyclic compound binding. The top 3 GO terms related to cellular components were cytosolic ribosome, ribosomal subunit, and cytosolic part. The top three KEGG pathways were pancreatic cancer pathway, pathways in cancer, and renal cell carcinoma pathway.
Table 2. Top 3 GO terms according to biological process, molecular function, and cellular component, and top 3 KEGG pathways related to the top 10 hub genes in module.
|(A) Biological processes|
|GO:0071456||Cellular response to hypoxia||5||2.6E-7||VEGFA, CREBBP, TP53, SIRT1, RPS27A|
|GO:0036294||Cellular response to decreased oxygen levels||5||3.2E-7||VEGFA, CREBBP, TP53, SIRT1, RPS27A|
|GO:0071453||Cellular response to oxygen levels||5||4.2E-7||VEGFA, CREBBP, TP53, SIRT1, RPS27A|
|(B) Molecular functions|
|GO:0002039||P53 binding||3||6.5E-4||CREBBP, TP53, SIRT1|
|GO:1901363||Heterocyclic compound binding||9||2.9E-3||CDC42, RPL35A, RPL30, CREBBP, TP53, MAPK8, SIRT1, RPS27A, SNRPG|
|GO:0097159||Organic cyclic compound binding||9||3.2E-3||CDC42, RPL35A, RPL30, CREBBP, TP53, MAPK8, SIRT1, RPS27A, SNRPG|
|(C) Cellular component|
|GO:0022626||Cytosolic ribosome||3||2.0E-3||RPL35A, RPL30, RPS27A|
|GO:0044391||Ribosomal subunit||3||4.5E-3||RPL35A, RPL30, RPS27A|
|GO:0044445||Cytosolic part||3||7.1E-3||RPL35A, RPL30, RPS27A|
|(D) KEGG pathway|
|hsa05212||Pancreatic cancer||4||4.4E-5||CDC42, VEGFA, TP53, MAPK8|
|Hsa05200||Pathways in cancer||5||6.1E-4||CDC42, VEGFA, CREBBP, TP53, MAPK8|
|hsa05211||Renal cell carcinoma||3||2.4E-3||CDC42, VEGFA, CREBBP|
|KEGG, Kyoto Encyclopedia of Genes and Genomes; Top 3 terms were selected according to P-value.|
Figure 4. Functional and signaling pathway enrichment analyses of the top 10 genes in the protein-protein interaction network. The GO and KEGG pathway analysis of the top 10 genes identified by the PPI network are shown. The top 3 GO terms related to biological process are cellular response to hypoxia (GO:0071456), cellular response to decreased oxygen levels (GO:0036294), and cellular response to oxygen levels (GO:0071453). The top 3 GO terms related to molecular functions are P53 binding (GO:0002039), heterocyclic compound binding (GO:1901363), and Organic cyclic compound binding (GO:0097159). The top 3 GO terms related to cellular component are cytosolic ribosome (GO:0022626), ribosomal subunit (GO:0044391), and cytosolic part (GO:0044445). The top 3 KEGG pathways are pancreatic cancer (hsa:05212), pathways in cancer (hsa:05200), and renal cell carcinoma (hsa:05211).
Construction of the network of hub genes and related microRNAs
Next, we used TargetScanHuman  and DIANA  software to identify miRNAs that target the hub genes. As shown in Figure 5, hsa-miR-34a-5p potentially targets TP53 and VEGFA, and hsa-miR-15a-5p potentially targets TP53 and RPS27A. VEGFA is highly expressed in osteoblast precursors  and induces osteogenic differentiation [25, 26], which suggests its potential role in bone turnover. Previous reports suggest that miR-16-5p, miR-20b-3p, miR-15a-5p, and miR-34a-5p regulate osteogenesis [27–30]. QRT-PCR analysis of serum samples showed significantly high miR-16-5p levels and significantly reduced VEGFA levels in the osteoporosis patients compared with the healthy subjects (Supplementary Figure 1 and Figure 7A). Therefore, we selected the miR-16-5p/VEGFA axis to validate further.
MiR-16-5p is upregulated in osteoporosis patients and exerts negative effect on osteogenic differentiation
QRT-PCR analyses showed that miR-16-5p levels were significantly higher in the post-menopausal osteoporosis patients compared to postmenopausal women without osteoporosis (n=10 each; Figure 6A). Furthermore, we analyzed miR-16-5p and VEGFA levels in control hMSCs and hMSCs transfected with agomiR-NC, agomiR-16-5p, antagomiR-NC, and antagomiR-16-5p. The miR-16-5p levels inversely correlated with VEGFA protein expression compared to the corresponding controls (Figure 6B and Figure 7D). Moreover, the expression of osteogenic genes, ALP, OCN, and RUNX2 were significantly increased in the antagomiR-16-5p group, but significantly reduced in the agomiR-16-5p group (Figure 6C). This suggests that miR-16-5p suppresses osteogenesis. Furthermore, alizarin red staining showed increased calcium deposition in the antagomiR-16-5p group compared to the corresponding controls (Figure 6D–6E). This suggests that miR-16-5p expression decreases bone density.
Figure 6. High miR-16-5p levels inhibit osteogenic differentiation in osteoporosis patients and hMSCs. (A) QRT-PCR analysis of miR-16-5p levels in postmenopausal patients with or without osteoporosis (n=10 per group). (B) QRT-PCR analysis of miR-16-5p levels in control hMSCs and hMSCs transfected with agomiR-NC, agomiR-16-5p, antagomiR-NC, and antagomiR-16-5p is shown. (C) QRT-PCR analysis shows the relative expression of osteogenic marker genes, ALP, OCN and RUNX2, in control hMSCs and hMSCs transfected with agomiR-NC, agomiR-16-5p, antagomiR-NC, and antagomiR-16-5p. (D, E) Alizarin red staining shows calcium deposition in control hMSCs and hMSCs transfected with agomiR-NC, agomiR-16-5p, antagomiR-NC, and antagomiR-16-5p for 21 days. Scale bar = 10 mm. Note: The data are represented as means ± SD. *p < 0.05, **p < 0.01, ***p < 0.001.
Figure 7. miR-16-5p directly targets VEGFA. (A) QRT-PCR analysis of VEGFA levels in postmenopausal patients with or without osteoporosis (n=10 per group). (B, C) Dual luciferase reporter assay results show firefly luciferase activity relative to Renilla luciferase activity in hMSCs transfected with dual luciferase vectors containing VEGFA-WT-3’UTR or VEGFA-MUT-3’UTR and miR-16-5p agonist (agomiR-16-5p) or negative control (agomiR-NC). (D) QRT-PCR analysis shows relative VEGFA mRNA levels in control hMSCs and hMSCs transfected with agomiR-NC, agomiR-16-5p, antagomiR-NC, and antagomiR-16-5p. (E) Western blotting analysis shows relative VEGFA protein levels expression in hMSCs transfected with siRNA-NC, and siRNA-VEGFA. (F) QRT-PCR analysis shows relative expression of osteogenic marker genes, ALP, OCN and RUNX2, in hMSCs transfected with siRNA-NC, and siRNA-VEGFA. (G, H) Alizarin red staining shows calcium deposition in hMSCs transfected with siRNA-NC, and siRNA-VEGFA for 21 days. Scale bar = 10 mm. The data are means ± SD. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
MiR-16-5p directly targets VEGFA
QRT-PCR analysis showed that VEGFA mRNA levels were significantly reduced in osteoporosis patients compared to healthy subjects (n=10 per group; p < 0.001; Figure 7A). This suggests that VEGFA levels correlate with osteoporosis progression. Next, we performed the dual luciferase reporter assay to confirm if miR-16-5p specifically binds to the 3’UTR of VEGFA mRNA. The results show high luciferase activity in the miR-16-5p plus Mutant type (Mut)-3’UTR-containing luciferase plasmid group, but, luciferase activity is significantly reduced in the miR-16-5p plus Wild type (WT)-3’UTR-containing luciferase plasmid group (Figure 7B–7C). These results demonstrate that miR-16-5p specifically targets VEGFA. Next, we analyzed control hMSCs and hMSCs transfected with agomiR-NC, agomiR-16-5p, antagomiR-NC, and antagomiR-16-5p. QRT-PCR analysis shows that VEGFA mRNA levels were significantly reduced in the agomiR-16-5p group and significantly increased in the antagomiR-16-5p group (p < 0.001; Figure 7D). These results confirm that miR-16-5p promotes osteoporosis by suppressing VEGFA expression.
Next, we transfected hMSCs with empty vector control or VEGFA-specific siRNAs (siRNA-NC or siRNA-VEGFA). Western blot analysis showed that VEGFA protein expression was significantly reduced in the siRNA-VEGFA group compared to the siRNA-NC group (Figure 7E). QRT-PCR analysis showed that the expression of osteogenic genes, ALP, OCN and RUNX2, was significantly reduced in the siRNA-VEGFA group compared to the siRNA-NC group (Figure 7F). This suggests that VEGFA promotes osteogenesis. Moreover, Alizarin red staining showed significantly reduced calcium deposition in the siRNA-VEGFA group compared to the siRNA-NC group (p < 0.001 Figure 7G–7H). This further suggests that VEGFA positively regulates bone density by promoting osteogenesis.
In this study, we compared gene expression microarray data from 20 osteoporosis patients and 20 healthy subjects, and identified 638 DEGs, including 371 up-regulated genes and 267 down-regulated genes. We constructed a PPI network showing functional interactions between the DEGs and identified 10 hub genes, namely, TP53, RPS27A, VEGFA, MAPK8, CDC42, CREBBP, SIRT1, RPL35A, RPL30, and SNRPG, using CytoScaPe plugin of Cytoscape software.
The top three hub genes are TP53, RPS27A, and VEGFA. Jia et al. reported that pri-miR-34b/c rs4938723 and TP53 Arg72Pro polymorphisms contribute to the risk of osteoporosis . Xie et al. reported that hMSCs from osteoporosis patients show significantly higher TP53 expression; moreover, TP53, SP1 and CTNNB1 transcription factors regulate most of the upregulated DEGs . Felthaus et al. demonstrated that SP1 and TP53 regulate osteogenic differentiation in both dental follicle cells (DFCs) and stem cells from human exfoliated deciduous teeth . The tumor suppressor p53 binds to OSX and prevents its binding to DLX5, thereby repressing osteoblast differentiation via deregulation of the osteogenic transcriptional network . Several studies show that VEGFA plays an important role in osteoporosis [35, 36]. VEGFA is a pro-angiogenic factor that is upregulated in response to uniaxial cyclic tensile strain in human adipose-derived stem cells (hASCs) and hMSCs from osteoporotic donors . In this study, we demonstrate that VEGFA is one of the top hub genes in the PPI network.
Hsa-miR-16-5p is the mature miRNA generated from precursor miRNAs such as miR-16-1, miR-16-2 . Duan et al. reported that osteogenic differentiation is reduced by miR-16-2 upregulation and enhanced by downregulation of miR-16-2 . Castaño et al. showed that antagomiR-16 increased the expression of a key osteogenic transcription factor RunX2 and the levels of OCN, an osteogenic biomarker . Wang et al. showed that Hsa-miR-16 was differentially expressed in the sera of patients with osteogenesis imperfecta, a genetic bone disease . In this study, we used TargetScanHuman and DIANA software and identified hsa-miR-16-5p as a potential miRNA that targets VEGFA mRNA. We confirmed this observation using dual luciferase reporter assays. Furthermore, we showed that miR-16-5p levels were upregulated and VEGFA levels were downregulated in osteoporosis patients. Moreover, alizarin red staining of hMSCs showed that high miR-16-5p expression reduced bone mineralization whereas miR-16-5p knockdown increased bone density.
Our study has a few limitations. First, we did not perform in vivo experiments to determine the role of miR-16-5p/VEGFA in osteoporosis. Secondly, we did not analyze other subtypes of osteoporosis in this study. Therefore, our results are applicable only to postmenopausal osteoporosis, which is the most common clinical subtype of osteoporosis. Further investigations are necessary to identify molecular mechanisms involved in other subtypes of osteoporosis. Moreover, future studies need to investigate the role of other miRNAs that regulate osteoporosis in order to gain a comprehensive understanding of the miRNA-mRNA axis that regulates osteoporosis.
In conclusion our study shows that miR-16-5p suppresses osteogenic differentiation by down-regulating VEGFA expression. Therefore, miR-16-5p/VEGFA axis is a potential therapeutic target for postmenopause-related osteoporosis.
Materials and Methods
Data search and identification of deferentially expressed genes (DEGs)
The RNA microarray data of 20 postmenopausal women with osteoporosis and 20 healthy postmenopausal women in the GSE56815 dataset was retrieved from NCBI Gene Expression Omnibus (GEO) database. The raw gene expression data was log2 transformed and differentially expressed genes were identified using the GEO2R in-built function with a default setting (https://www.ncbi.nlm.nih.gov/geo/geo2r/). DEGs with a P < 0.05 were considered as statistically significant. A heatmap was constructed from the log2 mRNA expression data using the pheatmap R package.
PPI network of the DEGs
We imported the DEGs into the Search Tool for the Retrieval of Interacting Genes (STRING) database  and identified protein-protein interactions with a combined score of >0.5. PPI networks were constructed using the Cytoscape software (version 3.7.2) . Molecular Complex Detection (MCODE) plugin from Cytoscape was used to screen modules of the PPI network, and those with MCODE score ≥ 5 and number of nodes >5 were selected. CentiScaPe 2.2 plugin in Cytoscape was used to calculate the degree, betweenness, and closeness centralities of nodes in the PPI network. Node degree is a measure that represents the number of connections associated with a specific node in the network. Closeness centrality defines how close a node is to all other nodes in the network. Betweenness centrality is the number of times a node acts as bridge along the shortest path between two other nodes.
GO and KEGG pathway enrichment analysis
The function and pathway enrichment analyses of the candidate genes were performed using the DAVID. GO annotation was performed to identify top 3 enriched GO terms associated with biological processes, molecular functions and cellular components. Moreover, the top enriched KEGG pathways involved in osteoporosis were also analyzed using DAVID.
The top ten genes with the highest degrees (hub genes) were visualized using the ggplot2 software . GO plot was used to visualize the results of hub gene enrichment analysis . Circular Visualization in R  was used to visualize chromosomal positions in a circular representation and the degree connectivity of the top 100 genes.
Cell culture and transfection
The hMSCs were obtained from the Huazhong University of Science and Technology, Wuhan, China and grown in a specific media designed for human mesenchymal stem cells (#MUXMA-90011, Cyagen, Suzhou, China) at 37°C in a 5% CO2 incubator. They were maintained for a maximum of 3 passages for experiments. The agomiR-16-5p, antagomiR-16-5p, agomiR-NC, antagomiR-NC, si-NC and si-VEGFA constructs were obtained from GenePharma (Shanghai). Cell transfection experiments were performed using Lipofectamine 3000 (ThermoFisher Scientific) according to manufacturer’s instructions using miRNA constructs at 200 mM and siRNA constructs at 50 mM.
From May 2016 to June 2018, peripheral blood samples from patients in Shanghai Tongji Hospital (10 healthy volunteers, 10 postmenopausal osteoporosis patients) were collected 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.
Total RNA was extracted using Trizol (#15596018, Invitrogen, USA) and reverse transcribed using the VersoTM cDNA Synthesis Kit (#AB-1054/A, ThermoFisher Scientific) according to the manufacturer’s protocol. The exosomal miRNAs were isolated using the SeraMir Exosome RNA purification Kit (System Biosciences, USA), and reverse transcribed and quantified using the TaqMan microRNA assay kit (Applied Biosystems, USA) for cDNA synthesis and qRT-PCR. The qRT-PCR reactions were performed using the Thermal Cycler C-1000 Touch system (#10021377, Bio-Rad CFX Manager, USA). U6 was used as a control for miRNA quantification whereas GAPDH was used as the internal control to quantify mRNA. Osteogenesis was estimated by analyzing the expression of osteogenic marker genes, ALP, RUNX2 and OCN by qRT-PCR analysis. The data was expressed as fold change relative to the appropriate controls. All the primers used in this study are listed in Table 3.
Table 3. QRT-PCR primers used in the study.
|Gene name||Primer sequence|
|hsa - miR-16-5p - Forward||TGGGGTAGCAGCACGTAAA|
|hsa - miR-16-5p - Reverse||CTCAACTGGTGTCGTGGAGTC|
|hsa-ALP - Forward||ACCACCACGAGAGTGAACCA|
|hsa-ALP - Reverse||CGTTGTCTGAGTACCAGTCCC|
|hsa-OCN - Forward||CAAAGGTGCAGCCTTTGTGTC|
|hsa-OCN - Reverse||TCACAGTCCGGATTGAGCTCA|
|hsa-Runx2 - Forward||TGGTTACTGTCATGGCGGGTA|
|hsa-Runx2 - Reverse||TCTCAGATCGTTGAACCTTGCTA|
|hsa-GAPDH - Forward||CCGTTGAATTTGCCGTGA|
|hsa-GAPDH - Reverse||TGATGACCCTTTTGGCTCCC|
Total protein samples were prepared from cells and callus samples using Protein Lysis buffer (#AS1004, Aspen, South Africa) containing 1% protease inhibitor (#AS1008, Aspen). Equal amounts of protein samples were separated on SDS-PAGE and transferred onto nitrocellulose membranes (#IPVH00010, Millipore, USA). Then, the membranes were blocked by incubation with 5% nonfat milk for 1 h followed by overnight incubation at 4°C with primary antibodies against VEGFA (1:500, Sigma, USA, #HPA069116) and GAPDH (1:10,000, Abcam, USA, #ab37168). Then, the blots were incubated with HRP-conjugated secondary antibody (#AS1058, Aspen). The blots were developed using enhanced chemiluminescence detection system and the expression of VEGFA relative to GAPDH was determined for all samples. Each experiment was repeated three times.
Alizarin red staining
The hMSCs were grown in 6-well plates in specific media containing 100 nM dexamethasone, 50 mM ascorbic acid, and 10 mM b-glycerophosphate to promote osteogenesis (#HUXMA-90021, Cyagen, USA). Then, the cells were washed twice with 1X PBS and fixed in 10% formalin for 15 minutes. Subsequently, the cells were stained with 1 mL 0.5% alizarin red staining solution at room temperature for 15 minutes. After rinsing the cells with distilled water for 5 minutes, the cells were mounted on slides and analyzed for red mineralized nodules using the charge-coupled device microscope. The absorbance was measured at 570 nm. The experiments were repeated in triplicate.
Luciferase reporter assays
The hMSCs were grown in 24-well plates (2.5×105 cells/well) and transfected with dual-luciferase vectors containing VEGFA-WT-3’UTR and VEGFA-MUT-3’UTR plus miR-16-5p agonist (agomiR-16-5p) or negative control (agomiR-NC). A Quik Change Site-Directed Mutagenesis Kit (Strata gene) was used to mutate the miR-16-5p binding-region in the VEGFA-3’UTR. The dual luciferase reporter kit (Promega) was used to perform the luciferase assay according to the manufacturer’s instructions. A luminometer (Glomax, Promega) was used to quantify luminescence from the firefly luciferase and control Renilla luciferase constructs in each group. The values of firefly luciferase activity were normalized to the corresponding Renilla signal.
The data are presented as means ± SD. Binary groups were compared using Student's t-test,, whereas, multiple groups (more than two) were compared using one-way ANOVA with Tukey post hoc test. The statistical analyses were conducted using the Graphpad Prism 8.0 software (Graphpad software, San Diego, CA, USA). P<0.05 was considered statistically significant. The schema of the approach process in this study is shown in Figure 8.
YY and YX conceived and designed the study; GY, HZ and BL supervised the study; TY and XY performed the bioinformatics analysis and experiments; HZ, WH and ZL analyzed the data; JX and YZ provided advice and technical assistance; TY and XY wrote the manuscript. All authors approved the final manuscript.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
This study was supported by the National Science Foundation of China (Grant Nos.31600754 and 81472144).
- 1. Burge R, Dawson-Hughes B, Solomon DH, Wong JB, King A, Tosteson A. Incidence and economic burden of osteoporosis-related fractures in the United States, 2005-2025. J Bone Miner Res. 2007; 22:465–75. https://doi.org/10.1359/jbmr.061113 [PubMed]
- 2. Yan G, Huang Y, Cao H, Wu J, Jiang N, Cao X. Association of breastfeeding and postmenopausal osteoporosis in Chinese women: a community-based retrospective study. BMC Womens Health. 2019; 19:110. https://doi.org/10.1186/s12905-019-0808-0 [PubMed]
- 3. Rajasree Vijayakumar DB. Osteoporosis: an under-recognized public health problem. Journal of Local and Global Health Science. 2016; 2:1–13. https://doi.org/10.5339/jlghs.2016.2
- 4. McClung MR, O’Donoghue ML, Papapoulos SE, Bone H, Langdahl B, Saag KG, Reid IR, Kiel DP, Cavallari I, Bonaca MP, Wiviott SD, de Villiers T, Ling X, et al, and LOFT Investigators. Odanacatib for the treatment of postmenopausal osteoporosis: results of the LOFT multicentre, randomised, double-blind, placebo-controlled trial and LOFT Extension study. Lancet Diabetes Endocrinol. 2019; 7:899–911. https://doi.org/10.1016/S2213-8587(19)30346-8 [PubMed]
- 5. Ma M, Chen X, Lu L, Yuan F, Zeng W, Luo S, Yin F, Cai J. Identification of crucial genes related to postmenopausal osteoporosis using gene expression profiling. Aging Clin Exp Res. 2016; 28:1067–74. https://doi.org/10.1007/s40520-015-0509-y [PubMed]
- 6. Cai X, Yi X, Zhang Y, Zhang D, Zhi L, Liu H. Genetic susceptibility of postmenopausal osteoporosis on sulfide quinone reductase-like gene. Osteoporos Int. 2018; 29:2041–47. https://doi.org/10.1007/s00198-018-4575-9 [PubMed]
- 7. Qian GF, Yuan LS, Chen M, Ye D, Chen GP, Zhang Z, Li CJ, Vijayan V, Xiao Y. PPWD1 is associated with the occurrence of postmenopausal osteoporosis as determined by weighted gene co-expression network analysis. Mol Med Rep. 2019; 20:3202–14. https://doi.org/10.3892/mmr.2019.10570 [PubMed]
- 8. Hausser J, Zavolan M. Identification and consequences of miRNA-target interactions—beyond repression of gene expression. Nat Rev Genet. 2014; 15:599–612. https://doi.org/10.1038/nrg3765 [PubMed]
- 9. Henshall DC, Hamer HM, Pasterkamp RJ, Goldstein DB, Kjems J, Prehn JH, Schorge S, Lamottke K, Rosenow F. MicroRNAs in epilepsy: pathophysiology and clinical utility. Lancet Neurol. 2016; 15:1368–76. https://doi.org/10.1016/S1474-4422(16)30246-0 [PubMed]
- 10. Alberti C, Manzenreither RA, Sowemimo I, Burkard TR, Wang J, Mahofsky K, Ameres SL, Cochella L. Cell-type specific sequencing of microRNAs from complex animal tissues. Nat Methods. 2018; 15:283–89. https://doi.org/10.1038/nmeth.4610 [PubMed]
- 11. Li H, Wang Z, Fu Q, Zhang J. Plasma miRNA levels correlate with sensitivity to bone mineral density in postmenopausal osteoporosis patients. Biomarkers. 2014; 19:553–56. https://doi.org/10.3109/1354750X.2014.935957 [PubMed]
- 12. Kocijan R, Muschitz C, Geiger E, Skalicky S, Baierl A, Dormann R, Plachel F, Feichtinger X, Heimel P, Fahrleitner-Pammer A, Grillari J, Redl H, Resch H, Hackl M. Circulating microRNA Signatures in Patients With Idiopathic and Postmenopausal Osteoporosis and Fragility Fractures. J Clin Endocrinol Metab. 2016; 101:4125–34. https://doi.org/10.1210/jc.2016-2365 [PubMed]
- 13. Cao Z, Moore BT, Wang Y, Peng XH, Lappe JM, Recker RR, Xiao P. MiR-422a as a potential cellular microRNA biomarker for postmenopausal osteoporosis. PLoS One. 2014; 9:e97098. https://doi.org/10.1371/journal.pone.0097098 [PubMed]
- 14. Li Z, Zhang W, Huang Y. MiRNA-133a is involved in the regulation of postmenopausal osteoporosis through promoting osteoclast differentiation. Acta Biochim Biophys Sin (Shanghai). 2018; 50:273–80. https://doi.org/10.1093/abbs/gmy006 [PubMed]
- 15. Wang WW, Yang L, Wu J, Gao C, Zhu YX, Zhang D, Zhang HX. The function of miR-218 and miR-618 in postmenopausal osteoporosis. Eur Rev Med Pharmacol Sci. 2017; 21:5534–41. https://doi.org/10.26355/eurrev_201712_13989 [PubMed]
- 16. You L, Pan L, Chen L, Gu W, Chen J. MiR-27a is Essential for the Shift from Osteogenic Differentiation to Adipogenic Differentiation of Mesenchymal Stem Cells in Postmenopausal Osteoporosis. Cell Physiol Biochem. 2016; 39:253–65. https://doi.org/10.1159/000445621 [PubMed]
- 17. Qiu J, Huang G, Na N, Chen L. MicroRNA-214-5p/TGF-β/Smad2 signaling alters adipogenic differentiation of bone marrow stem cells in postmenopausal osteoporosis. Mol Med Rep. 2018; 17:6301–10. https://doi.org/10.3892/mmr.2018.8713 [PubMed]
- 18. Qiao L, Liu D, Li CG, Wang YJ. MiR-203 is essential for the shift from osteogenic differentiation to adipogenic differentiation of mesenchymal stem cells in postmenopausal osteoporosis. Eur Rev Med Pharmacol Sci. 2018; 22:5804–14. https://doi.org/10.26355/eurrev_201809_15906 [PubMed]
- 19. Mangul S, Martin LS, Langmead B, Sanchez-Galan JE, Toma I, Hormozdiari F, Pevzner P, Eskin E. How bioinformatics and open data can boost basic science in countries and universities with limited resources. Nat Biotechnol. 2019; 37:324–26. https://doi.org/10.1038/s41587-019-0053-y [PubMed]
- 20. Carbonell P, Currin A, Jervis AJ, Rattray NJ, Swainston N, Yan C, Takano E, Breitling R. Bioinformatics for the synthetic biology of natural products: integrating across the Design-Build-Test cycle. Nat Prod Rep. 2016; 33:925–32. https://doi.org/10.1039/C6NP00018E [PubMed]
- 21. Teufel A. Bioinformatics and database resources in hepatology. J Hepatol. 2015; 62:712–19. https://doi.org/10.1016/j.jhep.2014.10.036 [PubMed]
- 22. Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. eLife. 2015; 4:4. https://doi.org/10.7554/eLife.05005 [PubMed]
- 23. Paraskevopoulou MD, Vlachos IS, Karagkouni D, Georgakilas G, Kanellos I, Vergoulis T, Zagganas K, Tsanakas P, Floros E, Dalamagas T, Hatzigeorgiou AG. DIANA-LncBase v2: indexing microRNA targets on non-coding transcripts. Nucleic Acids Res. 2016; 44:D231–38. https://doi.org/10.1093/nar/gkv1270 [PubMed]
- 24. Liu Y, Berendsen AD, Jia S, Lotinun S, Baron R, Ferrara N, Olsen BR. Intracellular VEGF regulates the balance between osteoblast and adipocyte differentiation. J Clin Invest. 2012; 122:3101–13. https://doi.org/10.1172/JCI61209 [PubMed]
- 25. Behr B, Tang C, Germann G, Longaker MT, Quarto N. Locally applied vascular endothelial growth factor A increases the osteogenic healing capacity of human adipose-derived stem cells by promoting osteogenic and endothelial differentiation. Stem Cells. 2011; 29:286–96. https://doi.org/10.1002/stem.581 [PubMed]
- 26. Hu K, Olsen BR. Osteoblast-derived VEGF regulates osteoblast differentiation and bone formation during bone repair. J Clin Invest. 2016; 126:509–26. https://doi.org/10.1172/JCI82585 [PubMed]
- 27. Li GQ, Wang ZY. MiR-20b promotes osteocyte apoptosis in rats with steroid-induced necrosis of the femoral head through BMP signaling pathway. Eur Rev Med Pharmacol Sci. 2019; 23:4599–608. https://doi.org/10.26355/eurrev_201906_18037 [PubMed]
- 28. Sang S, Zhang Z, Qin S, Li C, Dong Y. MicroRNA-16-5p Inhibits Osteoclastogenesis in Giant Cell Tumor of Bone. BioMed Res Int. 2017; 2017:3173547. https://doi.org/10.1155/2017/3173547 [PubMed]
- 29. Hui S, Yang Y, Li J, Li N, Xu P, Li H, Zhang Y, Wang S, Lin G, Li S, Qiu G, Zhao RC, Zhang J, Zhuang Q. Differential miRNAs profile and bioinformatics analyses in bone marrow mesenchymal stem cells from adolescent idiopathic scoliosis patients. Spine J. 2019; 19:1584–96. https://doi.org/10.1016/j.spinee.2019.05.003 [PubMed]
- 30. Francis M, Grider A. MiRNA-target interactions in osteogenic signaling pathways involving zinc via the metal regulatory element. Biometals. 2019; 32:111–21. https://doi.org/10.1007/s10534-018-00162-4 [PubMed]
- 31. Jia F, Sun R, Li J, Li Q, Chen G, Fu W. Interactions of Pri-miRNA-34b/c and TP53 Polymorphisms on the Risk of Osteoporosis. Genet Test Mol Biomarkers. 2016; 20:398–401. https://doi.org/10.1089/gtmb.2015.0282 [PubMed]
- 32. Xie W, Ji L, Zhao T, Gao P. Identification of transcriptional factors and key genes in primary osteoporosis by DNA microarray. Med Sci Monit. 2015; 21:1333–44. https://doi.org/10.12659/MSM.894111 [PubMed]
- 33. Felthaus O, Viale-Bouroncle S, Driemel O, Reichert TE, Schmalz G, Morsczeck C, Morsczeck C. Transcription factors TP53 and SP1 and the osteogenic differentiation of dental stem cells. Differentiation. 2012; 83:10–16. https://doi.org/10.1016/j.diff.2011.08.008 [PubMed]
- 34. Artigas N, Gámez B, Cubillos-Rojas M, Sánchez-de Diego C, Valer JA, Pons G, Rosa JL, Ventura F. p53 inhibits SP7/Osterix activity in the transcriptional program of osteoblast differentiation. Cell Death Differ. 2017; 24:2022–31. https://doi.org/10.1038/cdd.2017.113 [PubMed]
- 35. Toti P, Sbordone C, Martuscelli R, Califano L, Ramaglia L, Sbordone L. Gene clustering analysis in human osteoporosis disease and modifications of the jawbone. Arch Oral Biol. 2013; 58:912–29. https://doi.org/10.1016/j.archoralbio.2013.02.013 [PubMed]
- 36. Lee KH, Kim SH, Kim CH, Min BJ, Kim GJ, Lim Y, Kim HS, Ahn KM, Kim JH. Identifying genetic variants underlying medication-induced osteonecrosis of the jaw in cancer and osteoporosis: a case control study. J Transl Med. 2019; 17:381. https://doi.org/10.1186/s12967-019-2129-3 [PubMed]
- 37. Charoenpanich A, Wall ME, Tucker CJ, Andrews DM, Lalush DS, Dirschl DR, Loboa EG. Cyclic tensile strain enhances osteogenesis and angiogenesis in mesenchymal stem cells from osteoporotic donors. Tissue Eng Part A. 2014; 20:67–78. https://doi.org/10.1089/ten.tea.2013.0006 [PubMed]
- 38. Cai B, Ma M, Chen B, Li Z, Abdalla BA, Nie Q, Zhang X. MiR-16-5p targets SESN1 to regulate the p53 signaling pathway, affecting myoblast proliferation and apoptosis, and is involved in myoblast differentiation. Cell Death Dis. 2018; 9:367. https://doi.org/10.1038/s41419-018-0403-6 [PubMed]
- 39. Duan L, Zhao H, Xiong Y, Tang X, Yang Y, Hu Z, Li C, Chen S, Yu X. miR-16-2* Interferes with WNT5A to Regulate Osteogenesis of Mesenchymal Stem Cells. Cell Physiol Biochem. 2018; 51:1087–102. https://doi.org/10.1159/000495489 [PubMed]
- 40. Mencía Castaño I, Curtin CM, Duffy GP, O’Brien FJ. Harnessing an Inhibitory Role of miR-16 in Osteogenesis by Human Mesenchymal Stem Cells for Advanced Scaffold-Based Bone Tissue Engineering. Tissue Eng Part A. 2019; 25:24–33. https://doi.org/10.1089/ten.tea.2017.0460 [PubMed]
- 41. Wang Z, Lu Y, Zhang X, Ren X, Wang Y, Li Z, Xu C, Han J. Serum microRNA is a promising biomarker for osteogenesis imperfecta. Intractable Rare Dis Res. 2012; 1:81–85. https://doi.org/10.5582/irdr.2012.v1.2.81 [PubMed]
- 42. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P, Jensen LJ, Mering CV. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019; 47:D607–D613. https://doi.org/10.1093/nar/gky1131 [PubMed]
- 43. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
- 44. Wickham H. ggplot2. Elegant Graphics for Data Analysis. Springer-Verlag, New York, 2016. https://doi.org/10.1007/978-3-319-24277-4_9
- 45. Walter W, Sánchez-Cabo F, Ricote M. GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics. 2015; 31:2912–14. https://doi.org/10.1093/bioinformatics/btv300 [PubMed]
- 46. Gu Z, Gu L, Eils R, Schlesner M, Brors B. circlize Implements and enhances circular visualization in R. Bioinformatics. 2014; 30:2811–12. https://doi.org/10.1093/bioinformatics/btu393 [PubMed]