This study was approved by the Academic Committee of the Second Affiliated Hospital of Fujian Medical University, and it was performed according to the principles expressed in the Declaration of Helsinki. The study received ethics approval, and all patients gave written informed consent. Additionally, the public datasets used in this study were retrieved from the published literature.

Gene expression profiling interactive analysis (GEPIA, http://gepia.cancer-pku.cn/) is a newly developed online analytical tool, which is based on the sequencing database containing of 9,736 tumors and 8,587 normal samples from The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) programs. It provides key interactive and customizable functions including tumor/normal expression profiling and differential expression analysis according to cancer types or pathological stages, patient survival analysis, similar gene detection, correlation analysis, and dimensionality reduction analysis (Tang et al., 2017).

UALCAN (http://ualcan.path.uab.edu) is a comprehensive web portal that facilitates in‐depth analysis of TCGA gene expression data (Chandrashekar et al., 2017). In this study, we analyzed the expression of AQP genes across normal and LUAD tissues based on different tumor stages. Statistical significance was tested using Student’s t test, with the significance accepted at p < 0.05.

Total RNA was extracted from 10 paired stage I LUAD tissues and paracancerous tissues using the RNeasy Mini Kit (Qiagen, Germany) according to the manufacturer’s protocol. Then ribosomal RNA (rRNA) was removed using the RiboZero rRNA removal kit. The rRNA-depleted RNA was fragmented and reverse-transcribed. mRNA sequencing libraries were prepared using the VAHTS total RNA-seq Library Prep kit for Illumina (Vazyme NR603, China) following manufacturer’s instructions. Differential expression analysis for mRNA was performed using DESeq2 R package (https://bioconductor.org/packages/release/bioc/html/DESeq2.html). Genes were considered differentially expressed and retained for further analysis with |Log2 (fold change)| (|log2FC|) ≥1 and statistical p value ≤ 0.05. The heat map was used to display the expressed pattern of AQP family members in LUAD patients.

The total RNA of 28 paired stage I LUAD (13 men and 15 women; mean age, 64.61 ± 8.83 years) and adjacent normal tissues was obtained using TRIzol as per the established protocol. The Takara PrimeScript™ RT reagent kit (Takara, Japan) was utilized to synthetize cDNA from 1,000 ng RNA. TB Green Mixture (Takara Bio, Japan) was used to conduct RT-PCR using QuantStudio™ 5 Real-Time PCR Systems (Applied Biosystems). Every sample was assessed in triplicate. GAPDH served as control.

The relative gene expression obtained from this normalization was evaluated using 2^−△△CT methods. Sequencing data were provided in Supplementary Table S1.

IHC staining of LUAD tissues was performed in 5-μm sections. Paraffin-embedded sections were dewaxed and rehydrated in a series of alcohol to PBS. Endogenous peroxidase was then inactivated with 3% hydrogen peroxide at room temperature for 20 min. Then the slides were soaked in 0.1 mol/L citrate buffer (pH 6.0) and placed in an autoclave at 121°C for 3 min for antigen retrieval. After washing with PBS (pH 7.4) for 3 times, the sections were blocked with 1% BSA diluted in PBS at 37°C for 30 min, followed by incubation with primary antibody against AQP1 (1:100; Bioss, China) and AQP3 (1:200; Bioss, China) overnight at 4°C. Then the sections were incubated with the HRP-conjugated goat anti-mouse/rabbit antibody at room temperature for 1 h, followed by staining with the DAB until the appearance of brown color. Finally, the sections were counterstained with hematoxylin and mounted.

Kaplan–Meier analysis of four AQPs in LUAD tissues on overall survival (OS), first progression survival (FP), and post-progression survival (PPS) was performed using the online platform Kaplan–Meier plotter (https://kmplot.com/analysis/). Besides, we further displayed the OS among different tumor stages of LUAD patients using Kaplan–Meier curves and performed the log-rank (LR) test.

The cBioPortal (http://cbioportal.org) is a comprehensive web-based database that could visualize and analyze multidimensional genomic data of multiple kinds of tumors (Gao et al., 2013). Mutation frequencies, copy number variation (CNV), and the summary of the gene types in LUAD tissues were evaluated according to the online instructions of cBioPortal. The relationship between patients’ prognosis and gene mutation was analyzed using the tool of cBioPortal based on TCGA database, with a p-value < 0.05 regarded as significant.

GeneMANIA (http://genemania.org/) is an online database, a predictive server for analyzing physical interactions, co-expression, and information strength of target genes. Using this database, we analyzed the correlations between AQP superfamily molecules and their interactive genes.

Metascape is a comprehensive tool for gene annotation and enrichment analysis. In the current study, we evaluated the functions of the four AQP and their co-expression genes. The threshold value was set as 0.01, and the enrichment factor of >1.5 and a minimum count of 3 were considered important. (Zhou et al., 2019).

To explore the specific associations of AQP family genes with immune cells, we utilized the TIMER database for analysis (https://cistrome.shinyapps.io/timer/), which is a database designed for analyzing immune cell infiltrates in multiple cancers. In this study, we analyzed AQP expression levels in LUAD, and their associations with tumor purity and infiltrating immune cells including CD4+T cells, CD8+T cells, B cells, neutrophils, macrophages, and dendritic cells. Furthermore, we employed the TISIDB database (http://cis.hku.hk/TISIDB/) for evaluating whether AQP genes was related to immune infiltration therapy and subtypes in LUAD patients.

The GEPIA database was applied to investigate the expression differences of AQP genes between tumor and normal tissues in LUAD patients at the mRNA level (tumor cases = 483, normal cases = 347). According to the GEPIA database, AQP1, AQP3, AQP4, and AQP9 were significantly expressed in LUAD tissues compared to normal lung tissues (p < 0.05), but some genes including AQP0, AQP5, AQP6, AQP7, AQP8, AQP10, and AQP11 were not differentially expressed between tumor and normal tissues. In addition, we found that AQP2, AQP12A, and AQP12B were not expressed either in cancerous or in normal lung tissues (Figure 1).

Previous studies have reported that four AQP family members (AQP1, AQP 3, AQP 4, and AQP 5) were expressed in lung tissues (Warth et al., 2011; Verkman, 2012; Wittekindt and Dietl, 2019); therefore, we further explored the expression levels of these genes with different tumor stages for LUAD patients using UALCAN databases. Interestingly, all these four genes were found to be aberrantly expressed in stage I LUAD patients when compared with normal tissues. As shown in Figure 2A, AQP3 was significantly increased (p < 0.001), while the expressions of AQP1, AQP4, and AQP5 were markedly decreased in patients with stage I LUAD (p < 0.001).

To further validate the expression of these four AQP family molecules in stage I LUAD, we conducted RNA sequencing analysis in 10 pairs of stage I LUAD and paracancer tissues. The heat map clearly showed an AQP-related cluster (Figure 2B), which exhibited a higher AQP3 expression and lower AQP1 expression than matched normal tissues (p < 0.001). Furthermore, total RNA was extracted from 28 paired stage I LUAD tissues to further validate the results of RNA sequencing and UALCAN databases by RT-PCR. The results indicated that AQP3 expression was significantly increased and the expression of AQP1 was markedly decreased in LUAD tissues compared with corresponding paracancerous tissues (p < 0.05) (Figure 2C). However, the expressions of AQP4 and AQP5 were not significantly different between paired tissues. Additionally, we further performed IHC staining to determine protein expression of differential AQP1 and AQP3. The results also suggested that AQP3 was highly expressed and AQP1 was poorly expressed in stage I LUAD tissues compared with paracancerous tissues (Figure 2D).

To identify the effects of AQP genes on the progression and prognosis of LUAD, the differential expression levels of these four AQP genes were associated with LUAD patients using Kaplan–Meier plotter. The Kaplan–Meier survival curve analysis revealed that higher expression levels of AQP1/4/5 were related to longer OS. Conversely, the overexpression level of AQP3 indicated a poorer OS in all LUAD patients. Beyond that, the high level of AQP5 showed a positive effect on FP. And we found that AQP1/3/4/5 expression levels had no significant correlation with PPS (Figure 3A). Next, we further analyzed the association of the expression of the four AQPs with prognosis in different tumor stages. We found that AQP1/3/4/5 expression levels showed significant statistical differences in OS among all patients with stage I LUAD. By contrast, there was no statistical significance of AQP1/3/4/5 expression linked with OS among patients with stage II and stage III LUAD (Figure 3B).

And then we evaluated the frequency of genetic alterations of four AQP genes in LUAD patients using the cBioPortal database. As shown in Figure 4A, more than 15% of LUAD patients exhibited obvious alterations in the four AQP genes, including mutation, amplification, deep deletion, mRNA high expression, and multiple alterations. Besides, the genetic alterations of AQP1, AQP3, AQP4, and AQP5 among LUAD patients account for 2.4, 1, 2.2, and 0.8%, respectively (Figure 4B). Furthermore, the association between the alterations of AQP gene expression and LUAD prognosis was performed using the cBioPortal database. The results indicated that alterations of four AQP genes in LUAD patients were not significantly associated with OS (Figure 4C, p = 0.257) or DFS (Figure 4D, p = 0.470).

In order to investigate potential mechanisms of the four AQP family genes in LUAD, the gene–gene interaction network between AQP genes and their functionally related genes was established using GeneMANIA. The results demonstrated that 20 genes, including TRAPPC8, PLD2, GK, TP53, AQP12A, AQP11, AQP12B, RP11-407P15.2, FP325317.1, RP5-877J2.1, MIP, AQP8, SEC61A1, AQP10, AQP9, ATF1, AQP6, TRPV4, AQP2, and AQP7, were mainly associated with regulatory functions of aberrantly expressed four AQP family members in LUAD patients (Figure 5A).

Next, we used Metascape for gene ontology (GO), Kyoto Encyclopedia of Genes and Genome (KEGG), and protein–protein interaction (PPI) enrichment analyses. The data revealed the eight most enriched terms, including passive transport by aquaporins, water homeostasis, protein tetramerization, multicellular organismal water homeostasis, carbon dioxide transport, bile secretion, kidney development, and small GTPase-mediated signal transduction (Figure 5B). We also constructed a network of enriched terms colored by p-values (Figure 5C).

We first analyzed the correlation between abundance of tumor-infiltrating lymphocytes (TILs) and expression of four AQP molecules using TISIDB database. Interestingly, we found that the expression levels of AQP3/4/5 were obviously increased inflammation of the LUAD immune subtype (Figure 6A, p < 0.001).

Next, we further evaluated the association between AQP gene expression and LUAD molecular subtypes using TIMER. Our results showed that high AQP1/4 expression and low AQP3 expression were strongly related to high infiltrating abundances of macrophages, neutrophils, and dendritic cells in LUAD, while AQP5 expression was not significantly correlated with infiltration abundances of these lymphocytes. Additionally, we found that high AQP1/5 expression was notably associated with high infiltrating abundances of B cells and CD4+ cells (Figure 6B ).