Construction and validation of renal cell carcinoma tumor cell differentiation-related prognostic classification (RCC-TCDC): an integrated bioinformatic analysis and clinical study.

Текст сохранён. Часть иллюстраций пока доступна только в первоисточнике; загрузка будет проверена повторно.

Introduction

Globally, kidney cancer accounts for approximately 430,000 new cases and over 170,000 deaths each year, with its prevalence continuing to rise [ 1 ]. Renal cell carcinoma (RCC) is a highly heterogeneous malignancy that originates from the renal tubular epithelial cells and encompasses a broad spectrum of histopathological subtypes [ 2 ]. The primary subtypes include kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), and kidney chromophobe (KICH), representing approximately 75%, 15%, and 5% of all kidney cancer cases, respectively [ 3–5 ]. Each subtype is characterized by a unique genetic and molecular landscape within the tumor microenvironment (TME) [ 6–8 ]. In current clinical practice, localized RCC is effectively treated by nephron-sparing surgeries (generally preferred for smaller tumors ≤ 7 cm in diameter) and radical nephrectomy (indicated for larger tumors > 7 cm in diameter) for most patients [ 9 , 10 ]. Advanced RCC has also witnessed remarkable progress in targeted therapies and immunotherapy over recent decades [ 9 ]. However, no actual treatments are generally curative for advanced RCC, and many patients experience only modest improvements in survival and quality of life [ 11 ].

This awkward circumstance is largely attributed to the tumor heterogeneity of RCC, which highlights the significant molecular, cellular, and functional variations of tumors between patients, within a single patient, and even within the same tumor entity [ 12 , 13 ]. Notably, multi-directional partial differentiation of tumor cells can lead to distinct subclonal populations, each harboring unique genetic mutations, epigenetic modifications, and gene expression patterns [ 14 , 15 ]. This differentiation discrepancy among tumor cells and the resulting aggressive communications with other cellular components in the TME may drive heterogeneous malignant biological behaviors, adverse prognostic outcomes, and contribute to suboptimal therapeutic responses [ 16 , 17 ].

Therefore, unraveling the heterogeneity of RCC tumors and constructing a robust molecular classification system for RCC patients based on tumor cell differentiation are crucial for facilitating risk stratification and personalized therapies. In recent years, several high-impact studies with large RCC patient cohorts have proposed molecular classifications with significant therapeutic implications in the recent years [ 18–21 ]. However, these studies investigate molecular biomarkers based on bulk RNA sequencing (bulk RNA-seq) profiles, which capture average gene expression levels across entire cell populations, masking potential intra-tumoral heterogeneity in RCC. In contrast, single-cell RNA sequencing (scRNA-seq) enables gene expression analysis at the single-cell level, providing deeper insights into the RCC TME and dynamic processes such as cellular differentiation trajectories.

Accordingly, the present study aimed to explore the differentiation trajectory of heterogeneous RCC tumor cells and establish a novel molecular classification system for RCC patients by integrating scRNA-seq and bulk RNA-seq analyses with validation using large-sample pathological specimens and a retrospective clinical cohort.

Methods

Ethics approval and data extraction

The Ethics Committee of Xinhua Hospital Affiliated to Shanghai Jiao Tong University School of Medicine approved this study (XHEC-C-2021-145-1). A total of 14 scRNA-seq samples, including six benign adjacent kidney (BA) samples and eight RCC samples, were downloaded from the Gene Expression Omnibus (GEO) database (accession no. GSE159115 ) [ 17 , 22 ]. The clinical information, genome profiles, bulk RNA-seq profiles, reverse phase protein array (RPPA) profiles, and assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) profiles of 537 KIRC samples were obtained from The Cancer Genome Atlas (TCGA) ( https://tcga-data.nci.nih.gov/tcga/ ) database. External validation data was retrieved from GSE22541 , TCGA-KIRP, and TCGA-KICH cohorts [ 23 ]. Additional gene expression profiles of 12 metastatic RCC samples were obtained from the MET500 ( https://met500.med.umich.edu/ ) data collection [ 24 ]. The transcription factor (TF) profiles were obtained from the Cistrome cancer database ( http://cistrome.org/ ) [ 25 ]. More importantly, 370 patients diagnosed with kidney neoplasms, who underwent curative surgery at Xinhua hospital from 2016 to 2018 were included to construct a substantial retrospective cohort. Written informed consent was obtained from all participants upon their admission to the hospital. For participants who were minors, written informed consent was obtained from their legal guardians. Additionally, assent was obtained from minors where appropriate. Comprehensive statistical data, including demographic characteristics, clinical and pathological features, and follow-up prognostic information, were all collected.

Integrated bioinformatic analysis

ScRNA-seq cluster analysis

The scRNA-seq profiles of six BA samples and eight RCC samples from GEO were obtained for comprehensive data process. During quality control, empty droplets, droplets containing damaged or dying cells, and doublets or multiplets were excluded. Feature selection was carried out by identifying the 2000 most highly variable genes, which were then subjected to principal component analysis (PCA) for linear dimensionality reduction and ranking of principal components (PCs). Uniform manifold approximation and projection (UMAP) was subsequently applied for non-linear dimensionality reduction and visualization using the identified 20 final PCs by the inflection point. Differentially expressed genes (DEGs) were identified with a false discovery rate (FDR) < 0.05 and |average log 2 (fold change) value| > 0.5. Cell annotation was then conducted based on cell markers in DEGs of each unsupervised cluster to identify specific cell types. In addition, cell cycle and cell communication analyses were performed. After completing first-step scRNA-seq cluster analysis, epithelial cell profiles were subsequently retrieved for subsequent second-step subclustering. Finally, marker genes of each cell type were visualized by feature plots, violin plots, and heatmaps, generated by the Seurat R package and SCANPY in the Python 3.6 environment [ 26 , 27 ].

Pseudotime analysis of epithelial cells

Pseudotime analysis of epithelial cells was then performed using the Monocle2 R package. Cell trajectory analysis, differentiation fate-related gene analysis, branched expression analysis modeling (BEAM), and clustering were sequentially performed to identify differentiation fate-related genes and potential epithelial cell differentiation trajectories. Specifically, low-quality cells and genes with low dispersion were filtered out, and the top 1000 genes for ordering were screened out. Subsequently, the ‘DDRtree’ method was employed to reduce dimensionality and identify branching points across differential epithelial cell fates. Finally, differentiation fate-related genes identified by DEG analysis ( p < 0.05) in epithelial cells were termed ‘monocle2-sig genes’, while genes selected for ordering in the differentiation trajectories were designated as ‘monocle2-order genes’.

Determination of the tumor cell differentiation-related prognostic genes (TCDGs)

Furthermore, bulk RNA-seq profiles and clinical information from the TCGA-KIRC cohort were extracted for prognostic analyses. Thereafter, we obtained statistically significant genes in Kaplan-Meier survival analysis ( p < 0.00001) as ‘KM-sig genes’ and in univariable Cox regression ( p < 0.00001) as ‘uniCox-sig genes’. A Venn diagram was then generated to visualize the intersection of the four defined gene sets as tumor cell differentiation-related prognostic genes (TCDGs). A co-expression network of the TCDGs was also constructed to explore the correlation patterns among them.

Consensus clustering analysis of TCGA-KIRC patients based on TCDGs

Using the TCDGs, the ‘ConsensusClusterPlus’ R package was implemented to categorize bulk RNA-seq profiles into 2 to 9 clusters over 50 iterations. The optimal number of clusters (k value) was determined based on the relative change in the area under the cumulative distribution function (CDF) curve. The specific criteria included a sharp bend in the CDF curve, absence of clusters with very few samples, low intra-group correlation, and high inter-group correlation. The final classification of RCC patients into k distinct subtypes was defined as RCC tumor cell differentiation-related prognostic classification (RCC-TCDC). K-M survival analysis was applied to reveal prognostic differences among these RCC-TCDC subtypes. What’s more, PCA was employed on the expression profiles of the TCDGs to provide extra validation with continuous variables. Using the principal components obtained from PCA, we computed the Euclidean distance (PCA score) for each KIRC patient as the PCA-KIRC-score, and the differences among the identified subtypes were then explored. Moreover, based on the optimal cutoff value of PCA-KIRC-score determined by the ‘surv_cutpoint’ function, patients were divided into high-PCA-KIRC-score and low-PCA-KIRC-score subgroups. Then an additional K-M survival analysis was utilized to examine prognostic differences between them. Moreover, external validation of the TCDGs was conducted in GSE22541 KIRC cohort, TCGA-KIRP, and TCGA-KICH cohorts. Finally, profiles of 15 kidney cell lines from the Cancer Cell Line Encyclopedia (CCLE) were extracted and then classified according to the previously determined RCC-TCDC subtypes. Profiles from the Genomics of Drug Sensitivity in Cancer (GDSC) were then retrieved to find potential therapeutic drugs sensitive to each subtype.

Multi-omics analysis of TCGA-KIRC patients

A multi-omics analysis of TCGA-KIRC patients was then implemented. Initially, genomic mutations with different patterns in high-PCA-KIRC-score and low-PCA-KIRC-score subgroups were identified. Moreover, the ­differences in tumor mutation burden (TMB) and ­microsatellite instability (MSI) between high-PCA-KIRC-score and low-PCA-KIRC-score subgroups were unraveled. In addition, correlation analysis was utilized to investigate the correlation between PCA score and TMB or MSI. Subsequently, K-M survival analysis was performed to explore the prognostic relevance of TMB and PCA scores. Additionally, the copy number variation (CNV) and chromosomal locations of TCDGs were analyzed. Subsequently, the CIBERSORT algorithm was adopted to reveal immune cell infiltration in the RCC-TCDC subtypes [ 28 ]. In addition, the differential expression of programmed death-ligand 1 (PD-L1) and the differential response to immune checkpoint blockade measured by immunophenotype score (IPS) were compared between high-PCA-KIRC-score and low-PCA-KIRC-score subgroups were examined [ 29 ]. Finally, correlation analysis was performed to conduct co-expression analyses between PCA-KIRC-score and immune cell infiltration.

Differential expression analysis between primary and metastatic tumors

Differential expression analysis was utilized to uncover significant differences of DEGs and TFs between primary and metastatic RCC tumors from TCGA and MET500, identifying metastasis-related TCDGs (mTCDGs) and metastasis-related TFs (mTFs). Additionally, gene set variation analysis (GSVA) was implemented to delineate the variation of differentially enriched cancer hallmarks. Over-representation analysis (ORA) was then utilized to explore gene sets from the Molecular Signatures Database (MSigDB) ( https://github.com/tomastokar/gsoap ) in which mTCDGs were enriched.

Construction of the multivariable prognostic model

To reduce feature selection bias and avoid overfitting of the multivariable prognostic model, the least absolute shrinkage and selection operator (LASSO) regression was employed to screen out a mTCDGs-based gene signature. In this regression model, the response variable was the overall survival of TCGA-KIRC patients, while the independent variable was the normalized expression matrix of candidate mTCDGs. The model’s penalty parameter (λ) was computed by tenfold cross-validation, with the value of λ corresponding to the lowest partial likelihood deviance. The risk score was calculated by the following formula: risk score x = ∑ i = 1 n β i × mTCDG i

In the formula, ‘x’ represented the number of patients, ‘n’, ‘ β i ’ and ‘mTCDG i ’ represented the number, corresponding regression coefficient, and normalized expression level of each mTCDG, respectively.

Validation of the multivariable prognostic model

Upon multivariable prognostic model construction, patients were categorized into high-risk and low-risk subgroups based on the median risk scores. Cross-validation with a 60: 40 train/test split was then conducted for validation. Differential expression and prognostic differences were visualized between the high-risk and low-risk subgroups. Additionally, receiver operating characteristic (ROC) curves and the area under the curve (AUC) were applied to assess the predictive accuracy of the model. Furthermore, independent prognostic analysis was performed using univariable and multivariable Cox regression analyses for risk score along with age, T stage, N stage, and M stage. Chi-square tests were then utilized to examine differences in OS sensor, age, gender, race, grade, M stage, N stage, T stage, and cancer status between the high-risk and low-risk subgroups. Moreover, gene set enrichment analysis (GSEA) was performed to uncover the top five enriched Gene Ontology (GO) biological process (BP), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and cancer hallmarks of high-risk and low-risk subgroups.

Immunological relevance of the multivariable prognostic model

The CIBERSORT algorithm was implemented to visualize the relative proportion of 22 infiltrating immune cells. Additionally, single-sample GSEA (ssGSEA) was adopted to unveil the enrichment scores of 29 immune functions. Following that, significant differences in the fraction of 22 infiltrating immune cells and the enrichment of 29 immune functions were explored between the high-risk and low-risk subgroups. In addition, K-M survival analysis was applied to explore the prognostic differences between high and low immune cell infiltration and immune function enrichment. Moreover, We further investigated the differences of the established six distinct immune subtypes, including wound healing (C1), IFN-γ dominant (C2), inflammatory (C3), lymphocyte depleted (C4), immunologically quiet (C5), and TGF- β dominant (C6) between high-risk and low-risk subgroups by Chi-square tests [ 30 ]. In addition, in order to explore the immunotherapy response characteristics, the Tumor Immune Dysfunction and Exclusion (TIDE) database ( http://tide.dfci.harvard.edu ) was employed to reveal the differences in TIDE, T cell dysfunction, T cell exclusion, CD274, MSI, tumor-associated macrophage M2, myeloid-derived suppressor cell (MDSC), cancer-associated fibroblast (CAF), CD8, IFN-γ, and Merck18 (tumor inflamed score) scores between high-risk and low-risk subgroups [ 31 , 32 ].

Regulatory network construction, ATAC-seq validation, and inhibitor prediction for each RCC-TCDC subtype

Correlation analysis was implemented to explore the possible correlation relationships between key mTCDGs and mTFs, metastasis-related cancer hallmarks, infiltrating immune cells, immune functions, and RPPA proteomics in three RCC-TCDC subtypes. The most significant pairs (with top high correlation coefficients and p < 0.05) were selected for the final construction and visualization of a regulatory network for each subtype. Moreover, co-expression heatmaps were generated for these pairs. Additionally, GO and KEGG functional enrichment analyses were performed in each subtype. To further validate the transcriptional mechanism of key mTCDGs, ATAC-seq profiles were utilized to detect chromatin accessibility at their genomic loci. Finally, the Connectivity Map (CMap) ( https://portals.broadinstitute.org/cmap/ ) was utilized to identify potential inhibitors for each subtype. The enrichment score (ranging from −1 to 1) of each small molecular compound was measured, and the top 10 candidate inhibitors with high connectivity and p < 0.05 were filtrated for each subtype.

Large-sample pathological specimen and retrospective clinical study validation

Immunohistochemical (IHC) staining and scoring

Based on the differential expression of mTCDGs in the three subtypes, we selected specific markers for each subtype (S100A11, SH3BGRL3, TSPAN7, AQP1, and NPR3) with significantly elevated expression levels, taking the regulation networks and substantial relevant literature consultation into account. Afterwards, the tumor tissues of the 370 patients with kidney neoplasms were gathered to undergo a series of rigorous IHC staining processes on tissue microarrays (TMAs). Initially, the tissues were fixed, embedded, and sliced into paraffin samples. Subsequently, dewaxing and rehydration procedures were carried out, followed by antigen retrieval to expose specific antigens (markers of each subtype) via EDTA (pH = 9.0) or citric acid (pH = 6.0) antigen-retrieving buffers. Followingly, in order to block endogenous peroxidase, the samples were incubated in 3% hydrogen peroxide solution and subsequently washed in phosphate buffer saline solution. Further, the prevention of non-specific antibody binding was achieved by treating the samples with 3% bovine serum albumin. Primary antibodies, encompassing S100A11 (Proteintech, 10237-1-AP, 1:50), SH3BGRL3 (Affinity, DF9905, 1:500), TSPAN7 (Proteintech, 18695-1-AP, 1:50), AQP1 (Proteintech, 20333-1-AP, 1:500), and NPR3 (Proteintech, 26706-1-AP, 1:100), were then introduced to interact with the antigens. Subsequently, the secondary antibodies marked with the enzyme HRP were applied to bind to the primary antibodies upon the former interaction. For color development, the specific substrate 3,3′ Diaminobenzidine was added to react with horseradish peroxidase, which in the presence of the target antigens, induced observable color changes. To enhance visualization, a counterstain using haematoxylin was employed, providing contrast and highlighting the nuclei. Following this, the samples underwent dehydration through a series of graded alcohols. All TMAs were eventually imaged utilizing a microarray scanner. Afterwards, two experienced pathologists were enlisted to independently read and score the tissue spots on these TMAs. The conflicted results would be carefully re-assessed by a third pathologist, ensuring a thorough and reliable evaluation of the obtained data. Specifically, the IHC score of each sample was determined by both staining intensity (‘I’) and density (‘D’). ‘I’ indicated the intensity of stained markers in renal tumor cells, which was categorized into negative (0), weak (1), moderate (2), and intense (3). Moreover, ‘D’ defined the percentage of positive-staining cells, which was classified into 0% (0), 1–10% (1), 11–40% (2), 41–70% (3), and 71–100% (4). Eventually, the IHC score was calculated by I*D (ranging from 0–12), which was further divided into not detected (0), low (1–4), medium (5–8), and high (9–12). Finally, we classified all patients with thorough IHC staining data into three RCC-TCDC subtypes according to the overall differential expression of the five markers.

Comprehensive statistical analysis for Xinhua cohort

In our retrospective cohort consisting of 370 kidney neoplasm patients who received curative surgery at Xinhua hospital between 2016 and 2018, thorough clinicopathological information was gathered. Moreover, we performed systematic follow-up in March, 2021 to obtain additional prognostic data. After meticulously cleaning and processing these data, we first conducted a detailed descriptive demonstration of all gathered variables, including age, gender, histopathological classification, Fuhrman nuclear stage, TNM stage, stage, overall survival (OS), progression-free survival (PFS), number of patients in each RCC-TCDC subtype, and IHC staining score of each marker gene in Xinhua cohort. Following this, the differential IHC scores of the five markers across three subtypes were illustrated by violin plots. Additionally, we utilized the ‘surv_cutpoint’ function to reveal the optimal cut-off value of the IHC staining scores of the five markers, based on which we divided the patients into high-expression and low-expression subgroups. Following that, K-M survival analysis was implemented to unveil the differences in OS and PFS between the two subgroups. Subsequently, we investigate the clinical relevance of RCC-TCDC by employing Chi-square tests between various clinicopathological variables and three subtypes. Moreover, we conducted K-M survival analysis to illustrate the differences in OS and PFS among three subtypes. Following that, multivariable Cox regression analysis was carried out to explore the independent prognostic value of RCC-TCDC on OS and PFS of kidney neoplasm patients. Nomograms were then generated accordingly. Eventually, comprehensive model diagnostics were performed by residuals analysis, Schoenfeld tests, calibration curves, ROC curves, and the calculation of risk scores.

Statistical analysis

All statistical analyses were conducted by R version 4.2.2 (Institute for Statistics and Mathematics, Vienna, Austria; www.r-project.org ), Python version 3.6 software ( https://www.python.org/ ), and Strawberry Perl version 5.30.0.1 software ( https://www.perl.org/ ). For descriptive statistics, the mean ± standard deviation (SD) was employed for continuous variables with a normal distribution, while the median (range) was adopted for variables with a non-normal distribution. Moreover, parametric tests were used for group means when variables had a normal distribution and homogeneity of variance, including the t-test, ANOVA, and the Pearson correlation. Non-parametric tests were used for group medians when variables didn’t conform a normal distribution or homogeneity of variance, such as the Wilcoxon test, the Mann-Whitney U test, the Kruskal-Wallis test, and the Spearman correlation. If not specified above, those with a two-tailed p < 0.05 and FDR < 0.05 were considered statistically significant.

Results

Overview of the study

The schematic representation of this integrated study was outlined in Figure 1A . Additionally, the entire analysis process was detailed in a flowchart ( Figure S1 ). Specifically, based on an integrated bioinformatics analysis, we identified 32 TCDGs and developed a novel molecular classification system, RCC-TCDC, classifying TCGA-KIRC patients into three distinct molecular subtypes. The prognostic and multi-omics relevance was then validated. Moreover, we utilized 10 mTCDGs to establish a multivariable prognostic model and ­constructed regulatory networks with inhibitor prediction for each RCC-TCDC subtype. Notably, we then employed a large-sample pathological specimen and a retrospective clinical study for robust validation of RCC-TCDC by clinicopathological and prognostic implications. Eventually, a diagnostic kit will be designed and applied for a patent in the near future.

Identification of epithelial cells by first-step scRNA-seq cluster analysis of all cells

In the first-step scRNA-seq cluster analysis of six BA and eight RCC samples, we obtained 23 Seurat clusters (0–22), six major cell types, including B cell, endothelial cell, epithelial cell, myeloid cell, NK/T cell, and vascular smooth muscle cell (vSMC), and 14 sample origins (BA1-6 and RCC1-8), visualizing in the UMAP dimension reduction plots ( Figure 1B ). Cell annotation of the six cell types was displayed by feature expression of various marker genes ( Figure S2A ). Additionally, the cell number and cell proportion of the six cell types in each sample were demonstrated ( Figure S2B ). In addition, the average expression levels of four classic marker genes, LYZ, CD3D, CD24, and ACTA2, in the six cell types were illustrated, validating that LYZ, CD3D, CD24, ACTA2 expression were significantly up-regulated in myeloid cells, NK/T cells, epithelial cells, and vSMC, respectively ( Figure S2C ). Moreover, significantly up- and down-regulated DEGs of the 23 Seurat clusters were displayed in Figure S2D . Furthermore, differential expression of the top five DEGs in each cell type was revealed ( Figure S3A ). What’s more, significant cell communication among the six cell types was exhibited in Figure S3B, C . Specifically, myeloid cell, epithelial cell, and endothelial cell had predominantly active cell communications with other cells. Eventually, the distribution of G1, G2M, and S stages in a cell cycle was depicted, and cell cycle scores of G2M and S in each cell type were calculated ( Figure S3D, E ).

Delineating the differentiation trajectories of epithelial cells by second-step scRNA-seq cluster analysis and monocle2 pseudotime analysis

Histopathologically, RCC originated from the renal tubular epithelial cells, whose malignant differentiation process and heterogeneity played a pivotal role in the initiation, development and progression of RCC. Subsequently, we extracted epithelial cells to proceed with the second-step scRNA-seq cluster analysis and obtained 11 distinct epithelial clusters, as illustrated in Figure 1C . Moreover, the average expression levels of four malignant markers, MKI67, CD44, CD24, and PROM1, along with differential expression of the top five DEGs in the 11 epithelial clusters were also exhibited ( Figure S4A, B ). Besides, the significant cell communications among the 11 clusters were unveiled, in which epithelial clusters 1, 4, 5, 6, and 9 were shown to take over the coring communication status ( Figure S4C, D ). Overall, the obtained results revealed substantial heterogeneity of epithelial cells in RCC tumor.

Based on the identified 11 epithelial clusters, Monocle2 pseudotime analysis was implemented to reveal putative tumor cell differentiation trajectories ( Figure 1D ). Specifically, five states of epithelial clusters were found to differentiate along the pseudotime trajectories with two bifurcation points, each branching into two different sub-trajectories ( Figure S5A–C ). It was shown that epithelial clusters 1 and 4 were mainly distributed in initiating differential state 1, while clusters 0, 3, 5, and 8 had dominant distribution in ending differential states 3, 4, and 5, meanwhile cluster 9 was scattered in all five states. What’s more, BEAM analysis demonstrated the potential association between the differential expression of differentiation-fate-related genes of the five states and tumor cell differentiation in Figure S5D . Ultimately, we identified 1576 differentiation-fate-related DEGs ( p < 0.05) as ‘monocle2-sig genes’, and 441 genes considered for ordering in the differentiation trajectories as ‘monocle2-order genes’. So far, we have delineated the differentiation trajectory of RCC tumor cells and revealed candidate tumor cell differentiation-fate-related biomarkers. However, the prognostic relevance of these biomarkers required further exploration.

Determination of the TCDGs by integrating analysis of bulk RNA-seq profiles and prognostic data of patients from TCGA-KIRC cohort

Furthermore, 537 patients from TCGA-KIRC cohort were included for prognostic analysis, the baseline information was displayed in Table S1 . We identified 84 and 32 significant prognostic genes through K-M survival analysis ( p < 0.00001) as ‘KM-sig genes’ and through univariable Cox regression ( p < 0.00001) as ‘uniCox-sig genes’, respectively. Subsequently, we exerted a Venn diagram to show the overlap of 32 key genes among the four gene sets as the TCDGs ( Figure 1E ). In addition, the K-M survival curves of them were visualized in Figure S6A . Moreover, the results of prognostic analyses were detailed in Table 1 . Specifically, ACAT1, ADGRF5, AQP1, BBOX1, CYS1, DDC, GLYAT, HAO2, NPR3, PDK4, PLPP1, RCAN2, RHOB, SGMS2, SULT1C4, TJP2, and TSPAN7 were protective factors for KIRC patients. In contrast, C4orf48, CCNB1, CDK1, CXCL1, CXCL2, MDK, MTHFD2, NFKBIZ, PSAT1, PTTG1, S100A11, SAA1, SH3BGRL3, SLC25A37, and TIMP1 served as risk factors. Besides, we constructed a co-expression network demonstrating significant correlation relationships ( p < 0.0001) among the 32 TCDGs of five monocle2 states ( Figure S6B ).

Construction of RCC-TCDC and identification of three distinct subtypes by consensus clustering analysis

Through consensus clustering analysis, the dynamic change of the CDF curve with increasing consensus index and the AUC with consensus cluster number k changing from 2 to 9 were revealed in Figure S7A, B . Importantly, the k value was determined as 3 ( Figure S7C ), and hence we eventually constructed the RCC-TCDC by classifying the TCGA-KIRC patients into three molecular subtypes based on the 32 TCDGs. Specifically, the differential expression patterns of the 32 TCDGs across five states within three molecular RCC-TCDC subtypes were visualized ( Figures 2A and S7D). Notably, we found that S100A11, CXCL2, NFKBIZ, SH3BGRL3, CXCL1, SLC25A37, TIMP1, SAA1, PSAT1, C4orf48, MTHFD2, PTTG1, MDK, CDK1, and CCNB1 displayed relatively high expression in subtype 1. These up-regulated TCDGs were predominant in monocle2 states 1 and 5, representing the initial and early differentiation states in the pseudotime trajectory. Whereas subtype 3 exhibited elevated expression levels of RHOB, PDK4, SULT1C4, PLPP1, AQP1, HAO2, DDC, BBOX1, GLYAT, NPR3, SGMS2, TJP2, CYS1, and TSPAN7. Most of them were predominant TCDGs for monocle2 states 2, 3, and 4, indicating the late and terminal differentiation states in the pseudotime trajectory. Moreover, subtype 2 showed over-expression of predominant TCDGs for both monocle2 states 4 and 5, including SGMS2, TJP2, ADGRF5, TSPAN7, and RCAN2. Taken together, we defined the three molecular RCC-TCDC subtypes 1, 2, and 3 as RCC-low differentiation (RCC-LD), RCC-moderate differentiation (RCC-MD), and RCC-high differentiation (RCC-HD), respectively.

Subsequent K-M survival analysis revealed significant prognostic differences ( p < 0.001) in OS among the three subtypes ( Figure 2B ), with RCC-LD showing the worst prognosis, RCC-HD exhibiting the most favorable prognosis, and RCC-MD falling in between. To further explore the prognostic relevance of RCC-TCDC and reduce bias caused by categorical variables, we conducted PCA and calculated PCA-KIRC-score. The discriminative distribution of three RCC-TCDC subtypes was demonstrated in the PCA dimension reduction projection map ( Figure S7E ). Additionally, a significantly high PCA-KIRC-score was observed in RCC-LD compared to both RCC-MD ( p < 0.001) and RCC-HD ( p < 0.001) in Figure 2C . Subsequently, upon stratifying patients into two subgroups by the median value of PCA-KIRC-score, an additional K-M survival curve unveiled a significantly ( p < 0.001) worse prognosis in high-PCA-KIRC-score subgroup ( Figure 2D ). This further reinforced the unfavorable prognostic status of RCC-LD patients. In addition, these close links were depicted in a Sankey diagram ( Figure S7F ). Consequently, a novel molecular prognostic classification system, RCC-TCDC was preliminarily validated.

Moreover, external validation by the GSE22541 KIRC cohort was conducted, results in Figure S8 showed the same prognostic significance of TCDGs as in the TCGA-KIRC cohort. To further explore the implications of TCDGs in other pathological subtypes of RCC, we utilized TCGA-KIRP and TCGA-KICH cohorts for external validation. Specifically, we successfully established another classification system on KIRP and KICH patients (clusters 1, 2, and 3) based on the 32 TCDGs ( Figure S9A–D ). Significantly differential expression of the 32 TCDGs and significantly prognostic differences ( p < 0.001) were also observed within the three clusters ( Figures S9E and S10A ). Scattered distribution of three clusters in PCA was shown in Figure S10B . Additionally, K-M survival analysis of TCDGs in TCGA-KIRP ( Figure S10C ) and TCGA-KICH ( Figure S11 ) supported the broad prognostic relevance of TCDGs in other RCC subtypes.

Last but not least, the CCLE database was utilized to delineate the characteristic expression of the 32 TCDGs across 15 distinct kidney cell lines ( Figure 2E ). Subsequently, leveraging profiles from the GDSC database, we identified noteworthy differences in the AUCs of eight chemotherapeutic drugs among three RCC-TCDC subtypes through Kruskal-Wallis tests. The results indicated that RCC-LD might be potentially sensitive to selumetinib ( p = 0.0084), SN-38 ( p = 0.012), ZM-44739 ( p = 0.016), temsirolimus ( p = 0.026), SNX-2112 ( p = 0.039), veliparib ( p = 0.042), and cisplatin ( p = 0.048), while RCC-MD could exhibit potential sensitivity to bleomycin ( p = 0.011) ( Figure 2F ).

Multi-omics landscapes of RCC-TCDC

In genomic analysis, we presented the top 20 mutated genes across eight mutation types in both high-PCA-KIRC-score and low-PCA-KIRC-score subgroups ( Figure S12A, B ). Higher mutation frequency was shown in low-PCA-KIRC-score subgroup (81.15%) than high-PCA-KIRC-score (77.94%), both with VHL mutation dominating. Moreover, insignificant differences of TMB and ­significant differences ( p = 0.016) of MSI between low- and high-PCA-KIRC-score subgroups were displayed ( Figure S12C ). Besides, significantly positive ( R = 0.13, p = 0.022) and negative (R = −0.15, p = 0.0057) correlation relationships between TMB, MSI and PCA-KIRC-score across the three subtypes were revealed, respectively ( Figure S12D ). Subsequently, after dividing patients into two subgroups by the median value of TMB, K-M survival analysis revealed a better prognosis ( p < 0.001) in high-TMB (H-TMB) subgroup than low-TMB (L-TMB) subgroup ( Figure S12E ). Additional K-M survival curve unveiled significant prognostic differences ( p < 0.001) among the four subgroups ( Figure S12F ). Specifically, the L-TMB+high-PCA-KIRC-score subgroup exhibited the worst prognosis, while H-TMB+low-PCA-KIRC-score subgroup exhibited the best prognosis. In addition, CNV frequencies of the 32 TCDGs were illustrated in Figure S12G , which indicated the notable net gains of chromosomal regions of PTTG1 and losses of chromosomal regions of SH3BGRL3. Finally, the chromosomal position of the 32 TCDGs, along with their CNV status, were visualized in a color-coded circular plot ( Figure S12H ).

In immunomic analysis, the differential immune infiltration of 23 immune cells across three RCC-TCDC subtypes was illustrated ( Figure S13A ). Specifically, RCC-LD was indicated to have a generally high immune infiltration of various immune cells, including activated B cell, activated CD4 + T cell, activated CD8 + T cell, activated dendritic cell, MDSC, macrophage, natural killer T cell, natural killer cell, neutrophil, regulatory T cell (Treg), and T follicular helper cell, while RCC-MD generally exhibited a low immune infiltration, with RCC-HD falling in the middle. What’s more, positive correlation relationships between PCA-KIRC-score and multiple immune cells, including activated CD4 + T cell, activated CD8 + T cell, activated dendritic cell, MDSC, macrophage, natural killer T cell, natural killer cell, Treg, T follicular helper cell, type 1 helper cell, and type 2 helper cell were revealed by co-expression analysis ( Figure S13B ). Additionally, low-PCA-KIRC-score subgroup showed a notable elevation ( p = 0.00042) of PD-L1 expression compared to high-PCA-KIRC-score subgroup ( Figure S13C ). However, differential response to immune checkpoint blockade (anti-CTLA4 and anti-PD1 treatment) measured by IPS between high-PCA-KIRC-score and low-PCA-KIRC-score subgroups showed no significance ( Figure S13D ). Taken together, we revealed remarkable multi-omics relevance of RCC-TCDC. Particularly, this seemingly contradictory immune infiltration features conferred consistence with established immunological studies and notable clinical implications for RCC.

Differential expression analysis between primary and metastatic RCC tumors

Furthermore, we were committed to delineate the whole malignant biological processes underlying tumor cell differentiation, from mutation, initiation, development, metastasis, to colonization and development in a new site, forming a vicious progression cycle. Hence, we proceeded with differential expression analysis between primary and metastatic RCC tumors. All DEGs were presented in Figure S14A, B . Subsequently, 32 TCDGs were all identified to be DEGs between primary and metastatic tumors and defined as mTCDGs ( Figure 3A ). Specifically, ACAT1, ADGRF5, BBOX1, CYS1, DDC, GLYAT, HAO2, NPR3, PDK4, PLPP1, RCAN2, SGMS2, and TSPAN7 were relatively low expressed, while C4orf48, CCNB1, CDK1, CXCL1, NFKBIZ, PSAT1, PTTG1, S100A11, SAA1, SLC25A37, and TIMP1 exhibited elevated expression levels in metastatic RCC tumors compared to primary RCC tumors. Additionally, differential expression of 44 mTFs was also visualized ( Figure S14C, D ). Eventually, we employed GSVA to reveal the differential enrichment of 50 cancer hallmarks ( Figure S14E, F ). Specifically, unfolded protein response and DNA repair were notably up-regulated, while KRAS signaling dn and pancreas beta cells were notably down-regulated in metastatic RCC tumors.

Construction and validation of the multivariable prognostic model

Followingly, we characterized the enrichment of the 32 mTCDGs in various malignant pathways concerning the occurrence and progression of RCC by ORA ( Figure S15A ). Moreover, we conducted univariable Cox regression analysis to show their prognostic relevance in a forest plot ( Figure S15B ). Subsequently, we exerted Lasso regression to determine the optimal value of the parameter λ = 10 ( Figure 3B ). Hence, 10 selected mTCDGs, including S100A11, RHOB, CYS1, NFKBIZ, CXCL2, NPR3, PDK4, SLC25A37, PSAT1, and CDK1, were included to construct the multivariable prognostic model. In the subsequent cross-validation, differential expression patterns of the 10 mTCDGs were illustrated in all, train, and test sets ( Figures 3C and S15C). Specifically, S100A11, NFKBIZ, CXCL2, SLC25A37, PSAT1, and CDK1 exhibited relatively high expression, while RHOB, CYS1, NPR3, and PDK4 were relatively low expressed in high-risk subgroup compared to low-risk subgroup. Moreover, the distribution of risk scores and survival status was also shown ( Figure S15D, E ). Besides, K-M survival curves all showed significantly worse prognosis ( p < 0.001) in high-risk subgroup than low-risk subgroup ( Figure S15F ). What’s more, ROC curves revealed a decent accuracy of the model in all (AUC = 0.745), train (AUC = 0.752), and test (AUC = 0.757) sets ( Figure 3D ).

The independent prognostic prediction efficacy of the risk score was then validated by univariable Cox regression ( Figure 3E , p < 0.001, hazard ratio (HR) = 39.495, 95% confidence interval (CI) = 15.818–98.626) and multivariable Cox regression ( Figure 3F , p < 0.001, HR = 15.160, 95% CI = 5.658–40.622) analyses, taking the influence of confounding factors including age, M stage, N stage, and T stage into account.

Subsequently, clinical characteristics were then utilized to show differences between high-risk and low-risk subgroups by Chi-square tests ( Figure 3G ). Significantly, high-risk subgroup might be related to a higher probability of death ( p < 0.001), higher levels of T ( p < 0.001), N ( p < 0.05), M ( p < 0.001) stages, and unfavorable cancer status ( p < 0.001). Particularly, the significant differences ( p < 0.001) across M0, M1, and MX stages between high-risk and low-risk subgroups were highlighted in a cross table ( Figure 3H ). High-risk subgroup appeared to have more M1 and less M0 in comparison to low-risk subgroup, validating the metastasis-related feature of the prognostic model. Overall, the multivariable prognostic model exhibited remarkable clinical implications. Last but not least, GSEA revealed the top five enriched GOBP pathways ( Figure S16A, B ), KEGG pathways ( Figure S16C, D ), and cancer hallmarks ( Figure S16E, F ) in high-risk and low-risk subgroups, showing their distinct biological features.

Immunological relevance of the multivariable prognostic model

Firstly, we utilized the CIBERSORT algorithm to reveal differential infiltration of 22 immune cell types between high-risk and low-risk subgroups ( Figure S17A, B ). Moreover, significant differences in OS between high-infiltration and low-infiltration of immune cells were explored by K-M survival analysis ( Figure S17C ). Taken together, high-risk subgroup was found to exhibit a higher infiltration fraction of plasma cells, T cells CD4 memory activated, T cells follicular helper, Tregs, macrophages M0, and mast cells activating, while a lower infiltration fraction of B cells naïve, T cells CD4 memory resting, NK cells resting, monocytes, macrophages M1, macrophages M2, and mast cells resting, indicating a poor prognosis.

Additionally, we conducted ssGSEA to unveil the differential enrichment of 29 immune functions between high-risk and low-risk subgroups ( Figure S18A ). Also, significant differences in OS between high-enrichment and low-enrichment of immune functions were explored by K-M survival analysis ( Figure S18B ). Collectively, high-risk subgroup exhibited a higher enrichment of most immune functions, including aDCs, APC co-inhibition, APC co-stimulation, B cells, CCR, CD8 + T cells, check-point, cytolytic activity, DCs, HLA, inflammation-promoting, macrophages, MHC class I, NK cells, parainflammation, pDCs, T cell co-inhibition, T cell co-stimulation, T helper cells, Tfh, Th1 cells, Th2 cells, TIL, and type I IFN response, predominantly suggesting an unfavorable prognosis.

In addition, significant differences ( p < 0.001) across C1–C6 immune subtypes between high-risk and low-risk subgroups were identified ( Figure S19A ). Intriguingly, high-risk subgroup might be associated with more subtypes C4 and C6, but fewer subtype C3 in comparison to low-risk subgroup, indicating a terrible prognosis. Finally, TIDE analysis was performed to investigate differences in immunotherapeutic response characteristics between high-risk and low-risk subgroups ( Figure S19B ). Results showed a higher TIDE score, dysfunction score, exclusion score, MDSC, and CAF in high-risk subgroup than low-risk subgroup, suggesting a potentially positive immunotherapeutic response. Taken together, the multivariable prognostic model showed significant immune and immunotherapeutic correlations.

Regulatory network construction and inhibitor prediction for three RCC-TCDC subtypes

Regulatory networks were constructed by integrating upstream mTFs, key mTCDGs, downstream cancer hallmarks, infiltrating immune cells, enriched immune functions, and functioning proteomics in RCC-TCDC subtypes, including RCC-LD, RCC-MD, and RCC-HD, respectively, delineating the underlying molecular mechanisms ( Figure 4A ). Moreover, all significantly correlated components within the regulatory networks of three RCC-TCDC subtypes were also displayed in Figure 4B . Additionally, comprehensive GO and KEGG functional enrichment results were shown in Figure S20A–D . Moreover, to provide preliminary insights into the potential transcriptional regulation relationships, we utilized ATAC-seq profiles to detect the chromatin accessibility at the loci of key mTCDGs of each subtype ( Figure 4C ). Intriguingly, multiple binding peaks were identified at the chromosomal positions of SH3BGRL3, S100A11, CYS1, TJP2, CDK1, and PTTG1. What’s more, we further predicted the top 10 candidate inhibitors for RCC-LD, RCC-MD, and RCC-HD, respectively ( Figure 4D ). However, the underlying molecular mechanisms and these small molecular compounds holding potential for tailored treatment of RCC patients in each RCC-TCDC subtype warranted further in-depth investigation.

IHC staining and scoring for classifying patients with kidney neoplasms in Xinhua cohort based on RCC-TCDC

Following the comprehensive bioinformatic analysis that constructed and validated RCC-TCDC, we conducted a robust clinical validation involving 370 patients diagnosed with kidney neoplasms at Xinhua hospital ( Figure 5A ). The baseline characteristics of Xinhua cohort were presented in Table 2 and Figure S21A . In detail, our cohort comprised 33 patients aged 21–40, 194 patients aged 41–60, and 143 patients aged 61–88. Moreover, there were 114 female and 256 male patients. Regarding histopathological classification, the majority was KIRC (328), following by KIRP (18), KICH (8), and others (16). In addition, excluding 71 patients lost in follow-up, 276 patients were alive at the time of follow-up, with 23 patients deceased. Furthermore, 42 patients developed progression or death, while 257 patients did not. Afterwards, significant prognostic relevance was found in PFS censor, histopathological classification, stages, T stages, and Fuhrman nuclear stages by K-M survival curves ( Figure S21B–G ).

Subsequently, the previously identified differential expression patterns of the five selected markers, S100A11, SH3BGRL3, TSPAN7, AQP1, and NPR3, across three RCC-TCDC subtypes were shown in Figure 5B . Thereafter, we performed IHC staining and scoring of these markers on TMAs of kidney neoplasms retrieved from the 370 patients. Following this, we illustrated the significantly differential IHC scores of the five markers across three RCC-TCDC subtypes in Figure S22A . More importantly, a series of representative IHC images visualizing the differential staining results of S100A11, SH3BGRL3, TSPAN7, AQP1, and NPR3 in each RCC-TCDC subtype were displayed in Figure 5C . Specifically, S100A11 and SH3BGRL3 exhibited relatively high staining in RCC-LD, while medium or low staining in RCC-MD and RCC-HD. Additionally, TSPAN7 demonstrated relatively high staining in RCC-MD, medium or low staining in RCC-LD and RCC-HD. Moreover, AQP1 and NPR3 showed relatively high staining in RCC-HD, medium, low, or not detected staining in RCC-LD and RCC-HD. Furthermore, high expression levels of S100A11, SH3BGRL3, and TSPAN7 all revealed significantly worse prognoses for both OS ( p < 0.001, p = 0.004, and p = 0.009) and PFS ( p < 0.001, p < 0.001, and p = 0.009) in K-M survival curves ( Figure S22B–D ). Moreover, up-regulation of AQP1 was linked with a favorable prognosis for both OS ( p < 0.001) and PFS ( p < 0.001) ( Figure S22E ). Last but not least, low expression levels of NPR3 exhibited modest significance for longer OS ( p = 0.043) ( Figure S22F ).

Comprehensive statistical analyses for validating the clinical relevance of RCC-TCDC in Xinhua cohort

Chi-square tests were utilized to preliminarily reveal the clinical relevance of the RCC-TCDC in Xinhua cohort ( Figure 6A ). Intriguingly, a variety of clinical characteristics exhibited notable correlations with different RCC-TCDC subtypes. Specifically, significant differences in gender ( p = 0.03, Figure 6B ), stages ( p = 0.006, Figure 6C ), both T stages and detailed T stages ( p = 0.009 and p = 0.03, Figure 6D ), along with KIRC and histopathological classification ( p = 0.03 and p = 0.02, Figure 6E ) were identified across three RCC-TCDC subtypes, indicating that RCC-LD might be associated with a higher proportion of male, stage 3, T3 stage, and non-KIRC patients, and a lower proportion of female, stage 1, T1 stage, and KIRC patients. Regarding survival data, significant differences were revealed in OS censor and classification ( p = 0.02 and p < 0.001, Figure 6F ), along with PFS censor and classification ( p = 0.009 and p < 0.001, Figure 6G ) across three RCC-TCDC subtypes, suggesting that RCC-LD was correlated with a higher probability of death, progression, and a shorter OS and PFS.

Subsequently, we exerted K-M survival analysis to reveal significant differences ( p < 0.001) for OS across three RCC-TCDC subtypes ( Figure 7A ). Similar to previous bioinformatic results, RCC-LD exhibited the worst prognosis, RCC-HD represented the most favorable prognosis, with RCC-MD falling in between. Furthermore, we integrated age, gender, histopathological classification, and stage with RCC-TCDC into multivariable Cox regression analysis to identify RCC-LD as a significant risk factor ( p = 0.015, HR = 14.0, 95%CI = 1.67–117.8) for OS in patients with RCC, taking RCC-HD as a reference ( Figure 7B ). Likewise, significant differences ( p < 0.001) for PFS across three RCC-TCDC subtypes were also determined, with RCC-LD, RCC-MD, and RCC-HD exhibiting the worst, intermediate, and best prognosis, respectively ( Figure 7C ). Moreover, multivariable Cox regression analysis also revealed RCC-LD as a significant risk factor ( p = 0.010, HR = 4.0, 95%CI = 1.39–11.7) for PFS in RCC patients ( Figure 7D ). Taken together, RCC-TCDC was identified as an independent prognostic factor for RCC patients with robust validation.

Following that, nomograms were established for visualization ( Figure S23A, B ). Comprehensive model diagnostics were then performed for both multivariable Cox regression models. Initially, we observed a random distribution of residuals in residual plots ( Figure 7E ). Additionally, the modest predictive efficacy of the models was manifested by the calibration curves ( Figure 7F ). Moreover, robust AUC values varying with time for the ROC of the multivariable Cox regression models for OS and PFS were further displayed ( Figure 7G ). Specifically, ROC curves for predicting 1-year (AUC = 0.983), 2-year (AUC = 0.754), 3-year (AUC = 0.826), 4-year (AUC = 0.947), and 5-year (AUC = 0.940) OS, as well as 1-year (AUC = 0.837), 2-year (AUC = 0.733), 3-year (AUC = 0.778), 4-year (AUC = 0.831), and 5-year years (AUC = 0.774) PFS, were presented. In addition, risk scores were then calculated ( Figure S23C, D ). Favorable outcomes were obtained by risk curves, scatter plots, and K-M survival curves. Last but not least, Schoenfeld tests were conducted to validate the proportional hazards assumption for the multivariable Cox regression models ( Figure S23E, F ), revealing Schoenfeld residuals changing independent of time for each co-variate. Overall, a schematic diagram outlined the potential applications and future explorations of RCC-TCDC by designing a robust diagnostic kit and investigating tailored therapeutic options in Figure 7H .

Discussion

RCC is a heterogeneous malignancy with diverse genetic and molecular landscapes of tumor cells. There is an urgent need to develop molecular classification for RCC patients based on tumor cell differentiation for risk stratification and tailored treatment. The current study was initiated by a comprehensive bioinformatic analysis. By integrated analysis of scRNA-seq and bulk RNA-seq profiles, we determined 32 TCDGs and constructed RCC-TCDC, characterizing three subtypes with distinct prognosis: RCC-LD (high risk), RCC-MD (medium risk), and RCC-HD (low risk). Further exploration delved into the multi-omics relevance of the RCC-TCDC. Moreover, we established a multivariable prognostic model using 10 mTCDGs. In addition, we constructed regulatory networks and predicted candidate inhibitors for three RCC-TCDC subtypes. More importantly, to bolster the study’s robustness, we employed a large-sample pathological specimen and a retrospective clinical study including 370 patients with kidney neoplasms in Xinhua cohort for validation. After IHC staining and scoring of five specific markers on tumor tissues, patients were classified into RCC-LD (S100A11 + SH3BGRL3+), RCC-MD (TSPAN7+), and RCC-HD (AQP1 + NPR3+) subtypes according to RCC-TCDC. Finally, we obtained a series of significant results through Chi-square test, K-M survival analysis, multivariable Cox regression analysis, nomogram construction, and model diagnostics. Specifically, RCC-TCDC was identified as an independent prognostic factor for RCC patients, and RCC-LD was revealed to be significantly associated with the worst OS and PFS. Consequently, the clinical and prognostic relevance of RCC-TCDC was comprehensively validated. In fact, the majority of existing molecular classifications or diagnostic biomarkers of RCC with significant therapeutic implications were established by studies with large cohorts and bulk RNA-seq profiles [ 18–21 , 33 ]. However, there had been no robust gene classification system based on scRNA-seq data, and our study bridged this gap by constructing a novel molecular classification based on tumor cell differentiation and prognosis for RCC patients.

The high-risk RCC-TCDC subtype, RCC-LD (S100A11 + SH3BGRL3+), exhibited the worse OS in the TCGA-KIRC cohort and a higher frequency of advanced stage and T stage, as well as a shorter OS and PFS in Xinhua cohort. Moreover, elevated S100A11 and SH3BGRL3 expression showed significantly unfavorable prognosis at both the transcriptomic level in TCGA-KIRC cohort (OS, p < 0.001) and at the proteomic level in Xinhua cohort (OS, p = 0.004 and PFS, p < 0.001), serving as two risk predictors. Notably, S100A11, or S100 Calcium Binding Protein A11, had been extensively explored as an oncogene and a prognostic marker in numerous types of tumors, including RCC, gastric cancer, colorectal cancer, cervical cancer, and intrahepatic cholangiocarcinoma [ 34–39 ]. Specifically, previous research identified S100A11 as a potential prognostic indicator for KIRC and a specific protein marker for KIRP [ 40 , 41 ]. Additionally, regarding molecular mechanisms, an early study highlighted the essential role of S100A11 in efficient plasma membrane repair and aggressive cancer cell survival [ 42 ]. Moreover, potential regulation of S100A11 was found on proliferation, migration, and invasion of RCC cancer cells through the EGFR/AKT signaling pathway and E-cadherin [ 43 ]. In addition, recent findings have indicated that LINC00997 contributed to metastasis by regulating S100A11 expression in KIRC [ 44 ].

The other specific marker, SH3BGRL3, or SH3 domain-binding glutamate-rich protein like 3, was discovered to be a tumor necrosis factor (TNF)-α inhibitory protein, which has the ability to protect tumor cells against TNF-mediated apoptosis [ 45 , 46 ]. Recent investigations have also shed light on its oncogenic properties across various cancers. It has been implicated in promoting breast cancer cell migration, facilitating glioblastoma tumorigenesis, and regulating the progression of acute myeloid leukemia [ 47–49 ]. Intriguingly, a remarkable study has highlighted SH3BGRL3 as a potential prognostic biomarker for urothelial carcinoma [ 50 ]. In this study, SH3BGRL3 was uncovered to be positively correlated with higher histological grading, muscle invasiveness of urothelial carcinoma, and increased progression risk in nonmuscle-invasive bladder cancer patients. Moreover, the in vitro assays further displayed SH3BGRL3’s role in driving the epithelial-mesenchymal transition (EMT), cell migration, and proliferation of urothelial carcinoma. More importantly, a notable study focusing on KIRC also identified SH3BGRL3 as an independent risk factor for predicting progression and metastasis [ 51 ]. Elevated expression levels of SH3BGRL3 were observed in tumor tissues compared to normal tissues, and its up-regulation was associated with multiple aggressive clinicopathological features. Subsequent in vitro and in vivo experiments unveiled its involvement in promoting the growth and metastasis of KIRC through the EGFR/AKT signaling pathway.

The medium-risk RCC-TCDC subtype, RCC-MD (TSPAN7+), showed moderate prognostic outcomes in both the TCGA-KIRC and Xinhua cohorts. Interestingly, the contradictory prognostic implications of up-regulated TSPAN7 expression were observed between transcriptomic and proteomic levels. While TSPAN7 over-expression was associated with significantly longer OS in the TCGA-KIRC cohort ( p < 0.001), it was correlated with shorter OS and PFS in Xinhua cohort (both p = 0.009). TSPAN7, or Tetraspanin 7, belongs to the 4-transmembrane protein superfamily that plays a pivotal role in regulating tumor initiation, promotion, and metastasis [ 52 ]. Indeed, it has been extensively investigated across various cancers as a notably double-edged sword with both pro-tumor and anti-tumor functions [ 53–55 ]. On the one hand, an early study conducted by D. Wuttig et al. first identified TSPAN7 as a prognostic biomarker for KIRC in 2012 [ 23 ]. Specifically, high expression of TSPAN7 in tumor tissues was significantly associated with favorable prognosis at the transcriptomic level by genome-wide expression analyses and qPCR. However, this association was not replicated at the proteomic level by IHC staining on TMA. Nonetheless, TSPAN7-positive vessels in both cores exhibited modest significance for favorable disease-free survival ( p = 0.050) and tumor-specific survival ( p = 0.042). Still, the limited number of samples ( n = 106) and time efficiency also warranted further consideration. Additionally, up-regulation of TSPAN7 was uncovered in patients with multiple myeloma, and it might inhibit tumor development [ 56 ]. Moreover, low expression of TSPAN7 was found in bladder cancer (BCa) tissues and cell lines, which was a positive indicator of advanced tumor grade, tumor stage, and reduced OS at both transcriptomic and proteomic levels [ 57 ]. Furthermore, researchers elucidated the inhibition of TSPAN7 on BCa cell growth via the PTEN/PI3K/AKT signaling pathway. On the other hand, TSPAN7 also exerted pro-tumor effects in multiple cancers. In non-small-cell lung cancer, TSPAN7 was unveiled to be an adverse prognostic predictor that might facilitate tumor cell migration and proliferation through EMT [ 58 ]. In addition, in osteosarcoma, TSPAN7 exhibited enhancement of tumor cell invasion and metastasis by induction of EMT and activation of the FAK-Src-Ras-ERK1/2 signaling pathway [ 59 ].

The low-risk RCC-TCDC subtype, RCC-HD (AQP1 + NPR3+), demonstrated the best prognosis in both the TCGA-KIRC and Xinhua cohorts. Notably, high expression of AQP1 and NPR3 both exhibited extended OS at the transcriptomic level (both p < 0.001). However, at the proteomic level, while over-expression of AQP1 showed significantly longer OS and PFS (both p < 0.001), NPR3 displayed modestly significant association with shorter OS ( p = 0.043) and no significance in PFS ( p = 0.063). AQP1, or aquaporin 1, is a member of the water channel proteins existing in almost all human organs and tissues [ 60 ]. Its crucial role has been implicated in the early detection, prognostic prediction, tumor growth, angiogenesis, invasion, and metastasis in various cancers, particularly those originating from organs with essential water permeability, such as the kidney [ 61–64 ]. More intriguingly, it serves as a differentiation marker for proximal renal tubule cells, the origin of KIRC tumor [ 65 ]. An early investigation in the last century indicated transcriptional level of AQP1 as a potential favorable prognostic predictor in RCC patients following nephrectomy [ 66 ]. A subsequent TMA study further revealed positive correlation between high AQP1 expression and lower tumor grade in KIRC at the proteomic level [ 67 ]. Following that, a remarkable study presented that AQP1 was significantly higher in KIRC patients presenting non-typical symptoms, or exhibiting smaller size, lower grade, or lower stage tumors, or absence of microvascular invasion [ 68 ]. More importantly, it also showed that high AQP1 expression was associated with better prognosis in cancer-specific survival (CSS) and cancer-free survival (CFS) across all RCC patients, as well as in CSS among advanced metastatic RCC patients. Notably, these results were further supported by a recent large-sample IHC study, which underscored a significant correlation between AQP over-expression and several favorable clinicopathological characteristics, including younger age, lower TNM stage and WHO/ISUP nuclear grade, as well as the absence of distant metastasis [ 69 ]. Furthermore, up-regulation of AQP1 also served as a significantly independent predictor of both longer OS and PFS in KIRC. Additionally, elevated AQP1 expression was found to suppress invasion, migration, EMT, proliferation, apoptosis, and pyroptosis in Wilm’s tumor cells [ 70 , 71 ].

The other specific marker, NPR3, or natriuretic peptide receptor 3, which mediates natriuretic peptide degradation, was reported to act as an anti-tumor factor in multiple cancers. Specifically, a previous study revealed that MRCCAT1, a long non-coding RNA (lncRNA), contributed to KIRC cell proliferation, migration, and invasion in vitro and KIRC metastasis in vivo , by the inhibition of NPR3 and subsequent activation of the p38-MAPK signaling pathway [ 72 ]. Likewise, in hepatocellular carcinoma (HCC), lncRNA FENDRR was revealed to inhibit the viability and induce the apoptosis of HCC cells by targeting miR-362-5p degradation to promote NPR3 and suppress the p38-MAPK signaling pathway [ 73 ]. Moreover, NPR3 was uncovered to serve as a novel tumor suppressor in osteosarcoma, which was regulated by TF POU2F1 to inhibit the PI3K/AKT signaling pathway, thus inhibiting tumor cell growth [ 74 ]. Last but not least, a recent bioinformatic study identified NPR3 as a sarcomatoid differentiation-related gene whose up-regulation predicted a better prognosis in KIRC patients [ 75 ].

Taken together, existing evidence exhibited consistent clinical significance and mechanistic relevance of the selected five markers for three distinct RCC-TCDC subtypes, RCC-LD (S100A11 + SH3BGRL3+, high risk), RCC-MD (TSPAN7+, medium risk), and RCC-HD (AQP1 + NPR3+, low risk) with our own results.

However, despite the robust construction of RCC-TCDC for RCC patients, our study had certain limitations. First, while we constructed the novel molecular classification system using publicly available scRNA-seq and bulk RNA-seq profiles, the inclusion of our own RCC samples for sequencing and subsequent in silico analysis could have strengthened our findings. Second, despite the inclusion of a substantial number of pathological samples and patients ( n = 370) in Xinhua cohort for validation, it remains a single-center study, thereby limiting its generalizability. Third, although Xinhua cohort provided valuable clinicopathological features and follow-up information, there is room for further enhancement of data integrity. Forth, comprehensive in vitro and in vivo experiments are warranted to further validate the underlying regulatory mechanisms and candidate inhibitors for three RCC-TCDC subtypes.

Conclusions

Based on an integrated bioinformatics analysis, we identified 32 TCDGs and developed a novel molecular classification system, RCC-TCDC, classifying TCGA-KIRC patients into RCC-LD (S100A11 + SH3BGRL3+, high risk), RCC-MD (TSPAN7+, medium risk), and RCC-HD (AQP1 + NPR3+, low risk). Notably, we then employed a retrospective clinical study validating RCC-LD as an independent risk factor for RCC patients by univariable and multivariable analyses. Specifically, ROC curves for predicting 1-year (AUC = 0.983), 2-year (AUC = 0.754), 3-year (AUC = 0.826), 4-year (AUC = 0.947), and 5-year (AUC = 0.940) OS, as well as 1-year (AUC = 0.837), 2-year (AUC = 0.733), 3-year (AUC = 0.778), 4-year (AUC = 0.831), and 5-year years (AUC = 0.774) PFS all showed remarkable accuracy. In the near future, we will design a diagnostic kit based on RCC-TCDC and apply for a patent. Hopefully, we will subsequently conduct comprehensive in vitro and in vivo experiments to explore the clinical implications of identified underlying regulatory mechanisms and candidate inhibitors for three RCC-TCDC subtypes.

Supplementary Material

Supplemental Material

Acknowledgments

We thank all databases accessed in this study, including the GEO and TCGA database for allowing us to use their data. Conception/design: Runzhi Huang, Xingang Cui; collection and/or assembly of data: Keqin Dong, Xiuwu Pan; data analysis and interpretation: Yifan Liu, Yuntao Yao, Bingnan Lu, Guo Ji, Runzhi Huang; manuscript writing: Yifan Liu, Yuntao Yao, Bingnan Lu, Lei Wang, Haoyu Zhang, Zihui Zhao, Xinyue Yang; manuscript revising: Runzhi Huang, Wang Zhou, Xiuwu Pan, Xingang Cui; final approval of manuscript: Yifan Liu, Keqin Dong, Yuntao Yao, Bingnan Lu, Lei Wang, Guo Ji, Haoyu Zhang, Zihui Zhao, Xinyue Yang, Runzhi Huang, Wang Zhou, Xiuwu Pan, Xingang Cui; we confirm that all authors have read and approved the final work.

Ethics approval and consent to participate

The study was approved by the Ethics Committee of Xinhua Hospital Affiliated to Shanghai Jiao Tong University School of Medicine (XHEC-C-2021-145-1). Informed consent to participate was obtained by patients. This study adheres fully to the principles outlined in the Declaration of Helsinki.

Consent for publication

Written informed consent from all participants upon their admission to the hospital. For participants who were minors, written informed consent was obtained from their legal guardians. Additionally, assent was obtained from minors where appropriate.

Disclosure statement

The authors declare that they have no competing interests.

Data availability statement

The authors agree to share the data and materials supporting the results or analyses presented in this paper upon reasonable request, unless restricted by ethical, privacy, or security concerns.

References

1 Sung H , Ferlay J , Siegel RL , et al. Global Cancer Statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries . CA Cancer J Clin . 2021 ; 71 ( 3 ): 209 – 249 . doi: 10.3322/caac.21660 . 33538338

2 Escudier B , Porta C , Schmidinger M , et al. Renal cell carcinoma: ESMO clinical practice guidelines for diagnosis, treatment and follow-up . Ann Oncol . 2016 ; 27 ( suppl 5 ): v58 – v68 . doi: 10.1093/annonc/mdw328 . 27664262

3 Patard J-J , Leray E , Rioux-Leclercq N , et al. Prognostic value of histologic subtypes in renal cell carcinoma: a multicenter experience . J Clin Oncol . 2005 ; 23 ( 12 ): 2763 – 2771 . doi: 10.1200/JCO.2005.07.055 . 15837991

4 Choueiri TK , Motzer RJ. Systemic therapy for metastatic renal-cell carcinoma . N Engl J Med . 2017 ; 376 ( 4 ): 354 – 366 . doi: 10.1056/NEJMra1601333 . 28121507

5 Hsieh JJ , Purdue MP , Signoretti S , et al. Renal cell carcinoma . Nat Rev Dis Primers . 2017 ; 3 : 17009 . doi: 10.1038/nrdp.2017.9 . 28276433 PMC5936048

6 Srigley JR , Delahunt B , Eble JN , et al. The International Society of Urological Pathology (ISUP) vancouver classification of renal neoplasia . Am J Surg Pathol . 2013 ; 37 ( 10 ): 1469 – 1489 . doi: 10.1097/PAS.0b013e318299f2d1 . 24025519

7 Turajlic S , Swanton C , Boshoff C. Kidney cancer: the next decade . J Exp Med . 2018 ; 215 ( 10 ): 2477 – 2479 . doi: 10.1084/jem.20181617 . 30217855 PMC6170181

8 Moch H , Amin MB , Berney DM , et al. The 2022 World Health Organization classification of tumours of the urinary system and male genital organs-part A: Renal, ­penile, and testicular tumours . Eur Urol . 2022 ; 82 ( 5 ): 458 – 468 . doi: 10.1016/j.eururo.2022.06.016 . 35853783

9 Ljungberg B , Albiges L , Abu-Ghanem Y , et al. European association of urology guidelines on renal cell carcinoma: the 2022 update . Eur Urol . 2022 ; 82 ( 4 ): 399 – 410 . doi: 10.1016/j.eururo.2022.03.006 . 35346519

10 Yao Y , et al. Tracing the evolving dynamics and research hotspots in the kidney neoplasm and nephron sparing surgery field from the past to the new era . Cancer Med . 2024 ; 13 ( 12 ): e7336 . 39651783 10.1002/cam4.7336 PMC11192648

11 Johansson ME , Axelson H. Kidney cancer . Semin Cancer Biol . 2013 ; 23 ( 1 ): 1 – 2 . doi: 10.1016/j.semcancer.2012.05.006 . 22664338

12 Gerlinger M , Rowan AJ , Horswell S , et al. Intratumor heterogeneity and branched evolution revealed by multiregion sequencing . N Engl J Med . 2012 ; 366 ( 10 ): 883 – 892 . doi: 10.1056/NEJMoa1113205 . 22397650 PMC4878653

13 McGranahan N , Swanton C. Biological and therapeutic impact of intratumor heterogeneity in cancer evolution . Cancer Cell . 2015 ; 27 ( 1 ): 15 – 26 . doi: 10.1016/j.ccell.2014.12.001 . 25584892

14 Hart IR , Easty D. Tumor cell progression and differentiation in metastasis . Semin Cancer Biol . 1991 ; 2 ( 2 ): 87 – 95 . 1912522

15 Vicente-Dueñas C , Hauer J , Ruiz-Roca L , et al. Tumoral stem cell reprogramming as a driver of cancer: theory, biological models, implications in cancer therapy . Semin Cancer Biol . 2015 ; 32 : 3 – 9 . doi: 10.1016/j.semcancer.2014.02.001 . 24530939

16 Dagogo-Jack I , Shaw AT. Tumour heterogeneity and ­resistance to cancer therapies . Nat Rev Clin Oncol . 2018 ; 15 ( 2 ): 81 – 94 . doi: 10.1038/nrclinonc.2017.166 . 29115304

17 Zhang Y , Narayanan SP , Mannan R , et al. Single-cell analyses of renal cell cancers reveal insights into tumor microenvironment, cell of origin, and therapy response . Proc Natl Acad Sci U S A . 2021 ; 118 ( 24 ):e2103240118. doi: 10.1073/pnas.2103240118 . PMC8214680 34099557

18 Hakimi AA , Voss MH , Kuo F , et al. Transcriptomic profiling of the tumor microenvironment reveals distinct subgroups of clear cell renal cell cancer: data from a randomized phase III trial . Cancer Discov . 2019 ; 9 ( 4 ): 510 – 525 . doi: 10.1158/2159-8290.CD-18-0957 . 30622105 PMC6697163

19 Motzer RJ , Banchereau R , Hamidi H , et al. Molecular subsets in renal cancer determine outcome to checkpoint and angiogenesis blockade . Cancer Cell . 2020 ; 38 ( 6 ): 803 – 817.e4 . doi: 10.1016/j.ccell.2020.10.011 . 33157048 PMC8436590

20 Vano Y-A , Elaidi R , Bennamoun M , et al. Nivolumab, nivolumab-ipilimumab, and VEGFR-tyrosine kinase ­inhibitors as first-line treatment for metastatic clear-cell renal cell carcinoma (BIONIKK): a biomarker-driven, open-label, non-comparative, randomised, phase 2 trial . Lancet Oncol . 2022 ; 23 ( 5 ): 612 – 624 . doi: 10.1016/S1470-2045(22)00128-0 . 35390339

21 Büttner FA , Winter S , Stühler V , et al. A novel molecular signature identifies mixed subtypes in renal cell carcinoma with poor prognosis and independent response to immunotherapy . Genome Med . 2022 ; 14 ( 1 ): 105 . doi: 10.1186/s13073-022-01105-y . 36109798 PMC9476269

22 Barrett T , Wilhite SE , Ledoux P , et al. NCBI GEO: archive for functional genomics data sets–update . Nucleic Acids Res . 2013 ; 41 ( Database issue ): D991 – 5 . doi: 10.1093/nar/gks1193 . 23193258 PMC3531084

23 Wuttig D , Zastrow S , Füssel S , et al. CD31, EDNRB and TSPAN7 are promising prognostic markers in clear-cell renal cell carcinoma revealed by genome-wide expression analyses of primary tumors and ­metastases . Int J Cancer . 2012 ; 131 ( 5 ): E693 – 704 . doi: 10.1002/ijc.27419 . 22213152

24 Robinson DR , Wu Y-M , Lonigro RJ , et al. Integrative clinical genomics of metastatic cancer . Nature . 2017 ; 548 ( 7667 ): 297 – 303 . doi: 10.1038/nature23306 . 28783718 PMC5995337

25 Zheng R , Wan C , Mei S , et al. Cistrome data browser: expanded datasets and new tools for gene regulatory analysis . Nucleic Acids Res . 2019 ; 47 ( D1 ): D729 – D735 . doi: 10.1093/nar/gky1094 . 30462313 PMC6324081

26 Butler A , Hoffman P , Smibert P , et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species . Nat Biotechnol . 2018 ; 36 ( 5 ): 411 – 420 . doi: 10.1038/nbt.4096 . 29608179 PMC6700744

27 Wolf FA , Angerer P , Theis FJ. SCANPY: large-scale single-cell gene expression data analysis . Genome Biol . 2018 ; 19 ( 1 ): 15 . doi: 10.1186/s13059-017-1382-0 . 29409532 PMC5802054

28 Newman AM , Liu CL , Green MR , et al. Robust enumeration of cell subsets from tissue expression profiles . Nat Methods . 2015 ; 12 ( 5 ): 453 – 457 . doi: 10.1038/nmeth.3337 . 25822800 PMC4739640

29 Charoentong P , Finotello F , Angelova M , et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade . Cell Rep . 2017 ; 18 ( 1 ): 248 – 262 . doi: 10.1016/j.celrep.2016.12.019 . 28052254

30 Thorsson V , Gibbs DL , Brown SD , et al. The immune landscape of cancer . Immunity . 2018 ; 48 ( 4 ): 812 – 830.e14 . doi: 10.1016/j.immuni.2018.03.023 . 29628290 PMC5982584

31 Fu J , Li K , Zhang W , et al. Large-scale public data reuse to model immunotherapy response and resistance . Genome Med . 2020 ; 12 ( 1 ): 21 . doi: 10.1186/s13073-020-0721-z . 32102694 PMC7045518

32 Jiang P , Gu S , Pan D , et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy ­response . Nat Med . 2018 ; 24 ( 10 ): 1550 – 1558 . doi: 10.1038/s41591-018-0136-1 . 30127393 PMC6487502

33 Milella M , Rutigliano M , Lasorsa F , et al. The role of MUC1 in renal cell carcinoma . Biomolecules . 2024 ; 14 ( 3 ): 315 . doi: 10.3390/biom14030315 . 38540735 PMC10967890

34 Niu Y , Shao Z , Wang H , et al. LASP1-S100A11 axis promotes colorectal cancer aggressiveness by modulating TGF β /Smad signaling . Sci Rep . 2016 ; 6 ( 1 ): 26112 . doi: 10.1038/srep26112 . 27181092 PMC4867635

35 Koh SA , Lee KH. HGF-mediated S100A11 overexpression enhances proliferation and invasion of gastric cancer . Am J Transl Res . 2018 ; 10 ( 11 ): 3385 – 3394 . 30662594 PMC6291695

36 Meng M , Sang L , Wang X. S100 calcium binding protein A11 (S100A11) promotes the proliferation, migration and invasion of cervical cancer cells, and activates Wnt/ β -catenin signaling . Onco Targets Ther . 2019 ; 12 : 8675 – 8685 . doi: 10.2147/OTT.S225248 . 31695426 PMC6815786

37 Zhang M-X , Gan W , Jing C-Y , et al. S100A11 promotes cell proliferation via P38/MAPK signaling pathway in intrahepatic cholangiocarcinoma . Mol Carcinog . 2019 ; 58 ( 1 ): 19 – 30 . doi: 10.1002/mc.22903 . 30182496 PMC6587853

38 Zhang L , Zhu T , Miao H , et al. The calcium binding protein S100A11 and its roles in diseases . Front Cell Dev Biol . 2021 ; 9 : 693262 . doi: 10.3389/fcell.2021.693262 . 34179021 PMC8226020

39 Ji X , Qin X , Huang X , et al. S100A11: a potential carcinogen and prognostic marker that correlates with the immunosuppressive microenvironment in pan-cancer . J Cancer . 2023 ; 14 ( 1 ): 88 – 98 . doi: 10.7150/jca.78011 . 36605485 PMC9809332

40 Na CH , Hong JH , Kim WS , et al. Identification of protein markers specific for papillary renal cell carcinoma using imaging mass spectrometry . Mol Cells . 2015 ; 38 ( 7 ): 624 – 629 . doi: 10.14348/molcells.2015.0013 . 26062552 PMC4507028

41 Gabril M , Girgis H , Scorilas A , et al. S100A11 is a potential prognostic marker for clear cell renal cell carcinoma . Clin Exp Metastasis . 2016 ; 33 ( 1 ): 63 – 71 . doi: 10.1007/s10585-015-9758-6 . 26472670

42 Jaiswal JK , Lauritzen SP , Scheffer L , et al. S100A11 is ­required for efficient plasma membrane repair and survival of invasive cancer cells . Nat Commun . 2014 ; 5 : 3795 . doi: 10.1038/ncomms4795 . 24806074 PMC4026250

43 Liu L , Miao L , Liu Y , et al. S100A11 regulates renal carcinoma cell proliferation, invasion, and migration via the EGFR/Akt signaling pathway and E-cadherin . Tumour Biol . 2017 ; 39 ( 5 ): 1010428317705337 . doi: 10.1177/1010428317705337 . 28513300

44 Chang Y , Li N , Yuan W , et al. LINC00997, a novel long noncoding RNA, contributes to metastasis via regulation of S100A11 in kidney renal clear cell carcinoma . Int J Biochem Cell Biol . 2019 ; 116 : 105590 . doi: 10.1016/j.biocel.2019.105590 . 31442606

45 Berleth ES , Nadadur S , Henn AD , et al. Identification, characterization, and cloning of TIP-B1, a novel protein inhibitor of tumor necrosis factor-induced lysis . Cancer Res . 1999 ; 59 ( 21 ): 5497 – 5506 . 10554026

46 Berleth ES , Masso-Welch PA , Kazim LA , et al. Expression, tissue distribution, and cellular localization of the antiapoptotic TIP-B1 protein . J Leukoc Biol . 2001 ; 69 ( 6 ): 995 – 1005 . doi: 10.1189/jlb.69.6.995 . 11404387

47 Di Pisa F , et al. SH3BGRL3 binds to myosin 1c in a calcium dependent manner and modulates migration in the MDA-MB-231 cell line . BMC Mol Cell Biol . 2021 ; 22 ( 1 ): 41 . 34380438 10.1186/s12860-021-00379-1 PMC8356473

48 Nie Z , Cheng D , Pan C , et al. SH3BGRL3, transcribed by STAT3, facilitates glioblastoma tumorigenesis by activating STAT3 signaling . Biochem Biophys Res Commun . 2021 ; 556 : 114 – 120 . doi: 10.1016/j.bbrc.2021.03.165 . 33839406

49 Yang X , Wang Y , Rong S , et al. Gene SH3BGRL3 regulates acute myeloid leukemia progression through circRNA_0010984 based on competitive endogenous RNA mechanism . Front Cell Dev Biol . 2023 ; 11 : 1173491 . doi: 10.3389/fcell.2023.1173491 . 37397256 PMC10313326

50 Chiang C-Y , Pan C-C , Chang H-Y , et al. SH3BGRL3 protein as a potential prognostic biomarker for urothelial carcinoma: a novel binding partner of epidermal growth factor receptor . Clin Cancer Res . 2015 ; 21 ( 24 ): 5601 – 5611 . doi: 10.1158/1078-0432.CCR-14-3308 . 26286913

51 Yin L , Gao S , Shi H , et al. TIP-B1 promotes kidney clear cell carcinoma growth and metastasis via EGFR/AKT signaling . Aging (Albany NY) . 2019 ; 11 ( 18 ): 7914 – 7937 . doi: 10.18632/aging.102298 . 31562290 PMC6782011

52 Hemler ME. Tetraspanin proteins promote multiple cancer stages . Nat Rev Cancer . 2014 ; 14 ( 1 ): 49 – 60 . doi: 10.1038/nrc3640 . 24505619

53 Perot BP , Ménager MM. Tetraspanin 7 and its closest paralog tetraspanin 6: membrane organizers with key functions in brain development, viral infection, innate immunity, diabetes and cancer . Med Microbiol Immunol . 2020 ; 209 ( 4 ): 427 – 436 . doi: 10.1007/s00430-020-00681-3 . 32468130

54 Huang R , Sun H , Lin R , et al. The role of tetraspanins pan-cancer . iScience . 2022 ; 25 ( 8 ): 104777 . doi: 10.1016/j.isci.2022.104777 . 35992081 PMC9385710

55 Lu B , Liu Y , Yao Y , et al. Unveiling the unique role of TSPAN7 across tumors: a pan-cancer study incorporating retrospective clinical research and bioinformatic analysis . Biol Direct . 2024 ; 19 ( 1 ): 72 . doi: 10.1186/s13062-024-00516-8 . 39175035 PMC11340126

56 Cheong CM , Chow AWS , Fitter S , et al. Tetraspanin 7 (TSPAN7) expression is upregulated in multiple myeloma patients and inhibits myeloma tumour development in vivo . Exp Cell Res . 2015 ; 332 ( 1 ): 24 – 38 . doi: 10.1016/j.yexcr.2015.01.006 . 25637218

57 Yu X , Li S , Pang M , et al. TSPAN7 exerts anti-tumor effects in bladder cancer through the PTEN/PI3K/AKT pathway . Front Oncol . 2020 ; 10 : 613869 . doi: 10.3389/fonc.2020.613869 . 33489923 PMC7821430

58 Wang X , Lin M , Zhao J , et al. TSPAN7 promotes the ­migration and proliferation of lung cancer cells via epithelial-to-mesenchymal transition . Onco Targets Ther . 2018 ; 11 : 8815 – 8822 . doi: 10.2147/OTT.S167902 . 30588007 PMC6300375

59 Shao S , Piao L , Guo L , et al. Tetraspanin 7 promotes ­osteosarcoma cell invasion and metastasis by inducing EMT and activating the FAK-Src-Ras-ERK1/2 signaling pathway . Cancer Cell Int . 2022 ; 22 ( 1 ): 183 . doi: 10.1186/s12935-022-02591-1 . 35524311 PMC9074275

60 Magni F , Chinello C , Raimondo F , et al. AQP1 expression analysis in human diseases: implications for proteomic characterization . Expert Rev Proteom . 2008 ; 5 ( 1 ): 29 – 43 . doi: 10.1586/14789450.5.1.29 . 18282122

61 King LS , Agre P. Pathophysiology of the aquaporin ­water channels . Annu Rev Physiol . 1996 ; 58 ( 1 ): 619 – 648 . doi: 10.1146/annurev.ph.58.030196.003155 . 8815812

62 Moosavi MS , Elham Y. Aquaporins 1, 3 and 5 in different tumors, their expression, prognosis value and role as new therapeutic targets . Pathol Oncol Res . 2020 ; 26 ( 2 ): 615 – 625 . doi: 10.1007/s12253-019-00646-9 . 30927206

63 Moon CS , Moon D , Kang SK. Aquaporins in cancer biology . Front Oncol . 2022 ; 12 : 782829 . doi: 10.3389/fonc.2022.782829 . 35847914 PMC9278817

64 Liu Y , Lyu D , Yao Y , et al. The comprehensive potential of AQP1 as a tumor biomarker: evidence from kidney neoplasm cohorts, cell experiments and pan-cancer analysis . Hum Genomics . 2025 ; 19 ( 1 ): 15 . doi: 10.1186/s40246-025-00726-9 . 39988693 PMC11849320

65 Baer PC , Bereiter-Hahn J , Schubert R , et al. Differentiation status of human renal proximal and distal tubular epithelial cells in vitro: differential expression of characteristic markers . Cells Tissues Organs . 2006 ; 184 ( 1 ): 16 – 22 . doi: 10.1159/000096947 . 17190976

66 Takenawa J , Kaneko Y , Kishishita M , et al. Transcript levels of aquaporin 1 and carbonic anhydrase IV as predictive indicators for prognosis of renal cell carcinoma patients after nephrectomy . Int J Cancer . 1998 ; 79 ( 1 ): 1 – 7 . doi: 10.1002/(SICI)1097-0215(19980220)79:1<1::AID-IJC1>3.0.CO;2-5 . 9495349

67 Mazal PR , Stichenwirth M , Koller A , et al. Expression of aquaporins and PAX-2 compared to CD10 and cytokeratin 7 in renal neoplasms: a tissue microarray study . Mod Pathol . 2005 ; 18 ( 4 ): 535 – 540 . doi: 10.1038/modpathol.3800320 . 15502805

68 Huang Y , Murakami T , Sano F , et al. Expression of aquaporin 1 in primary renal tumors: a prognostic indicator for clear-cell renal cell carcinoma . Eur Urol . 2009 ; 56 ( 4 ): 690 – 698 . doi: 10.1016/j.eururo.2008.10.014 . 18930582

69 Lee S , Kim B , Jung M , et al. Loss of aquaporin-1 expression is associated with worse clinical outcomes in clear cell renal cell carcinoma: an immunohistochemical study . J Pathol Transl Med . 2023 ; 57 ( 4 ): 232 – 237 . doi: 10.4132/jptm.2023.06.17 . 37460397 PMC10369134

70 Liu H , Jin C , Yang X , et al. Identification of key genes and validation of key gene aquaporin 1 on Wilms’ ­tumor metastasis . PeerJ . 2023 ; 11 : e16025 . doi: 10.7717/peerj.16025 . 37904849 PMC10613441

71 Liu H , Jin C , Xia N , et al. Overexpression of aquaporin-1 plays a vital role in proliferation, apoptosis, and pyroptosis of Wilms’ tumor cells . J Cancer Res Clin Oncol . 2024 ; 150 ( 2 ): 85 . doi: 10.1007/s00432-024-05616-6 . 38334883 PMC10858134

72 Li J-K , Chen C , Liu J-Y , et al. Long noncoding RNA MRCCAT1 promotes metastasis of clear cell renal cell carcinoma via inhibiting NPR3 and activating p38-MAPK signaling . Mol Cancer . 2017 ; 16 ( 1 ): 111 . doi: 10.1186/s12943-017-0681-0 . 28659173 PMC5490088

73 Qian G , Jin X , Zhang L. LncRNA FENDRR upregulation promotes hepatic carcinoma cells apoptosis by targeting miR-362-5p Via NPR3 and p38-MAPK pathway . Cancer Biother Radiopharm . 2020 ; 35 ( 9 ): 629 – 639 . doi: 10.1089/cbr.2019.3468 . 32251605

74 Li S , Guo R , Peng Z , et al. NPR3, transcriptionally regulated by POU2F1, inhibits osteosarcoma cell growth through blocking the PI3K/AKT pathway . Cell Signal . 2021 ; 86 : 110074 . doi: 10.1016/j.cellsig.2021.110074 . 34229087

75 Zuo Y , Fu S , Zhao Z , et al. Sarcomatoid-associated gene risk index for clear cell renal cell carcinoma . Front Genet . 2022 ; 13 : 985641 . doi: 10.3389/fgene.2022.985641 . 36159988 PMC9493111