Identification of Immune-related Cuproptosis Genes in the Prognosis of Thyroid Cancer

doi:10.21203/rs.3.rs-3127004/v1

Background: Thyroid cancer (THCA) has gradually become the highest incidence of endocrine tumors, but its immunotherapy is still in its infancy. Cuproptosis has been reported to be closely related to various tumors. However, the mechanism of cuproptosis in THCA remains obscure.

Methods: TCGA-THCA dataset was downloaded from TCGA GDC portal. We then conducted survival analysis based on the cuproptosis scores (CPs), and constructed the models of RF and LASSO. GSVA and ssGSEA was subsequently performed. CIBERSORT was used to assess immune infiltration. TIDE algorithm was used to predict the likelihood of response to immunotherapy. qRT-PCR was used to confirm the expression of cuproptosis-related differentially expressed genes (CRDEGs) in THCA cells.

Results: We extracted 14 cuproptosis-related genes (CRGs) from TCGA-THCA, and 11 CRGs (SLC25A3, DLD, LIAS, ATP7B, PDHA1, SLC31A1, FDX1, DBT, GCSH, LIPT1 and DLAT) were CRDEGs. Survival analysis showed that the THCA patients with low-CPs had a relatively good prognosis. No matter which model (RF model or LASSO model) was chosen, the survival rate of high-risk group was lower in both Training dataset and Testing dataset, and 5 CRDEGs was selected as potential prognostic genes for THCA. GSVA showed that there were four differentially expressed hallmark gene sets (“APICAL_JUNCTION”, “APOPTOSIS”, “ANGIOGENESIS” and “COAGULATION” ) between high- and low-risk group in RF model and LASSO model. ssGSEA showed that pathways including glycerophospholipid metabolism, carbon metabolism and so on was differentially expressed among THCA patients with high- and low-CPs. CIBERSORT showed that the infiltrations of CD8⁺ T cells, Tregs, M0 Macrophages and M2 Macrophages, etc. exhibited significant differences. The scores of TIDE was higher in the high-CPs group, indicating THCA patients with high-CPs might be more likely to response immune checkpoint inhibits. qRT-PCR showed the mRNA levels of the above 5 CRDEGs were significantly downregulated in FTC-133 cells and 8505C cells.

Conclusions: This study identified 5 CRDEGs (FDX1, LIAS, LIPT1, DLD and PDHA1) and some immune cells associated with THCA, suggesting the treatment feasibility combined with targeting cuproptosis and immunotherapy in THCA.

thyroid cancer

cuproptosis

immune infiltration

FDX1

LIAS

LIPT1

DLD

PDHA1

With the increasing use of high-resolution ultrasonography, thyroid cancer (THCA) is the most common malignant tumor in the entire endocrine system (1), with a gradual upward trend (2). Targeted therapy was gradually being applied to TCHA, two kinase inhibitors, sorafenib (3) and Lenvatinib (4), have emerged for differentiated thyroid cancer, and vandetanib (5) and cabozantinib (6) have been approved for medullogenous thyroid cancer. However, the clinical evidence was insufficient. Therefore, proving personalized treatment for THCA patients with various gene types was the direction of future development.

Copper is an essential cofactor for all living things, and abnormal accumulation of intracellular copper triggers the aggregation of mitochondrial fatty acylated proteins and inhibits the synthesis of Fe-S cluster, eventually leading to a unique type of programmed cell death termed cuproptosis (7). Cuproptosis induction holds promise as a new diagnostic and therapeutic approach, and previous studies have identified many genes that regulate cuproptosis, some positively regulating cuproptosis and others negatively (7). Cuproptosis has been confirmed being associated with various cancers (8-10). Voli et al. (11) and Liu et al. (12) successively reported that copper level affected cancer immunotherapy. Similarly, two others studies showed the feasibility of copper chelators for the treatment cancers (13,14). However, the role of cuproptosis in immunotherapy for THCA remained still unclear .

Based on data analysis of high-throughput sequencing, this study investigated 499 THCA samples and 57 paracancerous tissue samples, based on 11 cuproptosis-related differentially expressed genes (CRDEGs) (SLC25A3, DLD, LIAS, ATP7B, PDHA1, SLC31A1, FDX1, DBT, GCSH, LIPT1, DLAT) to construct and test prognostic risk models and infer how THCA cells affect immune cell infiltration through cuproptosis. Overall, we have demonstrated an association between cuproptosis and THCA, thus indicating the value of further studies targeting this cellular mechanism.

2.1 Data Collection

The Cancer Genome Atlas (TCGA)-THCA dataset and the somatic mutation data of THCA was downloaded from TCGA Genomic Data Commons (GDC) portal (https://gdc.cancer.gov/) with the data type selected was “fragments per kilobase million” (FPKM) and “Masked Somatic Mutation”, respectively. TCGA-THCA included 499 THCA samples and 57 paracancerous tissue samples and the clinical data of THCA patients (n=499). The clinical baseline data could be seen in Table 1. The package “maftools” in R (15) was applied to visualize the somatic mutation and obtain the tumor mutation burden (TMB) of all patients. To analyze copy number changes of key genes in TCGA-THCA patients, we downloaded the patients’ “Masked Copy Number Segment” data (n=500) by the R package “TCGAbiolinks” (16).

In this study, cuproptosis-related genes (CRGs) were obtained according to a recent study reported by Bian et al. (17) for subsequent phenotype correlation analysis. We obtained 14 CRGs in total: SLC25A3, DLD, DLST, LIAS, ATP7B, PDHA1, SLC31A1, FDX1, DBT, GCSH, LIPT1, DLAT, ATP7A and PDHB.

2.2 Expression Analysis of CRGs

We compared the differential expression of CRGs between THAC group and control group, and boxplots were drawn to show the results. The CRGs with P-values less than 0.05 were defined as CRDEGs. Then we computed the correlation between CRDEGs using the Pearson analysis, and the results were visualized in heat map.

2.3 Building and testing of prognostic model based on CRDEGs

Prognostic models were contrasted using random forest (MF) and least absolute shrinkage and selection operator (LASSO) to evaluate the prognostic impact of CRDEGs on THCA patients. LASSO eigengenes were selected using R package “glmnet” (18). We next used “randomForestSRC” package (19) to accomplish the screening of THCA prognostic genes in RF model. Finally, according to the risk regression coefficients of all genes in each model, we constructed a risk score calculation formula in the LASSO and RF prognostic models, respectively. The calculation formula is:

In this study, THCA samples were randomly grouped into Training group (60% data) and Testing group (40% data). In addition, The “maxstst” package in R (20) was applied to calculate the optimal value of cut-off, which was the basis for grouping of THCA patients. Kaplan-Meier (K-M) method was applied for survival analysis in Training group and Testing group using “survival” in R (21), respectively. The data was visualized by the “survminer” in R package.

2.4 GSVA of Prognostic Models

Gene Set Variation Analysis (GSVA) was conducted to analyze the common pathways and calculate the score of the differentially expressed pathway between groups by using “GSVA” package (22) and “limma” package (23), respectively. The reference gene set “h.all.v7.5.1.symbols” was downloaded from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/) (24). The parameters were set as P. adjust < 0.05 and |log₂FC| > 0.1, which finally displayed through a heat map.

2.5 Difference Analysis Based on Prognostic Model

Differentially expressed genes (DEGs) were identified using the “limma” package by setting the parameters as P. adjust < 0.05 and |log₂FC| > 0.3, and “ggplot2” package (25) was used to map the volcano plot. We then performed differential analysis in LASSO prognostic model and RF prognostic model, respectively. We finally took the intersection of the results of differential analysis to obtain the final DEGs for subsequent analysis.

2.6 GO and KEGG Functional enrichment analysis

The “clusterProfiler” package (26) was used to conduct the Gene Ontology (GO) annotation analysis (27) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (28) on DEGs with false discovery rate (FDR) < 0.05. The top 8 terms for GO and KEGG were shown in bar graphs.

2.7 Calculation of CPs Based on TCGA-THCA Gene Expression Matrix

First, we combined the prognostic genes screened by LASSO and RF to obtain the final prognostic genes of THCA. Based on the gene expression matrix, Single sample Gene Set Enrichment Analysis (ssGSEA) was conducted to calculate the cuproptosis scores (CPs) of all samples in TCGA-THCA dataset using the “GSVA” package. Subsequently, “maxstat” package was performed to calculate the optimal cut-off of CPs, then the TCHA patients were grouped according to the cut-off value of CPs. The R package “survival” was used to draw K-M survival curves.

2.8 GSEA between CPs Models

Gene Set Enrichment Analysis (GSEA) (29) was conducted to investigate the differences in biological processes among THCA patients with high-CPs and low-CPs using R package “clusterProfiler” with FDR < 0.05.

2.9 Immune Infiltration Analysis

Immune infiltration was assessed by CIBERSORT (https://cibersort.stanford.edu/)(30). The “ggplot2” package in R was used to draw histogram to display the distribution of 22 types of immune cells in all samples. We also compared the different infiltration of immune cells among THCA patients according to CPs.

In addition, Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu) (31) was used to predict the response to immune checkpoint inhibitors of THCA patients according to CPs.

2.10 Construction of Clinical Prediction Model Based on CPs

We assessed the predictive power of clinicopathological features of THCA patients for overall survival (OS) by using Cox analysis. In addition, R package “rms” (32) was performed to explore the relationship between clinical factors and prognosis, and draw calibration curves for 1-, 3- and 5-year survival, respectively.

2.11 qRT-PCR Analysis

Quantitative real-time PCR (qRT-PCR) was used to detect the mRNA expression of 5 CRDEGs (FDX1, LIAS, PDHA1, LIPT1 and DLD) in two THCA cells (FTC-133 and 8505C) and Nthy-ori 3-1. The total RNA were extracted by TRIpure Total RNA Extraction Reagent (ELK Biotechnology, China) according the manufacturer’s protocol. EntiLink™ 1st Strand cDNA Synthesis Super Mix was used to synthesize first-strand cDNA. qRT-PCR was conducted on QuantStudio 6 Flex System PCR system using EnTurbo™ SYBR Green PCR SuperMix (ELK Biotechnology, China). The mRNA expression of hub genes was represented at a ratio to GAPDH. The primers in this study were shown as Table 7.

2.12 Statistical Analysis

R software 4.1.2 (https://www.r-project.org/) was used to perform all statistical analysis. Independent t-test and Wilcoxon rank-sum test conducted between two groups. One-way ANOVA was conducted among three groups in qRT-PCR, followed with Tukey-HSD test for post-hoc comparisons. The “survival” package was performed for survival analysis, which was shown as K-M survival curve. The “survival” package was used to conduct Cox analyses. R package “glmnet” and R package “randomForestSRC” were used to construct the LASSO model and RF model, respectively. Two-side P-value < 0.05 was considered as statistically significant.

3.1 Analysis of CRGs Expression and Overall Mutation Level of THCA Patients

First, we extracted 14 CRGs from TCGA-THCA dataset, and 11 CRGs (SLC25A3, DLD, LIAS, ATP7B, PDHA1, SLC31A1, FDX1, DBT, GCSH, LIPT1 and DLAT) were differentially expressed (Figure 1A). Among them, SLC25A3 was up-regulated in THCA group, while DLD, LIAS, ATP7B, PDHA1, SLC31A3, FDX1, DBT, GCSH, LIPT1 and DLAT were down-regulated in THCA group. We then conducted principal component analysis using CRDEGs, the results showed that tumor samples and normal samples were clearly distinguished (Figure 1B), indicating the importance of CRGs in THCA.

Mutation data were downloaded from the TCGA database to assess the influence of genetic mutations on THCA patients. The results showed that missense accounted for majority of mutations, and the highest frequency of variants in THCA patients was single nucleotide polymorphism with C > T accounting for the vast majority. The top 10 genes in THCA patients sorted by mutation rate was BRAF (59%), FRG1B (59%), NRAS (8%), NBPF10 (6%), NBPF1 (6%), UBBP4 (5%), KRTAP4 −11 (5%), CHEK2 (5%), EP400 (5%) and POTEC (4%) (Figure 1C).

Subsequently, we analyzed the mutation of top 20 genes (BRAF, FGR1B, NRAS, NBPF10, NBPF1, UBBP, KRTAP4-11, CHEK2, EP400, POTEO, TG, TTN, HRAS, FAM47C, MACF1, ZNF814, MT-ND5, MUC4, ZFHX3 and RPTN) with the highest mutation frequency in THCA patients, the results showed that missense mutation accounted for the main part of mutations (Figure 2A). At the same time, we analyzed the mutations and copy number variations (CNVs) of the above 11 CRDEGs among THCA patients, and missense mutations also accounted for the vast majority of these mutations. In addition, ATP7B exhibited the highest mutation frequency was ATP7B, and the main type of CNVs was censoring (Figure 2B).

3.2 Construction of THCA Prognostic Model

Pearson’s correlation analysis showed that the expressions of most CRDEGs had positive correlation with other CRDEGs, among which the highest correlation (r = 0.75) was between DLD and DLAT (Figure 3A), indicating the similar expression patterns among CRDEGs in THCA patients.

In order to identify whether prognostic genes were of diagnostic value, a RF model was constructed based on 11 CRDEGs (SLC25A3, DLD, LIAS, ATP7B, PDHA1, SLC31A1, FDX1, DBT, GCSH, LIPT1 and DLAT). We found that the accuracy of model was highest when the included variable was 5 (Figure 3B), and 5 genes (LIAS, LIPT1, DLD, FDX1, PDHA1) in model with the highest accuracy was selected as potential prognostic genes for THCA. Based on the above 5 genes in Cox regression models. The samples were grouped into high-risk group (44 samples) and low-risk group (225 samples) according to risk score and cut-off identified by the “maxstat” package.

As depicted in Figure 3C, LASSO model was also constructed. The characteristic parameters decreased and the absolute value of the coefficient increased with increasing λ. After simulating and selecting the number of features, we got two models (the best model and the simplest model) (Figure 3D). We chose the simplest model to construct the THCA prognostic model and identified 6 genes (SLC25A3, LIAS, ATP7B, PDHA1, SLC31A1 and GCSH) as potential prognostic genes of THCA. We next calculated the risk scores of all samples in accordance with the Cox risk regression coefficient of the above 6 genes. All samples are grouped into high-risk group (31 samples) and low-risk group (268 samples), the basis of grouping was consistent with RF model.

We next performed survival analysis in RF model and LASSO model, respectively, according to the above groupings. The results showed that no matter which model (RF model or LASSO model) was chosen, the survival rate of high-risk group was lower in both Training dataset and Testing dataset (Figure 3E, 3F, 3G and 3H). Meanwhile, the univariate COX regression results showed that the high-risk group exhibited higher degree of tumor progression regardless of model was chosen (Figure 4A and 4B). Finally, we compared the expressions of 5 prognostic genes (LIAS, LIPT1, DLD, FDX1 and PDHA1) of RF prognostic model and 6 prognostic genes (SLC25A3, LIAS, ATP7B, PDHA1, SLC31A1 and GCSH) of LASSO prognostic model between high-risk group and low-risk group in the Training dataset. The results represented as heat maps (Figure 3I and 3J) showed that the above prognostic gene differentially expressed between two groups in both two prognostic models.

3.3 GSVA Based on Risk Grouping

GSVA showed there were four differentially expressed hallmark gene sets (“APICAL_JUNCTION”, “APOPTOSIS”, “ANGIOGENESIS” and “COAGULATION” ) between two groups in RF model and LASSO model (Figure 5A and 5B). Among them, “APOPTOSIS” and “ANGIOGENESIS” have been proven to be closely associated with various cancers (35), indicating that the above two prognostic models might be highly accurate in the diagnosis and treatment of THCA.

3.4 Risk Score Grouping Gene Differential Expression Analysis & Functional Enrichment Analysis

We found 1749 DEGs (440 upregulated and 1311 downregulated) in RF prognostic model (Figure 6A) and 754 DEGs (289 upregulated and 465 downregulated) in LASSO prognostic model (Figure 6B) by using R package “limma”. The DEGs in RF prognositic model was shown in Table 2. Then 608 DEGs shown as Figure 6C and Table 3 were obtained by taking the intersection of DEGs in RF model and LASSO model.

As depicted in Figure 6D, 6E and 6F, GO enrichment analysis based on the above 608 DEGs indicated that DEGs were associated with biological processes (BP) (hormone transport, hormone secretion, signal release, etc.), cellular components (CC) (ion channel complex, cation channel complex, transmembrane transporter complex, etc.), and molecular functions (MF) (passive transmembrane transporter activity, receptor ligand activity, signaling receptor activator activity, etc.). The detailed results were shown in Table 4. KEGG functional enrichment analysis represented in Figure 6G and Table 5 showed that the DEGs had a main impact on pathways including cytokine−cytokine receptor interaction, viral myocarditis, cell adhesion molecules, arachidonic acid metabolism, Th17 cell differentiation, rheumatoid arthritis and so on.

3.5 CPs-based survival analysis and GSEA

9 prognostic genes (SLC25A3, LIAS, ATP7B, PDHA1, SLC31A1, GCSH, LIPT1, DLD and FDX1) of THCA was obtained by taking intersection of prognostic genes between RF prognostic model and LASSO prognostic model. Based on the above 9 prognostic genes, we calculated the CPs of all THCA patients using ssGSEA algorithm. THCA patients then was divided into two groups according to the optimal cut-off value of prognostic analysis based on CPs. As shown in Figure 7A, the results of survival analysis showed that the THCA patients with low-CPs had a relatively good prognosis (log-rank test P = 0.022). GSEA based on the fold differences of genes were shown in Table 6. The top 8 pathways with the smallest P-value were depicted in Figure 7B, among which glycerophospholipid metabolism (P = 0.0216) (Figure 7C), carbon metabolism (P = 0.0195) (Figure 7E), insulin signaling pathway (P = 0.0251) (Figure 7H), mineral absorption (P = 0.0251) (Figure 7I) and metabolic pathways (P = 0.0314) (Figure 7J) were significantly upregulated, and cell adhesion molecules (P = 0.0398) (Figure 7D), Th17 cell differentiation (P = 0.0375) (Figure 7F) and transcriptional mis-regulation in cancer (P = 0.0437) (Figure 7G) were significantly downregulated.

3.6 Transcription levels of 5 hub genes in THCA cells

To further identified our conclusions, the mRNA levels of the 5 hub genes (FDX1, LIAS, PDHA1, LIPT1 and DLD) were detected in Nthy-ori 3-1 cells, FTC-133 cells and 8505C cells by qRT-PCR using GAPDH as internal control. As depicted in Figure 8, the mRNA levels of FDX1, LIAS, PDHA1, LIPT1 and DLD were significantly downregulated in FTC-133 cells and 8505C cells. These results were consistent with the above prognostic models.

3.7 Immune Infiltration Analysis

As shown in Figure 9A, a panorama of 21 immune cells infiltration in all samples was constructed via CIBERSORT by removing immune cells when the sampling fractions of which were 0. We then explored the immune cells with significant differences according to the fractions of infiltration. As depicted in Figure 9B, the infiltrations of CD8⁺ T cells (P = 9.9e-05), Tregs (P = 0.01137), M0 Macrophages (P = 0.01134), M2 Macrophages (P = 0.00026 ), etc. exhibited significant differences in both high-CPs group and low-CPs group. We next separately analyzed and visualized the correlations among 21 immune cells. As showed in Figure 9C and 9D. In the high-CPs group, NK cells resting had significant positive correlation with Neutrophils (r = 0.36, P = 1.7e-06), and the negative correlation between CD8⁺ T cells and Resting memory CD4⁺ T cells was highest (r = -0.47, P = 6e-11). In the low-CPs group, Follicular helper T cells had significant positive correlation with Macrophages M1 (r = 0.46, P < 2e-16), and the negative correlation between Follicular helper T cells and Mast cells resting was highest (r = -0.45, P < 2e-16).

3.8 CPs-based Immunotherapy & Prognostic Model Construction

We further analyzed the characteristics of immunotherapy in two groups. The results showed that the overall TIDE assessment in high-CPs group was lower (Figure 10A). Among them, the scores of TIDE (P = 0.032) (Figure 10B) were higher in the high-CPs group, indicating the THCA patients with high-CPs might be more likely to response to immune checkpoint inhibitors.

Univariate Cox regression analysis and multivariate Cox regression analysis were conducted according to TNM stage, sex, age and grouping, the results showed that the age exerted significant impact on OS (Figure 10C and 10D). We also created a nomogram (Figure 10E) and a calibration curve (Figure 10F, 10G and 10H) using the R package “rms” package to estimate the rates of 1-, 3- and 5-year survival of THCA patients, and the results show that the nomogram model for 1-year survival was most consistent with the ideal model, and other nomogram models were basically consistent with the ideal model, indicating that our model had relatively high accuracy.

The incidence of THCA continues to rise worldwide, but the death rate from THCA has changed little over the past 50 years. Accurate identification of more serious or high-risk THCA has gradually become a crucial link in the diagnosis and treatment of THCA, and is also an urgent issue needed to be solved. Copper is a trace element whose homeostasis is essential for the normal operation of various enzymatic reactions, organs and metabolic processes. However, previous report have found that excess and accumulation of copper may lead to oxidative stress and cytotoxicity (33). Multiple studies have shown that the imbalance of copper homeostasis plays an important impact on tumors (34-37), and copper exerts an essential role in tumor immunity and antineoplastic therapies (38). However, its potential therapeutic value in THCA remains obscure. Cuproptosis is a new mode of cell death relying on copper-dependent targeting lipoylated tricarboxylic acid (TCA) cycle proteins (39). Here, we detailed how the pattern of cuproptosis in THCA is clinically significant. In addition, a cuproptosis scoring system has been constructed to assess the cuproptosis of individual THCA patient for raising understanding of THCA, eventually helping clinicians diagnose THCA as early as possible and choose more effective and personalized treatment options wherever possible.

Thyroid is an organ with prominent vascularization. Angiogenesis plays a significant role in promoting THCA progression. Copper is capable of promoting the progression and angiogenesis of tumor. Currently, targeting angiogenesis is an effective way to treat THCA. Apoptosis is involved in the development and progression of THCA, and recent studies had demonstrated that many natural substances had a pro-apoptotic effect in THCA cells, such as capsaicin (40), berberine (41) and myricetin (42). The above studies informed that the prediction models based on CRGs had high accuracy.

In our study, we used CPs to quantity the degree of cuproptosis, and THCA patients with high-CPs was associated with better outcome, and the prediction models in our study were relatively strong. To our knowledge, we firstly conducted the scheme described above to evaluate the relationship between cuproptosis THCA. Furthermore, we identified 5 CRGs (FDX1, LIAS, LIPT1, PDHA1, DLD) to predict the prognosis for patients with THCA, which was confirmed in vitro. FDX1, a core molecular of cuproptosis reducing Cu²⁺ to Cu⁺, was downregulated in most kinds of cancers (43). Lipoic acid synthetase (LIAS), a key enzyme synthesizing lipoate, has been proven to be involved in several metabolic diseases, including cancer (44).

LIPT1, Pyruvate dehydrogenase E1 component subunit alpha(PDHA1) and DLD are key enzymes involved in TCA cycle, the above three genes has been proven to be correlated with the prognosis of cancer (45-47). Golub et al(10) found that knockout of FDX1, LIAS, LIPT1, PDHA1 and DLD reduced the killing effect of two structurally different copper ionophores (elesclomol and diethyldithiocarbamate).

M1 macrophages exhibit an anti-tumor role, whereas M2 macrophages promotes the growth and metastasis of tumor. Therefore, converting the M2 macrophages into M1 macrophages is a promising treatment approach for cancers. Liu et al. (48) found that low dose treatment with bleomycin reversed M2 macrophage polarization to M1, eventually exhibiting anti-cancer effects on THCA cells. Our study showed that patients with high CPs was related to better survival prognosis. However, the fraction of M2 in patients with high-CPs was higher. Hence, we speculated that combination copper ionophores with targeting TAM may showed better efficacy. Yang et al. (49) determined that CD8⁺ T lymphocytes infiltration improved the survival of THCA. Tregs promotes the progression of tumor via inhibiting antitumor immunity and promoting angiogenesis (50). In our study, the fraction of Treg cells in patients with low-CPs was higher, indicating that cuproptosis might increase immune response via inhibiting Tregs infiltration and increasing Tregs infiltration. In order to confirm the immunotherapy based on CPs, we predicted immunotherapy efficacy between two groups. High-CPs group had a higher TIDE score, indicating that THCA patients with high-CPs might be more suitable for immune checkpoint inhubitors.

In summary, this study demonstrated that cuproptosis was associated with a good prognosis of THCA, and we identified five genes (FDX1, LIAS, LIPT1, DLD and PDHA1)as significant biomarkers in THCA. In addition, our study suggested that the treatment feasibility combined with targeting cuproptosis and immunotherapy in THCA. There are several limitations of this study. First, this study was based on a secondary analysis of public database which had selection bias expected to affect the accuracy of our results. Therefore, more clinical evidences were urgently needed to verify our findings. In addition, our study did not explore the direct action mechanisms of identified CRDEGs and immune cells on THCA. Therefore, In vivo and in vitro studies combined with the latest techniques (For example, single-cell sequencing, spatial transcriptomics research and Olink Proteomics, etc.) were needed to detect how the above CRDEGs and immune cells acting on THCA.

Ethics Approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of date and materials

The data that support the findings of this study are openly available (https://gdc.cancer.gov/) .

Competing interest

The authors declare no competing interests.

Funding

Basic and Applied Basic Research Foundation of Guangdong Province (Grant No: 2020A1515010775), National Natural Science Foundation of China (Grant No: 81704002), Science and Technology Foundation of Shenzhen city (Grant No: JCYJ20210324111212033)

Author’s contributions

Lin Zeng, Huilin Li and Deliang Liu：Conceptualization, Hengxia Zhao, Huilin Li and Deliang Liu: Methodology, Lin Zeng, Xin Qu and Ximing Wang: Software, Lin Zeng and Ni Xiong: Formal analysis, Lin Zeng and Ni Xiong: Writting, Huilin Li and Deliang Liu: Project administration and Funding acquisition.

Siegel RL, Miller KD, Jemal A. Cancer statistics, 2020. CA: a cancer journal for clinicians (2020) 70(1), 7-30. Doi: 10.3322/caac.21590
Global Burden of Disease Cancer Collaboration, Fitzmaurice C, Abate D, Abbasi N, Abbastabar H, Abd-Allah F, et al. Global, Regional, and National Cancer Incidence, Mortality, Years of Life Lost, Years Lived With Disability, and Disability-Adjusted Life-Years for 29 Cancer Groups, 1990 to 2017: A Systematic Analysis for the Global Burden of Disease Study. JAMA oncology (2019) 5(12): 1749-1768. doi: 10.1001/jamaoncol.2019.2996
Brose MS, Nutting CM, Jarzab B, Elisei R, Siena S, Bastholt L, et al. Sorafenib in radioactive iodine-refractory, locally advanced or metastatic differentiated thyroid cancer: a randomised, double-blind, phase 3 trial. Lancet (2014) 384(9940): 319–328. doi: 10.1016/S0140-6736(14)60421-9
Schlumberger M, Tahara M, Wirth LJ, Robinson B, Brose MS, Elisei R, et al. Lenvatinib versus placebo in radioiodine-refractory thyroid cancer. The New England journal of medicine (2015) 372(7):621-630. doi: 10.1056/NEJMoa1406470
Wells SA Jr, Robinson BG, Gagel RF, Dralle H, Fagin JA, Santoro M, et al. Vandetanib in patients with locally advanced or metastatic medullary thyroid cancer: a randomized, double-blind phase III trial. Journal of clinical oncology (2012) 30(2):134-141. doi: 10.1200/JCO.2011.35.5040
Elisei R, Schlumberger MJ, Müller SP, Schöffski P, Brose MS, Shah MH, et al. Cabozantinib in progressive medullary thyroid cancer. Journal of clinical oncology (2013) 31(29):3639-3646. doi: 10.1200/JCO.2012.48.4659
Tsvetkov P, Coy S, Petrova B, Dreishpoon M, Verma A, Abdusamad M, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science (2022) 375(6586):1254-1261. doi: 10.1126/science.abf0529
Ge EJ, Bush AI, Casini A, Cobine PA, Cross JR, DeNicola GM, et al. Connecting copper and cancer: from transition metal signalling to metalloplasia. Nature reviews Cancer (2022) 22(2):102-113. doi: 10.1038/s41568-021-00417-2.
Shanbhag VC, Gudekar N, Jasmer K, Papageorgiou C, Singh K, Petris MJ. Copper metabolism as a unique vulnerability in cancer. Biochimica Biophysica Acta Mol Cell Res (2021) 1868(2):118893. doi: 10.1016/j.bbamcr.2020.118893
da Silva DA, De Luca A, Squitti R, Rongioletti M, Rossi L, Machado CML, et al. Copper in tumors and the use of copper-based compounds in cancer treatment. J Inorg Biochem (2022) 226:111634. doi: 10.1016/j.jinorgbio.2021.111634
Voli F, Valli E, Lerra L, Kimpton K, Saletta F, Giorgi FM, et al. Intratumoral Copper Modulates PD-L1 Expression and Influences Tumor Immune Evasion. Cancer Res (2020) 80(19):4129-4144. doi: 10.1158/0008-5472.CAN-20-0471
Liu YL, Bager CL, Willumsen N, Ramchandani D, Kornhauser N, Ling L, et al. Tetrathiomolybdate (TM)-associated copper depletion influences collagen remodeling and immune response in the pre-metastatic niche of breast cancer. NPJ Breast Cancer (2021) 7(1):108. doi: 10.1038/s41523-021-00313-w
Bian Z, Fan R, Xie L. A Novel Cuproptosis-Related Prognostic Gene Signature and Validation of Differential Expression in Clear Cell Renal Cell Carcinoma. Genes (2022) 13(5):851. doi: 10.3390/genes13050851
Zhou P, Qin J, Zhou C, Wan G, Liu Y, Zhang M, et al. Multifunctional nanoparticles based on a polymeric copper chelator for combination treatment of metastatic breast cancer. Biomaterials (2019) 195:86-99. doi: 10.1016/ j.biomaterials.2019.01.007
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol (2004) 5(10):R80. doi: 10.1186/gb-2004-5-10-r80
Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, et al. TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res (2016) 44(8):e71. doi: 10.1093/nar/gkv1507
Bian Z, Fan R, Xie L. A Novel Cuproptosis-Related Prognostic Gene Signature and Validation of Differential Expression in Clear Cell Renal Cell Carcinoma. Genes (2022) 13(5):851. doi: 10.3390/genes13050851
Friedman J, Hastie T, Tibshirani R, Narasimhan B, Tay K, Simon N, et al. glmnet: Lasso and Elastic-Net Regularized Generalized Liner Models Version 4.1-4. (2022) R package Version 3.6.0. URL https://cran.r-project.org/package=glmnet
Ishwaran H, Kogalur UB. randomForestSRC: Fast Unified Random Forests for Survival, Regression, and Classification (RF-SRC) Version 3.1.1. (2022) R package version 3.6.0, URL https://cran.r-project.org/package=randomForestSRC
Hothorn T. maxstat: Maximally Selected Rank Statistic Version 0.7-25. (2017) R package version 1.7.0, URL https://cran.r-project.org/package=maxstat
Therneau TM, Lumley T, Elizabeth Atkinson, Cynthia C. survival: Survival Analysis Version 3.1.1. (2022) R package version 3.6.0, URL https://cran.r-project.org/package=survival
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC bioinformatics (2013) 14:7. doi: 10.1186/1471-2105-14-7
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi:10.1093/nar/gkv007
Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics (2011) 27(12):1739–1740. doi: 10.1093/bioinformatics/btr260
Wickham H, Chang W, Henry L, Pedersen TL, Takahashi K, Wilke C, et al. ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics Version 3.4.0. (2022) R package version 3.3, URL https://cran.r-project.org/package=ggplot2
Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) (2021) 2(3):100141. doi:10.1016/j.xinn.2021.100141
Gene Ontology Consortium. Gene Ontology Consortium: going forward. Nucleic Acids Res (2015) 43(Database issue):D1049-D1056. doi:10.1093/nar/gku1179
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res (2000) 28(1):27-30. doi:10.1093/nar/28.1.27
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA (2005) 102(43): 15545-15550. doi:10.1073/pnas.0506580102
Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling Tumor Infiltrating Immune Cells with CIBERSORT. Methods Mol Biol (2018) 1711:243-259. doi:10.1007/978-1-4939-7493-1_12.
Wang Q, Li M, Yang M, Yang Y, Song F, Zhang W, et al. Analysis of immune-related signatures of lung adenocarcinoma identified two distinct subtypes: implications for immune checkpoint blockade therapy. Aging (2020) 12(4):3312-3339. doi:10.18632/aging.102814.
Harrell FE Jr. rms: Rrgressions Modeling Strategies Version 6.3-0. (2022) R package version 3.5.0, URL https://cran.r-project.org/package=rms
Saporito-Magriñá CM, Musacco-Sebio RN, Andrieux G, Kook L, Orrego MT, Tuttolomondo MV, et al. Copper-induced cell death and the protective role of glutathione: the implication of impaired protein folding rather than oxidative stress. Metallomics (2018) 10(12):1743-1754. doi:10.1039/c8mt00182k.
Jiang Y, Huo Z, Qi X, Zuo T, Wu Z. Copper-induced tumor cell death mechanisms and antitumor theragnostic applications of copper complexes. Nanomedicine(Lond) (2022) 17(5):303-324. doi:10.2217/nnm-2021-0374.
Kaur P, Johnson A, Northcote-Smith J, Lu C, Suntharalingam K. Immunogenic Cell Death of Breast Cancer Stem Cells Induced by an Endoplasmic Reticulum-Targeting Copper(II) Complex. Chembiochem (2020) 21(24):3618-3624. doi:10.1002/cbic.202000553.
Nagai M, Vo NH, Shin Ogawa L, Chimmanamada D, Inoue T, Chu J, et al. The oncology drug elesclomol selectively transports copper to the mitochondria to induce oxidative stress in cancer cells. Free Radic Biol Med (2021) 52(10):2142-2150. doi:10.1016/j.freeradbiomed.2012.03.017.
Gao W, Huang Z, Duan J, Nice EC, Lin J, Huang C. Elesclomol induces copper-dependent ferroptosis in colorectal cancer cells via degradation of ATP7A. Mol Oncol (2021) 15(12):3527-3544. doi:10.1002/1878-0261.13079.
Percival SS. Copper and immunity. Am J Clin Nutr (1998) 67(5 Suppl):1064S–1068S. doi: 10.1093/ajcn/67.5.1064S.
Tsvetkov P, Coy S, Petrova B, Dreishpoon M, Verma A, Abdusamad M, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science (2022) 375(6586):1254-1261. doi: 10.1126/science.abf0529
Xu S, Cheng X, Wu L, Zheng J, Wang X, Wu J, et al. Capsaicin induces mitochondrial dysfunction and apoptosis in anaplastic thyroid carcinoma cells via TRPV1-mediated mitochondrial calcium overload. Cell Signal (2020) 75:109733. doi:10.1016/j.cellsig.2020.109733
Li L, Wang X, Sharvan R, Gao J, Qu S. Berberine could inhibit thyroid carcinoma cells by inducing mitochondrial apoptosis, G0/G1 cell cycle arrest and suppressing migration via PI3K-AKT and MAPK signaling pathways. Biomed Pharmacother (2017) 95:1225-1231. doi:10.1016/j.biopha.2017.09.010
Ha TK, Jung I, Kim ME, Bae SK, Lee JS. Anti-cancer activity of myricetin against human papillary thyroid cancer cells involves mitochondrial dysfunction-mediated apoptosis. Biomed Pharmacother (2017) 91:378–384. doi: 10.1016/j.biopha.2017.04.100.
Xiao C, Yang L, Jin L, Lin W, Zhang F, Huang S, et al. Prognostic and immunological role of cuproptosis-related protein FDX1 in pan-cancer. Front Genet (2022) 13:962028. doi:10.3389/fgene.2022.962028.
Mayr JA, Feichtinger RG, Tort F, Ribes A, Sperl W. Lipoic acid biosynthesis defects. J Inherited Metab Dis (2014) 37(4):553-563. doi: 10.1007/s10545-014-9705-8.
Chen Y, Xu T, Xie F, Wang L, Liang Z, Li D, et al. Evaluating the Biological Functions of the Prognostic Genes Identified by the Pathology Atlas in Bladder Cancer. Oncol Rep (2021) 45(1):191-201. doi:10.3892/or.2020.7853.
Deng L, Jiang A, Zeng H, Peng X, Song L. Comprehensive analyses of PDHA1 that serves as a predictive biomarker for immunotherapy response in cancer. Front Pharmacol (2022) 13:947372. doi: 10.3389/fphar.2022.947372.
Wang S, Xing N, Meng X, Xiang L, Zhang Y. Comprehensive bioinformatics analysis to identify a novel cuproptosis-related prognostic signature and its ceRNA regulatory axis and candidate traditional Chinese medicine active ingredients in lung adenocarcinoma. Front Pharmacol (2022) 13:971867. doi.org/10.3389/fphar.2022.971867
Liu H, Dong H, Jiang L, Li Z, Ma X. Bleomycin inhibits proliferation and induces apoptosis in TPC-1 cells through reversing M2-macrophages polarization. Oncol Lett (2018) 16(3):3858-3866. doi: 10.3892/ol.2018.9103.
Ha TK, Jung I, Kim ME, Bae SK, Lee JS. Anti-cancer activity of myricetin against human papillary thyroid cancer cells involves mitochondrial dysfunction-mediated apoptosis. Biomed Pharmacother (2017) 91:378–384. doi:10.1016/j.biopha.2017.04.100.
Facciabene A, Motz GT, Coukos G. T-regulatory cells: key players in tumor immune escape and angiogenesis. Cancer Res (2012) 72(9):2162-2171. doi:10.1158/0008-5472.CAN-11-3687.

Table 1 THCA sample baseline data table in TCGA database

Characteristic	Number of cases
T_stage
T1	142(28%)
T2	164(33%)
T3	169(34%)
T4	24(5%)
N_stage N0 N1 N2 M_stage M0 M1 Overall survival	227(45%) 223(45%) 49(10%) 282(57%) 217(43%)
Dead	483(97%)
Alive	16(3%)
Gender
Female	365(73%)
Male	134(27%)

Table 7 The primer sequence of qRT-PCR

Primer mane	Primer information	Sequence（5`-3`）		Tm value	CG%	Product
H-GAPDH	NM_001256799.2	sense	CATCATCCCTGCCTCTACTGG	59.4	57.1	259
H-GAPDH	NM_001256799.2	antisense	GTGGGTGTCGCTGTTGAAGTC	60.1	57.1	259
H-LIAS	NM_006859	sense	CAAGATTTTGTATCTGGTGATCTT	55.2	33.3	113
H-LIAS	NM_006859	antisense	TCTGTCTTTAGCCATGGAGGTAGT	60.8	45.8	113
H-PDHA1	NM_000284	sense	TGAAAGCAGATCAGCTGTATAAAC	57.3	37.5	257
H-PDHA1	NM_000284	antisense	GTAGAAGTTCTTGGCATACATGTG	57.9	41.6	257
H-LIPT1	NM_015929	sense	TATTATGTAGTACTGATGGGACGTT	57.1	36	116
H-LIPT1	NM_015929	antisense	CTTTTCCAAAAGATTTTTCACTAAG	53.8	28	116
H-DLD	NM_000108	sense	CAATAGTGTCATCTACAGGTGCTTT	59.1	40	197
H-DLD	NM_000108	antisense	GATGCGTTGAAAGTTTTTAGATATC	55.4	32	197
H-FDX1	NM_004109	sense	GTTGGTGATTCTCTGCTAGATGTT	58.8	41.6	147
H-FDX1	NM_004109	antisense	CTCATCAGTGATTGCATCTAACTTC	58	40	147

Table 2 To 6 are available in Supplementary Files section.

No competing interests reported.

Identification of Immune-related Cuproptosis Genes in the Prognosis of Thyroid Cancer

Status:

Version 1

Abstract

Figures

Introduction

Materials And Methods

Results

Discussion

Declarations

References

Tables

Additional Declarations

Supplementary Files

Status:

Version 1