We collected key regulatory genes for 13 programmed cell death patterns as PCD-related genes from the GSEA gene set, KEGG, Hallmark, and review articles, and the final gene list is the tandem regulatory genes for 13 PCD patterns (4, 7–15). A final collection of 1217 PCD-related genes was included in the analysis ( Table S1 ). We downloaded transcriptome profiles, corresponding clinical information and mutation data of ccRCC samples from The Cancer Genome Atlas (TCGA-KIRC) (16) species. In addition, we obtained DNA microcohort data and clinical characterization of ccRCC as a validation cohort from ICGC and GEO databases (ID: GSE36895 and GSE53757).

We used univariate Cox regression to assess the association of PCD-related genes with survival status in ccRCC patients, adjusting the cutoff P-value to 0.001. candidate genes were further narrowed down by Lasso and multivariate Cox regression methods to construct the most suitable signatures.The programmed cell death index (PCDI) was finally obtained for each patient by the following formula: programmed Cell death index (PCDI) = Coef(Gene 1) × Expr(Gene 1) + Coef(Gene 2) × Expr(Gene 2) +…… + Coef(Gene n) × Expr(Gene n). where Coef(Gene) represents the risk regression coefficient of Gene and Expr(Gene) represents the expression of Gene. Based on the median value of PCDI, we categorized ccRCC patients into high and low PCDI groups. We further searched the Human Protein Atlas database (https://www.proteinatlas.org/) (17) to obtain YBX3, BID, TCIRG1, IDUA, KDR, PYCARD, IFNG, and LRRK2 between clear cell renal cell carcinoma tissues and normal kidney tissues at the protein level for Histological validation.

Based on PCD-related genes, we implemented consensus clustering using the R package “ConsensusClusterPlus” to identify ccRCC subtypes (18).

We obtained reference genomes (Hallmark, c5go and c2kegg) from the Molecular Signatures Database (MSigDB). Hallmark, c5go and c2kegg analyses were performed using the R package “clusterProfiler”. Screening conditions were |NES| > 1, NOM p-value < 0.05.

To determine the mutational load of ccRCC, we used the R package “TCGAbiolinks” to obtain ccRCC mutation data (19). We then used the “maftools” package to analyze the mutation data and obtain the TMB for each patient. the results were visualized using a waterfall Figure (20). Based on the median, TMB was categorized into two groups, high TMB and low TMB, and the mutation characteristics were compared between different risk groups. In addition, KM curves were utilized to detect the relationship between TMB, different risk classes of TMB and prognosis of ccRCC.

We used the R package “ESTIMATE” to perform tumor microenvironmental analysis to calculate tumor purity, ESTIMATE score, immune cell score and stroma score for each sample (21). The proportion of 22 types of immune infiltrating cells was calculated by the CIBERSORT algorithm (https://cibersort.stanford.edu/) (22). The single sample gene set enrichment analysis (ssGSEA) algorithm was used to quantify the relative proportions of immune cells and immune function.

We derived IC50 by ProPhetic algorithm to assess drug response to common chemotherapeutic treatments for ccRCC, comparing drug sensitivity to the above chemotherapeutic treatments in patients with high PCDI and low PCDI. We downloaded the gene expression data of cancer cells to different drugs at the Tumor Drug Sensitivity Multi-Organics (GDSC) database (https://www.cancerrxgene.org/) (23) and calculated IC50 values to assess patients’ responses to chemotherapeutic drugs. The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm can be used to infer patient response to immunotherapy (24). In addition we downloaded anti-PD-1 and anti-CTLA4 IPS scoring data from ccRCC via the TCIA database (https://tcia.at/home) (25) to assess patient response to immune checkpoint inhibitors.The IMvigor210 study cohort evaluated the efficacy of atezolizumab (a PD-L1 targeting antibody) in 210 patients with locally advanced or metastatic uroepithelial cancer (26).The Liu study cohort evaluated 121 patients with metastatic melanoma treated with an anti-PD-1 inhibitor and the Kim study cohort of patients with metastatic gastric cancer treated with pembrolizumab (an anti-PD-1 inhibitor) (27). In addition, we collected genes that have been currently reported to be positively associated with immuno-efficacy and negatively associated with immuno-efficacy and analyzed them in association with PCDI (28–30).

Differential analysis of ccRCC and surrounding normal tissue was performed using the R package “limma” with a cutoff value set to log2 fold change (logFC) >1 and adjusted false discovery rate (FDR) <0.001 (31). Heatmap visualization was performed using the R package “pheatmap”. The R packages “rms” and “regplot” were used for the plotting of column line graphs and calibration curves. We calculated the tumor mutational load (TMB) of each patient and compared it between the high-risk and low-risk groups, and then applied the R package “Maftools” to explore the somatic mutation information in ccRCC patients. Principal component analysis (PCA) was performed using the “stats” package.

ACHN and 786−O cells were purchased from The Cell Bank of Type Culture Collection of The Chinese Academy of Sciences (Shanghai, China) and were cultured in accordance with the cell culture manual of the manufacturer. Passage time was less than 6 months and mycoplasma were eliminated according to the test report. ACHN cells were cultured in MEM medium (HyClone, Logan, UT, USA) and 786-O cells were cultured in RPMI-1640 medium (HyClone, Logan, UT, USA), the culture media contained 10% fetal bovine serum (FBS, HyClone) according to cell culture manual. All RCC cells were cultured in a humidified atmosphere of 5% CO₂ at 37°C. When the cell fusion rate was 90%, the cells were digested with 1-ml trypsin for 5 minutes, then adding 1 ml of medium into the solution for neutralization. Collected by centrifugation, the cells were resuspended with culture media and passaged.

ACHN and 786-O cell lines were cultured in 6-well plates for the transfection of the siRNA. When the cells confluence rate was 60%, the transfection complexes with IDUA siRNA (JTS scientific, Wuhan, China) were added into the culture solution by means of Lipofectamine™3000 (Invitrogen, USA) on the basis of the procedure in the instruction manual. The siRNA sequences were as follows (5′-3′): siIDUA-1# (sense ACUUUGAGGACAAGCAGCAGGUGUU and anti-sense AACACCUGCUGCUUGUCCUCAAAGU), siIDUA-2# (sense CAGCAGGUGUUUGAGUGGAAGGACU and anti-sense AGUCCUUCCACUCAAACACCUGCUG).

Total cell lysates were extracted by adding the RIPA solution (RIPA: PMSF, 100:1; P0013B; Beyotime Institute of Biotechnology) containing protease and phosphatase inhibitors. The concentration of total cell protein was assessed with Enhanced BCA Protein Assay Kit (cat. no. P0010; Beyotime Institute of Biotechnology). Denatured protein (40 μg/lane) realized the protein separation on a 10% sodium dodecyl sulfate polyacrylamide gel electrophoresis (140 V, 50 minutes) gel, then target protein was transferred to the polyvinylidene fluoride membrane (350 mA, 90 minutes) by transfer. The membrane was soaked in a sealed container using 5% fat-free milk at 37°C for 1 hour to realize the protein blockade. The membrane was incubated overnight with the primary antibody in 5% fat-free milk at 4°C. The information of the primary antibodies in the WB were as follows: anti-IDUA (1:1000, 55158-1-AP, Proteintech) and anti-GAPDH (1:5000, 10494-1-AP, Proteintech). The membrane was washed three times (5 minutes each) using tris-buffered saline tween-20, then the membrane was incubated with the corresponding secondary antibody at 37°C for 1 hour. EasySee Western Blot kit (Beijing Transgen Biotech, Beijing, China) was used for detection following the experimental operation manual of manufacturer. Density measurements were carried out using ImageJ software (National Institute of Health, Bethesda, MD, USA), and the density of target protein bands were normalized with GAPDH.

Cell viability was measured by colorimetric assay by means of Cell Counting Kit-8 (CCK-8) (Bimake, USA). 3 × 10³ cells were seeded in each well of 96-well plates. After 24 hours of transfection, CCK-8 working solution was added to the solution and incubated at 37°C for 1 hour. The absorbance of culture solution was measured at 450 nm using a plate reader (Model 680; Bio-Rad Laboratories) at 24, 48 and 72 h.

ACHN and 786-O cells were transfected with IDUA siRNAs in 24-well plates for 48 h, with EdU (BeyoClick™, EDU-555, China) added to the medium (1:1000) to assess the cell proliferation. Cells were incubated for 2 hours at 37°C, after labeling, the culture medium was removed, and 1-ml fixative solution (4% paraformaldehyde) was added for fix cells at room temperature for 20 minutes. Cells were incubated at room temperature for 15 minutes with 1-ml permeate (0.3% Triton X-100) and the cells were treated with click reaction buffer at 37°C for 30 min, after which Hoechst solution was added for 15 min at room temperature. A fluorescence microscope (Olympus Corporation, Japan) was used to obtain the images to show the cell proliferation.

600 μL medium (10% FBS) was added to the bottom of the well, which is the lower chamber. 1 × 10⁵ cells were suspended in 200 μL serum-free medium and the suspension were seeded into the upper chamber (Corning, NY, USA). Then the chambers were placed into the 24-well plates. Cells were incubated at 37°C and 5% CO₂ for 24 hours. After removing the remaining cells in the upper chamber using the swab, 4% paraformaldehyde was used to fix the cells for 10 minutes, and crystal violet stain was added for staining lasting for 10 minutes at room temperature. Cell migration was observed by means of optical microscopy, and ImageJ software lent itself to calculate the migration efficiency.

When the confluence rate of cells grew to 100% in the 6-well plates, a scratch was made in the center by means of the tip of the 1 mL pipette. Then cells were cultured with serum-free medium at 37˚C for 48 h. Cell images were gained with an microscope (EVOS XL system; AMEX12000; Thermo Fisher Scientific, Inc.). ImageJ software were applied to calculate the wound area.

Survival curves were plotted using the Kaplan-Meier method to compare the difference in survival between the two groups. p-value ≤ 0.05 was considered statistically significant. All statistical analyses were performed by R.

We performed consensus clustering analysis to explore ccRCC subtypes by collecting PCD-related genes (1217 genes in total). We found that differences between subgroups were most significant when k=2, indicating that ccRCC patients could be well categorized into two clusters ( Figure 1A ). The results of principal component analysis (PCA) showed that patients with renal clear cell carcinoma could be categorized into two clusters based on PCD-related genes ( Figure 1B ). We found a significant difference in patient overall survival (OS) between the two clusters (p < 0.001), with C1 associated with patients with a good prognosis and C2 associated with patients with a poor prognosis ( Figure 1C ). By comparing the differences in clinicopathologic characteristics between the two clusters, including ETHNICITY, T, N, M, and whether lymph node metastasis occurred. We found that not Hispanic or latino, T, N, M staging and occurrence of lymph node metastasis were higher in C2 than in C1, suggesting a higher malignant phenotype in C2 patients ( Figure 1D ).

We then analyzed the immune profile between the two clusters. We found that the tumor purity of C2 was significantly lower than that of C1, while the ESTIMATEScore, immune score and stromal score of C2 were significantly higher than that of C1, suggesting a high level of immune infiltration in the C2 phenotype compared to C1 ( Figure 1E ). Analysis of the differences in the expression of common immune checkpoints between the two clusters showed that most of the immune checkpoints were expressed significantly higher in type C2 than in type C1 ( Figure 1F ). The infiltration abundance of immune cells was assessed by the CIBERSORT algorithm, and the results showed that the level of immune cell infiltration was higher in type C2 ( Figure 1G ). Finally, the effect of the two clusters on the response to immunotherapy was assessed by the IPS score, and we found that C2 was significantly more effective than C1 against anti-PD-1 and/or anti-CTLA-4 inhibitors ( Figure 1H ).

We identified 8 genes and thus constructs with programmed cell death and its derived PCD by one-way Cox regression, lasso and multifactorial Cox regression analyses ( Figure 2A ). Our model derived the programmed cell death index (PCDI) for each patient by the following equation. PCDI= (0.0413*YBX3 exp.) + (0.1458*BID exp.) + (0.0341*TCIRG1 exp.) + (0.0943*IDUA exp.) + (-0.0192*KDR exp.) + (-0.0485*PYCARD exp.) + (0.1625*IFNG exp.) + (-0.0373*LRRK2 exp) ( Figure 2B ). We then explored the expression of PCDI genes in unpaired and paired samples of tumor and normal tissues in the TCGA-KIRC cohort, and showed that all eight genes comprising PCDI were significantly upregulated and expressed in tumor tissues (p < 0.05) ( Figures 2C, D ). We validated these results by two GEO datasets (GSE36895 and GSE53757), and the results were consistent ( Figure 2E, F ). The expression of the eight genes comprising PCDI at the protein level was analyzed by immunohistochemical staining of the HPA database, and we found that compared with normal tissues, the eight genes comprising PCDI stained more deeply immunohistochemically in tumor tissues, suggesting that the content of these eight proteins was higher than that in normal tissues in ccRCC ( Figure 2G ).

We found that patients in the high PCDI group had a poorer prognosis and were more likely to have a higher mortality rate through Kaplan-Meier survival curves (p < 0.05, Figure 3A ). Subsequently, we used the ICGC dataset as a validation cohort, and Kaplan-Meier analysis showed that patients in the high PCDI group were more likely to have poorer overall survival ( Figure 3B ). The results of univariate and multivariate Cox regression analyses performed showed that the p-value of PCDI was less than 0.05, suggesting that PCDI can be an independent prognostic factor in patients with ccRCC ( Figures 3C, D ). To further validate the clinical significance of PCDI, we analyzed the differences in PCDI between different clinical characteristics and showed that PCDI was higher in patients with C2, male, high grading, high T stage, M1 and N1 stage, suggesting that the higher our PCDI, the closer ccRCC was to advanced stage ( Figure 3E ). Analysis of survival curves stratified by clinical characteristics showed that PCDI significantly differentiated the prognosis of each clinical subgroup, with patients in the high PCDI group having a worse prognosis ( Figure S1 ). We evaluated the area under the curve (AUC) values of the TCGA cohort and compared them with other clinical traits, and the results showed that PCDI had high accuracy in predicting 1-, 3-, and 5-year survival in patients with ccRCC and was a better predictor of patient survival than other clinical traits ( Figures 3F, G ). With the results of univariate and multivariate Cox regression analyses, we modeled column-line plots for ccRCC patients to assess patient prognosis, with age, gender, Laterality, and PCDI included in the models ( Figure 3H ). The calibration curves of the column-line plots showed that the predicted 1-, 3-, and 5-year survival rates were more consistent with the actual survival rates at the reference line, suggesting that the constructed column-line plots could predict patient prognosis well ( Figure 3I ).

To explore the generalization of PCDI to other cancers, we used the modeling equations for PCDI described above to calculate PCDI values for patients with the other 32 cancer types in TCGA and to plot Kaplan-Meier survival curves for the high/low PCDI group. For Overall Survival (OS), patients in the high PCDI group in ACC, UVM, LAML, THYM, and TGCT had a poorer prognosis, whereas patients in the low PCDI group in BLCA, BRCA, ESCA, STAD, and SKCM had a poorer prognosis ( Figure 4A ). For Disease Specific Survival (DSS), patients in the high PCDI group in KIRP, TGCT, and THYM had poor Disease Specific Survival, while patients in the low PCDI group in BLCA, SKCM, and STAD had poor Disease Specific Survival ( Figure 4B ). For Disease Free Interval (DFI), patients in the high PCDI group in KIRP had poor Disease Free Interval, whereas patients in the low PCDI group in COAD had poor Disease Free Interval ( Figure 4C ). For Progression Free Interval (PFI), patients in the high PCDI group in KIRP, PAAD, and THYM had poor Progression Free Interval, whereas patients in the low PCDI group in BLCA and STAD had poor Progression Free Interval ( Figure 4D ). The above results suggest that PCDI not only has a good effect in predicting the prognosis of ccRCC (KIRC), but also has a predictive prognostic value in other cancers.

In order to investigate the regulatory pathways of tumorigenesis in the high PCDI group, we performed GSEA analysis of tumors in the high PCDI group, which showed that the high PCDI group was significantly enriched in CELL Killing, CYTOLYSIS, Necrotic cell death, Pyroptosis, and Apoptosis (p < 0.05) ( Figure 5A ). In addition, tumors in the high PCDI group were significantly enriched (p < 0.05) in immune function pathways such as B cell mediated immunity, NK cell activation, T cell activation and T cell mediated immunity, suggesting that our high PCDI group and tumor immune microenvironment closely associated with the tumor immune microenvironment ( Figure 5B ).

Using the ESTIMATE algorithm, we found that the high PCDI group had lower tumor purity but higher immune scores and stroma scores ( Figure 5C ). ssGSEA algorithm results showed that patients in the high PCDI group had better immune cell infiltration and immune-related functions compared to the low PCDI group ( Figure 5D ). cIBERSORT algorithm showed that the high PCDI group The level of immunostimulatory CD8 T cells was significantly higher than that of the low PCDI group, while the level of immunosuppressive M2-type macrophages was significantly lower in the high PCDI group, which suggested that our high PCDI group had high immune infiltration characteristics ( Figure 5E ).

To further explore the relationship between PCDI and the immune microenvironment, we compared the differences in the expression levels of common immune checkpoints, MHC molecules, cytokines, and cytokine receptors between the high/low PCDI groups, and the results showed that the expression levels of immune checkpoints, MHC molecules, cytokines, and cytokine receptors of patients in the high PCDI group were significantly higher than those in the low PCDI group (p < 0.05, Figure 6A ). Then we analyzed the correlation between PCDI and common immune checkpoints, and the results showed that most of the immune checkpoints were significantly positively correlated with PCDI (p < 0.05, Figure 6B ). These results suggest that the tumors in the high PCDI group showed “hot” tumor characteristics, and the efficacy of immune checkpoint inhibitor therapy in the high PCDI group may be better than that in the low PCDI group. By analyzing the association between TMB, MSI and PCDI, we found that patients in the high PCDI group showed higher TMB and MSI ( Figures 6C, D ). In addition, the prognosis of patients in the high TMB group was worse than that of patients with low TMB. ccRCC patients were divided into four groups after the combination of TMB and PCDI, and we found that the prognosis of patients in the high TMB+high PCDI group, the high TMB+low PCDI group, and the low TMB+high PCDI group were all similarly worse than that of patients in the high TMB+high PCDI group ( Figures 6E, F ). We also examined the mutation rates in ccRCC in the cBioPortal database for the eight genes from which PCDI was constructed and found a higher mutation rate of 15% for LRRK2 and a lower mutation rate of 4% for PYCARD compared to the other six genes from which PCDI was constructed ( Figure S2 ).

Subsequently, we collected external immunotherapy datasets to validate PCDI to predict immunotherapy efficacy, and the results showed that PCDI was significantly higher in patients in the immunotherapy-responsive group than in the non-responsive group in Imvigor210 (anti-PD-L1), Kim cohort (anti-PD-1), and Liu cohort (anti-PD-1) patients ( Figures 6G–I ). We then explored the association between PCDI group and immunophenotypic core (IPS) by assessing immunogenicity to predict patient response to immune checkpoint blockade (anti-PD1 and/or anti-CTLA4), and we found that IPS scores were higher in the high PCDI group, suggesting that patients in the high PCDI group may be better responders to immunotherapy ( Figure 6J ). We also compared the association between positive and negative immune efficacy-related genes and PCDI, and showed that PCDI was significantly positively associated with positive immune efficacy-related genes (POLD1 and POLE) and significantly negatively associated with negative immune efficacy-related genes (CTNNB1, EGFR, JAK1 and PTEN) ( Figure 6K ). All these results indicated that PCDI could better predict the immunotherapy effect, and the high PCDI group had better effect on immunotherapy.

To explore the association between PCDI and drug sensitivity, we calculated half-maximal inhibitory concentration (IC50) values for common drugs in ccRCC samples and compared them with those between PCDI subgroups. We found that many common chemotherapeutic agents for renal cell carcinoma were significantly different between the high/low PCDI groups (p < 0.05), and the IC50 values for Docetaxel were significantly lower in the high PCDI group than in the low PCDI group, suggesting that patients in the high PCDI group may have a better response to Docetaxel-based chemotherapy.The IC50 values for Lapatinib were significantly lower in the low PCDI group was significantly lower than that in the high PCDI group, suggesting that patients in the low PCDI group had a better response to Lapatinib-based chemotherapy ( Figure 7A ). Then, we analyzed the relationship between drug sensitivity and mRNA expression of the eight genes constructing PCDI using the GDSC and CTRP databases; positive correlation indicated that gene expression was associated with drug resistance, while negative correlation indicated that gene expression was associated with drug sensitivity. The results indicated that BID and YBX3 gene expression was associated with most drug sensitivities ( Figures 7B, C ).

To confirm the biological function of IDUA in RCC, we knocked down IDUA using two siRNAs in ACHN and 786-O cells. Western blotting showed that IDUA could be effectively silenced by two independent siRNAs ( Figure S3 ). The proliferation of ACHN and 786-O cells was noticeably decreased after down-regulation of IDUA expression compared with the corresponding control cells in the CCK8 assay ( Figures 8A, B ). In addition, EDU incorporation analysis also showed that the proportion of EDU-positive ACHN and 786-O cells was significantly reduced in the IDUA-interfered group compared with the corresponding control cells ( Figures 8C, D ). The results of transwell assay showed that interference with IDUA crippled the migratory capacity of ACHN and 786-O cells ( Figures 8E, F ). Subsequently, in the wound-healing assay, we found that IDUA knocking-down decelerated wound closure speed in both ACHN and 786-O cells ( Figures 8G, H ). These results indicated that IDUA could promote renal cell carcinoma cell proliferation and migration in vitro. IDUA is expected to be a potential target for the treatment of renal cell carcinoma.