The Prognostic and Immunological Value of Guanylate-Binding Proteins in Lower-Grade Glioma: Potential Markers or Not?

Seven guanylate-binding proteins (GBPs, GBP1–7), identified as a subfamily of interferon-γ-induced guanosine triphosphate hydrolases (GTPases), has been reported to be closely associated with tumor progression, metastasis, and prognosis of cancer patients in recent years. However, the expression patterns, prognostic value, immune infiltration relevance, and biological functions of GBPs in lower-grade glioma (LGG) remain elusive. In this study, by analysis and verification through multiple public data platforms, we found that GBP1, 2, 3, 4 were significantly upregulated in LGG tissues vs normal brain tissue. Analysis based on the Cox proportional hazard ratio and Kaplan–Meier plots demonstrated that the high expressions of GBP 1, 2, 3, 4 were significantly correlated with the poor prognosis of LGG patients. Correlation analysis of clinical parameters of LGG patients indicated that the expressions of GBP 1, 2, 3, 4 were significantly associated with the histological subtype and tumor histological grade of LGG. Furthermore, the correlation analysis of immune infiltration showed that the expressions of GBP1, 2, 3, 4 were significantly and positively correlated with the level of tumor immune-infiltrating cells. In particular, GBP1, 2, 3, 4 expressions were strongly correlated with the infiltration levels of monocyte, TAM, and M1/M2 macrophage, revealing their potential to regulate the polarity of macrophages. Finally, we used the GSEA method to explore the signaling pathways potentially regulated by GBP1, 2, 3, 4 and found that they were all closely associated with immune-related signaling pathways. Collectively, these findings suggested that GBP1, 2, 3, 4 were potent biomarkers to determine the prognosis and immune cell infiltration of LGG patients.


INTRODUCTION
Glioma is derived from astrocytes and/or oligodendrocytes and is one of the most common primary central nervous system tumors (Jang and Kim, 2018). Lower-grade glioma (LGG), the crucial pathological type of glioma, comprises grade II and grade III gliomas defined by the World Health Organization (WHO), mainly including anaplastic astrocytomas, oligodendrogliomas, and oligoastrocytomas (Brat et al., 2015).
LGG has the characteristics of diffuse infiltration, metastasis, and easy progression to higher-grade gliomas, which seriously affects human survival, especially young adults who enjoy an active life (Mazurowski et al., 2017;Liu et al., 2019). In recent years, comprehensive treatments such as postoperative chemotherapy, radiotherapy, and immunotherapy have made great progress, but the survival rate of LGG patients is still unsatisfactory and unpredictable. Thus, the identification of novel prognostic biomarkers or molecular targets is imperative for a highly accurate prediction of survival or guidance for individualized treatment of LGG patients.
Guanylate-binding protein (GBP) is classified as a unique subfamily of interferon-γ-induced guanosine triphosphate hydrolases (GTPases), which can hydrolyze guanosine triphosphate (GTP) to both guanosine diphosphate (GDP) and guanosine monophosphate (GMP) (Ghosh et al., 2006). In humans, seven GBP proteins (GBP1-7) with molecular weights in the range of 67-73 kDa have been well identified (Tripal et al., 2007). Studies have shown that GBPs, such as GBP1 and GBP2, are closely related to host defense against pathogens, and have antiviral and antibacterial activities in the process of host anti-infection and antiinflammatory defense (Vestal and Jeyaratnam, 2011;Honkala et al., 2019). However, the roles of GBPs in cancer are diverse and complicated. GBP1 upregulation is reported to be associated with decreased disease progression and better overall survival in patients with breast and colorectal cancer (Naschberger et al., 2008;Lipnik et al., 2010), while it is connected with increased disease progression, metastasis, and treatment resistance in ovarian cancer and glioblastoma (Duan et al., 2006;De Donato et al., 2012;Ji et al., 2019). GBP2 can enhance the invasion of glioblastoma (Yu et al., 2020), but inhibit the invasion ability of breast cancer cells . Thus, the functions of different GBPs in multiple cancers need to be further clarified.
The immune microenvironment has been determined to play a vital role in tumor biology . Tumor-infiltrating immune cells, including T and B lymphocytes, macrophages, neutrophils, dendritic cells, etc., are very important elements of the tumor microenvironment, which directly or indirectly regulate the growth and development of tumor cells and further affect the prognosis of many cancer patients including LGG (Domingues et al., 2016;. Recently, many promising preclinical and clinical immunotherapies have been implemented in malignant glioma, including immune checkpoint inhibitors, active or passive immunotherapy, etc. (Simonelli et al., 2018;Vismara et al., 2019), indicating that the immune components in the tumor microenvironment are of great value as prognostic biomarkers or therapeutic targets in glioma. Therefore, further exploration of immune infiltration regulation in the tumor microenvironment may support the treatment of cancers.
However, there are relatively a few studies on GBPs in LGG, and the prognostic value, the regulation of immune infiltration, and biological functions of GBPs in LGG need to be further clarified. In this study, we used public databases and online platforms and conducted a comprehensive and detailed analysis of the expression patterns, prognostic value, immune infiltration regulation, and biological functions of GBPs in LGG.

Data Collection
RNAseq data and corresponding clinical data of 509 LGG tissue samples were downloaded from TCGA (The Cancer Genome Atlas, https://portal.gdc.cancer.gov/). RNA array dataset (GSE4290) (Sun et al., 2006) was downloaded from the NCBI/GEO database (https:// www.ncbi.nlm.nih.gov/gds/). In this study, those data were used to perform gene expression analysis, clinical correlation analysis, and gene set enrichment analysis (GSEA) in LGG.

Oncomine Database Analysis
As the largest oncogene database and integrated data-mining platform in the world, Oncomine (http://www.oncomine.org) is used to compare transcriptome data between tumors and corresponding normal tissues in different types of cancer (Rhodes et al., 2004). In this study, relevant data were obtained to evaluate the expression of GBP family genes in LGG. The p-value cutoff was 0.05. Statistical differences were determined by Student's t-test.

GEPIA Database Analysis
GEPIA (http://gepia.cancer-pku.cn/) is an interactive web application that analyzes RNA sequencing expression data for more than 9,000 tumors and 8,000 normal samples from The Cancer Genome Atlas (TCGA) and GTEx projects (Tang et al., 2017). In this study, we performed gene expression analysis and prognostic analysis of GBP genes both in pan-cancer and in LGG with GEPIA. Besides, gene correlation analysis was also evaluated with the Spearman correlation coefficient by GEPIA. The p-value cutoff was 0.05. Student's t test was used to generate a p-value for expression, and a Kaplan-Meier curve was used for prognostic analysis.
Tumor Immune Estimation Resource Database Analysis TIMER (Tumor Immune Estimation Resource, https://cistrome. shinyapps.io/timer/) is a database designed for systematic analysis of immune cell infiltrates across diverse cancer types . In our study, we evaluated the correlation between GBP gene levels and the infiltration of immune cells as well as the correlation among GBP gene expressions and marker gene expressions of the infiltration of immune cells and clinical outcome. Specifically, we analyzed the correlation between differentially expressed GBPs and macrophage polarity through the "correlation" module in the TIMER database. The results of Univariate Cox survival analysis in "Survival" module is shown in Figure 5E, and the results of Multivariate Cox survival analysis is shown in Table 3. The correlation map of differentially expressed GBPs and macrophage-related marker genes is shown in Figure 6. Spearman correlation coefficient was chosen for the correlation analysis.

cBioPortal Database Analysis
CBioportal (http://www.cbioportal.org/) is an open platform for visualization, analysis, and download of multidimensional cancer genomics data (Gao et al., 2013). Based on the TCGA/LGG Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 651348 dataset (TCGA, provisional), we analyzed the genetic alterations and prognostic analysis of GBP genes in LGG.

GSCALite
GSCALite is a user-friendly web server for dynamic analysis and visualization of gene sets in 32 cancer types from TCGA (Liu et al., 2018). In this study, GSCALite was used to analyze the miRNA regulatory network of GBP genes in LGG using the "TCGA KIRC" dataset.

Search Tool for the Retrieval of Interacting Genes Database Analysis
The STRING (Search Tool for the Retrieval of Interacting Genes, https://string-db.org/) database aims to collect, score, and integrate both experimental as well as predicted protein-protein interaction (PPI) information and further achieve a comprehensive and objective global network, including direct (physical) as well as indirect (functional) interactions (Szklarczyk et al., 2019). In this study, we conducted a PPI network analysis of each GBP gene to explore the interactions of GBP genes.

GeneMANIA Database Analysis
GeneMANIA (http://www.genemania.org) provides information for protein and genetic interactions, pathways, co-expression, colocalization, and protein domain similarity of submitted genes and helps researchers predict the functions behind gene sets by constructing a protein-protein interaction (PPI) network (Warde-Farley et al., 2010).

Gene Set Enrichment Analysis
The gene set enrichment analysis (GSEA) was performed to identify significantly enriched groups of genes (Subramanian et al., 2005). In this study, the GSEA v4.0.3 software was applied to analyze biological pathway divergences between

Statistical Methods
In this study, SPSS 20.0 and GraphPad Prism 6.0 software were used for statistical analysis. The differential expression levels of GBPs were compared and analyzed by the Students't-test. Survival curves were generated from Kaplan-Meier Plotter, and their differences are analyzed using log-rank test in GEPIA. Chi-square test was used to determine the correlation between expressions of GBPs and clinical parameters. The correlation between expressions of GBPs and immune infiltration level or other marker genes in LGG were evaluated by Spearman's correlation and statistical significance. For GSEA, p < 0.05 and FDR (false discovery rate) q < 0.05 were considered as threshold values to estimate statistical significance.

Guanylate-Binding Proteins 1/2/3/4 Were Upregulated in Lower-Grade Glioma Patients
To explore the expression patterns of GBP family genes in brain and nervous system tumors, especially LGG, we first conducted a comprehensive analysis of the expression patterns of different GBP members using the Oncomine database. As shown in Figure 1, GBP1, 2, 3, 4, 5 were upregulated in the brain and nervous system tumors vs normal brain tissue, which was confirmed by data from 13, 7, 5, 1, and 1 datasets, respectively. Among them, the analysis results from the four datasets ( Table 1) simultaneously confirmed that the expression of GBP1 and GBP2 was significantly increased in LGG vs normal brain tissue, and the results from two datasets and one dataset, respectively, demonstrated that the expression of GBP3 and GBP4 significantly increased in LGG. To further determine the expression differences of these four GBP genes, we selected the GEPIA database and one GEO dataset for verification. As shown in Figures 2A-H, the data from both GEPIA database and GSE4290 dataset demonstrated that GBP1, 2, 3, 4 were significantly upregulated in LGG.

Prognostic Value of Guanylate-Binding Proteins 1/2/3/4 in Lower-Grade Glioma Patients
We continued to explore the prognostic value of GBP1, 2, 3, 4 in LGG patients using GEPIA database by evaluating the effect of gene expression on overall survival and disease-free survival of tumor patients. The survival significance maps ( Figures 3A,B) of pan-cancer based on the Cox proportional hazard ratio (HR) showed that GBP1, 2, 3, 4 had better prognostic value in LGG vs other tumor types, and highly expressed GBP1, 2, 3, 4 were all significantly unfavorable for both overall survival and disease-free survival of LGG patients. The Kaplan-Meier plots further demonstrated that LGG patients with highly expressed GBP1, 2, 3, 4 had shorter overall survival and disease-free survival time (Figures 3C-J).

Correlations of Guanylate-Binding Proteins 1/2/3/4 With Clinicopathological characteristics in Lower-Grade Glioma
Next, we analyzed the correlations between GBP1, 2, 3, 4 expressions and clinicopathological characteristics. The expression data and clinical data of 509 LGG patients were extracted from the TCGA database, and clinical parameters mainly include age, sex, histological subtype, and tumor histological grade. Chi-square test was used to determine the correlation between GBP1, 2, 3, 4 expressions and clinical parameters. As shown in Table 2, the expressions of GBP1, 2, 3, 4 were significantly correlated with the histological subtype and histological grade of LGG patients, and only the expression of GBP4 had a correlation with the age and gender of the patient. Furthermore, we explored the expression differences of GBP1, 2, 3, 4 in different histological subtypes and histological grades of LGG patients. We found that the expressions of GBP1, 2, 3, 4 were significantly increased in astrocytoma vs oligodendroglioma, and GBP2, 3 were highly expressed in oligoastrocytoma vs oligodendroglioma ( Figures 4A-D).
Besides, we found that the expressions of GBP 1, 2, 3, 4 were all significantly increased in poor histological grade of LGG patients ( Figures 4E-H), which was consistent with GBP 1, 2, 3, 4, which might be unfavorable factors for LGG patients.

Guanylate-Binding Protein P1/2/3/4 Expressions Were correlated With Immune cell Infiltration Levels in Lower-Grade Glioma Patients
Immune cell infiltration in the tumor microenvironment is an important factor affecting tumor progression and prognosis of cancer patients (Domingues et al., 2016). In order to explore whether GBP1, 2, 3, 4 regulates the level of infiltrating immune cells in the LGG microenvironment, we analyzed the correlation of the expressions of GBP 1, 2, 3, 4 with immune infiltrating cells based on the TIME database ( Figures 5A-D). Interestingly, we found that the expressions of GBP1, 2, 3, 4 were strongly and positively correlated with these immune-infiltrating cells, including B cells, CD8 + T cells, CD4 + T cells, macrophages, neutrophils, and dendritic cells. Furthermore, we explored the effects of six immune cells and GBP 1, 2, 3, 4 expressions on the prognosis of LGG patients. Univariate Cox survival analysis showed that the high infiltration levels of six types of immune cells and the high expressions of GBP 1, 2, 3, 4 indicated poor prognosis of LGG patients ( Figures 5E). Multivariate Cox survival analysis showed that macrophages, GBP1, and GBP2 were independent prognostic indicators for LGG patients ( Table 3). These findings indicated that GBP1, 2, 3, 4 may potentially regulate the level of immune cell infiltration in LGG, and a high level of immune cell infiltration is not conducive to patient survival.

Correlation Analysis Between Guanylate-Binding Protein 1/2/3/4 and Immune Markers in Lower-Grade Glioma Patients
To further clarify the relationship of GBP 1, 2, 3, 4 with immune infiltration, we analyzed the correlation between the expressions of GBP 1, 2, 3, 4 and gene markers of a variety of immune cells (Table 4). We found that GBP1 had a significant correlation with most of the gene markers of infiltrating immune cells, excluding one gene marker (STAT4) of T-helper 1 (Th1) cell, two gene markers (FOXP3 and STAT5B) of regulatory T cell (Treg), and two gene markers (KIR2DL1 and KIR3DL3) of natural killer cell. GBP2 had a strong correlation with almost all markers of infiltrating immune cells, except for the two gene markers (KIR2DL1 and KIR3DL3) of natural killer cell. GBP3 also showed a significant correlation with most markers of infiltrating immune cells, excluding two markers (FOXP3 and STAT5B) of Treg and four markers (KIR2DL1, KIR3DL1,  (Table 4 and Figure 6), which also had been verified in the GEPIA database (Table 5). This indicated that GBP1, 2, 3, 4 may be involved in regulating macrophage polarity in LGG.
Gene Alterations, co-expression, Interaction Network Analysis of Guanylate-Binding Protein 1/2/3/4 in Lower-Grade Glioma Then we focused on the gene alterations of GBP1, 2, 3, 4 in LGG using the cBioPortal platform. A total of 518 LGG patients were selected for this analysis. The genetic alteration frequency of GBP1, 2, 3, 4 in LGG, including amplification, high mRNA, deep deletion, and mutation,   was 3.02, 1.89, 4.15, 3.58%, respectively ( Figure 7A). The total genetic alteration frequency of GBP1, 2, 3, 4 was 5.66%, and high mRNA was the most common type of gene alteration in these samples ( Figure 7B). We further evaluated the impact of genetic alteration of GBP 1, 2, 3, 4 on patient survival and found that LGG patients with the genetic alteration of GBP1, 2, 3, 4 have shorter overall survival and disease-free survival time ( Figures 7C,D). We also found there was an inframe mutation in GBP2, and a missense mutation in GBP3 ( Figure 7E). We continue to explore the co-expression and interaction network of GBP 1, 2, 3, 4. We found that there was a strong expression correlation among GBP1, 2, 3 and 4 ( Figure 7F), and the functions of GBP1, 2, 3 and 4 were potentially regulated by different miRNAs ( Figure 7G). Also, there was a close interaction relationship between GBP1, 2, and 3 ( Figure  7H), and the functions of GBP1, 2, 3, 4 are mainly related to cellular response to type I interferon, interferon-gammamediated signaling pathway, chemokine activity, positive regulation of cAMP-mediated signaling, etc ( Figure 7I).

Pathway Enrichment Analysis of Guanylate-Binding Protein 1/2/3/4 in Lower-Grade Glioma
GSEA is used to explore the signaling pathways that are potentially regulated by GBP1, 2, 3, 4 in LGG. We divided the samples into the high-expression group and the lowexpression group based on the mean value. Pathways with higher frequency enriched in phenotype high of GBP1, 2, 3, 4 are presented in Figures 8A-D. We found that functions of GBP1, 2, 3, 4 were closely linked: 1) They were all involved in regulating immune-related signaling pathways, such as intestinal immune network for IgA production, primary immunodeficiency, B/T cell receptor signaling pathway, natural killer cell-mediated cytotoxicity, etc. 2) They were all closely related to cancer and participated in the regulation of cancer-related signaling pathways, such as JAK-STAT signaling pathway, apoptosis, etc. 3) They also regulated Toll-like receptor pathway signaling, NOD-like receptor signaling pathway, and chemokine signaling pathway. Together, these results indicated that GBP1, 2, 3, 4 had the potential to become therapeutic targets in LGG.

DISCUSSION
Previous studies have identified a family of IFN-inducible GTPases, namely, guanylate-binding proteins (GBPs), as a major nexus of IFN-driven complex homeostatic defense networks, which function in host defense to viral, bacterial, and protozoan pathogens (Tretina et al., 2019). In recent years, an increasing number of studies have also confirmed that GBPs are not only involved in regulating host immune defense but also closely related to tumor development and metastasis (Mustafa et al., 2018;Zhao et al., 2019;Yu et al., 2020), and some GBPs, such as GBP1 and GBP2, have shown good prognostic value in certain tumors, for example, breast, oral and colorectal cancer (Yu et al., 2011;Godoy et al., 2014;Wang et al., 2016). Elevated GBP1 expression has also been shown to be associated with chemotherapy resistance in lung, breast, and ovarian cancer (Duan et al., 2006;Fekete and Győrffy, 2019;Cheng et al., 2020). However, the biological function and prognostic value of individual GBP in LGG remain elusive. By analysis and verification through multiple public data platforms, we found that GBP1, 2, 3, 4 were significantly upregulated in LGG tissues vs normal brain tissue. Consistently, we further found that highly expressed GBP1, 2, 3, 4 were all significantly unfavorable for both overall survival and disease-free survival of LGG patients, suggesting the potential of GBP1, 2, 3, 4 as prognostic markers in LGG. To further explore the clinical significance of GBPs, we analyzed the correlation between their expressions and the clinical parameters of LGG patients, and we found that the expressions of GBP 1, 2, 3, and 4 were significantly associated with tumor histological grade of LGG. As the tumor histological grade increased, the expressions of GBP 1, 2, 3, and 4 significantly increased. This was consistent with GBP 1, 2, 3, 4, which might be unfavorable factors for LGG patients. Besides, the expressions of GBP1, 2, 3, 4 were significantly correlated with the histological subtype of LGG, and the expression of GBP4 also had a correlation with the age and gender of the patient.
Infiltrating immune cells in the tumor microenvironment, mainly including tumor-infiltrating lymphocytes (B cells, CD8 + T cells, and CD4 + T cells) and other immune cells (macrophages, neutrophils, and dendritic cells), have become the focus of current tumor research. Studies have shown that immuneinfiltrating cells play an indispensable function in the tumor microenvironment as a double-edged sword to promote or inhibit tumor cell progression (Fridman et al., 2012). On the one hand, immune infiltrating cells play an anti-tumor effect by monitoring and destroying cancer cells (Morvan and Lanier, 2016). On the other hand, studies have shown that cancer cells can evade the surveillance of immune-infiltrating cells through a variety of mechanisms or further manipulate these immune infiltrating cells to create a microenvironment that promotes tumor progression (Mantovani et al., 2008). The dual effect of immune-infiltrating cells on tumor cells may depend on the type of immune cells, the state of immune cells, and the microenvironment of different tumors (Domingues et al., 2016). In our study, we found that the expressions of GBP1, 2, 3, 4 were significantly and positively   Multivariate analysis showed that GBP1, GBP2, and macrophage infiltration are independent prognostic factors for LGG patients. The correlation between the expression of GBP 1, 2, 3, 4 and the expressions of other immune cell marker genes were assessed, further confirming the close connection between GBP 1, 2, 3, 4 and tumor immune-infiltrating cells. In particular, the expression of GBP 1, 2, 3, 4 are significantly and positively correlated with the marker genes of monocytes, TAM, M1, and M2 macrophages, indicating that GBP 1, 2, 3, 4 may be involved in regulating the polarity of macrophages. However, how GBP 1, 2, 3, 4 participate in the regulation of tumor immune-infiltrating cells requires further research to clarify. Then, we focused on the genetic alterations of GBP1, 2, 3 and 4 in LGG, and we found that mRNA high is the most common type of genetic alterations in all LGG patient samples with genetic alterations. The genetic alterations of GBP1, 2, 3 and 4 indicated a poor prognosis of LGG patients. We also found that there was an inframe mutation in GBP2, and a missense mutation in GBP3. Co-expression and interaction network analysis further revealed the close functional connection among them.
Finally, we used the GSEA method to explore the signaling pathways that may be potentially regulated by GBP1, 2, 3, and 4 in LGG. We found that GBP1, 2, 3, 4 are closely related to immune-related signaling pathways, such as intestinal immune network for IgA production, primary immunodeficiency, B/T cell receptor signaling pathway, natural killer cellmediated cytotoxicity, etc., which was consistent with the association between GBP1, 2, 3, 4 and immune cell infiltration that we explored above. Besides, we found that they also potentially regulated Toll-like receptor pathway signaling, NOD-like receptor signaling pathway, and chemokine signaling pathway. Together, these results indicated that GBP1, 2, 3, 4 had the potential to become therapeutic targets in LGG.
Our study still has some limitations. The analysis of gene transcription levels based on public data platforms cannot fully reflect the changes in protein levels. Therefore, experiments in vivo and in vitro are needed to verify our findings and further promote the understanding of GBPs in LGG. Despite these limitations, our study may help guide further investigation of GBPs in LGG.
In conclusion, we systematically and comprehensively analyzed the expression pattern, prognostic value, correlation with clinical parameters, immune infiltration relevance of GBPs in LGG, and further explored their potential regulatory signaling pathways. Our results indicated that GBP1, 2, 3, 4 were potential biomarkers that can be used to predict prognosis and tumor immune infiltration of LGG patients. We hope that our results can help clinicians better predict the survival of LGG patients or choose appropriate treatment methods or therapeutic drugs, thereby improving the survival prognosis of cancer patients.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.