Construction of a novel cuproptosis-related gene signature for predicting prognosis and estimating tumor immune microenvironment status in papillary thyroid carcinoma
BMC Cancer volume 22, Article number: 1131 (2022)
Cuproptosis, a new form of programmed cell death, has been recently reported to be closely related to tumor progression. However, the significance of cuproptosis-related genes (CRGs) in papillary thyroid carcinoma (PTC) is still unclear. Therefore, this study aimed to investigate the role of the CRG signature in prognosis prediction and immunotherapeutic effect estimation in patients with PTC.
RNA-seq data and the corresponding clinical information of patients with PTC were obtained from the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. Comprehensive analyses, namely, consensus clustering, immune analyses, functional enrichment, least absolute shrinkage and selection operator-multivariate Cox regression, and nomogram analysis, were performed to identify new molecular subgroups, determine the tumor immune microenvironment (TIME) status of the identified subgroups, and construct a clinical model. Independent verification cohort data and quantitative real-time polymerase chain reaction (qPCR) was performed to validate the expression of specific prognosis-related and differentially expressed CRGs (P-DECRGs).
In the TCGA database, 476 patients with PTC who had complete clinical and follow-up information were included. Among 135 CRGs, 21 were identified as P-DECRGs. Two molecular subgroups with significantly different disease-free survival and TIME statuses were identified based on these 21 P-DECRGs. The differentially expressed genes between the two subgroups were mainly associated with immune regulation. The risk model and nomogram were constructed based on four specific P-DECRGs and validated as accurate prognostic predictions and TIME status estimation for PTC by TCGA and GEO verification cohorts. Finally, the qPCR results of 20 PTC and paracancerous thyroid tissues validated those in the TCGA database.
Four specific P-DECRGs in PTC were identified, and a clinical model based on them was established, which may be helpful for individualized immunotherapeutic strategies and prognostic prediction in patients with PTC.
The incidence of thyroid carcinoma, the most common endocrine malignancy, has continuously and rapidly increased worldwide in recent decades . However, thyroid carcinoma mortality remains relatively stable at low levels . In 2021, the estimated numbers of new thyroid cancer cases and deaths in the United States were 44,280 and 2,200, respectively . Similar thyroid carcinoma trends have been reported in China . The most frequent histological type of thyroid carcinoma is papillary thyroid carcinoma (PTC), which accounts for approximately 90% of all cases . Most PTCs are well-differentiated, indolent, and have excellent prognoses when conventional treatments are implemented. However, several clinicopathological features, including age, tumor diameter, extrathyroidal extension, and cervical lymph node and distant organ metastases, are considered significantly unfavorable prognostic factors and result in a relatively high recurrence rate . Our institution has reported the risk factors associated with cervical lymph node metastasis of PTC and provided some clinical suggestions regarding prophylactic central lymph node dissection . To date, the molecular mechanisms of recurrence and metastasis in PTC are not fully understood. Therefore, a deeper investigation of potential therapeutic targets and more efficient prognostic models for PTC are of great clinical significance, especially for facilitating personalized treatment strategies.
Copper is an indispensable trace element that plays an essential role in various biological processes. Dysregulation of copper homeostasis, i.e., the abnormal alteration of intracellular copper levels, may trigger cytotoxicity and is considered a novel hallmark of malignant tumor progression . Recently, cuproptosis has been identified as a unique type of cell death related to copper homeostasis disorder, proving that the dysregulation of copper homeostasis plays a pivotal role in tumor development and progression. Additionally, several genes involved in cuproptosis have been identified . However, the role of cuproptosis-related genes (CRGs) in PTC remains poorly understood.
In the present study, we aimed to comprehensively investigate the correlations of CRGs with different survival statuses of patients with PTC, as well as the underlying molecular mechanisms. Our study highlighted the regulatory functions of CRGs in PTC progression. Moreover, the results shed light on novel strategies to predict prognosis and lay a foundation for individual-based therapeutic application in patients with PTC.
Materials and methods
Gene expression profiles and clinical information on thyroid carcinoma from The Cancer Genome Atlas (TCGA) database were obtained from the cBioPortal for Cancer Genomics, including 502 tumor and 58 non-tumor cases [10, 11]. We enrolled patients with PTC who had complete information on age, sex, tumor/lymph node/metastasis (TNM) staging system, follow-up, and survival status (disease-free survival [DFS]). The clinicopathological features of PTC samples from TCGA are presented in Table 1. Samples acquired from seven Gene Expression Omnibus (GEO) databases (GSE3678, GSE6004, GSE29265, GSE33630, GSE35570, GSE53157, and GSE60542) were integrated and defined as an independent verification cohort, including 196 PTC and 160 normal thyroid cases. The matched mapping information of GeneSymbol and ENSG_ID was extracted from the gff3 files and downloaded from GENCODE. The median was obtained when multiple matches existed, and finally, the converted expression spectrum (convert_exp.txt) was obtained. To deal with missing data completion, both genes and cases with unavailable (NA) ratios greater than 50% were removed. To standardize the data, we performed a log2 (X + 1) conversion. Datasets of cuproptosis-related genes (CRGs) were obtained from the Molecular Signatures Database (MSigDB) and a previously published cuproptosis-related study [9, 12].
A total of 20 patients who underwent thyroid surgery and were diagnosed with PTC by postoperative pathology between January 2022 and May 2022 at Shengjing Hospital of China Medical University were included. The inclusion criteria for selecting patients were age > 18 years, a postoperative pathological diagnosis of PTC, and complete patient clinical information. In addition, the exclusion criteria were other treatments, such as radiotherapy, chemotherapy, or immunotherapy, received before surgery; patients with serious systemic diseases, such as severe heart, liver, or kidney diseases and other malignant tumors; women in the gestational and lactational period; patients taking drugs with potential interference; patients who refused or were unable to sign informed consent; other histological types of thyroid carcinoma combined with PTC in postoperative pathological diagnosis. PTC and paracancerous thyroid tissues at least 2 cm away from PTC areas and confirmed as being without PTC cells by two pathologists were collected following protocols approved by the Ethics Committee of Shengjing Hospital of China Medical University. Written informed consent was obtained from all patients enrolled in this study.
Identification of molecular subgroups related to CRGs
Differentially expressed CRGs (DECRGs) in PTC and non-tumor cases from the TCGA database were determined using the Wilcoxon rank-sum test. The association between these DECRGs and PTC prognosis was compared using univariate Cox regression analysis. Next, cluster analysis of these prognosis-related DECRGs (P-DECRGs) was performed by “ConsensusClusterPlus” using agglomerative pam clustering with 1-Pearson correlation distances and resampling 80% of the samples for 10 repetitions . The optimal number of clusters was determined using an empirical cumulative distribution function plot. The R software package “stats” (version 3.6.0) was used for principal component analysis (PCA). Briefly, we first determined the Z-score on the expression spectrum and then used the “prcomp” function for dimension reduction analysis to obtain the reduced dimension matrix.
We used the “survfit” function of the R software package “survival” to analyze the prognostic differences between different groups of samples and the log-rank test method to evaluate the significance.
Differentially expressed genes (DEGs) between molecular subgroups related to P-DECRGs were identified using the R package “linear models for microarray data (Limma, version 3.40.6)” . Specifically, for the expression profile data set we obtained, we used the “lmFit” function to perform multiple linear regression. Further, we used the “eBays” function to compute moderated t- and F-statistics and log-odds of differential expression by empirical Bayes moderation of the standard errors towards a common value. Finally, we obtained the significance of the differences of each gene. For Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, we downloaded the “c5.go.v7.4.symbols.gmt” subset from the molecular signatures database and used the KEGG rest API to obtain the latest gene annotations for KEGG pathway analysis [15,16,17]. Considering these as the background, we mapped the genes into the background set and used the R software package “clusterProfiler” (version 3.14.3) for enrichment analysis of the gene set. The minimum and maximum gene sets were set to 5 and 5,000, respectively. P < 0.05, and false discovery rate (FDR) < 0.1 were considered statistically significant.
The stromal and immune scores and tumor purity, defined as the quantitative proportions of immune and stromal components and malignant cells in tumor samples in different groups, were calculated using the estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE) algorithm . The scores of four different immunophenotypes, including major histocompatibility complex (MHC) molecules, effector cells, immunosuppressive cells, and immune checkpoints, were calculated using the immunophenoscore (IPS) algorithm . The human immune cell subsets were calculated using the cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT) and the quantification of the tumor immune contexture from human RNA-seq data (quanTIseq) algorithms [20, 21]. The responses of patients with PTC to immune checkpoint blockade (ICB) therapy were predicted using the immune cell abundance identifier (ImmuCellAI) algorithm .
Establishment and verification of the risk model related to specific P-DECRGs
The PTC sample dataset from the TCGA database was randomly divided into training and verification cohorts (1:1). We used the R software package “glmnet” to integrate survival time and status and gene expression data and used the least absolute shrinkage and selection operator (LASSO)-Cox method for regression analysis. In addition, a tenfold cross-validation was performed to obtain the minimum lambda, which was defined as the optimal value. In addition, a risk model related to specific P-DECRGs was established using a multivariate Cox regression analysis. We used the R software package “survival” to integrate survival time and status and gene expression data and the Cox method to evaluate the prognostic significance of each gene. Patients with PTC were divided into high- and low-risk groups according to the medium value of the risk score. Next, we used the R software package “pROC” (version 22.214.171.124) to perform a receiver operating characteristic (ROC) analysis to assess the predictive efficiency of the risk model. Specifically, we obtained patient follow-up times and risk scores and used the pROC function to perform ROC analysis at multiple time points and the CI function of pROC to evaluate the area under the curve (AUC) and confidence interval to obtain the final AUC results.
Establishment and verification of a nomogram scoring system related to risk score
By integrating the DFS time, survival status, clinicopathological features, and risk score data, a predictive nomogram was established using the package “rms,” and the prognostic significance of these characteristics was evaluated . In the nomogram scoring system, each variable was matched to a score. The total score of each sample was the sum of all variable scores. Time-dependent ROC curves and calibration plots were used to assess the nomogram-scoring system.
RNA extraction and quantitative real-time polymerase chain reaction (qPCR)
Total RNA was extracted using an RNAiso Plus Kit (Takara, Dalian, China) according to the manufacturer’s instructions. After reverse transcription, cDNA was synthesized, and qPCR analysis was performed using SYBR Premix Real-time PCR Reagent (Takara) in a Roche LightCycler 480 II system (Roche Diagnostics Corporation, Indianapolis, IN, USA) according to the manufacturer’s protocol. The internal control was glyceraldehyde-3-phosphate dehydrogenase (GAPDH). The primer sequences of the selected target genes and internal controls are listed in Table 2. The target gene expression levels were calculated using the 2 –ΔΔCt method, similar to that in one of our previously reported studies .
Enumeration data are presented as numbers/percentages, whereas measurement data are presented as mean ± standard deviation (SD). Student’s t-test was used for the statistical analysis of two groups, and one-way ANOVA was used for three or more groups. Statistical analyses were performed using R software (version 4.0.3; R Foundation for Statistical Computing, Vienna, Austria). Statistical significance was set at P < 0.05.
Identification of P-DECRGs in PTC
A flow chart of the data analysis in the present study is shown in Fig. 1. Among the 502 patients with thyroid carcinoma from the TCGA dataset, 476 patients with PTC who had complete clinical and survival information were selected. To investigate the potential roles of CRGs in PTC, 135 CRGs were obtained from 12 gene sets in MSigDB and a previously published cuproptosis-related study (Additional file 1: Table S1). We matched these 135 genes with the RNA-seq data of 476 PTC and 58 non-tumor samples from the TCGA database. Compared with non-tumor samples, 107 CRGs were differentially expressed in PTC. Among these 107 genes, only 21 P-DECRGs, generated from the univariable Cox analysis results, were selected for subsequent analysis (Additional file 2: Table S2).
Identification of two molecular subtypes based on P-DECRGs in PTC
The consensus clustering approach was used to divide the 476 patients with PTC into subgroups based on the 21 P-DECRGs. The clustering results showed that optimal clustering stability was identified when K = 2 (Fig. 2a and b). A total of 267 patients with PTC were clustered into cluster 1 (C1), and the remaining 209 patients with PTC were clustered into C2 (Fig. 2c). The PCA results revealed significant differences in the expression levels of P-DECRGs between the two molecular subtypes (Fig. 2d). Moreover, patients in C2 had better DFS than those in C1 (P < 0.05; Fig. 2e). These results demonstrate that the 21 P-DECRGs could divide patients with PTC into two molecular subtypes with different prognostic results.
Functional enrichment of differentially expressed genes between the two P-DECRG-based molecular subtypes
To explore the underlying signaling mechanisms between the two P-DECRG-based molecular subtypes, the DEGs between these two clusters were identified. Compared with C2, the results showed that a total of 1,351 DEGs were identified, of which 775 genes were upregulated, and 576 genes were downregulated in C1 (Fig. 2f). GO biological process enrichment analysis revealed that the DEGs were enriched in biological adhesion, cell migration, locomotion, leukocyte migration, and the positive regulation of immune system processes (Fig. 2g). GO cellular component enrichment analysis revealed that the DEGs were enriched in the cell surface, external encapsulating structure, collagen-containing extracellular matrix, side of membrane, and MHC protein complex (Fig. 2h). GO molecular function enrichment analysis revealed that the DEGs were enriched in signaling receptor, antigen, glycosaminoglycan, growth factor, and cell adhesion molecule binding (Fig. 2i). KEGG enrichment analysis revealed that the DEGs were enriched in many pathways, including cell adhesion molecules, viral myocarditis, phagosome, focal adhesion, rheumatoid arthritis, and extracellular matrix (ECM)-receptor interaction (Fig. 2j). All these results demonstrated that the divided molecular subtypes were correlated with the regulation of the immune system process, which may be involved in the DFS of patients with PTC.
Different tumor immune microenvironment (TIME) statuses in the two P-DECRG-based molecular subtypes
To explore the differences in TIME status between the two molecular subtypes, we performed multiple immune analyses. Compared with those in C1, the ESTIMATE algorithm results revealed that patients with PTC in C2 had significantly lower stromal scores (P < 0.05), immune scores (P < 0.0001), and tumor purity (P < 0.0001; Fig. 3a). The IPS algorithm results showed that patients with PTC in C2 had significantly lower MHC molecules (P < 0.0001) and effector cells (P < 0.0001), and higher immunosuppressive cells (P < 0.0001) and immune checkpoints (P < 0.0001; Fig. 3b) compared to those in C1. The CIBERSORT algorithm results indicated that the infiltration levels of naive B cells, CD4 memory activated T cells, regulatory T cells (Tregs), M0 macrophages, resting dendritic cells, activated dendritic cells, and activated mast cells were significantly higher in C1 than those in C2. In contrast, memory B cells, CD8 T cells, resting NK cells, activated NK cells, monocytes, M2 macrophages, and eosinophils had significantly lower infiltration in C1 than those in C2 (Fig. 3c). The quanTIseq algorithm results showed that the infiltration levels of M1 macrophages, neutrophils, NK cells, and Tregs were significantly higher in C1 than those in C2, while CD4 T cells and dendritic cells had significantly lower infiltration in C1 compared to those in C2 (Fig. 3d). The results of the ImmuCellAI algorithm revealed that ICB therapy scores were significantly lower in C1 than those in C2 (Fig. 3e). In addition, the investigation of the expression differences of 33 known types of immune checkpoint molecules showed that 26 types were differentially expressed between the two molecular subtypes (Fig. 3f). All these results demonstrated significant differences in TIME status between the C1 and C2 subtypes.
Establishment of a risk model based on specific P-DECRGs in the training cohort
To assess the prognostic value of P-DECRGs in PTC, a risk signature model based on P-DECRGs in the training cohort was constructed. First, the potential P-DECRGs for establishing the risk model were screened using LASSO regression analysis. The results showed that the optimal lambda value was 0.0164684245862992, and 10 genes were filtered (Fig. 4a and b). In addition, multivariate Cox analysis was performed based on the genes generated from LASSO analysis, and four specific P-DECRGs of PTC were identified: Ferredoxin 1 (FDX1), Cyclin Dependent Kinase 1 (CDK1), Purinergic Receptor P2X 4 (P2RX4), and Microtubule Associated Protein 1 Light Chain 3 Alpha (MAP1LC3A). A risk model of these four specific genes was constructed to calculate the risk score as follows:
Risk score = -2.011 × the expression level of FDX1 + 1.139 × the expression level of CDK1 + 1.662 × the expression level of P2RX4 − 0.948 × the expression level of MAP1LC3A.(1).
Patients with PTC in the training cohort were successfully classified into high- and low-risk groups based on the medium risk score of the established risk model. We analyzed the relationship between different risk scores and patient follow-up times, events, and changes in the expression of various genes. It was observed that, with an increase in risk scores, the survival rate of patients decreased, and the expression of FDX1 and MAP1LC3A showed a downward trend. However, the expression of CDK1 and P2RX4 was upregulated with an increase in risk score (Fig. 4c). Patients with PTC in the high-risk group had a worse DFS than those in the low-risk group (Fig. 4d). The AUCs of the ROC curve for 1-, 3-, and 5-year survival were 0.84, 0.86, and 0.82, respectively, according to the results of time-dependent ROC analysis in the training cohort (Fig. 4e). These results suggest that the constructed risk model based on the four specific P-DECRGs showed considerable potential for predicting the DFS of patients with PTC.
Validation of the risk model based on specific P-DECRGs in the PTC verification cohort of PTC
To validate the constructed risk model based on the four specific P-DECRGs, patients with PTC in the TCGA verification cohort and those in the independent GEO verification cohort were also divided into high- or low-risk groups. Survival analysis revealed that patients from the TCGA verification cohort in the low-risk group had a better prognosis than those in the high-risk group (Fig. 4f). The AUC of the ROC curve for 1-, 3-, and 5-year survival was 0.68, 0.63, and 0.61, respectively, which indicated that the constructed risk model exhibited predictive capacity, according to the results of time-dependent ROC analysis in the verification cohort (Fig. 4g). These results validated that the risk model based on the four specific P-DECRGs was well established and associated with predicting the DFS of patients with PTC in the verification cohort.
Different TIME statuses in the two risk groups based on the constructed risk model
To explore the association between risk score and TIME status, we performed multiple immune analyses. In the training cohort, the ESTIMATE algorithm results revealed that patients with PTC in the high-risk group had significantly higher stromal scores (P < 0.01), immune scores (P < 0.0001), and tumor purity (P < 0.0001) than those in the low-risk group (Fig. 5a). The IPS algorithm results showed that patients with PTC in the high-risk group had significantly higher MHC molecules (P < 0.0001) and effector cells (P < 0.0001), and lower immunosuppressive cells (P < 0.0001) and immune checkpoints (P < 0.0001) than those in the low-risk group (Fig. 5b). The CIBERSORT algorithm indicated that the infiltration levels of naive B cells, plasma cells, CD4 memory activated T cells, follicular helper T cells, Tregs, M0 macrophages, M1 macrophages, eosinophils, and neutrophils were significantly higher in the high-risk group than those in the low-risk group, while memory B cells, resting NK cells, activated NK cells, monocytes, M2 macrophages, and activated mast cells had significantly lower infiltration in the high-risk group than those in the low-risk group (Fig. 5c). The quanTIseq algorithm results showed that the infiltration levels of M1 macrophages and Tregs were significantly higher in the high-risk group than those in the low-risk group, while CD4 T cells and dendritic cells had significantly lower infiltration in the high-risk group compared to those in the low-risk group (Fig. 5d). The ImmuCellAI algorithm results revealed that ICB therapy scores were significantly lower in the high-risk group than those in the low-risk group (Fig. 5e). Among the 33 known types of immune checkpoint molecules, 26 were differentially expressed between the two risk groups (Fig. 5f). Moreover, all these TIME status results in the TCGA verification cohort and independent GEO verification cohort were similar, demonstrating that there were associations between the risk score and TIME status in patients with PTC (Figs. 5g-6f).
Construction and calibration of a nomogram integrating the constructed risk model
A nomogram integrating risk score and clinicopathological features was constructed to predict the DFS of patients with PTC more precisely. The constructed nomogram showed the contribution of the risk score and clinicopathological features to the probability of DFS (Fig. 6g). The C-index of the nomogram in the training cohort reached 0.85 (95% CI: 0.78–0.92; P < 0.05). In addition, the AUCs of the ROC curve for 1-, 3-, and 5-year survival in the training cohort were 0.85, 0.88, and 0.85, respectively (Fig. 6h-i), and similar results were observed in the TCGA verification cohort (Fig. 6j-k). These results demonstrate that the integrated nomogram could accurately predict the DFS of patients with PTC.
Expressional validation of four specific P-DECRGs
The expression levels of four specific P-DECRGs were compared between PTC and non-tumor thyroid tissues in the independent GEO verification cohort. In addition, the expression differences of four specific P-DECRGs between selected PTC and paracancerous thyroid tissues were also estimated by qPCR. The comparison results showed that the mRNA expression levels of CDK1 and P2RX4 in PTC tissues were significantly higher than those in non-tumor thyroid tissues. Moreover, the mRNA expression levels of MAP1LC3A and FDX1 in PTC tissues were significantly lower than those in non-tumor thyroid tissues (Fig. 6l-m). The comparison results were consistent with those of the TCGA datasets.
Most PTCs are indolent, and patients have relatively favorable prognoses with a 10-year survival rate > 90% . Because of the excellent overall survival (OS) of PTC, we paid more attention to another important endpoint (DFS) to evaluate its prognostic significance for patients with PTC in this study. Among the 476 patients with PTC from the TCGA database, the OS rate was 99.58% (474/476), while the DFS rate was 90.34% (430/476). For patients with PTC, DFS more accurately reflects disease status and the impact of the disease on the physical and mental states of patients than OS.
Tumor heterogeneity is a crucial phenomenon involved in the interactions among various factors, including certain intracellular genetic changes and tumor microenvironment influences . As a result, tumor heterogeneity leads to the complexity of tumor cells and diversity in the therapeutic response . Therefore, it is essential to administer specific treatment strategies based on the molecular subtype classification of PTC, in line with the concept of precision medicine . Consensus clustering is a reliable approach in which several different clusters can be obtained to aggregate the clustering results and obtain a better clustering solution .
Accordingly, in this study, patients with PTC were divided into two different clusters, corresponding to two different molecular subtypes, using the consensus clustering method. Recent studies have suggested that copper is closely related to the occurrence and development of tumors . Molecular subtype classification based on CRGs may reveal the exact molecular mechanisms and specific signaling pathways involved in tumor progression, TIME characteristics, and patient prognosis. A prognostic model for cuproptosis-related subtypes has been developed for glioma, providing new insights into tumor prognosis assessment, neoplasm-immune interactions, and potential drug targets . Identifying cuproptosis-related subtypes and developing prognosis models in breast cancer are also useful predictors of prognosis and the tumor microenvironment . In hepatocellular carcinoma, the CRG scoring model may inspire new approaches for both clinically predicting prognosis and developing treatment strategies . However, the relationship between cuproptosis, PTC, and TIME status remains elusive. In the present study, patients with PTC were divided into clusters corresponding to two different molecular subtypes based on P-DECRGs. Further research indicated the diversity of prognosis and revealed the distinction of molecular mechanisms, especially those in immune regulatory mechanisms, between the two molecular subtypes. In addition, the constructed risk model was proven to be a well-established scoring system for predicting prognosis and evaluating the TIME status of each patient with PTC, which helped specify individualized treatment and follow-up strategy. Consequently, we speculate that cuproptosis plays an important role in tumor development, regulation of the TIME status, and prognosis in PTC.
Further bioinformatics analysis revealed the enrichment of many invasion-migration-associated pathways. These pathways, including biological adhesion, cell migration, locomotion, leukocyte migration, cell adhesion molecules, focal adhesion, and ECM-receptor interaction, have been widely confirmed to be associated with the prognosis of patients with PTC . In addition, some immune-related pathways, such as the positive regulation of immune system processes and MHC protein complexes, were also enriched. Therefore, cuproptosis may inhibit the progression and improve the prognosis of PTC by regulating invasion-migration-associated pathways and mediating these immune regulation processes.
The tumor microenvironment consists primarily of several different types of immune cells, extracellular matrix, and stromal cells, which provide tumor cells with support and nourishment. Tumor biology highly depends on the tumor microenvironment, especially the immune microenvironment . Copper, which is essential for the immune response, plays a vital role in both cellular and humoral immunity [36, 37]. Various immune cells, including B cells, T helper cells, macrophages, and NK cells, can be manipulated by copper to activate and maintain the immune system . In the present study, multiple immune analysis algorithms were applied to comprehensively evaluate the relationship between TIME status and the P-DECRG-based subtypes and the risk model. The results showed that the TIME features, as well as the abundance of many immune infiltrating cells, varied substantially across the two distinct P-DECRG-based molecular subtypes and various P-DECRG risk scores. Although there were some inconsistencies between the results of different immune algorithms, the significantly different trend of Tregs was always observed. Tumor-infiltrating Tregs inhibit antitumor immunity and promote cancer progression, with poor clinical outcomes . Recent studies have demonstrated that the CRG signature correlates with the infiltration of Tregs in osteosarcoma  and melanoma . Therefore, cuproptosis may improve the response to immunotherapy by inhibiting Treg infiltration into the tumor microenvironment of PTC. Based on current knowledge, immunotherapy, which has witnessed breakthrough progress in recent years, has provided a new strategy for tumor treatment. Among the various types of immunotherapies, immune checkpoint blockade therapy aims to inhibit or block the interaction of certain immune checkpoints with their ligands . Experimental and clinical studies of immune checkpoint inhibitors, including those targeting programmed cell death protein-1, programmed cell death-Ligand 1 (PD-L1), and cytotoxic T lymphocyte-associated antigen-4, have demonstrated improved treatment outcomes. Intra-tumoral copper supplementation enhances PD-L1 expression, affects the number of infiltrated CD8+ T and NK cells, and regulates key cancer immune evasion signaling pathways driven by PD-L1 . In the present study, significant differences were observed in the receipt of immune checkpoint blockade therapy between these two subtypes, and immune checkpoint blockade therapy could be a favorable choice for treating certain PTC subtypes.
Genome-wide CRISPR-Cas9 loss-of-function screens have been reported, and 10 genes closely related to cuproptosis have been identified . In the current study, we included a large number of genes, including 135 genes associated with copper homeostasis and transport, to comprehensively evaluate the influence of copper on the progression and prognosis of PTC. After stepwise screening and analysis, four specific P-DECRGs (CDK1, FDX1, MAP1LC3A, and P2RX4) related to the prognosis and TIME status of PTC were identified, which was not completely consistent with the results of a previously published study . Several explanations, including differences in tumor types, may explain this discrepancy.
FDX1 encodes a small iron-sulfur (Fe-S) protein that is present in the matrix of human mitochondria. In addition to its function in mitochondrial Fe-S cluster biogenesis, FDX1 is considered a versatile electron mediator involved in various physiological and pathological processes [45,46,47]. For example, in lung adenocarcinoma, FDX1 has been reported to be closely involved in fatty acid oxidation, glucose, and amino acid metabolism . However, in PTC, the functions and mechanisms of FDX1 have not been established. Nevertheless, the results of this study indicate that FDX1 may be negatively associated with the progression and prognosis of PTC through the mediation of cuproptosis.
CDK1 encodes a serine/threonine kinase that functions as a cell cycle checkpoint protein and plays a key role in controlling eukaryotic cell cycle progression . Binding to different cyclin proteins, CDK1 is considered the most critical cell cycle element and is sufficient to regulate numerous steps in cell proliferation and organ development . Recently, the cell cycle-independent function of CDK1 in enhancing overall protein synthesis has been uncovered . An integrative human pan-cancer analysis showed that dysregulation of CDK1 was identified in more than 20 human tumors and was significantly correlated with the development, progression, and microenvironment of tumor cells . Low-dose radiation exposure alters the expression of CDK1, which is associated with post-Chernobyl PTC development . In addition, CDK1 mediates the proliferation of PTC cells induced by high iodine . Our results in this study were consistent with those of other studies. Currently, the role of CDK1 in PTC progression remains unclear. We speculated that CDK1-related cuproptosis might be associated with energy metabolism reprogramming in PTC cells due to the vital function of CDK1 in regulating mitochondrial bioenergetics .
P2RX4 encodes a receptor protein of the P2X4 ATP-gated nonselective ionotropic channel, which has emerged as a key molecule controlling the release of neuronal brain-derived neurotrophic factor and is involved in pain processing . Emerging evidence has uncovered the role of P2RX4 in cancer biology, including cancer pain . In prostate cancer, P2RX4 is involved in enhancing tumor formation, mobility, TGFβ-1 induced invasiveness, and epithelial-to-mesenchymal transition [58, 59]. In hepatocellular carcinoma, P2RX4 is highly expressed and positively related to the growth and proliferation of cancer cells . Moreover, P2RX4 may play a vital role in regulating tumor development through inflammation and immune responses in the microenvironment . In the present study, P2RX4 expression was positively associated with the progression and prognosis of PTC by regulating the ionotropic channel related to cuproptosis.
Based on current knowledge, autophagy plays a crucial role in intracellular copper transportation . Microtubule-associated protein 1A/1B light chain 3 (LC3) is a well-known biomarker for evaluating autophagy. Three LC3 proteins exist in humans: LC3A, LC3B, and LC3C . Among these three members, MAP1LC3A encodes LC3A, which functions as a structural protein involved in mediating physical interactions between microtubules and other cytoskeleton elements. MAP1LC3A expression is downregulated in multiple tumors , and in the current study, MAP1LC3A expression was also suppressed in PTC, suggesting that its decreased expression may be related to the development of PTC cells. However, the exact molecular roles of MAP1LC3A in PTC pathogenesis require further study since the interactions between copper metabolism and autophagy are mostly unexplored.
The current study has some limitations. In the current study, no prognostic information on the selected patients with PTC was obtained due to the limited follow-up time. The number of patients with PTC selected for validation was small. Although the signatures of the four specific P-DECRGs were observed to be potential therapeutic targets or predictive biomarkers for PTC, the exact molecular mechanisms underlying their functions in cuproptosis have not yet been clarified. Therefore, further studies should include an expanded enrollment of cases and long-term follow-up to construct an independent validation cohort and confirm the predictive effect of the constructed risk model. In addition, further research is necessary to investigate the functions and molecular mechanisms of the four specific P-DECRGs in PTC.
In summary, in the present study, we identified two PTC molecular subtypes based on CRGs via consensus clustering. We also constructed a prognostic related predictive risk model using the signature of the four specific P-DECRGs, which was closely associated with the TIME status of patients with PTC. Our findings may provide valuable clinical guidance for either prognostic prediction or the development of precision treatment strategies for patients with PTC.
Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.
Papillary thyroid carcinoma
The Cancer Genome Atlas
Tumor immune microenvironment
Quantitative real-time polymerase chain reaction
Molecular Signatures Database
Principal component analysis
Differentially expressed genes
Kyoto Encyclopedia of Genes and Genomes
False discovery rate
Immune checkpoint blockade
Least absolute shrinkage and selection operator
Receiver operating characteristic
Area under the curve
Regulatory T cells
Major histocompatibility complex
Li M, Dal Maso L, Vaccarella S. Global trends in thyroid cancer incidence and the impact of overdiagnosis. Lancet Diabetes Endocrinol. 2020;8(6):468–70.
Li M, Brito JP, Vaccarella S. Long-Term Declines of Thyroid Cancer Mortality: An International Age-Period-Cohort Analysis. Thyroid. 2020;30(6):838–46.
Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin. 2021;71(1):7–33.
Wang J, Yu F, Shang Y, Ping Z, Liu L. Thyroid cancer: incidence and mortality trends in China, 2005–2015. Endocrine. 2020;68(1):163–73.
Megwalu UC, Moon PK. Thyroid Cancer Incidence and Mortality Trends in the United States: 2000–2018. Thyroid. 2022;32(5):560–70.
Ito Y, Miyauchi A, Kihara M, Fukushima M, Higashiyama T, Miya A. Overall Survival of Papillary Thyroid Carcinoma Patients: A Single-Institution Long-Term Follow-Up of 5897 Patients. World J Surg. 2018;42(3):615–22.
Liu C, Xiao C, Chen J, Li X, Feng Z, Gao Q, Liu Z. Risk factor analysis for predicting cervical lymph node metastasis in papillary thyroid carcinoma: a study of 966 patients. BMC Cancer. 2019;19(1):622.
Babak MV, Ahn D. Modulation of Intracellular Copper Levels as the Mechanism of Action of Anticancer Copper Complexes: Clinical Relevance. Biomedicines. 2021;9(8):852.
Tsvetkov P, Coy S, Petrova B, Dreishpoon M, Verma A, Abdusamad M, Rossen J, Joesch-Cohen L, Humeidi R, Spangler RD, et al. Copper induces cell death by targeting lipoylated TCA cycle proteins. Science. 2022;375(6586):1254–61.
Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA, Jacobsen A, Byrne CJ, Heuer ML, Larsson E, et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2012;2(5):401–4.
Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, Sun Y, Jacobsen A, Sinha R, Larsson E, et al. Integrative analysis of complexcancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013;6(269):pl1.
Liberzon A, Subramania A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27(12):1739–40.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010;26(12):1572–3.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7): e47.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.
Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–51.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, Trevino V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017;18(1):248–62.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.
Finotello F, Mayer C, Plattner C, Laschober G, Rieder D, Hackl H, Krogsdam A, Loncova Z, Posch W, Wilflingseder D, et al. Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. Genome Med. 2019;11(1):34.
Miao YR, Zhang Q, Lei Q, Luo M, Xie GY, Wang H, Guo AY. ImmuCellAI: A Unique Method for Comprehensive T-Cell Subsets Abundance Prediction and its Application in Cancer Immunotherapy. Adv Sci (Weinh). 2020;7(7):1902880.
Iasonos A, Schrag D, Raj GV, Panageas KS. How to build and interpret a nomogram for cancer prognosis. J Clin Oncol. 2008;26(8):1364–70.
Wang L, Huang Y, Liu C, Guo M, Ma Z, He J, Wang A, Sun X, Liu Z. Deltex3 inhibits Epithelial Mesenchymal Transition in Papillary Thyroid Carcinoma via promoting ubiquitination of XRCC5 to regulate the AKT signal pathway. J Cancer. 2021;12(3):860–73.
Haugen BR, Alexander EK, Bible KC, Doherty GM, Mandel SJ, Nikiforov YE, Pacini F, Randolph GW, Sawka AM, Schlumberger M, et al. 2015 American Thyroid Association Management Guidelines for Adult Patients with Thyroid Nodules and Differentiated Thyroid Cancer: The American Thyroid Association Guidelines Task Force on Thyroid Nodules and Differentiated Thyroid Cancer. Thyroid. 2016;26(1):1–133.
Pe’er D, Ogawa S, Elhanani O, Keren L, Oliver TG, Wedge D. Tumor heterogeneity. Cancer Cell. 2021;39(8):1015–7.
Dagogo-Jack I, Shaw AT. Tumour heterogeneity and resistance to cancer therapies. Nat Rev Clin Oncol. 2018;15(2):81–94.
Gawin M, Kurczyk A, Stobiecka E, Fratczak K, Polanska J, Pietrowska M, Widlak P. Molecular Heterogeneity of Papillary Thyroid Cancer: Comparison of Primary Tumors and Synchronous Metastases in Regional Lymph Nodes by Mass Spectrometry Imaging. Endocr Pathol. 2019;30(4):250–61.
Li J, Xie L, Xie Y, Wang F. Bregmannian consensus clustering for cancer subtypes analysis. Comput Methods Programs Biomed. 2020;189:105337.
Oliveri V. Selective Targeting of Cancer Cells by Copper Ionophores: An Overview. Front Mol Biosci. 2022;9:841814.
Chen B, Zhou X, Yang L, Zhou H, Meng M, Zhang L, Li J. A Cuproptosis Activation Scoring model predicts neoplasm-immunity interactions and personalized treatments in glioma. Comput Biol Med. 2022;148:105924.
Li Z, Zhang H, Wang X, Wang Q, Xue J, Shi Y, Wang M, Wang G, Zhang J. Identification of cuproptosis-related subtypes, characterization of tumor microenvironment infiltration, and development of a prognosis model in breast cancer. Front Immunol. 2022;13:996836.
Zhang Z, Zeng X, Wu Y, Liu Y, Zhang X, Song Z. Cuproptosis-Related Risk Score Predicts Prognosis and Characterizes the Tumor Microenvironment in Hepatocellular Carcinoma. Front Immunol. 2022;13:925618.
Reyes I, Reyes N, Suriano R, Iacob C, Suslina N, Policastro A, Moscatello A, Schantz S, Tiwari RK, Geliebter J. Gene expression profiling identifies potential molecular markers of papillary thyroid carcinoma. Cancer Biomark. 2019;24(1):71–83.
Pansy K, Uhl B, Krstic J, Szmyra M, Fechter K, Santiso A, Thuminger L, Greinix H, Kargl J, Prochazka K, et al. Immune Regulatory Processes of the Tumor Microenvironment under Malignant Conditions. Int J Mol Sci. 2021;22(24):13311.
Djoko KY, Ong CL, Walker MJ, McEwan AG. The Role of Copper and Zinc Toxicity in Innate Immune Defense against Bacterial Pathogens. J Biol Chem. 2015;290(31):18954–61.
Hackler J, Heller RA, Sun Q, Schwarzer M, Diegmann J, Bachmann M, Moghaddam A, Schomburg L. Relation of Serum Copper Status to Survival in COVID-19. Nutrients. 2021;13(6):1898.
Su Y, Zhang X, Li S, Xie W, Guo J. Emerging Roles of the Copper-CTR1 Axis in Tumorigenesis. Mol Cancer Res. 2022;20(9):1339–53.
Goswami TK, Singh M, Dhawan M, Mitra S, Emran TB, Rabaan AA, Mutair AA, Alawi ZA, Alhumaid S, Dhama K. Regulatory T cells (Tregs) and their therapeutic potential against autoimmune disorders - Advances and challenges. Hum Vaccin Immunother. 2022;18(1):2035117.
Yang M, Zheng H, Xu K, Yuan Q, Aihaiti Y, Cai Y, Xu P. A novel signature to guide osteosarcoma prognosis and immune microenvironment: Cuproptosis-related lncRNA. Front Immunol. 2022;13:919231.
Lv H, Liu X, Zeng X, Liu Y, Zhang C, Zhang Q, Xu J. Comprehensive Analysis of Cuproptosis-Related Genes in Immune Infiltration and Prognosis in Melanoma. Front Pharmacol. 2022;13:930041.
Petitprez F, Meylan M, de Reynies A, Sautes-Fridman C, Fridman WH. The Tumor Microenvironment in the Response to Immune Checkpoint Blockade Therapies. Front Immunol. 2020;11:784.
Voli F, Valli E, Lerra L, Kimpton K, Saletta F, Giorgi FM, Mercatelli D, Rouaen JRC, Shen S, Murray JE, et al. Intratumoral Copper Modulates PD-L1 Expression and Influences Tumor Immune Evasion. Cancer Res. 2020;80(19):4129–44.
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 (Basel). 2022;13(5):851.
Cai K, Tonelli M, Frederick RO, Markley JL. Human Mitochondrial Ferredoxin 1 (FDX1) and Ferredoxin 2 (FDX2) Both Bind Cysteine Desulfurase and Donate Electrons for Iron-Sulfur Cluster Biosynthesis. Biochemistry. 2017;56(3):487–99.
Wang Z, Dong H, Yang L, Yi P, Wang Q, Huang D. The role of FDX1 in granulosa cell of Polycystic ovary syndrome (PCOS). BMC Endocr Disord. 2021;21(1):119.
Tsvetkov P, Detappe A, Cai K, Keys HR, Brune Z, Ying W, Thiru P, Reidy M, Kugener G, Rossen J, et al. Mitochondrial metabolism promotes adaptation to proteotoxic stress. Nat Chem Biol. 2019;15(7):681–9.
Zhang Z, Ma Y, Guo X, Du Y, Zhu Q, Wang X, Duan C. FDX1 can Impact the Prognosis and Mediate the Metabolism of Lung Adenocarcinoma. Front Pharmacol. 2021;12:749134.
Xiang C, Sun WH, Ke Y, Yu X, Wang Y. CDCA8 Contributes to the Development and Progression of Thyroid Cancer through Regulating CDK1. J Cancer. 2022;13(7):2322–35.
Bhattacharyya N, Gupta S, Sharma S, Soni A, Bagabir SA, Bhattacharyya M, Mukherjee A, Almalki AH, Alkhanani MF, Haque S, et al. CDK1 and HSP90AA1 Appear as the Novel Regulatory Genes in Non-Small Cell Lung Cancer: A Bioinformatics Approach. J Pers Med. 2022;12(3):393.
Haneke K, Schott J, Lindner D, Hollensen AK, Damgaard CK, Mongis C, Knop M, Palm W, Ruggieri A, Stoecklin G. CDK1 couples proliferation with protein synthesis. J Cell Biol. 2020;219(3):e201906147.
Liu X, Wu H, Liu Z. An Integrative Human Pan-Cancer Analysis of Cyclin-Dependent Kinase 1 (CDK1). Cancers (Basel). 2022;14(11):2658.
Handkiewicz-Junak D, Swierniak M, Rusinek D, Oczko-Wojciechowska M, Dom G, Maenhaut C, Unger K, Detours V, Bogdanova T, Thomas G, et al. Gene signature of the post-Chernobyl papillary thyroid cancer. Eur J Nucl Med Mol Imaging. 2016;43(7):1267–77.
Lv C, Gao Y, Yao J, Li Y, Lou Q, Zhang M, Tian Q, Yang Y, Sun D. High Iodine Induces the Proliferation of Papillary and Anaplastic Thyroid Cancer Cells via AKT/Wee1/CDK1 Axis. Front Oncol. 2021;11:622085.
Xie B, Wang S, Jiang N, Li JJ. Cyclin B1/CDK1-regulated mitochondrial bioenergetics in cell cycle progression and tumor resistance. Cancer Lett. 2019;443:56–66.
Lalisse S, Hua J, Lenoir M, Linck N, Rassendren F, Ulmann L. Sensory neuronal P2RX4 receptors controls BDNF signaling in inflammatory pain. Sci Rep. 2018;8(1):964.
Zhang WJ, Luo C, Pu FQ, Zhu JF, Zhu Z. The role and pharmacological characteristics of ATP-gated ionotropic receptor P2X in cancer pain. Pharmacol Res. 2020;161:105106.
He J, Zhou Y, Arredondo Carrera HM, Sprules A, Neagu R, Zarkesh SA, Eaton C, Luo J, Gartland A, Wang N. Inhibiting the P2X4 Receptor Suppresses Prostate Cancer Growth In Vitro and In Vivo, Suggesting a Potential Clinical Target. Cells. 2020;9(11):2511.
Ghalali A, Ye ZW, Hogberg J, Stenius U. PTEN and PHLPP crosstalk in cancer cells and in TGFbeta-activated stem cells. Biomed Pharmacother. 2020;127:110112.
Asif A, Khalid M, Manzoor S, Ahmad H, Rehman AU. Role of purinergic receptors in hepatobiliary carcinoma in Pakistani population: an approach towards proinflammatory role of P2X4 and P2X7 receptors. Purinergic Signal. 2019;15(3):367–74.
Khalid M, Manzoor S, Ahmad H, Asif A, Bangash TA, Latif A, Jaleel S. Purinoceptor expression in hepatocellular virus (HCV)-induced and non-HCV hepatocellular carcinoma: an insight into the proviral role of the P2X4 receptor. Mol Biol Rep. 2018;45(6):2625–30.
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–24.
Baeken MW, Weckmann K, Diefenthaler P, Schulte J, Yusifli K, Moosmann B, Behl C, Hajieva P. Novel Insights into the Cellular Localization and Regulation of the Autophagosomal Proteins LC3A, LC3B and LC3C. Cells. 2020;9(10):2315.
Wang Z, Gao L, Guo X, Feng C, Lian W, Deng K, Xing B. Development and validation of a nomogram with an autophagy-related gene signature for predicting survival in patients with glioblastoma. Aging (Albany NY). 2019;11(24):12246–69.
This work was supported by the 345 Talent Project of Shengjing Hospital of China Medical University (30).
Ethics approval and consent to participate
The research protocol was approved by the Ethics Committee of the Shengjing Hospital of China Medical University (No. 2022PS547K). All the experimental protocols were performed in accordance with the Declaration of Helsinki. Written informed consent was obtained from all patients enrolled in this study. The datasets analyzed in this study were publicly available and, therefore, did not require ethics approval.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Wang, L., Yao, B., Yang, J. et al. Construction of a novel cuproptosis-related gene signature for predicting prognosis and estimating tumor immune microenvironment status in papillary thyroid carcinoma. BMC Cancer 22, 1131 (2022). https://doi.org/10.1186/s12885-022-10175-5