A Novel Cuproptosis-Related Prognostic Gene Signature and Validation of Differential Expression in Clear Cell Renal Cell Carcinoma

Clear cell renal cell carcinoma (ccRCC) is the most prevalent subtype of renal cell carcinoma, which is characterized by metabolic reprogramming. Cuproptosis, a novel form of cell death, is highly linked to mitochondrial metabolism and mediated by protein lipoylation. However, the clinical impacts of cuproptosis-related genes (CRGs) in ccRCC largely remain unclear. In the current study, we systematically evaluated the genetic alterations of cuproptosis-related genes in ccRCC. Our results revealed that CDKN2A, DLAT, DLD, FDX1, GLS, PDHA1 and PDHB exhibited differential expression between ccRCC and normal tissues (|log2(fold change)| > 2/3 and p < 0.05). Utilizing an iterative sure independence screening (SIS) method, we separately constructed the prognostic signature of CRGs for predicting the overall survival (OS) and progression-free survival (PFS) in ccRCC patients. The prognostic score of CRGs yielded an area under the curve (AUC) of 0.658 and 0.682 for the prediction of 5-year OS and PFS, respectively. In the Kaplan−Meier survival analysis of OS, a higher risk score of cuproptosis-related gene signature was significantly correlated with worse overall survival (HR = 2.72 (2.01–3.68), log-rank p = 1.76 × 10−7). Patients with a higher risk had a significantly shorter PFS (HR = 2.83 (2.08–3.85), log-rank p = 3.66 × 10−7). Two independent validation datasets (GSE40435 (N = 101), GSE53757 (N = 72)) were collected for meta-analysis, suggesting that CDKN2A (log2(fold change) = 1.46, 95%CI: 1.75–2.35) showed significantly higher expression in ccRCC tissues while DLAT (log2(fold change) = −0.54, 95%CI: −0.93–−0.15) and FDX1 (log2(fold change) = −1.01, 95%CI: −1.61–−0.42) were lowly expressed. The expression of CDKN2A and FDX1 in ccRCC was also significantly associated with immune infiltration levels and programmed cell death protein 1 (PD-1) expression (CDKN2A: r = 0.24, p = 2.14 × 10−8; FDX1: r = −0.17, p = 1.37 × 10−4). In conclusion, the cuproptosis-related gene signature could serve as a potential prognostic predictor for ccRCC patients and may offer novel insights into the cancer treatment.


Introduction
Renal cell carcinoma (RCC) is one of the most common cancer types in the urinary system affecting more than 430,000 individuals in 2020 worldwide [1]. Clear cell renal cell carcinoma (ccRCC) is the most prevalent and aggressive subtype accounting for about 70% of all RCCs [2]. Clinically, approximately a third of ccRCC patients will present with metastasis at the initial diagnosis, and a quarter of patients with localized disease will have a relapsing metastasis after curative surgical resection. ccRCC in the metastatic form is always associated with high mortality [3,4]. Therefore, given the substantial incidence and mortality of ccRCC, there is an urgent need to develop more efficient prognostic models.
Copper is an indispensable trace element involved in various biological processes. Recent studies showed that the copper levels of cancer patients are significantly elevated both in serum and tumor tissues compared to healthy counterparts [5][6][7]. While dysregulation of copper homeostasis may trigger cytotoxicity, alterations in intracellular copper levels may influence the development and progression of cancer [8]. Based on this mechanism, copper ionophores (disulfiram, dithiocarbamates, elesclomol, etc.) and copper chelators (trientine, tetrathiomolybdate, etc.) have been applied in anticancer treatment [9][10][11][12]. Recently, attention has been brought to a novel cell death pathway termed cuproptosis, and it has been proven that copper binds directly to the lipoylated components of the tricarboxylic acid (TCA) cycle, leading to toxic protein stress and, ultimately, cell death [13]. ccRCC is generally accompanied by a reprogramming of the tricarboxylic acid (TCA) cycle, downregulating the energy production through a TCA cycle and enabling tumor cells to survive in conditions of nutrient depletion and hypoxia and escape from the immune system [14][15][16]. Several genes involved in copper-induced cell death were identified, which may offer novel strategies to predict the prognosis of ccRCC patients.
In the present study, we intended to comprehensively investigate the molecular alterations and clinical relevance of cuproptosis-related genes (CRGs) in ccRCC. Our analysis highlights the importance of CRGs in ccRCC development and lays a foundation for the therapeutic application of cuproptosis regulators in ccRCC.

Multiomics Data Source and Preprocessing
We obtained data of cuproptosis-related genes and clinical information of ccRCC patients from the The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer. gov//, accessed on 15 March 2022). A total of 524 ccRCC patients and 72 adjunct nontumor samples were involved in the present study. Gene expression was measured with STAR as raw read counts, which subsequently transformed into transcripts per million (TPM). All gene features (e.g., chromosome positions, gene types, Ensembl IDs and official symbols) were annotated by the GENCODE project (v22) [17]. Clinical covariates, including the overall survival outcome, age, gender, tumor stage and histological grade, were attained from the previous related resource [18]. Only ccRCC patients with survival information were included in this study.
Somatic mutation data of ccRCC from whole exome/genome sequencing (WXS/WGS) were downloaded from the GDC TCGA-ccRCC project on the UCSC Xena server [19]. The MuTect2 algorithm, which assigns higher levels of confidence to somatic variants [20], was employed to identify additional germline mutations. Oncoplot was drawn according to the descending order of mutations using the R package "maftools" [21].
Digital focal-level copy number variation (CNV) values were calculated from tumor aliquots using a "masked copy number fragment" file by GISTIC2 [22] at the item level and then cut by a noise threshold of 0.3. A Cleveland dot plot was drawn to visualize the frequency of CNV by the R package "ggpubr".

Differential Expression Analysis and Validation
We investigated the differential expression levels of CRGs between tumor and normal samples. Under the criterion of |log 2 (fold change) | > 2/3 and p < 0.05, we considered it as statistically significant.
For validation, we collected 2 datasets (GSE40435 [23] and GSE53757 [24]) of 173 ccRCC samples from Gene Expression Omnibus (GEO). These data were generated using the platform of Affymetrix Human Genome U133 Plus 2.0 Array (GPL570). Box plots were adopted to compare the expression of CRGs in various datasets using the R package "ggplot2".
After computing the log 2 fold change and 95% confidence intervals, we performed a meta-analysis of the results of differential expression to improve the statistical power of our study. We used the Q test (I 2 statistics) as the assessment of the heterogeneity between multiple datasets. If there was no obvious heterogeneity (I 2 < 50%, p > 0.05), a fixed effects model was chosen. Otherwise, we selected a random effects model. A forest plot was utilized to show fold change and related 95% CI of CRGs using the R package "forestplot".

Gene Network and Enrichment Analysis of CRGs
To analyze the potential interactions of these genes, we performed a gene network analysis with the GENEMANIA website [25]. Furthermore, we implemented pathway enrichment analysis of CRGs with the Metascape [26] website. Gene Ontology (GO) as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) were used as references, and enrichment analysis was performed by the R package "clusterProfiler" [27]. We applied the Benjamini−Hochberg method for the multiple correction, and a false discovery rate (FDR) < 0.05 was considered to be of significance.

Construct Prognostic Signature of Cuproptosis-Related Genes with Penalized Regression
To quantify cuproptosis-related genes at the individual level, we developed a signature based on iterative sure independence screening (SIS) [28]. We used SIS and least absolute shrinkage and selection operator (LASSO)−penalized Cox regression to screen for CRGs associated with survival using the R package "SIS" [29].
We then calculated the risk score using the regression coefficients of the identified prognostic signature of CRGs for OS and PFS, respectively. Subsequently, we classified patients into the high-risk and low-risk groups according to the median value of risk scores. The Kaplan−Meier survival curve was plotted to compare the OS or PFS between high-risk and low-risk groups by the R package "ggsurvplot". Utilizing the signature, we also computed the 1-year survival, 3-year survival and 5-year survival based on the nearest neighbor nstimation (NNE) method [30]. Receiver operating characteristic (ROC) curves were computed for presenting the prediction ability using the R package "survivalROC". In addition, we investigated potential differences among subgroups stratified by age, gender and tumor stage. Enhanced regression nomograms of CRG scores and other clinical covariates of ccRCC patients were constructed by the R package "regplot". The calibration curves of the CRGs' scores and other clinical covariates of ccRCC patients were estimated using 1000 bootstrapping to determine bias-corrected estimates of predicted versus observed values, which was analyzed using the R package "rms".

Analysis of Correlation with Immune Infiltration
A Tumor Immune Estimation Resource (TIMER; cistrome.shinyapps.io/timer) [31] was used to investigate the relationship between the expression of CRGs and the abundance of six immune cells (CD4 + T cells, CD8 + T cells, B cells, neutrophils, dendritic cells and macrophages).
We also examined three important immune checkpoints (ICKs), including PD-1, PD-L1 and TIM-3 since the expression level of immune checkpoint-related genes is related to the treatment response of immune checkpoint inhibitors [32].
The Pearson correlation analysis was used to examine the association between TME or ICKs and CRGs.

Statistical Analysis
First, we performed a descriptive statistical analysis of ccRCC patients in TCGA. Continuous variables were described as mean ± standard deviation, and categorical variables were described by frequency and proportion. A Kruskal−Wallis rank sum test [33] was applied to examine the difference of CRG expression in various classifications of pathologic stage and histological grade of ccRCC patients.
We utilized the Benjamini−Hochberg method for multiple correction, and a FDR < 0.05 was the standard of statistical significance [34]. All statistical analyses were performed using R version 4.1.1 (The R Foundation). p values were two-sided, and we considered a level of p value less than 0.05 to be statistically significant.

Functional Enrichment and Protein-Protein Interaction Analysis of CRGs
To demonstrate the biological functions of CRGs, relevant pathways were analyzed by GO and KEGG databases. The biological processes of the 10 CRGs mainly involved in the GO analysis were the acetyl-CoA biosynthetic process from pyruvate, acetyl-CoA biosynthetic process, tricarboxylic acid cycle, acetyl-CoA metabolic process, thioester biosynthetic process, mitochondrial matrix, oxidoreductase complex, mitochondrial protein-containing complex, dihydrolipoyl dehydrogenase complex, tricarboxylic acid cycle enzyme complex, oxidoreductase activity, iron−sulfur cluster binding and metal cluster binding ( Figure 2A). Additionally, in the KEGG pathway enrichment analysis, the 10 CRGs were largely related to the TCA cycle, pyruvate metabolism, glycolysis/gluconeogenesis, carbon metabolism, lipoic acid metabolism, central carbon metabolism in cancer, biosynthesis of cofactors, glucagon signaling pathway, HIF-1 signaling pathway and D-Amino acid metabolism ( Figure 2B). A Protein-Protein Interaction (PPI) analysis was performed to explore the interactions of CRGs, showing that DLD, PDHB, DLAT and PDHA1 were hub genes ( Figure S1).

Construction of the Prognostic Signature of Cuproptosis-Related Genes in ccRCC
We further evaluated the association between the expression of CRGs and the prognosis in ccRCC. We found that all these genes except GLS were highly correlated with the overall survival (OS) in the univariate Cox proportional hazard regression model after the adjustment of age, gender, race and pathologic stage (Table S2) (Table S2).
Then, we separately constructed the prognostic signature of CRGs for OS and PFS in ccRCC utilizing SIS. For OS outcomes of ccRCC patients, three genes were selected to construct the prognostic score using their regression coefficients: Score os = −0.44 × FDX1−0.28 × DLAT + 0.23 × CDKN2A. A higher risk score of CRGs' signature was significantly associated with poor OS (HR = 2.72(2.01-3.68), log-rank p = 1.76 × 10 −7 , Figure 3A,B). We also applied a weighted risk score incorporating all related genes to estimate 1-, 3-and 5-year OS ( Figure 3C,F). The prediction accuracy evaluated by AUCs was reported to be 0.652, 0.633 and 0.658 in the 1-year, 3-year and 5-year ROC curves, respectively. For PFS, we built a risk score of cuproptosis-related gene signature using the following formula: risk score = ∑ n i=1 Coe f (Gene) × Expr(Gene) where Coef (Gene) was the coefficient of genes (i.e., FDX1, DLAT and CDKN2A) correlated with PFS, and Expr (Gene) was the expression signature of corresponding genes. The same analyses were conducted for the PFS outcome. Patients with a higher risk score had significantly shorter PFS (HR = 2.83 (2.08-3.85), log-rank p = 3.66 × 10 −7 , Figure 3G,H), and the reported AUCs of 1-year, 3-year and 5-year ROC curves for the prediction of PFS were 0.622, 0.634 and 0.682, respectively ( Figure 3G,I). In different subgroups of age, gender and pathologic stage, the CRG risk score also had a good performance (Table S3). In short, our constructed individual-level cuproptosis-related risk signature showed a significant association with survival of ccRCC.

Nomogram Development and Validation for ccRCC
To facilitate the clinical application of the prediction model, we integrated clinical information and gene features of patients from TCGA and performed the multivariable Cox regression model to develop the nomogram. Discrimination and calibration methods were applied in both OS and PFS outcomes (Figure 4). The c-index was calculated to be 0.77 for OS and 0.824 for PFS, reflecting a relatively excellent predictive performance of the nomogram. Meanwhile, calibration plots demonstrated favorable concordance between the predicted OS or PFS and the observed OS or PFS at 1, 3 and 5 years of survival ( Figure 4C,F).

Correlation between Expression of CRGs and Immune Infiltration Levels in ccRCC
It is uncertain whether CRGs would influence immune cell recruitment in the tumor microenvironment and, therefore, affect the prognosis of ccRCC. Thus, we performed an analysis to examine the relationships between CDKN2A, DLAT, FDX1 and LIAS and immune infiltration in ccRCC. The expression level of CDKN2A was positively associated with the immune infiltration level of CD8 + T cells (p = 2.89 × 10 −2 ) and negatively correlated with macrophages (p = 2.89 × 10 −2 ) ( Figure 6A). The DLAT expression level was positively correlated with the immune infiltration level of B cells (p = 1.40 × 10 −6 ), macrophages (p = 2.93E × 10 −13 ), neutrophils (p = 1.89 × 10 −6 ) and dendritic cells (p = 1.02 × 10 −4 ) ( Figure 6B). The FDX1 expression level was positively associated with the abundance of B cells (p = 2.33 × 10 −3 ) and macrophages (p = 1.73 × 10 −2 ) ( Figure 6C). Figure 6D presents the positive association between the LIAS expression and the abundance of CD8 + T cells (p = 3.86 × 10 −2 ), macrophages (p = 1.12 × 10 −5 ) and neutrophils (p = 1.15 × 10 −3 ). Our results also showed that CDKN2A expression in ccRCC had a positive correlation with PDCD1 expression levels (r = 0.24, p = 2.14 × 10 −8 ).   Figure 8, the expression levels of four genes (CDKN2A, DLAT, FDX1 and LIAS) varied in different pathologic stages and histological grades of ccRCC, with the exception of DLD, which did not show any statistical difference across histological grades (p = 0.063, Figure 8B). Specifically, there was a downward trend of FDX1 and LIAS expression and an upward trend of CDKN2A expression regardless of tumor stage or histological grade. These results suggested that the expression of CRGs may be correlated with disease grade and the presence of necrosis of ccRCC.

Discussion
In the present study, we explored the expression signature of 10 CRGs in ccRCC tissues and examined their relationships with OS and PFS. A novel cuproptosis-related prognostic score was constructed for the first time. In addition, functional analyses exhibited that pathways related to TCA cycle were enriched, and the CRGs were also proven to be associated with the grading and staging of ccRCC.
To the best of our knowledge, there have been no previous studies examining the correlations between CRGs and the development of ccRCC. Surprisingly, most cuproptosisrelated genes were differentially expressed between tumor and normal tissues, and all these genes were significantly associated with OS and PFS, suggesting a potential role of cuproptosis in the prognosis of ccRCC and the predictive value of this score in the prediction of ccRCC survivorship.
The prognostic score constructed in this study consisted of four cuproptosis-related genes (CDKN2A, DLAT, FDX1 and LIAS). FDX1 encodes a small iron−sulfur protein and was involved in the reduction of Cu 2+ to Cu 1+ . In addition, FDX1 played the role of upstream regulator in the process of protein lipoylation in the TCA cycle and key regulator of cuproptosis [13,35,36]. Dihydrolipoamide S-acetyltransferase (DLAT) was one of the components of the pyruvate dehydrogenase (PDH) complex. The oligomerization of DLAT was due to the integration of copper and lipoylated proteins in the TCA cycle [13]. Lipoic acid synthase (LIAS) encoded components of the lipoic acid pathway and synthesized a potent antioxidant termed α-Lipoic acid (LA) in mitochondria [35]. The CDKN2A expression functioned in the cell cycle control and was strongly correlated with the origin of a variety of tumors [36][37][38][39]. Furthermore, CDKN2A was shown to be absent in 76% of metastatic ccRCC samples, according to a meta-analysis [39].
Mutations leading to the overload of copper could trigger severe consequences. Nevertheless, it is feasible to manage intracellular copper levels within a certain range to selectively kill tumor cells [6]. Cuproptosis, an unconventional mechanism of cell death concerning the protein lipoylation in TCA cycle, might indicate novel insights to exploit copper toxicity to treat tumors [40]. Moreover, there have been some new findings about the pathologies of the urinary tract. Evidence has shown that there is an association between microbiota and the renal function, which is involved in the progression of several kidney diseases [41]. It has also been found that copper may contribute to the attenuation of bacterial colonization in the urethra, which may shed light on the potential therapeutic target [42]. Of interest, plant extracts contributing to the regulation and synthetization of copper-related metabolism may present as an alternative nonpharmacological intervention (NPI) strategy for cancer prevention and treatment [43,44]. For example, oregano extract could be one of the potential anticancer NPIs by invoking cell death through the mitochondrial and DNA damage pathways related to cuproptosis [44]. Additionally, plant extracts, such as the extract of the Camellia sinensis leaf, may be used to synthesize copper-related preparation and, thus, has the potential to play a part in the cancer treatment [45].
Our study has multiple advantages. Firstly, this study is the first to develop a prognostic model on the basis of CRGs. Cell death, an intense area of tumor research, has been proven fundamental to cancer origin and development [46]. Cuproptosis is a novel type of cell death that differs from any known mechanisms of cell death and is reliant on mitochondrial respiration. Such an unusual mechanism may bring up new solutions for the treatment of cancer. Moreover, compared to majority of the prognostic models which mainly focused on the overall survival, our study also took progression-free survival into consideration. PFS gain may indicate the reduction in tumor burden and symptom relief, especially for advanced cancer patients [47][48][49]. Additionally, we employed SIS to identify CRGs and adjusted covariates in both screening and testing procedures to make the results more robust.
There were several limitations in our study. First, although the validation of the expression of screened prognostic genes in TCGA and GEO demonstrated, to some extent, that three genes (CDKN2A, DLAT and FDX1) showed a robust prognostic effect, regrettably, no dataset with a sufficiently large sample size (n > 50) and clinical prognostic information was available for further validation, which is urgently warranted in future research. Second, although the prognostic score focusing on CRGs' expression signature showed a good performance in the prediction of ccRCC survivorship, some other significant genes with predictive values were not considered in this study. Third, given that the prognostic signature was built and validated by exploiting data from public databases, further biological evidence is needed apart from the statistical evidence we offered.
In summary, this study systematically analyzed the landscape of molecular alterations and interactive genes of cuproptosis in ccRCC. Our study demonstrated that these cuproptosis-related genes may play a crucial role in ccRCC outcomes. The prognostic risk score based on the expression signature of CRGs showed a good performance for the prediction of OS and PFS of ccRCC patients and was significantly associated with immune infiltration levels and PD1 expression. Our results would also provide novel insights in developing pharmacological and nonpharmacological therapeutic strategies related to cuproptosis for cancer prevention and treatment.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes13050851/s1, Figure S1: Protein-Protein Interaction (PPI) of cuproptosis-related genes (CRGs); Table S1: Results of differential expression of cuproptosis-related genes (CRGs) in various datasets; Table S2: Association results for cuproptosis-related genes derived from univariate and multivariate Cox proportional hazards model. Table S3: Association results of the cuproptosis-related gene signature of subgroup (age, gender and pathologic stage) derived from Cox proportional hazards model.