Transcriptomic analyses reveal the underlying pro-malignant functions of PTHR1 for osteosarcoma via activation of Wnt and angiogenesis pathways

Background Increasing evidence has indicated parathyroid hormone type 1 receptor (PTHR1) plays important roles for the development and progression of osteosarcoma (OS). However, its function mechanisms remain unclear. The goal of this study was to further illuminate the roles of PTHR1 in OS using microarray data. Methods Microarray data were available from the Gene Expression Omnibus database under the accession number GSE46861, including six tumors from mice with PTHR1 knockdown (PTHR1.358) and six tumors from mice with control knockdown (Ren.1309). Differentially expressed genes (DEGs) between PTHR1.358 and Ren.1309 were identified using the LIMMA method, and then, protein–protein interaction (PPI) network was constructed using data from STRING database to screen crucial genes associated with PTHR1. KEGG pathway enrichment analysis was performed to investigate the underlying functions of DEGs using DAVID tool. Results A total of 1163 genes were identified as DEGs, including 617 downregulated (Lef1, lymphoid enhancer-binding factor 1) and 546 upregulated genes (Dkk1, Dickkopf-related protein 1). KEGG enrichment analysis indicated upregulated DEGs were involved in Renin-angiotensin system (e.g., Agt, angiotensinogen) and Wnt signaling pathway (e.g., Dkk1), while downregulated DEGs participated in Basal cell carcinoma (e.g., Lef1). A PPI network (534 nodes and 2830 edges) was constructed, in which Agt gene was demonstrated to be the hub gene and its interactive genes (e.g., CCR3, CC chemokine receptor 3; and CCL9, chemokine CC chemokine ligand 9) were inflammation related. Conclusions Our present study preliminarily reveals the pro-malignant effects of PTHR1 in OS cells may be mediated by activating Wnt, angiogenesis, and inflammation pathways via changing the expressions of the crucial enriched genes (Dkk1, Lef1, Agt-CCR3, and Agt-CCL9). Electronic supplementary material The online version of this article (10.1186/s13018-017-0664-2) contains supplementary material, which is available to authorized users.


Background
Osteosarcoma (OS) is the most frequent primary malignant bone tumor developed in the metaphyses of long bones during childhood and adolescence, with an estimated incidence of approximately 3.5 per million [1,2]. Despite intensive multiagent chemotherapy and surgical resection have dramatically increased the 5-year survival rate to 70%, death still occurs in about 30% of patients with OS due to recurrence and metastasis (specially to the lung) [1,2]. Thus, improving understanding of the mechanisms of OS progression and exploiting underlying strategies for malignancy suppression has justifiably attracted a great deal of attention.
Recently, accumulating evidence has indicated parathyroid hormone type 1 receptor (PTHR1), a G-proteincoupled receptor, may play important roles in the pathogenesis of OS. PTHR1 is found to be highly expressed in OS cells and tissues (especially in metastatic or relapsed samples) [3][4][5][6][7]. Over-expression of PHTR1 promotes proliferation, motility, and invasion of OS cells, which can be reversed by shRNA-mediated gene silencing [3,7,8]. Further studies suggest PHTR1 may exert the tumorpromoting effects through being activated by its ligands, including parathyroid hormone (PTH) and parathyroid hormone-related peptide (PTHrP) [9,10]. Upon activation by PTH/PTHrP, PTHR1 induces the generation of cyclic AMP (cAMP) from ATP through adenylyl cyclase followed by the release of cAMP-dependent protein kinases [9][10][11]. Active protein kinases (PKA, PKC, or ERK) move to the nucleus and phosphorylates transcription factors, such as cAMP-response element-binding protein (CREB) and runt-related transcription factor 2 (Runx-2) which ultimately lead to the development of OS through regulating the expression their target genes (TGF-b1, transforming growth factor b1; CTGF, connective tissue growth factor; FGF-2, fibroblast growth factor; HAS2, HA-synthase-2 [3,12,13]). However, the functions of PTHR1 in OS remain not fully understood.
The goal of this study was to further illuminate the mechanisms of PTHR1 by analyzing the microarray data of OS [8]. Differentially expressed genes (DEGs) between OS tissues with and without PTHR1 knockdown were identified and then protein-protein interaction (PPI) network was constructed to screen crucial genes associated with PTHR1, which was not performed in the study of Ho et al. [8]. Our studies may provide new insights into the mechanisms of PTHR1 in OS and reveal some potential targets for treatment of OS.

Microarray data
The microarray data of OS were extracted from the Gene Expression Omnibus (GEO) database (http:// www.ncbi.nlm.nih.gov/geo/) under the accession number GSE46861 [8], which contained six tumors with shRNA PTHR1 knockdown and six tumors with shRNA control knockdown. The tumor tissues were obtained from Balb/c nu/nu mice undergoing mouse OS80 cell line injection into the back flank and grown for 4 weeks. Mouse OS80 was transfected with either renilla luciferase shRNA control (Ren.1309) or a shRNA specific for PTHR1 (PTHR1.358). Thus, PTHR1.358 and Ren.13096 cell samples were used to descript these two groups in the following analysis.

Data normalization and DEG identification
The raw data (CEL files) downloaded from the Affymetrix Mouse Gene 1.0 ST Array platform GPL6246 were preprocessed and normalized using the Robust Multichip Average (RMA) algorithm [14] as implemented in the Bioconductor R package (http://www.bioconductor.org/ packages/release/bioc/html/affy.html). The DEGs between PTHR1.358 and Ren.13096 cell samples were identified using the Linear Models for Microarray data (LIMMA) method [15] in the Bioconductor R package (http:// www.bioconductor.org/packages/release/bioc/html/limma. html). After the t test, the p value was corrected with the Benjamini-Hochberg (BH) algorithm [16]. Genes with an adjusted p < 0.05 and |logFC(fold change)| > 0.5 were considered differentially expressed.

PPI network construction
To screen crucial genes associated with PTHR1, the DEGs were imported into the PPI data that were collected from acknowledged STRING 10.0 (Search Tool for the Retrieval of Interacting Genes; http://stringdb.org/) database [17]. The PPIs with combined scores > 0.7 were selected to construct the PPI network which was visualized using Cytoscape software 2.8 (www.cytoscape.org/) [18]. Three topological properties, including degree [the number of interactions per node (protein)], betweenness (the number of shortest paths that pass through each node), and closeness centrality (the average length of the shortest paths to access all other proteins in the network) were calculated using the CytoNCA plugin in cytoscape software (http://apps.cytoscape.org/ apps/cytonca) [19] to rank the nodes in the PPI network. In addition, the Molecular Complex Detection (MCODE) plugin of Cytoscape software was also employed to identify functionally related and highly interconnected clusters from the PPI network with a degree cutoff of 2, node score cutoff of 0.2, k-core of 2, and maximum depth of 100 (http://baderlab.org/Software/MCODE) [20]. Significant modules were identified with MCODE score ≥ 4 and nodes ≥ 6.

Function enrichment analysis
Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses were performed to investigate the underlying functions of all DEGs and the DEGs in PPI network using The Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 online tool (http://david.abcc.ncifcrf.gov). A modified Fisher Exact p value < 0.05 was chosen as the cutoff point for GO and KEGG analyses.

Function enrichment analysis
The above differential genes were subjected to the online tool DAVID for function enrichment analysis with the mouse genome as background and p < 0.05 as the cutoff point. As a result, 24 KEGG pathways were enriched for upregulated DEGs, including Renin-angiotensin system (e.g., Agt, angiotensinogen), Renin secretion (e.g., Agt), and Wnt signaling pathway (e.g., Dkk1), while 16 pathways were for downregulated DEGs, including basal cell carcinoma (e.g., Lef1) ( Table 2).
In addition, several GO terms, including 909 biological process (GO-BP), 63 cellular component (GO-CC), and 109 molecular function (GO-MF) categories were also enriched for upregulated DEGs, while 578 GO-BP, 42 GO-CC, and 48 GO-MF categories were for downregulated DEGs. To simplify the results, only the GO terms containing PTHR1 gene was displayed in this study ( Fig. 1) because no KEGG pathway was obtained for PTHR1 gene. As expected, PTHR1 was found to be involved in cell proliferation process.

PPI network construction
A PPI network, including 534 nodes and 2830 edges (interaction relationships), was constructed after mapping the DEGs into the PPI data ( Fig. 2; Additional file 1). By calculating the degree, betweenness, and closeness centrality, Agt gene was found to be the most key hub gene (Table 3). More interestingly, Agt was shown to interact with PTHR1 in PPI network, further indicating PTHR1 may promote the development of OS by influencing the expression of this gene. The importance of this gene was also confirmed in the module analysis (Fig. 3). Five modules were screened according to the given parameters (Table 4), among which module 1 (including Agt) was considered as the most significant with MCODE score = 8 and nodes = 17. Function enrichment analysis of module 1 ( Table 2) showed chemokine-and cytokine-related inflammation pathways may be crucial, in which all enriched genes (CCR5, CXCL13, GNAI1, CCR3, CCR2, CCL9) could interact with Agt gene (Fig. 3; Additional file 1), indirectly illustrating the important role of Agt in OS.

Discussion
Using the microarray data of OS provided by Ho et al. [8], we found PTHR1 knockdown could induce the upregulation of Dkk1, but the downregulation of Lef1. Dkk1 is thought to act as a soluble inhibitor for Wnt signaling [21], while transcription factor Lef1 mediates Wnt signaling pathway by binding with its co-activator β-catenin [22]. Several studies have demonstrated that activation of Wnt signaling promotes OS cell proliferation and invasion [23], but contrast results can be obtained after its inhibition [24,25]. Accordingly, we believe Wnt pathway genes (Dkk1 and Lef1) may be an important downstream targets for PTHR1 to participate in the proliferation and invasion of OS, which was also identified in the study of Ho et al. [8]. Although accumulating evidence has confirmed the high expression of Lef1 regulates cell proliferation, migration, invasion, and cancer stem-like cell self-renewal, leading to poor prognosis of patients [26,27], their roles in OS remain rarely reported, and thus, this gene may be a new target for further exploration. The role of Dkk1 in OS remains still  controversial. In contrast to the theoretical expectation [28,29], as well as our result (lower expression in OS), some scholars recently have identified the elevated expression of Dkk1 in OS tissues and cells [30,31] and blockage of Dkk1 via a monoclonal antibody inhibits OS metastasis [32]. This indicates DKK1 represents a class of Janus-faced molecules with dichotomous roles in OS.
We hypothesize the underlying mechanisms may be related with the status of p53 in OS. It has been reported that Dkk-1 can be induced by wild-type p53, but not by mutant p53 (R249S) [33]. Thus, the downregulation of p53 in OS with wild-type p53 may lead to the lower expression of Dkk-1, while Dkk-1 may be increased in a p53-independent manner for OS initiation and maintenance when p53 mutant occurs, which is similar to the regulatory mechanism between p53 and p21 in cancer [34]. Also, a recent study indicates exogenous introduction of p53 and Dkk1 could obviously inhibit the growth of OS cells, cause the cell cycle arrest at G0/G1 phase and apoptosis of OS cells compared with Dkk1 and p53 alone [35], further predicting a synergic relation between p53 and Dkk1. Zhang et al. further found the antiproliferative effects of ursolic acid in OS cells may be mediated by upregulating p53 and then inhibit Wnt/βcatenin signaling [25]. In addition, p53 loss is observed to activate PTHrP-cAMP-CREB1 signaling [11] which is the downstream molecule of PTHR1 and thus may downregulate Dkk1 for OS as our study reported. However, further studies are also needed to confirm this mechanism of p53-PTHR1-Dkk1 in OS.  Furthermore, Agt was also shown to be upregulated after PTHR1 silencing. More interestingly, it could interact with PTHR1 in PPI network, indicating the change in its expression may be a crucial mechanism for explaining the roles of PTHR1 in OS, which was first identified in our study. Function enrichment analysis proved Agt may be involved in Renin-angiotensin system. As is well known, angiogenesis is an indispensable process for tumor growth and metastatic dissemination via providing essential oxygen and nutrients to proliferating cells and then a route for metastasis delivery [36]. Thus, targeted inhibition of angiogenesis may be potential approaches for prevention of OS progression. AGT, encoded by Agt gene, is a 452-amino-acid-residue protein that can be cleaved by renin to generate angiotensin I (AngI) which has been demonstrated to exert antiangiogenic properties in vitro and in vivo [37], suggesting the underlying anti-tumor activity of AGT. This conclusion has been further verified by recent studies. For example, Bouquet et al. showed adenovirus-mediated Agt overexpression inhibited tumor growth in preestablished human MDA-MB-231 mammary carcinomas in nude mice compared to controls and blocked tumorigenicity and pulmonary metastases of MDA-MB-231 and murine melanoma B16F10 cells when they were injected into C57BL/6 mice [38]. Vincent et al. revealed mice with bitransgenic HCC (hepatocellular carcinoma)/Hu-AGT-TG exhibited a significantly longer survival time than the HCC-TG mice and a decrease in both tumor growth and blood flow velocities of the liver through reducing of endothelial arterial markers (active Notch4, Delta-like 4 ligand and ephrin B2) [39]. However, the mechanism of Agt gene in OS remains unclear. In this  Degree the number of interactions per node (protein), betweenness the number of shortest paths that pass through each node, closeness centrality the average length of the shortest paths to access all other proteins in the network Fig. 3 The most significant module extracted from protein-protein interaction network. The red and green nodes represent the upregulated and downregulated genes, respectively study, we also predicted Agt might exert anti-malignancy effects by interacting with inflammation-related genes (such as CCR3, CC chemokine receptor 3; and CCL9, chemokine CC chemokine ligand 9). The relationship between inflammation genes and cancer development has been extensively studied. For example, it has been reported that CCR3 is highly expressed in breast cancer samples, especially luminal-like subtype [40]. Knockdown of CCR3 inhibited cellular proliferation, invasion, and migration, which was ERK signaling pathway-dependent [41,42]. CCL9 was also shown to be highly induced in Gr-1 + CD11b + immature myeloid cells and premetastatic lung of tumor-bearing mice. Knockdown of CCL9 in myeloid cells reduced tumor cell proliferation and metastasis [43].

Conclusion
Our present study preliminarily reveals PTHR1 may play important roles in the development and progression of OS by activating Wnt (Dkk1 and Lef1), angiogenesis, and inflammation pathways (Agt-CCR3 and Agt-CCL9). Lef1, Agt, CCR3, and CCL9 are all underlying new targets because no studies focused on them in OS. Thus, further in vitro and in vivo experimental studies were necessary to confirm the above findings. In addition, the role of Dkk1 is controversial and whether its expression is dependent on p53 status in OS also needs further investigation.