Unveiling the comorbidity hub: WT1 drives renal cancer progression in chronic kidney disease and confers sirolimus vulnerability

Introduction

Clear cell renal cell carcinoma (RCC; KIRC) accounts for the majority of RCC cases, comprising 75-85% of all RCC subtypes, and remains a leading cause of cancer-related mortality worldwide, with an estimated ~175,000 deaths annually (1,2). Global estimates suggest that >400,000 new RCC diagnoses are made each year and ~175,000 RCC-related deaths, an upward trend that has persisted for over a decade (3,4). Despite advances in targeted therapies and immune checkpoint inhibitors (5,6), outcomes for metastatic disease remain poor (7). A key reason lies in the marked biological heterogeneity of the tumor and its ability to rapidly develop resistance to treatment (8). What complicates matters further is that a number of patients with RCC also suffer from comorbid conditions, particularly chronic kidney disease (CKD). Studies report that 30-40% of patients with RCC have concomitant CKD (9), a combination that not only accelerates renal function decline but also limits tolerance to standard oncologic regimens (10). Clinically, this creates a difficult balancing act as aggressive immunotherapies may control tumor growth but risk worsening kidney damage, directly impacting quality of life and survival. Notably, patients with both KIRC and CKD experience a median survival reduction of 2-3 years compared with those with either condition alone (11,12).

At the molecular level, the interplay between RCC and CKD is poorly understood. Both diseases involve dysregulated immune pathways and metabolic disturbances (13), yet most prior research has treated them as separate entities. This siloed approach overlooks potential shared mechanisms, especially those involving immunogenic cell death (ICD), a form of regulated cell death that can stimulate antitumor immunity (14,15). ICD has received little attention in the context of RCC-CKD comorbidity, despite its theoretical relevance (16).

The aim of the present study was to address this gap. Rather than analyzing RCC or CKD in isolation, multi-omics data, from bulk and single-cell RNA sequencing to epigenetic and proteomic profiles, was integrated using machine learning to identify convergent molecular signatures. The present study aimed to identify shared molecular mechanisms linking CKD and RCC, with a focus on exploring the role of Wilms tumor 1 (WT1) in immune-metabolic crosstalk. Ultimately, we hope that the present study will move beyond descriptive associations and lay the groundwork for a comorbidity-aware framework in RCC management, a framework where treatment decisions account not just for tumor biology, but also for the underlying kidney health of patients. Fig. 1 outlines the integrated workflow of the present study, from data integration to therapeutic hypothesis generation.

Figure 1 Study flow chart. CKD, chronic kidney disease; RCC, renal cell carcinoma.

Materials and methods

Data sources and preprocessing

RNA-sequencing (seq) data from The Cancer Genome Atlas (TCGA)-KIRC dataset (https://portal.gdc.cancer.gov/), which included 539 tumor samples, 72 adjacent non-tumor samples and clinical records, were analyzed. For CKD, the GSE98603 dataset (GPL13497 platform; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98603) from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo)(), comprising renal cortical biopsies from 9 patients with CKD and 9 healthy controls, was used. An ICD-related gene list was compiled from GeneCards by querying ‘immunogenic cell death’ and retaining genes with a relevance score above the median score of 19.16. Raw counts from TCGA were converted to fragments per kilobase of transcript per million mapped reads (FPKM) and log2(FPKM + 1) transformed.

Meta-analysis

A systematic review and meta-analysis were conducted to assess whether CKD may be associated with an increased risk of incident RCC in cancer-free adults. PubMed (https://pubmed.ncbi.nlm.nih.gov/), Web of Science (www.webofscience.com), CNKI (https://www.cnki.net/), WanFang (https://www.wanfangdata.com.cn/index.html), VIP (https://www.cqvip.com/) and CBM (https://www.sinomed.ac.cn/zh/index.jsp?type=wx) were searched from inception to 2025, using predefined key words related to CKD and RCC. The initial search yielded 7,743 records; after removing duplicates, 5,218 unique citations underwent title/abstract screening, followed by full-text assessment of 42 potentially eligible articles. In total, 7 studies met the inclusion criteria (Table SI) (17-23), of which 4 reported sufficient binary data for meta-analysis. Study quality was assessed using the Newcastle-Ottawa Scale (NOS) (24) and all were rated high quality (NOS≥7). Given anticipated heterogeneity, a random-effects model was used to calculate the pooled odds ratio (OR). Publication bias was evaluated using funnel plots, in which the log[OR] was plotted against its standard error for each included study. Funnel plot asymmetry was visually inspected for the overall meta-analysis and for the East Asian subgroup analysis. Egger's regression test was also performed to quantitatively assess funnel plot asymmetry, with P<0.05 considered indicative of significant publication bias. All analyses were performed using RevMan software (version 5.4; The Cochrane Collaboration).

Inclusion criteria

The 7 included studies met the following criteria: i) Cohort, case-control or cross-sectional design assessing the association between CKD and incident RCC; ii) adult general population without baseline RCC; iii) CKD defined by estimated glomerular filtration rate (eGFR) <60 ml/min/1.73 m2, albuminuria (albumin-to-creatinine ratio ≥30 mg/g) or International Classification of Diseases diagnostic codes (Fifth Edition, 2016) (25); iv) RCC confirmed pathologically, radiologically or by cancer registry; v) reported ORs, hazard ratios (HRs) or relative risks with 95% confidence intervals (CIs), or provided raw data to calculate them; four studies contributed to the binary-variable meta-analysis, and the other three reported consistent positive associations without being pooled; vi) original peer-reviewed articles in English or Chinese, excluding reviews, case reports or preclinical studies; vii) NOS score ≥7, indicating high quality.

Identification and functional analysis of key CDK-associated genes (CKDGs)

Differentially expressed genes (DEGs) in the TCGA-KIRC cohort were identified by comparing tumor samples (n=539) with adjacent non-tumor tissues (n=72) using the edgeR package (v4.2.1) in R (https://www.r-project.org/). Genes with |log2[fold change (FC)]|>1 and false discovery rate (FDR) <0.05 were considered significantly dysregulated, following the standard practice for TCGA RNA-seq analysis (26). For the CKD GSE98603 dataset, differential expression analysis was performed via GEO2R-2025, applying thresholds of |log2FC|>0.5 and adjusted P<0.05 to account for the smaller effect sizes typically observed in non-neoplastic kidney disease (27) (Table SII).

The intersection of the KIRC DEGs, CKD DEGs and the ICD-related gene set from GeneCards was taken to identify common key genes (28,29), which were designated CKDGs. To interpret the biological roles and pathway involvements of the identified CKDGs, Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment were performed with the clusterProfiler R package (version 4.12.6), adopting a significance threshold of FDR<0.05 (30,31). The prognostic value of each CKDG was assessed using univariate Cox proportional hazards regression on overall survival (OS) data from the TCGA-KIRC cohort, estimating HRs and corresponding P-values (32). To explore functional relationships among the CKDGs, a protein-protein interaction (PPI) network was built using the STRING database (version 12.0; https://cn.string-db.org/), limited to Homo sapiens and retaining interactions with a confidence score ≥0.4. The resulting network was imported into Cytoscape (v3.9.1; https://cytoscape.org/) for visualization and hub gene identification based on degree centrality.

Single-cell RNA-seq data analysis

The single-cell RNA-seq GSE304466 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE304466) dataset from the GEO was analyzed to characterize the tumor immune microenvironment in KIRC. Raw count matrices were processed using Seurat (v5.1.0) in R (33). Cells were retained if they expressed ≥200 genes and had <20% mitochondrial reads; genes detected in <10 cells were excluded. Data were normalized using the LogNormalize method, scaled and subjected to principal component analysis. Clusters were identified via the Louvain algorithm on a k-nearest neighbor graph (resolution=0.6) and cell types were annotated based on canonical markers (34). To evaluate CKDG activity across cell populations, a CKDG signature score was calculated for each cell using Seurat's AddModuleScore function (35). For intercellular communication, cells were stratified into high- and low-CKDG-score groups (median split) and CellChat (v2.2.0) analysis was performed to compare ligand-receptor interaction patterns between groups (36,37). Pseudotime trajectory analysis of T-cell subsets was performed using the Monocle3 R package (version 1.0; https://github.com/cole-trapnell-lab/monocle3) within the R statistical computing environment (version 4.1.0). The analysis followed the standard Monocle3 workflow: Dimensionality reduction using Uniform Manifold Approximation and Projection, cell clustering, graph learning and pseudotime calculation. Naive CD4+ T cells were designated as the root of the trajectory.

Identification of CKD-associated molecular subtypes and immune landscape analysis

Unsupervised consensus clustering of TCGA-KIRC patients based on the expression of the 17 CKDGs was performed using the ConsensusClusterPlus R package (v1.64) (38). After evaluating consensus cumulative distribution functions and the proportion of ambiguous clustering across k=2-6, k=2 was selected as the optimal partition, yielding CKDGs-Cluster 1 and Cluster 2 (Table SIII). Kaplan-Meier survival curves and log-rank tests were performed using the survival R package (version 3.8-3; https://cran.r-project.org/web/packages/survival/index.html) to compare OS between subtypes (Table SIV). Heatmaps were generated using the pheatmap R package (version 1.0.13; https://cran.r-project.org/web/packages/pheatmap/index.html) to display CKDG expression patterns alongside key clinicopathological variables. To identify enriched biological pathways, Gene Set Enrichment Analysis (GSEA; v4.2.3) was performed using hallmark gene sets from MSigDB (https://www.gsea-msigdb.org/), considering results significant at FDR<0.25(39). The tumor immune microenvironment was assessed using two complementary approaches. ESTIMATE scores (stromal, immune and combined) were calculated using the estimate R package (version 1.0.13; https://github.com/oicr-gsi/ESTIMATE) to infer non-tumor content. CIBERSORTx (accessed via https://cibersortx.stanford.edu/) was applied to estimate the relative proportions of 22 immune cell types from bulk RNA-seq profiles (40). Differences in immune composition between clusters were tested using Wilcoxon rank-sum tests. Tumor immune dysfunction and exclusion (TIDE) scores were calculated for each sample using the TIDE web platform (http://tide.dfci.harvard.edu/). Gene expression profiles were submitted to the ‘Response Prediction’ module of the TIDE website, which returns TIDE scores (a higher score indicates greater tumor immune dysfunction/exclusion potential and a lower score predicts immunotherapy response).

Construction of a machine learning-based prognostic signature

A CKDG-based prognostic signature was constructed using TCGA-KIRC patients randomly split into training (70%) and testing (30%) cohorts (set.seed=123). In the training cohort, CKDGs were first screened for association with OS via univariate Cox regression (P<0.05). The significant candidates were then further prioritized using the Random Survival Forest algorithm as implemented in the randomForestSRC R package (version 3.2.0; https://cran.r-project.org/web/packages/randomForestSRC/index.html), retaining genes with positive minimal depth-based variable importance. The remaining genes were subsequently subjected to Least Absolute Shrinkage and Selection Operator (LASSO)-Cox regression (R package glmnet; v4.1-7; https://cran.r-project.org/package=glmnet) with 10-fold cross-validation; the optimal penalty parameter λ was selected at the value yielding the minimum partial likelihood deviance (λ.min) (41). Genes retained by LASSO were included in a multivariate Cox proportional hazards model to derive a risk score:

where βi is the regression coefficient and Expri is the expression of gene i. The median risk score from the training cohort was used to dichotomize patients into high- and low-risk groups in both cohorts. Model performance was evaluated by Kaplan-Meier survival analysis (log-rank test) and time-dependent receiver operating characteristic (ROC) curves (R package timeROC; version 0.4.1; https://cran.r-project.org/web/packages/timeROC/index.html).

Phenome-wide association study (PheWAS)

To explore the clinical relevance of CKDGs (particularly the top-priority gene, WT1) a phenome-wide scan was performed using publicly available PheWAS results from the genome-wide association study (GWAS) catalog (accessed June 2025; https://www.ebi.ac.uk/gwas/) and the UKBiobank TOPMed-imputed PheWeb (https://pheweb.org/UKB-TOPMed). Gene-trait associations derived from gene-based tests (such as MAGMA or VEGAS2) that mapped single nucleotide polymorphisms to WT1 and other CKDGs were focused on. Given the large number of tested phenotypes (>1,000), a Bonferroni-corrected significance threshold (P<0.05/N) was applied and associations at FDR<0.1 were reported. Renal carcinoma and related kidney disease phenotypes were highlighted a priori based on the focus of the present study (42).

MD simulations and molecular docking

The crystal structure of the WT1 zinc finger domain (PDB ID: 1A4Y; residues 270-340; https://www.rcsb.org/) was retrieved from the RCSB PDB (43) Virtual screening was performed using AutoDock Vina software (v1.1.2; http://vina.scripps.edu/) by docking compounds from the zinc small molecule database into the active site of WT1. Top-ranking candidate compounds were selected based on binding free energy (ΔG, kcal/mol) and binding mode analysis. MD simulations were carried out on the most promising protein-ligand complexes using Gromacs (2021) software (version 2021.1; https://www.gromacs.org/) (44,45). Following the MD simulations, complex stability was analyzed by computing the root mean square deviation (RMSD), reflecting atomic positional drift, and the radius of gyration (Rg), indicating overall compactness, across the production trajectory. Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) binding free energy calculations were performed using the gmx_MMPBSA tool (version 1.6.0; https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/) based on the molecular dynamics trajectories generated with Gromacs 2021.

Cell culture

Under standard culture conditions (37˚C, 5% CO2, humidified), UOK276 cells (kindly provided by Professor Zhang, Zhejiang Provincial Key Laboratory of Precision Diagnosis and Therapy for Urological Diseases, Hangzhou, China) were propagated in DMEM (Procell Life Science & Technology Co., Ltd.) supplemented with 10% fetal bovine serum (FBS) (Procell Life Science & Technology Co., Ltd.).

Cell counting kit-8 (CCK-8) assay

A dose-response assay was conducted to determine the IC50 of sirolimus (MedChemExpress) against UOK276 cells using the CCK-8 kit (Dojindo Laboratories, Inc.). Briefly, cells in the log phase were seeded (5x103/well) in 96-well plates. After adherence, the cells were exposed to a gradient of sirolimus concentrations (0, 1.25, 2.5, 5, 10, 20, 40 and 80 µg/ml) in quintuplicate for 48 h at 37˚C in a humidified atmosphere containing 5% CO2. Post-treatment, each well was incubated with 10 µl of CCK-8 solution for 2 h. The absorbance at 450 nm was then recorded. Cell viability (%) was calculated relative to the untreated control and the IC50 value was generated by non-linear curve fitting in GraphPad Prism 9.0 (Dotmatics).

Colony formation assay

A colony formation assay was employed to measure the long-term proliferative capacity of UOK276 cells following sirolimus treatment. Cells were seeded at a low density (500 cells/well) in 6-well plates and, after attachment, exposed to sirolimus at 0, 10, 15 or 20 µg/ml (triplicates per concentration). The culture was continued for 10 days with refreshment of drug-containing medium every third day. Following the appearance of macroscopic colonies in the control wells, the assay was terminated. Subsequent to PBS washing, cells were fixed with 4% paraformaldehyde (25˚C, 15 min) and stained with 0.1% crystal violet (25˚C, 30 min). After washing and drying, colonies (>50 cells) were counted using ImageJ (version 1.53t; National Institutes of Health). Three independent experimental repeats were performed.

Cell invasion assay (Transwell)

To determine whether sirolimus impairs the invasive potential of UOK276 cells, a Transwell invasion assay was performed. Chambers (8 µm pore; Corning, Inc.; cat. no. 3428) were pre-coated with a layer of Matrigel (BD Biosciences; cat. no. 356234) at 37˚C for 30 min. Cells, pre-exposed to sirolimus (0, 10, 15 or 20 µg/ml) for 48 h, were placed in the serum-free upper chamber (1x105 cells/ml in 200 µl), while the lower chamber contained medium with 20% FBS to induce chemotaxis. After 48 h, cells that had not invaded were removed from the upper side of the membrane. The cells that had invaded to the lower side were then fixed with 4% paraformaldehyde (25˚C, 30 min), stained with 0.1% crystal violet (25˚C, 30 min) and quantified by counting five random fields per membrane under an inverted microscope. Three independent experimental replicates were carried out.

Apoptosis detection: Terminal deoxynucleotidyl transferase dUTP nick-end labeling (TUNEL) staining

Sirolimus-induced apoptosis in UOK276 cells was assessed using the TUNEL assay. Cells were treated with sirolimus (0, 10, 15 or 20 µg/ml) for 48 h at 25˚C, then fixed with 4% paraformaldehyde (25˚C, 30 min) and permeabilized with 0.1% Triton X-100 (25˚C, 30 min). DNA fragmentation was detected using the In Situ Cell Death Detection Kit, Fluorescein (Roche Diagnostics; cat. no. 12156792910), according to the manufacturer's protocol. Nuclei were counterstained with DAPI (1 µg/ml at 25˚C for 10 min. Images were acquired using a Zeiss LSM 880 confocal microscope (63x oil objective; excitation 488 nm for FITC, 405 nm for DAPI). TUNEL-positive nuclei (green fluorescence) were quantified in five random fields per sample using ImageJ (version 1.53t). Data represent three independent biological replicates.

Statistical analysis

Statistical analyses were performed using R (v4.4.1) for bioinformatics and GraphPad Prism 9.0 for experimental data. All cell-based experiments were conducted in three independent biological replicates and data are presented as the mean ± SD. Normality was assessed using the Shapiro-Wilk test and homogeneity of variance by Levene's test. For comparisons across more than two groups, one-way ANOVA followed by Tukey's post hoc test was used for pairwise comparisons. P<0.05 was considered to indicate a statistically significant difference.

Results

Confirmation of a significant association between CKD and RCC risk

The present meta-analysis identified a significant association between pre-existing CKD and an increased risk of new-onset RCC, with a pooled OR of 2.55 (95% CI, 1.80-3.59; P<0.001; Fig. 2A). Notable heterogeneity was observed across studies (I²=97%). Subgroup analysis by geographical region revealed an OR of 3.26 (95% CI: 2.96-3.59; I²=0%) in North American studies and 2.11 (95% CI: 1.51-2.95; I²=95%) in East Asian studies (Fig. 2B). The between-group difference was statistically significant (P=0.01). Within the East Asian subgroup, marked methodological differences were noted between Park et al (18) and Park et al (23), including CKD diagnostic criteria (eGFR-only vs. eGFR + albuminuria), control selection (population-based vs. hospital-based) and adjustment for hypertension and diabetes, which may account for the residual heterogeneity. Notably, 3 additional studies (17,19,22) reported qualitatively consistent positive associations between CKD and RCC risk but were excluded from quantitative synthesis due to incompatible data formats. Visual inspection of the funnel plot for the overall analysis (Fig. 2C) revealed a roughly symmetrical distribution of effect estimates around the pooled log[OR], suggesting no substantial publication bias. Similarly, the funnel plot for the East Asian subgroup (Fig. 2D) showed generally symmetrical scatter, indicating that publication bias was unlikely to have distorted the subgroup estimate.

Figure 2 Forest plots of the association between chronic kidney disease and the risk of renal cell carcinoma. (A) Forest plot of the pooled analysis using a random-effects model (Mantel-Haenszel method) for all included studies. The overall OR was 2.55 (95% CI: 1.80-3.59; P<0.00001), with substantial heterogeneity (I²=97%). The size of each square is proportional to the study weight, and horizontal lines represent the 95% CI. (B) Forest plot of the subgroup analysis restricted to studies conducted in East Asia, using the same random-effects model. The pooled OR was 2.11 (95% CI: 1.51-2.95), with high heterogeneity (I2=95%). Funnel plot for publication bias assessment of the overall analysis. The scatter plot shows the log OR against its SE. Symmetry of the distribution was visually inspected. (C) Funnel plot for publication bias assessment of the overall analysis. The scatter plot shows the log OR against its SE. Symmetry of the distribution was visually inspected. (D) Funnel plot for the East Asian subgroup analysis, used to assess potential publication bias. CI, confidence interval; OR, odds ratio; SE, standard error.

Identification and characterization of key genes in RCC and CKD comorbidity

To identify shared molecular features between KIRC and CKD, DEGs from the TCGA-KIRC cohort and the CKD GSE98603 dataset were integrated. Differential expression analysis identified 10,826 DEGs in TCGA-KIRC (|log2FC|>1, FDR <0.05; Fig. 3A) and 15,918 DEGs in GSE98603 (|log2FC|>0.5, adj. P<0.05). Intersecting these DEGs with a curated set of 99 ICD-related genes (GeneCards relevance score ≥10) yielded 17 overlapping genes (Fig. 3B). These co-dysregulated, ICD-associated genes were designated CKDGs. KEGG pathway enrichment of the CKDGs highlighted involvement in oncogenic signaling, including the MAPK and PI3K-Akt pathways, as well as microRNAs in cancer (Fig. 3C). GO analysis showed enrichment in biological processes, such as epithelial tube morphogenesis and positive regulation of transferase activity, and in molecular functions including protein tyrosine kinase activity (Fig. 3D). Univariate Cox regression in the TCGA-KIRC cohort identified 4 CKDGs significantly associated with OS (P<0.05; Fig. 3E): 3 risk genes (HR>1) and 1 protective gene (HR<1). A PPI network constructed via STRING (confidence ≥0.7) identified WT1, adrenomedullin and C-C chemokine receptor 7 as hub genes based on degree centrality (Fig. 3F), prioritizing them for downstream validation.

Figure 3 Identification and characterization of CKDGs shared between KIRC and CKD. (A) Volcano map of differential gene expression in the TCGA-KIRC cohort. (B) Venn diagrams of CKDGs from the TCGA, GEO and GeneCards databases. (C) KEGG analysis and (D) Gene Ontology analysis based on CKDGs. (E) Prognostic analysis based on CKDGs. (F) Protein-protein interaction network based on the CKDGs. TCGA, The Cancer Genome Atlas; FDR, false discovery rate; KIRC, kidney renal clear cell carcinoma; GEO, Gene Expression Omnibus; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; MF, molecular function; CKDG, chronic kidney disease-associated gene; HR, hazard ratio; CI, confidence interval.

Single-cell RNA sequencing (scRNA-seq) reveals an association between the CKD signature score and the immune microenvironment

To investigate the cellular localization and function of CKDGs at single-cell resolution, the GSE304466 KIRC scRNA-seq dataset was analyzed. After rigorous quality control and unsupervised clustering, cells were annotated into nine major types (Fig. 4A-C): Cluster of differentiation 4 (CD4)+ T cells (42.6%), CD8+ T cells (39.2%), natural killer cells (4.9%), B cells (4%), macrophages (3.5%), cancer-associated fibroblasts (2.1%), endothelial cells (0.4%), regulatory T cells (Tregs) (0.2%) and dendritic cells (DCs; 0.1%). Expression profiling of the 17 CKDGs across these cell types showed a distinct pattern, with higher expression levels and proportions particularly in Tregs, CD4+ T cells and CD8+ T cells (Fig. 4E). A CKDG signature score was calculated for each cell using Seurat's AddModuleScore function (default parameters). The score varied across cell types (Fig. 4D), with Tregs showing the highest median score (Fig. 4F). When stratifying T cells by median CKDG score (high vs. low), the high-score group exhibited significantly higher cytotoxicity scores (based on granzyme B, perforin-1 and interferon-γ expression) and higher T-cell exhaustion scores [based on programmed cell death protein 1 (PDCD1), cytotoxic T-lymphocyte-associated protein 4 (CTLA4), lymphocyte-activation gene 3 (LAG3) and T-cell immunoreceptor with Ig and ITIM domains (TIGIT)], compared with the low-score group (P<2x10-16, Wilcoxon rank-sum test; Fig. 4G). CellChat analysis revealed a greater number of significant ligand-receptor interactions in the high CKDG score group (Fig. 5A). Differential communication analysis identified enhanced signaling between CD4+ T cells and Tregs, particularly involving the major histocompatibility complex class I (MHC-I) and platelet-derived growth factor subunit A (PDGFA) pathways (Fig. 5B). Pathway-level comparison showed significantly higher communication probability for the MHC-I, Annexin and PDGFA pathways in the high-score group (P<0.05, permutation test; Fig. 5C). Pseudotime trajectory analysis of T cell subsets, performed using Monocle3 with naive CD4+ T cells as the root, reconstructed their developmental path (Fig. 5D). Stratifying cells into ‘early’ and ‘late’ states based on the median pseudotime revealed a significant distribution difference (P<0.05, χ2 test). Cells in the ‘early’ state predominantly had low CKDG scores, whereas ‘late’-state T cells possessed significantly higher scores (Fig. 5E and F).

Figure 4 Single-cell landscape of CKDG expression in the GSE304466 dataset. (A) Heat maps, (B) cell type maps and (C) cell proportion maps of the single-cell dataset GSE304466 after strict quality control and cell annotation. (D) Projection of CKDG expression scores onto the cell type maps of GSE304466. (E) Bubble plot showing the expression levels of CKDGs across nine distinct cell types. (F) Violin plots quantifying the differences in CKDG expression among cell subpopulations. (G) The cytotoxicity score and the T-cell exhaustion score. CKDG, chronic kidney disease-associated gene.

Figure 5 CKDG-enriched T-cell communication and pseudotime trajectory analysis. (A) Differential interaction network between CKD score subgroups, identified by permutation test (false discovery rate <0.05). (B) Cell communication diagram showing ligand-receptor interactions among T-cell subsets. (C) Ranking chart of the differences in signal pathway activity among each group. (D) Pseudotime trajectory of T-cell differentiation, reconstructed using Monocle3 (Q<0.01). (E) Pseudotime trajectory colored by CKD score. (F) CKD score progression along pseudotime. CKDG, chronic kidney disease-associated gene.

CKDG-based molecular subtyping of bulk RNA data revealed KIRC subtypes with distinct immune features

Unsupervised molecular subtyping of the TCGA-KIRC cohort based on the 17 CKDGs using non-negative matrix factorization identified two stable subtypes at k=2, as supported by consensus clustering (Fig. 6A). The cohort was stratified into CKDG-Cluster 1 (n=546) and CKDG-Cluster 2 (n=698). Kaplan-Meier analysis revealed a significantly improved OS in Cluster 1 (P=0.0089, log-rank; Fig. 6B). A heatmap confirmed higher expression of multiple CKDGs in Cluster 2 (Fig. 6C). GSEA revealed enrichment in ‘Cytokine-cytokine receptor interaction’ and ‘Metabolic pathways’ in Cluster 2 (Fig. 6D). Clinicopathological comparison indicated a significant difference in T-stage distribution (P=0.003), but not in age, N-stage or M-stage (Fig. 6E). Immune deconvolution using CIBERSORTx revealed a higher infiltration of M0 macrophages in Cluster 1, whereas Cluster 2 had elevated CD8+ T cells, M2 macrophages and activated DCs (Fig. 7A). Monocyte chemoattractant protein-counter confirmed higher levels of T cells and cytotoxic lymphocytes in Cluster 2 (P<0.05; Fig. 7B). Consistently, ESTIMATE-derived stromal, immune and composite scores were all significantly higher in Cluster 2 (P<0.05; Fig. 7C). The expression of immune checkpoint genes, including CD40, CD48, CTLA4, LAG3, PDCD1 (PD-1) and TIGIT, was significantly upregulated in Cluster 2 (Fig. 7D). Of note, the TIDE score was also significantly higher in this cluster (P<0.001; Fig. 7E).

Figure 6 CKDG-based molecular subtyping of kidney renal clear cell carcinoma in The Cancer Genome Atlas cohort. (A) Consensus matrix heatmap showing two stable subtypes (k=2) identified by non-negative matrix factorization based on the 17 CKDGs. (B) Kaplan-Meier survival curves comparing overall survival between CKDG-Cluster 1 (n=546) and CKDG-Cluster 2 (n=698); log-rank P=0.0089. (C) Heat map showing the differences in gene expression of 17 chronic kidney disease-associated genes. (D) Gene Set Enrichment Analysis of the two subtypes. (E) Clinical characteristic analysis diagram of the two subtypes. CKDG, chronic kidney disease-associated gene.

Figure 7 Tumor microenvironment features of the chronic kidney disease-associated gene subtypes. (A) CIBERSORTx immune fractions. (B) MCP-counter scores. (C) ESTIMATE scores. (D) Immune checkpoint expression. (E) TIDE scores. Statistical analysis by Wilcoxon rank-sum test; *P<0.05, **P<0.01, ***P<0.001, ****P<0.001. MCP, Monocyte chemoattractant protein.

Construction and validation of a CKD-associated prognostic signature

To develop a clinically applicable risk stratification tool, a prognostic signature was constructed based on the 17 CKDGs. These candidates were further prioritized using the Random Survival Forest algorithm, retaining those with positive minimal depth-based variable importance (Fig. 8A and B). Subsequent LASSO-Cox regression identified an optimal 8-gene model at the penalty parameter (λ) corresponding to the minimum cross-validated partial likelihood deviance (Fig. 8D. A multivariate Cox proportional hazards model was then fitted using these 8 genes to compute a risk score for each patient in the TCGA-KIRC cohort (Fig. 8E).

Figure 8 Construction and validation of the 8-gene chronic kidney disease-associated gene-based prognostic signature in TCGA-KIRC. (A) Feature selection using random forest regression (ntree=500, mtry=15) to identify candidate genes. (B) Variable importance rankings from Random Survival Forest analysis used to prioritize candidate genes. (C) Distribution of risk scores, patient survival status and follow-up time across the cohort. LASSO-Cox regression model selection including (D) L1 penalty parameter (λ) trajectory and (E) cross-validated partial likelihood deviance versus log(λ). (F) Kaplan-Meier overall survival curves stratified by median risk score (log-rank P<0.0001). (G) Time-dependent receiver operating characteristic curves for 250-, 500- and 750-day survival prediction. TCGA, The Cancer Genome Atlas; KIRC, kidney renal clear cell carcinoma; AUC, area under the curve; TPR, true positive rate; FPR, false positive rate.

Patients were dichotomized into high- and low-risk groups using the median risk score. The risk score distribution, aligned with survival status and follow-up time, showed a higher proportion of death events in the high-risk group (Fig. 8C). Kaplan-Meier analysis confirmed a significantly poorer OS in the high-risk group (log-rank P<0.0001; Fig. 8F). Time-dependent ROC analysis yielded area under the curve values of 0.665, 0.602 and 0.563 for 250-, 500- and 750-day survival prediction, respectively (Fig. 8G), indicating a modest short-term but limited long-term discriminative ability.

Confirmation of WT1 as a key biomarker and computational screening for targeted therapeutics

To assess the broader clinical relevance of the 8-gene prognostic signature, a PheWAS was performed using publicly available GWAS summary statistics. Among the 8 CKDGs, only WT1 exhibited a nominal association with renal carcinoma (P=0.039; Fig. 9). No significant associations were observed for the other 7 genes. Given the prior biological plausibility of WT1 in kidney pathophysiology, in silico drug repurposing screening was conducted. Molecular docking of WT1 against a panel of Food and Drug Administration (FDA)-approved compounds identified sirolimus (-7.2 kcal/mol), liposomal daunorubicin (-7.5 kcal/mol) and tretinoin (-6.6 kcal/mol) as top candidates based on predicted binding affinity (Fig. 10A-C). To evaluate the stability of these interactions, 100-ns MD simulations were performed. Analysis of RMSD, Rg, solvent-accessible surface area and hydrogen bonding indicated that all three complexes remained structurally stable over the simulation period (Figs. 11 and 12A-C). Binding free energies estimated by MM/PBSA suggested favorable interactions, with sirolimus exhibiting the lowest (-122.88 kJ/mol) ΔG among the tested compounds (Fig. 12D). In combination, these computational analyses prioritized sirolimus as a candidate for further experimental validation of WT1 modulation.

Figure 9 Phenome-wide association study analysis chart of Wilms tumor 1.

Figure 10 Molecular docking map of WT1 with potential drug targets. (A) WT1 and sirolimus. (B) WT1 and daunorubicin liposomal. (C) WT1 and tretinoin. (D) WT1 and curcumin. (E) WT1 and cytarabine. (F) WT1 and halofuginone. (G) WT1 and deferoxamine. (H) WT1 and dimethyl sulfoxide. WT1, Wilms tumor 1.

Figure 11 Molecular dynamics simulation diagrams. (A) WT1 and sirolimus. (B) WT1 and daunorubicin liposomal. (C) WT1 and tretinoin. (D) WT1 and curcumin. (E) WT1 and cytarabine. WT1, Wilms tumor 1.

Figure 12 Molecular dynamics simulation diagrams. (A) WT1 and halofuginone. (B) WT1 and deferoxamine. (C) WT1 and dimethyl sulfoxide. (D) Energy Components diagram of the molecular dynamics simulation. WT1, Wilms tumor 1.

In vitro validation of the antitumor efficacy of sirolimus. Following the computational identification of sirolimus as a potential modulator of WT1, its antitumor effects in the renal carcinoma cell line UOK276 were evaluated. The CCK-8 assay yielded an IC50 of 15.09 µg/ml (≈16.7 µM) for sirolimus (Fig. 13A). Based on this, concentrations of 10, 15 and 20 µg/ml were selected for functional assays. Sirolimus treatment significantly reduced colony formation in a dose-dependent manner (P<0.05; Fig. 13B and C), indicating suppression of long-term proliferative capacity. Transwell invasion assays showed a concentration-dependent decrease in Matrigel®-penetrating cells (P<0.05; Fig. 13D and F), suggesting impaired invasive potential. Furthermore, TUNEL staining revealed increased apoptosis with higher drug concentrations (P<0.05; Fig. 13E and G). These results demonstrated that sirolimus exerts anti-proliferative, anti-invasive and pro-apoptotic effects in UOK276 cells under the tested conditions.

Figure 13 In vitro antitumor effects of sirolimus in UOK276 renal carcinoma cells. (A) Dose-response curve from the Cell Counting Kit-8 assay used to determine the IC50. (B) Bar graph showing the quantification of colony formation, indicating reduced proliferative capacity after sirolimus treatment. (C) Representative images of colony formation assays demonstrating decreased colony number and size following sirolimus exposure. (D) Representative images of Transwell Matrigel invasion assays showing reduced cell invasiveness after sirolimus treatment. (E) Representative fluorescence images of TUNEL staining indicating increased apoptosis following drug exposure. (F) Bar graph quantifying the Transwell invasion assay results, confirming a significant reduction in invaded cells. (G) Bar graph quantifying the TUNEL-positive cells, demonstrating a statistically significant increase in apoptosis. Statistical analysis by one-way ANOVA with Tukey's post hoc test; *P<0.05, **P<0.01, ***P<0.001. ns, no statistical significance; IC50, half-maximal inhibitory concentration; TUNEL, terminal deoxynucleotidyl transferase dUTP nick end labeling.

Discussion

The present integrative analysis identified 17 CKDGs that were dysregulated in both KIRC and CKD, with enrichment in MAPK and PI3K-Akt signaling, pathways known to drive proliferation and survival in stressed renal tissue. While epidemiological studies have long noted the RCC-CKD comorbidity (46,47), the present study provided a molecular framework linking these conditions through shared transcriptional programs. Of note, single-cell profiling revealed that CKDG expression was enriched in CD4+ T cells and associated with features of functional maturity, challenging the assumption that comorbidity uniformly suppresses antitumor immunity. However, the high-expression CKDG cluster exhibited poor prognosis despite an ‘immune-hot’ phenotype, possibly reflecting concurrent upregulation of immune checkpoints, a pattern that may indicate responsiveness to immunotherapy, as suggested by the elevated TIDE scores.

It is clinically necessary to explore the link between CKD and KIRC/RCC; CKD affects ~11% of adults in the United States (48) and is associated with an elevated risk of kidney cancer, even in the moderate stage, according to a previous meta-analysis (48). Identifying shared pathogenic mechanisms between CKD and RCC is essential for early detection and therapeutic intervention in this high-risk population.

In the present study, among the CKDGs, WT1 showed a nominal genetic association with renal carcinoma in PheWAS (P=0.039), although this did not survive multiple-testing correction. Computational docking and MD simulations predicted stable binding between WT1 and sirolimus and in vitro assays confirmed the antitumor effects of sirolimus in UOK276 cells. However, given the primary action of sirolimus as an mechanistic target of rapamycin (mTOR) inhibitor and the absence of direct evidence linking its efficacy to WT1 modulation, these results should not be interpreted as direct validation of WT1 as a druggable target. Future studies employing WT1 knockdown or CRISPR-based perturbation are needed to dissect its causal role. WT1 encodes a zinc finger DNA-binding protein that acts as a transcriptional activator or repressor depending on cellular context. It is required for normal genitourinary system formation and regulates mesenchymal-to-epithelial transition during kidney development (49). In normal and malignant hematopoiesis, WT1 has been implicated in the regulation of apoptosis, proliferation and differentiation (50).

Despite the integrative multi-omics approach and experimental validation performed in the present study, the present study has several limitations. Of note, the present findings are primarily based on in silico analyses and in vitro cellular models, lacking in vivo confirmation in animal models of RCC and CKD. In addition, the translational potential of the identified WT1-sirolimus axis remains to be evaluated in clinical settings. Future work will focus on two key directions: i) Leveraging WT1 as a molecular target to screen and validate novel therapeutic candidates for RCC through high-throughput drug discovery platforms; and ii) elucidating previously uncharacterized mechanisms of sirolimus action in RCC beyond canonical mTOR inhibition, particularly in the context of CKD-RCC comorbidity. Ultimately, advancing these findings into preclinical in vivo models and early-phase clinical studies will be essential to assess their therapeutic relevance and feasibility.

Beyond cancer, WT1 mutations underlie several congenital disorders. Congenital WT1 mutations cause WT, Wilms tumor, Aniridia, Genitourinary anomalies and Range of developmental delays syndrome, Denys-Drash syndrome and Frasier syndrome (51). Somatic WT1 mutations occur in acute and chronic myeloid leukemia, myelodysplastic syndrome and other blood neoplasms, while increased WT1 expression without mutation is observed in leukemia and solid tumors (52). These findings underscore the pleiotropic roles of WT1 across renal, hematological and developmental diseases. WT1 functions both as a tumor suppressor gene and as an oncogene, with this duality arising from isoform diversity and cellular context (53).

In conclusion, in the present study, a multi-omics framework that links CKD and KIRC was presented through a shared gene signature, offering new insights into their comorbid pathogenesis. The identification of WT1 as a candidate hub gene, coupled with preliminary evidence of sirolimus activity, has generated testable hypotheses for mechanistic and translational research. While the present findings do not yet support clinical application, they provide a rationale for exploring WT1-related pathways in the context of RCC-CKD comorbidity and highlight the need for models that recapitulate both tumor and renal dysfunction.

Finally, it is informative to put the in vitro observation of sirolimus activity in UOK276 cells into clinical context. Sirolimus was first approved by the FDA in 1999 for the prophylaxis of organ rejection in patients aged ≥13 years receiving renal transplants. It is a mTOR inhibitor with both immunosuppressive and antitumor activities (54). Sirolimus has also been used in drug-eluting stents to prevent restenosis following angioplasty. Regarding its safety profile, sirolimus is associated with increased susceptibility to infection, and possible development of lymphoma and other malignancies due to immunosuppression. The use of sirolimus in combination with tacrolimus has been associated with excess mortality and graft loss in patients following de novo liver transplant (55). However, conversion to sirolimus in renal transplant recipients with a history of cancer appears safe regarding renal function and graft survival, with patient survival largely dependent on tumor entity (56). These established clinical safety and efficacy profiles support the exploration of sirolimus as a repurposed agent in selected cancer populations.

Supplementary Material

Overall included data for the meta-analysis.

CKD-related gene data obtained from the GEO database.

FPKM file of the TCGA-KIRC queue.

Survival data file of the TCGA-KIRC queue.

Acknowledgements

Not applicable.

Funding

Funding: The present study was supported by Special Fund Project for Clinical Medical Research of Zhejiang Medical Association (grant no. 2025ZYC-A13), Science and Technology Co construction Project of China Traditional Chinese Medicine Comprehensive Reform Demonstration Zone (grant no. GZY-KJS-ZJ-2026-115) and Joint TCM Science & Technology Projects of National Demonstration Zones for Comprehensive TCM Reform (grant no. GZY-KJS-ZJ-115).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

WZ, PT, SW and ZS conceived and designed the experiments and analysis. WZ and ZS performed the experiments and statistical analysis. WZ, PT, SW and ZS provided technical support and wrote the paper. WZ and ZS provided funding support and supervision. All authors read and approved the final version of the manuscript. ZS and WZ confirm the authenticity of all the raw data.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References