Identification and Experimental Validation of Prognostic Signature and Peroxisome-Related Key Genes in Clear Cell Renal Cell Carcinoma

Introduction

Presently, renal cell carcinoma (RCC) is considered a prevalent tumor in adults, with increasing morbidity and mortality.^1,2 Clear cell RCC (ccRCC) makes up 70–75% of RCC cases and is responsible for the majority of kidney cancer-related deaths.^3,4 Although there has been notable progress in the treatment of ccRCC, the outcomes are still unsatisfactory.⁵ In the case of advanced ccRCC, conventional radiotherapy or chemotherapy generally fails to produce satisfactory therapeutic results.⁶ Additionally, the efficacy of targeted therapy for ccRCC is still limited.⁷ Hence, exploring novel potential signatures is essential for improving the prognostic prediction and individualized treatment of ccRCC.

Peroxisomes are widely present in eukaryotic cells and participate in various biological processes.^8,9 In recent years, peroxisomes have become key regulators of inflammation and immune responses.^10,11 Previously, increasing studies have demonstrated that peroxisome dysregulation contributes to cancer biology, including oxidative stress, lipid metabolism, and immune modulation. Catalase, a member of peroxisome, can neutralize reactive oxygen species (ROS).¹² The dysfunctions of peroxisomes disrupt the balance between ROS production and detoxification, leading to oxidative stress, which influences cancer progression.¹² Moreover, peroxisomes play a role in lipid metabolism.¹³ Dysregulation of peroxisomal lipid metabolism may result in the accumulation of lipid metabolic intermediates, which drives the metabolic reprogramming of cancer.¹⁴ Additionally, peroxisomes have been found to be involved in the development and function of immune cells in cancer.¹⁵ Peroxisome dysregulation may foster an immunosuppressive tumor microenvironment, promoting tumor growth and evasion of immune surveillance.¹⁶ Mounting evidence reveals that peroxisome is correlated with the progression of cancers. In hepatocellular carcinoma, PPAR-γ expression is significantly correlated with tumor size, metastatic recurrence, and tumor invasion.¹⁷ Dong et al found that PEX13 (a member of peroxisome) knockdown could attenuate the proliferation, migration, and invasive ability of pancreatic adenocarcinoma cells.¹⁸ In addition, it has been reported that PPAR-α is involved in metabolic reprogramming of RCC.¹⁹ However, the prognostic role of peroxisome-related genes (PGs) in ccRCC has not been clarified. Developing peroxisome-related risk models may provide individualized treatment for ccRCC patients.

In this study, a risk model, including 9 PGs, with good prognostic prediction ability was established. The predict value of prognosis was validated across multiple cohorts (TCGA, GSE167573 and local cohorts). Furthermore, we have assessed the association between risk model and various factors (survival outcome, pathological features, immune status, and chemotherapy). Additionally, the expression patterns of 9 PGs between tumor and normal tissues were examined via local clinical samples. Overall, the 9 PGs risk model was correlated with immune response and may be used independently to predict prognosis for ccRCC patients.

Materials and Methods

Data Preparation

As the training cohort (TCGA), transcriptome data and associated clinical data were obtained from the Cancer Genome Atlas database. GSE167573 cohort was downloaded from GEO databases. The 68 samples, collected from the First Affiliated Hospital of Wenzhou Medical University, were used as the external validation cohort (FAHWMU cohort). The Molecular Signatures Database (MsigDB) provides the 104 PGs. Differentially expressed genes (DEGs) were identified through the “limma” (package version 3.38.3) and follow the significant threshold as Log₂|Fold change| > 1.5 and p < 0.05.

Establishment of Risk Model

Univariate Cox regression analysis was used to identify the peroxisome-related DEGs, which are related to prognosis. The least absolute shrinkage selection operator (LASSO) Cox regression analysis was further performed to construct risk model. The risk score was calculated by summing the product of biased regression coefficients and transcriptional values for each gene (). According to the median value of the risk score, all samples with ccRCC were classified as either high- or low-risk. Dimensionality reduction cluster analysis, comprising Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis, was carried out using the “Rtsne” R package. The analysis of the Kaplan-Meier survival curve was performed by the “survival” R package. Through the use of the R package “survival ROC”, the receiver operating characteristic (ROC) curves were determined.

Evaluation of Immune Status

The relative infiltration levels of immune cell types were explored via Single sample gene set enrichment analysis (ssGSEA). Immune subtypes, used to assess immune infiltration within the tumor microenvironment (TME), were explored via R package “Immune Subtype Classifier”.²⁰ The relationships between the risk score and the genes related to immune checkpoint were assessed by Pearson correlation analysis. The immune checkpoint gene expression levels in high- and low-risk groups were compared using the Wilcoxon test.

Functional Enrichment Analyses

Enrichment analyses based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Hallmark pathway in GSEA software 4.1. Statistical significance is defined as P < 0.05 and a false discovery rate < 0.25.

Drug Sensitivity

The IC50 data for a total of 746 small molecules in cell lines were collected from the Genomics of Therapeutics Response Portal (CTRP) and the Genomics of Drug Sensitivity in Cancer (GDSC). The association between 9 PGs and drug sensitivity was investigated using Pearson correlation analysis.

Quantitative Real-Time PCR (qRT-PCR)

Briefly, total RNA was extracted with the Trizol reagent (Servicebio). Then cDNA was synthesized using the cDNA Synthesis Kit. RT-PCR was conducted based on the utility of 7500 Fast Real-Time PCR System.

Immunohistochemistry (IHC)

Sections of tissue were fixed in 4% formalin, followed by degreasing and rehydration. After blocking with 10% bovine serum albumin, specimens were incubated overnight at 4°C with the primary antibody. Nuclei were counterstained by hematoxylin. Further, sections were visualized by DAB color development kits. Images were captured by a Leica DM4B microscope.

Statistical Analyses

All statistical analyses were performed in R software or GraphPad Prism. The Wilcoxon test was applied to compare the expression levels of 9 PGs expression among tumor and normal tissues. Unpaired two-tailed Student’s t-test was used for comparison between two groups. Among these, P < 0.05 was considered to be statistically significant.

Results

Establishment of Risk Model

20 differentially expressed peroxisome-related genes (DEPGs) were identified (Figure 1A). For ccRCC, 13 DEPGs were found to be linked to the prognosis according to univariate Cox regression analysis (Figure 1B). The correlation network depicted the interactions among the 13 candidate DEPGs (Figure 1C). Subsequently, a risk model, encompassing 9 PGs, was constructed via LASSO analysis. The heatmap showed the expression patterns of the 9 PGs between ccRCC and adjacent normal tissues (Figure 1D). The coefficient for each gene included in the risk model was provided in Supplementary Table S1. Additionally, all 9 PGs included in the risk model were found to be linked to the prognosis of ccRCC (Supplementary Figure S1). Using the median risk score of the TCGA cohort as the threshold, patients were categorized into the high- and low-risk groups (Figure 2A). Furthermore, it was found that patients with high risk score exhibited more unsatisfied survival time than individuals with low risk score (Figure 2B). Additionally, dimensionality reduction and clustering techniques, such as PCA and t-SNE, showed distinct separation between high-risk and low-risk groups (Figure 2C and D). The Kaplan-Meier survival curve indicated better prognosis for ccRCC patients with low risk score (Figure 2E). With AUC values of 0.739, 0.709, and 0.712 for 1, 2, and 3 years, respectively, time-dependent ROC analysis further certificated the risk model’s strong predictive accuracy (Figure 2F). Moreover, as depicted in Figure 3A, univariate Cox regression analysis demonstrated that the risk score was significantly linked to the prognosis within the context of ccRCC. Obviously, the risk score served as an independent prognostic parameter in the TCGA cohort via multivariate Cox regression analysis (Figure 3B).

Figure 1 Identification of candidate peroxisome-related genes in the TCGA cohort. (A) Venn diagram to identify DEPGs between ccRCC and adjacent normal tissues. (B) Univariate Cox regression analysis identifying 13 DEPGs. (C) The correlation network of candidate genes. (D) Expression of the 9 PGs between ccRCC and adjacent normal tissues.

Figure 2 Prognostic analysis of the risk model in the TCGA cohort. (A) The distribution and median value of the risk scores. (B) Distributions of the overall survival (OS) status. (C) PCA plot. (D) t-SNE analysis. (E) Kaplan-Meier curves for OS of patients in the high-risk group and low-risk group. (F) Time-dependent ROC curves for OS.

Figure 3 The univariate and multivariate Cox regression analyses regarding TCGA cohort. (A) Univariate Cox regression analysis. (B) Multivariate Cox regression analysis.

Validation of the Prognostic Model in the External Cohorts

To assess the robustness of the prognostic model, independent external cohorts (GSE167573 and FAHWMU cohorts) were utilized for validation. Following the median value, patients with ccRCC in the GSE167573 and FAHWMU cohorts were divided into two groups (high-risk group or low-risk group) (Supplementary Figure S2A and S2C). Consistent with the aforementioned findings in the TCGA cohort, the high-risk group exhibited significantly shorter survival time compared to the low-risk group in both the GSE167573 and FAHWMU cohorts (Supplementary Figure S2B and S2D). PCA and t-SNE analyses revealed a distinct dichotomous distribution of samples into different separate risk groups within the external cohorts (Supplementary Figure S2E–H). In both the GSE167573 and FAHWMU cohorts, ccRCC patients in the high-risk group had a worse outcome (Supplementary Figure S2I and S2K). The AUC values at 1, 2, and 3 years in the GSE167573 cohort were 0.794, 0.766, and 0.823, respectively (Supplementary Figure S2J). In the FAHWMU cohorts, the AUC values at 1, 2, and 3 years were 0.801, 0.739, and 0.737, respectively (Supplementary Figure S2L). Overall, the risk model demonstrated a stable predictive capability in the external cohorts.

Risk Score in Different Clinical Characteristics Groups

The differences in risk score among various clinicopathological characteristic groups were investigated in the TCGA cohort (Figure 4A–G). The ccRCC patients with unfavorable pathological features (worse tumor stage, higher histological grade, and a raised TNM stage) were shown to have a higher risk score. These findings suggest that the risk score may be positively associated with clinicopathological status. Additionally, the ROC analysis supported that the 9 PGs risk model showed promising predictive capability (AUC = 0.738) in comparison with other clinical factors (Figure 4H).

Figure 4 The risk score in different groups stratified by clinical characteristics. (A) Age. (B) Gender. (C) Tumor grade. (D) Tumor stage. (E) Primary tumor size and extent stage. (F) Tumor lymph node stages. (G) Tumor distant metastasis stages. (H) The multifactor AUC. *P < 0.05, **P < 0.01, ****P < 0.0001, - P > 0.05.

Relationship Between Risk Score and Immune Status

The infiltration landscape of immune cells in the high- and low-risk groups was analyzed to investigate the impact of the risk score on the TME. There were variations in the abundance of multiple immune cell types, such as CD8⁺ T cells, Macrophages, Mast cells, Tumor infiltrating lymphocytes (TILs), and Regulatory T cells (Tregs), between the different risk groups in the TCGA cohort. Specifically, the high-risk group exhibited enrichment of CD8⁺ T cells, Macrophages, pDCs, T helper cells, Tfh, Th1 cells, Th2 cells, TILs, and Tregs (Figure 5A). In terms of immune function, there was a significantly stronger type II IFN response in the low-risk group compared to the high-risk group (Figure 5B). Furthermore, six immune subtypes, denoted as C1 (wound healing), C2 (INF-γ dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGF-β dominant),²¹ were defined to assess immune infiltration. As depicted in Figure 5C, the high-risk score exhibited a significant correlation with C1, C2, and C6, while the low-risk score was correlated with C3 and C4.

Figure 5 Association between risk score and immune status. (A) The scores of 16 immune cells. (B) The boxplots showing the 13 immune-related functions. (C) Comparison of the risk scores between different immune infiltrate subtypes. *P < 0.05, **P < 0.01, ****P < 0.0001, - P > 0.05.

Association Between Risk Score and Immune Checkpoints

Meanwhile, associations between immune checkpoint gene expression levels across risk categories were investigated. As for TCGA cohort, it was observed that CTLA4, LAG3, PDCD1, and TIGIT exhibited higher expression levels in the high-risk group compared to the low-risk group (Figure 6A). Moreover, the expression levels of CTLA4, LAG3, PDCD1, and TIGIT demonstrated a positive correlation with the risk score (Figure 6B–E). Accordingly, these results were also observed in the GSE167573 and FAHWMU cohorts (Supplementary Figure S3). Collectively, these findings indicate that the risk model effectively distinguishes the expression levels of immune checkpoint genes.

Figure 6 Association between risk score and immune checkpoint genes. (A) The expression of immune checkpoint genes in high- and low- risk groups. (B) Correlation between risk score and CTLA4. (C) Correlation between risk score and LAG3. (D) Correlation between risk score and PDCD1. (E) Correlation between risk score and TIGIT. ****P < 0.0001.

Functional Enrichment Analysis

To investigate the potential correlation between the risk model and biological functions as well as pathways, we conducted GSEA enrichment analyses on the KEGG and HALLMARK gene sets within the TCGA cohort. For KEGG enrichment gene sets, the p53 signaling pathway and the primary immunodeficiency pathway were shown to be considerably enriched in the high-risk group (Supplementary Figure S4A). As for the HALLMARK gene sets, the G2M checkpoint pathway, KRAS signaling pathway, and IL6-JAK-STAT3 signaling pathway were enriched in the group with high risk (Supplementary Figure S4B). Notably, the G2M checkpoint pathway plays a pivotal role in the regulation of cell cycle, which is closely linked to cancer progression. Hence, the association between G2M checkpoint pathway-related genes and risk score was analyzed. As depicted in Supplementary Figure S5, CDC25A, CDK1, and CDK2 exhibited elevated expression levels in the high-risk group and showed a positive correlation with the risk score.

Drug Sensitivity

In order to discover more effective and advanced medications for the treatment of ccRCC patients, as well as to develop additional drugs to improve the current limited situation, the correlations among the expressions of 9 prognostic PGs and the drug sensitivity of anti-tumor medications were assessed. As illustrated in Figure 7, the expressions of the prognostic 9 genes exhibited a correlation with the sensitivity of several medications, such as Lapatinib, Panobinostat, Mitomycin, and Docetaxel. Based on these results, the risk model may be used as a chemosensitivity predictor.

Figure 7 Relationship analysis between mRNA expression of genes and drug sensitivity. (A) CTRP database. (B) GDSC databases.

Validation of the Expressions of 9 PGs

The mRNA and protein expression levels of the 9 prognostic genes were assessed in ccRCC and the surrounding normal tissues. Figure 8A demonstrated that CEL, CRABP2, HSD3B7, SULT2B1, and TOP2A exhibited high expression in ccRCC tissues, whereas ABCB4, EPHX2, HAO2, and HSD11B2 showed low expression in ccRCC tissues. The IHC images further confirmed the increased expression of CEL, CRABP2, HSD3B7, SULT2B1, and TOP2A in ccRCC tissues (Figure 8B). Because of its high Hazard ratio value, TOP2A was chosen for validation. si-TOP2A#1 was applied in the subsequent experiments due to its better knockdown efficiency than si-TOP2A#2 (Figure 9A). The sequences of si-TOP2A and primers are listed in Supplementary Table S2. Obviously, the silence of TOP2A attenuated the proliferation, migration, and invasion ability of RCC cells via CCK-8 assay, wound healing assay and transwell invasion assay (Figure 9B–D). In addition, the expression levels of the prognostic genes were evaluated across various cancer types (Supplementary Figure S6). Notably, CEL, CRABP2, HSD3B7, TOP2A, HSD11B2, and SULT2B1 exhibited significantly higher levels in most cancer types, while ABCB4, EPHX2, and HAO2 displayed lower expression in majority cancers. The result further supports the prospective use of this model in other tumor investigations by indicating that most prognostic gene expression patterns in other malignancies are similar to those found in ccRCC.

Figure 8 Verification of the expression level of prognostic genes between ccRCC and adjacent normal tissues. (A) The mRNA expression level of prognostic genes in ccRCC and adjacent normal tissue detected by qRT-PCR. (B) Representative IHC images of prognostic genes in tumor and adjacent normal tissues. **P < 0.01, ***P < 0.001.

Figure 9 Effects of TOP2A knockdown on RCC cells. (A) mRNA expression level of TOP2A after treatment of si-TOP2A#1/#2 and si-NC. (B) CCK-8 assay. (C) Wound healing assay. (D) Transwell assay. ***P < 0.001.

Discussion

Among all urologic malignancies, ccRCC has the highest fatality rate.²² Advanced stage and metastasized ccRCC cancers have an unsatisfactory prognosis.²³ Hence, it is imperative to discover novel prognostic markers to attain prognostic prediction.

In this study, a risk model comprising 9 PGs (ABCB4, CEL, CRABP2, EPHX2, HAO2, HSD3B7, HSD11B2, SULT2B1, and TOP2A) was identified. Among these 9 prognostic genes, ABCB4, a member of the adenosine triphosphate-binding cassette, plays a regulatory role in chemoresistance in colorectal cancer.²⁴ The presence of CRABP2 (a member of the lipid calcium-binding protein/cytosolic fatty acid binding protein family) promotes the activation of pro-oncogenic pathways in various tumors such as breast cancer, neuroblastoma, and ovarian cancer.^25–27 EPHX2, which encodes soluble epoxide hydrolase (sEH) in vivo, is linked to the prognostic prediction of endometrial cancer and neuroblastoma.^28,29 sEH inhibitors are being considered as a potential treatment modality for renal-associated diseases.³⁰ Xiao et al reported that HAO2 is involved in inhibiting malignancy in ccRCC by promoting the lipid catabolic process.³¹ Recently, HSD3B7 was found to be associated with tumorigenesis in ccRCC by participating in bile acid biosynthesis.³² HSD11B2 emerges as a prognostic factor in colon cancer.³³ Recent research has revealed that up-regulated expression of SULT2B1 enhances the capability of colorectal carcinoma cells’ proliferation and invasion in vitro.³⁴ TOP2A has been reported to catalyze DNA double-strand breaks and be associated with a variety of dynamic processes.³⁵ Its aberrant expression contributes to the carcinogenesis of various malignancies, including ccRCC.^36–38 These 9 prognostic genes take part in the development of multiple cancers, suggesting that they could be useful as prognostic indicators for patients with ccRCC.

Herein, in diverse cohorts, survival analysis revealed that ccRCC patients with high risk score had poorer survival outcomes, indicating its excellent performance in prognostic prediction. Notably, the expression patterns of the 9 PGs were examined in clinical samples obtained from a local hospital. Additionally, a substantial correlation was found between the risk score and the clinicopathological status. Moreover, positive correlations were observed between G2M checkpoint-related genes (CDC25A, CDK1, and CDK2) and the risk score. Dysregulation of the G2M checkpoint (a critical component of the cell cycle) implicates RCC development.³⁹ The aberrant expressions of CDC25A, CDK1, and CDK2 are involved in the tumorigenesis and cancer progression.^40–43 The Above findings imply that the risk score may be associated with an abnormal cell cycle in ccRCC.

Presently, tumor immunology is of great interest. The TME is closely associated with the development of ccRCC.⁴⁴ It has been reported that the combinations of immunoscore and tumor grading as well as tumor staging contribute to the prognosis prediction of ccRCC patients.⁴⁵ Meanwhile, it was found that the predictive role of immunoscore differed in the non-ccRCC subtype versus ccRCC.⁴⁶ Besides, immune escape and medication resistance are significantly influenced by the immune microenvironment.⁴⁷ Differences were observed between the high- and low-risk groups in terms of immune cell subtypes, including CD8⁺ T cells, macrophages, and Treg cells. Recent studies have identified RCC as a highly infiltrated tumor by CD8⁺ T cells.⁴⁸ However, tumor-associated macrophages (TAMs), particularly M2 macrophages, tend to undermine the efficacy of immunotherapy.⁴⁹ The infiltration of TAMs subtypes and DC cells is associated with poor prognosis in RCC.^50–52 Additionally, immune checkpoint blockade therapy is a crucial method for RCC treatment.⁵³ Patients at high risk exhibited elevated levels of CD274, CTLA4, and LAG3. It implies that the risk score could potentially enhance the efficacy of personalized therapeutic approaches for ccRCC. Combining the above results, we reveal that the risk model may predict immunotherapy response and provide suggestions for individual treatment. Also, it was found that the expressions of the 9 prognostic PGs were closely associated with the drug sensitivity of multiple drugs. Among them, Lapatinib, Panobinostat, Mitomycin, and Docetaxel have been reported to be involved in the inhibition of RCC.^54–56 Insightfully, Docetaxel suppresses the development of ccRCC via inhibiting G2/M cell cycle.⁵⁶ Hence, the results of drug sensitivity suggest that our risk model may guide drug selection strategies and drug development for ccRCC patients. In sum, the risk model exhibits its excellent potential in clinical application.

In previous studies, other prognostic markers or gene signatures have been developed to predict the prognosis of ccRCC. Wei et al found that upregulated NIPAL4 is associated with a poor prognosis for ccRCC patients.⁵⁷ Notably, the role of NIPAL4 in prognosis is only validated in a public dataset. Additionally, Lyu et al constructed an immune-stromal score signature, having the AUC predictive values for 1-, 3-, and 5-year survival standing at 0.617, 0.669, and 0.68, respectively.⁵⁸ Compared with the previous studies, the prognostic value of our risk model was validated across three cohorts. Our risk model demonstrated robust performance, achieving AUC values of 0.739, 0.709, and 0.712 for 1-, 3-, and 5-year survival, respectively, suggesting the satisfactory predictive ability of this risk model. Meanwhile, the risk model was closely associated with immune-related factors as well as drug sensitivity. According to the experiments in vitro, we revealed that knockdown of TOP2A diminished the ability of proliferation, migration, and invasion of RCC cells. This finding could support the future development of therapeutic agents targeting TOP2A or the use of TOP2A as a diagnostic marker for ccRCC.

Our study on risk model related to peroxisomes in ccRCC highlights the importance of peroxisomes in cancer prognosis. We found that the 9 PGs are differentially expressed in pan-cancer. Therefore, it may provide insight that future research can explore peroxisomal functions in other cancers, identifying key enzymes and pathways involved as well as determining how their dysregulation contributes to tumorigenesis, metastasis, and treatment resistance.

In conclusion, our peroxisome-related risk model, exhibiting excellent predictive performance, may serve as an independent prognostic indicator for ccRCC, thus offering valuable insights into personalized treatment approaches for ccRCC.