To Explore the Mechanism of the GRM4 Gene in Osteosarcoma by RNA Sequencing and Bioinformatics Approach.

BACKGROUND Glutamate metabotropic receptor 4 (GRM4) has been correlated with the pathogenesis of osteosarcoma. The objective of this study was to explore the underlying molecular mechanism of GRM4 in osteosarcoma. MATERIAL AND METHODS The expression levels of GRM4 in four human osteosarcoma cell lines and hFOB1.19 cells were examined by real-time quantitative PCR (RT-qPCR). The U2OS cells of the highest GRM4 expression were transfected with lentivirus-mediated small interfering RNA (siRNA). The differentially expressed genes (DEGs) after GRM4 gene silencing were screened through RNA sequencing, and analyzed by bioinformatics. Additionally, the transcription factors (TFs) targeting GRM4 were predicted and the downstream protein-protein interaction (PPI) network was constructed using the bioinformatics approach. RESULTS A total of 51 significant DEGs were obtained, including 14 upregulated and 37 downregulated DEGs. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs indicated that four significant enrichment pathways were obtained. A total of six TFs that could be involved in the transcriptional regulation of GRM4 were detected. The results showed that 182 genes in the PPI network were significantly enriched in 14 pathways. The chemokines and chemokine receptors were found to be significantly enriched in three pathways. CONCLUSIONS The DEGs in the four significant enrichment pathways might participate in the development and progression of osteosarcoma through GRM4. The results revealed that EGR1 and CTCF are probably involved in the transcriptional regulation of GRM4, which participates in the progress of osteosarcoma by interacting with chemokines and their receptors.

Background

Osteosarcoma (OS) is the most prevalent primary malignant tumor of mesenchymal origin, arising in the metaphysis of long bone primarily affecting children and adolescents [1]. Neoadjuvant chemotherapy and surgical techniques have improved the 5-year survival rate to up to 70% in OS patients, however, the overall 5-year survival rate is around 20% to 30% in OS patients with metastatic disease at diagnosis, and the treatment protocol is controversial as the exact molecular mechanism of the development and progression of OS is not entirely clear [2,3]. Therefore, it is necessary to explore the molecular mechanisms of OS in order to improve therapeutic treatments.

The mGluRs are G-protein-coupled receptors, activated by the ligands that stimulate secondary messengers, including inositol triphosphate, diacylglycerol, and cyclic AMP (cAMP), to generate the signaling effect [4] and are mainly involved in maintaining the stability of the intracellular environment of the central nervous system. According to sequence similarity, pharmacology, and signal transduction mechanism, the mGluRs can be subdivided into three groups [5]. Group I (mGluR1 and mGluR5) has been shown to activate the phospholipase C pathway, leading to the hydrolysis of phosphatidylinositol and mobilization of intracellular Ca2+. Group II (mGluR2 and mGluR3) and group III (mGluR4, mGluR6, mGluR7, and mGluR8) couple negatively with adenylyl cyclase, suppressing the formation of cAMP, and further inhibiting protein kinase A [5,6].

Mounting evidence has revealed that GRM4 is correlated with poor prognosis in various cancers such as colorectal cancer [7], malignant glioma [8], rhabdomyosarcoma, and multiple myeloma [9]. Recently, a genome-wide association study suggested that mGluR4 may be involved in the pathogenesis of OS [10]. Moreover, in our previous study, we identified that the mGluR4 expression level was higher than that of normal bone tissues, and high expression of mGluR4 was positively correlated with higher Enneking stage, tumor metastasis, and poor prognosis in OS patients [11]. However, another study reported that the expression of mGluR4 was positively correlated with good prognosis in OS [12]. Unfortunately, the oncogenic details of the mechanism of GRM4 remains largely unknown. Therefore, in this study, we first selected the human OS cells with the highest expression of GRM4 from four human OS cell lines and hFOB1.19. Then the GRM4 gene was silenced via lentivirus-mediated siRNA. Next, the DEGs between U2OS cells and GRM4-siRNA cells were observed by RNA sequencing. Finally, we explored the underlying molecular mechanism of GRM4 based on bioinformatics approach.

Material and Methods

Cell culture

The human OS cell lines HOS, MG-63, U2OS, and Saos-2, the human osteoblast cell lines hFOB1.19, and human renal epithelial 293T cells were obtained from the Shanghai Cell Bank (Shanghai, China), and cultured in Dulbecco’s modified Eagle’s medium (DMEM, NY, USA) supplemented with 10% fetal bovine serum (FBS), 100 μg/mL streptomycin and 100 U/mL penicillin. All the cells were incubated at 37°C and 5% CO2 in a humidified atmosphere.

Real-time quantitative PCR

Total RNA was extracted from cultured cells using BiooPure™ RNA Isolation reagent (Bioo Scientific, TX, USA) for reverse transcription and cDNA was synthesized using the K1622 RevertAidTMFirst Strand cDNA Synthesis Kit (Fermentas) according to the manufacturer’s instruction. Real-time PCR was performed using Maxima SYBR Green/Rox qPCR master Mix (2x) (Fermentas) on an ABI Stepone Plus (USA) using the following primers:

GRM4-F, 5′-TGAGGGTGCTGTCACGATCC-3′;

GRM4-R, 5′-ACGTGGCTGCCCTTCTTGAG-3′;

β-actin-F, 5′-ACAGAGCCTCGCCTTTGCCGAT-3′;

β-actin-R, 5′-CTTGCACATGCCGGAGCCGTT-3′.

After 10 minutes at 95°C, DNA polymerase was activated, the cycling parameters were set as follows: one cycle of 95°C for 10 minutes, 40 cycles consisting of denaturation at 95°C for 15 seconds, annealing at 60°C for 30 seconds, and extension at 72°C for 15 seconds. β-actin was used as an endogenous reference for normalization. Gene expression was calculated using the 2−ΔΔCt method.

Design and synthesis of sequences of siRNAs and selection for the best silencing effect sequence of siRNA

Three coding regions corresponding to targeting the GRM4 were designed as siRNA target sequences under the guide of siRNA designing software offered by Genscript. These interfering sequences targeting the GRM4 mRNA (GeneBank accession: NM_000841) were designed and named siRNA-1 (GACCTTCAATGAGAATGGA), siRNA-2 (CACCTTAGAATAGAGCGGA) and siRNA-3 (CCCATCATCAAGCTTGAGT), whose coding regions corresponded directly to the GRM4 mRNA starting at 1889, 1998, and 2229 respectively. And a negative control sequence (NC-siRNA) was also designed. Based on the sequence of siRNAs and NC-siRNA (Table 1), the transcription templates of siRNAs specific for GRM4 were chemically synthesized by Sagene Co., Ltd. (Guangzhou, China). Following that, transfections were performed with Lipofectamine™ 2000 Reagent (Invitrogen) according to the manufacturer’s instruction. In order to detect the silencing effects of different siRNA sequences, the expression of GRM4 mRNA was determined using RT-qPCR.

Table 1

Sequences of specific siRNAs and negative control siRNA targeting GRM4 gene.

Name	Sense strand	Sequence
siRNA-1	Positive-sense strand	5′-GACCUUCAAUGAGAAUGGATT-3′
	Anti-sense strand	5′-UCCAUUCUCAUUGAAGGUCTT-3′
siRNA-2	Positive-sense strand	5′-CACCUUAGAAUAGAGCGGATT-3′
	Anti-sense strand	5′-UCCGCUCUAUUCUAAGGUGTT-3′
siRNA-3	Positive-sense strand	5′-CCCAUCAUCAAGCUUGAGUTT-3′
	Anti-sense strand	5′-ACUCAAGCUUGAUGAUGGGTT-3′
NC-siRNA	Positive-sense strand	5′-UUCUCCGAACGUGUCACGUTT-3′
	Anti-sense strand	5′-ACGUGACACGUUCGGAGAATT-3′

Construction of lentivirus vectors, transfection and selection for stable transfectants

The siRNA (GACCTTCAATGAGAATGGA) targeting the GRM4 mRNA (GeneBank accession: NM_000841) was designed, and scrambled shRNA sequence (TTCTCCGAACGTGTCACGTTTC) was used as a negative control (NC-siRNA). Based on the sequence of siRNA, the short hairpin RNA (shRNA) specific for GRM4 was designed and chemically synthesized by Sagene Co., Ltd. (Guangzhou, China). The sequence of shRNA was inserted into the BamHI and EcoRI sites of the pGLV3/H1/GFP+Puro Vector (Hyclone, Shanghai, China). The accuracy of the recombinant vectors was confirmed by restriction enzyme analysis and sequencing, namely GRM4-siRNA. Recombinant lentivirus vectors and negative control lentivirus vectors were produced by co-transfecting with the lenti-X™ HTX Packaging System (Catalog No. 631247, Clontech, USA) in 293T cells according to the manufacturer’s instruction. The lentiviruses were harvested at 48 hours post-transfection, then centrifuged and filtered through 0.45 μm cellulose acetate filters. The infectious titer was determined by hole-by-dilution titer assay. For transfection, U2OS cells were seeded into six-well plates at a density of 5.0×105 per well supplemented with 5 μg/mL Polybrene (Sagene Co., Ltd., Guangzhou, China), and transfected with recombinant GRM4-siRNA and NC-siRNA lentivirus at a multiplicity of infection (MOI) of 3, respectively. At 72 hours post-transfection, stable GRM4 knockdown cell clones were obtained by 5 μg/mL puromycin (Amersoco) selection. The individual puromycin-resistant clones were picked and then identified for GRM4 expression by RT-qPCR and western blot.

Western blot

The cultured cells were lysed with pre-cooling 4°C lysis buffer (0.004% bromophenol blue, 125 mM Tris-HCl (pH 6.8), 10% 2-mercaptoethanol, 4% sodium dodecyl sulfate, and 20% glycerol) for 60 minutes on ice. The protein concentration was detected using the bicinchoninic acid assay (BCA) method, and then placed in boiling water for 10 minutes. The same amount of protein (20 μg) from each sample per well underwent 12% SDS-PAGE and then was transferred onto a polyvinylidene difluoride (PVDF) membrane and blocked with TBST containing 5% dried skim milk for 60 minutes at room temperature. The PVDF membranes were incubated with rabbit polyclonal antibody of GRM4 (1: 1,000) overnight at 4°C and then with horseradish peroxidase-labeled goat anti-rabbit IgG (1: 5,000) for 60 minutes at room temperature. The result was observed after exposure of x-ray radiographs, development and fixation. Detection was performed with Uvipro image analysis system. All values were normalized against GAPDH values.

Preparation of cDNA library and RNA sequencing

Total RNA was isolated from U2OS cells and U2OS cells with GRM4 stable knockdown with TRIzol Reagent according to the manufacturer’s instruction (Invitrogen). Libraries were prepared according to the standard Illumina mRNA library preparation (Illumina). Purified mRNA was fragmented in fragmentation buffer. Using these short fragments as templates, random hexamer primers were used to synthesize the first-strand cDNA. The second-strand cDNA was synthesized by using buffer, dNTPs, DNA polymerase I, and RNase H. The short double-stranded cDNAs were purified with QiaQuick PCR extraction kit (Qiagen) and resolved with EB buffer for end repair, and addition of an “A” base, and the short fragments were ligated to Illumina sequencing adaptors. Next, DNA fragments with selected size were gel-purified with QiaQuick PCR extraction kit (Qiagen) and amplified by PCR. Constructed libraries were sequenced using Illumina HiSeq™ 2000 sequencing machine and the sequencing strategy was paired-end sequencing.

Processing of RNA sequencing data

The raw reads were filtered to generate clean reads and then saved in the high-quality fastq format. The raw reads were filtered on the basis of the following criteria: Reads with sequencing adaptors were removed, nucleotides with low-quality score were removed, artificial reads were discarded, and reads shorter than 50 were removed. The clean reads were used for subsequent analysis, and then aligned to the reference transcriptome sequences and reference genome by Bowtie2 and Tophat 2 [13]; reads matched with reference rRNA sequences were also mapped and discarded.

Testing of differentially expressed genes and pathway enrichment analysis

The expression of genes was measured and normalized using fragments per kilobase of exon per mapped million fragments (FPKM). The DEGs between U2OS cells and U2OS cells with GRM4 stable knockdown were identified by edger [14]. The selection criteria are the p-value ≤0.05 and fold change ≥2, and biological coefficient of variation (BCV)=0.4 according to the edgeR Users Guide. The p-value ≤0.05 and DEGs count ≥2 were used as cutoff criterion to identify the significantly enriched KEGG pathways by the DAVID bioinformatics tool.

Prediction of transcription factors and construction of PPI network

Based on the human TF-gene pairs obtained from the UCSC database, we used the bioinformatics method to select TFs that targeted the regulation of GRM4. PPIs were extracted via retrieving high-confidence protein interactions with GRM4 from the Search Tool for the Retrieval of Interacting Genes database (STRING) [15]. The combined score ≥0.9 was set as the threshold for construction of PPI network, and the corresponding network was visualized using Cytoscape software [16].

Functional annotation and pathway analysis of the genes identified by the PPI network

The DAVID bioinformatics tool was used to perform Gene Ontology (GO) and KEGG pathway analysis of the genes identified via the PPI network. The screening criterions to select the significant enrichment GO terms or KEGG pathways were set as p-value <0.05 and enriched genes count ≥2.

Statistical analysis

Statistical analysis was performed with SPSS 19.0 software (SPSS Inc., USA). All experiments were repeated in triplicate and data are expressed as mean ± standard deviation. The t-test was used for comparison between two groups. A value of p<0.05 was considered statistically significant.

Results

Screening the cell lines with the highest expression of GRM4 from OS cell lines and the human osteoblast hFOB1.19 cells and expressions of GRM4 after GRM4 silencing in U2OS cells

RT-qPCR showed that relative high expression of GRM4 was observed in OS cell lines MG-63, U2OS, HOS, and Saos-2 compared with human osteoblasts hFOB1.19 cell lines, and the expression of GRM4 in U2OS cells was the highest (Figure 1A). Compared with those cells transfected with NC-siRNA, the GRM4 mRNA level in U2OS cells that were transfected with siRNA-1, siRNA-2, and siRNA-3 were significantly reduced by 75%, 50%, and 68% respectively (Figure 1B). Therefore, the silencing effect of siRNA-1 was the best among the three sequences. Stable transformants were obtained by puromycin selection, and then identified by RT-qPCR and western blot. The results indicated that the GRM4 mRNA level in U2OS cells that transfected with GRM4-siRNA lentivirus were significantly reduced by 80%, compared with those transfected with NC-siRNA lentivirus (p=0.000, Figure 1C). The protein level of GRM4 was also significantly inhibited by GRM4-siRNA lentivirus transfection (Figure 1D).

Figure 1

Expressions of GRM4 mRNA in OS cell lines and the human osteoblast hFOB1.19 cells and expressions of GRM4 after GRM4 silencing in U2OS cells. (A) GRM4 mRNA in each cell line was detected by RT-qPCR. (B) Expressions of GRM4 mRNA after GRM4 silencing by three siRNAs in U2OS cells. (C) The mRNA level of GRM4 was significantly reduced in GRM4-siRNA treated group compared with NC-siRNA treated group, * p=0.000. (D) Western blot revealing that the protein levels of GRM4 was significantly inhibited by GRM4-siRNA lentivirus transfection.

Characterization of RNA sequencing data and analysis of differently expressed genes

U2OS cells and GRM4-siRNA cells were used for RNA sequencing. After mapping the data to reference genome, 40.3 million and 26.3 million reads pairs were obtained in U2OS cells and GRM4-siRNA cells respectively. The normalized expression level of each gene was measured by FPKM and the gene expression between the samples were highly correlated with Pearson correlation coefficient with a value of 0.8018. A total of 51 significant DEGs were obtained, including 14 upregulated and 37 downregulated DEGs. The clustering analysis showed that the DEGs of U2OS cells were distinct from that of GRM4-siRNA cells (Figure 2).

Figure 2

Heat map and unsupervised hierarchical clustering analysis.

Pathway enrichment analysis of differentially expressed genes

In total, there were four significant enrichment pathways (Table 2), including transcriptional misregulation in cancers, PPAR signaling pathway, complement and coagulation cascades, and rheumatoid arthritis.

Table 2

KEGG pathway analysis of differentially expressed genes.

Term	Pathway	Count	P-value	Genes
hsa05202	Transcriptional misregulation in cancers	3	0.021946	TRAF1, MMP3, HMGA2
hsa03320	PPAR signaling pathway	2	0.024251	MMP1, HMGCS2
hsa04610	Complement and coagulation cascades	2	0.028817	CFI, C4B
hsa05323	Rheumatoid arthritis	2	0.037374	MMP1, MMP3

TFs involved in the regulation of GRM4 and functional enrichment analysis of the genes identified by PPI network

A total of six TFs that could be involved in the transcriptional regulation of GRM4 were detected using the UCSC database. The results showed that 182 genes in PPI network were significantly enriched in 14 pathways (Figure 3, Table 3). The 182 genes were then classified according to the relevant biological functions of cellular component (CC), molecular function (MF), and biological progress (BP) (Figure 4). In our present study, the chemokines and chemokine receptors (CXCL1–3, 5–6, 9–13, 16; CXCR1–6; CCR1–10, 25, 27–28; and CCL5, 16, 19–21, 25, 27–28) were found to be significantly enriched in chemokine signaling pathway, cytokine-cytokine receptor interaction, and intestinal immune network for IgA production. Moreover, in the GO category of MF, the chemokines and chemokine receptors were mainly enriched in the chemokine receptor binding, chemokine activity, chemokine receptor activity, chemokine binding, C-C chemokine receptor activity, C-C chemokine binding and cytokine activity.

Figure 3

Integrated network of TF-target and PPI. The TFs implicated in the transcriptional regulation of GRM4 and the proteins interacted with GRM4 are presented in green nodes and pink nodes, respectively.

Table 3

KEGG pathway significantly enriched by the genes associated with GRM4.

Term	Pathway	Count	P-value
hsa04080	Neuroactive ligand-receptor interaction	71	5.47E-55
hsa04062	Chemokine signaling pathway*	50	3.96E-36
hsa04060	Cytokine-cytokine receptor interaction*	37	2.96E-16
hsa04742	Taste transduction	17	3.21E-13
hsa04540	Gap junction	15	1.40E-07
hsa04916	Melanogenesis	12	9.53E-05
hsa04914	Progesterone-mediated oocyte maturation	11	1.36E-04
hsa04020	Calcium signaling pathway	14	1.31E-03
hsa04672	Intestinal immune network for IgA production*	7	0.00238142
hsa05414	Dilated cardiomyopathy	9	0.004256388
hsa04912	GnRH signaling pathway	9	0.006242908
hsa04114	Oocyte meiosis	9	0.012242439
hsa04610	Complement and coagulation cascades	7	0.012870673
hsa04270	Vascular smooth muscle contraction	9	0.013554099

Pathways significantly enriched by chemokines and chemokine receptors.

Figure 4

GO analysis of the 182 genes identifed in PPI network. The enriched genes associated with GRM4 were categorized according to the relevant biological functions of molecular function, cellular component and biological process.

Discussion

The combined effects of environmental and genetic factors contribute to the development of OS, especially within a genetic background. The molecular mechanism leading to tumorigenesis or development of OS is unclear. In the past, studies of GRM4 primarily focused on the functions of the central nervous system. Up until now, accumulating evidence has indicated that GRM4 plays a crucial role in OS tumorigenesis or progression. Hence, we analyzed the DEGs after GRM4 silencing in U2OS cells by RNA sequencing.

Four significant enrichment pathways were obtained in our present study. Past research has revealed that TRAF1 was involved in the regulation of differentiation, stress responses, cell survival, and proliferation [17]. The function of TRAF1 is believed to have a disparate role in cancer cell apoptosis and survival [18]. In previous studies, the functions of MMP1 and MMP3 have been reported. MMP1 has been identified as a prognostic candidate marker for various types of malignancies. Furthermore, high levels of MMP1 expression has been indicated to have important roles in OS and metastasis to the lung [19], and have poor prognosis in OS patients [20]. MMP3 expression can be detected in various types of malignant cells. High expression of MMP3 was implicated in migration, invasion and metastasis in malignant tumors, including OS [21]. Nowadays, accumulating evidence has revealed that HMGA2 is involved in the development and progression of tumor by regulating the stemness of cancer stem cell, tumor cell metastasis, and epithelial-mesenchymal transition (EMT) [22]. Additionally, as a tumor oncogene, the expressions of HMGA2 are associated with Enneking stage, tumor size, lung metastasis, proliferation, and apoptosis in OS cells [23,24]. The expression of HMGCS2 is downregulated in Myc-dependent colon and rectum cancers [25]. Even so, high expression of HMGCS2 was correlated with advanced tumor, which indicated worse prognosis for local recurrence-free survival, metastasis-free survival, and disease-free survival in rectal cancer [26]. The coagulation and complement cascades were involved in the development of tumor [27]. Nevertheless, until now, the role of coagulation-related genes CFI and C4B in tumors remained unclear. The DEGs in the four significant enrichment pathways might be involved in the development and progression of OS via GRM4.

As we all know, TFs are essential for the regulation of gene expression, which can provide better clues on the underlying regulatory mechanisms for the targeting genes. Next, we predicted TFs that regulated GRM4 and constructed PPIs based on bioinformatics methods. The six identified TFs, EGR1 and CTCF in particular, were involved in the transcriptional regulation of GRM4. As transcription factors, EGR1 and CTCF are frequently involved in the regulation of cell differentiation, development, proliferation, and apoptosis [28,29]. Expression of EGR1 has been shown to be downregulated in several tumor cells, and forced expression of EGR1 inhibited migration and invasion but did not prevent growth in OS cell [30]. As a tumor suppressor candidate, CTCF may contribute to carcinogenesis directly or indirectly [31].

Furthermore, in the present study, 14 significant pathways were obtained based on the DAVID tool from 182 genes in the PPI network. Interestingly, based on the GO and pathway enrichment analysis of the 182 genes in the PPI network, we found that GRM4 is probably involved in tumorigenesis or progression of OS through the interactions with the chemokines and chemokine receptors. Approximately 50 chemokines and 20 chemokine receptors have been identified in humans. According to the pattern of cysteine residues, chemokines and chemokine receptors are divided into four families: C, CC, CXC, and CX3C, where C represents the cysteine and X represents noncysteine amino acids. Actually, several researches have indicated that chemokines and chemokine receptors were related to the progression, migration, invasion, and metastasis in various tumors [32]. On the contrary, the role of most of the remaining identified genes (such as ADCY1–9, and TAS2R1, 3–5, 16, 39–43, and 60) in tumors remained unclear. For instance, several research studies have indicated that the interaction of chemokines and chemokine receptors were involved in tumorigenesis and progression of cancers. CCL5 contributes to OS angiogenesis by stimulating the expression of VEGF [33]. The expression of CXCL5 was upregulated in OS and high expression of CXCL5 was implicated in migration and invasion in OS [34]. The chemokine-chemokine receptor axis was involved in the control of tumor angiogenesis, growth, invasion, and metastasis. The interaction of CCL5/CCR5 facilitated migration via αvβ3 integrin in OS [35]. Reportedly, CXCL12/CXCR4 axis was implicated in the pulmonary metastasis of OS [36,37], and facilitated OS cells migration via the Wnt/β-catenin and AKT signaling pathway [38]. The IL-8/CXCR1/Akt signaling pathway might participate in the process to improve chemotherapy sensitivity in OS after CXCR1 silencing [39]. Taken together, most of the chemokines and chemokine receptors were involved in the progression, migration, invasion, and metastasis of various tumors. The expression of chemokines and chemokine receptors is changed in many malignant tumors by either activating the oncogenes or inactivating the tumor suppressor genes. Any such change can subsequently lead to anomalous chemokine receptor signaling and thus affect the normal regulation of the chemokines [32]. This is in conformity with our finding that the chemokines and chemokine receptors were mainly enriched in the chemokine receptor binding, chemokine activity, chemokine receptor activity, chemokine binding, C-C chemokine receptor activity, C-C chemokine binding and cytokine activity in the GO category of MF. Unfortunately, the expression of the chemokines and chemokine receptors in our findings has no significant change after GRM4 silencing. This is might be due to the fact that the GRM4 gene plays a different role in different OS cell lines or different developmental stage.

Conclusions

This study identified the DEGs and related significant enrichment pathways in U2OS cells after GRM4 silencing. Moreover, as a role of tumor suppressor, EGR1 and CTCF may participate in the transcriptional regulation of GRM4, which is involved in the tumorigenesis or progression of OS by interacting with the downstream chemokines and their receptors. This research may provide new insight into the underlying molecular mechanisms of GRM4 in OS. However, further experiments are warranted to confirm these results.