We obtained the clinical information and raw fragment per kilobase (FPKM) values of THCA in The Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) datasets from the UCSC XENA database. We directly downloaded the series matrix files of the Affymetrix microarray profiles for GSE29265, GSE33630, GSE35570, and GSE60542 from the Gene Expression Omnibus (GEO) database. All the information about the public datasets is available in Supplementary Table S1. We used the “Combat” algorithm from the R package “sva” to eliminate the batch effects from different GEO datasets (Johnson et al., 2007), and we used the principal component analysis (PCA) to validate the impact of data normalization further. We named the aggregated four GEO datasets as GEO merge data groups for further research.

For biological function enrichment analysis, we employed a non-parametric and unsupervised method, ssGSEA, to estimate the activation of specific pathways based on gene signature “c5.all.v6.2.symbols” from the MSigDB database (Subramanian et al., 2005). Similarly, the relative abundance of immune cell infiltration levels in the tumor microenvironment (TME) was also conducted using the ssGSEA algorithm based on immune cell signature (Supplementary Table S2) (Bindea et al., 2013). Furthermore, the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 (Huang et al., 2009) was also employed to conduct GO analysis, consisting of biological processes, cellular component, molecular function, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis.

To explore the underlying biological function of cuproptosis regulators in THCA, we clustered the THCA samples into several modules based on the expression patterns. Firstly, we used sample clustering to detect outliers in TCGA and GEO merge sets. The soft thresholding power was set as 9 for subsequent co-expression module establishment in two datasets. Module–trait associations were applied to set up a relationship between modules and the expression of cuproptosis. The functional enrichment analysis was conducted for GO and KEGG analysis for each expression module, respectively. All the intersection analyses were performed online, and the WGCNA algorithm (Langfelder and Horvath, 2008) was screened using R.

Genes were first filtered by correlation with both FDX1 expression and the predicted IC50 values of the cuproptosis inducer elesclomol in TCGA and GEO cohorts (|r| > 0.3, P < 0.05). Only genes satisfying these criteria in both datasets were retained and defined as candidate cuproptosis-related genes. These intersected genes were subsequently subjected to univariate Cox regression analysis, and genes with P < 0.05 were selected for feature selection. A Sliding Windows Sequential Forward Feature Selection (SWSFS) framework was then implemented using a Ranger-based random survival forest model (Zhang et al., 2019). Briefly, features were sequentially added in a forward manner, and at each step, model performance was evaluated using the out-of-bag (OOB) error as an internal cross-validation metric. The Ranger model was constructed with 1,000 trees, log-rank splitting rules, and a minimum node size of 15. The number of features yielding the lowest OOB error was considered optimal. Based on this procedure, six genes were selected to construct the final cuproptosis-related score (RS). The RS was calculated as: RS = Σ (Coefᵢ × Expᵢ), where Coefᵢ represents the regression coefficient and Expᵢ denotes the expression level of each gene. The median RS was used as the cutoff to stratify patients into high- and low-score groups. The cuproptosis RS could be calculated using the formula: cuproptosis RS = Σ (Coef i × Exp i), where i is the members involved in the gene signature. We performed stratification analysis to test whether the cuproptosis RS was an independent prognostic factor in THCA. Based on the members involved in cuproptosis RS, we built a prognostic nomogram using the “rms” R package to predict 1- and 3- years PFS of THCA, and the predictive accuracy of this nomogram was assessed using the calibration curve.

Normal thyroid cell line, Nthy-ori 3-1, and PTC cell lines, TPC-1 and BCPAP, were purchased from the National Collection of Authenticated Cell Cultures (Shanghai, China). All cell lines are maintained in RPMI1640 media supplemented with 10% FBS and 1% ampicillin/streptomycin and cultured at 5% CO₂.

Small interfering RNA was purchased from GenePharma (Suzhou, China) and transfected using Lipofectamine 3000 (Invitrogen/Thermo Fisher Scientific). Total RNAs from cells were extracted using the Tiangen DNA kit (Tiangen Biotech, Beijing, China) and measured total RNA by SYBR Green One-Step qRT-PCR kit (Invitrogen, 11736059). The specific details of primers and siRNA sequences in this study are shown in Supplementary Table S3.

Approximately 1000 cells were seeded into six-well plates in triplicate and incubated for 5 days. Then cells were washed with PBS and fixed with 4% paraformaldehyde for 15 min, followed by 30 min incubation with 0.1% crystal violet. After enough washing by PBS, cell colonies containing more than 50 cells were counted and photographed. Cells were plated in 96-well plates (2000 cells/well), and the cell counting kit-8 (CCK-8) was performed by the manufacturer’s protocol (Biosharp, China).

The comparison of different groups was estimated using the Mann–Whitney U test (Wilcoxon rank-sum test). Spearman’s correlation analysis was performed to calculate the correlation coefficient between the two factors. Receiver operator characteristic curve (ROC) analysis was used to calculate the elements’ area under curve (AUC) value. The optimal cutoff point for a factor was determined using the “survminer” R package, and the “surv-cutpoint” function was used to repeat all potential cutoff points to obtain the best separation groups. Survival curves for prognostic analysis were employed using the Kaplan-Meier method, and significant differences were determined using the log-rank test. They represent immunohistochemical (IHC) images of cuproptosis regulators that were downloaded from the Human Protein Atlas (HPA) (https://www.proteinatlas.org). To predict the half-maximal inhibitory concentration (IC50) of a cuproptosis inducer, elesclomol, in a single sample, the Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/) was used by R package “pRRophetic” following the instructions described previously (Yang et al., 2013). The Tumor Immune Estimation Resource (TIMER) analysis was conducted and visualized in Sangerbox database (http://sangerbox.com/). The asterisks represent the statistical P-value (*P < 0.05; **P < 0.01; ***P < 0.001).

The main flowchart is shown in Figure 1. To explore the intrinsic relationship among the cuproptosis regulators in THCA, correlation analysis was first displayed among ten cuproptosis regulators in the TCGA cohort (Figure 2A). CDKN2A, a cuproptosis negative regulator, was negatively correlated with other cuproptosis positive regulators, indicating that a competitive relationship existed among the cuproptosis regulators. To elucidate the importance of cuproptosis regulators in THCA, differential analysis was conducted in GTEx and TCGA datasets. Most cuproptosis positive regulators, such as DLD, FDX1, LIAS, LIPT1, and PDHA1, were significantly lower expressed in THCA than in normal thyroid tissues (Figure 2B). Furthermore, the matched-pair analysis showed that all cuproptosis positive regulators were substantially lower expressed in THCA, while CDKN2A was observably higher in THCA than in normal THCA tissues (Figure 2C). ROC analysis further elucidated that the expression of cuproptosis regulators could well reflect the differences between the THCA and para-carcinoma tissues, the AUC of most of them was more than 0.7 (Figure 2D). For validation, we combined four GEO datasets as a whole cohort. PCA showed that the normalization effects were pretty good (Figure 2E). Correlation analysis showed that three cuproptosis negative regulators, CDKN2A, MTF1, and GLS, were negatively correlated with cuproptosis positive regulators (Figure 2F), which was consistent with the regulatory patterns of cuproptosis. Differentially expressed analysis further confirmed that all cuproptosis positive regulators were lower in THCA tissues, while two cuproptosis negative regulators, CDKN2A and GLS, were highly expressed in THCA tissues (Figure 2G). ROC analysis confirmed again that FDX1 and CDKN2A were representative cuproptosis regulators in THCA, with AUC >0.8 (Figure 2H). Survival analysis showed that THCA patients with higher expression of FDX1, LIPT1, and PDHB owned better progression-free survival rates (Figure 3A), and THCA patients with higher expression of FDX1, LIAS, LIPT1 showed better disease-free survival rates (Figure 3B). The higher expression of CDKN2A simultaneously showed worse progression-free and disease-free survival probability in THCA. From the perspective of clinical traits, we found that CDKN2A and GLS were higher expressed in T4, N1, stage IV, and BRAF mutation groups. At the same time, DLD, FDX1, LIPT1, PDHA1, and PDHB were significantly lower expressed in T4, N1, stage IV, and BRAF mutation groups (Figure 3C). We validated again at the protein level with the help of an IHC image from the HPA database. GLS was higher in THCA tissues, and FDX1 was negative in THCA tissues but positive in normal thyroid tissues (Figure 3D). Moreover, we could see that DLD, LIAS, and LIPT1 were higher in normal thyroid tissues, CDKN2A was significantly higher in THCA, and others showed no apparent differences by qualitative analysis (Supplementary Figure S1).

For systematically analyzing the biological function of cuproptosis regulators in THCA, we employed a weighted gene co-expression network analysis (WGCNA) to construct a co-expression network. After removing the outlier samples, we analyzed the protein-coding RNAs with WGCNA to identify the modules containing highly correlated genes in TCGA and GEO merge datasets (Figures 4A,B). A soft threshold (β = 9) was used to guarantee a scale‐free network (Figures 4C,D), which identified eight modules and fourteen modules, respectively, from TCGA and GEO merge datasets (Figures 4E,F). Among these modules in the TCGA cohort, the blue module was significantly negatively correlated with cuproptosis positive regulators and positively correlated with cuproptosis negative regulators, while yellow, green, and brown modules were on the contrary (Figure 4G). As for the GEO merge cohort, brown, green-yellow, magenta and turquoise modules were also highly associated with cuproptosis regulators (Figure 4H). KEGG analysis showed that genes involved in the blue module were mainly enriched in cytokine-cytokine receptor interaction and chemokine signaling pathways, indicating that cuproptosis positive regulators might be highly negatively associated with immune-related pathways, such as T, B cells receptor and NK, TH17 cells differentiation pathways (Figure 5A). As for the other three modules positively correlated with cuproptosis positive regulators, Rap1, PI3K-AKT, Hippo, and MAPK signaling pathways associated with tumorigenesis and progression processes were highly enriched (Figure 5A). KEGG analysis in the GEO merge cohort verified the results of enrichment analysis from TCGA (Figure 5B), indicating that cuproptosis regulators might be involved in initiating and progressing tumorigenesis in THCA. For further research, we employed GSVA for validation in TCGA and GEO merge cohorts, respectively (Figures 5C,D). We found that most cuproptosis positive regulators positively correlated with protein secretion, oxidative phosphorylation, MYC, MTORC1, DNA repair, and adipogenesis while negatively correlated with immune-related pathways, such as interferon response and inflammatory pathways.

As FDX1 exhibited the most robust statistical associations among the ten cuproptosis-related genes, the above analyses highlight its potential value as a prognostic biomarker in THCA rather than implying a definitive functional dominance. We next chose FDX1 as a core analytic target. Correlation analysis showed that the expression level of FDX1 was significantly negatively correlated with the predicted IC50 of elesclomol, a well-known inducer of cuproptosis, by GDSC analysis in TCGA and GEO merge cohorts respectively (Figures 6A,B). Correlation analysis was conducted in TCGA and GEO merge cohorts, respectively, among the expression of protein-coding RNAs and the expression of FDX1, and the predicted IC50 of elesclomol. After the intersection among these groups, 637 intersected genes were selected as cuproptosis-related genes for further analysis (Figure 6C). These genes were satisfied with the criterion of | correlation coefficients | > 0.3 and P–value <0.05 (Supplementary Table S4). GO analysis showed that most of these genes were located at mitochondrial and associated with cell metabolism processes. KEGG analysis showed that reactive oxygen species and oxidative phosphorylation pathways were highly enriched (Figure 6D). The STRING database depicted the specific protein interaction relationship among these genes, and the correlation networks between genes and pathways were also shown in Figure 6E. The metabolic pathways and mitochondrial-related phenotypes were mainly involved here. Then, univariate COX analysis was conducted for these genes to explore the prognostic values for THCA (Figure 6F; Supplementary Table S5). With the criterion of P–value <0.05, 105 genes were selected for further analysis.

Then, as shown in Figure 7A, we employed Ranger, a weighted version of random forest based on the SWSFS algorithm, to evaluate the importance of each cuproptosis-related gene. After processing with the SWSFS algorithm, the results showed that when the number of genes was 6, the out-of-bag (OOB) error rate was the lowest, indicating that its predictive ability was the strongest. Thus, the importance of these 6 genes is in the upper right of Figure 7A. And the coefficients of these 6 genes were shown in Figure 7B; Supplementary Table S6, and the cuproptosis-related score (RS) = −0.78044 * FDX1 - 0.36267 * RBPMS2 - 0.12573 * NMB +0.08012 * MAP1LC3A+ 0.35546 * GINM1 + 0.03373 * FAM229B. Survival analysis showed that THCA samples in the high cuproptosis RS group possessed a worse survival rate than those in low cuproptosis RS (Figure 7C). Moreover, the time-dependent AUC values of the cuproptosis RS for predicting 1-, 2- and 3-year progression-free survival rates in THCA exceeded 0.6 (Figure 7D). Compared with a single gene involved in cuproptosis RS as an aspect of PFS prediction, cuproptosis RS showed excellent predicted efficacy than all of them (Figure 7E). As for the DFS rate, the cuproptosis RS could also predict well (Figure 7F), and the AUC values exceeded 0.68 for the prediction of 1-, 2- and 3-year disease-free survival rates in THCA (Figure 7G). Stratified analysis showed that a high cuproptosis RS was correlated with dramatically worse PFS, regardless of whether the patient exhibited early- or advanced-stage THCA (Figure 7H). Likewise, we also found that regardless of THCA belonging to BRAF wild-type (WT) or BRAF mutation (MUT), the cuproptosis RS provided statistically significant PFS stratification (Figure 7H). To further facilitate the clinical use of the cuproptosis RS, a nomogram capable of predicting the 1- or 3-year survival probability of THCA patients was finally constructed (Figure 7I). The calibration curves at 1 and 3 years indicated good consistency between the prediction by the nomogram and actual PFS (Figure 7J).

As the cuproptosis RS displayed potent prognosis prediction efficacy, we wonder if there was any underlying biological function in the process of THCA. Firstly, a correlation analysis was conducted between the GSVA enrichment scores and cuproptosis RS (Figure 8A). We found that PI3K/AKT, oxidative phosphorylation, adipogenesis, MTORC1, and DNA repair pathways were significantly negatively correlated with cuproptosis RS. In contrast, inflammatory, Hedgehog, and KRAS signaling pathways are positively associated with cuproptosis RS. It indicated that the high cuproptosis RS group displayed higher growth and inflammatory signaling, but the low cuproptosis RS group possessed higher oxidative stress levels.

Because cuproptosis RS is strongly associated with immune-related pathways, immune cell infiltration levels were assessed by ssGSEA in THCA. We found that most immune cells correlated with cuproptosis RS (Figure 8B). Survival analysis showed that effector memory T cell (Tem) and Th2 cell were a risk factor for THCA, but Th17 cell was a protective factor for THCA (Figure 8C). Correlation analysis showed that cuproptosis RS was positively correlated with Tem and Th2 cells but negatively correlated with Th17 cells (Figure 8D). To validate the relationship between FDX1 and immune cells infiltration level, TIMER showed that FDX1 was significantly associated with CD8⁺ T cells and macrophages (Figure 8E).

To further elucidate the importance of genes involved in cuproptosis RS signature, we found that FAM229B, GINM1, and NMB were significantly higher expressed in THCA. At the same time, the expression of MAP1LC3A and RBPMS2 were markedly lower in THCA (Figure 9A). Using qRT-PCR, we found that the expression level of MAP1LC3A and RBPMS2 were dramatically upregulated in a normal thyroid cell line in comparison with. Still, the expression level of GINM1 was significantly higher in THCA cell lines (Figure 9B). The expression level of NMB and FAM229B was uncertain between THCA and the normal cell line. With the help of HPA, we confirmed that the expression of GINM1 was remarkably higher in THCA tissues, and MAP1LC3A was more elevated in normal tissues (Figure 9C). However, the expression level of NMB was both strong in THCA or normal thyroid tissues, and the IHC images of RBPMS2 and FAM229B remained uncertain. To investigate whether the biological function of FDX1 is related to the growth of THCA, the knockdown of FDX1 in cell lines was conducted using small interference RNA. Through qPCR, we found that the siRNAs significantly inhibited the expression of FDX1 (Figure 9D). Surprisingly, we found that the expression level of cuproptosis negative regulator GLS was elevated, while MAP1LC3A, RBPMS2, and GINM1 were downregulated (Figure 9D). Unexpectedly, the results of the colony-formation experiment and CCK-8 assays suggested that the knockdown of FDX1 could not affect the proliferation ability of thyroid cells (Figures 9E,F).