Transcriptome-based network analysis related to regulatory T cells infiltration identified RCN1 as a potential biomarker for prognosis in clear cell renal cell carcinoma

Abstract

Background

Regulatory T cells (Tregs) play a critical role in shaping the immunosuppressive microenvironment within tumors. Investigating the role of Tregs in Clear cell renal cell carcinoma (ccRCC) is crucial for identifying prognostic markers and therapeutic targets for ccRCC.

Methods

Weighted gene co-expression network analysis (WGCNA) was utilized to pinpoint modules related to Treg infiltration in TCGA-KIRC samples. Following this, consensus clustering was employed to derive two clusters associated with Treg infiltration in ccRCC. A prognostic model was then developed using the gene module associated with Treg infiltration. We then evaluated the ability of the prognostic model to predict ccRCC overall survival and demonstrated that RCN1 can be used as a target to predict ccRCC prognosis.

Results

We deduce that the two clusters associated with Treg infiltration exhibit distinct compositions of the immune microenvironment, pathway activations, prognosis, and drug sensitivities commonly utilized in ccRCC treatment. Furthermore, a 7-gene model risk score, developed based on ccRCC Treg infiltration, proved to be a reliable prognostic marker in both training and validation cohorts. Additionally, survival analysis indicated that RCN1 serves as a reliable prognostic factor for ccRCC. Single-cell sequencing analysis revealed that RCN1 is predominantly expressed in tumor cells. A pan-cancer analysis highlighted that RCN1 is linked with poor prognosis and the activation of inflammatory response pathways across various cancers.

Conclusion

We developed a prognostic model associated with Treg infiltration, which facilitates the clinical categorization of ccRCC progression. Moreover, our findings underscore the significant potential of RCN1 as a ccRCC biomarker.

Introduction

Renal cell carcinoma (RCC) accounts for 2% of all global cancer diagnoses, with its incidence rate gradually increasing. RCC is characterized by a variety of histological subtypes, each with distinctive pathological features. Clear cell renal cell carcinoma (ccRCC), which represents around 75% of all RCC instances, arises from the proximal tubular cells within the nephron and is the most common subtype [1,2,3]. Surgical removal at an early stage is considered highly advantageous for patients with ccRCC, yet about 30% of patients may face recurrence or metastasis after the tumor is removed. Antiangiogenic treatments, such as Sunitinib and Pazopanib, have shown effectiveness in managing metastatic RCC [4, 5]. however, challenges like resistance and recurrence persist, with some patients demonstrating innate resistance to these targeted therapies [6, 7]. Thus, the identification of early prognostic biomarkers is crucial for improving the treatment outcomes for ccRCC.

In the dynamic interplay between cancer progression and host immunity, regulatory T cells (Tregs) emerge as pivotal players, especially in the context of RCC. Characterized by their immunosuppressive capabilities, Tregs contribute significantly to the tumor microenvironment by promoting immune tolerance and enabling tumor escape from immune surveillance [8, 9]. In RCC, a malignancy marked by its resistance to traditional chemotherapy and its intricate interactions with the immune system, the role of Tregs has garnered considerable attention. These cells not only facilitate an immunosuppressive milieu conducive to cancer growth and dissemination but also impact the efficacy of emerging immunotherapeutic strategies [10, 11]. Understanding the mechanisms through which Tregs exert their influence in RCC is crucial for unveiling potential therapeutic targets. This involves exploring the balance between Treg-mediated immunosuppression and tumor-directed immune responses, the impact of Tregs on the effectiveness of current RCC treatments, and the potential for therapeutic modulation of Treg activity as a means to enhance anti-tumor immunity [12, 13]. Given the increasing incidence of RCC and the urgent need for more effective treatment modalities, dissecting the complex role of Tregs in renal carcinogenesis and progression holds promise for improving patient outcomes in this challenging malignancy.

In the realm of renal cell carcinoma (RCC), the search for reliable prognostic models is paramount for advancing patient management and tailoring personalized treatment strategies. RCC, with its diverse histological subtypes and complex biological behavior, poses significant challenges in predicting outcomes and selecting optimal therapeutic approaches [14,15,16]. Despite advances in diagnostic techniques and the development of targeted therapies, the prognosis for RCC patients, particularly those with advanced or metastatic disease, remains unpredictable. The advent of prognostic models based on genetic, molecular, and clinical parameters has opened new avenues for the early identification of high-risk patients, guiding treatment decisions, and improving survival outcomes [17,18,19,20]. These models integrate various prognostic factors to generate a more accurate prediction of disease progression, recurrence, and patient survival, thereby facilitating a more stratified approach to RCC treatment [21,22,23,24]. As the understanding of RCC's underlying molecular mechanisms continues to expand, the development and validation of these prognostic models become increasingly crucial for enhancing clinical decision-making and advancing the overall management of RCC.

Reticulocalbin-1 (RCN1), a protein found in the endoplasmic reticulum, is integral to maintaining calcium balance and preventing apoptosis triggered by ER stress [25]. Elevated RCN1 expression, linked to cancer development and invasion, has been noted across various cancers, including glioblastoma, non-small cell lung cancer, renal cell carcinoma, and oral squamous cell carcinoma [26,27,28,29]. Furthermore, increased RCN1 levels are correlated with resistance to sorafenib in hepatocellular carcinoma and resistance to doxorubicin in uterine cancer cells [30, 31]. In contrast, reducing RCN1 expression has been shown to curb cell growth and induce apoptosis through the activation of AKT and PTEN pathways in prostate cancer cells [32], highlighting its potential as a therapeutic target.

In our study, the CIBERSORT algorithm was applied and revealed a notable correlation between increased Treg presence and adverse patient outcomes in TCGA-KIRC cohort. Next, we used weighted gene co-expression network analysis (WGCNA) to identify gene modules that are closely related to Treg infiltration in ccRCC. Building upon these findings, our study ventured into stratifying ccRCC patients based on genes related to Treg cells infiltration, unveiling two distinct clusters with differing prognostic outcomes. We then found that there were distinct differences in immune infiltration, pathway activation, and predictability of drug effectiveness between the two clusters. Moreover, we constructed a prognostic model based on a gene module related to Treg cell infiltration, and this model's risk score demonstrated good performance in predicting the prognosis of ccRCC patients in both the training and validation sets. Subsequently, RCN1 was identified among the 7 genes used to construct the prognostic model and could predict overall survival in ccRCC. Furthermore, A comprehensive analysis across multiple cancer types revealed that RCN1 is associated with unfavorable outcomes and the activation of pathways related to inflammatory responses in a wide range of cancers. Accordingly, a prognostic model based on genes associated with Treg cell infiltration was developed and validated. This model offers potential prognostic utility for clinical applications, highlighting the possible importance of RCN1 in cancer progression.

Methods

Dataset source and data pre-processing

The analyses involved patients from two ccRCC cohorts (E-MTAB-1980 [33], Braun ccRCC 2020 [34]) and TCGA-KIRC. Patients without survival information and RNA sequencing (RNA-seq) data were excluded from the analysis. The clinical data and transcriptome expression data of E-MTAB-1980 and Braun ccRCC 2020 were downloaded from ArrayExpress (https://www.ebi.ac.uk/biostudies/arrayexpress/studies) and related original publication. Transcriptome data, measured as TPM (transcripts per million) values, along with clinical information, were obtained from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/) through the use of the 'TCGAbiolinks' package in R. The raw data for the single-cell RNA sequencing (scRNA-seq) study, identified as phs002065.v1.p1 [35], were retrieved from the Single Cell Portal (https://singlecell.broadinstitute.org/single_cell/study/SCP1288/tumor-and-immune-reprogramming-during-immunotherapy-in-advanced-renal-cell-carcinoma#study-summary). Processing of the scRNA-seq data was conducted using the 'Seurat' package in R, following the guidelines provided in the accompanying tutorial. In summary, cells exhibiting gene expression levels below 300 or above 6500, as well as those with mitochondrial gene expression exceeding 10%, were filtered out. The normalization and scaling of raw data were accomplished through the SCTransform function, preceding the execution of principal component analysis (PCA). To mitigate batch effects in the isolated scRNA-seq raw data, the "Harmony" package in R was employed. Distinct cell clusters within each scRNA-seq dataset were discerned through unsupervised clustering analysis and the unified manifold approximation and projection (UMAP) technique. Subsequently, each identified cell cluster was annotated utilizing known markers specific to cell types.

Weighted gene co-expression network analysis (WGCNA)

Gene co-expression networks were developed utilizing the TCGA-KIRC dataset through the 'WGCNA' package in R. Pearson’s correlation coefficient was calculated for each gene pair to create a similarity matrix. The 'WGCNA' package includes a power function capable of transforming this similarity matrix into an adjacency matrix. For all soft thresholds (b) yielding an R^2 greater than 0.9, we selected the automatically determined value of b (b = 5) as suggested by the WGCNA's pickSoftThreshold function. In line with the 'WGCNA' guide, we set the network merging criterion at a height of 0.25. We adhered to the default parameters provided by WGCNA for subsequent analyses unless stated otherwise.

Gene set variation analysis (GSVA) and assessment of immune cell infiltration

The GSVA analysis was conducted utilizing the 'GSVA' package in R to derive gene set enrichment scores. These scores facilitated the comparison of pathway enrichments, distinguishing between upregulated and downregulated pathways in the group with high-risk scores compared to those with low-risk scores. The gene sets for GSVA were acquired from the Molecular Signatures Database (MSigDB). Furthermore, the quantification of immune cell infiltration within TCGA-KIRC datasets was achieved through the application of the CIBERSORT algorithm, which was based on normalized gene expression data.

Acquisition of Treg cell-related clusters

A gene module linked to Treg cell infiltration was identified, followed by a univariate Cox regression analysis on the genes within this module. Out of these, 779 genes demonstrating a significant association with survival (P < 0.05) in the univariate analysis were selected for further analysis using the 'ConsensusClusterPlus' package in R, aimed at clustering patients from the TCGA-KIRC cohort. The determination of the optimal number of clusters (K) was based on the consensus value and the cumulative distribution function analysis, resulting in an optimal K value of 2. The Nearest Template Prediction (NTP) technique offers an efficient approach for making classification predictions based on a signature gene list and a test dataset, allowing for the assessment of predictive confidence in gene expression data for each patient. In our research, the top 50 upregulated genes from each cluster within TCGA-KIRC cohort were utilized to predict cluster membership in the E-MTAB-1980 cohort using the 'MOVICs' package in R.

Establishment of Treg-related prognostic model

A prognostic model associated with Treg cells infiltration was developed using two clear cell renal cell carcinoma (ccRCC) cohorts: TCGA-KIRC served as the training dataset, while E-MTAB-1980 was utilized for validation. From the previously mentioned 79 prognosis-related genes, key genes pertinent to Treg cell infiltration and patient prognosis were identified using two machine learning approaches: the least absolute shrinkage and selection operator (LASSO) logistic regression and the random forest method. Initially, the 'glmnet' package in R was utilized to select 21 prognostic genes (including CIB1, LGALS2, FIS1, EIF4EBP1, AUP1, ISG15, NDUFV1, GAMT, TRAPPC6A, NME4, NUTF2, B3GAT3, STAP2, HSBP1, PAXX, RCN1, TRAPPC2L, NUDT14, YDJC, PLAC9, and INMT) using LASSO regression with tenfold cross-validation. Subsequently, the random forest algorithm, implemented through the 'randomForestSRC' package in R, was employed to refine the selection to 10 features with prognostic significance. These features include ISG15, TRAPPC6A, HSBP1, STAP2, S100A11, RCN1, LGALS2, CIB1, BOP1, and PSENEN. Finally, seven common genes (CIB1, LGALS2, ISG15, TRAPPC6A, STAP2, HSBP1, and RCN1) were obtained to build the multivariate Cox regression models (both using stepwise regression). The risk score was calculated as follows: 0.411 × (CIB1 expression)-0.110 × (LGALS2 expression) + 0.273 × (ISG15 expression)-0.258 × (TRAPPC6A expression)-0.272 × (STAP2 expression)-0.450 × (HSBP1 expression) + 0.376 × (RCN1 expression).The risk score for both the TCGA-KIRC and validation cohorts was calculated using the identical model score threshold. Patients were categorized into low- and high-risk groups based on the median value of the risk scores. The R 'survival' package was then utilized to compare the overall survival (OS) differences between these groups.

Nomogram construction

A nomogram was constructed to incorporate both the risk score, age and tumor grade, utilizing the regplot function within the 'rms' package in R. To evaluate the predictive accuracy of our model, a receiver operating characteristic (ROC) curve was generated. Furthermore, to depict the variance between the predicted outcomes of our model and the actual observed patient survival, both calibration curves and decision curve analysis (DCA) were plotted.

Chemotherapeutic response prediction

Predictions of chemotherapeutic responses for each cluster were made utilizing the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/), the largest available pharmacogenomics database. Four chemotherapeutic agents frequently used in ccRCC treatment—sorafenib, axitinib, suunitinib, and Gemcitabine—were chosen for detailed analysis. These predictions were carried out using the 'pRRophetic' package in R. Ridge regression was employed to estimate the half-maximal inhibitory concentration (IC50) for the samples, and the prediction's accuracy was verified through tenfold cross-validation against the GDSC's training dataset. The default settings were used for all parameters, with the exception of specifying the tissue type as 'kidney'.

Pan-cancer analysis

The normalized mRNA expression data and clinical details for the TCGA pan-cancer cohorts (referenced in Supplementary Table 1) were acquired from the UCSC Xena Browser (https://xenabrowser.net/datapages/). The prognostic significance of RCN1 across different types of cancer prognosis was evaluated using univariate Cox regression analysis and Kaplan–Meier modeling. RCN1 expression data, treated as a continuous variable, was incorporated into the univariate Cox regression. Additionally, RCN1 expression was dichotomized for Kaplan–Meier survival analysis, with the division threshold determined by the "surv-cutpoint" function from the "survminer" R package (version 0.4.9). Both the log-rank p-value from the Kaplan–Meier method and the hazard ratio (HR) along with a 95% confidence interval (95%CI) were calculated. The results were visually summarized in a heatmap format.

Statistical analysis

Survival differences across groups were analyzed with Kaplan–Meier curves and the log-rank test. Pearson and Spearman analyses were utilized to compute correlation coefficients. For normally distributed continuous data, Student’s t-test was applied, while the Mann–Whitney U test was employed for data not following a normal distribution. When conducting comparisons across more than two groups, nonparametric data were analyzed using the Kruskal–Wallis test, and parametric data were examined with one-way ANOVA. These statistical procedures were carried out in R software (version 4.3.1), with P values < 0.05 deemed to indicate statistical significance.

Results

Tregs infiltration was associated with poor prognosis in ccRCC

Initially, the CIBERSORT algorithm was employed to evaluate the proportion of immune cell infiltration in patients. Within TCGA-KIRC cohort, patients exhibiting a higher level ofTregs infiltration demonstrated a poorer prognosis (Fig. 1A). Moreover, patients with advanced tumor stage and higher tumor grade exhibit higher Treg cells infiltration (Fig. 1B-C). Given the observation that an increased presence of Treg cells was linked to a poorer prognosis, we conducted a weighted gene co-expression network analysis (WGCNA) to identify the module associated with Treg cells infiltration (Fig. 1D). A soft threshold power of β = 5, achieving a scale-free fit index (R^2) of 0.90, was chosen to establish a scale-free network (Figure S1). Moreover, the correlation heatmap indicates that, within TCGA-KIRC, the yellow module showed a negative correlation with patient’s overall survival time (r = -0.14, P = 0.001) and a positive correlation with Treg cells infiltration (r = 0.19, P = 7e-06) (Fig. 1E).

Fig. 1

Tregs infiltration was associated with poor prognosis in ccRCC. A The Kaplan–Meier analysis demonstrated the relationship between the infiltration of Tregs and overall survival (OS) within the TCGA KIRC cohorts. Patients were categorized into "high" and "low" groups according to the median score for Tregs infiltration, derived from CIBERSORT analysis. Differences in Treg cell infiltration in different clinical stages (B) and different clinical grades (C) in TCGA-KIRC cohort. D Dendrogram of cluster modules analyzed by WGCNA results, where each color represents a different co-expression module and the top branches represent genes. E Correlation analysis between various phenotypes and co-expression modules revealed that genes within the yellow module showed a positive correlation with Treg cell infiltration and negative correlation with overall survival

Stratification of ccRCC based on Tregs infiltration

Utilizing the genes from the yellow module along with survival data from the TCGA-KIRC dataset, a univariate Cox regression analysis was conducted, identifying 79 genes associated with overall survival (OS) in TCGA-KIRC (Figure S2A). Through the use of the R ConsensusClusterPlus package for consistent clustering within the TCGA-KIRC dataset based on these 79 prognostic genes, two distinct clusters were identified: Cluster 1, comprising 354 cases, and Cluster2, consisting of 169 cases (Fig. 2A and Figure S2B). PCA analysis also showed that the gene expression pattern in these two clusters were distinct (Fig. 2B). Survival analysis further revealed that individuals in the Cluster2 group exhibited a worse prognosis compared to those in the Cluster1 group (Fig. 2C). The heatmap illustrated variations in the expression patterns of the 79 genes between Cluster1 and Cluster2 (Fig. 2D). Additionally, the chi-squared test revealed that patients in Cluster2 had higher tumor grades, more advanced stages, and poorer survival outcomes compared to those in Cluster1 (Fig. 2E). Next, we verified the existence of these two clusters in the E-MTAB-1980 cohort using the NTP algorithm (Fig. 2F). Patients in Cluster2 still has a worse prognosis than patients in Cluster1 (Fig. 2G).

Fig. 2

Stratification of ccRCC based on Tregs infiltration. A The consensus clustering plot shows that the TCGA-KIRC samples are divided into two clusters. B Principal component analysis of two clusters. C Kaplan–Meier survival analysis shows the difference in OS between the two clusters in TCGA-KIRC cohort. D Heatmap of expression patterns of 79 genes in two clusters in the TCGA-KIRC cohort. E Donut plot of survival status, tumor grade, and tumor stage for the two patient groups. F Classification predictions for E-MTAB-1980 cohort were carried out utilizing upregulated genes specific to TCGA-derived clusters, employing the Nearest Template Prediction (NTP) algorithm. G Kaplan–Meier survival analysis shows the difference in OS between the two clusters in E-MTAB-1980 cohort

Immune infiltration and pathway enrichment differences in different clusters

The diversity of immune cell infiltration across these two clusters was analyzed using the ssGSEA, MCPcounter and CIBERSORT algorithms. A heatmap displayed the profile of tumor-infiltrating immune cells in TCGA-KIRC patients (Fig. 3A). According to all three algorithms, Treg cells infiltration in Cluster2 was significantly increased compared to Cluster1; Moreover, analysis by three of these algorithms revealed that CD8 + T cell infiltrations were notably lower in Cluster2. Taken together, these findings indicate that Cluster2 exhibits a distinct immune phenotype compared to Clusters1, characterized by higher Treg cells infiltration and lower levels of immune activation. To illustrate the activation of signaling pathways within each cluster, GSVA enrichment scores were computed using gene sets from the HALLMARK pathways in MSigDB. Notably, Cluster2, in contrast to Cluster1, was distinguished by the enrichment of tumor-promoting pathways, such as Hypoxia, EMT and PI3K AKT Mtor SIGNALING pathway. While, Cluster1 was characterized by enrichement of metabolic pathways, including fatty acid metabolism and bile acid metabolism (Fig. 3B).

Fig. 3

Immune infiltration and pathway enrichment differences in different clusters. A A heatmap illustrates the disparities in immune cell infiltration among different subtypes, analyzed using the ssGSEA, MCPcounter, and Cibersort algorithms. Statistical discrepancies were evaluated using the Kruskal-Walli’s test. B Differences in HALLMARK pathway activities scored by GSVA between Cluster1 and Cluster2 groups. The red histogram represents pathways upregulated by Cluster2 and the blue represents pathways upregulated by Cluster1

Differences in immunotherapy efficacy and drug sensitivity between the two clusters

Then, the expression levels of HLA family genes and immune checkpoint markers across these two subtypes were examined within both TCGA-KIRC and E-MTAB-1980 datasets. The expression of HLA family genes and immune checkpoint genes in Cluster2 is lower than that in Cluster1 (Fig. 4A-B). Next, we verified the existence of these two clusters in the Braun ccRCC 2020 cohort using the NTP algorithm (Fig. 4C). Additionally, we categorized treatment responses into a binary model and observed that the proportion of patients exhibiting stable or progressive disease in Cluster2 was greater compared to Cluster1 (Fig. 4D). Moreover, Survival analysis further revealed that individuals in the Cluster2 group also exhibited a worse prognosis compared to those in the Cluster1 group (Fig. 4E). Given the potential for patients with ccRCC to develop resistance to various drugs, we assessed the responsiveness of the two clusters to four commonly used ccRCC therapeutic agents: gemcitabine, sorafenib, axitinib, and sunitinib. We developed a predictive model utilizing ridge regression based on the GDSC cell line dataset and verified its predictive reliability through tenfold cross-validation. We then calculated the half-maximal inhibitory concentration (IC50) for samples within the TCGA-KIRC dataset using the predictive models for these drugs. The findings indicated that patients in Cluster 1 exhibited a higher sensitivity to these therapeutic treatments (Fig. 4F).

Fig. 4

Differences in immunotherapy efficacy and drug sensitivity between the two clusters. Boxplots depict variations in the expression of immune-related and immune checkpoint genes across different subtypes within the TCGA-KIRC (A) and E-MTAB-1980 cohort (B). C Classification predictions for Braun ccRCC 2020 cohort were carried out utilizing upregulated genes specific to TCGA-derived clusters, employing the Nearest Template Prediction (NTP) algorithm. D The percentage of patients with response to immunotherapy in different subtypes in Braun ccRCC 2020 cohort. SD, stable disease; PD, progressive disease; CR, complete response; PR, partial response. E Kaplan–Meier survival analysis shows the difference in OS between the two clusters in Braun ccRCC 2020 cohort. F Differences in chemotherapy responsiveness between the two clusters in the TCGA-KIRC cohort

Construction and validation of Treg infiltration-related prognostic model

For the purpose of enhancing the clinical relevance of genes associated with Treg infiltration in prognosis assessment, we utilized two machine learning methods to pinpoint key genes among all 79 Treg infiltration-related prognostic genes. LASSO and random forest algorithms identified 21 and 10 key genes, respectively (Figure S3A-B). Subsequently, seven genes that were common to both algorithms: CIB1, LGALS2, ISG15, TRAPPC6A, STAP2, HSBP1, and RCN1 were chosen to develop a multivariate Cox regression mode (Figure S3C). Furthermore, a seven-gene risk model was employed to compute the risk score for each patient, subsequently categorizing patients into high-risk and low-risk subgroups. Patients with higher risk score had poorer prognosis in TCGA-KIRC cohort (Fig. 5A). Simultaneously, the AUC (Area Under the Curve) values at 1, 2, 3, and 5 years indicated that the risk score was significantly predictive of overall survival (OS) in TCGA-KIRC cohort (Fig. 5B). Figure 5C displays the risk score distributions, overall survival time associated mRNA expression profiles for the TCGA-KIRC cohort. Protective mRNAs (HSBP1, STAP2, LGALS2 and TRAPPC6A) exhibited higher expression levels in the low-risk group, whereas the other mRNAs (RCN1, CIB1 and ISG15) showed higher expression in the high-risk group. Additionally, the high-risk group have more deaths compared to the low-risk group (Fig. 5C). Moreover, patients with advanced tumor stage and higher tumor grade exhibit higher model risk score (Fig. 5D-E). Patients with higher risk score has more Treg cells infiltration (Fig. 5F). To confirm the prognostic importance of the risk score, the same calculation was applied to derive the Treg infiltration-related risk score in a validation cohort (E-MTAB-1980). The risk score demonstrated a comparable prognostic value in these cohorts, along with strong predictive accuracy for overall survival (Fig. 5G-H). And the risk score distribution plot also shows that the E-MTAB-1980cohort has the same expression pattern as TCGA-KIRC cohort (Fig. 5I).

Fig. 5

Construction and validation of Treg infiltration-related prognostic model. A Kaplan–Meier survival analysis shows the difference in OS between high-risk and low-risk score groups in the TCGA-KIRC cohort. B AUC values of risk score for 1-, 2-, 3- and 5-year OS in TCGA-KIRRC cohort. C The distribution of risk score, vital status and the prognostic gene expression patterns in the TCGA-KIRC cohort. Differences in model risk score in different clinical stages (D) and different clinical grades (E) in TCGA-KIRC cohort. F Differences in Treg cell infiltration across risk score groups in the TCGA-KIRC cohort. G Kaplan–Meier survival analysis shows the difference in OS between high-risk and low-risk score groups in the E-MTAB-1980 cohort. H AUC values of risk score for 1-, 2-, 3- and 5-year OS in E-MTAB-1980 cohort. I The distribution of risk score, vital status and the prognostic gene expression patterns in the E-MTAB-1980 cohort

Establishment of the nomogram model based on risk score

To improve the predictive capability of the aforementioned risk scores, a nomogram model was developed by integrating the risk score, age and tumor grade through multivariable Cox regression analysis (Fig. 6A). Calibration curves for disease-specific survival (DSS) at 3 and 5 years indicated a strong concordance between the predicted survival probabilities and the actual outcomes, highlighting the nomogram's reliability in forecasting survival. Furthermore, decision curve analysis (DCA) was conducted, revealing that the nomogram's prognostic accuracy surpassed that of the individual variables (Fig. 6B-C). Additionally, the outcomes of our investigation revealed that the nomogram's predicted AUC outperformed the risk score's AUC in both the TCGA-KIRC and E-MTAB-1980 cohorts (Fig. 6D-E). We applied two distinct methods to pinpoint drug candidates showing increased sensitivity in patients with high Treg infiltration-related risk scores, utilizing data on drug responses derived from CTRP and PRISM. Initially, a differential drug response analysis was conducted to compare high vs. low Treg infiltration-related risk score groups, aiming to discover compounds with notably lower estimated AUC values in the high-risk group (log2FC > 0.10). Subsequently, we identified compounds exhibiting a negative correlation between AUC values and Treg infiltration-related risk scores by assessing the Spearman correlation coefficient (Spearman’s r < -0.40 for CTRP and -0.35 for PRISM). This approach identified 11 compounds from PRISM (including rigosertib, MK-2461, vindesine, vinblastine, dolastain-10, talazoparib, verubulin, NVP-AUY922, topotecan, rubitecan and echinomycin) and four from CTRP (including uprosertib, leptomycin B, paclitaxel, topotecan, CR-1-31B, SB-743921 and BI-2536), all of which displayed lower estimated AUC values in the high-risk group and a negative correlation with the Treg infiltration-related risk score (Fig. 6F).

Fig. 6

Establishment of the nomogram model based on risk score. A Nomogram based on Treg infiltration-related risk score, age and tumor grade. B Disease-specific survival calibration curves at 3, and 5 years. C Nomogram decision curve analysis, Treg infiltration-related risk score, age and tumor grade. D Kaplan–Meier survival analysis shows the difference in OS between high-nomo and low-nomo score groups in the TCGA-KIRC cohort (left). AUC values of nomo score for 1-, 2-, 3- and 5-year OS in TCGA-KIRRC cohort (right). E Kaplan–Meier survival analysis shows the difference in OS between high-nomo and low-nomo score groups in the E-MTAB-1980 cohort (left). AUC values of nomo score for 1-, 2-, 3- and 5-year OS in E-MTAB-1980 cohort (right). F Analyses of Spearman correlation and differential responses for six compounds derived from CRTP, including comparisons of AUC values between high- and low-nomo score groups in response to these compounds. Additionally, Spearman correlation and differential response analyses for eleven compounds sourced from PRISM, accompanied by evaluations of AUC value differences between high- and low-nomo score groups in response to these four compounds. ‘*’ indicates P-value ≤ 0.05, ‘**’ indicates P-value ≤ 0.01, ‘***’ indicates P-value ≤ 0.001

Single-cell sequencing analysis based on RCN1 expression

Next, we performed survival analysis on the seven genes of the risk model in TCGA and E-MTAB-1980 cohort. Interestingly, only high expression of RCN1 indicated worse prognosis both in TCGA-KIRC and E-MTAB-1980 cohort (Figure S4). Therefore, we selected RCN1 for further analysis. We collected a single-cell sequencing data set of ccRCC and found that RCN1 is mainly expressed in tumor cells (Fig. 7A-B). Moreover, we divided patients into high RCN1 group and low RCN1 group based on the median RCN1 expression value of tumor cells. It was found that tumors with high expression of RCN1 had lower infiltration of tumor killer cells, such as NK cells, NKT cells and CD8 + T cells, but have more tumor-promoting cells, such as tumor-associated macrophages and Treg cells (Fig. 7C). Then we used Cellchat to analyze the differences in ligand receptors of various cell types in patients in the high and low RCN1 groups. Overall, tumor cells with high RCN1 expression have more ligand receptor interactions with other cell types in the tumor microenvironment than tumor cells with low RCN1 expression. Interestingly, compared to tumor cells with low RCN1 expression, tumor cells with high RCN1 expression have more ligand receptor interactions with Treg cells, such as SPP1-CD44, MIF-(CD74 + CXCR4), MIF-(CD74 + CD44), FN1-CD44 and APP-CD74 (Fig. 7D).

Fig. 7

Single-cell sequencing analysis based on RCN1 expression. A UMAP visualization maps of various cell clusters in phs002065.v1.p1 scRNA-seq cohort. B Feature plots show the expression of RCN1 in phs002065.v1.p1 scRNA-seq cohort. C Histogram of compositional differences among cell types from patients with high and low RCN1 expression. D Cell interaction analysis using CellChat. Bubble plots show all significant ligand-receptor pairs that contribute to the signaling sending from high and low-RCN1 malignant cellsto other cell types in phs002065.v1.p1 scRNA-seq cohort

Pan-cancer analysis of RCN1 in TCGA cohorts

Next, the expression of RCN1 in tumor versus normal tissues for 20 cancer types within the TCGA cohort was also investigated, revealing that RCN1 was upregulated in 60% of the cancers, such as BLCA, UCEC, HNSC, KIRP, COAD, LUSC, KIRC, LIHC, BRCA, KICH, LUAD, ESCA, and STAD (Fig. 8A). Moreover, we examined the association between RCN1 and the Hallmark pathways across the TCGA pan-cancer cohort, finding that RCN1 positively correlated with inflammatory pathway and EMT, including TNFA signaling pathway and inflammatory response signaling pathways, across various cancer types (Fig. 8B). Furthermore, the study assessed the impact of RCN1 expression on the survival prognosis in 32 cancer types, indicating that high RCN1 expression correlated with poorer survival outcomes in over 7 cancer types, including GBM, HNSC, KIRC, KIRP, LGG, LUAD and MESO (Fig. 8C).

Fig. 8

Pan-cancer analysis of RCN1 in TCGA cohorts. A Variations in RCN1 expression between tumor and normal tissues across 20 cancer types within the pan-cancer TCGA cohort. ‘*’ indicates P-value ≤ 0.05, ‘***’ indicates P-value ≤ 0.001, ‘****’ indicates P-value ≤ 0.0001. B Enrichment analysis of HALLMARK pathways between tumor tissues with high and low RCN1 across 33 cancer types in TCGA cohorts, NES normalized enrichment score in the GSEA algorithm, FDR false discovery rates. C Overview of how RCN1 expression correlates with overall survival (OS) among 32 cancer types in the TCGA pan-cancer cohort

Discussion

The incidence of renal cell carcinoma (RCC) is on the rise globally, with a mortality rate of approximately 20%. RCC is primarily categorized into three main histological subtypes: clear cell RCC (ccRCC), papillary RCC, and chromophobe RCC, with ccRCC being the most common (accounting for 70–80% of all cases) and the deadliest form [36]. Despite technological advancements leading to the development of new diagnostic methods and therapeutic approaches that have improved early-stage ccRCC patient outcomes, the overall survival rate remains unsatisfactory [37]. Therefore, identifying novel biomarkers for the prognosis and therapeutic targeting of ccRCC is crucial.

Recent research has highlighted the prognostic significance of Treg cells across different cancer types. The predominance of Treg cells is linked not just to adverse outcomes in a range of tumors but also to the creation of an immunosuppressive tumor microenvironment [38, 39]. Increased presence of Treg cells in tumor sites is associated with poorer prognosis in cancer patients [40, 41]. Remarkably, Tregs derived from tumors possess stronger suppressive abilities than their naturally occurring equivalents, making them a strategic target for improving outcomes in cancer therapy [42,43,44]. Despite the challenges in manipulating Treg cells to accurately regulate immune responses, there has been a growing focus on investigating Treg cell modifications in clear cell renal cell carcinoma (ccRCC) recently. Nonetheless, there is a significant lack of research on the unique roles and diagnostic potential of Treg cell in ccRCC. This deficiency highlights the urgent need for more comprehensive studies on the Treg subpopulations within ccRCC to uncover their roles in tumor diagnostics and treatment strategies. Our study discovered that an increased infiltration of Treg cells, estimated through the CIBERSORT algorithm, correlates with a poorer prognosis in patients with ccRCC. Furthermore, the tumor grading and clinical staging of patients are positively associated with Treg cell infiltration. These results are consistent with results from several previous studies.

Consensus clustering, a sophisticated method integrating multiple clustering algorithms to identify stable and consistent cluster structures, has significantly impacted the study of cancer heterogeneity [45]. In cancer research, the application of consensus clustering has enabled the stratification of patients into groups with different prognostic outcomes, sensitivity to treatments, and risk factors, thus offering a more nuanced understanding of tumor biology. Moreover, it facilitates the discovery of novel biomarkers and therapeutic targets by delineating the molecular profiles associated with specific cancer subtypes [46]. In ccRCC, A study indicates that ccRCC can be divided into four subtypes based on gene expression patterns, and these four subtypes are able to differentiate patients with varying prognoses and sensitivities to targeted drugs [47]. In another study, Grigory Andreevich Puzanov identified the most aggressive ccRCC subtype associated with metastasis using genes related to coagulation (FGA, FGG) and genes associated with changes in tumor immune characteristics (ENAM, IGFBP1, IL6) [48]. In our study, Given the significant role of Treg cells in determining clinical outcomes and their contribution to the immunosuppressive TME, it was deduced that a gene module indicative of Treg cells presence in ccRCC could be utilized to formulate a prognostic model. This model holds potential in forecasting clinical outcomes in ccRCC cases. The result of analysis indicated that a high Treg cells infiltration correlates with decreased survival rates in the TCGA-KIRC cohort. To further assess the efficacy of Treg cells as a prognostic indicator for ccRCC, patients were divided into two groups based on the expression of genes within an Treg cells-related module. This division revealed disparities in overall survival (OS) and clinical characteristics between the groups.

In our study, we explored the differential expression of HLA family genes and immune checkpoint markers across the two identified subtypes within both the TCGA-KIRC and E-MTAB-1980 datasets. Our analyses revealed that Cluster2 exhibited lower expression levels of these genes compared to Cluster1, suggesting a diminished immune response capacity that may influence treatment efficacy and disease progression. The validation of these clusters in the Braun ccRCC 2020 cohort, employing the NTP algorithm, further solidifies our findings and underscores the robustness of our subclassification. The investigation into treatment responses highlighted a significant difference in outcomes between the two clusters. Specifically, Cluster2 demonstrated a higher proportion of patients with stable or progressive disease, indicating a potential for poorer treatment efficacy. This observation was corroborated by survival analysis, where individuals in Cluster2 were found to have a worse prognosis than those in Cluster1. Such findings emphasize the importance of cluster-based stratification in predicting patient outcomes and guiding treatment decisions. Our exploration into drug resistance and responsiveness to common ccRCC therapeutic agents—gemcitabine, sorafenib, axitinib, and sunitinib—provides critical insights into the potential for personalized treatment approaches. By developing a predictive model based on the GDSC cell line dataset and validating its accuracy through tenfold cross-validation, we established a framework for estimating drug sensitivity. The calculated half-maximal inhibitory concentrations (IC50) for samples in the TCGA-KIRC dataset revealed that Cluster 1 patients are more susceptible to these treatments, suggesting a potential for tailored therapeutic strategies that could improve patient outcomes.

Moreover, seven genes were identified and used to develop a multivariate Cox regression model. Analysis across two ccRCC cohorts indicated that patients were stratified into high- and low-risk categories based on Treg infiltration-related risk scores. Patients in the high-risk category showed poorer overall survival than those in the low-risk category. Tregs, known for inhibiting CD8 + T cell cytotoxicity, facilitating B-cell proliferation, and encouraging tumor growth, are implicated in the adverse outcomes seen in the high-risk group. This underlines the critical impact of Treg infiltration and the consequential suppression of immune responses as key factors driving the poorer prognosis in high-risk patients compared to their low-risk counterparts. Therefore, the risk score derived from the prognostic model related to Treg cells serves as a potential predictor of overall survival in ccRCC.

Many studies have indicated that CIB1, LGALS2, ISG15, TRAPPC6A, STAP2, HSBP1, and RCN1 play significant roles in the progression of cancer, sensitivity to drug treatment, and the survival prognosis of cancer patients [31, 49,50,51,52,53,54]. In our study, high expression of RCN1 was found to be associated with poorer prognosis in two ccRCC cohorts, hence we selected RCN1 for further investigation. Reticulocalbin 1 (RCN1) is a multifunctional protein residing in the endoplasmic reticulum (ER), belonging to the family of calcium-binding proteins. The expression of RCN1 has been documented across multiple tissues, highlighting its essential role in normal cellular physiology [55]. However, its involvement in disease processes, particularly in cancer, has garnered increasing attention. Research has demonstrated that RCN1 expression levels are altered in various types of cancer, suggesting a potential role in tumorigenesis and progression [29, 31, 56]. Elevated RCN1 expression has been associated with poor prognosis in several cancers, including colorectal cancer, prostate cancer, and ccRCC, implicating it as a potential biomarker for cancer diagnosis and prognosis [28, 32, 57]. In our study, single-cell sequencing analysis also showed that RCN1 is mainly expressed on tumor cells and patients with high RCN1 expression have higher Treg cell infiltration. Moreover, compared with tumor cells with low expression of RCN1, cellchat analysis showed that tumor cells with high expression of RCN1 were interacted with Treg cells through SPP1-CD44, MIF-(CD74 + CXCR4), and MIF-(CD74 + CD44). Previous studied showed that macrophage migration inhibitory factor (MIF) promotes tumor growth by increasing Tregs generation [58]. These results also reveal to a certain extent the relationship between high RCN1 expression and increased Treg cell infiltration, and further indicated that RCN1 may regulate Treg cell infiltration in ccRCC.

The investigation into the expression of Reticulocalbin 1 (RCN1) across various cancer types within the TCGA cohort has revealed significant insights into its role in oncogenesis and tumor progression. Our analysis showed that RCN1 was upregulated in a majority of cancers analyzed, indicating its pervasive role in tumoral environments. Specifically, an increased expression of RCN1 was observed in 60% of the cancers examined, including but not limited to BLCA, UCEC, HNSC, and KIRP. This upregulation suggests that RCN1 may play a fundamental role in the development and progression of these cancers, acting possibly as a facilitator of tumor growth and survival. Moreover, our examination of the association between RCN1 and various Hallmark pathways unveiled a positive correlation with pathways known to enhance tumor aggressiveness, such as the inflammatory pathway and epithelial-mesenchymal transition (EMT). The correlation with TNFα signaling and other inflammatory response pathways across multiple cancer types underscores the potential of RCN1 to modulate the tumor microenvironment, promoting conditions favorable for cancer progression. Furthermore, the impact of RCN1 expression on survival outcomes across 32 cancer types revealed a stark correlation between high RCN1 expression and poorer survival outcomes in over seven cancer types, including GBM and LUAD. This association highlights the prognostic potential of RCN1, suggesting that its expression level could serve as a biomarker for survival outcomes in a subset of cancers. The findings from this comprehensive analysis shed light on the multifaceted role of RCN1 in cancer biology. These insights necessitate further investigation into the mechanistic underpinnings of RCN1's role in cancer and its potential utility in prognostic assessments and therapeutic interventions. Given its significant association with both cancer progression and patient prognosis, RCN1 represents a promising target for the development of novel cancer therapies and management strategies.

Conclusion

Our research has developed a prognostic model related to Treg cells for forecasting OS in ccRCC, while also investigating the predictive significance of RCN1 in the context of ccRCC prognosis. We aim to enrich the existing knowledge on the impact of Treg cells on ccRCC's biology and its prognosis. Furthermore, we propose that RCN1 could serve as an innovative biomarker for predicting clinical outcomes in ccRCC patients.