This article is mentioned in:

Introduction

GBM, a grade IV glioma according to the World Health Organization (WHO) classification, is the most malignant primary brain tumor with strong invasive potential and rapid recurrence after surgery. Although aggressive treatment strategies have been used to prolong overall survival (OS) in patients with GBM, the median 5-year OS rate is <20% (1). With the development of molecular biology techniques in recent years, high-throughput sequencing and microarray techniques have been used to improve our understanding of the molecular events in GBM (2). A combination of histological and molecular characteristics of glioma, including the presence of isocitrate dehydrogenase-1 (IDH-1) mutations and 1p/19q codeletion, are associated with GBM according to the 2016 WHO classification (2). The identification of this combination has improved the accuracy of glioma diagnosis, treatment and prognostic evaluation (2). Based on data from The Cancer Genome Atlas (TCGA) Research Network, GBM is divided into four subtypes, mesenchymal, classical, neural and proneural, with different cellular features and genetic contexts (3). The molecular characterization of glioma has improved patient stratification and provided insight into novel strategies for treating GBM. However, exploring the molecular events involved in GBM is still necessary for developing a targeted therapy (4).

Vacuolar (V)-ATPase is a complex consisting of multiple subunits that play roles in multiple biological processes in mammalian cells, intracellular membrane-associated V-ATPase can acidifies lysosomes, endosomes and secretory vesicles, and influences numerous processes associated with these organelles, including vesicular trafficking, endocytosis, autophagy, receptor recycling and protein degradation (5). V-ATPase is comprised of two functional domains: The peripheral V1 domain, which consists of eight subunits (A-H) and is mainly responsible for ATP hydrolysis, and the V0 domain, which consists of five subunits (a, c, c', d and e) and is responsible for proton translocation. The main function of V-ATPase is to transport protons into intracellular compartments, promote the acidification of endosomes and lysosomes and excrete intracellular protons into the extracellular space, thus maintaining H+ homeostasis (5). Multiple studies have demonstrated that V-ATPase plays critical roles in several cancer types, especially in tumor invasion and migration (5–10). For example, V-ATPase is overexpressed at the plasma membrane of invasive MB231 human breast cancer cells (7), and inhibiting V-ATPase using proton pump inhibitors or small interfering RNAs can suppress cancer cell line proliferation and metastasis in animal models (6,9,10). Overexpression of V-ATPase at the plasma membrane of invasive cancer cells facilitates the activation of proteinases under low pH conditions and further modifies components of the extracellular matrix, such as matrix metalloproteinases (8). An acidic pH in the tumor microenvironment can induce VEGF expression and recruit and polarize macrophages (11).

Subunit ‘a’ is comprised of four isoforms (a1-a4) and is located at the V0 domain, which is responsible for the subcellular localization of V-ATPase (5). The delivery of V-ATPase specifically to the plasma membrane of breast cancer cells relies on the overexpression of the V0a3 subunit, which is also known as T cell immune regulator 1 (TCIRG1). TCIRG1 normally localizes to lysosomes in osteoclasts and insulin-containing secretory vesicles in pancreatic β cells, suggesting that distinct subunit compositions support cancer-specific functions (5–7). Overexpression of TCIRG1 has also been reported in hepatocellular carcinoma, melanoma and breast cancer (7,9,10). The role of the V-ATPase a3 subunit in mouse cytotoxic T lymphocytes is responsible for the acidification of cytotoxic granules and contributes to the maturation and efficient transport of cytotoxic granules to the immune synapse (12). However, the expression profiles of TCIRG1 and its function in GBM prognosis, as well as the underlying mechanisms, remain to be understood (13,14). In the present study, data mining was used to characterize the profiles of TCIRG1 expression in glioma and other malignant tumors, and investigated the association between TCIRG1 expression levels and OS in patients with different GBM subtypes. Furthermore, the correlation between TCIRG1 expression and tumor-infiltrating immune cells in the environment surrounding GBM was explored using the Tumor Immune Estimation Resource (TIMER) and Gene Expression Profiling Interactive Analysis (GEPIA). The study revealed the important function of TCIRG1 in GBM and explained the potential relationship and mechanism of the interaction between TCIRG1 and tumor immunity.

Materials and methods

Gene expression array and online tool analysis

TCIRG1 mRNA expression in diverse cancer types was evaluated with the Oncomine database (https://www.oncomine.org/resource/main.html) with cut-off values for the probability and fold-change were 0.05 and 1.5, respectively. TCIRG1 mRNA expression data were further analyzed using the TIMER database (https://cistrome.shinyapps.io/timer/) (15), which contains >10,000 samples from 32 different kinds of cancer from TCGA for estimating tumor purity and the abundance of immune infiltrates, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells (DCs).

A gene expression array (AffyU133a) of TCGA-GBM was obtained by using the University of California Santa Cruz Xena browser (https://xenabrowser.net), which contains 529 GBM tissues and 10 non-tumor tissues. Gene expression microarray data from the GSE16011 dataset [eight normal brain tissues, 159 GBM tissues and 117 low-grade glioma (LGG) tissues] (16) and RNA sequencing data (443 cases of LGG, 249 cases of GBM and one case without a defined WHO grade) from Chinese Glioma Genome Atlas (CGGA) were downloaded through the GlioVis portal (http://gliovis.bioinfo.cnio.es/) (17) and CGGA portal (http://www.cgga.org.cn/), respectively. The clinical data included the demographic and survival data, IDH-1 status, GBM subtypes (3), and chemotherapy and radiation therapy information. GBM subtypes classification were performed by Verhaak methods (3), and the specificity and sensitivity of TCIRG1 expression with mesenchymal subtype GBM were evaluated by receiver operating characteristic (ROC) curves and the value of area under the curve (AUC). Genes that showed a strong correlation (correlation |r|≥0.5) with TCIRG1 in GBM samples from the TCGA-GBM and GSE16011 cohorts were calculated using Pearson's coefficient analysis (18), data derived from the GlioVis portal (17) and the R2: Genomics Analysis and Visualization platform (http://r2.amc.nl). GO category analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 (https://david.ncifcrf.gov/) (19). The protein-protein interaction (PPI) networks were constructed using the STRING online tool (version 11.0; http://string-bd.org/) (combined score ≥0.4 was used as the cut-off criterion) and further visualized using Cytoscape software (version 3.7.1) (20). The connectivity degree value of each node in the PPI network was calculated by the CytoHubba plugin, and the module analyses of the PPI network were performed by the MCODE plugin (node score cut-off=0.2; K-core=2; maximum depth=100; degree cut-off=2), in which modules contained >20 nodes were regarded as the key modules in the PPI network. The relationship between TCIRG1 and immune cell infiltration in GBM was evaluated using the TIMER database, and the immune and stromal scores were calculated using Estimation of STromal and Immune Cells in MAlignant Tumor Tissues using Expression data (ESTIMATE) (https://bioinformatics.mdanderson.org/estimate/) (15,21), which is widely used to evaluate immune scores and stromal scores in cancer (22–24). The immune gene markers that were significantly correlated (Spearman's correlation with P-value <0.05) with TCIRG1 expression in the TIMER database were further identified using the Gene Expression Profiling Interactive Analysis (GEPIA) online database (http://gepia.cancerpku.cn/index.html) (25). External validation of TCIRG1 expression based on immune signatures in GBM was performed using GSE16011 with the GlioVis online tool by Spearman's correlation coefficient analysis (16,17).

Clinical sample preparation

Fourteen tumor resection and peritumoral samples (fresh-frozen) were obtained from the Department of Neurosurgery of the First Hospital of China Medical University (Shenyang, China) from September 2018 to October 2019, including four samples of peritumoral tissues as a control group, four samples of LGG and six samples of GBM as a glioma group (patients age ranged from 40 to 65, median age 55, including six females and eight males). Tumor classification was based on the WHO grade classification, and the experiments were approved by The Ethics Committee of the First Hospital of China Medical University (Shenyang, China).

Western blotting (WB) experiment

Total protein was extracted from the aforementioned tissues using radioimmunoprecipitation buffer (Beyotime Institute of Biotechnology) and determined the protein concentration using an Enhanced BCA Protein Assay kit (Beyotime Institute of Biotechnology). The proteins in per lane were 120 µg loaded onto 12% gels, resolved using SDS-PAGE and transferred to a polyvinylidene fluoride membrane. After blocking in 5% bovine serum albumin (Beijing Solarbio Science & Technology Co., Ltd.) for 2 h at room temperature, the membranes were incubated overnight at 4°C with primary antibodies against TCIRG1 (1:1,000; cat. no. A15382; ABclonal Biotech Co., Ltd.) and GAPDH (1:4,000; cat. no. AC027; ABclonal Biotech Co., Ltd.). The secondary antibody, horseradish peroxidase-conjugated goat anti-rabbit IgG (1:10,000, cat. no. AS0014; ABclonal Biotech Co., Ltd.), was then applied for 2 h at 37°C. GAPDH served as the internal control for all WB experiments. Densitometry for WB was performed by ImageJ software (version 18.0; National Institutes of Health) and GraphPad Prism 7 (GraphPad Software).

Statistical analysis

High and low expression groups in each dataset were defined by TCIRG1 median expression value and the data were presented as mean ± standard deviation. GraphPad Prism 7 (GraphPad Software) was used to compare the differences in TCIRG1 expression among different groups. Kaplan-Meier survival analysis using log-rank tests were performed to compare the OS data between the two cohorts. Column analysis was performed using GraphPad Prism 7 using unpaired t-tests (two groups) and one-way ANOVA with Tukey's post hoc test. A multivariate Cox regression analysis was used to evaluate the prognostic value of TCIRG1 for the OS of patients with GBM using SPSS 21.0 (IBM Corp). Correlation analysis was performed by Pearson's coefficient analysis or Spearman analysis. All experiments were repeated three times. P<0.05 was considered to indicate a statistically significant difference.

Results

TCIRG1 expression in glioma and other types of cancer

Using the Oncomine database, the difference in TCIRG1 mRNA expression in different tumors and adjacent normal tissues was analyzed. TCIRG1 expression levels were increased in bladder, brain and central nervous system, breast, cervical, colorectal, esophageal, gastric, head and neck and kidney cancer, leukemia, liver and lung cancer, lymphoma, melanoma, myeloma and pancreatic cancer and sarcoma (P<0.05). In addition, in some datasets, TCIRG1 expression in breast cancer, leukemia, lung cancer, melanoma, pancreatic and prostate cancer and sarcoma was lower compared with that in normal tissues (Fig. 1A).

Figure 1. TCIRG1 expression in different types of cancer. (A) Expression levels of TCIRG1 in cancer and normal tissues in the Oncomine database. (B) TCIRG1 expression data from TCGA were compared using the TIMER database. Increased TCIRG1 expression was found in GBM according to data from (C) CGGA, (D) GSE16011 and (E) TCGA-GBM. (F and G) Protein levels of TCIRG1 were validated using western blot experiments, and increased expression was found in the glioma group. ****P<0.0001 vs. WHO III or LGG, **P<0.01 and *P<0.05 vs. respective control. TCIRG1, T cell immune regulator 1; ns, no statistical significance; CGGA, Chinese Glioma Genome Atlas; TCGA, The Cancer Genome Atlas, GBM, glioblastoma; WHO, World Health Organization; TPM, transcripts per million.

To further evaluate TCIRG1 expression in different human tumors, TCGA datasets were searched using the online tool TIMER (15). The results showed that TCIRG1 expression was higher in bladder cancer, breast cancer, cholangiocarcinoma, colon adenocarcinoma, esophageal carcinoma, head and neck cancer, kidney renal clear cell carcinoma, kidney renal papillary cell carcinoma, liver hepatocellular carcinoma, lung adenocarcinoma, rectum adenocarcinoma, stomach adenocarcinoma, thyroid carcinoma and uterine corpus endometrial carcinoma tissues compared with that in the corresponding normal tissues (Fig. 1B). However, low expression was found in lung squamous cell carcinoma (Fig. 1B), which was similar to the results in Oncomine.

An analysis of the mRNA levels of TCIRG1 in gliomas of different WHO grades showed that TCIRG1 expression was higher in GBM tissues compared with that in LGG and normal brain tissues (Fig. 1D and E). Through WB experiments, it was demonstrated that TCIRG1 was also more highly expressed in glioma tissues compared with that in control tissues (Fig. 1F and G). These results indicated that TCIRG1 expression could serve as a biomarker of glioma malignancy.

Increased TCIRG1 expression predicts a mesenchymal subtype of GBM

To explore the expression profiles of TCIRG1 in different subtypes of GBM, Verhaak classification methods were used (3). In both the TCGA-GBM and GSE16011 datasets, the highest TCIRG1 expression was observed in the mesenchymal subtype, and the proneural subtype had the lowest TCIRG1 expression (Fig. 2A and C). Furthermore, TCIRG1 expression and mesenchymal subtype data were used to generate ROC curves. In the TCGA-GBM cohort, the AUC was 0.867 (Fig. 2B). In the GSE16011 cohort, the AUC was 0.848 (Fig. 2D). These results indicated that TCIRG1 expression might serve as an indicator to predict the mesenchymal subtype of GBM.

Figure 2. Association between TCIRG1 expression and GBM subtypes. Mesenchymal subtype of GBM had the highest TCIRG1 expression in (A) TCGA-GBM and (B) receiver operating characteristic curves showed high sensitivity in predicting the mesenchymal subtype, with an AUC value of 0.867. Similar results were also found in (C) GSE16011, mesenchymal subtype GBM had the highest TCIRG1 expression and (D) the value of AUC was 0.848. ****P<0.0001 vs. mesenchymal or classical; ***P<0.001 vs. mesenchymal. TCIRG1, T cell immune regulator 1; TCGA, The Cancer Genome Atlas, GBM, glioblastoma multiforme; AUC, area under the curve.

High TCIRG1 expression is associated with poor OS in patients with GBM

The effect of TCIRG1 gene expression on the OS of patients with GBM was further analyzed. In TCGA-GBM and GSE16011 datasets, high TCIRG1 expression predicted poor OS in GBM patients (Fig. 3A and B). Multivariate Cox regression analysis was then used to investigate the independent role of TCIRG1 expression in the OS of patients with GBM. The results showed that after adjusting for age, sex, IDH-1 status, mesenchymal subtype, radiotherapy and chemotherapy, TCIRG1 still affected the OS of patients with GBM (Table I). These results indicated that TCIRG1 expression together with IDH-1 status, radiotherapy and chemotherapy are risk factors for the poor prognosis in GBM.

Figure 3. High TCIRG1 expression is associated with unfavorable overall survival for GBM. (A and B) Kaplan-Meier curve analyses showed that patients with GBM with high TCIRG1 expression were more likely to have poor OS in the TCGA-GBM and GSE16011 datasets. TCIRG1, T cell immune regulator 1; TCGA, The Cancer Genome Atlas, GBM, glioblastoma multiforme.

Table I. Multivariate Cox regression analysis in The Cancer Genome Atlas-Glioblastoma dataset.

Table I.

Multivariate Cox regression analysis in The Cancer Genome Atlas-Glioblastoma dataset.

	Univariate				Multivariate
Variables	HR	P-value	95% CI		HR	P-value	95% CI
TCIRG1 expression	1.583	0.003	1.173	2.138	1.525	0.007	1.124	2.068
Radiation therapy	5.521	<0.001	3.719	8.197	3.853	<0.001	2.486	5.972
Chemotherapy	2.916	<0.001	2.308	4.173	1.907	0.002	1.274	2.854
Mesenchymal subtype	1.267	0.146	0.921	1.743
Sex	0.884	0.450	0.643	1.216
Age	1.170	0.304	0.867	1.579
IDH-1 expression	3.685	0.004	1.511	8.991	2.761	0.027	1.124	6.781

[i] HR, hazard ratio; CI, confidence interval; TCIRG1, T cell immune regulator 1; IDH-1, isocitrate dehydrogenase-1.

GO analysis of TCIRG1-correlated genes in GBM

Pearson's correlation coefficient analysis was used to identify genes that showed strong positive and negative correlations (correlation ≥0.5) with TCIRG1 in the TCGA-GBM cohort using the GlioVis portal (17). In total, there were 244 genes with positive correlations and 76 genes with negative correlations (Table SI). Through the R2 platform, 217 genes were negatively correlated and 716 genes were positively correlated with TCIRG1 in GBM in the GSE16011 dataset (Table SII). Then all correlated genes were uploaded to the DAVID online tool to identify GO categories.

The top 10 items ranked according to the P-value for the GO categories are presented in Fig. 4. Genes that were positively correlated with TCIRG1 were mainly involved in the immune response. Regarding the biological processes category, genes that showed positive correlations with TCIRG1 were significantly enriched in the ‘immune response’, ‘immune effector process’, ‘defense response’, ‘regulation of the immune system process’ and ‘positive regulation of the immune system process’ (Fig. 4A). GO cellular component analysis revealed that genes that showed positive correlations with TCIRG1 were highly enriched in ‘membrane-bound vesicles’, ‘lysosomes’, ‘lytic vacuoles’, ‘extracellular region’ and ‘extracellular exosomes’ (Fig. 4B). GO molecular function analysis indicated that genes that showed positive correlations with TCIRG1 were enriched in ‘protein complex binding’, ‘cytokine receptor activity’, ‘cell adhesion molecule binding’, ‘receptor binding’ and ‘signal transducer activity’ (Fig. 4C). Moreover, genes that showed negative correlations were mainly involved in physiological functions, such as ‘nervous system development’ and ‘neuron development’ (Fig. 4D). For cellular component and molecular function, negative correlation genes were mainly enriched in microtuble, spindle, tubulin and microtubule binding (Fig. 4E and F). Similarly, GO category analyses of GSE16011 showed consistent results with the results for TCGA-GBM (Fig. S1). These results indicated that TCIRG1 mainly plays important roles in the host immune system in GBM.

Figure 4. GO analysis of genes that showed positive and negative correlations with TCIRG1 in The Cancer Genome Atlas-Glioblastoma cohort. Terms for (A) GO-BP, (B) GO-CC and (C) GO-MF of genes that showed positive correlations with TCIRG1. Terms for (D) GO-BP, (E) GO-CC and (F) GO-MF of genes that showed negative correlations with TCIRG1. Red nodes represent gene counts, and black bars represent -Log10 (P-value). GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function.

PPI network construction and module analyses of TCIRG1-associated genes

The correlated genes acquired from TCGA-GBM were used to create a PPI network. A total of 272 nodes and 1,584 edges were generated from the PPI network (Fig. 5A). The connectivity degree value of each node in the PPI network was calculated using the CytoHubba plugin (Table SIII). The genes in the PPI network with top 10 degree value were integrin subunit αM (ITGAM), integrin subunit β2 (ITGB2), toll-like receptor 2 (TLR2), cytochrome b-245 β chain (CYBB), CD44 molecule (CD44), intracellular adhesion molecule 1 (ICAM1), interleukin 10 receptor subunit α (IL10RA), spleen-associated tyrosine kinase (SYK), pleckstrin (PLEK) and C-C motif chemokine ligand 5 (CCL5). The MCODE plugin was used to select significant modules in the PPI network, and three modules were selected as the key modules in the PPI network (Fig. 5B-D). GO biological process analyses of these key modules were performed, and the top 10 items ranked by P-value are listed in Table II, which suggested that genes in the key modules were significantly associated with the immune response.

Figure 5. PPI network construction and key module selection in the PPI network. (A) PPI network generated with the genes correlated with T cell immune regulator 1 expression containing 272 nodes and 1,584 edges. (B-D) Three key modules were selected by the MCODE plugin with >20 nodes in each module. Node size was defined by the connectivity degree (calculated by the CytoHubba plugin, Table SIII), and color indicates Pearson's r value. PPI, protein-protein interaction.

Table II. GO biological process of key modules in the protein-protein interaction network.

Table II.

GO biological process of key modules in the protein-protein interaction network.

A, Module 1
Term	Gene count	P-value
GO:0006955~immune response	19	5.26×10−13
GO:0006952~defense response	16	1.41×10−9
GO:0006954~inflammatory response	12	1.76×10−9
GO:0031347~regulation of defense response	12	3.28×10−9
GO:0032103~positive regulation of response to external stimulus	9	6.72×10−9
GO:0031349~positive regulation of defense response	10	8.13×10−9
GO:0030595~leukocyte chemotaxis	8	1.91×10−8
GO:0050729~positive regulation of inflammatory response	7	2.78×10−8
GO:0070887~cellular response to chemical stimulus	18	4.29×10−8
GO:0007166~cell surface receptor signaling pathway	18	4.84×10−8
B, Module 2
Term	Gene count	P-value
GO:0006955~immune response	14	8.44×10−9
GO:0051707~response to other organism	10	4.98×10−7
GO:0043207~response to external biotic stimulus	10	4.98×10−7
GO:0009607~response to biotic stimulus	10	7.70×10−7
GO:0009617~response to bacterium	8	4.25×10−6
GO:0006952~defense response	11	8.94×10−6
GO:0009605~response to external stimulus	12	2.05×10−5
GO:0032496~response to lipopolysaccharide	6	3.88×10−5
GO:0071222~cellular response to lipopolysaccharide	5	3.97×10−5
GO:0071219~c10llular response to molecule of bacterial origin	5	4.64×10−5
C, Module 3
Term	Gene count	P-value
GO:0002684~positive regulation of immune system process	18	1.73×10−16
GO:0006955~immune response	20	6.16×10−16
GO:0002682~regulation of immune system process	19	2.53×10−15
GO:0048584~positive regulation of response to stimulus	20	9.26×10−14
GO:0050778~positive regulation of immune response	14	1.40×10−12
GO:0050776~regulation of immune response	15	2.33×10−12
GO:0016337~single organismal cell-cell adhesion	14	2.92×10−12
GO:0002252~immune effector process	14	3.91×10−12
GO:0098602~single organism cell adhesion	14	7.22×10−12
GO:0070489~T cell aggregation	12	8.22×10−12

[i] GO, Gene Ontology.

TCIRG1 expression is associated with the level of immune cell infiltration in GBM

Tumor purity and infiltrating lymphocytes are important predictors of prognosis in patients with glioma (26). Therefore, whether TCIRG1 expression was correlated with immune infiltration levels in GBM was investigated. The correlations among TCIRG1 expression and the levels of immune infiltration were determined according to data from the TIMER database (15). The results showed that TCIRG1 expression was negatively correlated with tumor purity (coefficient=−0.408, P=3.39×10−18) and weakly significantly corelated with CD8+ T cell infiltration (coefficient=−0.285, P=2.96×10−9) but positively correlated with the CD4+ T cell (coefficient=0.122, P=0.0127) and dendritic cell infiltration levels (coefficient=0.572, P=1.15×10−37) (Fig. 6A). The correlation coefficients between TCIRG1 expression and CD8+ T cells (coefficient=−0.285) and CD4+ T cells (coefficient=0.122) infiltration were weak. No statistical significance was found for B cells, macrophages or neutrophils. Although some correlation coefficients were small, the results still indicated that the TCIRG1 expression level is a potentially valuable factor for the immune infiltration of GBM, particularly that of CD4+ T cells and dendritic cells. Immune and stromal scores calculated by the ESTIMATE algorithm in GBM can predict distinct molecular subtypes and overall survival (22), and the relationship between the immune and stromal scores and TCIRG1 expression was further investigated. The results indicated that there was a significant difference in immune and stromal scores of patients with GBM with high TCIRG1 compared with low TCIRG1 expression (Fig. 6B-D). These results showed the potentially important role of TCIRG1 expression in the local immune response in GBM.

Figure 6. Correlation between TCIRG1 expression and immune cell infiltration and immune scores in GBM. (A) TCIRG1 expression was negatively significantly correlated with tumor purity and CD8+ T cell infiltration but positively significantly correlated with CD4+ T cell and dendritic cell infiltration levels. No statistical significance was found in B cells, macrophages or neutrophils. Increased TCIRG1 expression predicted (B) high stromal scores (B), (C) high immune scores and (D) high ESTIMATE scores in GBM samples. ****P<0.0001. TCIRG1, T cell immune regulator 1; cor, correlation.

Association between TCIRG1 expression and distinct immune markers

To clarify the TCIRG1-related immunological processes in GBM and their association with several types of infiltrating immune cells, a list of immune markers, including CD8+ T, T and B cells, monocytes, tumor-associated macrophages (TAMs), M1/M2 macrophages, neutrophils, natural killer cells and DCs, were selected from the TIMER and GEPIA datasets (15,25). The roles of T cells with different functions as described by previous studies (27–29) were also examined, such as T helper (h)1, Th2, T follicular (f)h, Th17, T regulatory (reg) and exhausted T cells. Based on TIMER, after adjustment for tumor purity, the TCIRG1 expression level was significantly correlated with gene markers of tumor-associated macrophage (CD68, coefficient=0.573, P=2.4×10−13), neutrophil (ITGAM, coefficient=0.692, P=7.76×10-21), th2 (STAT6, coefficient=0.602, P=7.40×10−15; STAT5A, coefficient=0.561, P=1.01×10−12) and T regulatory cell (TGFB1, coefficient=0.548, P=4.28×10−12). in total, 36 positively correlated genes (P<0.05) of 57 immune cell gene markers were identified in GBM (Table III).

Table III. Correlation analysis of TCIRG1 expression with related gene markers in TIMER database.

Table III.

Correlation analysis of TCIRG1 expression with related gene markers in TIMER database.

A, CD8+ T cell
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
CD8A	0.18	3.05×10−02	0.09	3.00×10−1
CD8B	0.15	6.49×10−02	0.03	6.96×10−1
B, T cell
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
CD3D	0.23	3.62×10−03	0.07	4.49×10−1
CD3E	0.40	2.95×10−07	0.32	1.66×10−4
CD2	0.31	1.10×10−04	0.17	4.43×10−2
C, B cell
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
CD19	0.18	2.53×10−02	0.16	6.70×10−2
CD79A	0.11	1.62×10−01	0.07	4.30×10−1
D, Monocyte
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
CD86	0.41	2.67×10−07	0.31	1.83×10−4
CSF1R	0.53	0.00×10°	0.47	9.20×10−9
E, Tumor-associated macrophage
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
CCL2	0.37	3.37×10−06	0.24	4.81×10−3
CD68	0.61	0.00×10°	0.57	2.40×10−13
IL10	0.33	2.68×10−05	0.21	1.28×10−2
F, M1 macrophage
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
NOS2	0.09	2.71×10−01	0.13	1.39×10−1
IRF5	0.52	0.00×10°	0.43	1.19×10−7
PTGS2	0.39	8.71×10−07	0.30	3.53×10−4
G, M2 macrophage
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
CD163	0.52	0.00×10°	0.45	4.39×10−8
VSIG4	0.43	5.81×10−08	0.35	2.60×10−5
MS4A4A	0.38	1.38×10−06	0.29	5.68×10−4
H, Neutrophil
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
CEACAM8	0.08	3.33×10−01	0.01	4.15×10−1
ITGAM	0.70	0.00×10°	0.69	7.76×10−21
CCR7	0.36	7.49×10−06	0.28	1.17×10−3
I, Natural killer cell
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
KIR2DL1	0.13	1.16×10−01	0.11	2.12×10−1
KIR2DL3	0.12	1.42×10−01	0.06	5.01×10−1
KIR2DL4	0.30	1.83×10−04	0.28	8.84×10−4
KIR3DL1	0.13	9.79×10−02	0.13	1.27×10−1
KIR3DL2	0.10	2.36×10−01	0.10	2.34×10−1
KIR3DL3	0.19	1.73×10−02	0.19	2.30×10−2
KIR2DS4	0.21	9.45×10−03	0.17	4.51×10−2
J, Dendritic cell
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
HLA-DPB1	0.42	8.08×10−08	0.34	5.74×10−5
HLA-DQB1	0.22	5.50×10−03	0.16	6.99×10−2
HLA-DRA	0.33	4.18×10−05	0.20	1.89×10−2
HLA-DPA1	0.35	1.00×10−05	0.28	1.07×10−3
CD1C	0.14	7.59×10−02	−0.02	8.37×10−1
NRP1	0.55	0.00×10°	0.52	9.65×10−11
ITGAX	0.54	0.00×10°	0.49	1.19×10−9
K, T helper 1
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
TBX21	0.13	1.14×10−01	0.16	5.98×10−2
STAT4	0.24	2.55×10−03	0.15	8.08×10−2
STAT1	0.18	2.56×10−02	0.24	5.40×10−3
IFNG	0.13	1.04×10−01	0.09	3.12×10−1
TNF	0.20	1.50×10−02	0.08	3.41×10−1
L, T helper 2
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
GATA3	0.40	3.30×10−07	0.38	3.88×10−6
STAT6	0.67	0.00×10°	0.60	7.40×10−15
STAT5A	0.61	0.00×10°	0.56	1.01×10−12
IL13	−0.02	8.01×10−01	0.04	6.10×10−1
M, T follicular helper
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
BCL6	0.28	3.86×10−04	0.29	6.30×10−4
IL21	−0.12	1.29×10−01	−0.15	9.07×10−2
N, T helper 17
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
STAT3	0.33	4.30×10−05	0.37	7.78×10−6
IL17A	−0.01	9.54×10−01	−0.05	6.02×10−1
O, T regulatory cell
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
FOXP3	0.49	1.62×10−10	0.48	4.51×10−9
CCR8	0.34	2.01×10−05	0.27	1.54×10−3
STAT5B	0.03	6.85×10−01	0.15	6.90×10−2
TGFB1	0.62	0.00×10°	0.55	4.28×10−12
P, T cell exhaustion
	None		Purity
Gene markers	Cor	P-value	Cor	P-value
PDCD1	0.35	8.12×10−06	0.32	1.49×10−4
CTLA4	0.41	2.05×10−07	0.32	1.31×10−4
LAG3	0.24	2.86×10−03	0.31	2.59×10−4
HAVCR2	0.40	3.07×10−07	0.31	1.89×10−4
GZMB	0.30	1.48×10−04	0.18	3.27×10−2

[i] Cor, Spearman's correlation coefficient; None, correlation without adjustment; Purity, correlation adjusted by purity.

It was observed that the levels of expression of the majority of monocytes, TAMs, M1/M2 macrophages, dendritic cells, Th2 cells, Tregs and T cell exhaustion immunomarkers were positively correlated with TCIRG1 expression, which were similar to the results found in GEPIA, CD68 (coefficient =0.600, P=3.50×10−17), TGFB1 (coefficient=0.690, P=4.3×10−24), STAT6 (coefficient =0.710, P=3.80×10−27) and STAT5A (coefficient =0.620, P=5.80×10−19) (Table IV).

Table IV. Correlation between TCIRG1 expression with immune gene markers in GEPIA database.

Table IV.

Correlation between TCIRG1 expression with immune gene markers in GEPIA database.

		Glioblastoma		Normal
Immune cells	Gene markers	Cor	P-value	Cor	P-value
M1 macrophage	IFR5	0.56	9.40×10−15	0.50	4.30×10−8
	PTGS2	0.42	3.30×10−8	0.33	4.90×10−4
M2 macrophage	CD163	0.47	3.30×10−10	0.58	6.10×10−11
	VSIG4	0.41	7.30×10−8	0.62	2.80×10−12
	MS4A4A	0.39	2.90×10−7	0.59	2.50×10−11
Tumor-associated macrophage	CCL2	0.39	3.40×10−7	0.47	5.50×10−7
	CD68	0.60	3.50×10−17	0.44	3.10×10−6
	IL10	0.30	1.10×10−4	0.34	3.20×10−4
Monocyte	CD86	0.39	2.00×10−7	0.44	2.20×10−6
	CSF1R	0.56	7.90×10−15	0.44	2.70×10−6
Regulatory T cells	FOXP3	0.54	1.10×10−13	−0.22	2.70×10−2
	CCR8	0.32	2.40×10−5	−0.01	0.90
	TGFB1	0.69	4.30×10−24	0.82	3.50×10−26
T cell exhaustion	PDCD1	0.39	3.70×10−7	0.37	1.20×10-4
	CTLA4	0.38	4.00×10−7	0.10	0.32
	LAG3	0.29	1.80×10−4	0.21	3.50×10−2
	HAVCR2	0.43	1.20×10−8	0.47	5.20×10−7
	GZMB	0.33	2.20×10−5	0.16	9.40×10−2
Dendritic cell	HLA-DPB1	0.42	1.70×10−8	0.69	6.10×10−16
	HLA-DRA	0.32	3.70×10−5	0.63	4.40×10−13
	HLA-DPA1	0.37	1.30×10−6	0.60	1.40×10−11
	NRP1	0.57	3.20×10−15	0.30	1.80×10−3
	ITGAX	0.63	4.60×10−19	0.40	2.20×10−5
T helper 2	GATA3	0.44	3.50×10−9	0.13	0.20
	STAT6	0.72	3.80×10−27	0.35	2.80×10−4
	STAT5A	0.62	5.80×10−19	0.70	6.70×10−17

[i] Cor, Spearman's correlation coefficient.

Furthermore, the GBM cases in the GSE16011 dataset were used to validate the correlations between TCIRG1 and immune signatures using Spearman's correlation coefficient through the GlioVis online tool (17). The results showed that TCIRG1 expression was significantly correlated with most immune signatures in the GSE16011 dataset (Table SIV).

High TCIRG1 expression was related to dense DC infiltration in GBM. The expression of the DC markers HLA-DPB1, HLA-DRA, HLA-DPA1, NRP1 and ITGAX were significantly correlated with TCIRG1 expression (Table SIV). These results further demonstrated the close relationship between TCIRG1 expression and infiltration of DCs.

In addition, a correlation between TCIRG1 expression and markers of Tregs and T cell exhaustion was reported. FOXP3, CCR8, TGFB1, PDCD1, CTLA-4, LAG3, HAVCR2 and GZMB were all positively correlated with TCIRG1 expression (Table III and IV), which further supports the potential of TCIRG1 as an important regulator of the immune response and local immune tolerance in GBM. Hence, these findings suggested that TCIRG1 expression is significantly related to immune infiltration in GBM, suggesting that TCIRG1 may play an important role in immune escape in the GBM microenvironment.

Discussion

GBM is the most highly malignant primary brain tumor with the shortest 5-year OS rate (1). Tumor-related immune responses and immunotherapy have been used in the treatment of malignant cancer, and the successful application of the CTLA-4 and PD1 checkpoints inhibitor in melanoma and non-small cell lung cancer has generated interest in glioma immunotherapy (30,31). However, as a heterogeneous cancer, GBM results in different prognoses in different patients who receive the same treatment (2). Therefore, mining data on the molecular mechanisms in the GBM microenvironment remains necessary.

V-ATPase is a macromolecular complex that is overexpressed in cancer cells, and subunit ‘a’ affects the subcellular localization of V-ATPase (6,7,9,10,14). TCIRG1 is an a3 subunit of V-ATPase that is normally expressed in osteoclasts and insulin-containing secretory vesicles in pancreatic β cells; however, aberrantly overexpressed TCIRG1 in cancer cells promotes tumor metastasis and migration potential, indicating that TCIRG1 might be a specific marker for tumor malignancy (5,7,8).

In hepatocellular carcinoma, melanoma and breast cancer, overexpression of TCIRG1 has been reported to be correlated with tumor malignancy (7,9,10). Increased TCIRG1 expression has been found in highly invasive cell lines, such as MDA-MB-231 breast cancer cells (7). In contrast, in non-invasive MCF7 breast cancer cells, TCIRG1 expression is lower compared with that in MDA-MB-231 cells, and treatment with a V-ATPase inhibitor does not reduce the invasion potential of MCF7 cells, as found in MDA-MB-231 cells (7,32). In addition, TCIRG1 is aberrantly overexpressed in the highly metastatic B16-F10 mouse melanoma cell line compared with the poorly metastatic B16 mouse melanoma cell line (9). TCIRG1-knockdown or treatment with a specific inhibitor in B16-F10 cells inhibits metastasis in an in vivo mouse model (9). However, the expression profile of TCIRG1 and its function in GBM remain to be fully understood (13,14,33). In the present study, by analyzing transcriptional data in public datasets, it was reported that TCIRG1 was highly expressed in several types of malignant tumors, which was consistent with previous studies (6,10). However, in some datasets from the Oncomine database, the results showed slight differences, as TCIRG1 expression was lower in tumors compared with in adjacent normal tissues. Discrepancies in TCIRG1 expression may result from data assembly methods and fundamental mechanisms applicable to distinct biological characteristics.

In malignant tumors of the central nervous system, the highest TCIRG1 mRNA expression was found in GBM tissues compared with LGG and normal brain tissues. Protein level validation was performed using WB analysis, and although the number of clinical samples was limited, a significant difference between normal brain tissues and glioma tissues was observed. Upon further investigation of the effect of TCIRG1 expression on GBM subtypes based on Verhaak methods (3), it was reported that TCIRG1 was highly expressed in the mesenchymal subtype, which is the most malignant subtype and has the worst prognosis of all GBM subtypes (3,4). Therefore, the present study offers novel insight into the function of TCIRG1 in glioma malignancy evaluation and highlights the potential of TCIRG1 as a biomarker of mesenchymal GBM.

Kaplan-Meier curve analyses from two independent datasets (TCGA-GBM and GSE16011) demonstrated the prognostic value of TCIRG1 expression in GBM. Elevated TCIRG1 expression significantly reduced OS in patients with GBM from both cohorts. Notably, after adjusting for IDH-1 status, age, sex, mesenchymal subtype, chemotherapy and radiation therapy information, TCIRG1 still played an important role in the prognosis of GBM. Therefore, the present study provides insight into the potential use of TCIRG1 as a prognostic marker in GBM, and TCIRG1 may become a novel marker in targeted therapy.

To further explore the biological function of TCIRG1 in GBM, genes were selected that showed positive correlations with TCIRG1 in the TCGA-GBM and GSE16011 datasets by performing Pearson's correlation coefficient analysis. The GO category analysis revealed that these genes were enriched in the ‘immune response’, ‘immune effector process’, ‘defense response’, ‘regulation of the immune system process’ and ‘positive regulation of the immune system process’. However, this is potentially due to the tumorigenesis function of TCIRG1 (7,9,10), and the negatively correlated genes were mainly enriched in physiological functions. Similar results have also been reported by Wang et al (18), in which programmed cell death 1 ligand 2 is a glioma malignancy-related gene, and its positively correlated genes are enriched in the immune response. In contrast, the negatively correlated genes mainly participate in physiological processes, such as nervous system development and neuron development. The present study demonstrated that TCIRG1 was involved in malignant glioma phenotypes and might participate in processes of the host immune system in glioma. Based on the PPI network analysis, the biological function of key modules also supported the aforementioned results.

Another important finding in the present study was that TCIRG1 expression was correlated with immune infiltration levels in GBM. As GBM has long been recognized as an immunosuppressive neoplasm that is characterized by the activation of various immune escape mechanisms (34), and to clarify the importance of TCIRG1 in the immune response, the TIMER database (15) and ESTIMATE algorithm (21) were used to compare the expression of this gene with immune cell infiltration and local immune scores. The results indicated that TCIRG1 expression was positively correlated with CD4+ T cells, DC infiltration and high immune scores but negatively correlated with tumor purity and CD8+ T cell infiltration. A correlation between the genetic markers of Tregs as well as T cell exhaustion and TCIRG1 expression was found in the TIMER database. For other immune cells, TCIRG1 expression was positively correlated with markers of M2 macrophages, monocytes and TAMs and some markers in M1 macrophages, Th2 cells, DCs, natural killer cells and neutrophils. A less significant association was reported between TCIRG1 expression and Th1, Tfh and Th17 cells. However, the mechanism of TCIRG1 expression involved in regulating the GBM local immune response remains to be elucidated.

Tumor cells can influence immune cells by producing lactic acid to change the local pH in the tumor microenvironment (11). The main function of V-ATPase is to transport protons into intracellular compartments, promote the acidification of endosomes and lysosomes and excrete intracellular protons into the extracellular space, thus maintaining H+ homeostasis (5). TCIRG1 regulates V-ATPase localization and functions in mammalian cells, and promotes acidification and maturation of cytotoxic granules in cytotoxic T lymphocytes (12). The effects of V-ATPase on regulating DC maturation, antigen presentation as well as the regulation of tumor-associated macrophages have also been reported (8,35,36). Therefore, it is reasonable to speculate that local immune infiltration in GBM may affect TCIRG1-mediated activities in the tumor microenvironment. The present results indicated the novel function of TCIRG1 in the immune response in addition to its function in V-ATPase regulation. Since immune regression in GBM is an important poor prognostic factor (30), further studies on TCIRG1 regulation of the local immune response in GBM are required.

Despite the fact that the present study used transcriptional data across different databases, several limitations still exist. First, the use of different microarrays and sequencing data from public databases inevitably introduced systematic bias (28), and the bias could be found in Oncomine data mining in the present study. Second, as a potential prognostic and diagnostic biomarker in GBM, the molecular mechanisms underlying the function of TCIRG1 in survival of patients with GBM remain uncertain, and the results did not conclusively demonstrate that TCIRG1 affects patient survival through immune infiltration. The roles of TCIRG1 in glioma cell differentiation, proliferation, invasion and migration have not been explored. Third, the present study only conducted a bioinformatics analysis of TCIRG1 expression and immune cell infiltration and patient survival in distinct public databases, but functional experiments were not performed. Further testing using in vivo and in vitro models, such as the effect of TCIRG1 on glioma cell invasion and proliferation and its effect on chemotherapy, will therefore be necessary before determining whether TCIRG1 inhibition could be an effective method for treating GBM.

In conclusion, TCIRG1 expression is associated with glioma malignancy and may be a marker of unfavorable prognosis in patients with GBM, and it could be regarded as a prognostic biomarker and an indicator of immune infiltration in GBM.

Supplementary Material

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The current study work was funded by the Liaoning Province Special Professor Project (grant no. 3110517003).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

CQ and SO were responsible for conception and design of the study. CQ acquired the data, performed the statistical analyses and wrote of the manuscript. LL, JH, GW and JL participated in data collection and helped analyze the data. SO also analyzed the data. All authors contributed to drafting and revising the article, gave final approval of the version to be published and agree to be accountable for all aspects of the work.

Ethics approval and consent to participate

The study was approved by The Ethics Committee of the First Hospital of China Medical University (Shenyang, China). All patients provided written informed consent.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References