Identification of key genes and validation of key gene aquaporin 1 on Wilms’ tumor metastasis

Materials

The gene expression profile data sets with accession numbers GSE57269, GSE73209 and GSE66405 were found from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) (Edgar, Domrachev & Lash, 2002). Previously, Shukrun et al. (2014) submitted chip data (GEO accession GSE57269) to Platform GPL15207, Karlsson et al. (2016) submitted chip data (GEO accession GSE73209) to Platform GPL10558, and Ludwig et al. (2016) submitted chip data (GEO accession GSE66405) to Platform GPL17077. The GSE57269 data set contains 7 Wilms’ tumor samples (Registration number: GSM1378119∼GSM1378125) and 1 normal kidney sample (Registration number: GSM1378118). The GSE73209 data set contains 32 Wilms’ tumor samples (Registration number: GSM1888569∼GSM1888600) and 6 normal kidney samples (Registration number: GSM1888603∼GSM1888608). The GSE66405 data set contains 28 Wilms’ tumor samples (Registration number: GSM1621548∼GSM1621575) and 4 normal kidney samples (Registration number: GSM1621576∼GSM1621579).

Methods

GO and KEGG pathway enrichment analysis

Enrichment analysis is a commonly used method to analyze gene sets derived from genomic studies and identify their associated biological functions and signaling pathways. For the analysis of gene function and pathway enrichment, we utilized the Database for Annotation, Visualization, and Integrated Discovery (DAVID), (https://david.ncifcrf.gov), a comprehensive web-based resource (Huang da, Sherman & Lempicki, 2009). This resource facilitated Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses in our study. Significance was determined based on an adjusted p value threshold of < 0.05.

PPI network and key genes were selected

The protein–protein interaction network (PPI) was investigated using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/) (Szklarczyk et al., 2019). The DEGs were subjected to analysis in STRING, leading to the construction of the PPI network. To ensure reliability, a minimum interaction score of 0.7, denoting high confidence, was employed as the cutoff threshold. In the PPI network, high confidence score was used as indicators to filter and determine the reliability and repeatability of protein interactions. Setting 0.7 as the minimum value was done to ensure that the resulting PPI network has a higher level of confidence and biological significance. The PPI network was visualized using Cytoscape v3.7.2 software (https://cytoscape.org/) (Shannon et al., 2003). To identify significant modules within the network, we utilized the Molecular Complex Detection (MCODE) plugin, a tool integrated with Cytoscape. The module identification process involved applying specific criteria, including a degree cutoff of 2, a node score cutoff of 0.2, a K-core value of 2, and a maximum depth of 100 (Bader & Hogue, 2003). By applying the aforementioned criteria, it is possible to filter and identify modules that exhibit high connectivity, dense interactions, and biological relevance, enabling a better understanding of protein interactions and biological functions. We employed the CytoHubba plugin in Cytoscape (Chin et al., 2014) to rank and explore crucial genes within the PPI network modules. For the top 5 key genes, we conducted enrichment analyses using DAVID, applying a significance criterion of an adjusted p value < 0.05, for both GO and KEGG enrichment analyses.

Gene regulatory network analysis

We utilized the eXpression2Kinases (X2K) web-based application (Clarke et al., 2018) (https://amp.pharm.mssm.edu/X2K/) to examine the participation of upstream cellular signaling pathways and establish the regulatory connections among transcription factors (TFs), protein kinases, and their target genes. By employing default parameters, X2K was utilized to examine the upstream regulatory network governing the gene expression regulation of the upregulated genes. The inferred network, which consisted of co-expressed differentially expressed genes, was constructed and visualized.

Cell culture and tissues

WT Wit-49 cells and human renal tubular epithelial cells (HK-2) were obtained from the Cell Bank of the Chinese Academy of Sciences. These cells were cultured in a medium supplemented with 10% fetal bovine serum (FBS) under standard conditions: 37 °C and 5% CO2. From May 2013 to May 2021, WT and matched normal renal tissues from 32 children with pathologically confirmed WT were collected from the Affiliated Hospital of Qingdao University. None of these patients had received adjuvant radiotherapy or chemotherapy prior to surgical intervention. This study was conducted in accordance with the ethical guidelines set by the Ethics Committee of the Affiliated Hospital of Qingdao University (Registration No. QYFY WZLL 27451). Prior written informed consent was obtained from the parents or legal guardians of all the participating children.

RNA extraction and quantitative real-time PCR (qRT–PCR) analysis

Total RNA from the gastrocnemius was extracted using TRIzol reagent (Invitrogen, Grand Island, NY, USA), and reverse transcription of RNA into complementary DNA (cDNA) was performed using a reverse transcription kit (TaKaRa, Shiga, Japan) and oligo (dT) primers. qRT–PCR (Ying et al., 2017) was conducted using gene-specific primers and SYBR Green qPCR mix (Thermo Fisher, Waltham, MA, USA). The primer pairs were as follows: Renin forward: 5′-CTCTACACTGCCTGTGTGTATC-3′, Renin reverse: 5′-CACTGACTGTCCCTGTTGAATA-3′. NPHS2 forward: 5′-CTGTGAGTGGCTTC-TTGTCCTC-3′, NPHS2 reverse: 5′-CCTTTG GCTCTTCCAGGAAGCA-3′. SLC12A3 forward: 5′-TGGACGACCATTTC CTACCTGG-3′, SLC12A3 reverse: 5′-CACTCGGTGAAG TTCCAGCCAT-3′. SLC12A1 forward: 5′-AGGCTCTTTCCTACGTGAGTGC-3′, SLC12A1 reverse: 5′-GCCACTGTTCTTGGTAAAGGCG-3′. AQP1 forward: 5′-CTGTGGGATTAA CCCTGCTCG-3′, AQP1 reverse: 5′-GAAG -CTCCTGGAGTTGATGTCG-3′. GAPDH was used as a reference for each sample.

Western blotting

Western blotting (WB) was conducted following a previously described protocol (Barik et al., 2014). Total protein was extracted using RIPA Lysis Buffer (Beyotime, China), and the protein concentration was determined using the Bicinchoninic Acid (BCA) Protein Assay (Beyotime, China). SDS-PAGE was performed to separate the protein samples, which were then transferred to a nitrocellulose membrane, blocked, and probed with primary antibodies, including anti-AQP1, anti-Vimentin (1:1000; Abcam, Cambridge, UK), E-Cadherin (1:1000; CST, USA), N-Cadherin (1:1000; CST, USA), and anti-β-actin (1:1000; Abcam, Cambridge, UK). The membranes were subsequently incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies (1:1000; Invitrogen, Waltham, MA, USA). Immunoreactive protein bands were visualized using an HRP chemiluminescence detection reagent (ECL; Thermo Fisher Scientific), and the blots were imaged using a ChemiDoc MP imaging system (Bio-Rad, Hercules, CA, USA). The grayscale values were quantified by ImageJ.

Immunohistochemical (IHC) staining

Paraffin sections of WT and matched normal kidney tissues (more than five cm from the edge of the cancer focus) were stained with an anti-AQP1 antibody. Before staining, the paraffin-embedded tissue sample was cut into 4 µm sections, dewaxed, rinsed in anhydrous xylene, and rehydrated by gradient ethanol for antigen repair (Al-Dhohorah et al., 2016). To block endogenous peroxidase activity, the sliced samples were exposed to 3% hydrogen peroxide, followed by three washes with phosphate-buffered saline (PBS). Subsequently, the samples were sealed with a 5% FBS-PBS solution for 15 min and then incubated overnight at 4 °C with the primary antibody. On the following day, the samples were incubated with the secondary antibody for 60 min. Finally, 3,3′-diaminobenzidine (DAB) was detected by immunohistochemical Avidin-Biotin Complex (ABC) staining. Under a light microscope, the cytoplasm or nucleus was stained yellow or brown for positive staining, and samples without staining were negative. The images were taken at 100 × and 200 × magnification. ImageJ software was used to analyze the optical density of the immunohistochemical photos. Three pathological sections of each specimen were taken to detect the cumulative integrated optical density (IOD), and the average value was calculated.

Transfection

Wit-49 cells were maintained in a 37 °C incubator with 5% CO2 until they reached approximately 70% confluence. The AQP1 cDNA (NM_198098) plasmid or a negative control plasmid was obtained from GeneChem (Shanghai, China) and transfected into Wit-49 cells using Lipofectamine 3000 (Invitrogen, CA, USA). After 24 h of incubation, the medium was replaced, and the corresponding functional assays were performed 24 h later. Protein extraction was conducted 48 h after transfection, and the overexpression of AQP1 was confirmed by Western blotting.

Wound healing to detect cell migration

Confluent monolayers of AQP1-overexpressing cells and negative control cells were subjected to a defined wound by manually removing a 300–500 µm strip of cells using a standard 200 µl pipette tip. The wounded monolayers were washed twice to eliminate nonadherent cells. The speed of wound healing was quantified by measuring the average linear migration rate of the wound edges over a 24-hour period. Cell migration was observed under a light microscope, and cell mobility was calculated using ImageJ analysis software as follows: Cell mobility = (0 h scratch area-24 h scratch area)/0 h scratch area. The average was calculated from 3 independent experiments.

Transwell assays to detect cell invasion

For the analysis of invasion ability, chambers coated with Matrigel (Corning, NY, USA, 354480) were utilized to culture AQP1-overexpressing cells and negative control cells. After incubating for 24 h, the cells that migrated through the Transwell chambers were fixed using methanol and stained with Giemsa. The number of migrated cells was quantified by examining five randomly chosen fields per well under a microscope. The average value was calculated based on data obtained from three independent experiments.

Phalloidin-iFlour 488 staining of the cytoskeleton

Cells were seeded onto cell climbing tablets and cultured for 24 h. After that, they were fixed using 4% paraformaldehyde for 20 min. To enhance cell permeability, a 0.1% Triton X-100 solution was applied for 15 min. For fluorescence staining of F-actin, coverslips were treated with Phalloidin-iFlour 488 reagent (ab176753) (1:1000; Abcam, Cambridge, UK) at room temperature for 45 min, while being protected from light. After washing the cells three times with PBS, nuclear staining was performed using DAPI in the fluorescence assay. Finally, all cells were observed using a Laser Scanning Confocal Microscope (LSCM) (Nikon, Tokyo, Japan) following mounting. Phalloidin-iFlour 488 staining of F-actin revealed the cytoskeleton, which appeared as green fluorescence, and DAPI staining was used to stain the nuclei. Then, we overlapped the two parts.

Statistical analysis

The mean and standard deviation (SD) were used to present the data. Statistical analysis was performed using SPSS version 22.0 software (SPSS, Chicago, IL, USA). T tests were conducted for statistical evaluation, and a p value < 0.05 was considered indicative of significant differences.

Identification of WT-related genes

The research flowchart of this study is shown in Fig. 1. GSE57269 contains 7 WT samples and 1 normal kidney sample, GSE73209 contains 32 WT samples and 6 normal kidney samples, and GSE66405 contains 28 WT samples and 4 normal kidney samples. A total of 5128 DEGs in WT were screened, including 1284 upregulated genes and 3844 downregulated genes (—log2fold-change (FC)— ≥1.5 and adjusted p value < 0.05) (Figs. 2A, 2B, 2C). The intersection of the three groups of DEGs contained 73 genes (co-expressed DEGs) (Fig. 2D).

GO and KEGG pathway enrichment analysis of the co-expressed DEGs

In order to gain a deeper understanding of the co-expressed DEGs, we conducted enrichment analyses for GO and KEGG pathways. The GO enrichment analysis included categories such as biological process (BP), cell component (CC), and molecular function (MF). The results revealed that the co-expressed DEGs were primarily associated with BP terms such as “sodium ion homeostasis”, “sodium ion transmembrane transport”, and “multicellular organismal water homeostasis.” In terms of CC, they were enriched in terms such as “extracellular exosome” and “apical plasma membrane”. As for MF, the DEGs were mainly involved in “protein binding” and “sodium: potassium: chloride symporter activity” (Fig. 3A). KEGG pathway analysis identified significant enrichment of the co-expressed DEGs in processes such as “Metabolic pathways”, “Proximal tubule bicarbonate reclamation”, and ”Aldosterone-regulated sodium reabsorption” (Fig. 3B). Taken together, these findings highlight that the co-expressed DEGs may be associated with homeostatic and transmembrane transport of intracellular water and ions, mainly involved in metabolism-related signaling pathway.

Key genes were selected by the PPI network

The STRING database was applied to construct a protein-protein interaction (PPI) network comprising 73 co-expressed DEGs to explore core modules and key genes, both of which could play a central role in the network. The PPI network, including 48 nodes (genes) and 93 edges (interactions), was further analyzed by two built-in algorithms within Cytoscape (Fig. 4A). MCODE identified highly interconnected protein subgraphs from PPI network, which often represent functionally relevant protein complexes. 2 modules were retrieved from the PPI network constructed using co-expressed DEGs (Figs. 4B–4C). Among them, the module 1 included seven nodes (genes) and 11 edges (interactions), and the clustering score (density multiplied by the number of members) was 3.667; the module 2 included five nodes (genes) and seven edges (interactions), and the clustering score (density multiplied by the number of members) was 3.5. To rank the important genes within the co-expressed DEGs, the CytoHubba plug-in in Cytoscape was employed based on the degree centrality (Fig. 4E). The top five key genes, identified based on their node degree, were REN, NPHS2, SLC12A3, SLC12A1, and AQP1.

X2K construction and visualization of the upstream regulatory network

We used X2K to construct and visualize the upstream regulatory network of the co-expressed DEGs (Fig. 4D). We computationally predicted the involvement of upstream cell signaling pathways and identified the regulatory associations between TFs, protein kinases, and their target genes. E2F Transcription Factor 1 (E2F1) and CAMP Responsive Element Modulator (CREM) were the most important TFs that regulated the co-expressed DEGs. In addition, Mitogen-Activated Protein Kinase 14 (MAPK14) and cyclin-dependent kinase 4 (CDK4) were identified as the most important kinases controlling the co-expressed DEGs. Together, the prediction of the upstream regulatory network controlling the co-expressed DEGs could help us to conduct in-depth mechanism research.

GO enrichment analysis of the top five key genes

The functional annotations of top five key genes were performed by DAVID, GO analysis results revealed of BP, CC and MF with p value (Table 1).

Validation of top five key genes expression

In the GEO dataset analysis, REN, NPHS2, SLC12A3, SLC12A1, AQP1 were significantly downregulated in WT tissues. To validate this finding, the expression levels of the top 5 key genes in WT Wit-49 cells and HK12 normal kidney cells, we performed RT-qPCR. The results revealed that the mRNA expression levels of top five key genes were all downregulated in the WT Wit-49 cells. No statistically significant relationship existed of REN and SLC12A1 (Fig. 5A).

To further investigate the expression levels of AQP1 in renal tubules and WT tissues. Western blotting showed that the expression of AQP1 protein in the WT Wit-49 cells was significantly lower than that in normal renal tubular epithelial HK2 cells (Fig. 5B). We used immunohistochemical (IHC) staining to detect the expression levels in normal renal tissues and in WT tissues. AQP1 was expressed in different parts of the tissues, showing a brown color, and was located in the cytoplasm and cell membrane. AQP1 was expressed in glomeruli, renal tubules, and cancerous tissues and was clearly stained in normal renal tissues. The AQP1 level in pathological sections was quantified by ImageJ software, and the average optical density (AOD) was obtained. The results showed that the normal renal tissues had a higher AOD value than the WT tissues. The protein levels of AQP1 in WT tissues, which were detected by IHC, also confirmed that the protein-level expression of AQP1 is lower in WT (Fig. 5C).

Overexpression of AQP1 suppresses the migration and invasion of WT cells in vitro

The functional evaluation of AQP1 in WT cells was performed using an overexpression plasmid, which resulted in a substantial increase in AQP1 expression in the AQP1 overexpression group (Fig. 6A). Wound healing and Transwell assays were conducted to investigate the impact of AQP1 on WT cell migration and invasiveness. The findings revealed that AQP1 overexpression significantly reduced the migratory and invasive capacity of Wit-49 cells (Fig. 6B). Notably, the scratch healing rate in the AQP1 overexpression group (32.90 ± 5.00)% was significantly lower than that in the negative control group (78.70 ± 4.70)% (t = 6.709, P < 0.01), indicating that AQP1 overexpression significantly inhibited the migration ability of WT cells. The results of Transwell assays showed that the number of perforated cells was significantly lower in the AQP1 overexpression group (179.00 ± 37.74) than in the negative control group (487.67 ± 41.72) (t = 5.486, P < 0.01), suggesting that AQP1 overexpression significantly inhibited the invasion capacity of WT cells. To investigate the potential impact of AQP1 overexpression on tumor metastasis, the influence of AQP1 on epithelial-mesenchymal transition (EMT) was assessed using Western blotting. The results demonstrated that the AQP1 overexpression group exhibited decreased levels of Vimentin and N-Cadherin, while the expression of E-Cadherin showed an opposite trend compared to the negative control group (Fig. 6A). Collectively, these findings suggest that AQP1 overexpression suppresses the migration, invasion, and EMT of WT cells.

The expression of F-actin in the negative control group increased significantly and metastasized to the pericellular periphery. Moreover, the cell volume increased, and lamellipodia, filopodia and stress fibers appeared. Together, F-actin expression was lower, and lamellipodia and filopodia formation decreased in the AQP1 overexpression group. By regulating the remodeling and dynamic changes of the cytoskeleton, cells were able to alter their shape, generate crawling pseudopods, and reorganize and rearrange microfilaments and microtubules, which could guide directional cell migration. The cell volume increased, and lamellipodia, filopodia, and stress fibers appeared, these changes collectively indicate an increase in cell migration capability. Overall, AQP1 overexpression decreased the formation of lamellipodia and filopodia, thereby inhibiting cell migration (Fig. 6C).