Prognostic and Immune Infiltration Signatures of GIMAP Family Genes in Clear Cell Renal Cell Carcinoma

Background: Clear cell renal cell carcinoma (ccRCC) is a common malignant tumor of the urinary system characterized by abundant immunocytes infiltration. The impact of guanosine triphosphatases (GTPases) of immunity-associated proteins (GIMAPs) on the tumor immune microenvironment (TIME) and prognosis of ccRCC is unclear. Methods: The expression of GIMAPs in ccRCC was determined through multiple datasets (ONCOMINE, TCGA and UALCAN). The relationship between GIMAP family members was analyzed through Spearman correlation analysis. The interaction among the GIMAPs protein was analyzed using STRING. Prognostic values of GIMAPs were evaluated by Survival analysis, Lasso and Cox regression analysis; Prognostic risk model and nomogram were constructed. The correlation between GIMAPs and TIME was explored using TIMER, Cibersort and Pearson correlation analysis. Gene set enrichment analysis (GSEA) was performed to discuss their function and mechanism in ccRCC. Results: GIMAPs were over-expressed in ccRCC and significantly related to overall survival (OS) of the patients. GIMAPs were positively correlated with each other, the risk model based on GIMAPs had good prognostic value in ccRCC. GIMAPs mainly expressed in TIME and were associated with abundant immunocytic infiltration in ccRCC, the risk model also had close correlation with TIME. Our results showed GIMAPs may affect the development of ccRCC by regulating the amount and antitumor activity of immunocytes in TIME. Conclusions: GIMAPs were over-expressed in ccRCC, and their expression levels were significantly related to the OS of patients and immunocytic infiltration in TIME. GIMAPs are potential therapeutic targets and prognostic biomarkers for ccRCC.

Keywords

GIMAPs

clear cell renal cell carcinoma

tumor immune microenvironment

prognosis

biomarker

1. Introduction

Renal cell carcinoma is an aggressive urinary system malignancy, endangering human health. In 2020, over 4,000,000 new cases and about 1,800,000 deaths had been reported globally [1]. Clear cell renal cell carcinoma (ccRCC) is the most common histological subtype of renal cell carcinoma, accounted for approximately eighty percent of all cases [2]. Patients with early ccRCC can achieve favorable outcome by radical surgery, but patients with advanced ccRCC and postoperative metastasis have poor prognosis [3, 4]. CcRCC is insensitive to chemotherapeutic drugs and radiotherapy. Interleukin-2 and $\alpha{}$ -interferon have some benefits in the treatment of advanced ccRCC, but they only work for a fraction of patients and sometimes can result in severe side effects [5]. In recent years, anti-angiogenic tyrosine kinase inhibitors and immune checkpoint inhibitors have showed effectiveness in the treatment of advance-stage ccRCC patients, however some patients do not benefit from this therapy [6, 7, 8]. In The Cancer Genome Atlas (TCGA), ccRCC had higher infiltration scores of T cells, Th1, cytotoxic cells, dendritic cells (DCs) and neutrophils, and lower scores of Th2 and regulatory T cells (Tregs), in comparison with other 18 epithelial cancers [9]. It has been shown that tumor immune microenvironment (TIME) correlate with the effect of immunotherapy and clinical outcomes [10, 11, 12]. TIME rich in Th2 and Tregs was associated with tumor mutation burden and inhibited immune response [13]. The decrease of tumor-associated macrophages (TAMs) can enhance CD8 ${}^{+}$ T cells infiltration and their capacity to migrate in to tumor cells [14]. Therefore, finding the key factors that affecting TIME for predicting the efficacy of immunotherapy and establishing new targets in ccRCC is of great importance.

Human guanosine triphosphatases of immunity-associated proteins (GIMAP) family genes are located on chromosome 7, spanning about 500 KB, and include seven functional genes (GIMAP1, GIMAP2, GIMAP4, GIMAP5, GIMAP6, GIMAP7, GIMAP8) and a pseudogene [15]. The GIMAP proteins are similar in N-end sequence and contain guanine nucleotide binding domain called GTPase [15, 16]. Most of these proteins participate in the maintenance and development of lymphocytes. The deficiency of GIMAP5 leads the decrease of peripheral T, B cells and natural killer (NK) cells in mice [17, 18]. GIMAP1 is important for the maintenance of T cells’ proliferation and mature the function of B cells [19, 20]. GIMAP4 may promote T cell apoptosis [21]. Knocking out GIMAP6 makes Jurkat T cells more susceptible to apoptosis inducers [22]. Recent studies have found the dysregulation of GIMAPs in a number of tumor types, including hepatocellular cancer, endometrial cancer and non-small cell lung cancer, GIMAPs were significantly associated with the prognosis of patients with lung adenocarcinoma and endometrial cancer [23, 24, 25, 26]. CcRCC is characterized by the infiltration of abundant immunocytes, nonetheless, the influence of the GIMAP family on the immunological microenvironment of ccRCC and the prognosis of ccRCC patients is unclear.

In this study, we sought to investigate the impact of GIMAPs on the immunological microenvironment and prognosis of ccRCC. We identified GIMAPs’ expression in ccRCC from multiple databases. Then, we discussed the impact of GIMAPs on outcomes in ccRCC patients. We then elucidated the interrelationship between GIMAPs and the immunocytes in ccRCC and discussed the molecular mechanism of GIMAPs affecting the development of ccRCC. Our results showed that GIMAPs could regulate TIME, affect the prognosis of these patients and are potential therapeutic targets for ccRCC.

2. Materials and Methods

2.1 Identification of GIMAPs Expression in ccRCC

We contrasted the expression of GIMAPs mRNA in normal renal tissues and ccRCC using ONCOMINE (https://www.oncomine.org, accessed on 14 October 2021), which can provide powerful and reliable function to analyze multiple expression of gene characteristics of tumors in the Gene Expression Omnibus, TCGA and published literature [27]. Then, we downloaded RNA-seq and clinical data of patients with ccRCC from TCGA (https://portal.gdc.cancer.gov/, accessed on 16 December 2021) and UCSC Xexa (https://xenabrowser.net/datapages/, accessed on 31 December 2021). We were able to analyze clinical data for 529 patients and RNA sequencing data for 530 tumor samples and 72 normal samples, the expression of genes are normalized as log2 (TPM + 1). The expression of GIMAPs mRNA between them were compared to verify the results in ONCOMINE. GIMAPs’ effect was evaluated by receiver operating characteristic (ROC) curves and multivariate logistic regression analysis for the diagnosis of ccRCC. The protein expression of GIMAPs in normal renal tissues and ccRCC were compared by UALCAN (http://ualcan.path.uab.edu, accessed on 21 October 2021). UALCAN is an online tool that can analyze tumor patients’ transcriptome data and clinical parameter in TCGA, also the expression of proteins of tumors in the Clinical Proteomic Tumor Analysis Consortium (CPTAC) [28, 29].

2.2 Protein-Protein Interaction

STRING (version: 11.5, https://cn.string-db.org/, accessed on 3 November 2021) was utilized to evaluate the interreaction among the GIMAPs protein and construct the network diagram. STRING is a database containing many types of protein-protein interactions [30].

2.3 Survival Analysis, Lasso and Cox Regression Analysis

Patients’ inclusion criteria: the patients with RNA-seq data and complete survival data (OS and survival status). 527 patients were included. On the basis of the median expression of GIMAPs, the patients were divided into Low and High groups. The effect of GIMAPs on OS of patients with ccRCC in TCGA was assessed using the Kaplan Meier (KM) curve. To further analyze the effect of GIMAPs and clinical parameters (tumor stage, tumor grade, tumor longest dimension, gender and age) on ccRCC patients’ prognosis, patients were randomly split into two sets (a training set: 322 patients (61%) and a test set: 205 patients (39%)). In the training set, Lasso regression analysis was used to select the genes in GIMAPs to construct the prognostic risk model, then risk scores of patients both in the training set and test set were calculated. Cox regression analysis was utilized to determine the independent prognostic factors in risk scores and clinical parameters, then the independent prognostic factors were used to construct a nomogram. The performance of the risk model and nomogram were comprehensively evaluated using the KM curve, calibration curve, concordance index, and ROC curve.

2.4 Immune Infiltration Analysis

The interrelationship among GIMAPs family members and between GIMAPs expression and immunocytes in ccRCC was evaluated using TIMER (https://cistrome.shinyapps.io/timer/, accessed on 21 October 2021). The amount of infiltrating immunocytes were estimated by the TIMER algorithm. The correlation was identified by Spearman correlation analysis. TIMER is online tool that can analyze immunocytes infiltration of different cancers from TCGA [31]. TME scores of ccRCC in TCGA were downloaded from ESTIMATE (https://bioinformatics.mdanderson.org/estimate/, accessed on 18 October 2022), immunocytes infiltration scores of ccRCC were calculated using CIBERSPRT algorithms [32].

2.5 Co-Expression Analysis and Function Enrichment Analysis

The co-expression analysis was carried out using LinkedOmics (http://linkedomics.org, accessed on 2 November 2021) to obtain the genes co-expressed with GIMAPs in ccRCC. The ccRCC patients’ data were from TCGA, and the statistical method was Pearson correlation analysis. LinkedOmics is a web service that can analyze various tumor patients’ clinical parameters and gene expressions in CPTAC and TCGA. Webgestalt (http://www.webgestalt.org, accessed on 2 November 2021) was used for gene set enrichment analysis (GSEA) of the genes co-expressed with GIMAPs. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were chosen for GSEA, Rank Criteria: Pearson correlation coefficient (PCC); Minimum Number of Genes: 3; Simulations: 1000. WebGestalt is an online tool that can perform functional enrichment analysis of genes of multiple species [33].

2.6 Statistical Analyses

The Student’s t-test in SPSS (25.0.0.0, IBM Corp., Chicago, IL, USA) was used to compare the GIMAPs mRNA expression in the ccRCC tumor and normal renal tissues. In evaluating the value of GIMAPs for ccRCC’s diagnosis, “plotROC” package was used to plot ROC curves and calculated Area Under Curve (AUC), and the “glm” function was used for the multivariate logistic regression analysis. In evaluating the impact of GIMAPs on the prognosis of ccRCC patients, “glmnet” package was used to perform Lasso regression analysis, and “survival”, “survminer” and “survivalROC” packages were used to perform Cox regression analysis, plot Kaplan Meier survival curves, nomogram, the ROC curves and calibration curves, and calculate concordance index. R (v.4.2.1) was used to perform the statistical analysis, p $<$ 0.05 signified statistical significance.

3. Results

3.1 GIMAPs Were Overexpressed in ccRCC Tumor Tissues

In ONCOMINE, Yusenko et al. [34] and Lenburg et al. [35] Renal dataset both showed higher mRNA expression of GIMAP2, GIMAP4, GIMAP5, GIMAP6 and GIMAP7 in ccRCC tumor tissues than normal renal tissues (p $<$ 0.05), The mRNA expression of GIMAP8 in the two tissues did not differ significantly, and the Yusenko Renal dataset showed that ccRCC tumor tissues had higher GIMAP1 mRNA expression than normal renal tissues (p $<$ 0.05), Gumz et al. [36] and Jones et al. [37] Renal dataset both showed that ccRCC tumor tissues had higher mRNA expression of GIMAP4, GIMAP5, and GIMAP6 than normal renal tissues (Table 1)[34, 35, 36, 37].

Table 1.The difference of GIMAPs’ mRNA expression between normal renal tissues and ccRCC tumors.

GIMAPs	ccRCC tumor cases	Normal renal cases	Fold Change	p-value	t-score	Dataset
GIMAP1	26	5	2.622	4.29 × 10 ${}^{–4}$ *	4.554	Yusenko Renal
GIMAP2	9	9	2.536	9.24 × 10 ${}^{–4}$ *	4.326	Lenburg Renal
GIMAP2	26	5	4.974	3.98 × 10 ${}^{–7}$ *	12.909	Yusenko Renal
GIMAP4	10	10	3.300	8.75 × 10 ${}^{–5}$ *	5.008	Gumz Renal
GIMAP4	23	23	2.317	1.93 × 10 ${}^{–13}$ *	10.214	Jones Renal
GIMAP4	9	9	2.112	2.32 × 10 ${}^{–5}$	6.106	Lenburg Renal
GIMAP4	26	5	3.300	1.40 × 10 ${}^{–8}$ *	9.996	Yusenko Renal
GIMAP5	10	10	2.119	6.05 × 10 ${}^{–5}$ *	5.432	Gumz Renal
GIMAP5	23	23	1.440	5.38 × 10 ${}^{–7}$ *	6.170	Jones Renal
GIMAP5	9	9	1.927	6.12 × 10 ${}^{–5}$ *	5.496	Lenburg Renal
GIMAP5	26	5	2.599	0.001*	4.924	Yusenko Renal
GIMAP6	10	10	4.199	5.73 × 10 ${}^{–7}$ *	7.705	Gumz Renal
GIMAP6	23	23	2.325	6.32 × 10 ${}^{–8}$ *	6.307	Jones Renal
GIMAP6	9	9	2.374	0.002*	3.667	Lenburg Renal
GIMAP6	26	5	3.321	3.07 × 10 ${}^{–9}$ *	8.988	Yusenko Renal
GIMAP7	9	9	1.814	0.019*	2.314	Lenburg Renal
GIMAP7	26	5	3.364	1.63 × 10 ${}^{–4}$ *	5.671	Yusenko Renal
GIMAP8	9	9	–1.037	0.758	–0.719	Lenburg Renal
GIMAP8	26	5	1.033	0.469	0.081	Yusenko Renal

* p $<$ 0.05. ccRCC, clear cell renal cell carcinoma; GIMAP, guanosine triphosphatases of immunity-associated proteins.

The GIMAPs mRNA expression pattern was verified in TCGA. CcRCC tumor tissues showed higher mRNA expression of GIMAP1 compared with normal renal tissues (p $<$ 0.001) (Fig. 1A), GIMAP2 (p $<$ 0.01) (Fig. 1B), GIMAP4 (p $<$ 0.01) (Fig. 1C), GIMAP5 (p $<$ 0.001) (Fig. 1D), GIMAP6 (p $<$ 0.001) (Fig. 1E), GIMAP7 (p $<$ 0.001) (Fig. 1F) and GIMAP8 (p $<$ 0.001) (Fig. 1G) (Supplementary Table 1). To evaluate the ability of GIMAPs expression level to predict the type of samples (normal kidney tissue or ccRCC tissue), we plotted ROC curves, the results showed AUC of all ROC curves were more than 0.7 (Fig. 1H). We performed multivariate logistic regression analysis with 7 genes in the GMIPA family to determine their predictive ability for the diagnosis of ccRCC, the results indicated GIMAP1, GIMAP2, GIMAP5, GIMAP6 and GIMAP8 were independent diagnostic factors (p $<$ 0.05) (Table 2).

Fig. 1.

The expression of GIMAPs mRNA in ccRCC tumor and normal renal tissues. (A–G) Compared to normal renal tissues, the ccRCC tumor had higher mRNA expression of GIMAPs. (H) ROC curves showed high sensitivity and specificity of these GIMAPs in distinction between normal renal tissues and ccRCC tumors. ** p $<$ 0.01, *** p $<$ 0.001. ccRCC, clear cell renal cell carcinoma; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Table 2.Diagnostic value of GIMAPs in ccRCC.

Gene symbol	Logistic regression analysis
Gene symbol	OR	95% CI	p‑value
GIMAP1	0.0269	0.0036–0.1673	$<$ 0.001*
GIMAP2	1.1541	0.4149–3.3329	$<$ 0.001*
GIMAP4	4113.901	414.5919–70,187.91	0.787
GIMAP5	827.8691	80.2729–15,301.57	$<$ 0.001*
GIMAP6	0.0024	0.0003–0.0136	$<$ 0.001*
GIMAP7	4.1545	1.0478–19.7302	0.054
GIMAP8	0.0675	0.0158–0.2911	$<$ 0.001*

* p $<$ 0.05. OR, odds ratio.

We then studied the GIMAPs protein expression of ccRCC in CPTAC. There was no total protein but only phosphoprotein data about GIMAP5. The expression of GIMAP5 phosphoprotein in normal renal tissues and ccRCC was not significantly different (Fig. 2D), but GIMAP1 (Fig. 2A), GIMAP2 (Fig. 2B), GIMAP4 (Fig. 2C), GIMAP6 (Fig. 2E), GIMAP7 (Fig. 2F), and GIMAP8 (Fig. 2G) total proteins were expressed at lower levels in normal renal tissues than ccRCC (p $<$ 0.001). The results showed that all GIMAPs mRNA and 6 GIMAPs (GIMAP1, GIMAP2, GIMAP4, GIMAP6, GIMAP7 and GIMAP8) protein were over-expressed in tumor tissues of ccRCC.

Fig. 2.

GIMAPs protein expression in normal renal tissues and ccRCC tumors. (A–C,E–G) Most GIMAPs (GIMAP1, GIMAP2, GIMAP4, GIMAP6, GIMAP7, and GIMAP8) had higher total protein expression in ccRCC tumors than in normal renal tissues. (D) It had no significant difference between GIMAP5 phosphoprotein expression of normal renal tissues and ccRCC tumors. *** p $<$ 0.001.

3.2 Relationship among GIMAP Family Members

First, we used STRING to detect the interaction between GIMAPs protein, and found that all combined scores except with GIMAP2 were more than 0.6 (Fig. 3A, Supplementary Table 2). Then, the relationship among GIMAP family members in ccRCC was analyzed. All Spearman correlation coefficients (SCCs) between GIMAP members were more than 0.5 (p $<$ 0.001), and all SCCs except with GIMAP2 were more than 0.8 (p $<$ 0.001) (Fig. 3B). It can be seen that GIMAP family members were closely related to each other.

Fig. 3.

Relationship between GIMAP family members. (A) Protein-protein interaction network diagram of GIMAP family members. (B) Heatmap of spearman correlation coefficients between GIMAP family members in ccRCC.

3.3 Impact of GIMAPs on ccRCC Patients’ Prognosis

There were 322 cases in the training set and 205 cases in the test set (Table 3). Patients with higher GIMAP1 (Fig. 4A), GIMAP2 (Fig. 4B), GIMAP4 (Fig. 4C), GIMAP5 (Fig. 4D), GIMAP6 (Fig. 4E), GIMAP7 (Fig. 4F), and GIMAP8 (Fig. 4G) expression levels had longer OS than those with lower expression levels, according to survival analysis (p $<$ 0.05).

Table 3.Clinical traits of patients with ccRCC in the training and test sets.

Clinical factor		Training set (n = 322)	Test set (n = 205)	Overall (n = 527)
Survival time (day)		1354.97 $\pm$ 1006.87	1313.23 $\pm$ 940.19	1338.74 $\pm$ 980.78
Status
	Live	219	136	355
	Dead	103	69	172
Gender
	Female	114	69	183
	Male	208	136	344
Age (year)
	$<$ 50	65	42	107
	50–59	83	55	138
	60–69	96	54	150
	$\geq$ 70	78	54	132
Tumor longest dimension (cm)
	0–0.9	17	12	29
	1–1.9	167	106	273
	2–2.9	49	35	84
	$>$ 3	21	15	36
	Missing	69	36	105
Tumor grade
	1	9	4	13
	2	133	93	226
	3	128	77	205
	4	48	27	75
	Missing	4	4	8
Tumor Stage
	I	162	101	263
	II	35	22	57
	III	74	48	122
	IV	50	32	82
	Missing	1	2	3

Fig. 4.

Kaplan Meier curves showed the impact of GIMAPs on OS of ccRCC patients. (A–G) Patients with higher expression levels of all GIMAPs had longer OS than those with lower expression levels. OS, overall survival.

Because of the close correlation between GIMAP family members, a risk model base on GIMAPs was constructed by Lasso regression to comprehensively evaluate the impact of GIMAPs on prognosis of patients with ccRCC. In the training set, when the model lambda ( $\lambda{}$ ) was at minimum value, the optimal prediction model based on 6 genes was constructed (Fig. 5A,B). Ln (risk score) = 0.1477192 $\times{}$ GIMAP2 + 0.5610788 $\times{}$ GIMAP4 + 0.3883533 $\times{}$ GIMAP5 – 0.6306035 $\times{}$ GIMAP6 + 0.1594977 $\times{}$ GIMAP7 – 0.6300014 $\times{}$ GIMAP8. In order to analyze the prognostic value of the risk model, the patients were separated into two groups by the median risk scores (High risk and Low risk). In the training set, the Low risk group had more patients alive and higher GIMAPs expression than the High risk group (Fig. 5C), similar population distribution and GIMAPs expression were found in the test set (Fig. 5D). In the training set, compared to the patients in High risk group, the patients in Low risk group had longer OS (p $<$ 0.0001) (Fig. 5E), ROC curve showed the sensitivity and specificity of the model for prognostic prediction. AUC of 1 year, 2 year and 3 year survival were 0.652, 0.705 and 0.72 respectively (Fig. 5F). In the test set, the patients in the Low risk group also had longer OS than those in the High risk group (p $<$ 0.0001) (Fig. 5G), AUC of 1 year, 2 year, and 3 year survival were 0.654, 0.653 and 0.665 respectively (Fig. 5H).

Fig. 5.

GIMAPs-based risk model for ccRCC patients’ prognosis. (A) The coefficients of GIMAPs when the model lambda ( $\lambda{}$ ) was minimum value. (B) Partial likelihood deviation plot. (C,D) Population distribution and GIMAPs’ expression of Low and High risk groups in the test and train sets. (E,G) In both the train and test sets, the patients in Low risk group had longer OS than those in High risk group. (F,H) In train and test sets, ROC curves both showed high sensitivity and specificity of the risk score for prognostic prediction of patients with ccRCC.

High risk score, old age, high tumor grade and stage were found to be unfavorable for OS of ccRCC patients according to the Cox regression analysis based on risk score and clinical parameters (p $<$ 0.05), and these four parameters were independent prognostic factors in ccRCC (p $<$ 0.05) (Table 4). To further assess the risk models’ prognostic value, we analyzed the impact of risk scores on patients’ OS in different clinical subgroups, in all clinical subgroups (Grade 1–2, Grade 3–4, Age $<$ 60, Age $\geq$ 60, Stage I–II, Stage III–IV) and in total patients. Patients in Low risk groups had longer OS than patients in High risk groups (p $<$ 0.05) (Fig. 6A–G). Then these independent prognostic factors were used to construct a nomogram (Fig. 6H). The ROC curves indicated that the nomogram was sensitive to predict patients’ prognosis (AUC of 1, 3 and 5 year survival: 0.858, 0.812 and 0.776) (Fig. 6I). Excellent consistency of the nomogram predictions were demonstrated in the calibration curves (Fig. 6J). The concordance index of the nomogram was 0.776 (95% CI: 0.758–0.793), indicating good accuracy of the prognostic model.

Table 4.Univariate and multivariate Cox regression analysis with OS of ccRCC.

		Univariate analysis			Multivariate analysis
Covariates		HR	95% CI	p-value	HR	95% CI	p-value
Gender
	(male vs. female)	0.946	0.693–1.291	0.724
Age (year)
	(51–59 vs. $<$ 50)	1.569	0.915–2.691	0.101	1.296	0.751–2.237	0.351
	(60–69 vs. $<$ 50)	1.863	1.107–3.135	0.019*	1.282	0.752–2.184	0.361
	( $\geq$ 70 vs. $<$ 50)	3.076	1.862–5.083	$<$ 0.001*	2.569	1.543–4.277	$<$ 0.001*
Stage
	(II vs. I)	1.221	0.657–2.267	0.528	1.192	0.638–2.227	0.583
	(III vs. I)	2.611	1.734–3.931	$<$ 0.001*	1.877	1.220–2.887	0.004*
	(IV vs. I)	6.467	4.412–9.478	$<$ 0.001*	4.471	2.856–6.999	$<$ 0.001*
Tumor longest dimension		1.242	0.996–1.548	0.054
Tumor grade		2.317	1.890– 2.842	$<$ 0.001*	1.395	1.106–1.760	0.005*
Risk score		3.199	2.350–4.355	$<$ 0.001*	1.914	1.363–2.686	$<$ 0.001*

* p $<$ 0.05. HR, hazard ratio; OS, overall survival.

Fig. 6.

The impact of risk scores on OS of ccRCC patients in TCGA. (A–G) In all clinical subgroups and total patients, the patients in Low risk groups had longer OS than patients in High risk groups. (H) Predictive nomogram base on risk score and clinical parameters for OS in ccRCC patients. (I) ROC curve respecting nomogram’s prognostic prediction’s specificity and sensitivity. (J) Calibration curve respecting the accuracy of nomogram for predicting overall survival (OS) at 1 year, 3 years and 5 years.

3.4 GIMAPs Were Closely Correlated with Immunocytes Infiltration in ccRCC

All GIMAPs’ expression had negative correlation with the purity of tumor (p $<$ 0.001), which suggested that GIMAPs were mainly expressed in the tumor microenvironment. The SCCs between GIMAPs expression and B cell infiltration were less than 0.5 (p $<$ 0.001). All GIMAPs except of GIMAP8 had high correlation with CD8 ${}^{+}$ T cell (SCC $>$ 0.5, p $<$ 0.001). GIMAP1 and GIMAP8 had high correlation with CD4 ${}^{+}$ T cell infiltration (SCC $>$ 0.5, p $<$ 0.001). GIMAP2 had high correlation with macrophages (SCC $>$ 0.5, p $<$ 0.001). GIMAP2, GIMAP4 and GIMAP6 had high correlation with neutrophils (SCC $>$ 0.5, p $<$ 0.001). GIMAP2 and GIMAP4 had high correlation with DCs (SCC $>$ 0.5, p $<$ 0.001) (Fig. 7).

Fig. 7.

Correlation between GIMAPs and immunocytes infiltration in ccRCC. All GIMAPs mRNA expression were negative to the purity of tumor, indicated GIMAPs were mainly expressed in tumor immune microenvironment. The degree of immunocyte infiltration (B cell, CD4 ${}^{+}$ T cell, CD8 ${}^{+}$ T cell, macrophage, dendritic cell, neutrophil, and dendritic cell) was closely correlated with the expression of GIMAPs.

The risk model based on GIMAPs was closely correlated with the TIME of ccRCC, the risk score was highly associated with immune score (SSC = 0.23, p $<$ 0.001) (Fig. 8A) and negatively correlated with stromal score (SSC = –0.33, p $<$ 0.001) (Fig. 8B). No significant correlation was found between the risk score and the ESTIMATE score (Fig. 8C). The main immunocytes infiltrating in ccRCC included M2 macrophages, resting memory CD4 ${}^{+}$ T cells, M1 macrophages, CD8 ${}^{+}$ T cells, activated NK cells and so on. Comparing with the High risk group, there were more M2 macrophages, M1 macrophages, resting memory CD4 ${}^{+}$ T cells and naive B cells, less CD8 ${}^{+}$ T cells, Tregs, gamma delta T cells, plasma cells and M0 macrophages in Low risk group (p $<$ 0.05) (Fig. 8D,E). In the common immune checkpoints, compared with the High risk group, the Low risk group had lower expression of PD1 and CTLA4 (p $<$ 0.05) (Fig. 8F,H), higher expression of PDL1 (p $<$ 0.05) (Fig. 8G) in the Low risk group. These data demonstrated that GIMAPs have an important impact on immune cell infiltration in ccRCC.

Fig. 8.

Relationship between the immune microenvironment of ccRCC and the risk score. (A–C) Risk scores’ relationship with stromal score, immune score and ESTIMATE score. (D) Heatmap showed immunocytes infiltration in ccRCC of Low and High risk groups. (E) The comparison of immunocytes infiltration in the Low and High risk groups. (F–H) The comparison between the immune checkpoints (CTLA4, PD1 and PDL1) expression in High and Low risk groups. * p $<$ 0.05, ** p $<$ 0.01, *** p $<$ 0.001, *** p $<$ 0.0001, ns p $>$ 0.05.

3.5 Biological Function and Pathways of GIMAPs in ccRCC

We used co-expression analysis to identify genes that correlated with GIMAPs in ccRCC, and then performed GSEA of genes co-expressed with each GIMAP. The results of co-expression analysis are shown in Supplementary Table 3. There were total 93 KEGG pathways (p $<$ 0.05, FDR $<$ 0.05) obtained by GSEA (Fig. 9, Supplementary Table 4). All GIMAPs were strongly associated with the differentiation of helper T cells, the cytotoxicity associated with natural killer cell and so on. Pathways that negatively correlated with all GIMAPs include oxidative phosphorylation, citrate cycle and so on. It was demonstrated that lots of pathways correlated with GIMAPs play important roles in TIME.

Fig. 9.

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways correlated with GIMAPs in ccRCC. (A) Heatmap of pathways correlated with GIMAP family members in ccRCC. (B) The common pathways positively correlated with GIMAP family members. NES, normalized enrichment score; FDR, false discovery rate. *** p $<$ 0.001.

4. Discussion

Recent research showed that TIME was crucial to the occurrence and development of multiple tumors [38]. CcRCC is considered as an immunogenic tumor characterized of abundant immunocytes infiltration, in particular T cell infiltration [39]. Different types and functional states of immunocytes have different effects on ccRCC [40]. Therefore, identification of key factors affecting TIME is beneficial to predict the prognosis and explore new treatments for ccRCC patients. In this study, we found that GIMAPs were over-expressed in TIME of ccRCC, and they had close correlation with clinical outcomes and infiltration of immunocytes in this tumor.

The GIMAP family proteins are associated with immunity, which all have binding domains for GDP/GTP. GIMAPs can regulate biological functions and the states of a variety of immunocytes [17, 18, 19, 20, 21, 22]. GIMAPs are closely related to the autoimmune regulation of diabetes, allergy and asthma [41]. The dysregulation of GIMAPs is not only related to immune related diseases, but also various tumors. The blood and tumor tissues of hepatocellular carcinoma patients showed down-regulated levels of GIMAP6 and GIMAP5 expression, suggesting GIMAP6 and GIMAP5 possibly participated in the pathogenic mechanism of hepatocellular carcinoma [24]. It has been identified that the majority of GIMAPs (GIMAP1, GIMAP4, GIAMP6, GIMAP7 and GIMPA8) were down-regulated in endometrial cancer, and that low GIMAPs expression were associated with a poor prognosis and closely linked to immunocytes infiltration [25]. Studies indicated that GIMAPs showed low expression in lung adenocarcinoma tissues, and low GIMAPs expression were closely related to poor clinical outcomes [23, 26, 42]. Similarly, our study showed that the prognosis for patients with low GIMAP expression was poor, and in different clinical subgroups, the patients in the Low risk groups had longer OS than those in High risk groups, suggesting that the GIMAPs-based risk model has a good capacity for predicting the prognosis of patients with ccRCC. GIMAP family members were negatively correlated with tumor purity of lung adenocarcinoma [26]. In our research, GIMAP family members had the similar correlation with tumor purity of ccRCC, indicating that GIMAPs were mainly located in TIME, and that the risk score had a positive relationship with the immune score and a negative relationship with the stromal score. We found that GIMAPs’ high expression reflected the characteristics of high immunocytes infiltration in ccRCC. However, different from other tumors, GIMAPs’ expression was notably up-regulated in ccRCC. High immunocytes infiltration in ccRCC dose not always means favorable prognosis for patients. The types and states of infiltrating immunocytes are key factors for prognosis of ccRCC patients [43]. For this reason, we analyzed the biologic function of GIMAPs in TIME.

CcRCC is infiltrated by abundant immunocytes, mainly DCs, macrophages, NK cells, CD4 ${}^{+}$ T cells, and CD8 ${}^{+}$ T cells [39]. Our study yielded similar results, the main immunocytes infiltrated in ccRCC were CD4 ${}^{+}$ T cells, CD8 ${}^{+}$ T cells, macrophages, NK cells and Tregs. As the most important antitumor immune cell, the amount and antitumor activity CD8 ${}^{+}$ T cells are significantly related to the immunotherapeutic effect and clinical prognosis in tumors [44]. Studies have reported that GIMAP1, GIMAP5 and GIMAP6 are important to maintain the amount of CD8 ${}^{+}$ T cells [19, 22, 45]. GIMAP1, GIMAP4, GIAMP6, GIMAP7 and GIMPA8 are negatively related to CD8 ${}^{+}$ T cells infiltration in endometrial cancer [25]. However, all GIMAP family members are a positively related to CD8 ${}^{+}$ T cells infiltration in lung cancer [26], GIMAP7 also has positive correlation with CD8 ${}^{+}$ T cells infiltration in pancreatic adenocarcinoma [46]. Our study indicated that all GIMAPs had positive correlation with CD8 ${}^{+}$ T cells infiltration. It can be seen that GIMAPs had different relationship with CD8 ${}^{+}$ T cells infiltration in different cancers, and that GIMAPs are important to regulate the amount of CD8 ${}^{+}$ T cells in TIME of ccRCC. Multiple cells regulated CD8 ${}^{+}$ T cells’ antitumor activity, which included DCs, helper T cells, Tregs, macrophages and so on. CcRCC is reported to secrete cytokines, which affect differentiation of DCs, resulting in decrease or loss of CD8 ${}^{+}$ T cells’ antitumor activity [47]. Mature DCs are linked to the activation of CD8 ${}^{+}$ T cells and ccRCC’s favorable prognosis [40, 48]. M1 macrophages are able to increase CD8 ${}^{+}$ T cell cytotoxicity by expressing IL-12, IFN- $\gamma{}$ and TNF [49], M2 macrophages can cooperate with Tregs to prevent CD8 ${}^{+}$ T cell migration to tumor cells [50]. TAMs have high heterogeneity, and targeting specific TAM subgroups may be a good treatment strategy for ccRCC [43]. In the treatment of ccRCC, bevacizumab not only could decrease the microvascular density of tumors, but also could reduce the the number of CD68 ${}^{+}$ macrophages [51]. Tregs can suppress antitumor activity of CD8 ${}^{+}$ T cells by releasing TGF- $\beta{}$ [52]. In this study, we found that GIMAPs were positively correlated with DCs and macrophages infiltration. Compared to the Low risk group, there were less M1 and M2 macrophages and more Tregs infiltration in the High risk group. It was likely that low M1 macrophages and high Tregs infiltration would decrease CD8 ${}^{+}$ T cells’ antitumor activity, resulting in poor outcome of patients with ccRCC. Studies indicated Th1 cells produced IL-2, IFN- $\gamma{}$ and LT- $\alpha{}$ to participate in the cell-mediated immune response, while the increase of Th2 and Tregs inhibited the immune response and were related to poor prognosis in ccRCC [13, 53]. It has been reported that the lack of GIMAP5 leaded to abnormal differentiation of helper T cells, resulting in an increase of pathogenic Th2 and Th17 cells to promote allergic airway disease [54]. Our study showed GIMAPs had a significant positive correlation with the differentiation of Th1, Th2, and Th17 cells as well as the chemokine signaling pathway in ccRCC. Adhesion molecules were required for antigen-presenting cells to activate T cells [55, 56]. It has been reported that adhesion molecules can regulate the polarization of Th1/Th2, and the maintenance and development of Tregs and T cells [57, 58, 59]. In our study GIMAPs were closely related to cell adhesion molecules. It is high possible that GIMAPs can regulate CD8 ${}^{+}$ T cells’ antitumor activity in TIME of ccRCC.

NK cells are another important antitumor immunocytes in TIME. Studies showed that IL-2 inhibited development of ccRCC by increasing proliferation and cytotoxicity of NK cells [60]. A study indicated that a high proportion of NK cells in TIME had correlated with favorable prognosis of ccRCC patients [61]. As important immunosuppressive cells, Tregs can suppress the cytotoxicity of NK cells by releasing TGF- $\beta{}$ . Our results showed NK cells were one of the main immunocytes infiltrating in ccRCC. There were no significant differences in activated NK cells infiltration between High and Low risk groups. The High risk group had more Tregs than the Low risk group, Tregs might inhibit the NK cells’ anti-tumor activity. GIMAP5 is important to maintain NK cells in peripheral blood and lymphoid organs [18]. In our study, all GIMAPs were positively related to NK cell mediated cytotoxicity. We found that GIMAPs are able to inhibit the development of ccRCC by regulating antitumor activity of NK cells.

Immunotherapy has become an important method for the treatment of ccRCC [62, 63], but only small proportion of patients benefit from it [64]. Studies have shown that there are a large number of attenuated and functional defective CD8 ${}^{+}$ cell infiltration in ccRCC [39]. PD1 and CTLA4 can suppress CD8 ${}^{+}$ T cells’ activation and increase the depletion of T cell [65, 66]. A poor prognosis is associated with high PD1 expression and more CD8 ${}^{+}$ T cells infiltration in ccRCC [67]. In our study, compared with the Low risk group, there were more CD8 ${}^{+}$ cell infiltration, higher expression of PD1 and CDLA4 expression in High risk group, the patients in the High risk group had shorter OS. These results indicated that the model based on GIMAPs had good prognosis prediction capacity, and activating the anti-tumor activity of CD8 ${}^{+}$ T cells by regulation of GIMAPs may be a new treatment strategy for ccRCC.

Metabolic reprogramming is an important signature of ccRCC, including increased aerobic Glycolysis, decreased mitochondrial Oxidative phosphorylation, increased Fatty acid metabolism and so on [68, 69, 70]. It can provide sufficient energy and substances for tumor growth, meanwhile, it is beneficial for tumors to adapt to hypoxic environments, resist oxidative stress, and evade host immune surveillance [71]. Study showed that ccRCC patients with high levels of glycolytic enzymes had lower progression free survival and cancer specific survival than patients with low levels of glycolytic enzymes [72]. Inhibition of aerobic glycolysis by 2-DG could reduce the proliferation and activity of low-grade ccRCC, and promotion of fatty acid oxidation by Etomoxir could inhibit the proliferation and activity of high-grade ccRCC [73]. Studying tumor metabolism reprogramming are important to find new strategies for the diagnosis and treatment of ccRCC. Inhibition of NDUFA4L2 could reduce the vitality of ccRCC cells, increase mitochondrial mass, and induce ROS production during hypoxia [74]. Depletion of MUC1 could inhibit the migration and proliferation of ccRCC cells [75]. A study about lung cancer showed that serum metabolomic fingerprints could serve as a “collective” biomarker for predicting immune checkpoint inhibitor responses, capable of predicting individual treatment outcomes with an accuracy of $>$ 80% [76]. Our research showed that GIMAPs were negative correlatted with a large number of energy metabolism pathways, including “Oxidative phonology”, “Propanoate metabolism”, “Pyruvate metabolism”, “Cysteine and methionine metabolism” and so on. GIMAPs may affect the growth and metastasis of ccRCC by regulating its metabolism.

In this study, we comprehensively explored the role of GIMAPs in ccRCC. Previous studies and our own researches all showed that GIMAPs were important to regulate TIME. As an immunogenic tumor, ccRCC was significantly affected by GIMAPs, which could inhibit the development of tumor by increasing the amount and antitumor activity of infiltrating immunocytes.

There are some limitations in this research. All data is derived from databases, and the total samples sizes are limited. More samples are needed to verify our results, and further experiments are essential to reaffirm the role of GIAMPs in ccRCC.

5. Conclusions

Our research indicated that GIMAPs were over-expressed in ccRCC, and that the GIMAPs’ expression was closely related to the prognosis of ccRCC patients. In addition, GIMAPs were highly correlated with the immunocytes infiltration in TIME. We suggest that GIMAPs are potential prognostic biomarkers and therapeutic targets of ccRCC.

Availability of Data and Materials

All data produced and detailed in this article is accessible from open databases. On reasonable request, the corresponding author will provide additional information.

Author Contributions

XML and MZ conceived and designed the project. MJZ and XZ collected the data. JH and XML performed the interpretation of data. MZ and MJZ performed the statistical analysis. MJZ and MZ wrote the manuscript. XML, JH and MZ revised the article. All authors read and approved the final manuscript.

Ethics Approval and Consent to Participate

Not applicable.

Acknowledgment

We appreciate the generosity of the researchers from ONCOMINE, TCGA and CPTAC for sharing the huge amount of data.

Funding

This research received no external funding.

Conflict of Interest

The authors declare no conflict of interest.

Supplementary Material

Supplementary material.rar

References

[1]

Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: a Cancer Journal for Clinicians. 2021; 71: 209–249.