Fibrogenesis-driven tumor progression in clear cell renal cell carcinoma: prognostic, therapeutic implications and the dual role of neuropilin-1

Abstract

Background

Clear cell renal cell carcinoma (ccRCC) is the predominant subtype of renal cancer, with a poor prognosis driven by therapy resistance and a propensity for recurrence. Tumor microenvironment (TME)-associated fibrosis accelerates disease progression by fostering immune evasion. Neuropilin-1 (NRP1), a key mediator in fibrotic signaling and cancer biology, has been implicated in these processes. However, the genetic correlation between fibrogenesis and ccRCC remains largely unexplored, necessitating a focused analysis of fibrogenesis-related genes (FRGs) to identify novel prognostic markers and therapeutic strategies.

Methods

This study utilized an integrative bioinformatics framework to identify prognosis-associated fibrogenesis-related genes (pFRGs) and applied non-negative matrix factorization (NMF) to stratify ccRCC patients based on fibrotic signatures. A machine learning-derived prognostic model was developed to categorize patients into high-risk and low-risk groups, with tumor microenvironment (TME) features analyzed across these subgroups. The pro-tumorigenic role of NRP1 via the TGF-β/SMAD signaling pathway was validated in vitro and in vivo.

Results

Twelve pFRGs were identified, with elevated expression correlating with reduced survival. NMF revealed two ccRCC subtypes with different fibrotic and immune profiles. The high-fibrosis subtype showed worse survival and a pro-tumorigenic TME. The risk model demonstrated robust predictive performance (AUCs: 0.738, 0.731, 0.711 for 1-, 2-, and 3-year survival). High-risk patients, marked by immune dysfunction, exhibited worse survival but greater immunotherapy sensitivity. Among the pFRGs, NRP1 was upregulated in ccRCC, and paradoxically associated with favorable prognosis in TCGA, primarily due to stromal enrichment. In vitro and in vivo experiments confirmed that NRP1 promotes ccRCC proliferation, migration, and invasion by enhancing TGF-β/SMAD-driven epithelial-mesenchymal transition (EMT).

Conclusion

Fibrosis is a critical driver of ccRCC progression, linking fibrogenesis-related genes to poor prognosis, immune suppression, and tumor aggressiveness. NRP1 was identified as a central regulator of fibrosis-induced tumor progression through the TGF-β/SMAD signaling pathway. Combining NRP1 inhibition with anti-fibrotic therapies presents a potential strategy for enhancing therapeutic outcomes in ccRCC.

Graphical abstract

Introduction

Clear cell renal cell carcinoma (ccRCC) is the most common subtype of kidney cancer, with surgical resection as the primary treatment due to its resistance to radiotherapy and chemotherapy. However, approximately 30% of patients experience recurrence or metastasis even after early surgical intervention. For advanced ccRCC, the 5-year survival rate is alarmingly low at just 23% [1, 2]. This underscores the urgent need to elucidate the molecular mechanisms underlying ccRCC’s aggressive proliferation and invasion, with the goal of identifying early diagnostic and prognostic markers and developing effective therapeutic targets.

Fibrosis, typically a response to chronic injury, involves the excessive buildup of extracellular matrix (ECM), driven by fibroblast and immune cell activation [3]. This response is similar to the tumor microenvironment in cancer, where fibroblasts, inflammatory cells, and ECM changes foster an environment that supports cancer progression [4]. Fibrotic responses are common across various malignancies, such as lung, pancreatic and breast cancers [5,6,7,8]. In the tumor microenvironment (TME), fibrogenesis involves the activation of cancer-associated fibroblasts (CAFs), excessive extracellular matrix (ECM) deposition, and intratumoral fibrosis (ITF) [4, 9,10,11,12]. These elements have proven pivotal in driving tumorigenesis, evading the immune response, promoting metastasis, and fostering resistance to drugs across various solid tumors [4, 13]. In ccRCC, fibrosis and cancer share key pathways, including TGF-β signaling, epithelial-to-mesenchymal transition (EMT), and increased ECM stiffness. These factors collectively foster an immunosuppressive and tumor-promoting microenvironment [14]. Innovative therapeutic strategies that target these shared pathways, with a special interest in agents that may simultaneously mitigate fibrosis and cancer progression, offering a promising approach for treating ccRCC [15].

The advent of multi-omics has positioned transcriptome profiling as a critical tool in ccRCC research, enabling biomarker discovery, tumor heterogeneity analysis, and investigations into metastasis and therapy resistance. Additionally, exploring fibrosis in ccRCC progression has highlighted new therapeutic opportunities through in-depth mining of patient-derived transcriptomic datasets. Recent single-cell RNA sequencing (scRNA-seq) analysis of human RCC revealed that pro-fibrotic signaling pathways, such as TGF-β, WNT, and mTOR, play critical roles in driving tumorigenesis and enhancing tumor aggressiveness [16]. Despite this, the relationship between fibrogenesis and ccRCC lacks genetic-level research, and the potential for cancer-associated fibrosis to serve as a prognostic indicator in ccRCC warrants exploration.

This study leveraged scRNA-seq and spatial transcriptomics data from TCGA and GEO databases to uncover prognostic fibrogenesis-related genes (pFRGs). By applying non-negative matrix factorization (NMF), ccRCC samples were divided into two subtypes characterized by distinct fibrotic signatures, clinical features, and immune profiles. A machine learning-derived prognostic model demonstrated reliable predictive accuracy (AUC: 0.738 for 1-year, 0.731 for 2-year, and 0.711 for 3-year survival) and was validated using an independent dataset, effectively stratifying patients into risk categories. This stratification supports personalized prognoses and facilitates the development of immunotherapy strategies. Additionally, in vitro and in vivo experiments confirmed NRP1’s role in enhancing ccRCC progression. Notably, our findings provide the first evidence that NRP1 promotes tumor aggressiveness by regulating EMT via the TGF-β/SMAD pathway, thereby linking NRP1 activity to ccRCC advancement.

Materials and methods

Dataset acquisition and processing

We accessed transcriptional data for 537 ccRCC samples from the TCGA-KIRC cohort and 100 normal samples from the Genotype-Tissue Expression (GTEx) database through the UCSC Xena platform (https://xenabrowser.net/datapages/). A total of 78 fibrogenesis-related genes (FRGs) were curated from the Harmonizome database (Table S1) [17]. To enhance our analysis, we retrieved additional datasets from the Gene Expression Omnibus (GEO) and ArrayExpress database, including GSE29609, GSE53757, GSE126964 and E-MTAB-3267 for validation purposes.

Identification of prognosis-related FRGs

Differential expression analysis was conducted using the “limma-voom” method, applying an adjusted p-value threshold of < 0.05 and |logFC|> 1 [18]. Principal component analysis (PCA) was performed on the identified differentially expressed genes (DEGs). By intersecting DEGs with the FRGs, we obtained a set of DEFRGs between tumor specimens and control specimens. Prognosis-associated DEFRGs (pFRGs) were identified through univariate Cox regression analysis using the “coxph” function, with p < 0.05 as the significance threshold. Kaplan–Meier (K–M) survival analysis for pFRGs was conducted through the GEPIA2 database, applying 75% and 25% as cutoff thresholds [19]. The"Rcircos"package was employed to visualize the chromosomal location of pFRGs [20]. Gene alterations were visualized using the"Maftools"package [21].

Transcriptome analysis at the single-cell level

Single-cell RNA sequencing (scRNA-seq) data were obtained from two sources. The GSE171306 dataset was accessed via the Tumor Immune Single-cell Hub 2 (TISCH2) platform and analyzed through the platform’s integrated MAESTRO pipeline, which includes quality control, batch effect correction, cell clustering, and automatic cell type annotation [22]. Additionally, the GSE159115 dataset was downloaded from the GEO database and processed locally using the Seurat (v5.0.2) R package. Standard preprocessing steps were applied, including filtering, normalization, scaling, dimensionality reduction, clustering, and annotation based on canonical marker genes. The BayesPrism R package was used to perform deconvolution by integrating cell-type proportion information from scRNA-seq into TCGA bulk RNA-seq data [23].

Validation of spatial transcriptomics and protein expression

We utilized the STOmicsDB portal and SpatialTME platform for ccRCC spatial transcriptomic data acquisition, analysis, and visualization of the expressions of pFRGs in clinical samples [24, 25]. Additionally, we consulted the online Human Protein Atlas (HPA) database for protein expression data [26].

Immune infiltration analysis

Immune infiltration of pFRGs was evaluated using the single-sample gene set enrichment analysis (ssGSEA) algorithm. STRING database was utilized to generate Protein–protein interaction (PPI) networks [27].

Non-negative matrix factorization (NMF) algorithm

Using the “NMF” R package, we performed NMF clustering analysis based on pFRG expression to identify potential ccRCC molecular subtypes [28]. The “brunet” criterion was applied with 500 iterations, exploring cluster numbers (k) from two to eight, requiring a minimum of ten samples per cluster. Optimal k was determined through cophenetic, dispersion, and silhouette metrics, resulting in two distinct molecular clusters. To assess clustering robustness, we conducted tenfold cross-validation based on reconstruction error and bootstrap resampling (n = 10) with pairwise Adjusted Rand Index (ARI) calculation. Immune cell abundance and responses within each cluster were estimated using ssGSEA. Clinical feature differences between clusters were analyzed via the Chi-square test and visualized as stacked bar plots.

Functional enrichment analysis

Gene Set Variation Analysis (GSVA) was performed with the “GSVA” R package to investigate the biological characteristics of each cluster [29]. To explore fibrosis-associated biological differences between molecular subtypes, we also curated 57 fibrosis-related gene sets from the Molecular Signatures Database (MsigDB). Gene sets were retrieved by performing keyword searches (“fibrosis OR fibrogenesis”) using the MsigDB portal’s “Search” function. Retrieved gene sets were then manually reviewed, and those directly associated with extracellular matrix remodeling, fibroblast activation, TGF-β signaling, and other fibrosis-relevant processes were retained for GSVA Analysis. We also utilized HALLMARK Predefined gene datasets from the MSigDB database [30]. Enrichment scores were analyzed via the"limma"package for differential analysis [31].

Development of the prognostic stratified model

To construct the prognostic stratified model, patients were randomly divided into a training set and a testing set in a 7:3 ratio. The training set served to build the model, while the testing set was used for validation. Using the “glmnet” R package [32], LASSO regression with 1000-fold cross-validation was applied to select key FRGs. These selected FRGs were further analyzed through multivariable Cox regression to identify the optimal PR-FRGs. The model was defined as a linear combination of expression levels and Cox regression coefficients for each optimal PR-FRG, calculated as:

$$Risk\, score = \,\sum\limits_{i = 1}^{n} {\left( {coef_{i} \times Exp_{i} } \right)}$$

Coef represents the Cox hazard ratio coefficient and Exp represents the expression level of the corresponding PR-FRG.

Assessment of the prognostic stratification model

CcRCC patients were stratified into high- and low-risk groups using the median risk score as the cutoff. Kaplan–Meier (K-M) survival analysis was performed to compare survival outcomes between the two groups. The model’s predictive accuracy was evaluated using time-dependent ROC curves generated with the “timeROC” package [33]. To visualize survival trends, risk curves and heatmaps were constructed for both groups. A nomogram combining clinical characteristics was developed with the “rms” R package, and its performance was validated using calibration plots.

Comprehensive TME analysis

Using the IOBR package [34], we employed various immune infiltration algorithms to assess immune cell composition in tumor samples. The ssGSEA algorithm quantified immune cell infiltration and analyzed immune-related pathway activity. To evaluate immune escape potential and predict immunotherapy response, the tumor immune dysfunction and rejection (TIDE) score was calculated. Additionally, immunophenoscore (IPS) data from the Cancer Immunome Atlas (https://tcia.at/home) were utilized to identify ccRCC patients likely to benefit from immunotherapy.

TMB and drug sensitivity anaysis

Tumor mutational burden (TMB) for each patient was calculated using TCGA somatic mutation data with the R package “maftools.” Drug sensitivity to chemotherapeutic agents was assessed by estimating IC50 values through the “oncoPredict” package [35] and the Genomics of Drug Sensitivity in Cancer database. Ridge regression models were employed to predict individual drug responses.

Clinical samples

A total of 20 tissue sections, including ccRCC tumors and adjacent non-cancerous tissues, were collected from surgical patients at Linyi People’s Hospital. These samples were formalin-fixed and paraffin-embedded to analyze NRP1 expression. Pathological diagnoses of ccRCC specimens were independently confirmed by two pathologists. Ethical approval was granted by the Science and Technology Ethics Committee of Linyi People’s Hospital (approval number: 202409-H-027), and all samples complied with defined standards for patient data and specimen quality. Informed consent was obtained from all participants.

Cell culture and NRP1 knockdown

Human ccRCC cell lines (ACHN, RRID:CVCL_1067; 786-O, RRID:CVCL_1051; 769-P, RRID:CVCL_1050; and Caki-1, RRID:CVCL_0234) and the human renal proximal tubular epithelial cell line HK2 were sourced from FuHeng Biology. Cells were cultured in the following media: Minimum Essential Medium (MEM, Gibco) for ACHN, RPMI-1640 (Invitrogen-Gibco) for 786-O and 769-P, McCoy’s 5 A (Gibco) for Caki-1, and DMEM (Gibco) for HK-2. All media were supplemented with 10% fetal calf serum (FCS) and 100 U/mL penicillin/streptomycin. Cell cultures were maintained at 37 °C in a humidified atmosphere with 5% CO₂.

NRP1 knockdown was performed using small interfering RNAs (siRNAs) targeting NRP1 (sense: 5’-CCGACAGCGCGAUAGCAAATT-3’; antisense: 5’-UUUGCUAUCGCGCUGUCGGTT-3’) from RiboBio Co., Ltd. (Guangzhou) and transfected with Lipofectamine 2000 (Thermo Fisher Scientific). Transfected cells were incubated for 48 h prior to subsequent experiments.

Cell proliferation and apoptosis assays

786-O and Caki-1 cells were trypsinized, adjusted to 5 × 10⁴ cells/mL, and seeded at 100 µL per well in 96-well plates. The following morning, cells were transfected according to experimental groups. The medium was refreshed with complete medium 5 h after transfection, and cells were further cultured. t 24, 48, and 72 h, 10 µL of CCK-8 solution was added to each well (final volume: 100 µL). Absorbance at 450 nm was measured after 4 h of incubation using a microplate reader.

For EdU assays, the BeyoClick™ EdU Cell Proliferation Kit (Alexa Fluor 488; Beyotime, Shanghai, China) was used to label proliferating cell nuclei with green fluorescence.

Cell cycle analysis involved treating, trypsinizing, and washing cells with PBS, followed by fixation in 70% ethanol at 4 °C for 4 h or overnight. Cells were resuspended in staining solution, incubated in the dark at 37 °C for 20 min, and analyzed by flow cytometry (BD Accuri C6; BD Biosciences, USA) using 488 nm excitation. Cells were resuspended in Binding Buffer, stained with 5 μL Annexin V-FITC, and incubated in the dark for 15 min to detect apoptosis. PI was added immediately before flow cytometry to assess apoptotic stages.

Western blot analysis

RIPA buffer with protease and phosphatase inhibitors was used to lyse the cells, and protein concentrations were measured using a BCA assay. Equal protein amounts (20–30 μg per lane) were loaded onto SDS-PAGE gels for separation and then transferred onto PVDF membranes. To block non-specific binding, membranes were treated with 5% non-fat milk in TBST for 1 h at room temperature. The expression of key proteins involved in apoptosis, cell cycle regulation, EMT, and the TGF-β/SMAD pathway in 786-O and Caki-1 cells was analyzed using primary antibodies. These included NRP1 (Proteintech Cat# 60067-1-Ig, RRID:AB_2150840, 1:1000), BAX (Abcam Cat# ab32503, RRID:AB_725631, 1:1000), BCL2 (Abcam Cat# ab182858, RRID:AB_2715467, 1:2000), cleaved-caspase 3 (Abcam Cat# ab32042, RRID:AB_725947, 1:500), CDK4 (Proteintech Cat# 11026-1-AP, RRID:AB_2078702, 1:1000), PCNA (Proteintech Cat# 10205-2-AP, RRID:AB_2160330, 1:5000), E-cadherin (Proteintech Cat# 20874-1-AP, RRID:AB_10697811, 1:2000), N-cadherin (Proteintech Cat# 22018-1-AP, RRID:AB_2813891, 1:2000), Vimentin (Abcam Cat# ab92547, RRID:AB_10562134, 1:1000), SMAD2 (Abcam Cat# ab228765, 1:500), SMAD3 (Abcam Cat# ab40854, RRID:AB_777979, 1:1000), p-SMAD2 (Abcam Cat# ab188334, RRID:AB_2732791, 1:1000), and p-SMAD3 (Cell Signaling Technology Cat# 9520, RRID:AB_2193207, 1:1000). GAPDH (Proteintech Cat# 60004-1-Ig, RRID:AB_2107436, 1:10000) served as the loading control. Primary antibody incubation was performed overnight at 4 °C, followed by TBST washing and a 1-h incubation with HRP-conjugated secondary antibodies at room temperature. Enhanced chemiluminescence (ECL) was used for protein detection. Bands were visualized and quantified using ImageJ, with expression levels normalized to GAPDH.

Enzyme-linked immunosorbent assay (ELISA)

Culture supernatants from 786-O and Caki-1 cells were collected after transfection with siRNA (si-NC, siNRP1, or siNRP1 + TGF-β) for 48 h and subjected to ELISA to quantify the concentration of total TGF-β. The TGF-β1 ELISA Kit (Multisciences Cat# EK981-48) was used to measure TGF-β levels.

Wound healing and migration and invasion experiments

Wound healing: 786-O and Caki-1 cells transfected with si-NC or si-NRP1 were seeded in 6-well plates. After 48 h, a 10 μL Eppendorf Tip was used to create a scratch in the cell monolayer. The medium was replaced with serum-free medium, and images were captured at 0 and 48 h using an inverted microscope (IX73, Olympus, Germany).

Transwell experiment: Migration assays were conducted by seeding 786-O and Caki-1 cells, adjusted to 2 × 10⁵ cells/mL, into the upper chambers of Transwell inserts (Costar 3422, Corning, USA) at 10,000 cells per well in serum-free medium. Complete medium containing 10% FBS was added to the lower chambers as a chemoattractant. After 48 h, cells that migrated to the underside of the membrane were fixed with paraformaldehyde, stained with crystal violet, and counted under a microscope. For invasion assays, Matrigel-coated Transwell inserts were used, following the same procedure as described for the migration assay.

Xenograft model

5- to 6-week-old male BALB/c (RRID:IMSR_APB:4790) nude mice were purchased from SPF Biotechnology Co., Ltd (Beijing). All animal experiments were approved by the Animal Experiments and Experimental Animal Welfare Committee of Capital Medical University and conducted according to animal welfare guidelines (approval number: AEEI-2024–288). Mice were subcutaneously injected with 5 × 10⁶ NRP1-knockdown or control ccRCC cells suspended in 100 µL PBS in the right hind flank. Tumor volumes were measured twice weekly with calipers and calculated using the formula V = (L × W²)/2, where L represents tumor length and W is the width.

Immunohistochemistry and TUNEL staining

At the endpoint, tumor tissues were harvested, fixed, dehydrated, and embedded in paraffin. Hematoxylin and eosin (H&E) staining was utilized to evaluate general tissue morphology. For immunohistochemistry, sections were incubated with primary antibodies against Ki67 (Abcam Cat# ab16667, RRID:AB_302459, 1:500), CD31 (Abcam Cat# ab182981, RRID:AB_2920881, 1:500), NRP1 (Abcam Cat# ab81321, RRID:AB_1640739, 1:200), and Vimentin (Abways Cat# CY5134, 1:200). Ki67 was used to measure cell proliferation, CD31 to assess angiogenesis and blood vessel density, NRP1 to evaluate its expression in tumor tissues, and Vimentin to identify mesenchymal marker expression indicative of EMT. Tumor apoptosis was detected via TUNEL staining following the manufacturer’s protocol (P0017, PINUOFEI, Wuhan, China).

TGF-β/SMAD pathway analysis

The activity of the TGF-β/SMAD pathway in NRP1-knockdown cells was evaluated by analyzing SMAD2/3 phosphorylation levels via Western blot. To investigate potential rescue of TGF-β/SMAD pathway activity, exogenous TGF-β (10 ng/mL) was added to NRP1-knockdown cells. The effects on cell proliferation, migration, and invasion were further assessed, confirming that NRP1 promotes tumor progression through the TGF-β/SMAD pathway.

Statistical analysis

R software (v4.3.1) and GraphPad Prism (v9.5, GraphPad Software, USA) served as the primary tools for statistical analyses and data visualization. Continuous variables were displayed as mean ± SD. For two-group comparisons, normally distributed data were analyzed using a two-tailed t-test, while non-normally distributed data were evaluated with the Mann–Whitney U test. One-way ANOVA was employed for analyzing differences among multiple groups. Proportions or rates were used to express categorical data, with the chi-square test or its continuity correction applied for comparisons between high- and low-NRP1 expression groups. Relationships between variables were evaluated using Spearman’s correlation test. Kaplan–Meier curves illustrated survival outcomes, and independent prognostic factors were identified through univariate and multivariate Cox proportional hazard models. A P < 0.05 threshold defined statistical significance, reported as follows: ns (not significant), * (P < 0.05), ** (P < 0.01), *** (P < 0.001).

Results

Identification of pFRGs

In this study, an analysis was conducted on cohorts comprising 537 cases of clear cell renal cell carcinoma and 100 control cases (Table S2). Utilizing a dataset encompassing 78 fibrosis-associated genes, a total of 28 DEFRGs were identified (Figure S1 A-D, Table S3). Subsequently, a univariate Cox analysis was employed to further discern the pFRGs, resulting in the identification of 12 pFRGs (ADAM12, AGER, COL1 A1, ENTPD1, GREM1, IFNG, KCNN4, LRAT, NRP1, TGFB1, TGFBI, and VTN) (Fig. 1A). The prognosis analysis revealed poorer overall survival (OS) and disease-free survival (DFS) among patients with elevated pFRG signature expression, as illustrated by K–M curves (Fig. 1B, C). Moreover, an assessment of the individual impact of each pFRG on prognosis indicated that, apart from ENTPD1 and NRP1, elevated expression of the remaining 10 genes seemed to imply a worsened prognosis (Figure S2).

Fig. 1

Identification of pFRGs. A Twelve pFRGs selected by univariate Cox PH Model. B, C Survival analysis for OS (B) and DFS (C) in high (upper quartile) and low (lower quartile) prognostic signature expression groups. D CNV Frequencies among 12 pFRGs. E PCA plot of the prognostic signature expression profile between ccRCC and normal tissues. F Expression levels of 12 pFRGs in tumor and normal samples. G Heatmap displaying correlations of pFRG expression. H PPI network for pFRGs. r, correlation coefficient. CNV, copy number variation; Dim, dimension; HR, hazards ratio; TPM, transcripts per million

Genomic and transcriptomic characteristics of pFRGs

Chromosomal circular plots were employed to denote the chromosomal positions of pFRGs, predominantly located on chromosomes 10, 17, and 19 (Figure S1E). In the analysis of somatic copy-number variations (SCNVs), the gene TGFBI displayed the highest frequency of amplifications, whereas LRAT exhibited the highest frequency of deletions (Fig. 1D). A mere 2.43% of samples (nine out of 370 samples) exhibited mutations within 6 pFRGs (Figure S1 F). Through expression profiling analysis, a principal component analysis (PCA) plot revealed 12 pFRGs to discriminate between tumor and non-tumor tissues (Fig. 1E). Remarkably differential expression was observed across all pFRGs (Fig. 1F), and further validation of pFRG expression was undertaken using 2 external validation sets (Figure S1 G-H). Furthermore, we provide a comprehensive overview of interactions between 12 pFRGs in a heatmap (Fig. 1G). Most pFRGs demonstrate positive correlations with each other, with COL1 A1 exhibiting the strongest correlation with ADAM12 (r = 0.76). It is noteworthy that AGER displayed a negative correlation with a significant number of pFRGs, suggesting it may have an opposing or distinct role within the fibrotic or immune landscape of ccRCC. Additionally, a PPI network was constructed and visualized using the STRING database to illustrate potential mechanistic crosstalk among the 12 pFRGs (Fig. 1H).

Expression profile of pFRGs in ccRCC across cellular, spatial and protein levels

Using the TISCH 2.0 databse, a total of 11,427 cells from GSE171306 scRNA-seq dataset were clustered into 25 groups and annotated into 10 distinct cell types (Figure S3 A). Marker genes between each cluster were calculated and represented in Figure S3B. 6 pFRGs (COL1 A1, ENTPD1, IFNG, NRP1, TGFB1, and TGFBI) were selected based on their appreciable and heterogeneous expression across cell clusters at the single-cell level. Notably, 5 of them (excluding IFNG) exhibited predominant expression in fibroblasts (Fig. 2A). Besides, we conducted a spatial transcriptomic analysis utilizing the STOMICS database to reinforce our understanding of the spatial expression patterns of pFRGs within the ccRCC TME (Figure S3 C).

Fig. 2

Single cell data analysis and protein level of pFRGs. A Stacked violin plots illustrating expression levels of pFRGs in different cell clusters. B UMAP projections of the expression patterns of 6 highly expressed pFRGs between cell types. C Immunohistochemistry images showing the expression of optimal pFRGs in tumor and normal specimens. The protein expression levels are classified as high, low, and ND. Macro macrophage, Mono monocyte, NK natural killer cell, Tconv conventional T cell, Tprolif proliferative T cell, UMAP uniform manifold approximation and projection, ND not detected

Using the HPA database, we assessed the protein levels of six selected pFRGs (Fig. 3B). Tumor tissues demonstrated significantly increased expression of ENTPD1, TGFBI, and TGFB1 relative to normal tissues, as illustrated in the figure. TGFBI and TGFB1, both central mediators of fibrogenesis, likely promote extracellular matrix remodeling and immune suppression within the TME. These findings corroborate transcriptomic data, further substantiating the importance of these pFRGs in ccRCC progression and fibrotic TME regulation.

Fig. 3

NMF clustering of ccRCC patients. A Consensus matrix generated from NMF clustering with k = 2. B The UMAP plot demonstrating a difference in gene expression between the two clusters in the entire genome. C Differential expression of 12 pFRGs between the two clusters. D Immune infiltration patterns of two clusters. E Overall survival differences between the two clusters. F Stacked bar graph displaying differences of clinicopathological characteristics between two ccRCC clusters. BCR B-cell receptor, HR hazards ratio, MDSC myeloid-derived suppressor cell, TCR T-cell receptor, TPM transcripts per million

Immune infiltration and interactions associated with pFRGs

To gain insights into the immune infiltration status of pFRGs, we conducted a ssGSEA analysis to assess their expression within immune cells and immune pathways (Figure S3D). The heatmap reveals that, whereas most of the remaining genes exhibit a positive correlation, VTN, LRAT, and AGER negatively correlate with immune responses, suggesting these genes may contribute to immune evasion or suppression within the ccRCC microenvironment. This network sheds light on potential protein interactions among the pFRGs, further elucidating their functional relationships within the context of ccRCC and the immune microenvironment.

Stratification of ccRCC patients based on pFRGs

To further investigate the relationship between fibrosis and ccRCC, we conducted NMF analysis using the expression profiles of the 12 pFRGs in the ccRCC cohort. The cophenetic, dispersion, and silhouette metrics identified k = 2 as the optimal cluster number (Fig. 3A, Figure S4). All ccRCC samples were subsequently divided into two distinct clusters, designated as the C1 cluster and the C2 cluster. The UMAP plot demonstrates significant differences in whole gene expression profiles between the two clusters (Fig. 3B). To evaluate the stability of the NMF clustering, we performed bootstrap resampling (n = 100) and k-fold cross-validation, which yielded a high average Adjusted Rand Index (ARI = 0.838) and a consistently narrow reconstruction error range (1.61–1.99), supporting the robustness of the factorization results. In Fig. 3C, we present the expression landscape of the 12 pFRGs between these two clusters. Notably, except for AGER, which was upregulated in the C2 cluster, all remaining differentially expressed pFRGs showed upregulation in the C1 cluster. Subsequent ssGSEA analysis reveals that the C1 cluster exhibits an overall higher level of immune infiltration and immune responses (Fig. 3D). However, this immune activity likely reflects a dysfunctional or exhausted immune response, consistent with the immunosuppressive effects of fibrosis-associated pathways enriched in Cluster C1. In terms of clinical characteristics, the C1 cluster displayed significantly higher histologic grade, pathological stage, and TMN stage compared to the C2 cluster (Table 1, Fig. 3F), aligning with their worse overall survival as shown by Kaplan–Meier analysis (HR = 1.86, P < 0.001, Fig. 5E).

Table 1 The patient characteristics of the two clear cell renal cell carcinoma clusters

Fibrotic and functional landscapes of ccRCC clusters

To investigate the connection between the two subtypes and fibrosis further, we manually curated 57 fibrosis-related datasets from MsigDB for GSVA analysis (Table S4). Significant differences between the two subtypes were observed across nearly all fibrosis-related pathways, such as fibroblast growth factor signaling, EMT, and extracellular matrix organization. The bar graphs on the right side quantify these differences with Cluster C1 showing higher enrichment scores compared to Cluster C2, demonstrating that Cluster C1 has a more severe degree of fibrosis (Fig. 4A). Moreover, the hallmark gene set (h.all.v2023.1.Hs.symbols.gmt) served as the foundation for GSVA analysis., indicating distinct biological behaviors and pathway activities associated with each subtype (Fig. 4B). As illustrated in the box plots, pathways linked to inflammation, EMT, and angiogenesis were predominantly enriched in Cluster C1, emphasizing its greater tumor progression risk. Cluster C2, however, showed enhanced activity in metabolic pathways, underscoring distinct biological functions and clinical implications.

Fig. 4

GSVA pathway enrichment analysis of the two ccRCC clusters. A Enriched fibrosis-related pathways of two clusters based on GSEA database. B Box plot displaying HALLMARK pathways that are upregulated, downregulated or not correlated between two clusters

To validate the stability and generalizability of the NMF-based clustering approach, we further applied the same NMF procedure to an independent ccRCC dataset (E-MTAB-3267), which again identified k = 2 as the optimal cluster number. The resulting subgroups exhibited similar clinical and fibrotic features to those in the original cohort, supporting the robustness and transferability of our model (Figure S5).

In summary, the C1 cluster is characterized by higher expression levels of pFRGs, significant enrichment in fibrosis-related and inflammatory pathways, and poorer clinical outcomes, suggesting a more aggressive tumor phenotype driven by fibrosis and immune dysregulation. Conversely, Cluster C2 displays a metabolic-driven phenotype with less fibrotic activity, potentially indicating a less aggressive tumor microenvironment. These findings underscore the potential utility of targeting fibrotic and inflammatory pathways in Cluster C1 patients to improve prognosis and guide personalized treatment strategies.

Development and assessment of the prognostic stratified model

Ten-fold cross-validated LASSO regression identified 6 optimal pFRGs (VTN, KCNN4, IFNG, ENTPD1, COL1 A1, and AGER) by minimizing partial likelihood deviance (Fig. 5A, B). We then applied multivariate Cox regression analysis to these six genes to construct the final prognostic model. A risk score for PR-FRGs was calculated for each sample based on their expression levels and coefficients (Fig. 5C). Samples were divided into high-risk and low-risk groups based on the median risk score for downstream analyses. (Fig. 5J). The prognostic performance of the PR-FRG model was validated via ROC curves, which demonstrated its predictive capability over 1-, 3-, and 5-year timeframes. The AUC values for the training set, internal testing set, and full dataset ranged from 0.673 to 0.755 (Fig. 5D–F), indicating robust performance. The high-risk group demonstrated a poorer survival prognosis than the low-risk group, as shown by K-M survival analysis (Fig. 5G). Notably, high-risk samples exhibited lower expression levels of ENTPD1, while the remaining five PR-FRGs (VTN, KCNN4, IFNG, COL1 A1, and AGER) were upregulated (Fig. 5J). We further validated the prognostic value of our model using an external GEO dataset (GSE29609), as shown in Fig. 5H, I. The risk score was integrated with clinical factors such as histologic grade, pathological stage, and TMN classification to construct a nomogram for enhancing the PR-FRG model’s clinical relevance (Fig. 5K). Predicted survival probabilities demonstrated good consistency with actual survival outcomes at 1-, 3-, and 5-year intervals, as shown by the calibration plots (Fig. 5L). This comprehensive analysis underscores the potential clinical utility of our PR-FRG model in predicting ccRCC prognosis.

Fig. 5

Construction and validation of fibrogenesis-related prognostic model. A, B 6 candidate pFRGs identified by LASSO regression with ten-fold cross-validation, using the minimum partial likelihood deviance to determine the optimal lambda. C Coefficient profiles of 6 optimal pFRGs selected for the model. D-F ROC curves demonstrating the predictive performance of the prognostic model in the training (D), test (E) and whole (F) sets. G K–M survival showing a different prognosis in high- and low-risk groups. H, I Model performance validation using an external GEO dataset with ROC (H) and K-M (I) curves. J Risk plot of survival status of ccRCC patients and expression profiles of 6 pFRGs between risk groups. K A nomogram for estimation of survival rates for ccRCC patients. L Calibration plots of the nomogram at 1, 3, and 5 years. AUC area under curve, HR hazards ratio

Insights into differential immunotherapy efficacy through variations in immune infiltration between risk groups

Renal Clear Cell Carcinoma is characterized by significant heterogeneity in tumor progression and the TME, necessitating risk stratification to guide prognosis and personalized treatment strategies. Significant disparities in tumor microenvironment composition and immunotherapy efficacy were revealed through risk score analysis of ccRCC patients. As shown in Fig. 6A, high-risk patients had elevated Grade, advanced Stage, and poorer prognosis compared to low-risk patients (p < 0.001), with N classification showing no notable differences. However, T and M classifications are markedly higher in the high-risk group (p < 0.001). These results collectively indicate that high-risk patients have more advanced disease and poorer clinical outcomes, demonstrating the predictive validity of the risk score for patient stratification. Differences in TME cell composition between high- and low-risk groups are further depicted in the heatmap, with most immune cells enriched in the high-risk group except for neutrophils, which were predominant in the low-risk group. Figure 6B from ESTIMATE analysis supports these findings, showing elevated immune and ESTIMATE scores, reduced tumor purity, and no significant stromal score differences in the high-risk group. These observations highlight increased immune infiltration in high-risk patients.

Fig. 6

TME, immunotherapy and drug sensitivity analysis between high- and low- risk groups. A Immune landscape obtained by various algorithms with clinical features. B Comparison of the distribution of stroma, immune, ESTIMATE score and tumor purity. C Higher abundance of fibroblasts in high-risk group. D TIDE, Dysfintion, Exclusion and MDSC scoring by TIDE algorithm. E Differences in immune subtype between risk groups. F Immunophenoscore (IPS) comparison between different risk groups. G Classical immune checkpoint genes expression level. H Mutation landscapes. I Potential drug screening for different risk groups and correlation between IC50 and risk scores

Notably, both MCPcounter and xCell algorithms reveal that the high-risk group has a higher abundance of fibroblasts and CAFs (Fig. 6C), consistent with increased fibrosis, which is another evidence of fibrosis contributing to worse prognosis.

To further elucidate the complexity of the TME, TIDE scores were calculated for both groups (Fig. 6D). The high-risk ccRCC patients demonstrated significantly higher TIDE scores, reflecting increased immune dysfunction, as evidenced by higher Dysfunction scores and enriched MDSC levels, despite no differences in Exclusion scores. This suggests that while the high-risk group exhibits greater immune cell infiltration, the effectiveness of these immune cells is severely compromised due to dysfunction and exhaustion, limiting their ability to effectively eliminate tumor cells.

Interestingly, IPS predictions indicate greater immunotherapy sensitivity in the high-risk group (Fig. 6E), which aligns with elevated expression of classical immune checkpoint genes such as PDCD1 (PD-1), CTLA-4, and LAG3. (Fig. 6F). These findings suggest that, despite the immune dysfunction indicated by high TIDE scores, Patients in the high-risk group, characterized by abundant immune cell infiltration and increased immune checkpoint expression, may respond more effectively to immune checkpoint inhibitors targeting PD-1/PD-L1 and CTLA-4 pathways..

Immune subtype analysis further stratified patients based on the immune landscape (Fig. 6G). The high-risk group exhibited greater proportions of C1 (wound healing), C2 (IFN-gamma dominant), and C6 (TGF-beta dominant) subtypes compared to the low-risk group. Subtype C1 is often associated with high stromal content and angiogenesis, while C2 represents a pro-inflammatory TME with high immune cell infiltration. Conversely, C6, dominated by TGF-beta signaling, reflects a highly immunosuppressive TME enriched with fibrosis and immune exclusion. These subtype distributions align with the high fibrosis levels and immunosuppressive TME features in the high-risk group.

TMB analysis reveals that while the overall mutation frequency is similar between the two groups, the high-risk group shows a slightly higher frequency of disruptive mutations, such as frame-shift (Fig. 6H). This suggests greater genomic instability in the high-risk group, which may contribute to increased neoantigen load.

Lastly, drug sensitivity analysis was conducted for key targeted therapies used in advanced ccRCC (Fig. 6I). Axitinib, Sorafenib, and Erlotinib are used targeted therapies for advanced renal cell carcinoma, particularly ccRCC [36, 37]. Sensitivity to VEGFR-targeted therapies, such as Axitinib and Sorafenib, was significantly higher in the low-risk group, as indicated by lower IC50 values. However, the high-risk group demonstrated greater sensitivity to Erlotinib, an EGFR inhibitor, with reduced IC50 values. Risk scores were negatively correlated with IC50 values for Axitinib and Sorafenib but positively correlated with IC50 values for Erlotinib. These results indicate that high-risk patients are less likely to benefit from VEGFR-targeted therapies but may respond more effectively to EGFR-targeted agents, emphasizing the role of risk scores in personalized therapy.

In all, our risk stratification effectively captures the clinical, molecular, and immune heterogeneity in ccRCC, providing valuable insights into prognosis, tumor microenvironment characteristics, and guiding personalized therapeutic strategies.

Context-dependent expression and prognostic association of NRP1 in ccRCC

Neuropilin-1 (NRP1) has emerged as a multifunctional molecule critically involved in fibrosis, immune regulation, and tumor progression, making it a compelling target in both cancer biology and fibrotic disease [38]. In our single-cell RNA-seq analysis, NRP1 was found to be expressed across multiple cell types, including endothelial cells, fibroblasts, and tumor epithelial cells (Fig. 2A), highlighting its diverse roles within the ccRCC TME. These findings prompted us to further investigate role of NRP1 in ccRCC.

To explore the clinical significance of NRP1 in ccRCC, we first performed a pan-cancer analysis using TCGA data to assess NRP1 expression and its association with prognosis across various cancer types. NRP1 is overexpressed in various tumor types compared to normal tissues (Figure S6 A), with particularly high expression observed in breast cancer (BRCA), glioblastoma (GBM), and stomach adenocarcinoma (STAD). NRP1 showed strong diagnostic performance in ccRCC, with an AUC of 0.829 (CI 0.782–0.876), as demonstrated by ROC analysis (Figure S6B). Interestingly, high NRP1 expression is associated with improved prognosis in only ccRCC (KIRC) (Figure S6 C), which contrasts with previous reports of its pro-tumorigenic role in ccRCC, prompting our further investigation [39,40,41].

To clarify the paradox between the favorable prognostic association of NRP1 in TCGA and its reported tumor-promoting roles, we performed multi-omics analyses incorporating spatial transcriptomics and single-cell RNA-seq (GSE159115). These results confirmed and extended our initial findings from the GSE171306 dataset, showing that NRP1 is predominantly expressed in non-malignant stromal cells—including endothelial cells, pericytes, and fibroblasts—while its expression in tumor epithelial cells remains relatively low (Fig. 7A–H).

Fig. 7

Multi-omics characterization of NRP1 expression patterns and its context-dependent prognostic value in ccRCC. A H&E-stained tissue section from ccRCC patient for spatial transcriptome. B Cell type annotation of spatial transcriptomics (ST) data using deconvolution analysis. C Spatial expression map of NRP1 overlaid on the ST data. D, E UMAP plots of the ST dataset displaying cell clusters (D) and NRP1 expression levels across cell types (E). F,G UMAP plots of an independent scRNA-seq dataset (GSE171306) showing annotated cell types in ccRCC (F) and NRP1 expression (G). H Violin plots showing the distribution of NRP1 expression across different cell types. I Correlation analysis between NRP1 expression and the inferred proportion of major cell types across TCGA samples, based on BayesPrism deconvolution. J Scatter plot with marginal density curves showing a strong inverse correlation between NRP1 expression and tumor cell proportion (Spearman R = −0.640, p < 0.001). K Box plot comparing NRP1 expression levels between samples with high versus low tumor cell proportion. L, M Kaplan–Meier survival analysis stratified by NRP1 expression in Low-Tumor (L) and High-Tumor (M) groups

Since bulk RNA-seq data represent aggregated expression from both tumor and stromal compartments, we hypothesized that the favorable prognostic signal in TCGA may be primarily driven by NRP1 expression in the stroma. To test this, we performed BayesPrism-based deconvolution by integrating cell-type composition inferred from single-cell data into the TCGA bulk RNA-seq cohort. The results revealed that NRP1 expression was positively associated with stromal cell proportions, particularly pericytes and endothelial cells, and negatively correlated with tumor cell content Fig. 7I, J). Moreover, NRP1 levels were significantly lower in samples with high tumor content (Fig. 7K). Importantly, Kaplan–Meier survival analysis stratified by tumor content demonstrated that NRP1 exhibited a favorable prognostic effect in TME-rich (Low-Tumor) samples (p = 0.00026), while this association disappeared in tumor-rich (High-Tumor) samples (p = 0.54) (Fig. 7L, M). These findings suggest that the prognostic value of NRP1 is context-dependent and likely reflects its stromal expression, rather than tumor-intrinsic function, prompting us to further investigate the role of NRP1 within tumor cells through functional assays.

Clinical validation of NRP1 expression

We examined NRP1 expression in a clinical cohort of 20 ccRCC patients. Tumor tissues exhibited significantly higher NRP1 expression than adjacent normal tissues, as revealed by IHC staining (Fig. 8A, Figure S6D), further confirmed by quantitative mean optical density (OD) values (Fig. 8B). We also analyzed NRP1 expression across different tumor grades and found an increasing trend from Grade 2 to Grade 4 (Fig. 8C). A statistically significant difference in expression was observed between Grade 2 and Grade 4 tumors (p < 0.01). The corresponding representative IHC images are shown in Figure S6E. These findings indicate that NRP1 has clinical significance in ccRCC pathogenesis.

Fig. 8

Effects of NRP1 knockdown on ccRCC cell proliferation and apoptosis. A, B IHC analysis of ccRCC patient samples showing significantly higher NRP1 expression in tumor tissues. C NRP1 expression across different tumor grades in ccRCC. D NRP1 expression across various ccRCC cell lines, with selection of high-expression models 786-O and Caki-1. E Validation of siRNA knockdown efficiency, identifying siNRP1-3 as the optimal sequence. F CCK-8 assay results showing reduced cell viability post-NRP1 knockdown. G EdU staining indicating decreased cell proliferation. H Flow cytometry results showing G1 phase arrest. I Apoptosis analysis revealing increased apoptotic rates in siNRP1 cells. J Western blot analysis of apoptosis-related proteins

Effects of NRP1 knockdown on cell proliferation and apoptosis in ccRCC cells

To investigate NRP1’s function in ccRCC, we selected two cell lines with high NRP1 expression, 786-O and Caki-1, for further study (Fig. 8D). Among three siRNA sequences targeting NRP1, siNRP1-3 achieved the highest knockdown efficiency in both cell lines and was selected for subsequent experiments (Fig. 8E).

We evaluated cell proliferation via CCK-8 and EdU assays. Figure 8F demonstrates a significant reduction in cell viability in 786-O and Caki-1 cells transfected with siNRP1 over 24, 48, and 72 h compared to controls (p < 0.01). EdU staining corroborates these findings, showing decreased EdU incorporation in NRP1 knockdown cells, indicating reduced proliferation rates (Fig. 8G).

Flow cytometry was conducted to evaluate the effects of NRP1 on cell cycle distribution and apoptosis. NRP1 knockdown induced G1 phase arrest in both cell lines (Fig. 8H), suggesting that NRP1 may facilitate cell cycle progression. Additionally, Fig. 8I demonstrates a significant increase in apoptosis rates following NRP1 knockdown. Western blot analysis reveals that NRP1 knockdown upregulates the expression of the pro-apoptotic protein BAX and downregulates the anti-apoptotic protein BCL2, supporting its role in promoting cell survival. Additionally, the increase in cleaved caspase-3 and the reduction in CDK4 and PCNA expression indicate enhanced apoptosis and inhibited cell proliferation upon NRP1 knockdown (Fig. 8J). In conclusion, The knockdown of NRP1 in ccRCC cell lines suppresses proliferation, causes G1 phase arrest, and triggers apoptosis, highlighting its value as a therapeutic target.

Inhibition of NRP1 reduces cell migration, invasion, and EMT

Using Transwell and wound-healing assays, we examined the effects of NRP1 on ccRCC cell migration and invasion. NRP1 knockdown significantly reduced the migratory and invasive capabilities of 786-O and Caki-1 cells (Fig. 9A, B). Quantification of cell migration and invasion confirms these findings (Fig. 9C, D).

Fig. 9

NRP1 Knockdown suppresses ccRCC cell migration, invasion, EMT, and in vivo tumor growth. A, B Transwell migration and invasion assays showing significantly reduced migratory and invasive capacities post-NRP1 knockdown. C, D Quantitative analysis of migration and invasion. E, F Wound-healing assay demonstrating reduced cell motility with NRP1 knockdown. G GSEA plot showing the association between NRP1 and EMT. H Western blot analysis indicating EMT inhibition, with increased E-cadherin and decreased N-cadherin and Vimentin expression. I Representative images of tumors from sh-NC and shNRP1 xenograft models, showing reduced tumor sizes in the shNRP1 group. J Tumor growth curves showing significantly lower tumor volume in the NRP1 knockdown group. K Tumor weight measurements indicating reduced tumor mass post-knockdown. L H&E staining and IHC staining results showing reduced NRP1, Ki67 and CD31 expression, decreased Vimentin, and increased TUNEL-positive cells

Wound-healing assays revealed a significant reduction in migration rates over 48 h in NRP1-knockdown cells (Fig. 9E, F), suggesting a role for NRP1 in enhancing ccRCC cell motility. Subsequently, GSEA analysis of the TCGA cohort further showed a positive association between NRP1 and EMT (Fig. 9G). Consistently, NRP1 knockdown increased E-cadherin expression while reducing N-cadherin and Vimentin levels (Fig. 9H), indicating its role in promoting EMT and invasive behavior in ccRCC cells.

In vivo effects of NRP1 knockdown on tumor growth

To validate these findings in vivo, Xenograft assays in nude mice were performed using 786-O and Caki-1 cells with stable NRP1 knockdown. BALB/c-nu mice were injected with either 786-O/Caki-1 LV-control cells or 786-O/Caki-1 LV-sh-NRP1 cells to establish xenograft models. Figure 9I shows representative images of tumors from the NRP1 knockdown (shNRP1) and control groups (sh-NC) in both cell lines, with significant reductions in tumor volume observed over 21 days (Fig. 9J). Tumor weight was also significantly lower in the shNRP1 group, indicating that NRP1 knockdown effectively inhibits tumor growth in vivo (Fig. 9K).

Histological analysis of tumor tissues with H&E staining in Fig. 9L revealed that tumors from the sh-NC group exhibited dense cellular architecture with large, irregular nuclei and high cell density, consistent with aggressive tumor characteristics. In contrast, shNRP1 tumors displayed a looser structure and reduced cellularity, indicating decreased proliferative activity. IHC analysis further confirmed effective NRP1 knockdown, as shown by diminished NRP1 staining in shNRP1 tumors compared to controls. Additionally, Ki67 and CD31 expression levels were significantly lower in the NRP1 knockdown group, suggesting reduced proliferation and angiogenesis. Vimentin expression was also reduced, while TUNEL-positive cells were increased, indicating enhanced apoptosis and suppression of EMT in the absence of NRP1. These results collectively highlight NRP1’s role in promoting ccRCC tumor growth, angiogenesis, and cellular invasion in vivo.

NRP1 promotes ccRCC progression via the TGF-β/SMAD pathway

Considering that NRP1, TGFB1, and TGFBI genes exhibit high expression abundance in single-cell analyses of the ccRCC tumor microenvironment (Fig. 2A, B), and given the established link between NRP1 and the TGF-β pathway in other cancers [42, 43], we investigated the correlation between NRP1 and TGF-β signaling components in ccRCC. GSEA results of TCGA cohort demonstrate NRP1 positively correlating with the TGF-β (-SMAD) signaling pathway (Fig. 10A).

Fig. 10

NRP1 promotes ccRCC progression through the TGF-β/SMAD signaling pathway. A GSEA plot of the association between NRP1 and the TGF-β in HALLMARK, KEGG genesets. B The protein levels of TGF-β in the cell supernatants of ccRCC cells with or without NRP1 knockdown. C Western blot analysis of SMAD2 and SMAD3 phosphorylation levels showing reduction post-NRP1 knockdown. D, E Transwell migration and EdU proliferation assays showing partial recovery of migration, invasion, and proliferation capacities with TGF-β treatment in NRP1 knockdown cells. F, G Wound-healing and viability assays confirming partial restoration of motility and viability. H Flow cytometry analysis showing reduced apoptosis with TGF-β treatment. I ELISA results for TGF-β levels. J Western blot indicating partial restoration of SMAD2 and SMAD3 phosphorylation levels following TGF-β treatment. K Proposed mechanism of NRP1 in promoting ccRCC progression

ELISA analysis demonstrates that NRP1 knockdown significantly reduces TGF-β protein levels, along with Western blot analysis showing decreased phosphorylation of SMAD2 and SMAD3 (Fig. 10B, C), suggesting that NRP1 may promote ccRCC progression via TGF-β/SMAD pathway activation.

To confirm the involvement of TGF-β signaling in NRP1’s pro-tumorigenic effects, we treated NRP1-knockdown cells with exogenous TGF-β. TGF-β treatment partially restored the migratory, invasive, and proliferative capacities that were reduced in siNRP1 cells, as assessed by Transwell and EdU assays (Fig. 10D, E). Wound-healing and cell viability assays (Fig. 10F, G) further supported these findings, showing that TGF-β treatment increased cell motility and viability in NRP1-knockdown cells. Additionally, flow cytometry reveals that TGF-β treatment reduced apoptosis rates in siNRP1 cells, indicating a role for TGF-β in enhancing cell survival (Fig. 10H). ELISA results for TGF-β levels, along with Western blot analysis of SMAD pathway components (Fig. 10I, J) show that TGF-β treatment partially restored phosphorylation of SMAD2 and SMAD3, supporting the hypothesis that NRP1 exerts its tumorigenic effects in ccRCC through activation of the TGF-β/SMAD pathway (Fig. 10K).

Discussion

The role of fibrosis in cancer is highly intricate and remains insufficiently elucidated. In ccRCC, fibrosis manifests as intratumoral fibrosis (ITF), extracellular matrix (ECM) remodeling, cancer-associated fibroblasts (CAFs), and epithelial-mesenchymal transition (EMT) [16]. CAFs, which originate from various sources, contribute to RCC progression by interacting with tumor cells, promoting immunosuppression, and enhancing ECM deposition. They are linked to poor prognosis in ccRCC [44]. Although ITF is not directly linked to prognosis, it correlates with factors like Fuhrman grade and lymphovascular invasion [45]. Interestingly, early-stage protective fibrosis, especially in peritumoral capsules (PCs), may restrict tumor spread and improve prognosis by acting as a barrier to tumor-immune interactions, with lower PC fibrosis associated with higher recurrence and metastasis risk [46]. Given the complex role of fibrosis in ccRCC, a deeper understanding of fibrosis-associated genes is crucial for uncovering the mechanisms driving ccRCC progression and identifying potential therapeutic targets.

To address this, we employed an integrated bioinformatics approach, combining bulk RNA sequencing, single-cell RNA sequencing (scRNA-seq), and spatial transcriptomics to explore fibrosisgenesis-related features in ccRCC. Through Cox regression analysis, we initially identified a set of pFRGs significantly correlated with patient prognosis, revealing key players involved in the fibrotic landscape of ccRCC. To further elucidate the cellular and spatial distribution of these genes, scRNA-seq and spatial transcriptomics were used to map gene expression within specific cell populations across the tumor and surrounding tissues. This analysis highlighted distinct fibrotic cell types and regions that contribute to disease progression, underscoring the spatial complexity of the fibrotic response in ccRCC.

Using non-negative matrix factorization (NMF), we stratified ccRCC patients into two molecular subtypes based on pFRG expression. These subtypes displayed distinct fibrotic profiles, immune infiltration patterns, and clinical outcomes. The C1 cluster, characterized by increased fibrotic activity, higher immune infiltration, and enrichment in inflammatory and EMT-related pathways, exhibited significantly worse overall survival compared to the C2 cluster. Conversely, the C2 cluster displayed higher enrichment of metabolic pathways, suggesting a less aggressive tumor phenotype. GSVA analysis of fibrosis-related pathways further highlighted that C1 had higher enrichment scores for fibroblast growth factor signaling, ECM organization, and EMT, corroborating its pro-fibrotic and tumor-promoting characteristics.

Furthermore, leveraging machine learning, we identified a set of 12 pFRGs significantly correlated with patient outcomes. A robust prognostic model that effectively categorized ccRCC patients into high-risk and low-risk group, demonstrating high predictive accuracy. High-risk patients demonstrated significantly worse clinical outcomes, characterized by higher histologic grades, advanced pathological stages, and worse TMN classifications. Notably, the distinct TME features and increased fibrotic activity observed in the high-risk group underscore the prognostic and mechanistic relevance of pFRG-based stratification compared to the low-risk group. TME analysis results suggest that, despite the high TIDE scores indicating immune suppression, the abundant immune infiltration and elevated immune checkpoint expression in high-risk patients may render them more responsive to checkpoint inhibitors, particularly combination therapies targeting PD-1/PD-L1 and CTLA-4 pathways. However, effective immunotherapy may require strategies to alleviate immune exhaustion and modulate the fibrotic TME. This bioinformatics-driven framework offers a comprehensive approach to understanding the complexity of TME in ccRCC and identifying patient-specific targets for therapeutic intervention.

Among the pFRGs identified, Neuropilin-1 (NRP1) stands out due to its significant role in fibrosis and tumor progression. Originally known for its function in axonal guidance through Semaphorin signaling, NRP1 has been recognized for its broader contributions to vascular development, immune regulation, and cancer [38]. By interacting with ligands such as VEGF165 and TGF-β, NRP1 facilitates angiogenesis, immune modulation, and fibrotic processes, including hepatic cirrhosis and pulmonary fibrosis [47,48,49]; targeting NRP1 in renal tubular epithelial cells (TECs) has recently been shown to reduce kidney injury and fibrosis via TGF-β and TNF-α signaling [50]. In cancer, NRP1 expression generally correlates with aggressive behavior and poor prognosis, promoting immune evasion by enhancing Treg function in the TME, supporting angiogenesis, and increasing invasiveness via interactions with ligands like Semaphorins [38]. Our pan-cancer analysis corroborates these findings, highlighting NRP1 as a potential immunotherapy target. In ccRCC specifically, NRP1 enhances dedifferentiation, proliferation, invasion, and metastasis through pathways like VEGF and Sonic Hedgehog (Shh) [41]. On CD8 + T cells, NRP1 binds tumor-derived Semaphorin-3 A (SEMA3 A), disrupting cytoskeletal function and inhibiting T cell mobility, thereby limiting tumor infiltration [51]. Additionally, NRP1 promotes radiation resistance and immune suppression in RCC cells, underscoring its therapeutic potential for enhancing immune-based therapies [40].

Paradoxically, our TCGA analysis showed that high NRP1 expression correlated with improved survival in ccRCC. To investigate this contradiction, we integrated spatial and single-cell transcriptomics with BayesPrism-based deconvolution. These analyses revealed that NRP1 is mainly expressed in stromal cells—including endothelial cells, fibroblasts, and pericytes—rather than tumor epithelial cells. Its expression was positively correlated with stromal content and negatively correlated with tumor content. Moreover, survival analyses showed that NRP1’s favorable prognostic effect was limited to TME-rich samples and absent in tumor-rich cases. This suggests that the prognostic signal in bulk datasets primarily reflects stromal, not tumor-intrinsic, NRP1 expression. This compartment-specific function is supported by Morin et al. [52], who reported that perivascular NRP1 expression was independently associated with improved survival in RCC, potentially through VEGFR2-mediated endothelial-stromal crosstalk.

In contrast, our in vitro and in vivo experiments, which specifically target tumor-cell–intrinsic NRP1, demonstrated that NRP1 promotes proliferation, migration, and invasion of ccRCC cells. Furthermore, NRP1 knockdown in xenograft models significantly reduced tumor growth, supporting its oncogenic role in tumor cells. These findings align with previous reports identifying NRP1 as a key promoter of RCC progression [39,40,41]. Notably, xenografts, while partially mimicking the TME, lack human stromal and immune compartments. Thus, the tumor-suppressive phenotype following NRP1 knockdown reflects the loss of tumor-intrinsic NRP1 function, not stromal effects. These findings collectively reinforce the context-dependent nature of NRP1, with divergent roles in tumor cells versus the surrounding microenvironment.

TGF-β/SMAD is a well-established driver of fibrosis and plays a dual role in ccRCC by promoting both fibrotic ECM remodeling and EMT [42]. Mechanistically, we found that NRP1 modulates the TGF-β/SMAD pathway—a key driver of fibrosis and EMT in ccRCC. NRP1 knockdown reduced SMAD2/3 phosphorylation, suggesting its role as a TGF-β co-receptor that enhances fibrotic signaling. EMT, tightly regulated by TGF-β, facilitates tumor cell motility and invasiveness by shifting epithelial cells toward a mesenchymal phenotype [53]. In our study, NRP1 knockdown restored E-cadherin expression and reduced N-cadherin and Vimentin levels, indicating a mesenchymal-to-epithelial transition and loss of invasive potential.

Beyond EMT, the interplay between TGF-β signaling, fibrosis, and immune suppression defines a tumor-promoting microenvironment. Fibrotic ECM stiffening impairs immune cell infiltration, while CAFs secrete immunosuppressive cytokines and recruit regulatory T cells [54, 55]. Our single-cell and spatial analyses further demonstrated that fibrosis-related genes—including NRP1—are differentially expressed across CAFs, immune cells, and endothelial populations, suggesting that fibrosis and immune suppression are co-regulated processes contributing to ccRCC aggressiveness.

Given its distinct expression patterns across tumor and stromal compartments and its dual roles in promoting tumor progression and shaping the microenvironment, NRP1 holds promise as both a context-dependent prognostic biomarker and a potential therapeutic target in ccRCC. Its tumor-intrinsic expression is associated with enhanced proliferation, migration, and EMT via TGF-β/SMAD signaling, whereas its stromal enrichment appears linked to favorable prognosis in TME-rich contexts. This functional divergence highlights the need for cell-type–specific interpretation and intervention strategies.

Moving forward, further studies are warranted to delineate the molecular mechanisms regulating NRP1 in different cellular contexts, identify its upstream activators and downstream effectors, and explore the therapeutic efficacy of selectively targeting NRP1 in tumor cells. Advances in delivery systems—such as nanoparticle-based approaches, tumor-specific promoters, or ligand-guided antibodies—may help achieve selective inhibition of tumor-cell NRP1 while preserving stromal function. In parallel, integration of NRP1 expression into prognostic models stratified by tumor-stromal composition could enhance patient risk classification and guide precision therapy. Together, these directions support the translational potential of NRP1 as a biomarker-driven therapeutic node in fibrotic and immune-modulated ccRCC.

Despite the comprehensive nature of our study, several limitations should be acknowledged. First, reliance on public databases like TCGA limits generalizability, underscoring the need for validation using larger, multi-center cohorts—especially those incorporating immunohistochemical (IHC) data from different institutions—to confirm the reproducibility and clinical relevance of our observations. Additionally, while in vitro and in vivo models elucidated NRP1’s role, these models do not fully replicate the complexity of the human tumor microenvironment, particularly immune interactions, making patient-derived models or clinical studies essential for further confirmation. Furthermore, while we focused on the dual role of NRP1 in tumor cells versus stromal compartments, its cell-type–specific function across various TME components—including immune cells, fibroblasts, and endothelial cells—requires further dissection through single-cell resolution and lineage-tracing approaches. Lastly, while NRP1 is proposed as a therapeutic target, we did not explore combinatorial strategies such as co-targeting TGF-β/SMAD pathways, and preclinical studies are needed to evaluate potential therapeutic synergies in ccRCC treatment.

Conclusions

Our study highlights the prognostic and therapeutic value of fibrogenesis-related genes in ccRCC, focusing specifically on the dual role of NRP1 in fibrosis and tumor progression. Through an integrative bioinformatics analysis, we identified key pFRGs and underscored the critical role of fibrosis in shaping the TME and driving ccRCC progression. The identification of molecular subtypes based on pFRG expression highlights the heterogeneity of fibrotic responses in ccRCC and offers a framework for patient stratification and personalized therapy. To enhance prognostic accuracy, We developed a robust machine learning-based model for risk classification, which accurately reflected clinical aggressiveness and immune contexture. Among these genes, NRP1 emerged as a central mediator, promoting proliferation, migration, and invasion of ccRCC cells via EMT and activation of the TGF-β/SMAD signaling pathway. In contrast, its expression in stromal components correlated with favorable prognosis in bulk transcriptomic data. These findings collectively reveal the context-dependent function of NRP1 and underscore its value as both a prognostic indicator and a therapeutic target in ccRCC.