Integrated bioinformatics analysis of the NEDD4 family reveals a prognostic value of NEDD4L in clear-cell renal cell cancer

Data acquisition and processing

The kidney clear cell Carcinoma (KIRC) project of the TCGA (TCGA_KIRC) data set and GEO data sets (GSE40435 and GSE53757) were acquired and processed using the method described in our previous study (Zhang et al., 2020). Briefly, normalized RNA-seq data (607 cases, HTSeq-FPKM), phenotype information (985 cases), and survival data (979 cases) of the TCGA_KIRC were downloaded from the GDC hub of UCSC xena website (http://xena.ucsc.edu/public). Gene expression data of the GSE40435 and GSE53757 data sets (series matrix file) was downloaded from the Gene Expression Omnibus (GEO) database through the GEOquery package in the R software (version 3.6.2). After data processing, 526 tumor samples with survival data of the TCGA_KIRC data set, 202 samples of the GSE40435, and 144 samples of the GSE53757 were used for further analysis. All data for this study were acquired from public databases and did not require approval from the Ethics Committee.

Online database analysis

The transcriptional level of NEDD4L in renal cell carcinoma was verified in the Oncomine database (https://www.oncomine.org/resource/login.html). The expression data of 5 groups of renal cell carcinoma and normal tissues (transformed by log2) were retrieved and compared statistically. The default threshold was as follows: P <1E⁻⁴, multiple change >2, and the top 10% of genes. The staining intensity of NEDD4L in RCC samples and normal kidneys were downloaded from the Human Protein Atlas database (https://www.proteinatlas.org/).

Functional analysis and enrichment analysis

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were carried out using the clusterProfiler package in R version 3.6.2 (Hsing et al., 2020). Gene set enrichment analysis (GSEA) was used to investigate pathways enriched in the high- and low-risk subgroups. C2.cp.kegg.v7.1.symbols.gmt was chosen as the gene set database. These pathways were considered to be significantly enriched when the following criteria were met: nominal p-value < 0.05, false discovery rate q-value < 0.25, and absolute normalized enrichment score >1.

Development of the prognostic model

The least absolute shrinkage and selection operator (LASSO) Cox regression analysis was performed using the glmnet package in R. The analysis generated key gene signatures, and their corresponding coefficients were obtained by multi-variate cox analysis. A new score was calculated by multiplying the gene expression value of each gene and their corresponding coeffecient, as follows:

Score = − 0.06919*SOWAHB − 0.18041*WDR72 − 0.22124*EPB41L4A-DT − 0.07547*C4orf19 − 0.22609*CDKL2 − 0.05303*IRF6 − 0.16642*DHRS7 − 0.31446*NUPR2.

The risk score was calculated by subtracting the minimum score of the cohort from this score, and dividing the difference by the absolute value of the maximum score of the cohort, namely riskscore = (Score-min(Score))/abs(max(Score)).

Development and evaluation of the nomogram

The independent risk factors determined by multivariate Cox regression analysis were chosen to develop a nomogram for prediction of total overall survival (OS) probability. The consistency index (C-index) was used to evaluate the consistency between the occurrence frequency of the actual results and the prediction probability of the model. The nomogram was generated using the rms package in R software (Li et al., 2019a; Li et al., 2019b).

Statistical analysis

The data collected were analyzed by default as described using web resources. The rest of the data was analyzed using R software (version 3.6.2). The patients were divided into two subgroups by the median value of gene expression or risk score. Briefly, the correlation between the expression of NEDD4L and the clinicopathological features of ccRCC patients was examined by chi-square test with SPSS software (IBM Corp., Armonk, NY, USA). The expression of NEDD4L and other genes was analyzed by Spearman correlation analysis with R software (version 3.6.2). Students’ t-test was used to compare the level of gene expression. Univariate/multivariate Cox regression were performed by using the ‘survminer’ package in R. Kaplan–Meier survival analysis was used to compare the OS, relapse-free survival (RFS), disease specific survival (DSS) and progression free interval (PFI) of high- and low-CIFI groups Time-dependent receiver operator characteristic (ROC) analyses and subsequent calculation of the area under the curve (AUC) (Le, 2019; Le et al., 2019) were performed using the ‘timeROC’ package in R. Packages in R used for data analysis and graph plotting included ggstatsplot, ggplot2, ggpubr, limma, survminer, survival, tidyverse, dplyr, and plyr. A value of p < 0.05 was considered to be statistically significant (*, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001).

Patient characteristics

The data were collected from the TCGA_KIRC data set, including 72 normal subjects and 530 patients with primary tumors with clinical and gene expression data. The 526 patients with survival information were used for follow-up analysis (Table S1). There were 183 females (34.79%) and 343 males (65.21%). The onset age ranged from 26 to 88 years old, and the median age was 60 years old, 245 (46.58%) cases were younger than 60 years, and 281 (53.42%) cases were older than 60 years. The classification of the tumors according to the TNM (tumor, node, metastasis) staging system was as follows: 267 (50.76%) cases with T1 stage and 259 (49.24%) cases with T2-4 stage; 238 (45.25%) cases with N0 stage and 288 (54.75%) cases with N1/X stage; 418 (79.47%) cases with M0 stage, 106 (20.15%) cases with N1/X stage, and 2 (0.38%) cases of unknown stage; 239 (45.44%) cases with G1/2 stage, 279 (53.04%) cases with G3/4, and 8 (1.52%) cases of unknown stage. The classification of the tumors according to the Overall Stage Grouping (OSG) was as follows: 318 (60.46%) cases with stage I/II, 205 (38.97%) cases with stage III/IV, and 3 (0.57%) cases of unknown stage; 330 (60.46%) cases with stage I/II, 147 (38.97%) cases with stage III/IV, and 49 (0.57%) cases of unknown stage. The GSE40435 (101 normal subjects and 101 primary tumors), and GSE53757 (70 normal subjects and 70 primary tumors) data sets from the GEO database were used for verification.

NEDD4L is downregulated in ccRCC

We screened the transcription levels of NEDD4 family members in the TCGA_KIRC. The differentially expressed genes (DEGs) between normal and tumor tissues were analyzed using the following criteria: fold change (FC) >2 and FDR< 0.05. The results revealed that only NEDD4L was meaningfully down-regulated in tumors compared with that in normal tissues (Fig. 1A). Subsequently, two data sets from the GEO database were used for verification, and the results confirmed that NEDD4L was significantly down-regulated in ccRCC carcinoma tissues (Figs. 1B–1C). Data from the Oncomine database made up of five separate queues also confirmed that NEDD4L was meaningfully down-regulated in ccRCC (Fig. 1D). Immunohistochemical results showed that the expression level of NEDD4L in renal cell carcinoma was significantly decreased than that in normal kindey (Fig. S1).

Low expression of NEDD4L correlated to dismal prognosis in ccRCC

The TCGA_KIRC data set was used to investigate the prognostic value of NEDD4L in ccRCC. After data filtering, a total of 526 samples with gene expression and survival data were included. The median of the gene expression was set as the cut-off value. Kaplan–Meier survival analysis showed that patients with low NEDD4L expression had shorter OS (P < 0.0001), PFS (P < 0.0001) and DSS (P < 0.0001) (Figs. 1E–1G). Patients were divided into different subgroups according to sex, age, laterality, cancer status, histological grade, TNM stage, and OSG tumor stage. Analysis results of all but the N1 and M1 subgroups showed that patients with high NEDD4L expression had longer OS (Figs. 2 and 3, Figs. S2 & S3). Univariate and multivariate analysis were used to validate the prognostic value of NEDD4L in ccRCC. The univariate analysis results demonstated that low NEDD4L expression was associated with shorter OS (hazard ratio [HR]: 2.323; 95% CI [1.683−3.208]; P = 0.001) (Table 1). Multivariate COX regression analysis revealed that low expression of NEDD4L was still an independent factor associated with the deterioration of OS (HR: 1.905; CI [1.33–2.73], P < 0.001) (Table 1). These results suggested that NEDD4L may be an independent prognostic factor for ccRCC.

NEDD4L downregulation correlated with worse clinicopathological features of ccRCC

To investigate the relationship between the expression of NEDD4L and clinicopathological features in renal cell carcinoma, patients were divided into low and high expression subgroup on the base of the median expression of NEDD4L. Chi-square test showed that lower levels of NEDD4L were significantly associated with adverse prognostic features, including ccRCC with higher T stages (OR = 0.456 (0.322−0.646) for T1 vs. T2/3/4, P < 0.001), M stage (OR = 0.530 (0.342−0.821) for M0 vs. M1/X, P = 0.004), OSG tumor stage (OR = 0.429 (0.299−0.615) for stage I/II vs. stage III/IV, P < 0.001), histological grade (OR = 0.456 (0.321−0.647) for G1/2 vs. G3/4, P < 0.001) (Table 2).

Co-expression and enrichment analysis

Since NEDD4L was downregulated in ccRCC and correlated with worse prognosis of patients with ccRCC, we undertook to elucidate its potential role in this disease. In TCGA_KIRC and GSE40435 data sets, the correlation between the expression of NEDD4L and other genes in tumor samples was analyzed. A total of 494 genes in the TCGA_KIRC and 515 genes in the GSE40435 were identified as NEDD4L-related genes, using the following criteria: coefficient >0.5 and p value <  0.05. Among them, 206 NEDD4L-related genes were shared between the two datasets (Fig. S4, Table S2 ) and were subjected to GO and KEGG analysis. GO analysis revealed that the proteins encoded by these genes were mainly located in mitochondrial membrane and matrix, involved in the catabolism of many metabolites, such as lipid oxidation, fatty acid oxidation, fatty acid catabolism, monocarboxylic acid catabolism, and performed various molecular functions, such as DNA binding and ATP enzyme activity (Fig. 4A, Table S3). Accordingly, KEGG analysis showed that these genes play a major role in amino acid metabolism, oxidative phosphorylation and citrate cycle (Fig. 4B). These results suggest that NEDD4L is closely related to ccRCC cell metabolism.

Development of a NEDD4L-related prognostic model

To develop a NEDD4L-related prognostic model for ccRCC, we first performed univariate Cox regression analysis on the aforementioned 206 NEDD4L-related genes, which identified 183 genes with significant prognostic relevance (p < 0.05, Table S2). Patients from TCGA_KIRC (n = 526) were randomly and equally divided into two groups, one is to establish a LASSO Cox regression model, and the other is to evaluate the effectiveness of the model. The 183 NEDD4L-related genes with significant prognostic relevance were input into a LASSO Cox regression model, which generated eight key genes (Figs. 5A–5B), namely DHRS7, NUPR2, C4orf19, CDKL2, SOWAHB, WDR72, EPB41L4A-DT, and IRF6 (Fig. 5C). The risk score was calculated by inputting the selected signature genes into the above formula. The median risk score was set to a cutoff value, dividing patients into low-risk subgroups and high-risk subgroups. Prognostic analysis with Kaplan–Meier method showed that patients with ccRCC in the low risk subgroup had significantly longer DSS (Figs. 6A–6C). The analysis of the time-related receiver operating characteristic (ROC) in the TCGA training data set shows that the risk score had a favorable predictive value (Fig. 6D, area under the curve (AUC) at 1 year = 0.72, AUC at 3 years = 0.75, AUC at 5 years = 0.8). The predictions of the established risk scores were further confirmed by the TCGA test data set and the entire TCGA data set. The ROC analysis of the TCGA testing dataset showed that the AUC was 0.76 at 1 year, 0.69 at 3 years, and 0.68 at 5 years (Fig. 6E), while that of the TCGA entire dataset showed that the AUC was 0.74 at 1 year, 0.71 at 3 years, and 0.74 at 5 years (Fig. 6F). As reflected by the AUC values at 1 year, 3 year and 5 years, the signatures developed in our work had relatively better performance in predicting outcome of ccRCC patients (Fig. 6F, Fig. S5) (Zhang, Qiu & Yang, 2021; Fei, Chen & Xu, 2021; Chu et al., 2021; Xing et al., 2021).

Identification of NEDD4L-related signaling pathways by GSEA

Various biological processes are recognized to be involved in the occurrence and development of ccRCC. Accordingly, we hypothesized that the worse survival related to low expression of NEDD4L may be associated with some activated signal pathways in ccRCC. To evaluate this hypothesis, we performed a GSEA of the differences in the NEDD4L-related low risk group in TCGA_KIRC to verify the key signaling pathways related to NEDD4L expression. There were significant differences in the richness of MSigDB sets (c2.cp.v6.2 symbols) of multiple pathways (normalized P < 0.05, FDR < 0.25). According to the normalized enrichment score (NES), the signaling pathways most significantly associated with NEDD4L expression are shown in Fig. 7 and Table 3. In particular, NEDD4L expression was associated with biosynthesis of unsaturated fatty acids, insulin signaling pathway, and fatty acid metabolism, supporting a possible role of NEDD4L in the regulation of various metabolic processes

Establishment of a nomogram based on NEDD4L

Univariate and multivariate Cox regression analysis were achieved on the TCGA_KIRC ccRCC data set to investigate whether the risk score was an independent predictor of prognosis in patients with ccRCC. The adjustment of conventional clinical models, including age, sex, laterality, TNM stage, cancer status, histological grade, and tumor stage, indicated that the risk score was also an independent prognostic factor, which confirmed it had a strong ability to predict OS in patients with ccRCC (OR = 1.9329 (1.3429–2.782), P < 0.001, Figs. 8A–8B)

To establish a quantitative method for ccRCC prognosis, we used independent clinical risk factors to develop a nomogram (Fig. 8C). The point scale of nomogram is used to assign points to each variable based on multivariate Cox analysis. We drew a horizontal line to determine the points of each variable, then added the points of all variables together to calculate the total score of each patient, and then normalize it to a distribution of 0 to 45. The estimated 1 -, 3-and 5-year survival rates of ccRCC patients were calculated by drawing a vertical line between the total point axis and each prognostic axis (Li et al., 2019a; Li et al., 2019b). The C-index of the nomogram 0.7627102 (0.7255857–0.7998348) suggesting a good prediction effect. Thus, our results prove that nomogram was the best model for predicting the prognosis of renal cell carcinoma compared with a single risk factor.