Jinghong Yuan1†
Zhao Yuan2†Aifang Ye3Tianlong Wu4
Jingyu Jia5Jia Guo6
Jian Zhang1Tao Li1
Xigao Cheng1,4,5*- 1Department of Orthopaedics, The Second Affiliated Hospital of Nanchang University, Nanchang, China
- 2Clinical Research Center, The Second Affiliated Hospital of Nanchang University, Nanchang, China
- 3Department of Otorhinolaryngology, Jiangxi Provincial Children’s Hospital, Nanchang, China
- 4Institute of Orthopaedics of Jiangxi Province, Nanchang, China
- 5Institute of Minimally Invasive Orthopaedics of Nanchang University, Nanchang University, Nanchang, China
- 6Department of Orthopaedics, Jiangxi Provincial People’s Hospital Affiliated to Nanchang University, Nanchang, China
Background: G protein subunit gamma 12 (GNG12) is observed in some types of cancer, but its role in osteosarcoma is unknown. This study hypothesized that GNG12 may be a potential biomarker and therapeutic target. We aimed to identify an association between GNG12 and osteosarcoma based on the Gene Expression Omnibus and the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) databases.
Methods: Osteosarcoma samples in GSE42352 and TARGET database were selected as the test cohorts. As the external validation cohort, 78 osteosarcoma specimens from The Second Affiliated Hospital of Nanchang University were collected. Patients with osteosarcoma were divided into high and low GNG12 mRNA-expression groups; differentially expressed genes were identified as GNG12-related genes. The biological function of GNG12 was annotated using Gene Ontology, Kyoto Encyclopedia of Genes and Genomes, gene set enrichment analysis, and immune infiltration analysis. Gene expression correlation analysis and competing endogenous RNA regulatory network construction were used to determine potential biological regulatory relationships of GNG12. Overall survival, Kaplan–Meier analysis, and log-rank tests were calculated to determine GNG12 reliability in predicting survival prognosis.
Results: GNG12 expression decreased in osteosarcoma samples. GNG12 was a highly effective biomarker for osteosarcoma [area under the receiver operating characteristic (ROC) curve (AUC) = 0.920], and the results of our Kaplan–Meier analysis indicated that overall survival and progression-free survival differed significantly between low and high GNG-expression group (p < 0.05). Functional analyses indicated that GNG12 may promote osteosarcoma through regulating the endoplasmic reticulum. Expression correlation analysis and competing endogenous RNA network construction showed that HOTTIP/miR-27a-3p may regulate GNG12 expression. Furthermore, the subunit suppresses adaptive immunity via inhibiting M1 and M2 macrophage infiltration. GNG12 was inhibited in metastatic osteosarcoma compared with non-metastatic osteosarcoma, and its expression predicted survival of patients (1, 3, and 5-year AUCs were 0.961, 0.826, and 0.808, respectively).
Conclusion: This study identified GNG12 as a potential biomarker for osteosarcoma prognosis, highlighting its potential as an immunotherapy target.
1 Introduction
Osteosarcoma is among the most common primary solid malignant bone tumors in adolescents and young adults. Usually originating in the metaphyses of long bones, the cancer is characterized by heterogeneous presentation, high mortality, and an annual incidence of 8–11 million among 15–19-year-olds (1). After comprehensive treatment with methods such as radiotherapy and chemotherapy, the 5-year survival rate is 65%–70% (2, 3), but these improvements have been mainly limited to patients with non-metastatic disease. Moreover, the progress of research aiming to raise osteosarcoma survival has stalled over the past 30 years (4). Therefore, it is critical to determine reliable predictors related to osteosarcoma metastasis and prognosis and to make available novel targets for therapy and prognosis prediction. However, due to the complex molecular mechanisms of osteosarcoma, the predictive capability of traditional clinical information is limited. Thus, it is crucial to find novel prognostic biomarkers to predict survival and metastases more accurately in osteosarcoma. Although biomarkers such as FAT10 and MYC have been associated with osteosarcoma in recent studies, their reliability requires further investigation (5, 6).
G protein subunit gamma 12 (GNG12) is a protein-coding gene located on chromosome 1, first reported in 1995 (7). In the GeneCards database (www.genecards.org), related pathways include “Translation Translation regulation by Alpha-1 adrenergic receptors” and “Sweet Taste Signaling”. Gene Ontology (GO) annotations related to this gene include obsolete signal transducer activity and phosphate ion binding. Previous research (8) suggested that GNG12 knockdown in BV-2 cells increased nitric oxide levels and tumor necrosis factor alpha (TNF-α) expression in response to lipopolysaccharide (LPS) stimulation, indicating that GNG12 was an essential negative regulator of inflammation. Transforming growth factor beta (TGF-β)-induced chondrogenic differentiation of mesenchymal stem cells is promoted when C-type natriuretic peptide/natriuretic peptide receptor-B reduces GNG12 expression (9). The first association of GNG12 with cancer was a study showing that low GNG12 expression increases the proliferation of endometrial cancer (10). Subsequently, HOXA13 upregulation was found to be a promoter of lung squamous cancer progression through reducing GNG12 expression (11). Although these studies suggest an important role for GNG12 in cancer progression, we currently know little about the underlying mechanisms and function of this protein in osteosarcoma progression and immunology.
This study thus aimed to determine the role of GNG12 in patients with osteosarcoma. We obtained two test cohort from the Gene Expression Omnibus (GEO) database and the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database. To clarify biological function, we compared expression matrices of high and low GNG12-expression groups to identify potential GNG12-related genes. These genes were screened using functional enrichment analyses, gene set enrichment analysis (GSEA) regulation networks, immune infiltration analysis, and pan-cancer analysis. We also collected osteosarcoma samples to evaluate GNG12 predictive ability on survival and prognosis, using immunohistochemistry, histochemistry score (H-Score) analysis, and receiver operating characteristic (ROC) analysis.
Our results demonstrated that low GNG12 expression is a potential biomarker of osteosarcoma, linked to poor prognosis. GSEA revealed that GNG12 expression was associated with “extracellular matrix organization”, “core matrisome”, “cytoplasmic ribosomal proteins”, “cell adhesion molecules cams”, “class A 1 rhodopsin-like receptors”, “GPCR ligand binding”, and “diseases of metabolism”, and “Matrisome”. Immune infiltration analysis demonstrated that GNG12 regulated the proportions of macrophages (M0, M1, and M2) and mast cells to influence the tumor microenvironment. Finally, we confirmed the feasibility of GNG12 as a biomarker in our collected osteosarcoma samples.
2 Materials and Methods
2.1 Patients and Identification of Differentially Expressed Genes
The GEO (http://www.ncbi.nlm.nih.gov/geo) database is a free public gene-expression data repository containing microarray and high-throughput sequencing data. Gene-expression datasets from GSE42352 (12, 13), including 103 osteosarcoma cases and 15 normal controls, were collected from GEO as the test cohort. Using the median value of GNG12 expression as the cutoff, the 103 patients were divided into high (n = 52) and low (n = 51) GNG12 mRNA expression groups. Significant differentially expressed genes (DEGs) between high and low GNG12-expression groups were screened with the limma package (http://www.bioconductor.org/packages/release/bioc/html/limma.html) in R version 3.6.3 (http://www.R-project.org/) (14). Genes were considered differentially expressed with an adjusted p < 0.5 and |log fold change| (|logFC|) ≥ 1. Heat maps were constructed using R package “pheatmap”, version 1.0.12 (https://cran.rproject.org/web/packages/pheatmap/index.html), and volcano plots were generated with “ggplot2” version 3.3.3 (https://cran.r-project.org/web/packages/ggplot2/index.html). In addition, gene expression profiles and clinical data of osteosarcoma cases were collected from TARGET (https://ocg.cancer.gov/programs/target). The TARGET cohort (n = 101) were divided into high/low groups using the same cutoff as the GEO cohort (high, 51 cases; low, 50 cases). The same method as described above was used to determine DEGs between TARGET high/low groups. Common DEGs between the GSE42352 and TARGET datasets were represented with a Venn diagram made in ggplot2.
The Second Affiliated Hospital of Nanchang University provided 78 osteosarcoma specimens from January 2012 to December 2016 as an external validation cohort. Samples were provided as formalin fixed paraffin-embedded blocks.
This research was approved by the Ethics Committee of The Second Affiliated Hospital of Nanchang University [Review (2020) No. (086)]. All participants provided informed consent.
2.2 Functional Enrichment Analysis
Functional enrichment analyses (Gene Ontology, GO; https://www.geneontology.org and Kyoto Encyclopedia of Genes and Genomes, KEGG; https://www.genome.jp/kegg/) were performed on DEGs between high and low GNG12-expression groups in the test cohort, using clusterProfiler version 3.14.3 (http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html). The results were visualized with the “ggplot2” package (https://cran.r-project.org/web/packages/ggplot2/index.html). Xiantao Academic online (https://www.xiantao.love) was used to generate a diagram of the functional enrichment network analysis.
Significant differences in function and pathways between high/low GNG12-expression groups were determined using GSEA in clusterProfiler and Xiantao Academic online. For each analysis, gene set permutations were performed 1,000 times to obtain a normalized enrichment score (NES). Enrichment was considered significant with adjusted p < 0.05, false discovery rate (FDR) q < 0.25, and |NES| > 1. The selected reference gene set was c2.cp.v7.2.symbols.gmt (Curated), and results were visualized in ggplot2.
2.3 Protein–Protein Interaction Network Construction and Hub-Gene Extraction
A protein–protein interaction (PPI) network of DEGs was constructed using Metascape online tools (https://metascape.org) with the following parameters: min network size = 3 and max network size = 500. Crucial proteins in this PPI network were screened using the Molecular Complex Detection (MCODE, http://apps.cytoscape.org/apps/mcode), a plug-in of Cytoscape version 3.7.2 (https://cytoscape.org/).
2.4 Immune Infiltration Analysis
CIBERSORT (https://cibersort.stanford.edu/) is an algorithm developed to characterize the cellular composition of complex tissues through gene expression profiles compared with a signature matrix (LM22). LM22 comprises 547 genes that define 22 immunization cell subtypes. This analysis was used to calculate infiltration abundance of 22 immune cells (T cells, B cells, plasma cells, natural killer cells, myeloid subgroups) in 103 osteosarcoma samples across high and low GNG12-expression groups. We also analyzed any correlations between the two GNG12 groups in terms of expression distribution among the 22 infiltrating immune cells; the correlation heat map was then plotted. Red represents a positive correlation, and blue represents a negative correlation, with darker colors and values closer to 1 indicating stronger correlations.
2.5 GNG12 Expression Analysis
The UCSC XENA platform (https://xenabrowser.net/datapages/) is a pan-cancer HTSeq-TPM database of TCGA and GTEx, processed uniformly by Toil. Based on XENA data, differential GNG12 mRNA expression between tumor and normal tissues was determined in R and visualized with ggplot2. Spearman’s correlations were calculated between GNG12 expression in the TARGET cohort with the expression of nine DEGs (GREM1, CAV1, PTGDS, PRB2, PRB1, NDUFB9, AMBN, FAT3, and IBSP) common across GSE42352 and TARGET datasets. Again, the analysis was performed in R and the results visualized in ggplot2.
2.5.1 Construction of Competing Endogenous RNA Network
A differential long noncoding RNA (lncRNA) expression matrix was obtained from TARGET using the limma package. The miRcode database (http://www.mircode.org/) was used to predict a series of highly conserved microRNAs (miRNAs) linked to differentially expressed lncRNA. The miRNA target genes were predicted using miRDB (http://mirdb.org/), miRTarBase (http://miRTarBase.cuhk.edu.cn/), and TargetScan (http://www.targetscan.org/mamm_31/) databases. Predicted mRNAs were associated with the six GNG12-related genes (GREM1, CAV1, NDUFB9, AMBN, FAT3, and IBSP) to establish a competing endogenous RNA (ceRNA) network.
2.6 Survival Analysis
Kaplan–Meier (KM) analyses of the TARGET database were performed using R packages survival (https://www.rdocumentation.org/packages/survival/versions/2.42-3) and survminer (https://cran.rstudio.com/web/packages/survminer/index.html), to determine whether GNG12 expression accurately predicted overall survival time (OS) and progression-free survival (PFS). The timeROC package (https://cran.r-project.org/web/packages/timeROC/index.html) was used to analyze 1-, 3-, and 5-year survival prognosis of GNG12-related lncRNAs, miRNAs, and genes. All outcomes were visualized in ggplot2.
2.7 Assessing Diagnostic Performance
An ROC curve was generated using the pROC package (https://cran.r-project.org/web/packages/pROC/). The ordinate was drawn with the true positive rate (sensitivity), and the abscissa was the false positive rate (1 − specificity). This analysis determined whether GNG12 expression can distinguish between 15 standard samples and 103 osteosarcoma samples. Also determined was the cutoff value that produces the highest likelihood ratio for assessing the recognition threshold of GNG12 as a biomarker for distinguishing osteosarcoma from normal tissue.
2.8 Immunohistochemical Staining and Evaluation
To further verify the suitability of GNG12 expression for predicting survival, immunohistochemistry was performed on paraffin sections following the standard protocol (GNG12, Ab204757, 1:100). All slides were observed and photographed under a XSP-C204 microscope (CIC). GNG12 expression was assessed using the H-score.
2.9 Validation of GNG12 as a Prognostic Predictor of Osteosarcoma
Known survival outcomes of osteosarcoma patients with high and low H-scores were subjected to KM analysis using the R package “survival” and visualized in “surviminer”. The accuracy of GNG12 expression in predicting survival prognosis (1, 3, and 5 years) was verified with timeROC and visualized in ggplot2. An AUC >0.5 directly reflects the prognostic value of the biomarker. Additionally, an AUC closer to 1 and to the (0, 1) point indicates greater authenticity of the prognostic prediction.
2.10 Statistical Analysis
Statistics were performed in Microsoft Excel and R version 3.6.3. The log-rank test was used to perform the KM survival analysis. Continuous variables that violated the normality assumption according to the Shapiro–Wilk normality test were analyzed with the Wilcoxon rank sum test. Categorical variables were analyzed with chi-square tests. Significance was set at p < 0.05 (ns, p > 0.05; *p < 0.05; **p < 0.01; ***p < 0.001).
3 Results
3.1 Identification of DEGs
GNG12 expression was significantly different between the standard samples (n = 15) and osteosarcoma samples (n = 103) (Figure 1A). After screening for DEGs between high and low GNG12-expression groups (Figure 1B), we found 210 coexpressed genes (165 upregulated and 45 downregulated, Supplementary Table 1) that we have visualized in a volcano plot (Figure 1D) and a heat map (Figure 1E). Using the same method to divide the TARGET cohort (n = 101) into high/low GNG12-expression groups (Figure 1C), we obtained 347 DEGs (178 upregulated and 169 downregulated, Supplementary Table 2 and Figures 1F, G). We then screened out nine common DEGs (six upregulated and three downregulated) across the two DEG sets (Figures 1H, I).
Figure 1 GNG12 and coexpressed gene difference analysis. (A) Compared with normal tissues, GNG12 is significantly lower expressed in tumor tissues. (B) In GSE42352 cohort, the expression level of GNG12 in low/high GNG12-expression groups. (C) In TARGET cohort, the expression level of GNG12 in low/high GNG12 expression groups. (D) The median expression of GNG12 in 103 osteosarcoma samples in the GSE42352 data set was divided into high and low expression groups, and the significantly different genes between the two groups were displayed in the form of a volcano graph. (E) A heat map showing the significantly different genes between the two groups of the GSE42352 data set. (F) GNG12 was used in the median expression values of 101 osteosarcoma samples in the TARGET data set divided into high and low expression groups, showing the significant difference genes between the two groups in the form of volcano maps. (G) The TARGET data set shows the significant difference genes between the two groups in the form of a heat map. (H) Venn diagram of the intersection of upregulated DEGs in the GEO and TARGET cohort. (I) Venn diagram of the intersection of downregulated DEGs in the GEO and TARGET cohort. ***p < 0.001.
3.2 Functional Enrichment Analysis and GSEA
The GO and KEGG enrichment results for the GEO cohort (Figures 2A, B) indicated that the following terms were enriched: “protein targeting to endoplasmic reticulum (ER)”, “establishment of protein localization to ER”, “focal adhesion”, and “cell adhesion molecule binding” (Figure 2C). For the TARGET cohort, enriched terms were “ossification”, “cartilage development”, “integrin binding”, “PI3K-Akt signaling pathway”, and “Focal adhesion” (Figure 3).
Figure 2 GO and KEGG enrichment analysis in GEO cohort. (A) The bubble graph of clusterProfiler package for GO and KEGG enrichment analysis; the bubble size represents the number of gene enrichment, and color represents significance in GSE42352 data set. (B) The bar graph of clusterProfiler package for GO and KEGG enrichment analysis; length represents significance in GSE42352 data set. (C) A function enrichment network in GSE42352 data set.
Figure 3 GO and KEGG enrichment analysis in TARGET cohort. (A) The bubble graph of clusterProfiler package for GO and KEGG enrichment analysis; the bubble size represents the number of gene enrichment, and color represents significance in TARGET data set. (B) The bar graph of clusterProfiler package for GO and KEGG enrichment analysis; length represents significance in TARGET data set. (C) A function enrichment network in TARGET data set.
We selected msigdb.v7.0.entrez.gmt as the reference gene set for our GSEA. For the GEO cohort, DEGs were involved in “extracellular matrix organization” from the Reactome database, “core matrisome” from NABA, “cytoplasmic ribosomal proteins” from WP, and “cell adhesion molecules cams” from KEGG (Figure 4A). For the TARGET cohort, DEGs were involved in “class A 1 rhodopsin like receptors”, “GPCR ligand binding”, and “diseases of metabolism” from the Reactome database, and “Matrisome” from NABA (Figure 4B).
Figure 4 GSEA enrichment analysis. (A) GSEA enrichment analysis in GEO cohort. (B) GSEA enrichment analysis in TARGET cohort.
3.3 Construction of PPI Networks and Hub-Gene Screening
Based on the Metascape online tool, we constructed two PPI networks of DEGs from the GEO and TARGET cohorts (Figures 5A, C). We then screened hub-gene clusters using the MCODE clustering algorithm. In the GEO cohort, screened hub genes were linked in the following networks: MCODE_1 (cytoplasmic translation, cytosolic ribosome, and ribosome; EIF3A, RPL23, RPL14, RPL7, RPL27A, RPS28, RPLP1, RPS3A, RPL10A, and RPN2), MCODE_2 (NLS-bearing protein import into the nucleus, protein import into the nucleus, and import into the nucleus; KPNA2, RGPD6, RGPD3, RGPD8, SUMO2, ACTR3, WIPF1, ARPC5, and GMFG), and MCODE_3 (ER, growth factor activity, and extracellular matrix structural constituent; VCAN, CCN1, AMBN, CDH2, VGF, and IL6; Figure 5B). Hub genes of the TARGET cohort were involved in more networks. These were MCODE_1 (extracellular matrix structural constituent conferring tensile strength, collagen trimer, and protein digestion and absorption; COL8A2, COL9A1, COL11A1, COL11A2, COL15A1, COL22A1, COL24A1, P4HA3, COL14A1, COL27A1, COL2A1, COL3A1, COL6A1, and COL6A2), MCODE_2 (in utero embryonic development, peptidyl-tyrosine phosphorylation, and peptidyl-tyrosine modification; FGFR4, ISLR, FGF2, IGF2, TGFB3, CAV1, TP53, ACTA2, ALDH1A1, PYGM, GPI, GFPT2, and H2AW), MCODE_3 (G-protein beta-subunit binding, neuropeptide hormone activity, and heterotrimeric G-protein complex; GNG12, GAL, CORT, NPB, GNG4, GNGT2, HCRT, HRH1, PTGFR, ADRA1D, and TBXA2R), MCODE_4 (ER lumen, biomineral tissue development, and biomineralization; MELTF, PENK, MEPE, SPP2, MOTUM, DMP1, FBN1, AMBN, AMELX, and GPC3), MCODE_5 (anchored component of membrane; NTM, ALPL, NTNG1, CD109, and LY6K), MCODE_6 (lamellar body, mutivesicular body and respiratory gaseous exchange by respiratory system; SFTPB, SFTPC, SFTPA1, and SFTPA2), and MCODE_7 (Wnt signaling pathway; CCND2, TCF7L2, and FOSL1; Figure 5D).
Figure 5 PPI and hub gene clusters network construction. (A) A PPI network from GSE42352 data set. (B) Three hub gene clusters obtained by the MCODE clustering algorithm from GSE42352 data set. (C) A PPI network from TARGET data set. (D) Seven hub gene clusters obtained by the MCODE clustering algorithm from TARGET data set.
3.4 Correlation Between GNG12 Expression and Immune Infiltration
For further immune infiltration analysis, we used CIBERSORT package and LM22 algorithm to calculate the infiltration abundance of 22 kinds of immune cells between high/low expression level of GNG12 groups in test and validation cohort, including different B cells, T cells, natural killer cells, plasma cells, and different myeloid subsets. The results were visualized as two violin maps (Figures 6A, B) and two correlation heat maps (Figures 6C, D) to reveal the differences in the expression of immune cells between groups.
Figure 6 Immune cell infiltration difference and correlation analysis. (A) Using the CIBERSORT package, the LM22 algorithm was used to calculate the difference in expression of 22 immune cells in the GNG12 high and low expression groups in the 103 osteosarcoma samples in the GEO set. Blue represents the low expression group, and red represents the high expression group. (B) Violin diagram of the difference in expression of 22 immune cells infiltration between the high and low expression groups of GNG12 in 101 osteosarcoma samples in the TARGET set. Blue represents the low expression group; red represents the high expression group. (C) Correlation heat map of 22 immune cells in the GEO set. (D) Correlation heat map of 22 immune cells in the TARGET set. The p-value represents the significance of the difference; p < 0.05 is considered a significant difference.
3.5 Expression Analysis of GNG12-Related Genes
Evaluation of differential GNG12 expression in pan-cancer revealed that GNG12 was downregulated in 12 tumor types (i.e., ACC, BLCA, BRCA, CESC, KIRC, LAML, OV, PCPG, PRAD, TGCT, UCEC, and UCS) and upregulated in 11 (CHOL, COAD, DLBC, ESCA, GBM, LGG, LIHC, PAAD, READ, STAD, and THYM; Figure 7A). When we used a gene-expression heat map to examine the nine common DEGs across GEO and TARGET cohorts, we found that GREM1 and CAV1 expression were positively correlated with GNG12 expression in the TARGET cohort. Additionally, NDUFB9, AMBN, FAT3, and IBSP expression were negatively correlated (Figure 7B). We then generated six scatterplots of correlations between the expression of these six genes and of GNG12 (Figures 7C–H).
Figure 7 Expression analysis of GNG12-related genes. (A) Differential expression of GNG12 in pan-cancer. (B) A gene expression-related heat map of the nine common DEGs between GEO and TARGET sets. (C–H) Scatter plots of the expression correlation between GREM1, CAV1, NDUFB9, AMBN, FAT3, and IBSP and GNG12. Significance was set at p < 0.05 (ns, p < 0.05; *p < 0.05; **p < 0.01; ***p < 0.001).
3.6 Competitive Endogenous RNA Network
Using limma, we obtained 22 differentially expressed lncRNAs (DELs) from the TARGET database and visualized them on a heat map (Figure 8A). We then used the miRcode website to predict miRNAs that were highly conserved with DELs. Furthermore, we identified mRNAs targeted by the miRNAs through a series of databases and then constructed a ceRNA network with DEGs (Figure 8B).
Figure 8 GNG12-related ceRNA network construction and survival analysis. (A) The TARGET data set shows the significant difference lncRNAs between the two groups in the form of a heat map. (B) A ceRNA network: blue diamond, lncRNAs; green triangle, miRNAs; red round, mRNAs. (C) The effect of low GNG12 expression on the prognosis of osteosarcoma overall survival (OS) is statistically significant. (D) The effect of high miR-27a-3p expression on the prognosis of osteosarcoma overall survival (OS) is statistically significant. (E) The effect of low HOTTIP expression on the prognosis of osteosarcoma overall survival (OS) is statistically significant. (F) The effect of low GNG12 expression on the prognosis of osteosarcoma progress-free survival (PFS) is statistically significant. (G) Time-dependent ROC curves, 1 year (AUC = 0.692), 3 years (AUC = 0.658), and 5 years (AUC = 0.603). (H) The area under the ROC curve (AUC) in the GEO set assesses the performance for distinguishing osteosarcoma from normal tissue of GNG12 (AUC = 0.92).
3.7 Survival Analysis of GNG12 in TARGET Cohort
The results of our KM analysis indicated that OS and PFS differed significantly between low and high GNG-expression groups (Figures 8C, F). In addition, the ceRNA network revealed that low HOTTIP and GNG12 expression led to high risk. Elevated hsa-miR-27a-3p expression also resulted in high risk (Figures 8D, E). Through time-dependent ROC curves, we found that GNG12 was a reliable predictor of patient prognosis (Figure 8G).
We used an ROC curve to analyze whether GNG12 expression could distinguish between 15 normal samples and 103 tumor samples from the GEO dataset. We found an AUC of 0.920, indicating that GNG12 expression was an excellent biomarker for distinguishing osteosarcoma from normal tissue. (Figure 8H).
3.8 Validating Prognostic Value of GNG12 Expression in Osteosarcoma
We used 78 patients with osteosarcoma as a validation cohort to confirm the prognostic reliability of GNG12 (Table 1 and Supplementary Table S3). GNG12 expression levels (high vs. low) of these patients were confirmed using immuno-histochemistry (IHC) (Figures 9A, B). The results of KM analysis showed that patients with low GNG12 expression tend to have poor survival prognosis (Figure 9C). Analysis of time-dependent ROC curves revealed that 1-/3-/5-year AUCs were 0.961, 0.826, and 0.808 (Figure 9F). In addition, GNG12 expression was lower in metastatic than in non-metastatic osteosarcoma (Figure 9D); patients with metastatic osteosarcoma have a poor survival prognosis (Figure 9E).
Table 1 Baseline clinical characteristics of validation cohort.
Figure 9 Validation the GNG12 expression and prognostic value of GNG12. (A, B) Low/high H-score of GNG12 ICH images. (C) The effect of low GNG12 expression on the prognosis of osteosarcoma overall survival (OS) is statistically significant in the validation cohort. (D) In the validation cohort, the H-score level of GNG12 between metastasis and non-metastasis groups. (E) The effect of metastasis on the prognosis of osteosarcoma overall survival (OS) is statistically significant in the validation cohort. (F) Time-dependent ROC curves: 1 year (AUC = 0.961), 3 years (AUC = 0.826), and 5 years (AUC = 0.808). ***p < 0.001.
4 Discussion
The transducer and transmembrane-signal regulator GNG12 plays an important role in cancer progression through enhancing cell proliferation (15, 16). It is indispensable for guanosine triphosphatase (GTPase) activity, which functions in GTP-to-GDP catalysis (17). GTPase is associated with proliferation, invasion, and migration of tumors in multiple cancers, including osteosarcoma (18, 19). Despite this link, few studies have explored the role of GNG12 in osteosarcoma. Therefore, in this study, we hypothesized that GNG12 expression may influence osteosarcoma progression. Furthermore, GNG12 may have prognostic and therapeutic value for this cancer.
Analysis of data from GEO indicated that GNG12 expression was significantly downregulated in osteosarcoma tissues compared with normal samples. This reduced expression was positively correlated with poor prognosis in the TARGET cohort. Our functional analyses corroborated previous studies, showing that GNG12 was involved in ossification (20–22), protein targeting to ER (23, 24), extracellular-related terms (25, 26), focal adhesion (25, 27), PI3K-Akt signaling pathway (28, 29), GPCR ligand binding (30), and the matrisome. Through PPI networks, we identified the ER lumen as the location of hub gene clusters present in both GEO and TARGET cohorts. Regulating ER function plays a major role in the treatment of osteosarcoma (23, 31). Specifically, H2S-releasing doxorubicins caused misfolding of ER-related proteins to enhance ER-dependent apoptosis, which is effective against doxorubicin-resistant osteosarcoma in vitro (23). Additionally, lexibulin induces autophagy and apoptosis of osteosarcoma in vitro through triggering mutually enhanced reactive oxygen species and ER stress, suggesting that it could be an effective drug for osteosarcoma therapy (31). Taken together, our research and previous studies provided insight on GNG12’s role in cancer pathogenesis and demonstrated that the protein is a potential biomarker of osteosarcoma.
Integrated landscapes depicting the osteosarcoma microenvironment may help explain how the cancer responds to immunotherapy and help the development of new treatment strategies (32, 33). Our immune infiltration analysis of the GEO cohort showed that the high and low GNG12-expression groups differed significantly in macrophages (M0, M1, and M2), activated CD4+ memory T cells, mast cells, activated dendritic cells, and eosinophils. The same analysis on the TARGET cohort revealed that the three macrophage subtypes and mast cells were significantly different. Combining the two results indicated that macrophages and mast cells play a crucial role in the immune microenvironment of osteosarcoma. Moreover, osteosarcoma occurrence and development are probably related to inflammation and metabolic pathways. Thus, novel osteosarcoma therapy could potentially target macrophages and mast cells to improve their distribution in patients (34, 35). Findings from a mouse tail metastasis model indicated that M2 macrophages enhanced metastasis of osteosarcoma cells (K7M2 WT), thus identifying M2-polarized tumor-associated macrophages (TAMs) as a therapeutic target. The same study also found that all-trans retinoic acid (ATRA) inhibited M2 polarization of TAMs via downregulating MMP12 expression, thereby limiting osteosarcoma (36). Our study thus examined the correlation between GNG12 expression and MMP12. We found that MMP12 expression is negatively correlated with GNG12 expression, and high MMP12 expression can also predict poor osteosarcoma PFS (Supplementary Figure S1A, B). Therefore, GNG12 may inhibit M2 polarization of TAMs, which negatively regulates MMP12 expression and thus suppresses osteosarcoma metastasis.
Most studies have shown that ceRNA plays an important role in osteosarcoma occurrence and development (37, 38). We constructed a ceRNA network based on TARGET data and found that lncRNA HOTTIP (a low risk factor) can function as a ceRNA, decoying miR-27a-3p (a high risk factor) to promote GNG12-mediated metastasis. Targeting HOTTIP may therefore prove to be beneficial for osteosarcoma treatment (39). Liu et al. reported that miR-27a-3p were upregulated in osteosarcoma to promote the proliferation and invasion of osteosarcoma cells via inhibiting expression of TET1 (40). Although these results suggest that HOTTIP may regulate GNG12 through ceRNA and is involved in osteosarcoma progress, further experimental studies are needed to confirm these conclusions. Nevertheless, HOTTIP, miR-27a-3p, and GNG12 are all potential therapeutic targets for osteosarcoma.
This study has several limitations. First, the sample size of our healthy controls was much smaller than the sample size of patients with osteosarcoma. Future studies would benefit from balancing the number of participants in each group. Second, our outcomes were validated in a set of 78 patients with osteosarcoma, using an H-score to evaluate GNG12 expression in osteosarcoma tissue. Although the results confirmed that patients with low GNG12 expression had poor prognosis, this was a retrospective analysis. Thus, future research should employ a prospective methodology to avoid analysis bias. Finally, we did not verify our findings using in vitro or in vivo experiments, meaning that the exact mechanisms of GNG12 involvement in osteosarcoma remain unclear. This is an important topic for further investigation. Therefore, we highlight several areas in which further work is needed to deepen our understanding. First, as the response of GNG12, it would be interesting to examine the basic expression of these predicted DEGs and hub genes with Western blot, quantitative PCR (qPCR), IHC, immuno-fluorescence (IF) assays, and so on. Second, to clarify the function of GNG12, DEGs, and hub genes in osteosarcoma, a clean loss-of-function and gain-on-function study with tissue- and cell-type specificities remains warranted.
In conclusion, GNG12 mRNA and protein expression were downregulated in osteosarcoma, a change that was related to poor prognosis. We found that GNG12 may promote osteosarcoma through regulating ER function. HOTTIP/miR-27a-3p are candidates for regulating GNG12 expression to inhibit macrophage infiltration and suppress adaptive immunity. We provided important evidence supporting GNG12 as a biomarker for osteosarcoma prognosis and highlighting its potential as a target for immunotherapy.
Data Availability Statement
The datasets analyzed for this study can be found in the TARGET repository (https://target-data.nci.nih.gov/Public/OS/mRNA-seq/) and Gene Expression Omnibus (GEO, https://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42352/matrix/). In addition, the validation cohort presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics Statement
The studies involving human participants were reviewed and approved by Ethics Committee of The Second Affiliated Hospital of Nanchang University [Review (2020) No. (086)]. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.
Author Contributions
JY, ZY, and XC contributed to the conception and design of this study. JY and ZY collected the data sets from the database. JY, TW, JJ, JG, JZ, and TL performed the bioinformatics and statistical analysis. JY, ZY, and AY wrote the first draft of the manuscript. All authors contributed to the article and approved the submitted version.
Funding
The present study was supported by the National Natural Science Foundation of China (grant nos. 8166090137 and 8186090165), the National Natural Science Foundation of Jiangxi Province (No. 911084415069), and General Plan of Jiangxi Provincial Health Commission (No. 20195568).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.758845/full#supplementary-material
References
1. Ritter J, Bielack SS. Osteosarcoma. Ann Oncol (2010) 21 Suppl 7:vii320–325. doi: 10.1093/annonc/mdq276
2. Isakoff MS, Bielack SS, Meltzer P, Gorlick R. Osteosarcoma: Current Treatment and a Collaborative Pathway to Success. J Clin Oncol (2015) 33(27):3029–35. doi: 10.1200/JCO.2014.59.4895
3. Rickel K, Fang F, Tao J. Molecular Genetics of Osteosarcoma. Bone (2017) 102:69–79. doi: 10.1016/j.bone.2016.10.017
4. Bishop MW, Janeway KA, Gorlick R. Future Directions in the Treatment of Osteosarcoma. Curr Opin Pediatr (2016) 28(1):26–33. doi: 10.1097/MOP.0000000000000298
5. Yi X, Deng X, Zhao Y, Deng B, Deng J, Fan H, et al. Ubiquitin-Like Protein FAT10 Promotes Osteosarcoma Growth by Modifying the Ubiquitination and Degradation of YAP1. Exp Cell Res (2020) 387(2):111804. doi: 10.1016/j.yexcr.2019.111804
6. Feng W, Dean DC, Hornicek FJ, Spentzos D, Hoffman RM, Shi H, et al. Myc Is a Prognostic Biomarker and Potential Therapeutic Target in Osteosarcoma. Ther Adv Med Oncol (2020) 12:1758835920922055. doi: 10.1177/1758835920922055
7. Morishita R, Nakayama H, Isobe T, Matsuda T, Hashimoto Y, Okano T, et al. Primary Structure of a Gamma Subunit of G Protein, Gamma 12, and Its Phosphorylation by Protein Kinase C. J Biol Chem (1995) 270(49):29469–75. doi: 10.1074/jbc.270.49.29469
8. Larson KC, Lipko M, Dabrowski M, Draper MP. Gng12 Is a Novel Negative Regulator of LPS-Induced Inflammation in the Microglial Cell Line BV-2. Inflammation Res (2010) 59(1):15–22. doi: 10.1007/s00011-009-0062-2
9. Babadagli ME, Tezcan B, Yilmaz ST, Tufan AC. Matrilin-3 as a Putative Effector of C-Type Natriuretic Peptide Signaling During TGF-β Induced Chondrogenic Differentiation of Mesenchymal Stem Cells. Mol Biol Rep (2014) 41(9):5549–55. doi: 10.1007/s11033-014-3448-3
10. Orchel J, Witek L, Kimsa M, Strzalka-Mrozik B, Kimsa M, Olejek A, et al. Expression Patterns of Kinin-Dependent Genes in Endometrial Cancer. Int J Gynecol Cancer (2012) 22(6):937–44. doi: 10.1097/IGC.0b013e318259d8da
11. Zhang R, Deng Y, Zhang Y, Zhai GQ, He RQ, Hu XH, et al. Upregulation of HOXA13 as a Potential Tumorigenesis and Progression Promoter of LUSC Based on qRT-PCR and Bioinformatics. Int J Clin Exp Pathol (2017) 10(10):10650–65.
12. Kuijjer ML, Peterse EF, van den Akker BE, Briaire-de Bruijn IH, Serra M, Meza-Zepeda LA, et al. IR/IGF1R Signaling as Potential Target for Treatment of High-Grade Osteosarcoma. BMC Cancer (2013) 13:245. doi: 10.1186/1471-2407-13-245
13. Kuijjer ML, van den Akker BE, Hilhorst R, Mommersteeg M, Buddingh EP, Serra M, et al. Kinome and mRNA Expression Profiling of High-Grade Osteosarcoma Cell Lines Implies Akt Signaling as Possible Target for Therapy. BMC Med Genomics (2014) 7:4. doi: 10.1186/1755-8794-7-4
14. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007
15. Luo C, Zhao S, Dai W, Zheng N, Wang J. Proteomic Analyses Reveal GNG12 Regulates Cell Growth and Casein Synthesis by Activating the Leu-Mediated Mtorc1 Signaling Pathway. Biochim Biophys Acta Proteins Proteom (2018) 1866(11):1092–101. doi: 10.1016/j.bbapap.2018.08.013
16. Li J, Jin C, Zou C, Qiao X, Ma P, Hu D, et al. GNG12 Regulates PD-L1 Expression by Activating NF-κb Signaling in Pancreatic Ductal Adenocarcinoma. FEBS Open Bio (2020) 10(2):278–87. doi: 10.1002/2211-5463.12784
17. Cai YD, Zhang S, Zhang YH, Pan X, Feng K, Chen L, et al. Identification of the Gene Expression Rules That Define the Subtypes in Glioma. J Clin Med (2018) 7(10):350. doi: 10.3390/jcm7100350
18. Lopez-Guerrero AM, Espinosa-Bermejo N, Sanchez-Lopez I, Macartney T, Pascual-Caro C, Orantos-Aguilera Y, et al. RAC1-Dependent ORAI1 Translocation to the Leading Edge Supports Lamellipodia Formation and Directional Persistence. Sci Rep (2020) 10(1):6580. doi: 10.1038/s41598-020-63353-5
19. Marko TA, Shamsan GA, Edwards EN, Hazelton PE, Rathe SK, Cornax I, et al. Slit-Robo GTPase-Activating Protein 2 as a Metastasis Suppressor in Osteosarcoma. Sci Rep (2016) 6:39059. doi: 10.1038/srep39059
20. Eckardt JJ, Ivins JC, Perry HO, Unni KK. Osteosarcoma Arising in Heterotopic Ossification of Dermatomyositis: Case Report and Review of the Literature. Cancer (1981) 48(5):1256–61. doi: 10.1002/1097-0142(19810901)48:5<1256::AID-CNCR2820480534>3.0.CO;2-4
21. Iannaci G, Luise R, Sapere P, Costanzo RM, Rossiello R. Extraskeletal Osteosarcoma: A Very Rare Case Report of Primary Tumor of the Colon-Rectum and Review of the Literature. Pathol Res Pract (2013) 209(6):393–6. doi: 10.1016/j.prp.2013.03.010
22. Hoshi M, Takada J, Oebisu N, Nakamura H. Osteosarcoma Arising From the Parapharyngeal Space: A Case Report. Oncol Lett (2014) 8(3):1240–2. doi: 10.3892/ol.2014.2265
23. Buondonno I, Gazzano E, Tavanti E, Chegaev K, Kopecka J, Fanelli M, et al. Endoplasmic Reticulum-Targeting Doxorubicin: A New Tool Effective Against Doxorubicin-Resistant Osteosarcoma. Cell Mol Life Sci (2019) 76(3):609–25. doi: 10.1007/s00018-018-2967-9
24. Ward BK, Rea SL, Magno AL, Pedersen B, Brown SJ, Mullin S, et al. The Endoplasmic Reticulum-Associated Protein, OS-9, Behaves as a Lectin in Targeting the Immature Calcium-Sensing Receptor. J Cell Physiol (2018) 233(1):38–56. doi: 10.1002/jcp.25957
25. Chen Y, Pasapera AM, Koretsky AP, Waterman CM. Orientation-Specific Responses to Sustained Uniaxial Stretching in Focal Adhesion Growth and Turnover. Proc Natl Acad Sci USA (2013) 110(26):E2352–61. doi: 10.1073/pnas.1221637110
26. Cui J, Dean D, Hornicek FJ, Chen Z, Duan Z. The Role of Extracelluar Matrix in Osteosarcoma Progression and Metastasis. J Exp Clin Cancer Res (2020) 39(1):178. doi: 10.1186/s13046-020-01685-w
27. Hu C, Chen X, Wen J, Gong L, Liu Z, Wang J, et al. Antitumor Effect of Focal Adhesion Kinase Inhibitor PF562271 Against Human Osteosarcoma In Vitro and In Vivo. Cancer Sci (2017) 108(7):1347–56. doi: 10.1111/cas.13256
28. Wang R, Liu W, Wang Q, Li G, Wan B, Sun Y, et al. Anti-Osteosarcoma Effect of Hydroxyapatite Nanoparticles Both In Vitro and In Vivo by Downregulating the FAK/PI3K/Akt Signaling Pathway. Biomater Sci (2020) 8(16):4426–37. doi: 10.1039/D0BM00898B
29. Ma Z, Yang J, Yang Y, Wang X, Chen G, Shi A, et al. Rosmarinic Acid Exerts an Anticancer Effect on Osteosarcoma Cells by Inhibiting DJ-1 via Regulation of the PTEN-PI3K-Akt Signaling Pathway. Phytomedicine (2020) 68:153186. doi: 10.1016/j.phymed.2020.153186
30. Luo J, Sun P, Siwko S, Liu M, Xiao J. The Role of GPCRs in Bone Diseases and Dysfunctions. Bone Res (2019) 7:19. doi: 10.1038/s41413-019-0059-6
31. Wang Z, Yin F, Xu J, Zhang T, Wang G, Mao M, et al. CYT997(Lexibulin) Induces Apoptosis and Autophagy Through the Activation of Mutually Reinforced ER Stress and ROS in Osteosarcoma. J Exp Clin Cancer Res (2019) 38(1):44. doi: 10.1186/s13046-019-1047-9
32. Zhou Y, Yang D, Yang Q, Lv X, Huang W, Zhou Z, et al. Single-Cell RNA Landscape of Intratumoral Heterogeneity and Immunosuppressive Microenvironment in Advanced Osteosarcoma. Nat Commun (2020) 11(1):6322. doi: 10.1038/s41467-020-20059-6
33. Hong W, Yuan H, Gu Y, Liu M, Ji Y, Huang Z, et al. Immune-Related Prognosis Biomarkers Associated With Osteosarcoma Microenvironment. Cancer Cell Int (2020) 20:83. doi: 10.1186/s12935-020-1165-7
34. Zhang C, Zheng JH, Lin ZH, Lv HY, Ye ZM, Chen YP, et al. Profiles of Immune Cell Infiltration and Immune-Related Genes in the Tumor Microenvironment of Osteosarcoma. Aging (Albany NY) (2020) 12(4):3486–501. doi: 10.18632/aging.102824
35. Heymann MF, Lézot F, Heymann D. The Contribution of Immune Infiltrates and the Local Microenvironment in the Pathogenesis of Osteosarcoma. Cell Immunol (2019) 343:103711. doi: 10.1016/j.cellimm.2017.10.011
36. Zhou Q, Xian M, Xiang S, Xiang D, Shao X, Wang J, et al. All-Trans Retinoic Acid Prevents Osteosarcoma Metastasis by Inhibiting M2 Polarization of Tumor-Associated Macrophages. Cancer Immunol Res (2017) 5(7):547–59. doi: 10.1158/2326-6066.CIR-16-0259
37. Wang Y, Zeng X, Wang N, Zhao W, Zhang X, Teng S, et al. Long Noncoding RNA DANCR, Working as a Competitive Endogenous RNA, Promotes ROCK1-Mediated Proliferation and Metastasis via Decoying of miR-335-5p and miR-1972 in Osteosarcoma. Mol Cancer (2018) 17(1):89. doi: 10.1186/s12943-018-0837-6
38. Zhou FC, Zhang YH, Liu HT, Song J, Shao J. LncRNA LINC00588 Suppresses the Progression of Osteosarcoma by Acting as a ceRNA for miRNA-1972. Front Pharmacol (2020) 11:255. doi: 10.3389/fphar.2020.00255
39. Liu K, Ni JD, Li WZ, Pan BQ, Yang YT, Xia Q, et al. The Sp1/FOXC1/HOTTIP/LATS2/ YAP/β-Catenin Cascade Promotes Malignant and Metastatic Progression of Osteosarcoma. Mol Oncol (2020) 14(10):2678–95. doi: 10.1002/1878-0261.12760
Keywords: osteosarcoma, GNG12, biomarker, therapeutic target, prognosis
Citation: Yuan J, Yuan Z, Ye A, Wu T, Jia J, Guo J, Zhang J, Li T and Cheng X (2021) Low GNG12 Expression Predicts Adverse Outcomes: A Potential Therapeutic Target for Osteosarcoma. Front. Immunol. 12:758845. doi: 10.3389/fimmu.2021.758845
Received: 15 August 2021; Accepted: 17 September 2021;
Published: 06 October 2021.
Edited by:
Nicholas Vitanza, Seattle Children’s Hospital, United StatesReviewed by:
Ashleigh Fordham, Children’s Cancer Institute Australia, AustraliaXiaoling Zhang, Zhejiang University, China
Copyright © 2021 Yuan, Yuan, Ye, Wu, Jia, Guo, Zhang, Li and Cheng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xigao Cheng, 228206846@qq.com
†These authors have contributed equally to this work