Integrated network analysis to explore the key genes regulated by parathyroid hormone receptor 1 in osteosarcoma

Background 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. Methods 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. Results 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. Conclusions ZFPM2, LEF1, NR4A2, HAS2, and RHOC might be potential targets of PTHR1 in OS.


Background
As an invasive malignant tumor, osteosarcoma (OS) often occurs in tubular long bones [1]. OS is the most common type of primary bone cancer, which has high occurrence rate in children and teenagers [2]. OS has rapid growth, high metastatic potential, and local aggressiveness; thus, it can result in high mortality [3]. In childhood cancers, OS accounts for about 2.4% of all malignant cancers and is the eighth most common tumor [4]. 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 [7]. 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 [10].
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 [13]. 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. 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 [14] 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) [15] package in R. The |log 2 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 [16]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database can be applied for revealing gene functions and connecting genomic information with functional information [17]. Using the Database for Annotation, Visualization and Integrated Discovery (DAVID, version: Fig. 1 The heatmap of the differentially expressed genes by using pheatmap package 6.8, parameter: Classification Stringency was set as Medium) online tool [18], 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 [19], 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 [20]. Based on the CytoNCA plugin (version 2.1.6, parameter: without weight) [21] 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 [22] in the PPI network were selected. Moreover, the MCODE plugin (version 1.2) [23] 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 [18].

Integrated network analysis
Using WEB-based gene set analysis toolkit (WebGestalt) tool [24], 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) [25] 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 [20].

DEG analysis and hierarchical cluster analysis
Compared with the control OS samples, there were a total of 871 DEGs (438 upregulated and 433  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  . 3 The most significant modules identified from the protein-protein interaction (PPI) networks for the upregulated genes (a) and the downregulated genes (b). The significant modules were analyzed by using MCODE plugin, and then module networks were visualized using Cytoscape software Table 2 The GO_BP terms and pathways enriched for the upregulated genes (A) and the downregulated genes (B) involved in the most significant modules.

GO, Gene
Ontology; BP, biological process

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.

Discussion
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 Fig. 4 The miRNA-gene regulatory network for the upregulated genes. The interactions of miRNAs-gene were predicted by WEB-based gene set analysis toolkit tool, and regulatory network was visualized using Cytoscape software. Red circles and white quadrangles represent upregulated genes and miRNAs, respectively 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 [26]. NURR1 and NURR77 may contribute to increase the migratory potential of fetal FBMSCs, which may mediate the local immune response specifically [27]. 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 [30]. These declared that PTHR1 might function in OS through targeting NR4A2.
Through regulating the expressions of versican, HAS2, and hyaluronan, transforming growth factor β2 (TGF-β2) may lead to the metastasis of OS cells [32]. HAS2 plays critical roles in osteoblast differentiation and development by mediating high molecular weight hyaluronan synthesis [33]. RHOC and MMP9 expression levels have close association, and their high expressions are closely correlated with the formation, development, invasion, and metastasis of OS [34]. 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 [35]. Therefore, HAS2 and RHOC might also be targets of PTHR1 in OS.
Via regulating the expression of ZFPM2, hypoxiainduced miR-429 contributes to the differentiation of osteoblastic cells [36]. LEF1 can delay osteoblast maturation and regulates the expression levels of some genes in osteoblasts [37]. LEF1 has an essential role Fig. 5 The miRNA-gene regulatory network for the downregulated genes visualized using Cytoscape software. The interactions of miRNAs-gene were predicted by WEB-based gene set analysis toolkit tool, and regulatory network was visualized using Cytoscape software. Green circles and white quadrangles represent downregulated genes and miRNAs, respectively in the activation of alpha 1 (XI) collagen (COL11A1), and COL11A1 inhibits the terminal differentiation of osteoblasts [38]. LEF1 serves as a transcriptional effector of the Wnt/β-catenin pathway and is critical for the tumor invasion induced by hepatocyte growth factor (HGF) [39]. 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.

Conclusions
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 Fig. 6 The integrated network for the downregulated genes visualized by Cytoscape software. The miRNA-gene regulatory relationships in the network were predicted by using WEB-based gene set analysis toolkit tool, whereas transcription factor-genes regulatory relationships in the network were predicted by using iRegulon plugin. Green circles, white quadrangles, and white triangles represent downregulated genes, miRNAs, and transcription factors, respectively should be conducted to confirm these findings obtained from bioinformatics analysis.