Research | Open | Published:
Integrated network analysis to explore the key genes regulated by parathyroid hormone receptor 1 in osteosarcoma
World Journal of Surgical Oncologyvolume 15, Article number: 177 (2017)
As an invasive malignant tumor, osteosarcoma (OS) has high mortality. Parathyroid hormone receptor 1 (PTHR1) contributes to maintaining proliferation and undifferentiated state of OS. This study is designed to reveal the action mechanisms of PTHR1 in OS.
Microarray dataset GSE46861, which included six PTHR1 knockdown OS samples and six control OS samples, was obtained from the Gene Expression Omnibus database. The differentially expressed genes (DEGs) were identified and then performed with enrichment analysis separately using the limma package and DAVID online tool. Then, protein-protein interaction (PPI) network and module analyses were conducted using Cytoscape software. Using the WebGestalt tool, microRNAs (miRNAs) were predicted for the DEGs involved in the PPI network. Following this, transcription factors (TFs) were predicted and an integrated network was constructed by Cytoscape software.
There were 871 DEGs in the PTHR1 knockdown OS samples compared with the control OS samples. Besides, upregulated ZFPM2 was involved in the miRNA-DEG regulatory network. Moreover, TF LEF1 was predicted for the miRNA-DEG regulatory network of the downregulated genes. In addition, LEF1, NR4A2, HAS2, and RHOC had higher degrees in the integrated network.
ZFPM2, LEF1, NR4A2, HAS2, and RHOC might be potential targets of PTHR1 in OS.
As an invasive malignant tumor, osteosarcoma (OS) often occurs in tubular long bones . OS is the most common type of primary bone cancer, which has high occurrence rate in children and teenagers . OS has rapid growth, high metastatic potential, and local aggressiveness; thus, it can result in high mortality . In childhood cancers, OS accounts for about 2.4% of all malignant cancers and is the eighth most common tumor . Therefore, revealing the molecular mechanisms of OS and developing novel therapies are of great importance.
Via upregulating matrix metalloproteinase 2 (MMP2), astrocyte elevated gene 1 (AEG1) functions in OS progression can be used for predicting the progression and prognosis of the disease [5, 6]. A previous study reports that miR-203 is a tumor suppressor by mediating RAB22A, member RAS oncogene family (RAB22A) expression, and has correlation with the progression and carcinogenesis of OS . Through nuclear factor kappa B (NF-κB) signaling pathway and mitochondria pathway, the inhibitor of growth 4 (ING4) plays a suppressive role in OS progression and serves as a potential target for treating OS [8, 9]. miR-24 inhibits the metastasis of OS via regulating activated Cdc42-associated kinase (ACK1) through AKT/MMP pathways, which may be applied for the diagnosis and therapy of OS . However, the pathogenesis of OS has not been comprehensively revealed.
Parathyroid hormone receptor 1 (PTHR1) signaling plays a critical role in keeping proliferation and undifferentiated state of OS; thus, PTHR1 suppression may be used to inhibit OS proliferation and promote differentiation [11, 12]. Overexpression of PTHR1 may contribute to OS progression through affecting aggressive phenotype and microenvironment . To explore the action mechanisms of PTHR1 in OS, we downloaded the microarray dataset GSE46861 which included both PTHR1 knockdown OS samples and control OS samples. Then, differentially expressed genes (DEGs) were identified and performed with enrichment analysis. In addition, protein-protein interaction (PPI) network and module analyses, as well as integrated network analysis, were conducted to further screen the key targets of PTHR1 in OS.
Data source and data preprocessing
Microarray dataset GSE46861, which was sequenced on the platform of GPL6246 [MoGene-1_0-st] Affymetrix Mouse Gene 1.0 ST Array [transcript (gene) version], was obtained from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) database. GSE46861 included six PTHR1 knockdown OS samples and six control OS samples. The raw data of GSE46861 was normalized using the Robust Multiarray Average (RMA) method  in the R package Affy (version: 1.52.0).
DEG screening and hierarchical cluster analysis
After the data normalization, the DEGs between PTHR1 knockdown OS samples and control OS samples were screened using the Linear Models for Microarray Analysis (limma, version: 3.30.3)  package in R. The |log2 fold change (FC)| > 0.58 and p value < 0.05 were considered as the thresholds. In addition, pheatmap package (version 1.0.8, https://cran.r-project.org/web/packages/pheatmap/index.html) in R was utilized to perform the hierarchical cluster analysis.
Functional and pathway enrichment analysis
Gene Ontology (GO) database can annotate genes and gene products from molecular function (MF), biological process (BP), and cellular component (CCo) aspects . The Kyoto Encyclopedia of Genes and Genomes (KEGG) database can be applied for revealing gene functions and connecting genomic information with functional information . Using the Database for Annotation, Visualization and Integrated Discovery (DAVID, version: 6.8, parameter: Classification Stringency was set as Medium) online tool , GO functional and KEGG pathway enrichment analyses were carried out for the DEGs, with the threshold of p value < 0.05.
PPI network and module analyses
Search Tool for the Retrieval of Interacting Genes (STRING, version: 10.0) database , which included PPI pairs of multiple organisms, was used to analyze the PPI pairs among the DEGs. The required confidence (combined score) > 0.4 was used as the cutoff criterion. Afterwards, PPI network was visualized using Cytoscape software . Based on the CytoNCA plugin (version 2.1.6, parameter: without weight)  in Cytoscape software, topological property analysis was performed for the nodes in the PPI network. According to closeness centrality (CCe), betweenness centrality (BC), and degree centrality (DC) scores, the hub nodes  in the PPI network were selected. Moreover, the MCODE plugin (version 1.2)  in Cytoscape software was used for identifying the significant modules in the PPI network. In addition, enrichment analysis for the genes involved in the most significant module was conducted using the DAVID online tool .
Integrated network analysis
Using WEB-based gene set analysis toolkit (WebGestalt) tool , miRNAs were predicted for the DEGs involved in the PPI network, with false discovery rate (FDR, that was adjusted p value) < 0.05 and number of target genes ≥ 5 as the thresholds. Based on the iRegulon plugin (version 1.3)  in Cytoscape software, transcription factors (TFs) were further predicted for the miRNA-DEG regulatory network. The Minimum NEScore > 3 and FDR on motif similarity < 0.001 were set as thresholds. Finally, the obtained transcription regulation relationships were merged into the miRNA-DEG regulatory network, and the integrated network was visualized by Cytoscape software .
DEG analysis and hierarchical cluster analysis
Compared with the control OS samples, there were a total of 871 DEGs (438 upregulated and 433 downregulated genes) in the PTHR1 knockdown OS samples. The heatmap of the DEGs is shown in Fig. 1, which indicated that the expression of identified DEGs could correctly distinguish the two kinds of samples.
Functional and pathway enrichment analysis
Enrichment analysis was conducted for the upregulated genes and the downregulated genes, respectively. A total of 121 GO_BP terms and 20 pathways were enriched for the upregulated genes. The top five GO_BP terms and pathways are shown in Fig. 2a, mainly including immune system process (GO_BP, p value = 4.07E−07) and Staphylococcus aureus infection (pathway, p value = 1.09E−08). Besides, a total of 65 GO_BP terms and five pathways were enriched for the downregulated genes. Similarly, the top five GO_BP terms (such as sterol biosynthetic process, p value = 1.02E−06) and pathways (such as Biosynthesis of antibiotics, p value = 6.96E−05) are shown in Fig. 2b.
PPI network and module analyses
The PPI network for the upregulated genes had 280 nodes and 1090 interactions. According to DC scores, the nodes with degrees larger than 30 are listed in Table 1. The most significant module (score = 25.786) identified from the PPI network for the upregulated genes had 29 nodes and 361 interactions (Fig. 3a). The upregulated genes involved in the most significant module were mainly enriched in immune system process (GO_BP, p value = 1.63E−07) and Staphylococcus aureus infection (pathway, p value = 2.59E−10) (Table 2 (A)).
Besides, the PPI network constructed for the downregulated genes included 262 nodes and 503 interactions. In the PPI network, discs, large homolog 4 (DLG4) was the only node with degree larger than 20. Module analysis for the PPI network of the downregulated genes showed that the most significant module (score = 10.364) had 12 nodes and 57 interactions (Fig. 3b). The downregulated genes involved in the most significant module were mainly enriched in sterol metabolic process (GO_BP, p value = 4.60E–14) and steroid biosynthesis (pathway, p value = 1.10E−08) (Table 2 (B)).
Integrated network analysis
The miRNA-DEG regulatory network constructed for the upregulated genes (such as zinc finger protein, multitype 2, ZFPM2) and the downregulated genes are separately shown in Figs. 4 and 5. There was no TF predicted for the miRNA-DEG regulatory network of the upregulated genes. Only two TFs (SRY (sex determining region Y)-box 12, SOX12, and lymphoid enhancer binding factor 1, LEF1) were predicted for the miRNA-DEG regulatory network of the downregulated genes. The integrated network for the downregulated genes had 278 nodes and 1144 edges (Fig. 6). Especially, LEF1, nuclear receptor subfamily 4 group A member 2 (NR4A2), hyaluronan synthase 2 (HAS2), and ras homolog family member C (RHOC) had higher degrees in the integrated network.
In this study, a total of 871 DEGs were identified in the PTHR1 knockdown OS samples compared with the control OS samples, including 438 upregulated and 433 downregulated genes. ZFPM2 was involved in the miRNA-DEG regulatory network constructed for the upregulated genes. There was no TF predicted for the miRNA-DEG regulatory network of the upregulated genes. Besides, TF LEF1 was predicted for the miRNA-DEG regulatory network of the downregulated genes. In the integrated network, LEF1, NR4A2, HAS2, and RHOC had higher degrees.
In rat bone marrow stromal cells (BMSCs), neurotrophic factors, NURR1 (also named NR4A2), tyrosine hydroxylase (TH), and nestin genes have spontaneous expression . NURR1 and NURR77 may contribute to increase the migratory potential of fetal FBMSCs, which may mediate the local immune response specifically . NURR1 and PPARγ coactivator-1α (PGC-1α) may play pivotal roles in mediating osteoblast function and cAMP-dependent osteoblast gene expression [28, 29]. NURR1 maintains cartilage homeostasis via selectively inhibiting the expression of MMP gene during inflammation . These declared that PTHR1 might function in OS through targeting NR4A2.
Hyaluronan synthesized by HAS2 affects the proliferation, invasion, and motility of MG-63 OS cells . Through regulating the expressions of versican, HAS2, and hyaluronan, transforming growth factor β2 (TGF-β2) may lead to the metastasis of OS cells . HAS2 plays critical roles in osteoblast differentiation and development by mediating high molecular weight hyaluronan synthesis . RHOC and MMP9 expression levels have close association, and their high expressions are closely correlated with the formation, development, invasion, and metastasis of OS . RHOC has different expressions in SOSP-9607E10 and SOSP-9607H9 OS cell lines, and its overexpression functions in the invasion and metastasis of OS through inducing cell migration . Therefore, HAS2 and RHOC might also be targets of PTHR1 in OS.
Via regulating the expression of ZFPM2, hypoxia-induced miR-429 contributes to the differentiation of osteoblastic cells . LEF1 can delay osteoblast maturation and regulates the expression levels of some genes in osteoblasts . LEF1 has an essential role in the activation of alpha 1 (XI) collagen (COL11A1), and COL11A1 inhibits the terminal differentiation of osteoblasts . LEF1 serves as a transcriptional effector of the Wnt/β-catenin pathway and is critical for the tumor invasion induced by hepatocyte growth factor (HGF) . LEF1 mediates bone density, osteoblast differentiation, and skeletal strength, and Lef1ΔN regulates the terminal differentiation in osseous cells [40, 41]. Thus, ZFPM2 and LEF1 might be targets of PTHR1 in OS.
In conclusion, a total of 871 DEGs were screened from the PTHR1 knockdown OS samples. Besides, ZFPM2, LEF1, NR4A2, HAS2, and RHOC might be targets of PTHR1 in OS. However, more experimental researches should be conducted to confirm these findings obtained from bioinformatics analysis.
- ACK1 :
Activated Cdc42-associated kinase
- AEG1 :
Elevated gene 1
Bone marrow stromal cells
Database for Annotation, Visualization and Integrated Discovery
Differentially expressed genes
Discs, large homolog 4
False discovery rate
Gene Expression Omnibus
- HAS2 :
Hyaluronan synthase 2
- HGF :
Hepatocyte growth factor
- ING4 :
Inhibitor of growth 4
Kyoto Encyclopedia of Genes and Genomes
- MMP2 :
Matrix metalloproteinase 2
Nuclear factor kappa B
- NR4A2 :
Nuclear receptor subfamily 4 group A member 2
- PGC-1α :
Parathyroid hormone receptor 1
- RHOC :
Ras homolog family member C
Robust Multiarray Average
Search Tool for the Retrieval of Interacting Genes
- TGF-β2 :
Growth factor β2
Neurotrophic factors, tyrosine hydroxylase.
Luetke A, Meyers PA, Lewis I, Juergens H. Osteosarcoma treatment—where do we stand? A state of the art review. Cancer Treat Rev. 2014;40:523–32.
Ottaviani G, Jaffe N, Eftekhari F, Raymond AK, Yasko AW. Pediatric and adolescent osteosarcoma. Berlin: Springer; 2011.
Mohseny AB, Tieken C, van der Velden PA, Szuhai K, de Andrea C, Hogendoorn PC, Cleton-Jansen AM. Small deletions but not methylation underlie CDKN2A/p16 loss of expression in conventional osteosarcoma. Genes Chromosomes Cancer. 2010;49:1095–103.
Ottaviani G, Jaffe N. The epidemiology of osteosarcoma. Cancer Treat Res. 2009;152:3–13.
Wang F, Ke ZF, Sun SJ, Chen WF, Yang SC, Li SH, Mao XP, Wang L. Oncogenic roles of astrocyte elevated gene-1 (AEG-1) in osteosarcoma progression and prognosis. Cancer Biol Ther. 2011;12:539–48.
Wang F, Ke ZF, Wang R, Wang YF, Huang LL, Wang LT. Astrocyte elevated gene-1 (AEG-1) promotes osteosarcoma cell invasion through the JNK/c-Jun/MMP-2 pathway. Biochem Biophys Res Commun. 2014;452:933–9.
Yang D, Liu G, Wang K. miR-203 acts as a tumor suppressor gene in osteosarcoma by regulating RAB22A. PLoS One. 2015;10:e0132225.
Mei L, Ye Z, Zhang H, Li L, Peng H, Hong X, Yu Z, Mao C. Delivery of inhibitor of growth 4 (ING4) gene significantly inhibits proliferation and invasion and promotes apoptosis of human osteosarcoma cells. Sci Rep. 2014;4:7380.
Xu M, Xie Y, Sheng W, Miao J, Yang J. Adenovirus-mediated ING4 gene transfer in osteosarcoma suppresses tumor growth via induction of apoptosis and inhibition of tumor angiogenesis. Technol Cancer Res Treat. 2015;14:617–26.
Liu Z, Zhang Y, Li Y, Liu B, Zhang K. miR-24 represses metastasis of human osteosarcoma cells by targeting Ack1 via AKT/MMPs pathway. Biochem Biophys Res Commun. 2017;486:211–7.
Ho PW, Goradia A, Russell MR, Chalk AM, Milley KM, Baker EK, Danks JA, Slavin JL, Walia M, Crimeenirwin B. Knockdown of PTHR1 in osteosarcoma cells decreases invasion and growth and increases tumor differentiation in vivo. Oncogene. 2014;34:2922–33.
Ho PW, Russell M, Goradia A, Chalk A, Slavin J, Dickins R, Martin TJ, Walkley C. Abstract A54: PTHR1 signaling regulates the invasion and differentiation stage of osteosarcoma. Cancer Res. 2014;74:A54.
Yang R, Hoang BH, Kubo T, Kawano H, Chou A, Sowers R, Huvos AG, Meyers PA, Healey JH, Gorlick R. Over-expression of parathyroid hormone Type 1 receptor confers an aggressive phenotype in osteosarcoma. Int J Cancer. 2007;121:943–54.
Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003;4:249–64.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47–e47.
Consortium TGO. Gene Ontology Consortium: going forward. Nucleic Acids Res. 2015;43:1049–56.
Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2015;44:D457–D62.
Huang DW. DAVID Bioinformatics Resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. 2007;35:W169.
Franceschini A, Szklarczyk D, Frankild S, Kuhn M, Simonovic M, Roth A, Lin J, Minguez P, Bork P, von Mering C. STRING v9. 1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 2013;41:D808–D15.
Saito R, Smoot ME, Ono K, Ruscheinski J, Wang P-L, Lotia S, Pico AR, Bader GD, Ideker T. A travel guide to Cytoscape plugins. Nat Methods. 2012;9:1069–76.
Tang Y, Li M, Wang J, Pan Y, Wu FX. CytoNCA: a cytoscape plugin for centrality analysis and evaluation of protein interaction networks. Bio Systems. 2015;127:67–72.
He X, Zhang J. Why do hubs tend to be essential in protein networks? PLoS Genet. 2006;2:826–34.
Bader GD, Hogue CW. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. 2003;4:1–27.
Wang J, Duncan D, Shi Z, Zhang B. WEB-based GEne SeT AnaLysis Toolkit (WebGestalt): update 2013. Nucleic Acids Res. 2013;41:77–83.
Janky RS, Verfaillie A, Imrichová H, Sande BVD, Standaert L, Christiaens V, Hulselmans G, Herten K, Sanchez MN, Potier D. iRegulon: from a gene list to a gene regulatory network using large motif and track collections. PLoS Comput Biol. 2014;10:e1003731.
Fatemeh M, Maryam HGK, Taghi GM, Taghi L. Spontaneous expression of neurotrophic factors and TH, Nurr1, Nestin genes in long-term culture of bone marrow mesenchymal stem cells. Cell J. 2012;13:243.
Maijenburg MW, Gilissen C, Melief SM, Kleijer M, Weijer K, Ten BA, Roelofs H, Van Tiel CM, Veltman JA, de Vries CJ. Nuclear receptors Nur77 and Nurr1 modulate mesenchymal stromal cell migration. Stem Cells Dev. 2012;21:228–38.
Nervina JM, Magyar CE, Pirih FQ, Tetradis S. PGC-1alpha is induced by parathyroid hormone and coactivates Nurr1-mediated promoter activity in osteoblasts. Bone. 2006;39:1018–25.
Siyahian AJ. PTH and Nurr1 mediated gene expression in osteoblasts. Dissertations & Theses - Gradworks. Los Angeles: University of California; 2009.
Mix KS, Attur MG, Al-Mussawir H, Abramson SB, Brinckerhoff CE, Murphy EP. Transcriptional repression of matrix metalloproteinase gene expression by the orphan nuclear receptor NURR1 in cartilage. J Biol Chem. 2007;282:9492–504.
Nishida Y, Knudson W, Knudson CB, Ishiguro N. Antisense inhibition of hyaluronan synthase-2 in human osteosarcoma cells inhibits hyaluronan retention and tumorigenicity. Exp Cell Res. 2005;307:194–203.
Nikitovic D, Zafiropoulos A, Katonis P, Tsatsakis A, Theocharis AD, Karamanos NK, Tzanakakis GN. Transforming growth factor-beta as a key molecule triggering the expression of versican isoforms v0 and v1, hyaluronan synthase-2 and synthesis of hyaluronan in malignant osteosarcoma cells. IUBMB Life. 2006;58:47–53.
Falconi D, Aubin JE. LIF inhibits osteoblast differentiation at least in part by regulation of HAS2 and its product hyaluronan. J Bone Miner Res. 2007;22:1289–300.
Du HZ, Liu XM, Zhang YX. Expression of RhoC and MMP-9 in osteosarcoma and its clinicopathological significance. Mod Prev Med. 2010;37:1794–7.
Tang CY, Jiang-Dong NI. Molecular mechanisms of RhoC inducing invasion and metastasis of osteosarcoma. J Clin Res. 2007;24:1639–43.
Huang J, Peng J, Cao G, Lu S, Liu L, Li Z, Zhou M, Feng M, Shen H. Hypoxia-induced MicroRNA-429 promotes differentiation of MC3T3-E1 osteoblastic cells by mediating ZFPM2 expression. Cell Physiol Biochem. 2016;39:1177–86.
Kahler RA, Galindo M, Lian J, Stein GS, van Wijnen AJ, Westendorf JJ. Lymphocyte enhancer-binding factor 1 (Lef1) inhibits terminal differentiation of osteoblasts. J Cell Biochem. 2006;97:969–83.
Kahler RA, Yingst SM, Hoeppner LH, Jensen ED, Krawczak D, Oxford JT, Westendorf JJ. Collagen 11a1 is indirectly activated by lymphocyte enhancer-binding factor 1 (Lef1) and negatively regulates osteoblast maturation. Matrix Biol. 2008;27:330–8.
Huang FI, Chen YL, Chang CN, Yuan RH, Jeng YM. Hepatocyte growth factor activates Wnt pathway by transcriptional activation of LEF1 to facilitate tumor invasion. Carcinogenesis. 2012;33:1142–8.
Hoeppner LH, Secreto F, Jensen ED, Li X, Kahler RA, Westendorf JJ. Runx2 and bone morphogenic protein 2 regulate the expression of an alternative Lef1 transcript during osteoblast maturation. J Cell Physiol. 2009;221:480–9.
Galindo M, Kahler RA, Teplyuk NM, Stein JL, Lian JB, Stein GS, Westendorf JJ, van Wijnen AJ. Cell cycle related modulations in Runx2 protein levels are independent of lymphocyte enhancer-binding factor 1 (Lef1) in proliferating osteoblasts. J Mol Histol. 2007;38:501–6.
Availability of data and materials
All data generated or analyzed during this study are included in this published article.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Donghui Guan and Honglai Tian are co-first authors.
About this article
- Parathyroid hormone receptor 1
- Differentially expressed genes
- Protein-protein interaction network
- Integrated network