MiR-16-5p regulates postmenopausal osteoporosis by directly targeting VEGFA

Results

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
SIRT1	Sirtuin 1	33
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.

Figure 3. The degree, betweenness, and closeness centrality of hub genes. The analysis shows that TP53, RPS27A, and VEGFA were the top three hub genes with greatest degree, betweenness and closeness.

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
Term	Name	Count	P-value	Genes
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
Term	Name	Count	P-value	Genes
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
Term	Name	Count	P-value	Genes
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
Term	Name	Count	P-value	Genes
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 [22] and DIANA [23] 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 [24] 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.

Figure 5. The miRNA-mRNA targeting prediction network. The results of TargetScanHuman and DIANA software analyses show the miRNAs that potentially target the mRNAs of the top 3 hub genes, TP53, RPS27A, and VEGFA.

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.

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 [42] and identified protein-protein interactions with a combined score of >0.5. PPI networks were constructed using the Cytoscape software (version 3.7.2) [43]. 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.

Circular visualization

The top ten genes with the highest degrees (hub genes) were visualized using the ggplot2 software [44]. GO plot was used to visualize the results of hub gene enrichment analysis [45]. Circular Visualization in R [46] 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% CO₂ 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.

Blood collection

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.

QRT-PCR analysis

Total RNA was extracted using Trizol (#15596018, Invitrogen, USA) and reverse transcribed using the Verso^TM 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-U6-Forward	CTCGCTTCGGCAGCACA
hsa-U6-Reverse	AACGCTTCACGAATTTGCGT
hsa-VEGFA-Forward	AGGGCAGAATCATCACGAAGT
hsa-VEGFA-Reverse	AGGGTCTCGATTGGATGGCA
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

Western blotting

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×10⁵ 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.

Statistical analysis

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.

Figure 8. Schematic representation of the experimental strategy in the present study.

Author Contributions

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.

Funding

This study was supported by the National Science Foundation of China (Grant Nos.31600754 and 81472144).

References

• 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]