Comprehensive analysis of cuproptosis-related genes in prognosis, tumor microenvironment infiltration, and immunotherapy response in gastric cancer

Cuproptosis is the most recently identified copper-dependent cell death form that influences tricarboxylic acid (TCA) cycle. However, the relationship between cuproptosis and clinical prognosis, tumor microenvironment infiltration (TME), and response to immunotherapy remains unclear. Single-sample gene-set enrichment analysis (ssGSEA) was employed to construct cuproptosisScore (cpS) and 1378 gastric cancer (GC) patients from five independent public datasets were classified into high- or low-cpS groups according to the median of cpS. Then the impacts of cuproptosis on tumor microenvironment infiltration (TME), biological function, response to immunotherapy, and clinical prognosis of GC were evaluated. RiskScore and nomogram were constructed using Lasso Cox regression algorithm to validate its predictive capability in GC patients. Compared to patients with high cpS, patients with low cpS exhibited poorer prognosis, higher TNM stage, and stronger stromal activation. Meanwhile, the analysis of response to immunotherapy confirmed patients with high cpS could better benefit from immunotherapy and had a better susceptibility to chemotherapeutic drugs. Then, 9 prognosis-related signatures were collected based on differentially expressed genes (DEGs) of cpS groups. Finally, a riskScore model was constructed using the multivariate Cox (multi-Cox) regression coefficients of prognosis-related signatures and had an excellent capability of predicting 1-, 3-, and 5-year survival in GC patients. This study revealed the role of curproptosis in TME, response to immunotherapy, and clinical prognosis in GC, which highlighted the significant clinical implications of curproptosis and provided novel ideas for the therapeutic application of cuproptosis in GC.


Introduction
Gastric cancer (GC) is the fifth most commonly diagnosed cancer and the third leading cause of cancer-related deaths in the world (more than 1 million cases resulted in more than 768,000 deaths in 2020 worldwide) (Ajani et al. 2022). Currently, multidisciplinary treatments, such as chemotherapeutics and target therapeutics, are the only choice for patients who cannot be surgically treated, but they were frequently reported with treatment failure and side effects (Ajani et al. 2022). There is a rapid development in the field of immunotherapy in recent years, including the use of checkpoint inhibitors targeting programmed cell death ligand 1 (PD-L1), cytotoxic T lymphocyte antigen 4 (CTLA4), or adoptive T-cell therapy (ACT) (Chen and Mellman 2013). However, only about 30% of patients benefit from immunotherapy in most tumors (Tang et al. 2020). Therefore, it is an important direction to explore novel GC therapeutic approaches.
As a vital micronutrient, copper is involved in fundamental and conserved life processes throughout all forms of life (Denoyer et al. 2015). It serves as a catalytic cofactor for essential enzymes involved in the regulation of energy conversion, iron collection, and intracellular oxidative metabolism (Kim et al. 2008). The previous study also reported copper plays a key role in numerous biochemical processes, including synthesis of hemoglobin and redox processes (Mroczek-Sosnowska et al. 2015). It is necessary to maintain a strict level of copper concentration in cells. Too low or too high copper concentration can cause significant influence on cellular physiology, not only in bacteria but also in human cells (Lutsenko 2010). Intracellular copper overload induced by elesclomol can inhibit colorectal cancer with high efficiency (Gao et al. 2021). Moreover, increased copper concentration was reported to be associated with the inhibition of endometrial cancer and prostate cancer (Atakul et al. 2020). Thus, copper has a potentiality to be a new therapeutic target used to treatment of tumors by increasing intracellular copper accumulation (Ge et al. 2021). Of note, the previous study of Tsvetkov et al. indicated that intracellular copper accumulation could induce a novel cell death, termed cuproptosis, and it occurs by the direct binding of copper and lipoylated components of the tricarboxylic acid (TCA) cycle. Several cuproptosis-related genes (CRGs) were identified by Tsvetkov et al., which provides a new strategy for treatment and prognosis prediction of tumors (Tsvetkov et al. 2022). Considering there have been no previous research about the effect of cuproptosis to GC in TME and prognosis, the correlation between cuproptosis and GC was explored in this study.
Based on CRGs, we developed a cuproptosisScore (cpS) to explore the differences of tumor microenvironment infiltration (TME), biological function, response of immunotherapy, and clinical prognosis in different groups of GC. Then, prognosis-related signatures were identified to further validate the effect of cuproptosis for predicting prognosis. We suggested that this study will lay a foundation for the therapeutic application of cuproptosis in GC.

Collection and preprocessing of GC dataset
A flowchart of our study is shown in Fig. S1. In the study, five independent cohorts from patients with gastric cancer were collected: GSE62254/ACRG, GSE57303, GSE84437, GSE15459, and TCGA-STAD (The Cancer Genome Atlas-Stomach Adenocarcinoma) (Tab. S1). The raw data from the microarray data sets and the relevant clinicopathological data were downloaded from the Gene-Expression Omnibus (GEO; https:// www. ncbi. nlm. nih. gov/ geo/). The upper quartile fragments per kilobase of transcript per million mapped reads (FPKM-UQ) and the relevant prognostic data of STAD were downloaded from the UCSC Xena browser (https:// xenab rowser. net/ datap ages/). In addition, a normal gastric cohort of 174 samples was downloaded from Genotype-Tissue Expression Project (GTEx) (https:// www. gtexp ortal. org/). Then, we converted the FPKM value to the transcripts per kilobase million (TPM) value of TCGA-STAD by "data.table", "tibble", "dplyr", and "tidyr" R packages, and were believed to be identical to those from microarrays (Trapnell et al. 2010). The microarray data adopted a robust multiarray averaging method with the "affy" and "simpleaffy" R packages to perform background adjustment and quantile normalization. Batch effects from non-biological technical biases were corrected using the "ComBat" algorithm of "sva" R package. Finally, a total of 1378 GC patients were gathered in the subsequent analyses.

Single-sample gene-set enrichment analysis (ssGSEA) for CRGs
According to the previous study of Tsvetkov et al.,12 CRGs existing in all cohorts were selected in the subsequent analyses, including PDHA1, FDX1, SLC31A1, ATP7B, ATP7A, DBT, LIPT1, DLD, PDHB, DLAT, GCSH and LIAS (Tsvetkov et al. 2022). Single-sample gene-set enrichment analysis (ssGSEA) was applied to quantify CRGs in each sample of five GC cohorts with the "GSVA" R package and cuproptosisScore was constructed through ssGSEA score. The patients were divided into low-and high-score groups, we termed as low-cuproptosisScore and high-cuproptosisScore groups according to the median value of the ssGSEA score. Detailed information of two groups is shown in Tab. S2.

Prognosis and other clinical characteristics in cpS groups
To evaluate the clinical value of the two groups identified by ssGSEA, we performed the following analysis. Correlation between prognosis and CRGs was assessed by univariate Cox regression (uni-Cox) analysis. The differences in overall survival (OS) among different groups were assessed using Kaplan-Meier (KM) analysis by the "survival" and "survminer" R packages.
Considering comprehensive clinicopathological information in ACRG cohort (GSE62254), we explored the relationships between cpS groups and related clinical characteristics in ACRG cohort, including age, sex, survival status, ACRG subtype, TNM stage, and histology subtype.

Analysis of correlation between cpS groups and TME infiltration
We evaluated the immune, stromal, and ESTIMATE score of tumor samples by ESTIMATE algorithm with "ESTIMATE" R package, and quantified levels of 29 immune infiltration cells from the gene set of the previous study of Charoentong et al. in each sample by ssGSEA with "GSVA" R package (Charoentong et al. 2017).

Response of immunotherapy and drug susceptibility analysis
To estimate the effect of cpS predicting the response of immunotherapy, we performed a systematical collected for immunotherapy-related cohorts. Anti-PD-L1 cohort (GSE78220) and ACT cohort (GSE100797) were downloaded from GEO, and anti-CTLA4 cohort of previous study was downloaded from http:// tide. dfci. harva rd. edu (Yan and Richmond 2021). Then, stratified analysis of the response of immunotherapy and KM survival analysis were performed.
To explore differences in the therapeutic response of chemotherapeutic drugs of GC, we calculated the semi-inhibitory concentration (IC50) values of chemotherapeutic drugs commonly used to treat GC using the "pRRophetic" R package.

Differentially expressed genes (DEGs) associated with the cuproptosis phenotype
To further identify genes associated with the cuproptosis phenotype, the empirical Bayesian approach of "limma" R package was used to determine DEGs between two cpS groups based on ACRG cohort (Ritchie et al. 2015). The significance criteria for determining DEGs was set as adjusted p value < 0.05 and |log2(fold change) |> 0.6. Then, the patients were classified into high-and low-geneScore groups according to the median value of the ssGSEA score.

Construction of cuproptosis-related signatures and riskScore model
To verify prognostic value of cuproptosis in GC, the risk score was calculated to quantify the cpS groups of the individual tumors. The DEGs were subjected to uni-Cox analysis based on the clinical data of all GC cohorts to identify cuproptosis-related gene signatures. Then, 300 patients of ACRG cohort were randomly categorized into training and testing sets at a ratio of 6:4. Based on training set, prognosis-related gene signatures were identified using the Lasso Cox regression algorithm with the "glmnet" R package, and a p value was set at < 0.05 as well as run for 1000 cycles. A random stimulation was set at 1000 times to avoid overfitting in each cycle (Zhao et al. 2021). Finally, a riskScore model was constructed, and the risk score was calculated with the following formula: risk score = Σ(Expi * coefi). where Coefi and Expi denote the multivariate Cox regression (multi-Cox) coefficient and expression of each gene, respectively. Based on the cut-off value of risk score calculated by the "survminer" R package, a total of 180 patients in the training set of ACRG cohort were divided into low-riskScore and high-riskScore groups. KM survival analysis and the generation of receiver operating characteristic (ROC) curves were used to estimate the effect of predicting prognosis. Finally, the testing set and other GC cohorts (GSE15459, GSE57303, GSE84437, STAD) were divided into low-and high-riskScore groups to evaluate model stability.

Development and estimation of a nomogram
Potential differences between prognosis with riskScore and clinical information, such as age, sex, survival status, ACRG subtype, TNM stage, and histology subtype, were investigated by uni-Cox and multi-Cox analysis. Enhanced regression nomograms of riskScore and clinical covariates of ACRG samples were constructed by the "regplot" R package. The calibration curves based on 1000 bootstrapping were constructed using "rms" R package to determine biascorrected estimates of predicted versus observed values.

Analysis of functional and pathway enrichment
Gene set variation analysis (GSVA), an unsupervised and non-parametric algorithm, is commonly used to estimate the variation in pathway and biological process activity in expression dataset. (Hänzelmann et al. 2013). We performed GSVA analysis using "GSVA" R packages to explore the difference on biological process between cpS groups. The gene sets of "h.all.v7.5.1.symbols" downloaded from MSigDB database was selected to running GSVA analysis. |log2(fold change) |> 0.3 and adjusted p < 0.05 was considered as statistically significance.

Statistical analysis
All analysis was conducted using the R V.4.1.2 software and GSEA software. Wilcoxon or Kruskal-Wallis tests were used for comparison between two or multiple groups, respectively. When the data was non-normally distributed, the Mann-Whitney test or Dunn's test were used. KM survival curves were plotted using the "survival" R package and analyzed by the log-rank test. The correlation between cpS or other score and clinicopathologic features was analyzed using the χ 2 test. Statistical significance criteria were set at p < 0.05.

Landscape of genetic variation of CRGs in gastric cancer
A total of 12 CRGs were collected in this study (Tsvetkov et al. 2022). First, copy number variations (CNV) and somatic mutations of these genes in the TCGA-STAD cohort were evaluated. 67 (15.47%) were found had mutations in 433 samples, of which ATP7B had the highest mutation frequency, followed by ATP7A, DLAT, and DLD. CNV alteration frequency showed a prevalent CNV alteration in 12 genes (Fig. 1A). ATP7B and SLC31A1 had significant CNV increase, while DLAT, FDX1, DBT, and PDHB had significant CNV deletion (Fig. 1B). The location of CNV alteration of 12 genes on chromosomes is shown in Fig. 1C. Then, we further analyzed the relationship of expression level among 12 CRGs and noticed that there was a significant positive correlation between these genes ( Fig. 1D). To further explore whether these genes can distinguish between normal and disease samples in GC, we performed principal component analysis (PCA) and evaluated the relationship between the expression of 12 genes in tumor and normal samples. The results showed that these genes had a good power in distinguishing GC from normal tissue (Fig. 1E). The boxplot showed that the expression level of 12 CRGs was markedly upregulated in GC tumors than that in normal tissues (Fig. 1F). Taken together, these analyses suggested that the high heterogeneity and expression imbalance of CRGs played a crucial role in the GC occurrence and progression.

Construction and clinical characteristics of cpS groups in GC
To identify the expression pattern of CRGs involved in tumorigenesis, 1378 patients of five GC cohorts (GSE62254/ACRG, GSE57303, GSE84437, GSE15459, TCGA-STAD) were enrolled into a meta-cohort for further analysis. 1378 patients were divided into high-and low-cpS group according to the median value of the ssGSEA score. PCA showed that the distribution of the two groups can be markedly distinguished ( Fig. 2A). Compared to low-cuprotosisScore group (LCSG), CRGs demonstrated markedly higher expression in high-cuproptosisScore group (HCSG) ( Fig. 2B; Tab. S2). The distribution of cpS, the survival status, and survival time of GC patients in two cpS groups are shown in Fig. 2C. KM survival analysis further confirmed that the patients with higher cpS showed better OS (Fig. 2D). The results of uni-Cox analysis revealed the relationship of negative correlation between cpS with prognosis had statistical significance and most of CRGs were associated with worse prognosis (Fig. 2E).
Comprehensive clinicopathological data of ACRG cohort (GSE62254) can be used to explore the relationship between cpS and prognosis. Similar to meta-cohort, GC patients of ACRG cohort were also divided into HCSG and LCSG through ssGSEA analysis. (Fig. S2A, B). Then, GC CRGs. E PCA for the expression profiles of 12 CRGs to distinguish tumors from normal samples in a pooled cohort consisted of TCGA-STAD and GTEx. *p < 0.05, **p < 0.01, ***p < 0.001. F The expression of 12 CRGs between normal and tumor samples. CRGs, cuproptosis-related genes; GC, gastric cancer; CNV, copy number variations; PCA, principal component analysis; TCGA-STAD, The Cancer Genome Atlas-Stomach Adenocarcinoma. GTEx, Genotype-Tissue Expression Project patients were classified by age, gender, status, TMN stage, tumor stage, histological subtype, and ACRG subtype. The cpS was not significantly associated with age (p = 0.061) and gender (p = 0.29). Samples in LCSG showed markedly higher TNM (p < 0.05) and tumor stage (p = 0.022) than in HCSP. Of note, Samples with EMT molecular subtype and diffuse histological subtype were characterized by LCSG, while MSI subtype and intestinal histological subtype were characterized by the HCSG (Fig. 2F; Fig.  S2C). In GC, the EMT molecular subtype and diffuse histological type were markedly associated with poorer survival, while MSI and intestinal histological type were linked to a better clinical outcome (Menju and Date 2021). KM survival analysis also showed a poorer prognosis in LCSG (Fig. S2D). Therefore, a poorer prognosis of LCSG was correlated with stromal activation, rapid progression, and high malignancy.

Functional annotation and characteristics of the TME in cpS groups
To further explore the relationship between the biological behaviors and prognosis of cpS groups, GSVA enrichment analysis was performed. Of note, LCSG was markedly enriched in not only stromal activation pathways, such as angiogenesis, hypoxia, and EMT, but also immune activation pathways, such as inflammatory response, IL-6/JAK/ STAT3 signaling, and IL2/STAT5 signaling ( Fig. 3A; Tab. S3). Then subsequent analyses of ESTIMATE also indicated stromal, immune, and ESTIMATE score of the LCSG were statistically higher than HCSG (Fig. 3B). However, LCSG did not show a matching survival advantage (Fig. 2D). The previous study indicated cuproptosis was distinct from all other known mechanisms of regulated cell death occurs and the relationship between cuproptosis and TME was not clear (Tsvetkov et al. 2022). In addition, other studies suggested that tumors with immune excluded phenotype also showed the presence of abundant immune cells, while these immune cells remained in the stroma surrounding tumor Activated stroma in TME suppressed the effect of T-cell (Chen and Mellman 2017). To further verify this assumption, we analyzed the relationship between cpS and other immune-related behaviors (Tab. S4) (Mariathasan et al. 2018). It also showed that the stromal activity of LCSG was significantly increased, such as EMT, transforming growth factor-β (TGFβ), and angiogenic pathway activation (Fig.  S3). We then found the fold change of stromal activation pathways was significantly higher than immune activation Fig. 3 Functional annotation and characteristics of the TME in two cpS groups. A GSVA showing the activation status of biological pathways in two cpS groups. The heatmap was used to visualize these biological pathways, and yellow represented activated pathways and blue represented inhibited pathways. The gastric cancer cohorts were used as sample annotations. B The comparison of stromal, immune, and ESTIMATE score between two cpS groups. *p < 0.05, **p < 0.01, ***p < 0.001. C The fold change of stromal activa-tion pathways and immune activation pathways in LCSG vs HCSG. Red represented stromal activation pathways and blue represented immune activation pathways. D The difference of infiltration levels of 29 immune cells between two cpS groups. TME, tumor microenvironment infiltration; GSVA, gene-set variation analysis; cpS, cuprop-tosisScore; HCSG, high-cuprotosisScore group; LCSG, low-cuproto-sisScore group 1 3 Fig. 4 The comparison of response to immunotherapy and drug susceptibility between two cpS groups. A-D The upper part showed the difference of proportion of patients who respond to immunotherapy in two cpS groups. The lower part showed the survival analysis of two cpS groups. GSE78220 is anti-PD-L1 cohort. GSE100797 is ACT cohort. Nathanson2017 pre and Nathanson2017 post are anti-CTLA4 cohorts. E The prediction of drug susceptibility in cpS groups from five GC cohorts (GSE62254, GSE57303, GSE84437, GSE15459, and TCGA-STAD). ACT, adoptive T-cell therapy; cpS, cuproptosisScore; GC, gastric cancer pathways in LCSG (Fig. 3C). Finally, TME cell infiltration also indicated abundant immune cells presented in LCSG (Fig. 3D). The above results proved that the effect of TME cells in LCSG was suppressed due to the stromal activation.

Prediction the immunotherapy response of GC in cpS groups
Then, the cuproptosisScore was used to predict the response of immunotherapy of GC to explore the sensitivity and clinical prognosis of targeted medicine for cuproptosis. In the anti-PD-L1 cohort (GSE78220), ACT cohort, and anti-CTLA4 cohort, patients with higher cpS had statistically longer progression-free survival and better therapeutic response than those with lower cpS (Fig. 4A-D). Meanwhile, drug sensitivity analysis indicated HCSG had a higher IC50 of 10 chemotherapeutic drugs, such as elesclomol, midostaurin, and nilotibin ( Fig. 4E; Fig. S4). In summary, our study indicated that cpS could effectively predict the response of therapy of GC in the common frontier field of immunotherapy.

Construction and validation of the DEGs
To further explore the potential biological behavior of cpS groups, 1125 cuproptosis phenotype-related DEGs were selected using "limma" R packages based on ACRG cohort (Tab. S5). 300 GC samples of ACRG cohort were classified into high-geneScore group (HGSG) and low-geneScore group (LGSG) through the median value of the ssGSEA score. To our surprise, the expression of these DEGs exhibited nearly complete opposite clinicopathological trends compared with the cpS groups. Patients with alive status, MSI subtype, intestinal histological subtype, and lower tumor TNM stage were mainly concentrated in LGSG, while patients with dead status, EMT subtype, diffuse histological subtype, and higher TNM stage were mainly concentrated in HGSG. Consistent with above results, we found most patients with high geneScore had low cpS and vice versa (Fig. 5A). The expression level of 12 CRGs in LGSG was generally upregulated (Fig. 5B). Above findings suggested that these DEGs had a negative correlation with CRGs in the regulation of biological function.
We then investigated clinical prognosis and TME immune regulation of the cuproptosis-related DEGs. KM curves showed that patients in LGSG had a better prognosis than in HGSG (Fig. 5C). Consistent with GSVA results of cpS groups (Fig. 3A), GSEA analysis found pathways considered T-cell suppressive and immune activation pathways were associated with high geneScore (Tauriello et al. 2018) (Fig. 5D; Fig. S5A; Tab. S6). ESTIMATE analysis indicated that stromal, immune, and ESTIMATE score of the HGSG were statistically higher than LGSG (Fig. S6A). Examination of the known signatures depicted the fold change of stromal activation pathways was significantly higher than immune activation pathways in LGSG, which was consistent with the above results (Fig. 5E) (Charoentong et al. 2017). In addition, the collected signatures related to cytokine and chemokine were extracted from the previous study, including TGFβ/EMT pathway and immune activation (Zeng et al. 2019). It was found that the expression of all genes was significantly upregulated in HGSG (Fig. S6B-C). Finally, the alluvial diagram exhibited that HGSG with the EMT of ACRG molecular subtype was correlated with lower cpS and was associated with a poorer prognosis ( Fig. 5F; Tab. S7). The above results showed again that cuproptosis played a non-negligible regulation role in the functional regulation of TME.

Construction of the riskScore model of prognosis-related gene signatures
We further evaluated the association between cuproptosis and the prognosis in GC based on ACRG cohort. 731 cuproptosis-related signatures were found to affect prognosis based on uni-Cox analysis (Tab. S8). To avoid overfitting the prognostic signature, nine prognosis-related gene signatures with minimal lambda were extracted by the Lasso regression (Fig. 6A, B). RiskScore model was constructed using multi-Cox regression coefficients of nine prognosis-related signatures. The risk score was calculated as follows: risksco re = ECE1 × (0.215) + PLAT × (0.164) + FLOT1 × (0.339) + NOX4 × (0.119) + TMCC1 × (0.266) + LOXL4 × (0.162) + P REX2 × (0.143) + PDE12 × (-0.211) + PACRGL × (-0.272). Then, 300 patients of ACRG cohort were classified in high-riskScore group (HRSG) and low-riskScore group (LRSG) using the cut-off values obtained with the survminer R package. It was found that patients with high riskScore showed poorer OS than those with low riskScore (Fig. 6C). KM survival analysis also indicated HRSG had a poorer prognosis than LRSG (Fig. 6D). In addition, the 1-, 3-, and 5-year survival rates of riskScore were represented by area under the curve (AUC) values of 0.638, 0.683, and 0.681, respectively (Fig. 6E). To further validate the stability of the riskScore model, we further validated the prognostic value of several external GC cohorts. Our findings indicated that the riskScore model had a stable ability to predict prognosis in all GC cohorts (Fig. 6F-J).
To explore the correlation in riskScore between cpS and geneScore, we performed subsequent analyses. Patients with higher riskScore had lower cpS and higher geneScore (Fig. 6K, L); Tab. S7). According to the results of the TME in cpS groups, HRSG may be more related to stromal activation than immune activation. A 300 patients of ACRG cohort was classified into high-and low-geneScore groups using ssGSEA. The heatmap showed the differences in clinicopathological features and expression levels of DEGs between two geneScore groups. B The expression of 12 CRGs between two geneScore groups. *p < 0.05, **p < 0.01, ***p < 0.001, ns (p > 0.05) C Survival analysis showed OS time of two geneScore groups. D GSEA showing the activation status of biological pathways in two cpS groups. E The fold change of stromal activation pathways and immune activation pathways in LGSG vs HGSG. Red represented stromal activation pathways and blue represented immune activation pathways. F Alluvial diagram showed the changes of ACRG molecular subtypes, cpS, geneScore, and survival outcomes. DEGs, differentially expressed genes; ssGSEA, single-sample gene-set enrichment analysis; OS, overall survival; GSEA, gene-set enrichment analysis; cpS, cuproptosisScore; HGSG, high-geneScore group; LGSG, low-geneScore group

Nomogram development and estimation of riskScore
We further explored the clinical value of riskScore in GC of ACRG cohort. Through uni-Cox and multi-Cox regression analysis, we found that multiple factors including TNM stage, tumor stage, and riskScore were independent prognostic parameters (Fig. 7A, B). Then, a nomogram was developed with the results of multi-Cox regression to predict the 1-, 3-, and 5-year OS rates (Fig. 7C). The calibration plot showed that the nomogram performed well in predicting patient OS according to an ideal model (Fig. 7D).
Furthermore, the calibration plot constructed by using other GC cohorts also indicated an excellent ability to predict the prognosis of the nomogram (Fig. 7E-G).

Discussion
Copper is an indispensable trace element in life processes, imbalance of intracellular copper concentration can occur different effects on biological processes (Ge et al. 2021). Meanwhile, numerous research proved that copper imbalance is associated with the development and metastasis of a Fig. 6 Construction of prognosis-related signatures and riskScore model. A The LASSO coefficient profile of nine prognosis-related signatures. B Cross-validation for turning parameter selection through minimum criteria in the LASSO model. C Distribution of riskScore and survival status of GC patients. D Survival analysis showed OS time of two riskScore groups. E The 1-, 3-, and 5-year ROC curves of ACRG cohort. F-J Validation of the riskScore model in fifth external independent sets and a meta-cohort consisted of five GC cohorts (GSE62254, GSE57303, GSE84437, GSE15459, and TCGA-STAD). K Differences in riskScore among two cpS groups in ACRG cohort. L Differences in riskScore among two geneScore groups in ACRG cohort. GC, gastric cancer; OS, overall survival; ROC, receiver operating characteristic; cpS, cupropto-sisScore variety of tumors (Atakul et al. 2020). The previous study by Tsvetkov et al. suggested a copper-related cell death, termed cuproptosis, causes aggregation of lipid acylated proteins and loss of iron-sulfur (Fe-S) cluster proteins and increases proteotoxic stress through intracellular copper accumulation (Tsvetkov et al. 2022). Taken together, cuproptosis can provide a novel potential for tumor treatment. Thus, we explored the correlation between cuproptosis and GC in this study, and a novel cuproptosisScore about GC was constructed based on 12 CRGs for the first time. We then comprehensively investigated the relationship between cpS and the molecular alterations, TME, response of immunotherapy, and clinical prognosis in GC. Functional analyses revealed that stroma-and immune-related pathways were enriched.
12 CRGs was obtained from the previous study, including FDX1, LIAS, SLC31A1, ATP7B, ATP7A, DBT, LIPT1, DLD, PDHA1, PDHB, DLAT, and GCSH (Tsvetkov et al. 2022). Then, we collected some reports to explain the prognostic role of the gene signature. FDX1 and lipoylated proteins play a key role in cell death induced by copper ionophores. The decreased expression of FDX1 and LIAS inhibit the occurrence of cuproptosis. In addition, FDX1 was found to be associated with the glucose metabolism, fatty acid oxidation, and amino acid metabolism and can impact the prognosis of lung adenocarcinoma by mediating the metabolism (Zhang et al. 2021). SLC31A1 can enhance the effect of chemotherapeutic drugs on esophageal squamous cell carcinoma and osteosarcoma (Fujita et al. 2021). ATP7A and ATP7B are P-type Cu-transporting ATPases and retain copper homeostasis in the organism by transferring the copper ions across the cellular membranes (Lenartowicz and Krzeptowski 2010). The activity of dihydrolipoamide branched chain transacylase E2 (DBT) may be enhanced in kidney under hypoxia (Ahn et al. 2015). As an enzyme that activates TCA cycle-associated 2-ketoacid dehydrogenases,  (Solmonson et al. 2022). Dihydrolipoamide dehydrogenase (DLD) is the third catalytic enzyme of the pyruvate dehydrogenase complex (PDHC) and is involved in TCA cycle by converting pyruvate to acetyl-CoA. Deficiency of DLD is associated with metabolic diseases, such as liver failure, encephalopathy, and intellectual disability (Ambrus and Adam-Vizi 2018). PDHA1, PDHB, and DLAT are involved in glycolytic regulation, which was also reported to be close correlation with poor prognosis of GC patients (Song et al. 2019). Cleavage System Protein H (GCSH) involved in glycine metabolize is an enzymatically inactive cofactor of the glycine cleavage system, antisense regulation of GSH can determine viability of breast cancer cells (Adamus et al. 2018). In summary, CRGs are closely associated with the development and prognosis of tumors. Consistent with the above reports, the expression level of 12 CRGs in high-cpS patients with better prognosis was markedly upregulated in this study, which demonstrate CRGs will be important targets in study for treatment of tumors.
Construction of cpS was to quantify CRGs in each sample of GC and explore clinical characteristics, prognosis, and TME of cuproptosis in GC. In the study, HCSG and LCSG were classified according to the median of cpS. LCSG with the lower expression of CRGs had poorer OS than HCSG, which proved the effectiveness of copper-induced cell death to tumor cells (Tsvetkov et al. 2022). However, ESTIMATE analysis and TME cell infiltration indicated patients with low cpS were characterized by high immune activation. Above results are contrary to common studies that hyperimmune infiltration is accompanied by a better prognosis (Jochems and Schlom 2011). Surprisingly, through functional analysis, we found not only immune activation but also stromal activation pathways were enriched in LCSG. The previous studies suggested that there are tumors with immune excluded phenotype with abundant immune cells, while these cells remained in the stroma surrounding tumor cell nests rather than penetrate the parenchyma of tumor (Chen and Mellman 2017). The stroma can be limited to the tumor envelope or penetrate the tumor itself, which is similar to the immune cells inside the tumor (Gajewski 2015). In addition, clinicopathological features of ACRG cohort showed that LCSG was characterized by EMT subtype, diffuse histological subtype, and high TMN stage. The stromal activity of LCSG was stronger than the immune activity, including the highly expressed angiogenesis, EMT, and TGFβ pathways, which were considered to suppress T-cells. In summary, we suggest stromal activation suppresses activity of TME cells and plays a dominant role in GC of low cpS.
The therapeutic effect of tumors is an inevitable topic. In recent years, immunotherapies for tumors have undoubtedly emerged a great breakthrough. In this study, we found patients with high cpS exhibited a better response to immunotherapy and prognosis, including anti-PD-L1, ACT, and anti-CTLA4. Consisting with the results, HCSG had a better drug sensitivity. Although there were no studies about the relation between cuproptosis and immunotherapy, the previous studies reported intratumoral copper levels influence PD-L1 expression in cancer cells to enhanced immunotherapy (Voli et al. 2020). We suggest cuproptosis may enhance the effect of anti-PD-L1 through increasing intracellular copper accumulation. It is a pity that the mechanism of cuproptosis affecting the tumor therapy by ACT and anti-CTLA-4 has been unclear. However, as a representative of anti-CTLA-4, ipilimumab was approved by the US Food and Drug Administration for the treatment of metastatic melanoma, and ipilimumab plus nivolumab, an anti-PD-L1 antibody, was considered to improve outcomes of treatment (Larkin et al. 2019). Thus, the results of our study suggest cuproptosis may be involved in the mechanism by which combination therapy with targeted drugs enhances the effect of immunotherapy. Considering an excellent effect of predicting drug sensitivity by cpS, we suggest cuproptosis plays a non-negligible regulation role in immunotherapeutic mechanism. Of note, the previous study reported elesclomol can induce copper chelation to inhibit colorectal cancer (Gao et al. 2021). Tsvetkov et al. proved the sensitivity was restored by the addition of copper when cells grown in the absence of serum were resistant to elesclomol (Tsvetkov et al. 2022). These reports and our study suggest cuproptosis is the potential mechanism of tumor cells death induced by elesclomol. Therefore, we conclude cuproptosis is a novel and important study direction for immunotherapy of GC.
The predictive capability of cuproptosis for prognosis of tumors is very promising. At present, there were some studies constructing risk model based on CRGs, and these models have a good capability of predicting prognosis of melanoma, hepatocellular carcinoma, and kidney renal clear cell carcinoma (Lv et al. 2022). In addition, it was reported that the prognosis-related lncRNA signatures were used to develop risk model with excellent predictive capability in soft tissue sarcoma and head and neck squamous cell carcinoma (Yang et al. 2022). However, there have been no studies to explore cuproptosis phenotype-related gene signatures for predicting prognosis. Thus, we first developed geneScore based on 1125 DEGs to explore cuproptosisrelated gene signatures. GeneScore groups exhibited a complete opposite clinicopathological and TME trends compared with the cpS groups. Similarly, patients with high geneScore showed a poor prognosis, higher TMN stage, stronger stromal activation than immune activation. This suggests these signatures had a negative correlation with CRGs in the regulatory mechanism and cuproptosis really benefits GC patients with better prognosis. Finally, we selected 9 prognosis-related signatures from 731 cuproptosis-related gene signatures, and a riskScore model was constructed to facilitate the clinical application of predicting prognosis. Nomogram showed an excellent predicting effect of riskScore to 1-, 3-, and 5-year survival in GC. In summary, quantifying samples based on DEGs and cuproptosis-related gene signatures can comprehensively assess cuproptosis of individual tumor and can further identify tumor TME subtypes and survival status. We suggest our study will provide new targets to study cuproptosis applying for prognosis of GC.
Despite the results of this study generally meet the expectations, there are still several limitations. First, all analyses were based on data from public databases, and all patients enrolled in this study were obtained retrospectively. Therefore, further experiments in vivo and in vitro are needed to confirm our findings. Second, we did not further analyze the relationship between TME with the riskScore. Third, the tissue sample of the patients is inconsistent, including serum and surgically resected tissue, and there are a large number of clinical variables, such as surgery, neoadjuvant chemotherapy, and chemoradiotherapy, not to be further analyzed, which may influence the final analysis results.
In conclusion, this study about CRGs revealed an extensive regulatory mechanism of TME, response to immunotherapy, and clinical prognosis in GC. Construction of prognosis-related signatures and riskScore model pointed a further mechanism exploration of cuproptosis. These findings highlight the significant clinical implications of CRGs and provide novel ideas for the therapeutic application of cuproptosis in GC.

Author contributions
All authors contributed to the study conception and design. HN: collected all public datasets and was responsible for the main analysis, then wrote this original draft. HW and MZ: reviewed and edited this original draft. YN, XC and ZZ: were responsible for prognosis analysis of cpS groups. XH: was responsible for data integration. QZ, PC and JF: were responsible for funding acquisition, supervision of the study. FW: was the lead author of the study and guided selection of analysis. All authors read and approved the final manuscript. All authors read and approved the final manuscript. Data Availability All public datasets enrolled in this study could download from GEO database (https:// www. ncbi. nlm. nih. gov/ geo/) and the UCSC Xena browser (https:// xenab rowser. net/ datap ages/). All data generated or analyzed during this study are included in the supplementary information files of this article.