Gene expression data (transcripts per million, TPM) and the relevant ccRCC prognostic and clinical data were obtained from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) databases. GSE40435, GSE53757, GSE121636, and GSE156632 were obtained from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). RPPA protein expression data were obtained from ucsc xena (https://xena.ucsc.edu/). The E-MTAB-1980, E-MTAB-3267 and E-MTAB-3218 datasets were obtained from the ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/). CheckMate 025 was obtained from (CM-025; NCT01668784) and CheckMate-010 from (CM-010; NCT01354431). The IMvigor210 cohort was downloaded from IMvigor210CoreBiologies. A schematic diagram of the study design is shown in Supplementary Fig. 1.

Molecular subtypes based on pyroptosis-related genes were identified by the “ConsensusClusterPlus” package. Repetitions were performed 1000 times to ensure the stability of results [11].

GSVA enrichment was performed using the “GSVA” R package [12]. The “limma” package was used to identify differentially expressed genes (DEGs) [13]. The selection of differentially expressed genes was an adjusted p-value less than 0.05 and an absolute value of log${}_2$ Fold Change $\ge$0.1.

Immune cells in the tumor microenvironment (TME) were quantified using the CIBERSORTx algorithm and MCP-counter [14, 15]. Tumor purity scores were estimated with the “ESTIMATE” package [16]. The Tumor Immune Estimation Resource (TIMER) database was used to evaluate the TME [17] (Supplementary Table 1).

The cancer immune cycle reflects the anticancer immune response and is comprised of 7 steps [18] (Supplementary Table 2). The relative activities of these steps determine the fate of tumor cells. The activity of each step was analyzed by single-sample gene set enrichment analysis (ssGSEA) based on gene expression of each samples [19].

DEGs (log${}_2$ Fold Change $>$0.5) were identified using a protein interaction algorithm and STRING analysis tool. Thirty central nodes with optimal connectivity were identified using the CytoHubba’s Maximal Click Centrality (MCC) function, and a prediction model was established using LASSO-Cox analysis [20]. Kaplan-Meier analysis was used to evaluate patient survival, and a time-dependent receiver operating characteristic (ROC) curve was used to evaluate the accuracy of the model.

Barcodes with 5000 total UMIs (Unique Molecular Identifiers), 200 genes expressed, or mitochondrial genes were filtered out. Highly expressed genes were selected by genes that retained at least 3 UMIs in at least 200 cells [21]. Cells were separated by coarse clustering using the scatter v1.12.2 package [22]. Perform Principal Component Analysis (PCA) on normalized data and correct for batch effects on the top 200 PCs using the Mutual Nearest Neighbors (MNN) method [23].

The human normal cortex/proximal tubule epithelial cell line HK-2 and two renal clear cell carcinoma cell lines, 769-P and Caki-1, were from the Chinese Academy of Sciences Cell Bank (Shanghai, China). The HK-2, 769-P and Caki-1 cell lines were cultured in Dulbecco’s modified Eagle’s medium (biosharp, Hefei China), RPMI 1640 Medium (biosharp, China) and McCoy’s 5a Modified Medium (biosharp, China), respectively. This was added with 10% fetal bovine serum (VivaCell BIOSCIENCES, Hefei, China), 100 U/mL penicillin and 100 µg/mL streptomycin (Gibco, New York, NY, USA). All cell lines were found to be mycoplasma-free using the MycAway™ Plus-Color One-Step Mycoplasma Detection Kit (Yeasen Biotechnology, Shanghai, China) and were authenticated shortly before use with a PCR technique (Procell Corporation, Wuhan, China).

RNA was isolated from paracancerous and tumor tissues (16 pairs) using TRIZOL reagent (Invitrogen, Waltham, MA, USA) according to the manufacturer’s protocol. Tumor tissues and paracancerous were collected from 16 ccRCC patients in the Department of Urology, Fourth Affiliated Hospital of Harbin Medical University. Informed consent was obtained and signed by patients for all tissues.The use of tissues was approved by the Ethics Committee of the Fourth Affiliated Hospital, Harbin Medical University. The ReverTraAce Qpcr RT Kit (Toyobo, Tokyo, Japan) was used for qRT-PCR experiments. The following primers were used for qRT-PCR: GSDMB: 5${}^{\prime }$-AGTCTTTGGGTTCGGAGGAT-3${}^{\prime }$ (F), 5${}^{\prime }$-CTGTCTGGGTCCTCCATGT-3${}^{\prime }$ (R); ACTB: 5${}^{\prime }$-CTTCCTTCCTGGGCATGG-3${}^{\prime }$ (F), 5${}^{\prime }$-GCCGCCAGACAGCACTGT-3${}^{\prime }$ (R).

The protein expression level in cell lines was determined by Western blot analysis. Protease inhibitors were used to isolate and lyse samples in RIPA buffer (Beyotime, Shanghai, China). Western blot was performed using rabbit polyclonal antibody against GSDMB (12885-1-AP, ProteinTech Group, Chicago, IL, USA) and mouse monoclonal antibody against $\beta$-actin (66009-1-Ig, ProteinTech Group, Chicago, IL, USA).

A total of 35 tissue specimens were collected from 35 patients in the Department of Urology, Fourth Affiliated Hospital of Harbin Medical University. The use of tissues was approved by the Ethics Committee of the Fourth Affiliated Hospital, Harbin Medical University. GSDMB protein expression in low- and high-grade ccRCC and in normal tissues was assessed by IHC. Paraffin-embedded samples were deparaffinized, rehydrated, and placed in citrate buffer at 98 °C for 15 min for antigen retrieval. They were then incubated with anti-GSDMB antibody (Proteintech, 1:200 dilution). The expression was then examined by DAB kit (ORIGIN, Beijing, China).

Correlation coefficients were calculated using Spearman analysis and distance correlation analysis. The log-rank test was used to determine the significance of differences between survival curves. ROC curves and the area under the curve (AUC) were obtained using the “pROC” and “timeROC” packages. Clinical characteristics were compared by chi-square or Fisher’s exact test.

Many clinical and genomic studies have shown that ccRCC is a high immune-infiltrating tumor type [24]. In the present study, the immune score for ccRCC was significantly correlated with poor overall survival (OS) and with tumor grade and stage (Fig. 1a, Supplementary Fig. 2), in accordance with the findings of a recent study [25]. Cytokines produced by pyroptotic cells recruit immune cells that subsequently infiltrate the TME, thereby having a tumor-promoting or suppressive role [26]. Pyroptosis therefor plays a significant role in the progression of ccRCC. The differential expression of pyroptosis-related genes (PRGs) between ccRCC and normal tissue is shown in Supplementary Fig. 3. Univariate cox analysis was used to identify genes associated with overall survival in ccRCC. Twelve genes (AIM2, CASP4, GSDMB, GSDME, NLRP1, NLRP6, NOD2, NLRP1, NLRP6, PYCARD, SCAF11, TIRAP and HMGB1) were selected for subsequent analyses due to having p 0.05. However, only NLRP6, SCAF11, TIRAP, and HMGB1 were potentially protective genes (i.e., hazard ratio (HR) 1) for ccRCC patients, with the remaining genes being associated with increased risk (HR $>$1) (Fig. 1b). A pyroptosis network was created to demonstrate the comprehensive landscape of PRG interactions, regulations and connections, together with their prognostic value in ccRCC patients (Fig. 1c). The expression of these PRGs (log${}_2$ FC $>$1) in the TME, as determined by single-cell RNA sequencing, is presented in Fig. 1d–f and in Supplementary Fig. 4.

Fig. 1.

Overview of pyroptosis-related genes in clear cell renal cell carcinoma (ccRCC). (a) Relationship between the Immune Score and tumor grade and stage (*p 0.05; **p 0.01; ***p 0.001; ****p 0.0001). (b) Forest plot presents pyroptosis-related genes associated with overall survival in ccRCC in The Cancer Genome Atlas (TCGA) database. (c) Interactions among pyroptosis-related genes (PRGs) in ccRCC. (The lines connecting PRGs indicate the interactions between them, while the width of the lines indicates the strength of the connection between them. Blue represents negative correlation and purple represents positive correlation, and the size of the point represents the strength of the p-value. Purple dots represent risk factors and green dots represent protective factors). (d–f) Expression of pyroptosis-related genes in the microenvironment of renal clear cell carcinoma based on single-cell RNA sequencing data.

The molecular mechanisms underlying the heterogeneity of ccRCC were studied with the aim of optimizing personalized treatment strategies. Using unsupervised clustering, two different patterns were identified based on the expression profiles of PRGs with prognostic value (log${}_2$ Fold Change $>$1). These were comprised of 330 samples in Cluster 1, and 196 samples in Cluster 2 (Fig. 2a, Supplementary Fig. 5). Survival analysis showed that patients in Cluster 2 had worse over survival than patients of Cluster 1 (p 0.001) (Fig. 2b). Biological functions were assessed by GSVA enrichment analysis between the two subtypes. Cluster 1 was associated with oncogenic pathways such as mTOR, Notch, MAPK, Wnt signalling, metabolic pathways, cell cycle, and DNA damage response. Cluster 2 contains more immune-related pathways, including antigen processing and presentation, and cytokine-cytokine receptor interactions, implying this subgroup may be involved in the tumor immunological process (Fig. 2d). Furthermore, Cluster 2 had higher immune and stromal scores compared with cluster 1, cluster 2 had higher immune and stromal scores Cluster 1 (Fig. 2c), implying the clusters had distinctive TME cell-infiltrating properties. Further investigation of the stromal composition revealed that Cluster 2 had more cancer-associated fibroblasts, while Cluster 1 had more endothelial cells. These results demonstrate that Cluster 1 tumors had increased angiogenesis, similar to the result from the TIMER database (Fig. 2e). Moreover, evaluation of the subsets at the protein level further support the clustering results in RPPA database (Supplementary Table 3). The AKT pS473, AMPKALPHA, c-KIT, EGFR pY1068, EGFR pY1173, MAPK pT202Y204, MEK1, SMAD1, SMAD4, CD31, VEGFR2, and STAT3 pY705 proteins were significantly expressed in Cluster 1 (Fig. 2f). These proteins participate in various carcinogenic processes and are also associated with angiogenesis. For example, CD31 is widely used as a vascular marker to assess tumor angiogenesis [27, 28], while VEGFR2 promotes angiogenesis by affecting the proliferation and migration of vascular endothelial cells [29]. Therefore, Cluster 1 was referred to as the “angiogenesis subtype” [30].

Fig. 2.

Cluster 1 shows features of an angiogenesis phenotype. (a) 526 ccRCC patients were divided into two subtypes based on the consensus clustering matrix. (b) Analysis of overall survival in two subgroups. (c) Tumor microenvironment (TME) score in the two clusters (***p 0.001). (d) Identification of common pathways in the two subtypes at the transcriptome level using gene set variation analysis (GSVA). (e) Comparison of the TME between the two subtypes at the transcriptional level based on the Tumor Immune Estimation Resource (TIMER) database (**p 0.01; ***p 0.001). (f) Comparison of proteins associated with the angiogenic signalling pathway between the two subtypes based on the RPPA database (*p 0.05; **p 0.01; ***p 0.001).

A recent study reported that high immune infiltration in ccRCC was associated with worse patient outcome [31]. The immune cells in the microenvironment between the two subtypes were evaluated based on CIBERSORT. Cluster 1 showed abundant infiltration of activated innate immune cells, including activated dendritic cells, M1 macrophages, naïve B cells, resting NK cells, mast cells and resting T4 memory cells. Although Cluster 2 showed significant immune cell infiltration, this subgroup was enriched with regulatory T (Treg) cells and CD8${}^{+}$ T cells (Fig. 3b). KEGG pathway enrichment analysis revealed that DEGs were primarily enriched in immune-related pathways, with cytokines and cytokine-receptor interactions being the most significantly enriched. Consequently, we investigated chemokines, tumor necrosis factor and interleukin in the TME of the two subtypes. Cytokine expression was significantly different between the subtypes. IL-4, IL-10 and IL-6 are crucial for immunosuppression and can promote cancer metastasis [32]. IL-6 can also promote MDSC differentiation, thus promoting immunosuppression [33]. CXCL10, CCL4, CCL5, CCL8, CCL19 and CXCL13 have been associated with worse survival of ccRCC patients [34, 35]. CCL5 recruits Treg cells into tumor masses [36]. Markers associated with T cell exhaustion, such as TIGIT, LAG3, LYN, and TOX, were significantly expressed in Cluster 2 (Fig. 3a). Although there were more infiltrating CD8${}^{+}$ T cells, the T cells in Cluster 2 were exhausted or non-functional, thus confirming our hypothesis. Single-cell RNA sequencing revealed heterogeneity among infiltrating CD8${}^{+}$ T cells in the ccRCC microenvironment, with most showing high expression levels for markers of exhaustion (Fig. 3c,d). Moreover, exhausted CD8${}^{+}$ T cells were associated with tumor stage [37]. Therefore, Cluster 2 was referred to as the “inflamed subtype” [38]. Anti-angiogenic therapy and immune checkpoint inhibitors, alone or in combination, can significantly improve the clinical prognosis of patients with advanced ccRCC. This is due to the high vascular density and the large number of immune cells infiltration in the microenvironment of ccRCC [39]. VEGF, VEGFR, EGF, EGFR, and KIT are targets for first-line anti-angiogenesis therapy. The expression of these targets and of immune checkpoints in the two subgroups was confirmed here (Fig. 3e). The results suggest the “angiogenesis subtype” may be more sensitive to targeted therapy, while the “inflamed subtype” may be more sensitive to immunotherapy.

Fig. 3.

Cluster 2 shows an inflamed phenotype. (a) The thermogram shows different expression for chemokines, interleukins, tumor necrosis factor and T cell exhaustion factors between Clusters 1 and 2 at the transcriptome level (*p 0.05; **p 0.01; ***p 0.001). (b) The abundance of infiltrating immune cells in the tumor microenvironment was analyzed separately for the two pyroptosis-related clusters using CIBERSORTx. (c) Characteristics of CD8${}^{+}$ T cells in the tumor microenvironment of ccRCC. (d) Percentage of exhausted CD8${}^{+}$ T cells in the TME (P1 represents the first patient, P2 means the second patient, and P3 means the third patient). (e) Expression levels of angiogenic genes and immune checkpoint genes at the transcriptome level between the two subtypes (*p 0.05; **p 0.01; ***p 0.001).

The antitumor immune response requires a series of stepwise events termed the cancer immune cycle. Neoantigens produced by oncogenesis are first captured by dendritic cells (DCs) (step 1). The antigens captured by MHCI and MHCII molecules are then presented to T cells (step 2). Effector T cells are then primed and activated (step 3). Finally, the activated effector T cells move to the tumor bed (step 4) and subsequently infiltrate (step 5). As a result, they selectively recognize cancer cells through T-cell receptor (TCR) interactions (step 6) and kill the target cancer cell (step 7) [18]. Differences between the two subtypes were observed here for the cancer-immune cycle (Fig. 4a). Cluster 2 had a higher score for step 7 due to more infiltration of CD8${}^{+}$ T cell. However, a high number of exhausted CD8${}^{+}$ T cells limits the recognition of cancer cells, and immune checkpoint blockade can re-awaken exhausted CD8${}^{+}$ T cells to elicit antitumor immunity [40]. The tumor mutation burden (TMB) is a biomarker of the immune checkpoint blockade response [41]. Commonly mutated genes in ccRCC were evaluated here in the two clusters. BAP1 mutation was more frequent in Cluster 2, while ATM and PBRM1 mutations were more frequent in Cluster 1 (Fig. 4b). The TMB was not statistically different between the two clusters (Supplementary Fig. 7). Tumors with BAP1 mutation have been associated with considerably worse over survival rates [42, 43]. Furthermore, BAP1-mutant tumors are associated with mTORC1 activation [44]. The clinical features of the two clusters were investigated next. Tumors from Cluster 2 had higher grade, stage and immune score, thus explaining their worse prognosis (Fig. 4c). We further assessed genes that were previously related to specific biological processes [45]. Cluster 1 contains more oncogenic pathways such as epithelial-mesenchymal transition (EMT) and Wnt signaling, as well as the cell cycle and angiogenesis pathway (Fig. 4c). Cluster 2 was enriched in immune pathways, including CD8${}^{+}$ T effectors, immune checkpoint, and Pan-F-TBRS pathways. Personalized treatment for ccRCC could therefore be offered based on the characteristics of significant vascular infiltration in Cluster 1, and of immunosuppression in Cluster 2. Anti-angiogenic drugs that target VEGF signalling pathways, including sorafenib and sunitinib, have been widely used to treat advanced carcinomas [46, 47]. In accordance with our conjecture, the Cluster 1 subgroup in the E-MTAB-3267 cohort showed a higher objective response rate to sunitinib than Cluster 2 (Fig. 4d). Based on the same genes, tumor samples in the PD-1 inhibitor-treated cohorts E-MTAB-3218 and CheckMate-010 were divided into two distinct subgroups, thus demonstrating the robustness of our cluster analysis (Fig. 4e,f). As predicted, Patients in cluster 2 had a higher reactivity to nivolumab than Cluster 1 patients.

Fig. 4.

Comprehensive analysis of clinical and mutational features of the two clusters, as well as their response to treatment. (a) Cancer immune cycle analysis in the two clusters. (b) The top 20 mutations found in the two subtypes of ccRCC. (c) The clinical features, mutation and functional characteristics of the subtypes were comprehensively analyzed (*p 0.05; **p 0.01; ***p 0.001). (d) Response of the two subtypes to sunitinib therapy (****p 0.0001). (e) The response to nivolumab treatment was assessed between the two subtypes using progression-free survival (CR means complete response; PR means partial response; SD means stable disease; PD means progressive disease) (****p 0.0001). (f) Response of the two subtypes to nivolumab therapy in CheckMate-010 cohort (**p 0.01).

A prognostic signature was developed to predict prognosis and the response to immune checkpoint blockade therapy. Univariate Cox regression analysis was performed using the MCC function to identify the top 30 central genes in the two subgroups associated with OS, followed by LASSO penalty Cox regression analysis of the significant variables (p 0.05). A 9-gene signature was then constructed based on the optimum $\lambda$ value. The risk score for each patient was calculated using the following formula: (0.184 $\times$ HJURP exp) + (0.329 $\times$ NUF2 exp) + (–0.203 $\times$ KIF15 exp) + (–0.390 $\times$ MELK exp) + (–0.101 $\times$ TPX2 exp) + (0.250 $\times$ PLK1 exp) + (0.561 $\times$ CDCA3 exp) + (0.050 $\times$ CTLA4 exp) + (0.050 $\times$ FOXP3 exp). Patients in the TCGA dataset were then divided into high- and low-risk groups according to the median score. Survival analysis revealed the low-risk group had a significantly better clinical outcome than the high-risk group (Fig. 5a; p 0.01). Using the same algorithm, risk scores were calculated for ccRCC patients classified as low- and high-risk subgroups in the E-MTAB-1980 dataset to further validate the prognostic signature. Kaplan-Meier survival curves showed the OS results for ccRCC patient subroups was consistent with those from the TCGA dataset (Fig. 5c; p 0.001). To investigate whether this signature could also predict responsiveness to PD-1 therapy, risk scores were calculated using the same algorithm in the CheckMate-025 and CheckMate-010 datasets for metastatic ccRCC patients previously treated with VEGF-directed therapy. The signature predicted both progression-free survival (PFS) and OS of metastatic ccRCC patients treated with nivolumab. The low-risk group showed better response to nivolumab treatment than the high-risk group (Fig. 5e,f). Moreover, in the IMvigor210 cohort the risk score also reflected objective response rates and PFS in urothelial carcinoma patients treated with the PD-L1 inhibitor atezolizumab (Supplementary Fig. 8). Univariate and multivariate Cox analyses were also used to evaluate whether the signature was an independent prognostic factor for ccRCC patients. Univariate analysis indicated the signature was strongly associated with OS in ccRCC patients from the TCGA dataset [hazard ratio (HR): 3.08 95% CI: 2.46–3.85, p 0.001]. Multivariate analysis further showed the signature was an independent predictor of OS (HR: 2.10, 95% CI: 1.63–2.70, p 0.001; Fig. 5g). The test dataset confirmed the signature was an independent predictor of OS for ccRCC patients in univariate (HR: 4.16, 95% CI: 2.46–7.03, p 0.001) and multivariate (HR: 3.62, 95% CI: 1.89–6.95, p 0.001) analyses (Fig. 5h). These results indicate the 9-gene signature may be useful for the clinical evaluation of prognosis. AUC analysis of cases in the TCGA database revealed the signature showed good accuracy, with values of 0.73, 0.72 and 0.75 at 1, 5 and 10 years of follow-up, respectively (Fig. 5b). The AUC values calculated from the E-MTAB-1980 dataset were 0.89, 0.82 and 0.75 at 1, 5 and 10 years, respectively (Fig. 5d).

Fig. 5.

Construction and validation of a subtype-based prognostic signature. (a) Survival analysis of high- and low-risk groups in patients from the TCGA database. (b) Receiver operating characteristic (ROC) curves demonstrated the accuracy of the risk score in the TCGA database. (c) Overall survival (OS) analysis of the high- and low-risk groups in patients from the E-MTAB-1980 cohort. (d) ROC curves demonstrated the accuracy of the risk score in the E-MTAB-1980 dataset. (e) Progression-free survival (PFS) in the high- and low-risk groups of patients from the CheckMate-025 and CheckMate-010 cohorts. (f) OS in the high- and low-risk groups of patients from the CheckMate-025 and CheckMate-010 cohorts. (g) Univariate (left) and multivariate (right) analysis of the training set from the TCGA cohort. (h) Univariate analysis (left) and multivariate (right) analysis of the test set from the E-MTAB-1980 cohort.

The PRGs with prognostic significance were associated with immune cell infiltration, immune checkpoints (Fig. 6a), and activation of the epithelial-mesenchymal transition (EMT) (Fig. 6b). Multivariate cox regression analysis showed that GSDMB was the most significant prognostic factor for ccRCC compared with other pyroptosis-related genes (Fig. 6c). Moreover, GSDMB was overexpressed in ccRCC tissue in GSE40435 and GSE53757 from the GEO dataset (Supplementary Fig. 9). Kaplan-Meier analysis showed that GSDMB overexpression was associated with poor OS of ccRCC patients (Supplementary Fig. 10). The TCGA database showed that GSDMB overexpression was significantly correlated to the stage and grade of ccRCC (Fig. 6d–f). qPCR showed that GSDMB expression at the transcription level was associated with ccRCC grade (Supplementary Fig. 11). Differential expression of GSDMD at the protein level was found between ccRCC cells and normal kidney cells, confirming previous findings (Fig. 6g,h). Furthermore, IHC revealed that GSDMB expression was related to ccRCC grade (Fig. 6i,j). Therefore, the above results indicate that GSDMB may promote the occurrence and development of renal clear cell carcinoma, although the precise mechanisms remain to be determined.

Fig. 6.

Comprehensive analysis of clinical features associated with the two clusters. (a) Relationship between 8 prognostic genes and infiltrating immune cells in the microenvironment. (b) Relationship between expression of the 8 prognostic genes, angiogenesis-related genes, and epithelial-mesenchymal transition (EMT)-related genes. (c) Multivariate cox regression analysis of overall survival was performed for 8 pyroptosis-related genes. (d–f) Relationship between GSDMB and clinical features of clear cell renal cell carcinoma (*p 0.05; ***p 0.001; ****p 0.0001). (g,h) Relative level of GSDMB protein expression in the HK-2, Caki-1 and 769P cell lines (*p 0.05; **p 0.01). (i) Representative IHC staining for GSDMB in normal, low and high grade ccRCC. (j) GSDMB expression level in normal (n = 15), low grade (n = 15) and high grade (n = 5) ccRCC tissues (**p 0.01; ****p 0.0001).