The analysis procedure of this study is shown in Figure 1 . To estimate the fraction of tumor-infiltrating immune cells (TICs) and the amount of immune or matrix components in the KIRC samples, RNA sequencing (RNA-seq) data from 551 patients and corresponding clinical data of 537 cases were downloaded from The Cancer Genome Atlas (TCGA) database. The baseline data of enrolled cases are shown in Table 1 .

Using CIBERSORT and ESTIMATE calculation methods, the scores of the two main fractions of TICs were calculated. Subsequently, DEGs were drawn out by performing an overlap between Cox regression analysis and protein–protein interaction (PPI) network. AQP9 was identified as a latent predictor through the process. In subsequent analyses, the prognostic value of AQP9 was verified by survival analysis, differential expression, etc. Then, gene set enrichment analysis (GSEA) and correlation with TICs were performed to explore the potential mechanism. Finally, we conducted the preliminary verification through multiple immunohistochemistry (mIHC) and real-time quantitative PCR (RT-qPCR).

Kaplan–Meier survival analysis was adopted for the calculation of ImmuneScore, StromalScore, and ESTIMATEScore in order to figure out the relationship between estimated ratio and patients’ survival duration. Higher estimated scores indicate a larger number of corresponding components in the TME. ESTIMATEScore is the sum of the two ingredients and represents the comprehensive percentage of all immune-related ingredients.

As shown in Figure 2A , the immune components present a negative correlation with the overall survival rate of KIRC patients. However, the other components had no significant correlation with the overall survival rate ( Figures 2B, C ). These results implied that the immune components in TME were more suitable for indicating the prognosis of KIRC patients.

We analyzed the clinical information of KIRC patients in order to clarify the relationship between the percentage of immune and stromal components and pathological properties. As shown in Figure 2 , ImmuneScore was positively correlated with M and N classification of TNM stage ( Figure 2D ). StromalScore was not significantly correlated with TNM stage ( Figure 2E ), whereas ESTIMATEScore was strongly related with increasing N stage ( Figure 2F ; p = 0.042). It can be inferred from these results that the elevated proportion of the immune component was strongly associated with KIRC progression such as invasion and metastasis of tumor.

In order to elucidate the clear changes in gene distribution related to two immune-related ingredients in TME, the DEGs were divided into high and low score groups, and comparative analysis was performed subsequently. The result of intersection is shown in the Venn diagram, and it turns out that a total of 218 upregulated genes were shared by the ImmuneScore and StromalScore, as well as 68 downregulated genes shared by the low score group. These DEGs (a total of 286 genes) may be decisive factors of TME status ( Figures 3A–C ). It is shown by the results of Gene Ontology (GO) enrichment analysis that almost all DEGs were enriched in the functions related to immunity, such as humoral immunity and complement system ( Figures 3D–F ). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis also showed enrichment for primary immunodeficiency, cytokine–cytokine receptor interactions, and hematopoietic cell profiles ( Figures 3G–I ). Thus, the overall function of DEGs appeared to be related to immune-associated activities, suggesting that the presence of immune factors was a distinguishing part of TME in KIRC.

Further explorations focus on the detailed mechanisms of those genes in KIRC. PPI network was performed using Cytoscape software [National Institute of General Medical Sciences (NIGMS)] based on the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database. Interaction results of 286 genes are shown in Figure 4A , and the top 30 genes with the largest number of nodes in PPI network have been listed in the bar plot ( Figure 4B ). Uni-COX regression was subsequently performed to find out which gene is more important among the 286 DEGs ( Figure 4C ). Then, an overlap between the PPI network and uni-Cox regression was performed, and only four genes (AQP9, TNFRSF13B, IGLL5, HSD11B1) were presented eventually ( Figures 4A–D ).

Among the four candidates, TNF Superfamily Member 13b TNFRSF13B is known as a member of the TNF Superfamily, which proved to induce apoptosis by interacting with other members of TNF receptor family (19, 20). Also, studies have suggested that Immunoglobulin Lambda Like Polypeptide 5 IGLL5 may be involved in the regulation of TME (21). However, we believe such influence of Immunoglobulin Lambda Like Polypeptide 5 comes from a deeper layer, since it encodes the lambda polypeptides of immunoglobulin, which are often extensive and nonspecific. Therefore, we believe that this kind of connection is meaningless. As for Hydroxysteroid 11-Beta Dehydrogenase 1 HSD11B1, its translation product is a catalyst for the generation of corticosteroids. Current research believes that its function is mainly concentrated in endocrine (22, 23).

AQP9, which is an aquaporin, is believed to mainly focus on metabolism and white blood cell functions (18). The changes in the number of CD8+ T cells in the immune microenvironment are particularly related to lipid metabolism, and AQP9 happens to be connected with the regulation of this aspect (24, 25). Therefore, we chose AQP9 for further research.

In this study, all KIRC samples were divided into two groups with high or low AQP9 expression. The conclusion that the survival duration of KIRC patients with high AQP9 expression is shorter than that in the low expression group can be obtained from survival analysis ( Figure 4G ). Furthermore, it can be inferred from the results of Wilcoxon rank sum test that the expression of AQP9 in tumor samples is significantly higher than that in normal samples ( Figure 4E ). Similar results were observed in the paired analysis between normal and tumor tissue from the same patient ( Figure 4F ). The above results clearly show that the expression of AQP9 in TME was negatively correlated with survival outcomes. In particular, the expression of AQP9 increased with the progression of the TNM phase ( Figures 4H–K ).

Furthermore, we downloaded the clinical data ( Table 1 ) and RNA-seq data of 90 KIRC patients from the International Cancer Genome Consortium (ICGC) database (including tissues of 90 primary tumors and 46 adjacent tissues) in order to verify the results above. Consistent with the results of TCGA, AQP9 was significantly higher in tumor tissues (p < 0.001) ( Figure 5A ), and higher AQP9 levels were associated with shorter overall survival duration (p = 0.044) ( Figure 5B ).

GSEA results are applicable to both the high and the low expression groups because the AQP9 level is positively correlated with survival and TNM stage. As shown in Figure 5F , the genes in the high expression group of AQP9 present an association with immune-related activities [such as P53, Janus kinase (JAK)/signal transducer and activator of transcription (STAT), nod-like receptor signaling pathway, cytokine and its receptor interaction]. As for the low-AQP9 group, an apparent enrichment of these genes can be observed in fatty acid metabolism pathways ( Figure 5G ). Spliceosome is a complex composed of ribonucleoprotein (RNP) and protein, which is involved in the occurrence and development of various diseases by functioning as stimulus for the splicing of pre-mRNA (26). On the one hand, it has been reported that the spliceosome is involved in the regulation of some key genes in lipid metabolism. For instance, removal of exon 9 in fatty acid synthase (FASN) by splicing significantly reduced the production of the synthesis of fatty acid (27). On the other hand, as an aquaporin, AQP9 is indeed very likely to participate in the transport of fatty acids. Some previous studies have emphasized the relationship between fatty acid metabolism and polarization of TME (28). Therefore, it can be drawn from these results that AQP9 may be a potential indicator of TME status.

Immunotherapy, which is used in cancer treatment to block the immune checkpoints such as programmed cell death-1 (PD-1) and programmed cell death-ligand 1 (PD-L1), has shown promising outcomes in advanced cancers. However, it is not effective in all patients. In the subsequent analysis, we performed correlation analysis between AQP9 and PD-1, PD-L1 to clarify if AQP9 is able to predict the benefits of immunotherapy for patients. Furthermore, the patients who received anti-PD-L1 immunotherapy in the IMvigor210 cohort were assigned to response or non-response group, and expressions of AQP9 between the two groups were compared. However, none of these analyses indicated that AQP9 could be correlated to response to immunotherapy ( Figures 5C–E ).

The proportion of tumor-infiltrating immune-related subgroups was obtained by using the CIBERSORT algorithm. Therefore, the relationship between AQP9 and the immune microenvironment can be further derived, and 21 immune cell profiles were also constructed in KIRC samples. Among them, there were 11 kinds of immune cells presenting obviously different expressions between tumor and adjacent tissues ( Figure 6A ). Then, immune cell correlation was performed, and a total of nine immune cells were shown to be related to the expression of AQP9. After the intersection between differentially expressed immune cells and related cells, nine immune cells related to the TME were finally obtained ( Figure 6B ). Five of the immune cells were positively correlated with AQP9 expression, including regulatory T cells (Tregs), M2 macrophages, M0 macrophages, CD4 T cells, and neutrophils, while resting mast cells, resting natural killer (NK) cells, and CD8 T cells were negatively correlated with AQP9 levels. These results served as strong evidence to further confirm that AQP9 was involved in the regulation of TME immune reactivity ( Figures 6C–K ).

In order to further verify the potential correlation between the AQP9 level and tumor stage mentioned above, we applied this on six KIRC samples of different stages and a para-cancerous sample from the sample bank ( Table 2 ). Total RNA expression was extracted from these samples for RT-qPCR. As shown in Figure 7A , AQP9 expression was relatively higher in high-stage KIRC specimens. In the subsequent analysis, a simple linear regression was performed on the AQP9 levels of different stages, and the results verified that the AQP9 expression level was positively correlated with the tumor stage (R² = 0.5217, p < 0.001) ( Figure 7B ).

Since the potential connection between AQP9 and immune cells in the TME may be the key to its function, mIHC was performed on pathological specimens of renal clear cell carcinoma in order to further verify the possible relevance in clinical samples. It was observed in the tumor center that when the tumor tissue presented a high AQP9 level, the infiltration of CD8+ T cells was significantly blocked. Conversely, the recruitment of CD8+ T cells was relatively active in areas with low AQP9 level ( Figure 7B ). In addition, consistent with the above results, the quantity of M2 TAMs (labeled by CD163) and AQP9 showed obvious consistency. Therefore, the active AQP9 expression may be deeply involved in the polarization process of M2 TAMs ( Figures 7C, D ). These results preliminarily verified our inferences in the bioinformatics analysis that AQP9 was indeed related to M2 TAM levels and CD8+ T-cell recruitment.

All statistical analyses and graphing were performed by the R software (version 4.0.2) or GraphPad Prism (version 8.0.1). The Student’s t-test was used to compare the expressive difference of genes. The Spearman correlation analysis was used to evaluate the relationships between AQP9 level and the expression of various immune-related cells. Additionally, the survival analysis was based on the Kaplan–Meier method. A p-value <0.05 was considered significant.

RNA-seq data of 551 KIRC cases (normal samples, 54 cases; tumor samples, 497 cases) and corresponding data on clinical characteristics of 537 cases were downloaded from TCGA database, too (https://portal.gdc.cancer.gov/).

Clinical data ( Table 1 ) and RNA-seq data of 90 KIRC patients were downloaded from the ICGC database (including tissues of 90 primary tumors and 46 adjacent tissues), and project RECA-EU samples are limited to primary tumor (https://dcc.icgc.org/).

The ratio of the immune matrix components of each sample in TME has been expressed in the form of three scores: ImmuneScore, StromalScore, and ESTIMATEScore, which were achieved by the ESTIMATE algorithm (R language version 4.0.2). The higher the scores are, the larger the proportions of the corresponding components are in TME.

The Survminer package in R language is used to perform survival analysis. Among 497 cases, 458 tumor samples had detailed survival information and were capable of being included in the survival analysis. The survival curve was drawn by the Kaplan–Meier method, and the log rank was used as the statistical significance test. p < 0.05 was considered significant.

According to the median of ImmuneScore and StromalScore, 497 tumor samples were divided into high and low groups. Difference analysis is performed using the limma package, and DEG is generated by comparing high- and low-score samples. DEGs with fold change larger than 1 after transformation of log2 (high-score group/low-score group) and false discovery rate (FDR) <0.05 were considered significant.

GO and KEGG enrichment analyses using DEGs were performed by R language with the aid of clusterProfiler, enrichplot, and ggplot2 packages. p-values <0.05 were considered significantly enriched.

Heatmaps of DEGs were produced by R language with package pheatmap.

The clinical data corresponding to the KIRC sample were downloaded from TCGA. Wilcoxon rank sum or Kruskal–Wallis rank sum test is used as significance test according to the number of clinical stages compared in the R software.

Cytoscape version 3.8.1 was used to rebuild the PPI network that was constructed from the STRING database. Only nodes with a confidence level greater than 0.95 were included in the construction of the network.

R with package survival was used for univariate Cox regression. The order of p-values in univariate Cox is the main basis for gene sequencing.

GSEA was performed using the software gsea-4.1.0 (downloaded from Broad Institute), and Hallmark was used as target sets (downloaded from Molecular Signatures Database). All samples were included, and gene sets with NOM p < 0.05 were considered significant.

Correlation plots of AQP9 with the target checkpoint of immune therapy (PD-1 and PD-L1) were drawn according to the data from TCGA referred above.

IMvigor210 was performed to determine the predictive value of the AQP9. The IMvigor210 dataset was downloaded from http://research-pub.gene.com/IMvigor210CoreBiologies. A total of 298 urothelial cancer cases with complete clinical information were enrolled.

The CIBERSORT calculation method was adopted to estimate the TIC abundance distribution in all tumor samples, followed by quality filtering, and only samples with p < 0.05 were included in the analysis.

The datasets are available in TCGA database (https://portal.gdc.cancer.gov).

Total RNA was extracted using with EZNA Total RNA Kit (Omega Bio-tek) from six tumors and one para-carcinoma normal sample from the Sample Bank of Zhongda Hospital Affiliated to Southeast University. Primers were diluted in ddH₂O with SYBR Green PCR Master Mix (Applied Biosystems, Japan). Running platform is StepOne Plus Real-time PCR system (Applied Biosystems, Foster City, CA, USA). Transcriptional expression was determined as the fold change of AQP9 relative to Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH). Primer sequences for AQP9 and GAPDH are as follows:

The relative expression of AQP9 mRNA in KIRC was represented as 2^-ΔΔCt  = 2^-  (ΔCt (AQP9tumor) −  ΔCt (AQP9pan) , and a bar graph was drawn accordingly.

KIRC tissues were stained by brightfield and mIHC. AQP9 rabbit antibody was purchased from Abcam, catalog number ab8428. We used 4-μm tissue sections for both brightfield and fluorescence multiplex IHC. For brightfield, tissue sections were dewaxed and incubated in an autoclave for 5 min at 121°C in Tris-EDTA pH 7.8 antigen retrieval solution prior to blocking of endogenous peroxidase and incubation of the antibody. For multiplex fluorescence IHC, antibody staining was performed in subsequent steps including peroxidase blocking, application of the primary antibody, detection with a secondary horseradish peroxidase (HRP)-conjugated antibody, fluorescence dye detection, and microwave treatment to remove the bound antibodies. Three circulations to stain three different antibodies were performed in each experiment. Slides were subsequently counterstained with diamidino-2-phenylindole (DAPI) and mounted in antifade solution.