Expression profile data, phenotype data, and mutation data of GBM were downloaded from the Xena¹ database. The GBM data packet named “01A,” which including 114 samples (Table 1), was selected. The data itself has been log₂(data + 1). The information uses FPKM (Fragments Per Kilobase of transcript per Million mapped reads) without the normalization between arrays of the limma package. Twenty immune checkpoints and 122 immunomodulator-related genes were obtained from the document “Siglec15 shapes a non-inflamed tumor microenvironment and predict the molecular subtype in bladder cancer.” GBM data were obtained from CGGA (Chinese Glioma Genome Atlas) database. Data samples from two different data packets contain 693 cancer data and 325 cancer data, respectively. The number of GBM and rGBM in data packet of 693 data samples is 209, and the number of GBM and rGBM in data packet of 325 samples is 109. The raw data is RSEM, which has been converted to log2(data + 1).

The Hmisc package of R was used to calculate the correlation analysis of 20 immunological checkpoints in the data of the gene AQP4 and TCGA. Spearman’s technique was utilized for correlation analysis.

The R’s CIBERSORT was applied to assess the proportion of infiltrating 22 immune cells in cancer samples. The correlation between the gene AQP4 and the 22 immune cells implanted ratio was calculated. The same analysis was performed in low grade glioma (LGG), GBM, LGG + GBM, lung squamous cell carcinoma, lung adenocarcinoma, and breast cancer.

The limma package in R was used to do the difference analysis. The median expression of the AQP4 gene served as the grouping criterion; individuals with expression levels higher than the median were classified into high expression groups, while those with expression levels lower than the median were divided into low expression groups. The screening criteria for differential genes were P-value < 0.05, |log2FC| > 0.5853. The differential genes were intersected with 122 immunomodulator-related genes to obtain differential immunomodulatory genes, plotted by R’s heatmap package.

Immune circulation data for GBM were obtained from the TIP (Trip Medical Database) database and plotted using the ggplot2 package for R.

The differential genes of AQP4 high and low groups were analyzed for functional enrichment by the cluster profile package of R.

Heatmaps of tumor-associated immune cell effector genes were created using the R package pheatmap and data from the literature “Siglec15 shapes a non-inflamed tumor microenvironment and predict the molecular subtype in bladder cancer.”

The patients were grouped according to clinical characteristics, including age, sex, immunophenotype, molecular type, IDH mutation, and MGMT, and the expression of the AQP4 gene in different groups was compared.

Immune scores, stromal scores, and tumor purity were calculated through the estimate package for R. The R package maftools were used to calculate TMB. The R package Hmisc was then used to calculate the connection between the expression value of gene AQP4 and TMB, immunological score, stromal score, and tumor purity.

The R’s GSVA package was used to perform the pathway enrichment scores of AQP4 gene high and low groups, and the AQP4 gene expression in different immune groups in IMvigor210 (bladder cancer) and GSE91061 (melanoma) were analyzed. GSE126045 (NSCLC) was excluded from data analysis due to a lack of data.

The number of AQP4 mutations is low. Copy number variation and methylation data were obtained from the Xena database.²

The GenVisR package for R displayed high and low-grouped mutations in AQP4 gene expression.

Screening of cancer-type therapeutic targets was performed using data from the TTD (Therapeutic Target Database) database.³

As possible factors, the intersection of differential genes in high and low AQP4 gene groups, differential genes in high and low immunological score groups, and differential genes in high and low matrix score groups was used. The R package used is limma, and the filtering conditions are P-value < 0.05, |log2FC| > 0.585.

The intersection of the above differential genes was subjected to single factor cox analysis by R survival. For subsequent analysis, genes with P-values less than 0.05 were screened as prognostic factors.

Lasso analysis was performed by R’s glmnet package, followed by multivariate cox regression to obtain model scores.

The survival package in R was used to perform KM (Kaplan-Meier) survival analysis on the models. The survival ROC package in R was used to examine the AUC value of the model score and the different analyses of the genes in the high and low model score groups.

The above analysis was carried out using ICGC (The International Cancer Genome Consortium) data.

Combined with clinical information, univariate and multivariate cox regression analysis verified that risk scores were independent prognostic factors (binary classification method, at least stage, age, and gender were included). Univariate and multivariate separate predictive analyses were performed on risk scores.

The R’s rms package, foreign, and survival packages were used to construct a nomogram for clinical features and risk scores.

The human tissue microarrays were purchased from Outdo Biotech Company (Shanghai, China). A total of 77 glioma tissue samples and three normal samples were included. Slides were deparaffinized and rehydrated and were then immersed in target retrieval solution (pH 6) and boiled at medium heat three times for 10 min each in a microwave. After the slides were blocked with 3% BSA, the sections were incubated with primary antibody against AQP4 (Santa Cruz Biotechnology, Dallas, TX, USA) followed by HRP-labeled anti-rabbit IgG secondary antibody. The specimens were counterstained with hematoxylin. The negative control was obtained by replacing the primary antibody with a regular rabbit IgG. Target-positive cells were counted in 3-4 different fields and imaged using an Olympus microscope, and the immunoreactions were evaluated independently by two pathologists blinded to the clinicopathologic information in order to ensure unbiased assessment of tissue morphology. The antibody staining intensity was classified as follows: no staining, 0; weak staining, 1; moderate staining, 2; and strong staining, 3. A five-point scale was used to classify the percentage of cells stained: 0 (no positive cells), 1 (<25% positive cells), 2 (25–50% positive cells), 3 (50–75% positive cells), and 4 (>75% positive cells). The score for each tissue was calculated by multiplying the intensity index by the percentage index, yielding a score between 0 and 12. The median AQP4 score was employed to determine the cutoff value. Tumors with AQP4 scores less than or equal to the median were designated as having “low expression,” whereas those with scores greater than the median were designated as having “high expression.”

The Spearman test was utilized in every correlation analysis. The ns denotes that there is no significant difference, * means p-value < 0.05, ^** means p-value < 0.01, ^*** means p-value < 0.001, ^**** means p-value < 0.0001.

The expression of AQP4 protein in tumor samples from 77 glioma patients was evaluated by immunohistochemical (IHC) staining of a tissue microarray. The expression of AQP4 protein in each IHC sample was classified as either high (score > 6) or low (score ≤ 6), using the median of the IHC scores as the cutoff value. The correlations between AQP4 protein expression and the clinicopathologic variables of 77 glioma patients are shown in Figure 1A. The results indicated that low AQP4 expression was significantly correlated with lower clinical grade glioma (grade 1/2), while high AQP4 expression predicted higher clinical grade glioma (grade 3/4) in patients (Figure 1A) (X ² = 12.434, p < 0.001). In addition, younger age was significantly correlated with low expression of AQP4 (X ² = 4.812, p = 0.028), and female patients tended to have lower expression of AQP4 than male patients (X ² = 4.434, p = 0.035). More importantly, we also found that AQP4 was highly expressed in glioma tissues and that the expression level of AQP4 was negatively correlated with the glioma grade (Figures 1B,C) (Spearman’s rank correlation rs = 0.297, p = 0.013). Kaplan-Meier OS analyses revealed that the low AQP4 expression group (IHC staining indicating strong and moderate expression) had a better prognosis than the high AQP4 group (IHC staining indicating weak and negative expression) (p < 0.01) (Figure 1D). Overall, these data demonstrated the prognostic significance of AQP4 in glioma.

In the correlation analysis between the gene AQP4 and 20 immune checkpoints, only five immune checkpoints were significantly correlated. CD276, LAG3, VTCN1, and PVR were negatively correlated with the gene AQP4, and BTLA was negatively correlated with the gene AQP4 (Figures 2A,B and Supplementary Tables 1–3), and the correlations were all statistically significant. We also looked into the relationship between AQP4 expression and the number of immune invading cells (Figures 2C,D). The AQP4 gene was strongly connected with five immune cells, including B cells, out of 22 immune invading cells. Macrophages are immature. M0, plasma cells, T-cells, CD4 memory resting, and monocytes are some terms used to describe M0, plasma cells, T-cells, and monocytes (Supplementary Tables 4–6). The tumor types we included in the study mainly included LGG, GBM, LGG + GBM, lung squamous cell carcinoma, lung adenocarcinoma, and breast cancer. In Macrophages M0, all cancers except LGG were significantly associated with the gene AQP4. Only breast cancer had no significant correlation with the gene AQP4 in monocytes and T-cells, memory resting. In contrast, LGG, GBM, LGG + GBM, lung squamous cell carcinoma, and lung adenocarcinoma all had significant correlations (Supplementary Tables 7, 8).

Initially, we looked at how AQP4 gene expression affects immune cells’ ability to present and process antigens. A total of 122 genes were collected from the literature, 32 of which were substantially different in the AQP4 high and low expression groups, with the majority of the differential genes being low in the AQP4 high expression group. This indicated that antigen presentation and processing were attenuated in the high AQP4 group (Figure 3A and Supplementary Table 9). We further analyzed the anti-cancer immune process activity differences between high and low AQP4 gene expression groups (Figure 3B). In the high and low AQP4 expression groups, only the following five processes were significantly different between the high and low groups, including the release of cancer cell antigens, tumor antigen presentation, initiation and activation, eosinophil recruiting, and cancer cell clearance. The other 18 immunological mechanisms, on the other hand, did not differ significantly between the high and low groups (Figure 3B and Supplementary Table 10). Following that, we examined the relationships between the AQP4 gene, immunotherapy-related prediction pathway enrichment scores, and anti-cancer immune process activity (Figure 3C).

We obtained 26 immunotherapy-related predictive pathways from the literature. The three groups with the strongest correlation in the anti-cancer immune process activity correlation analysis are Step 3 Priming and activation, Step 7 killing cancer cells, and Step 4 eosinophil recruiting, the correlations are −0.373, −0.288, and −0.186. The enrichment score was obtained by GSVA among the immunotherapy-related prediction pathways. Finally, we found that 14 pathways were significantly correlated with AQP4. The three groups with the strongest correlation were FGFR3-coexpressed genes, EGFR ligands, and pyrimidine metabolism, with correlations of 0.432, 0.422, and −0.366, respectively (Figure 3C and Supplementary Table 11).

Among the biological process (BP) pathways, extracellular matrix organization, extracellular structure organization, external encapsulating structure organization, organic anion transport, and vascular process in the circulatory system have the most significant differences. The cellular component (CC) pathway showed substantial differences in the collagen-containing extracellular matrix, basolateral plasma membrane, basal plasma membrane, the basal section of the cell, and collagen trimer. In the MF pathway, the extracellular matrix structural constituents, monocarboxylic acid transmembrane transporter activity, carboxylic acid transmembrane transporter activity, organic acid transmembrane transporter activity, and growth factor binding had the most significant differences. Among the KEGG pathways, Proximal tubule bicarbonate reclamation, Bile secretion, Protein digestion and absorption, Gastric acid secretion, and AGE-RAGE signaling pathway in diabetic complications showed the most significant differences (Figures 4A–D and Supplementary Tables 12–15). The relationship between AQP4 gene levels and tumor-associated immune cell effector genes was also investigated (Figure 4E). Seven genes were significantly altered (p < 0.05) among the 35 immune cell effector genes. MARCO, SPON2, TBX21, MMP8, and SLAMF8 were among the genes with low expression in the AQP4 high expression group, while PRR5L and FGL2 were substantially expressed in the AQP4 high expression group (Figure 4E and Supplementary Table 16).

The AQP4 gene did not differ significantly in age, sex, IDH mutation, and MGMT. There were significant differences between C1 and C4 in immunophenotyping. There were significant differences in molecular typing between Classical and G-CIMP, Classical and Mesenchymal, Classical and Proneural, G-CIMP and Neural, Mesenchymal and Neural, and Neural and Proneural (Figure 5A). AQP4 gene expression was found to be substantially linked with microsatellite instability (MSI) but not with tumor mutational burden (TBI), immunological score, stromal score, or tumor purity (Figure 5B and Supplementary Table 17). In addition, we further analyzed the effect of AQP4 gene expression level on the enrichment score of immunotherapy prediction-related pathways, as well as differences in AQP4 gene expression among different immune response groups in the immune dataset (Figure 5C). According to the findings, only eight groups of immunotherapy predicted pathway enrichment scores were significantly different between the high and low expression groups of the AQP4 gene. Other immunological groups exhibited significant differences in progressive disease (PD) and partial response (PR)/complete response (CR) when GSE91061 was analyzed using the AQP4 gene (Figure 5C).

The relationship between methylation and expression of the AQP4 gene was also investigated. AQP4 has three methylation sites, three of which have a strong relationship with the gene’s expression (Figure 5D). We also analyzed the mutation status in the high and low AQP4 gene expression groups and displayed the top 20 genes with mutation rates in the AQP4 low and high expression groups (Figure 5E). In addition, we evaluated AQP4-related therapeutic targets in GBM. To extract data for analysis, we used the TDD database. The findings revealed that three genes, EGFR, ABAT, and PDGFRA, significantly differed between high and low AQP4 expression groups in GBM (Figure 5F).

As a preliminary step, we searched for possible factors in the intersection of differential genes in the AQP4 gene high and low expression groups, the immunological score high and low group, and the matrix score group. The AQP4 gene high and low expression group had 291 differential genes, the immunological score high and low group had 883 divergent genes, the matrix score group had 792 differential genes, and the intersection of the three groups had 42 difference genes (Figures 6A–C and Supplementary Tables 18–20). Subsequently, we performed a univariate cox analysis on the 42 differential genes and obtained eight prognostic factors with significant differences (Figure 6D and Supplementary Table 21). In addition, duplicate elements were eliminated using LASSO dimensionality reduction to develop a useful predictive risk scoring model. After lasso dimension reduction, we obtained a model consisting of three genes from the eight prognostic factors. Then perform multi-factor screening on these three genes to get the final scoring model constructed by the three genes (Figure 6E and Supplementary Table 22), that is, Risk score = RARRES1*(0.06121847) + SOCS3*(0.16404808) + TTYH1*(−0.07449025).

We further evaluated the predictive power of the prognostic model. The high and low groups are divided according to the median value of the risk score of the prediction model. In evaluating the predictive model from the TCGA data, there was a significant difference in the survival analysis, with an AUC value of 0.774 for the ROC curve in the third year (Figures 7A–F). Following that, we validate the dataset with an external source. The survival analysis was significantly different across the two CGGA data sets, with AUC values of 0.602 and 0.695 for the 3-year ROC curves, respectively (Figures 7G–J).

Initially, both univariate and multivariate independent prognostic assessments of risk scores revealed substantial differences (Figures 8A,B). Then, using a nomogram and clinical features, we created a clinical practice guide. In the constructed nomogram, risk score and age had a more significant impact on prognosis, and the combined AUC value of risk score and clinical features in the ROC curve was 0.777 (Figures 8C–G).