Gene Expression Omnibus (GEO), NCBI’s public functional genomics data repository, was thoroughly queried for datasets related to ccRCC. The inclusion criteria are as follows: (1) Datasets contain human ccRCC tissue samples. (2) Raw sequencing data are available for each sample. (3) Datasets contain normal kidney tissues as control. Gene expression data in each dataset were downloaded, imported, and analyzed using R (version 4.0.3). The limma package was used to identify differentially expressed genes in each study (Ritchie et al., 2015). And genes with absolute log₂FC > 1 and adjusted p-value < 0.05 were defined as differentially expressed genes compared to the control group in a dataset. The robust rank aggregation method, which is parameter-free and robust to outliers, noise, and errors, was implemented to screen shared differentially expressed genes across all datasets (Kolde et al., 2012).

To determine the biological processes, molecular functions, cellular components, and pathways relevant to the differentially expressed genes, Gene Ontology (GO) term enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were implemented using clusterProfiler package, which supports both hypergeometric test and gene set enrichment analysis for many ontologies and pathways (Yu et al., 2012).

STRING¹ is a database of known and predicted PPI that incorporates both physical and functional associations between proteins. We used STRING to build a PPI network from common differentially expressed genes across all datasets and Cytoscape² software to visualize the network. The MCODE plugin (Bader and Hogue, 2003) in Cytoscape was used to find the functional modules in the PPI network.

The Kaplan–Meier plotter is an online tool capable of assessing the survival effect of 54k genes on a wide spectrum of cancer types, which integrates cancer sources from GEO, European Genome-phenome Archive, and TCGA Program (Nagy et al., 2018). We estimated the survival effect on ccRCC patients of all selected hub genes using the Kaplan–Meier plotter, and log-rank test was applied to compare the survival distributions of different groups. Hazard ratio (HR) with 95% confidence interval was calculated when there was a difference.

Gene Expression Profiling Interactive Analysis (GEPIA) is a Web server for cancer and normal gene expression profiling based on TCGA and Genotype-Tissue Expression (GTEx) data (Tang et al., 2017). We used GEPIA to profile target gene expression in different cancer types. The Human Protein Atlas (HPA³) database provides a comprehensive map of all the human proteins in cells, tissues, and organs. We explored the expression of target gene at the protein level in different human tissues with the help of HPA database.

The Open Targets Platform is a robust data integration for access to and visualization of potential drug targets associated with numerous diseases (Carvalho-Silva et al., 2019). Associations between target gene and diseases were evaluated using the Open Target Platform.

The level 3 data in The Cancer Genome Atlas Kidney Renal Clear Cell Carcinoma (TCGA-KIRC) dataset were acquired by UCSC Xena browser (Goldman et al., 2020). Only ccRCC patients with gene expression data were included in this study. ccRCC patients with a history of neoadjuvant treatment and without clear pathologic stage and histologic grade were excluded from this study. Extracted patient data include gender, pathologic stage, histologic grade, lymph node invasion, metastasis, overall survival, overall survival time, and genomic data of each TCGA-KIRC case. The genomic data include RNA-seq data of total gene expression, exon expression, methylation status, and DNA copy number alterations (CNAs) of a target gene. Cox regression analysis was used to evaluate the effect of different factors on ccRCC patients’ overall survival.

Formalin-fixed, paraffin-embedded ccRCC tissues were subjected to immunohistochemical analysis to estimate the expression of AOX1. Briefly, 10-μm cross section was deparaffinized and rehydrated in a series of xylene and graded ethanol solutions. Heat-induced antigen retrieval was performed in Tris-EDTA buffer (pH 9.0, G1203; Servicebio, China) in a 99°C water bath for 15 min. Intrinsic peroxidase activity was deactivated by incubation with hydrogen peroxide for 30 min at room temperature. After blocking with 5% goat serum in phosphate-buffered saline (PBS) for 1 h at room temperature, ccRCC tissue sections were incubated with AOX1 antibody (1:100 dilution in PBS, 19495-1-AP; Proteintech, China) overnight at 4°C. Isotype control of primary antibody was used as a negative control to differentiate non-specific background signal. Tissue sections were then incubated with a goat anti-rabbit horseradish peroxidase (HRP)-conjugated secondary antibody for 1 h at room temperature, followed by colorimetric detection using 3,3’-diaminobenzidine (DAB) (SA1022; Boster, China) per manufacturer’s protocol. Finally, the sections were counterstained with hematoxylin, and ImageJ (version 2.1.0/1.53c⁴) was used to measure and analyze the staining intensity of AOX1 following the published protocol (Crowe and Yue, 2019).

All ccRCC clinical samples were obtained with signed informed consent from the Department of Urology, Renmin Hospital of Wuhan University. This study was approved by the ethical review board of Renmin Hospital of Wuhan University.

Wherever applicable, a two-tailed Student’s t-test, Mann–Whitney U test, or Kruskal–Wallis test was used to test the significance of differences between groups. Dunnett’s correction is used to prevent family wise error rates when performing multiple hypothesis tests. Survival analysis was performed with Kaplan–Meier curves using log-rank test. Data were plotted and analyzed using GraphPad Prism 8 (GraphPad Software, CA) and presented as mean ± SEM unless stated elsewhere. We used the receiver operating characteristic (ROC) curve and calculated the area under the curve (AUC) score to show the diagnostic ability of AOX1 expression in classification of patients. Effect size is expressed as an unstandardized mean difference or the percentage of one group mean to another in this study.

The workflow of this study is shown in Figure 1A. A comprehensive search in the GEO database yielded eight eligible datasets, namely, GSE15641, GSE68417, GSE46699, GSE53757, GSE40435, GSE36895, GSE16441, and GSE16449. Detailed information about them is shown in Table 1. Eight GEO ccRCC datasets contain 331 normal kidney samples and 374 ccRCC samples in total. Differentially expressed genes in each dataset were first identified using the limma package in R with a cutoff value of absolute log₂FC > 1 and adjusted p-value < 0.05. The differentially expressed genes from each dataset range from 1,000 to 3,500. All differentially expressed genes from the eight datasets were then analyzed using the RobustRankAggreg package (Kolde et al., 2012) to determine the shared differentially expressed genes across all datasets. Among them, 211 and 145 genes were significantly downregulated and upregulated in ccRCC tissues compared with normal kidney tissues. In addition, the top 20 upregulated and downregulated genes can be found in Figure 1B. The volcano plots for the differentially expressed genes are shown in Supplementary Figure 1.

To gain mechanistic and functional insight into the differentially expressed genes between ccRCC tissues and normal kidney tissues, GO term and KEGG enrichment analyses were performed using the clusterProfiler package in R (Yu et al., 2012). The GO hierarchy contains three sub-ontologies: biological process, molecular function, and cellular component. The GO term enrichment analysis reveals that the major biological processes and molecular functions associated with upregulated differentially expressed genes are response to low oxygen levels, extracellular structure organization, receptor–ligand activity, cytokine activity, and monosaccharide binding (Figure 2A). The KEGG pathway analysis shows that the upregulated differentially expressed genes are primarily enriched in the hypoxia-inducible factor (HIF)-1 signaling pathway, focal adhesion, peroxisome proliferator-activated receptor (PPAR) signaling pathway, extracellular matrix–receptor interaction, and cholesterol metabolism (Figure 2B). For downregulated differentially expressed genes, the top biological processes and molecular functions include small-molecule process, organic acid catabolic process, carboxylic acid catabolic process, kidney development, and various transmembrane transporter activities (Figure 2C). Additionally, the downregulated differentially expressed genes are mainly enriched in PPAR signaling pathway, carbon metabolism, collecting duct acid secretion, fatty acid degradation, and diverse amino acid metabolism (Figure 2D). These results suggest that ccRCC cancer cells might better adapt to low oxygen conditions, which confer on them a survival advantage in hypoxic microenvironments in cancer tissue, and altered carbon, organic acid, and various amino acid metabolisms characterize a metabolic reprogramming in ccRCC cells.

We used the STRING database to establish a PPI network for the shared differentially expressed genes, and the network was processed by Cytoscape software for visualization (Figure 3A). In a PPI network, hub genes are those with more interactions with other nodes, and they are more likely to be essential genes for cells than other genes (Jeong et al., 2001). We identified 42 hub genes in the PPI network of all differentially expressed genes with a cutoff degree of 8, and the top 10 hub genes are KNG1, ITGB2, C3, TYROBP, ALB, FTCD, VWF, PROC, HRG, TIMP1 (Figure 3A). In addition, five functional modules were identified in the established PPI network by MCODE (Bader and Hogue, 2003), and the genes in each module are listed in Table 2. GO term and KEGG pathway analyses were applied to determine the function of each module (Figures 3B,C for module 4; Supplementary Figure 2 for other modules). The significant biological processes in GO terms and pathways enriched in KEGG analysis for module 4 are related to HIF-1 signaling, glycolysis, amino acid metabolism, and fatty acid metabolism (Figures 3B,C). Module 1 is mainly associated with protein modification and peptide degradation (Supplementary Figure 2A). Module 2 is enriched in cell adhesion, Pap1 pathway, PI3K/Akt pathway, Ras pathway, and mitogen-activated protein kinase (MAPK) pathway (Supplementary Figures 2B,C). Functions of module 3 involve cell chemotaxis, chemokine signaling, tumor necrosis factor (TNF) signaling, interleukin (IL)-17 signaling, Toll-like receptor signaling, and cAMP signaling (Supplementary Figures 2D,E). Module 5 is enriched in extracellular structure modification, glycan catabolism, protein digestion, and PI3K/Akt pathway (Supplementary Figures 2F,G).

Based on the centrality-lethality rule, which implies that the loss of a hub gene is more lethal than the loss of a non-hub protein in a PPI network (He and Zhang, 2006), the hub genes are of great biological significance. We checked the impact of hub genes on survival of ccRCC patients using the Kaplan–Meier plotter (Nagy et al., 2018). Strikingly, 22 out of 42 identified hub genes were associated with better or worse survival in ccRCC patients (Figure 3D), which validates their essential roles as hub genes.

PLG, vWF, and AOX1 are the three most protective hub genes for ccRCC patients, and ccRCC patients with a higher expression of these genes have 50% more overall survival advantage over those with a lower expression (Figure 3D). While the roles of PLG and vWF in ccRCC had been well elaborated, AOX1 has not yet been reported. We validated AOX1 mRNA expression with TCGA-KIRC dataset and AOX1 protein expression with immunohistochemistry. AOX1 gene expression data from 72 normal kidney tissues and 533 ccRCC tissues in TCGA-KIRC dataset were extracted and analyzed using UCSC Xena (Goldman et al., 2020). The expression of AOX1 mRNA in ccRCC tissues is only 34.2% (95% CI: 23.6–50.0%) as much as that in normal kidney tissues, which confirmed our differential expression analysis (Figure 4A). The ROC curve demonstrates the potential ability of AOX1 expression level as an indicator to discriminate between normal and ccRCC tissues (AUC = 0.69, p < 0.0001) (Figure 4B). In addition, 10 ccRCC clinical samples and six paired normal kidney samples were stained and quantified for AOX1 protein expression. ccRCC patients’ information can be found in Supplementary Table 1. Consistent with the results obtained from differential expression analysis and TCGA-KIRC data mining, AOX1 protein expression in normal kidney tissues is 2.5 times (95% CI: 1.6–3.6) higher than that in ccRCC tissues (Figures 4C,D). We further compared the localization and level of AOX1 protein expression in normal kidney tissues, and the result shows that proximal tubules have the most AOX1 protein expression, followed by distal tubules, while the glomeruli have the least AOX1 protein expression (Figures 4E,F).

To understand the effect of AOX1 on the prognosis of ccRCC patients, 516 eligible ccRCC patients from TCGA-KIRC dataset were subgrouped into two groups based on AOX1 mRNA expression in UCSC Xena (Goldman et al., 2020). Kaplan–Meier curve analysis was performed to compare the overall survival between the two groups. The result demonstrates that ccRCC patients with low AOX1 mRNA expression levels had worse overall survival than those with high AOX1 mRNA expression levels [high expression vs. low expression hazard ratio (HR): 0.52, 95% CI: 0.38–0.72] (Figure 5A). To clarify the definite relationship of AOX1 with different clinical factors of ccRCC, we estimated the mRNA expression levels of AOX1 in different patient subgroups. ccRCC patients aged between 40 and 60 years and older than 60 years have 2.1 times (95% CI: 1.2–3.9) and 2.3 times (95% CI: 1.3–4.2) higher AOX1 mRNA expression than that in patients younger than 40 years, respectively (Figure 5B). In addition, male ccRCC patients express 2.4 times more AOX1 mRNA (95% CI: 1.9–3.1) than female patients (Figure 5C). Higher histologic grade and lymph node invasion often correlate with poor prognosis and relapse. Moreover, ccRCC patients with lymph node invasion also have lower AOX1 expression levels (Figure 5D). No difference was found in AOX1 expression levels in ccRCC patients with respect to pathologic stages, metastasis, smoking, or histologic grade (Supplementary Figure 3). To investigate the effect of different factors on ccRCC patients’ overall survival, we performed Cox regression analyses. Cox regression confirmed that higher AOX1 mRNA expression (HR: 0.89, 95% CI: 0.83–0.96) has a protective role. In addition, higher pathologic stage (HR: 1.68, 95% CI: 1.43–1.96), higher histologic grade (HR: 1.53, 95% CI: 1.20–1.94), and older patients (HR: 1.03, 95% CI: 1.02–1.05) are associated with worse overall survival (Table 3).

To understand the potential mechanism of AOX1 gene dysregulation in ccRCC, we examined the DNA CNAs of the AOX1 gene in ccRCC. DNA CNAs are somatic variations that lead to duplication or deletion of certain sections of the chromosome and are quite common in human cancers (Tang et al., 2017). In TCGA-KIRC database, DNA CNA data were available for 503 ccRCC patients. Among them, 24 have either AOX1 homozygous deletion (–2) or heterozygous deletion (–1), and 78 of them have either AOX1 low-level (1) or high-level (2) amplification (Figure 6A). We checked the correlation between AOX1 mRNA expression and its DNA CNAs. Although AOX1 DNA amplification (1/2) has no impact on its mRNA expression, AOX1 DNA deletion (–1/–2) is associated with significantly lower AOX1 mRNA expression (Figure 6B), with only 22.2% (95% CI: 12.1–80.8%) AOX1 mRNA expression as much as that of AOX1 DNA normal (0) group. Kaplan–Meier survival curve shows that the ccRCC patient subgroup without AOX1 DNA deletion (CNA ≥ 0) has a survival advantage over that with AOX1 DNA deletion (CNA < 0) (CNA ≥ 0 vs. CNA < 0 HR: 0.25, 95% CI: 0.11–0.56) (Figure 6C). DNA methylation is another epigenetic tool to regulate gene expression in eukaryotic cells; therefore, we compared the methylation status of 20 CpG sites of the AOX1 gene in ccRCC patients subgrouped by AOX1 expression. Four out of the 20 examined CpG sites have significantly higher methylation levels in ccRCC patients with lower AOX1 mRNA expression (Figure 6D), and the methylation levels of these four CpG sites are negatively correlated with AOX1 mRNA expression (Figure 6E). All in all, these results imply that hypermethylation and copy number deletion may cause the downregulation of AOX1 in ccRCC.

To acquire a general knowledge on AOX1 expression in normal tissue, we queried the HPA database for the AOX1 gene. The HPA database offers thorough mRNA and protein expression profiles of AOX1 in 37 normal tissue types, and most of the examined tissues have low to no AOX1 mRNA or protein expression (Figure 7A). Among all the listed tissue types, liver, adrenal gland, testis, and kidney have the most abundant AOX1 protein expression, while liver, adrenal gland, kidney, and pancreas are the four tissues with the highest AOX1 mRNA expression level. Next, we compared AOX1 mRNA expression in 31 different types of tumor samples and normal tissue samples using GEPIA (Tang et al., 2017). AOX1 mRNA expression is significantly dysregulated in 23 out of 31 types of tumors, and one of the 23 types of tumors has upregulated AOX1 level, while 22 of the 23 types of tumors have downregulated AOX1 levels, including ccRCC (Figure 7B). We explored the potential diseases associated with AOX1 gene on the Open Target Platform (Carvalho-Silva et al., 2019) and identified 11 different kinds of diseases related to AOX1 gene dysregulation (Figure 7C). The disease spectrum covers almost all the major systems in the human body, including urinary system diseases, genetic familial or congenital diseases, immune system diseases, endocrine system diseases, respiratory or thoracic diseases, gastrointestinal diseases, pancreas diseases, cell proliferation disorders, nervous system diseases, and reproductive system diseases. The identified diseases are mostly cancers, and ccRCC is among them (Figure 7C). To conclude, AOX1 is abundant in normal kidney tissue, and it is downregulated not only in ccRCC but also in many other cancers.