A novel risk model of three SUMOylation genes based on RNA expression for potential prognosis and treatment sensitivity prediction in kidney cancer

Introduction: Kidney cancer is one of the most common and lethal urological malignancies. Discovering a biomarker that can predict prognosis and potential drug treatment sensitivity is necessary for managing patients with kidney cancer. SUMOylation is a type of posttranslational modification that could impact many tumor-related pathways through the mediation of SUMOylation substrates. In addition, enzymes that participate in the process of SUMOylation can also influence tumorigenesis and development.

Methods: We analyzed the clinical and molecular data which were obtanied from three databases, The Cancer Genome Atlas (TCGA), the National Cancer Institute’s Clinical Proteomic Tumor Analysis Consortium (CPTAC), and ArrayExpress.

Results: Through analysis of differentially expressed RNA based on the total TCGA-KIRC cohort, it was found that 29 SUMOylation genes were abnormally expressed, of which 17 genes were upregulated and 12 genes were downregulated in kidney cancer tissues. A SUMOylation risk model was built based on the discovery TCGA cohort and then validated successfully in the validation TCGA cohort, total TCGA cohort, CPTAC cohort, and E-TMAB-1980 cohort. Furthermore, the SUMOylation risk score was analyzed as an independent risk factor in all five cohorts, and a nomogram was constructed. Tumor tissues in different SUMOylation risk groups showed different immune statuses and varying sensitivity to the targeted drug treatment.

Discussion: In conclusion, we examined the RNA expression status of SUMOylation genes in kidney cancer tissues and developed and validated a prognostic model for predicting kidney cancer outcomes using three databases and five cohorts. Furthermore, the SUMOylation model can serve as a biomarker for selecting appropriate therapeutic drugs for kidney cancer patients based on their RNA expression.

Introduction

Kidney cancer is one of the most common and lethal urological malignancies. Clear cell renal cell carcinoma (ccRCC) composes approximately 70%–85% of all kidney cancers and arises from the epithelium of the proximal tube of the kidney (Li et al., 2019a). Patients with ccRCC have a poor outcome, although targeted therapy and immune checkpoint blockade (ICB) agents have been recommended as first-line treatment for metastatic ccRCC (Albiges et al., 2019). More accurate prognosis prediction and personalized treatment plans are essential for the full-course management of patients with ccRCC. Today, many risk models that consist of clinical and molecular features have been proposed to predict the outcome in patients with kidney cancer (Wei et al., 2019; Ning et al., 2021; Li et al., 2022).

Posttranslational modification (PTM) is a process to modify polypeptide chains through chemical modifications, which makes protein functions more diverse. About 400 types of PTMs have been reported. SUMOylation is carried out through small ubiquitin-related modifier (SUMO) proteins and enzymes that connect SUMO proteins to a lysine residue in the target protein (Ramazi and Zahiri, 2021). SUMOylation has multiple physiological functions, such as chromatin remodeling, stemness, cell identity, and protein stability. Abnormal SUMOylation has been proven to be associated with cancer occurrence, neurodegeneration, and infection (Celen and Sahin, 2020). In cancers, SUMOylation of promyelocytic leukemia/retinoic acid receptor alpha (PML/RARA) is important for cellular transformation and essential for the pathogenesis of acute promyelocytic leukemia. Androgen receptor SUMOylation is associated with tumorigenesis of prostate cancer, while SUMOylation of proto-oncogene MYC occurs in B cell lymphoma (Sahin et al., 2022). SUMOylation of tumor suppressor gene p53 might have two opposite functions in cancer, cancer promotion, and cancer suppression (Bettermann et al., 2012). Some drugs that target SUMOylation have been discovered. Ginkgolic acid could reduce cellular SUMOylation and inhibit the invasion capacity of lung, colon, and liver cancers (Baek et al., 2017; Qiao et al., 2017; Li et al., 2019b; Sahin et al., 2022).

In kidney cancer, transcription factors, hypoxia-inducible factor 1A (HIF1A) and hypoxia-inducible factor 2A (HIF2A) are important for tumorigenesis and development of ccRCC through activating transcription of genes, such as vascular endothelial growth factor, which could impact the angiogenesis signaling pathway (Schodel et al., 2016). Previous studies have shown that SUMOylation could also take part in the carcinogenesis of kidney cancer. SUMOylation of HIF1A, E3 ligase hypoxia-associated factor could bind to HIF2A and enhance its transcriptional activity, which then promotes the metastasis of ccRCC cells (Koh et al., 2015). In addition, mutation of microphthalmia-associated transcription factor (MITF) could impair the SUMOylation of MITF, and then improve the transcriptional activity of HIF1A through binding to its promoter (Bertolotto et al., 2011). A study has also shown that SUMOylation-related genes could help predict the outcome of patients with bladder cancer (Guo et al., 2022). However, the role of SUMOylation-related genes have not been well discussed in kidney cancer.

In this study, we screened the differentially expressed SUMOylation genes at the RNA level between kidney cancer and normal kidney tissues. A SUMOylation risk model was constructed to predict prognosis and guide the drug selection for the treatment of patients with kidney cancer. We then found the different immune statuses between various SUMOylation risk groups and analyzed protein expression further. This study not only built a robust prognosis model but also screened the significant SUMOylation-related genes, which may be essential for the development of kidney cancer.

Methods

Data acquiring

Three databases were used in this study: The Cancer Genome Atlas (TCGA), the National Cancer Institute’s Clinical Proteomic Tumor Analysis Consortium (CPTAC), and ArrayExpress. Specifically, the Kidney Renal Clear Cell Carcinoma (KIRC) cohort in the TCGA database, PDC000127 cohort in the CPTAC database, and E-MTAB-1980 cohort in the ArrayExpress database were enrolled in our study (Sato et al., 2013; Clark et al., 2019). Clinical and RNA expression data of these three cohorts, and protein expression data of the PDC000127 cohort and Chinese Fudan cohort were downloaded and processed (Clark et al., 2019; Qu et al., 2022). Data of patients who had complete clinical data in the three cohorts were used for prognosis analysis. Detailed patient lists and clinical features are shown in the supplementary file and summarized in Supplementary Table S1. The SUMOylation genes were obtained from the molecular signature database (Supplementary Material).

Analysis of differential SUMOylation gene expression

RNA-seq data from 72 normal kidney tissues and 539 tumor tissues in the TCGA-KIRC cohort were used to screen the differential SUMOylation gene expression. The “limma” package was used, and genes with a fold change>1 and false discovery rate <0.05 were considered differentially expressed (Ritchie et al., 2015). The intersection of differential SUMOylation genes, and the TCGA-KIRC, CPTAC, and E-MTAB-1980 gene lists were used for further clinical prognosis analysis.

SUMOylation-related gene risk model construction

In total, 513 TCGA-KIRC patients, 89 CPTAC patients, and 99 E-MTAB patients were diagnosed with ccRCC by pathology, and with complete clinical data (age, sex, clinical stage, tumor stage, metastasis stage, node stage, grade, and survival time more than 1 month), were analyzed for building a prognosis risk model. First, the 513 TCGA-KIRC patients were randomly separated into a discovery TCGA cohort (256 patients) and a validation TCGA cohort (257 patients) in a 1:1 ratio, using the set.seed function of Base R, and the seed was set as 4443. The discovery TCGA cohort was used for constructing a risk model, and the differential SUMOylation genes were analyzed by univariate cox regression to find the prognosis-related genes. These genes were analyzed with Lasso regression analysis, and the selected genes were tested in multivariate cox regression analysis using the Akaike information criterion method. The coefficient and expression value of these genes were employed to compute a risk score and construct a risk model, which was named the SUMOylation risk model in this study.

The risk score of each patient in the five cohorts (discovery TCGA cohort, validation TCGA cohort, total TCGA cohort, CPTAC cohort, and E-MTAB-1980 cohort) was calculated. Patients in each cohort were allocated into high- and low-risk groups according to the median SUMOylation risk score as the cut-off value. Outcomes of the different groups of patients were evaluated by log-rank analysis and presented by Kaplan-Meier curve, and the receiver operating characteristic curve was conducted. The propensity-matched analysis (PSM), subgroup analysis, and interaction test was conducted to further test and validate the role of the risk group in outcome predicting. The PSM test was conducted using the “MatchIt” package and the caliper was set as 0.03. The role of SUMOylation risk score in prognosis prediction was investigated by univariate and multivariate cox regression sequentially. In addition, a nomogram was built, and then multi-ROC and calibration analyses were conducted to test the nomogram.

The correlation between the SUMOylation risk group and molecular features

Gene Set Variation Analysis (GSVA) analysis which was conducted by “c2. cp.kegg.v6.2. symbols.gmt” and “GSVA” package was used to clarify the molecular features in the tumor tissues of different risk groups. The infiltration immune cell proportion in tumor tissues of TCGA-KIRC patients was computed and assessed using CIBERSORT (https://cibersort.stanford.edu/) (Newman et al., 2015).

In addition, the immune subtype was raised according to the 160 immune expression signatures, and the all-type tumors in the TCGA database were clustered into six subtypes, C1 (wound healing), C2 (IFN-g dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-b dominant) (Thorsson et al., 2018). Subtype data of the TCGA-KIRC tumors were taken out.

A previous study proposed and computed the TIDE score of samples across tumor types in the TCGA database based on the expression of specific genes. The results are published at https://tide.dfci.harvard.edu/, and related information of TCGA-KIRC patients was extracted from there (Jiang et al., 2018).

Drug treatment sensitivity predicting information was packaged in the “pRRophetic” package (Yang et al., 2013; Geeleher et al., 2014). The targeted agent’s treatment sensitivity of TCGA-KIRC patients was calculated using the package and represented half maximal inhibitory concentration (IC50).

The SUMOylation risk model genes encoding protein expression data of CPTAC and Chinese Fudan cohort were obtained from the published data (Clark et al., 2019; Qu et al., 2022). The IHC data and images of renal cancer and kidney tissues were obtained and analyzed using the HPAnalyze package (Tran et al., 2019).

Immunohistochemistry (IHC) detection and IHC score

The tissue microarray (TMA) containing 30 kidney cancer tissues and 30 adjacent normal kidney tissues was purchased from SHANGHAI OUTDO BIOTECH CO., LTD (Lot: HKid-CRCC060PG-01). Primary antibodies targeting CDCA8 (PROTEINTECH: 12465-1-AP), CDH1 (PROTEINTECH: 20874-1-AP), and PPARA (PROTEINTECH: 66826-1-Ig) were used to detect the corresponding protein expression via IHC. The IHC procedure involved deparaffinization and rehydration of the TMA, followed by heating the slides in Tris-EDTA (PH 8.0) for antigen retrieval. Subsequently, the slides were incubated in 0.3% H₂O₂ for 30 min and blocked in 10% goat serum to prevent non-specific binding. The slides were then incubated with antibodies against CDCA8 at a 1:100 dilution, CDH1 at a 1:600 dilution, and PPARA at a 1:100 dilution overnight at 4°C. Finally, the slides were incubated with the secondary antibody and visualized using DAB chromogen.

We used ImageJ software to analyze the staining intensity of IHC. We measured the staining intensity for each tissue section and divided the samples into quartiles. Based on these quartiles, we classified the samples into three groups: low expression (staining intensity below the 25th percentile), moderate expression (staining intensity between the 50th and 75th percentiles), and high expression (staining intensity above the 75th percentile).

Statistical analysis

Continuous variables were analyzed by the t-test, while categorical variables were examined using the chi-square test. All statistical analyses were conducted using R software, and a two-sided p-value of < 0.05 was considered statistically significant.

Results

Differential SUMOylation genes in TCGA-KIRC dataset

Among these SUMOylation genes, 29 genes were abnormally expressed in kidney cancer tissues. Of these, 17 genes (CDKN2A, H4C5, AURKB, BIRC5, H4C11, TOP2A, cell division cycle-related protein 8 [CDCA8], H4C12, NR5A2, BLM, NFKB2, PML, H4C8, H4C9, RARA, NR1H3, and H4-16) were upregulated and 12 genes (CHD3, peroxisome proliferator-activated receptor alpha [PPARA], NR4A2, CDH1, VDR, PPARGC1A, THRB, NR3C2, TFAP2A, TFAP2C, NOS1, and TFAP2B) were downregulated than in the adjacent normal kidney tissues (Supplementary Material; Supplementary Figure S1).

SUMOylation gene risk model for predicting the survival of patients with kidney cancer

In these 29 genes, seven (AURKB, BIRC5, CDCA8, CDH1, NR1H3, PPARA, and TOP2A) were correlated with the prognosis of kidney cancer through univariate cox regression analysis using the discovery TCGA cohort, and five genes (AURKB, CDCA8, CDH1, NR1H3, and PPARA) were selected from Lasso regression (Supplementary Figure S2). In addition, three genes (CDCA8, CDH1, and PPARA) presented as independent factors of kidney cancer by multivariate cox regression (Supplementary Material). Finally, a SUMOylation risk model was constructed with the following equation: risk score = 0.413023054*CDCA8- 0.024059083*CDH1-0.12352751*PPARA. The patients’ risk scores in these five cohorts were computed, and patients in each cohort were divided into two risk groups. In the discover TCGA cohort, patients with high-risk scores often had a short survival time (Figure 1A), and log-rank analysis also showed that patients in the high-risk group had a significantly worse outcome (p < 0.001, Figure 2A). These findings were verified by the four other cohorts (validation TCGA cohort, total TCGA cohort, CPTAC cohort, and E-MTAB-1980 cohort) (Figures 1B–E; Figures 2B–E).

FIGURE 1

FIGURE 1. The survival and SUMOylation risk score distribution and SUMOylation risk model gene expression status. In each sub-figure, the upper panel shows the risk score distribution, the dashed horizontal line shows the median value of the risk score, and the dashed vertical line shows the patients with the median value of the risk score (green dots show the low-risk group, and red dots shows the high-risk group). The middle panel shows the survival status (red dots represent dead status, and green dots represent live status), and the bottom panel shows the heatmap of the SUMOylation risk model gene expression in each cohort. (A) Discovery TCGA cohort. (B) Validation TCGA cohort. (C) Total TCGA cohort. (D) CPTAC cohort. (E) E-MTAB-1980 cohort.

FIGURE 2

FIGURE 2. The different outcomes in various SUMOylation risk groups. In each figure, the red line represents the outcome of patients in the high SUMOylation risk group, the blue line shows the outcome of patients in the low SUMOylation risk group, and the box shows the live patient count during the entire follow-up period. Log-rank analysis was used to compare the outcomes of patients in different risk groups and presented by the Kaplan-Meier curve. (A) Discovery TCGA cohort. (B) Validation TCGA cohort. (C) Total TCGA cohort. (D) CPTAC cohort. (E) E-MTAB-1980 cohort. (F) Total TCGA cohort (After PSM). Propensity-matched analysis: PSM.

In the total TCGA cohort, the high SUMOylation risk group is significantly correlated with Male patients (p = 0.003), high tumor stage (p < 0.001), metastasis stage (p < 0.001), clinical stage (p < 0.001), and tumor grade (p < 0.001) (Table 1). After the PSM analysis, each group has been assigned 180 patients, and there is no significant difference in clinical features between the high and low-risk groups (Table 1; Figure 3A). And the high-risk group still presents a poorer prognosis than the low-risk group (Figure 2F).

TABLE 1

TABLE 1. The clinical feature in different risk group in total TCGA and total TCGA cohort after PSM.

FIGURE 3

FIGURE 3. The propensity-matched analysis (PSM), subgroup analysis, and interaction test of SUMOylation risk group in the total TCGA cohort. (A) The red line showed the SMD value of the original total TCGA cohort. The green lines showed the SMD value of each variable after the adjustment of the PSM analysis. (B) The forest figure showed the result of the subgroup analysis and interaction test, and low-risk group as the reference.

In addition, the subgroup and interaction analysis also show high-risk group associated with poor outcomes in the patients with the following features, ≤65 years, male, clinical stage I, clinical stage IV, grade III, grade IV, T1, T2, M0+Mx, M1 and N0+Nx. The contrary result in the N1 subgroup might be due to the limited number of patients (Figure 3B).

Although clinical features present in each cohort may vary in functions for predicting the survival status, the SUMOylation risk score could be an independent prognosis factor in all five cohorts and presented a good prognostic effect (Table 2; Figure 4).

TABLE 2

TABLE 2. Univariate and Multivariate Cox regression reveals riskscore is an independent risk factor in kidney cancer.

FIGURE 4

FIGURE 4. The SUMOylation risk score is an independent poor prognosis factor of kidney cancer patients. In each subfigure, the left panel shows the multivariate cox regression results, and the right panel presents the receiver operating characteristic curve of 1, 3, and 5 years (the 1 and 3 years of the CPTAC cohort). (A) Discovery TCGA cohort. (B) Validation TCGA cohort. (C) Total TCGA cohort. (D) CPTAC cohort. (E) E-MTAB-1980 cohort.

In the total TCGA cohort, age, clinical stage, Grade, and SUMOylation risk score presented as independent risk factors, and these four variables were used to build a nomogram to predict the survival rate of patients with kidney cancer (Figure 5).

FIGURE 5

FIGURE 5. A nomogram for predicting the survival rate of kidney cancer patients. (A) Using multivariate cox regression results of the total The Cancer Genome Atlas cohort, patient age, Grade, clinical stage, and SUMOylation risk score show independent risk factors of patient outcomes. Each feature was assigned a value, and the total value of patients corresponded to values of 1-, 3- and 5-year survival rates. (B) The nomogram shows the highest AUC value of in multi-ROC analysis. (C) Calibration curves of the nomogram for predicting the outcome of 1, 3, and 5 years in the total TCGA cohort.

Different SUMOylation subgroups present various molecular features

GSVA analysis revealed the proteasome, DNA replication, glycosaminoglycan biosynthesis chondroitin sulfate, glycosaminoglycan biosynthesis keratan sulfate pathways were significantly highly enriched in kidney cancer tissues with the high-risk group (Figure 6A).

FIGURE 6

FIGURE 6. The cell signaling and immune status in the different risk groups of kidney cancer tissues. (A) The GSVA analysis shows seven pathways were up-regulated in the low-risk group, while four pathways were up-regulated in the high-risk group. (B) The different types of immune cells infiltrated in the tumor tissues, *** indicates p < 0.001, ** indicates p < 0.01, and * indicates p < 0.05. (C) The patients’ immune subtypes in the different risk groups. C1 (wound healing), C2 (IFN-g dominant), C3 (inflammatory), C4 (lymphocyte depleted), and C6 (TGF-b dominant).

In the two groups (the high-risk and low-risk groups), 10 types of immune cell infiltration statuses were different. Among them, T cell CD8, T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs), T cells gamma delta, and macrophages M0, were higher in the tissues of the high-risk group, while T cells CD4 memory resting, monocytes, macrophages M2, and mast cells resting were higher in the low-risk group (Figure 6B). Tissues in the high-risk group presented higher ratios of C1, C2, C4, and C6 (total 20%) than tissues in the low-risk group (total 6%). Additionally, they had a lower C3 ratio (80%) than tissues in the low-risk group (94%) (Figure 6C).

SUMOylation subgroups may be a biomarker of kidney cancer treatment

Tissues in the high-risk group showed a high TIDE score and a lower MSI score (Figure 7A). In addition, tumor tissues in the high-risk group had a lower IC50 of axitinib (p = 0.047), sorafenib (p < 0.001), sunitinib (p < 0.001), pazopanib (p < 0.01), and temsirolimus (p < 0.001) (Figure 7B).

FIGURE 7

FIGURE 7. The TIDE score and targeted drug treatment sensitivity prediction in different risk groups. (A–C) TIDE, MSI, and exclusion scores of the two SUMOylation risk groups (*** indicates p < 0.001, ns indicates no significant difference). (D–H) Five targeted drugs half maximal inhibitory concentration prediction.

RNA and protein expression status of the SUMOylation risk model genes

At the RNA level, CDCA8 was overexpressed, while CDH1 and PPARA were downregulated in kidney cancer tissues, in TCGA-KIRC (Figure 8A). PPARA was downregulated at both the RNA (Figure 8A) and protein levels (Figure 8F). CDH1 was downregulated in cancer tissues at both the RNA (Figure 8A) and protein levels in the CPTAC (Figure 8B) and Chinese Fudan cohort which detecting by proteome sequencing (Figure 8C), and TMA tissues (Figure 8E). However, CDCA8 protein was downregulated in the CPTAC kidney cancer tissues (Figure 8B) and TMA tissues (Figure 8D), which is contrary to the RNA level (Figure 8A).

FIGURE 8

FIGURE 8. The RNA and protein expression of the SUMOylation risk model genes. (A) The RNA expression values of cell division cycle-related protein 8 (CDCA8) (p < 0.001), CDH1 (p < 0.001), and peroxisome proliferator-activated receptor alpha (PPARA) genes (p < 0.001), between normal and kidney cancer tissues from the TCGA-KIRC dataset. (B) CDCA8 and CDH1 protein expression of the normal and tumor tissues in the National Cancer Institute’s Clinical Proteomic Tumor Analysis Consortium cohort dataset. (C) CDH1 protein expression status in Chinese Fudan cohort. (D) The IHC detection of CDCA8 protein expression in TMA tissues. Left panel: the bar plot shows the CDCA8 expression status detected by IHC. Middle panel: The whole TMA slide of CDCA8 detection; the tissue with the blue rectangle means these tissues were excluded from statistics because the tissue was too little. Right panel: displayed detailed CDCA8 statistic data. (E) The IHC detection of CDH1 protein expression in TMA tissues. Left panel: the bar plot shows the CDH1 expression status detected by IHC. Middle panel: The whole TMA slide of CDH1 detection; the tissue with the blue rectangle means these tissues were excluded from statistics because the tissue was too little. Right panel: displayed detailed CDH1 statistic data. (F) The IHC detection of PPARA protein expression in TMA tissues. Left panel: the bar plot shows the PPARA expression status detected by IHC. Middle panel: The whole TMA slide of PPARA detection; the tissue with the blue rectangle means these tissues were excluded from statistics because the tissue was too little. Right panel: displayed detailed PPARA statistic data.

Discussion

Here, we analyzed the SUMOylation-related gene expression status at the RNA level in kidney cancer tissues and built a risk model using three SUMOylation genes, CDH1, CDCA8, and PPARA, to indicate the outcomes based on differentially expressed SUMOylation genes. The risk model also presented as a potential biomarker to reflect the sensitivity of targeted therapy and immune status in tumor tissues. Our study analyzed the protein expression data of these three SUMOylation risk model genes.

Few studies previously focused on the role of SUMOylation in kidney cancer. Overexpression of SUMO-Specific Protease 1 (SENP1), a SUMO protease, could reduce the SUMOylation and ubiquitination of HIF2α, and then increase the local invasion and metastasis capacity of ccRCC cells (Lee et al., 2022). While SENP1 deficiency could lead to deSUMOyation of HIF1α and then inhibit cell proliferation (Dong et al., 2016). In addition, the SUMOylation of hypoxia-associated factor (HAF) could mediate HIF1α degradation, and promoter HIF2α-dependent transcription, and then promote ccRCC development and morbidity (Koh et al., 2015). In our study, the GSVA analysis showed proteasome pathways were significantly highly enriched in kidney cancer tissues in the high-risk group (Figure 6A). Proteasomes take part in the degrading protein process, and it was involved in lots of cellular functions, the dysregulation of proteasomes could result in uncontrolled growth, immune escape, drug resistance, and EMT of cancer (Chen et al., 2017). The poor prognosis of kidney patients in the high-risk group might be associated with the abnormal proteasome in tumor tissue.

CDH1 is also known as E-Cadherin, and its mutation has proven to correlate with tumorigenesis, progression, and invasion of gastric, breast, and colorectal cancers. In addition, the CDH1 protein is a subunit of an E3 ubiquitin ligase, an anaphase-promoting complex or cyclostome (APC/C), and can participate in the mechanism of P53 regulating the G2 checkpoint through SUMOylation (Wang et al., 2020). Previous studies have shown that HIF can mediate suppression of CDH1 gene transcription, and result in the downregulation of E-cadherin in VHL-deficient kidney tumor tissues. In addition, one study showed that hypermethylation of the E-cadherin promoter also contributes to the downregulation of CDH1 (Dulaimi et al., 2004; Esteban et al., 2006). These results may partly explain the reason that CDH1 was downregulated both at the RNA and protein levels. Moreover, Zhang et al. (2017) revealed that reduced E-cadherin in kidney cancer could facilitate tumor progression by activating the WNT/β-catenin pathway. PPARA is a transcription factor, and it can modulate lipid, glucose, and amino acid metabolism, regulate inflammation, and display a tumor suppressor or oncogene role in many cancer types (Tan et al., 2021). PPARA can be SUMOylated and results in the function nuance regulation (Wadosky and Willis, 2012). Indeed, the SUMOylation of PPARA should result in a decrease in its transcriptional activity. In kidney cancer, inhibiting PPARA through an antagonist or siRNA could induce apoptosis and cell cycle arrest of kidney cancer cell lines, and it has been proven to be a diagnostic and prognostic marker for ccRCC (Abu Aboud et al., 2013). Although our study showed that PPARA is downregulated in kidney cancer tissues, Omran et al. revealed that PPARA protein expression was correlated with Grade, where it had a higher expression in grade 4 tissues than in grade 1 tissues (Abu Aboud et al., 2013). The CDCA8 gene encodes a protein that takes part in the composing of chromosomal passenger complex. CDCA8 is upregulated in bladder and prostate cancer tissues, and knockdown of CDCA8 inhibits tumor cell proliferation (Gao et al., 2020; Xu et al., 2022). In addition, CDCA8 expression status could predict the prognosis of prostate and liver cancers (Shuai et al., 2021; Xu et al., 2022). CDCA8 was identified as a core gene involved in the metastasis of ccRCC cells. (Peng et al., 2020). The function of CDCA8 in kidney cancers should be further investigated.

In this study, we determined that CDH1 (Figures 8A–C, E), and PPARA (Figures 8A, F) were downregulated both at the RNA and protein levels. However, CDCA8 protein was downregulated in CPTAC kidney cancer tissues assessed by proteome sequencing (Figure 8B) and in kidney cancer TMA tissues by the IHC method (Figure 8D), which is contrary to the RNA level (Figure 8A). Many reasons might contribute to the contradiction between the protein and RNA expression levels of CDCA8. Studies had revealed that MiR-133a-3p played crucial roles in the regulation of CDCA8 in oesophageal cancer (Wang et al., 2022), and miR-133b regulated the CDCA8 expression in lung adenocarcinoma (Hu et al., 2021) through a post-transcriptional regulation manner. However, the mechanism of the contrary expression level of RNA and protein in kidney cancer has not been well discussed, which should be studied in the future.

Our SUMOylation risk model indicated that patients in the high-risk group have a worse outcome. The immune subtype analysis revealed a significantly higher ratio of C1, C2, C4, and C6 subtypes and a lower C3 subtype in the high-risk model group. The previous study proposed that patients in the C3 subtype often had the best prognosis of all subtypes, while patients in the C6 subtype had the worst outcome. Different outcomes of the distinct SUMOylation risk group might be attributable to the varying immune statuses in tumor tissues, which were assessed by the immune subtype and immune cell analysis.

A high TIDE score was associated with the immune escape capacity of tumor cells, while a low IC50 might indicate a sensitivity to drug treatment. In this study, tissues of the high-risk group presented a high TIDE score and a low IC50 of the five targeted drugs, which implied that patients with high-risk scores might be more suitable for receiving the ICB plus targeted agent treatment as the first-line treatment.

In our study, we examined the RNA expression status of SUMOylation genes in kidney cancer tissues and constructed and verified a SUMOylation prognosis model for predicting the outcome of kidney cancer using three databases and five cohorts. In conclusion, the SUMOylation risk model may be a biomarker for the selection of treatment drugs for kidney cancer.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The studies involving human participants were reviewed and approved by the Shanghai Outdo Biotech Company. The patients/participants provided their written informed consent to participate in this study.

Author contributions

Conceived and designed the experiment: X-HN and S-CL; Performed the experiments and analyzed the data: X-HN, S-CL, Z-KJ, L-JY, and X-LW; Interpretation of the findings: X-HN, S-CL, Z-KJ, and J-JY; All the authors contributed to writing the manuscript.

Funding

This work was supported by the Postdoctoral Research Grant in Henan Province (grant number 1901004), the Henan Science and Technology Research Program (grant number 2018020142), The Natural Science Foundation of Henan Province (212300410265), and the Key scientific research projects of colleges and universities in Henan Province (23A320038). The funders had no role in the study design, data collection, analysis, the decision to publish, or the preparation of the manuscript. We thank Chen-yang Sima for his kind help in data analysis.

Acknowledgments

The results shown here are in whole or part based upon data generated by the TCGA Research Network (https://www.cancer.gov/tcga), ICGC database (https://dcc.icgc. org/) and CPTAC (https://proteomic.datacommons.cancer.gov/pdc/study/PDC000127).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2023.1038457/full#supplementary-material

References

Abu Aboud, O., Wettersten, H. I., and Weiss, R. H. (2013). Inhibition of PPARα induces cell cycle arrest and apoptosis, and synergizes with glycolysis inhibition in kidney cancer cells. PLoS One 8, e71115. doi:10.1371/journal.pone.0071115

PubMed Abstract | CrossRef Full Text | Google Scholar

Albiges, L., Powles, T., Staehler, M., Bensalah, K., Giles, R. H., Hora, M., et al. (2019). Updated European association of urology guidelines on renal cell carcinoma: Immune checkpoint inhibition is the new backbone in first-line treatment of metastatic clear-cell renal cell carcinoma. Eur. Urol. 76, 151–156. doi:10.1016/j.eururo.2019.05.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Baek, S. H., Ko, J. H., Lee, J. H., Kim, C., Lee, H., Nam, D., et al. (2017). Ginkgolic acid inhibits invasion and migration and TGF-beta-induced EMT of lung cancer cells through PI3K/Akt/mTOR inactivation. J. Cell. Physiol. 232, 346–354. doi:10.1002/jcp.25426

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertolotto, C., Lesueur, F., Giuliano, S., Strub, T., de Lichy, M., Bille, K., et al. (2011). A SUMOylation-defective MITF germline mutation predisposes to melanoma and renal carcinoma. Nature 480, 94–98. doi:10.1038/nature10539

PubMed Abstract | CrossRef Full Text | Google Scholar

Bettermann, K., Benesch, M., Weis, S., and Haybaeck, J. (2012). SUMOylation in carcinogenesis. Cancer Lett. 316, 113–125. doi:10.1016/j.canlet.2011.10.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Celen, A. B., and Sahin, U. (2020). Sumoylation on its 25th anniversary: Mechanisms, pathology, and emerging concepts. FEBS J. 287, 3110–3140. doi:10.1111/febs.15319

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Zhang, Y., and Guo, X. (2017). Proteasome dysregulation in human cancer: Implications for clinical therapies. Cancer Metastasis Rev. 36, 703–716. doi:10.1007/s10555-017-9704-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, D. J., Dhanasekaran, S. M., Petralia, F., Pan, J., Song, X., Hu, Y., et al. (2019). Integrated proteogenomic characterization of clear cell renal cell carcinoma. Cell. 179, 964–983 e31. doi:10.1016/j.cell.2019.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, B., Gao, Y., Kang, X., Gao, H., Zhang, J., Guo, H., et al. (2016). SENP1 promotes proliferation of clear cell renal cell carcinoma through activation of glycolysis. Oncotarget 7, 80435–80449. doi:10.18632/oncotarget.12606

PubMed Abstract | CrossRef Full Text | Google Scholar

Dulaimi, E., Ibanez de Caceres, I., Uzzo, R. G., Al-Saleem, T., Greenberg, R. E., Polascik, T. J., et al. (2004). Promoter hypermethylation profile of kidney cancer. Clin. Cancer Res. 10, 3972–3979. doi:10.1158/1078-0432.CCR-04-0175

PubMed Abstract | CrossRef Full Text | Google Scholar

Esteban, M. A., Tran, M. G., Harten, S. K., Hill, P., Castellanos, M. C., Chandra, A., et al. (2006). Regulation of E-cadherin expression by VHL and hypoxia-inducible factor. Cancer Res. 66, 3567–3575. doi:10.1158/0008-5472.CAN-05-2670

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, X., Wen, X., He, H., Zheng, L., Yang, Y., Yang, J., et al. (2020). Knockdown of CDCA8 inhibits the proliferation and enhances the apoptosis of bladder cancer cells. PeerJ 8, e9078. doi:10.7717/peerj.9078

PubMed Abstract | CrossRef Full Text | Google Scholar

Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9, e107468. doi:10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, Q., Xiao, X. Y., Wu, C. Y., Li, D., Chen, J. L., Ding, X. C., et al. (2022). Clinical roles of risk model based on differentially expressed genes in mesenchymal stem cells in prognosis and immunity of non-small cell lung cancer. Front. Genet. 13, 823075. doi:10.3389/fgene.2022.823075

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, C., Wu, J., Wang, L., Liu, X., Da, B., Liu, Y., et al. (2021). miR-133b inhibits cell proliferation, migration, and invasion of lung adenocarcinoma by targeting CDCA8. Pathol. Res. Pract. 223, 153459. doi:10.1016/j.prp.2021.153459

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24, 1550–1558. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Koh, M. Y., Nguyen, V., Lemos, R., Darnay, B. G., Kiriakova, G., Abdelmelek, M., et al. (2015). Hypoxia-induced SUMOylation of E3 ligase HAF determines specific activation of HIF2 in clear-cell renal cell carcinoma. Cancer Res. 75, 316–329. doi:10.1158/0008-5472.CAN-13-2190

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, M. H., Sung, K., Beebe, D., Huang, W., Shapiro, D., Miyamoto, S., et al. (2022). The SUMO protease SENP1 promotes aggressive behaviors of high HIF2α expressing renal cell carcinoma cells. Oncogenesis 11, 65. doi:10.1038/s41389-022-00440-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Meng, X., Zhang, D., Xu, X., Li, S., and Li, Y. (2019). Ginkgolic acid suppresses the invasion of HepG2 cells via downregulation of HGF/cMet signaling. Oncol. Rep. 41, 369–376. doi:10.3892/or.2018.6786

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q. K., Pavlovich, C. P., Zhang, H., Kinsinger, C. R., and Chan, D. W. (2019). Challenges and opportunities in the proteomic characterization of clear cell renal cell carcinoma (ccRCC): A critical step towards the personalized care of renal cancers. Semin. Cancer Biol. 55, 8–15. doi:10.1016/j.semcancer.2018.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S. C., Jia, Z. K., Yang, J. J., and Ning, X. H. (2022). Telomere-related gene risk model for prognosis and drug treatment efficiency prediction in kidney cancer. Front. Immunol. 13, 975057. doi:10.3389/fimmu.2022.975057

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi:10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Ning, X. H., Li, N. Y., Qi, Y. Y., Li, S. C., Jia, Z. K., and Yang, J. J. (2021). Identification of a hypoxia-related gene model for predicting the prognosis and formulating the treatment strategies in kidney renal clear cell carcinoma. Front. Oncol. 11, 806264. doi:10.3389/fonc.2021.806264

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, R., Wang, Y., Mao, L., Fang, F., and Guan, H. (2020). Identification of core genes involved in the metastasis of clear cell renal cell carcinoma. Cancer Manag. Res. 12, 13437–13449. doi:10.2147/CMAR.S276818

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiao, L., Zheng, J., Jin, X., Wei, G., Wang, G., Sun, X., et al. (2017). Ginkgolic acid inhibits the invasiveness of colon cancer cells through AMPK activation. Oncol. Lett. 14, 5831–5838. doi:10.3892/ol.2017.6967

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, Y., Feng, J., Wu, X., Bai, L., Xu, W., Zhu, L., et al. (2022). A proteogenomic analysis of clear cell renal cell carcinoma in a Chinese population. Nat. Commun. 13, 2052. doi:10.1038/s41467-022-29577-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramazi, S., and Zahiri, J. (2021). Posttranslational modifications in proteins: Resources, tools and prediction methods. Database (Oxford) 2021, baab012. doi:10.1093/database/baab012

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahin, U., de The, H., and Lallemand-Breitenbach, V. (2022). Sumoylation in physiology, pathology and therapy. Cells 11, 814. doi:10.3390/cells11050814

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, Y., Yoshizato, T., Shiraishi, Y., Maekawa, S., Okuno, Y., Kamura, T., et al. (2013). Integrated molecular analysis of clear-cell renal cell carcinoma. Nat. Genet. 45, 860–867. doi:10.1038/ng.2699

PubMed Abstract | CrossRef Full Text | Google Scholar

Schodel, J., Grampp, S., Maher, E. R., Moch, H., Ratcliffe, P. J., Russo, P., et al. (2016). Hypoxia, hypoxia-inducible transcription factors, and renal cancer. Eur. Urol. 69, 646–657. doi:10.1016/j.eururo.2015.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Shuai, Y., Fan, E., Zhong, Q., Chen, Q., Feng, G., Gou, X., et al. (2021). CDCA8 as an independent predictor for a poor prognosis in liver cancer. Cancer Cell. Int. 21, 159. doi:10.1186/s12935-021-01850-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, Y., Wang, M., Yang, K., Chi, T., Liao, Z., and Wei, P. (2021). PPAR-Α modulators as current and potential cancer treatments. Front. Oncol. 11, 599995. doi:10.3389/fonc.2021.599995

PubMed Abstract | CrossRef Full Text | Google Scholar

Thorsson, V., Gibbs, D. L., Brown, S. D., Wolf, D., Bortone, D. S., Ou Yang, T. H., et al. (2018). The immune landscape of cancer. Immunity 48, 812–830.e14. doi:10.1016/j.immuni.2018.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Tran, A. N., Dussaq, A. M., Kennell, T., Willey, C. D., and Hjelmeland, A. B. (2019). HPAanalyze: an R package that facilitates the retrieval and analysis of the Human Protein Atlas data. BMC Bioinforma. 20, 463. doi:10.1186/s12859-019-3059-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wadosky, K. M., and Willis, M. S. (2012). The story so far: Post-translational regulation of peroxisome proliferator-activated receptors by ubiquitination and SUMOylation. Am. J. Physiol. Heart Circ. Physiol. 302, H515–H526. doi:10.1152/ajpheart.00703.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Zhu, L., Lin, X., Huang, Y., and Lin, Z. (2022). MiR-133a-3p inhibits the malignant progression of oesophageal cancer by targeting CDCA8. J. Biochem. 170, 689–698. doi:10.1093/jb/mvab071

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Tian, J., Huang, C., Ma, J., Hu, G., Chen, Y., et al. (2020). P53 suppresses SENP3 phosphorylation to mediate G2 checkpoint. Cell. Discov. 6, 21. doi:10.1038/s41421-020-0154-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, J. H., Feng, Z. H., Cao, Y., Zhao, H. W., Chen, Z. H., Liao, B., et al. (2019). Predictive value of single-nucleotide polymorphism signature for recurrence in localised renal cell carcinoma: A retrospective analysis and multicentre validation study. Lancet Oncol. 20, 591–600. doi:10.1016/S1470-2045(18)30932-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, Y., Peng, Y., Shen, M., Liu, L., Lei, J., Gao, S., et al. (2022). Construction and validation of angiogenesis-related prognostic risk signature to facilitate survival prediction and biomarker excavation of breast cancer patients. J. Oncol. 2022, 1525245. doi:10.1155/2022/1525245

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, W., Soares, J., Greninger, P., Edelman, E. J., Lightfoot, H., Forbes, S., et al. (2013). Genomics of drug sensitivity in cancer (GDSC): A resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. 41, D955–D961. doi:10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Yang, M., Shi, H., Hu, J., Wang, Y., Sun, Z., et al. (2017). Reduced E-cadherin facilitates renal cell carcinoma progression by WNT/β-catenin signaling activation. Oncotarget 8, 19566–19576. doi:10.18632/oncotarget.15361

PubMed Abstract | CrossRef Full Text | Google Scholar