Molecular function validation and prognostic value analysis of the cuproptosis-related gene ferredoxin 1 in papillary thyroid carcinoma.

FDX1 functions as a key upstream regulator of the cuproptosis signaling, and its downregulation can lead to resistance against copper-induced cell death. To explore its expression status in PTC, Wilcoxon tests were used to ascertain the expression difference significance between primary PTC tissues and adjacent normal tissues. As shown in Fig.  1 A–D, PTC displayed a notably lower FDX1 expression compared to normal tissues in all datasets (TCGA-THCA, ENA_PRJEB11591, GSE33630 , GSE60542 , GSE50760 ). Besides, a single-cell data analysis of GSE191288 revealed that a total of 24,024 cells were grouped into 15 clusters. Cluster 0 was identified as thyroid follicular epithelial cells with high expression of biomarkers such as EPCAM, TSHR, and TG (Fig.  1 E,F). These thyroid follicular epithelial cells were further subdivided into 6 subclusters (Fig.  1 G). Among them, cluster 2 was identified as benignant cells with high expression of TPO, while the remaining clusters were identified as malignant cells characterized by high expression of MET and KRT19(Fig.  1 H,I). The feature plots indicated that FDX1 was significantly underexpressed in tumor cells when compared with normal thyroid follicular epithelial cells (Fig.  1 J). Besides, our own confirmed that the expression of FDX1 in PTC tumor tissues was significantly lower than that in normal thyroid tissues detected by PCR and immunohistochemistry (Fig.  6 A–C). Fig. 1 FDX1 was downregulated in PTC. ( A – D ) Boxplots showing FDX1 was significantly downregulated in PTC samples compared with normal samples. ( E , F ) Single-cell data analysis showed that a total of 24,024 cells were categorized into 15 clusters and cluster 0 was identified as thyroid follicular epithelial cells. ( G , H ) Thyroid follicular epithelial cells were further grouped into 7 subgroups by the UMAP algorithm. ( I , J ) Violin plots showing the expression level of the biomarkers. FDX1 was downregulated in PTC. ( A – D ) Boxplots showing FDX1 was significantly downregulated in PTC samples compared with normal samples. ( E , F ) Single-cell data analysis showed that a total of 24,024 cells were categorized into 15 clusters and cluster 0 was identified as thyroid follicular epithelial cells. ( G , H ) Thyroid follicular epithelial cells were further grouped into 7 subgroups by the UMAP algorithm. ( I , J ) Violin plots showing the expression level of the biomarkers. To explore the correlation of the expression level of FDX1 with the prognosis in PTC, KM survival analyses were conducted for DFS among PTC patients with complete survival data. As depicted in Fig.  2 E, the DFS analysis results from the TCGA-THCA dataset showed that patients with lower FDX1 expression had a higher recurrent rate. Additionally, Wilcoxon tests were performed to identify the relationship of FDX1 with tumor clinical features. The results showed that patients in the extrathyroidal extension (ETE) group or T3/T4 group had significantly lower expression of FDX1 (Fig.  2 A–D). Accordingly, we postulate that the dysregulation of the FDX1-related cuproptosis signaling pathway could significantly affect the prognosis of PTC patients. Fig. 2 Downregulated FDX1 indicated poor prognoses. ( A – C ) Boxplots showing FDX1 was significantly negatively correlated with tumor size, regional lymph node metastasis, and extrathyroidal extension. ( D ) Boxplots showing FDX1 was significantly negatively correlated with tumor stage. ( E ) The KM survival analysis showed that patients with lower FDX1 expression had a higher recurrent rate. Downregulated FDX1 indicated poor prognoses. ( A – C ) Boxplots showing FDX1 was significantly negatively correlated with tumor size, regional lymph node metastasis, and extrathyroidal extension. ( D ) Boxplots showing FDX1 was significantly negatively correlated with tumor stage. ( E ) The KM survival analysis showed that patients with lower FDX1 expression had a higher recurrent rate. Considering the vital role of FDX1 in the cuproptosis signaling pathway and its profound impact on the prognosis in PTC patients, we assumed that FDX1-related genes could potentially serve as indicators for assessing the status of cuproptosis and predicting PTC prognosis. To test this hypothesis, we first conducted Pearson correlation analyses across three datasets (TCGA-THCA, ENA_PRJEB11591, and GSE33630_GSE60542). Through these analyses, we identified 2876 FDX1-related genes based on a coefficient threshold of |r| > 0.1 and p  < 0.01 (Supplementary Table 1). Then, these genes were used for enrichment analysis, while the results of GO analysis revealed significant enrichment in the mitochondrial respiratory-related biological processes, such as cellular respiration, respiratory electron transport chain, and mitochondrial respiratory chain complex assembly (Fig.  3 A). Furthermore, the Reactome Gene Sets analysis showed that these genes were significantly enriched in the pathways of the citric acid (TCA) cycle and respiratory electron transport, respiratory electron transport (Fig.  3 B). Taken together, these findings strongly implicate the biological function of FDX1 and its related genes in the process of cuproptosis. Moreover, we identified 187 genes that are highly correlated with FDX1, based on a coefficient threshold of |r| > 0.5 and p  < 0.01 (Supplementary Table 2). Through univariate analyses, we discovered that 40 of these genes could statistically affect DFS prognosis (Fig.  3 C). Then, we incorporated these 40 genes into a LASSO regression model using the TCGA-THCA cohort data. This process led to the selection of six key genes: HIGD1A, SPCS2, PPIL1, CLCNKB, FHL1, and IQGAP3 (Fig.  3 D). As a result, these six genes were defined as cuproptosis-related genes in PTC and were selected to further construct the cuproptosis-related prognosis model. Fig. 3 Functional annotation of FDX1-related genes. ( A , B ) GO enrichment analysis and Reactome Gene Sets enrichment analysis of FDX1 and its related genes. ( C ) Univariate analyses showed that 40 out of FDX1-related genes could statistically affect DFS prognosis. ( D ) The LASSO Cox regression model was constructed with the 40 genes and six signature genes were eventually identified according to the best fit profile. Functional annotation of FDX1-related genes. ( A , B ) GO enrichment analysis and Reactome Gene Sets enrichment analysis of FDX1 and its related genes. ( C ) Univariate analyses showed that 40 out of FDX1-related genes could statistically affect DFS prognosis. ( D ) The LASSO Cox regression model was constructed with the 40 genes and six signature genes were eventually identified according to the best fit profile. The cuproptosis-related Cox proportional hazards model was established using the six carefully selected genes in the dataset TCGA-THCA. The CRRS of each patient was determined as outlined in the methods section. The KM survival analyses showed that patients with higher CRRS displayed a poorer prognosis (Fig.  4 A). The time-dependent ROC analyses indicated that the CRRS had a good prognosis prediction accuracy, surpassing that of clinical features, like age, sex, and tumor size within the TCGA-THCA cohort (Fig.  4 B). To test whether the CRRS statistically enhances the prognosis prediction accuracy of the traditional clinicopathological staging systems in PTC patients, we constructed a clinical prognosis model solely based on clinical features of the tumor stage. In addition, a combined prognosis model incorporated both clinical features and the CRRS. The risk score of each patient in both two models was calculated as described before (Fig.  4 B). These results from IDI analyses showed a statistically significant improvement in 1- and 3-year DFS prognosis prediction accuracy of the combined prognosis model compared to the clinical prognosis model (Supplementary Fig. 1). Significant prognosis differences were observed among the various risk stratification subgroups (Fig.  4 C). In the nomogram plot, weighted scores were derived from risk features, including tumor stage and CRRS, and these scores were applied to predict the 1–5-year DFS rate of PTC patients (Fig.  4 D). The calibration curve plot exhibited good performances of the nomogram plot in prognosis prediction accuracy when compared to an ideal model (Supplementary Fig. 2). Hence, our findings indicated that cuproptosis-related signatures could serve as a reliable prognosis predictor for PTC patients. Fig. 4 Development and validation of the CRRS model. ( A ) The KM survival analysis showed that patients with higher CRRS expression had a higher recurrent risk. ( B ) The time-dependent ROC analyses indicated that the CRRS had a good prognosis prediction accuracy for 1-, 3-, and 5-year recurrent rates compared with the clinical features. The combined prognosis model constructed with tumor stage and CRRS could statistically improve the prediction accuracy. ( C ) The KM survival analysis showed significant prognosis differences among the various risk stratification subgroups. ( D ) A Nomogram plot was developed with the two features of tumor stage and CRRS to predict the 1-, 3-, and 5-year recurrent rates of PTC patients. Development and validation of the CRRS model. ( A ) The KM survival analysis showed that patients with higher CRRS expression had a higher recurrent risk. ( B ) The time-dependent ROC analyses indicated that the CRRS had a good prognosis prediction accuracy for 1-, 3-, and 5-year recurrent rates compared with the clinical features. The combined prognosis model constructed with tumor stage and CRRS could statistically improve the prediction accuracy. ( C ) The KM survival analysis showed significant prognosis differences among the various risk stratification subgroups. ( D ) A Nomogram plot was developed with the two features of tumor stage and CRRS to predict the 1-, 3-, and 5-year recurrent rates of PTC patients. PTC patients were divided into two subgroups, high- and low-CRRS groups, based on the median of CRRS. To investigate differences in immune status within the tumor microenvironment between these subgroups, we initially performed the ESTIMATE algorithm to determine the stromal score and immune score of each patient. Our finding revealed that the high-CRRS group exhibited a notably higher proportion of immune infiltration and stroma components compared to the low-CRRS group in TCGA-THCA (Fig.  5 B). Moreover, ssGSEA results for 28 types of immune cells indicated the high-CRRS group had a higher proportion of immunosuppressive cells, such as Treg cells and MDSCs compared with the low-CRRS group in the TCGA cohort (Fig.  5 D). Furthermore, the boxplot of 12 types of immune checkpoints showed that 10 out of 12 were significantly overexpressed in the high-CRRS group (Fig.  5 C). Besides, GSEA was conducted on the TCGA-THCA cohort to identify the difference in biological processes between the two CRRS subgroups. The results indicated that the high-CRRS group was significantly enriched in the hallmarks of G2M checkpoint, mitotic spindle, epithelial-mesenchymal transition, and IL2 JAK STAT3 signaling (Fig.  5 A). These biological processes are well-recognized as intimate to tumorigenesis and progression. Therefore, our results suggested that the dysregulation of cuproptosis may induce alterations in the tumor microenvironment, thereby facilitating tumor progression. To evaluate the role of the six-genes signature in predicting prognosis and tumor microenvironment of other endocrine cancer contexts, we respectively performed the same bioinformatic analyses with breast carcinoma (BRCA) cohort and prostate adenocarcinoma (PRAD) cohort as we did in thyroid cancer (THCA). And, the results showed these six genes had similar expression differences between the tumor samples and the normal samples in THCA, BRCA, and PRAD. The risk score evaluated with these six genes had good performance in predicting prognosis in all three cohorts (Supplementary Fig. 3). However, the role of these six genes in the immune microenvironment of these three types of cancers had significant inconsistencies (Supplementary Fig. 4). Therefore, the findings in this study have certain specificity in THCA compared with BRCA and PRAD. Fig. 5 Immune landscape of CRRS-Based classification. ( A ) The GSEA plots showing the significantly enriched hallmarks of the high-CRRS group compared to the low-CRRS group. Violin plots showing the distribution of immune scores ( B ) and the expression difference of 12 types of immune checkpoints between the high- and low-CRRS subgroups. ( D ) Violin plots displaying the enrichment score distribution of 28 types of immune cells obtained by the method of ssGSEA. Immune landscape of CRRS-Based classification. ( A ) The GSEA plots showing the significantly enriched hallmarks of the high-CRRS group compared to the low-CRRS group. Violin plots showing the distribution of immune scores ( B ) and the expression difference of 12 types of immune checkpoints between the high- and low-CRRS subgroups. ( D ) Violin plots displaying the enrichment score distribution of 28 types of immune cells obtained by the method of ssGSEA. To explore the biological function of FDX1 in PTC, we constructed the FDX1-overexpressed (FDX1-OE) cell models utilizing the PTC cell lines of K1 and TPC-1. The results of the CCK8 assay showed that overexpression of FDX1 could significantly inhibit cell proliferation of K1 and TPC-1 (Fig.  6 D). In addition, the results of flow cytometry showed that overexpression of FDX1 could significantly promote cell death of K1 and TPC-1 (Fig.  6 E,F). The results of wound healing assays and transwell assays showed that overexpression of FDX1 could significantly inhibit the migratory and invasive capacities of K1 and TPC-1 cells (Supplementary Figs. 5 and 7 A-D). What’s more, overexpression of FDX1 also downregulated E-cadherin and upregulated N-cadherin and Vimentin, promoting epithelial-mesenchymal transition (EMT) (Fig.  7 E). Thus, our findings indicated that the downregulation of FDX1 may promote the progression and recurrence in PTC by inhibiting cell death and enhancing the viability of tumor cells. Fig. 6 Functional verification of FDX1 in vitro experiments. ( A – C ) The results of qPCR and IHC showed that FDX1 was significantly downregulated in PTC tissues compared with normal tissues. ( D ) CCK-8 assays showing the viability of K1 and TPC-1 cells with overexpressed FDX1. ( E , F ) Flow cytometry analysis with Annexin V-PI staining was performed to evaluate the percentage of apoptotic cells in K1 and TPC-1 with upregulated FDX1. Functional verification of FDX1 in vitro experiments. ( A – C ) The results of qPCR and IHC showed that FDX1 was significantly downregulated in PTC tissues compared with normal tissues. ( D ) CCK-8 assays showing the viability of K1 and TPC-1 cells with overexpressed FDX1. ( E , F ) Flow cytometry analysis with Annexin V-PI staining was performed to evaluate the percentage of apoptotic cells in K1 and TPC-1 with upregulated FDX1. Fig. 7 Overexpression of FDX1 could inhibit the viability of K1 and TPC-1 cells ( A , B ) Transwell assays showing the migration of K1 and TPC-1 cells with overexpressed FDX1. ( C , D ) Transwell assays showing the invasion of K1 and TPC-1 cells with overexpressed FDX1. ( E ) Western blot assay showed expression profiles of N-cadherin, E-cadherin, Vimentin, BRCA1, and RAD51, with GAPDH serving as a loading control, following FDX1-overexpression. Overexpression of FDX1 could inhibit the viability of K1 and TPC-1 cells ( A , B ) Transwell assays showing the migration of K1 and TPC-1 cells with overexpressed FDX1. ( C , D ) Transwell assays showing the invasion of K1 and TPC-1 cells with overexpressed FDX1. ( E ) Western blot assay showed expression profiles of N-cadherin, E-cadherin, Vimentin, BRCA1, and RAD51, with GAPDH serving as a loading control, following FDX1-overexpression. To investigate the regulatory mechanism of FDX1 in cuproptosis signaling, we utilized the online tool STRING to perform PPI network analysis, inputting FDX1 and its 187 most related genes. The PPI network comprised a total of 184 nodes and 281 edges (Fig.  8 A). Our finding indicated that the hub gene CYCS, with a maximum degree of 29, played a pivotal role in the regulatory network and exhibited a high interaction confidence score of 0.98 with FDX1 (Supplementary Table 3). Moreover, the results of Pearson correlation analysis showed a significantly positive association between FDX1 and CYCS (Fig.  8 B). Interestingly, CYCS, an electron carrier protein, holds an important role in apoptosis, aligning with our hypothesis regarding the molecular function of FDX1. We verified our findings in clinical samples, CYCS was significantly decreased in PTC tissues compared to normal ones (Fig.  8 C). In addition, the expression of CYCS was enhanced along with the overexpression of FDX1 in the cell lines of K1 and TPC-1 (Fig.  8 D). We observed that while FDX1 overexpression did not significantly alter the expression of the DNA repair gene BRCA1 , it markedly reduced RAD51 expression in TPC-1 cells (Fig.  7 E). This suggests that FDX1 may modulate DNA repair pathways, potentially implicating it in cuproptosis-related signaling. Fig. 8 Construction of the PPI network-based FDX1. ( A ) The PPI network comprised a total of 184 nodes and 281 edges and showed a high interaction confidence of FDX1 with CYCS. ( B ) The results of Pearson correlation analysis showed a significantly positive association between FDX1 and CYCS. ( C ) The results of qPCR showed that CYCS was significantly downregulated in PTC tissues compared with normal tissues. ( D ) The expression of CYCS was enhanced along with the overexpression of FDX1 in the cell lines of K1 and TPC-1. Construction of the PPI network-based FDX1. ( A ) The PPI network comprised a total of 184 nodes and 281 edges and showed a high interaction confidence of FDX1 with CYCS. ( B ) The results of Pearson correlation analysis showed a significantly positive association between FDX1 and CYCS. ( C ) The results of qPCR showed that CYCS was significantly downregulated in PTC tissues compared with normal tissues. ( D ) The expression of CYCS was enhanced along with the overexpression of FDX1 in the cell lines of K1 and TPC-1.

The mRNA expression profiles in Fragments Per Kilobase Million (FPKM) format along with corresponding clinical information for 497 primary thyroid cancer (THCA) patients, 1081 breast carcinoma (BRCA) patients, and 483 prostate adenocarcinoma (PRAD) patients were downloaded from The Cancer Genome Atlas TCGA-THCA program ( https://portal.gdc.cancer.gov/ ). Additionally, we obtained normalized gene expression matrixes of the microarray datasets ( GSE33630 , GSE60542 , and GSE83520 ) with their corresponding clinical information from Gene Expression Omnibus (GEO) ( https://www.ncbi.nlm.nih.gov/geo/ ). The raw sequencing data of project PRJEB11591, comprising 125 PTC samples and 65 normal samples, were retrieved from the European Nucleotide Archive (ENA). Furthermore, the normalized single-cell sequencing expression profiles of GSE191288 including 6 PTC samples and 1 normal sample were also downloaded from GEO. These data underwent additional quality control and filtering using the R package “Seurat”. Given the pivotal function of FDX1 in the cuproptosis signaling pathway, the Wilcoxon rank-sum test was applied to compare the expression levels of FDX1 between the primary PTC tissues and normal thyroid tissues across multiple datasets, including TCGA-THCA, GSE33630 , GSE60542 and PRJEB11591, respectively. Additionally, the single-cell sequencing dataset GSE191288 was utilized to assess the expression difference of FDX1 between malignant thyroid cells and normal thyroid cells. Subsequently, Pearson correlation analyses were performed on PTC samples from the cohorts of TCGA-THCA, GSE33630 , GSE60542 , and PRJEB11591, to pinpoint FDX1-related genes. Genes exhibiting correlation coefficients |r| > 0.5 and p  < 0.01, along with FDX1, were recognized as cuproptosis-related genes. To identify the cuproptosis-related genes that are significantly associated with disease-free survival (DFS) in PTC, univariate analyses were conducted, with a significance threshold set at p  < 0.05. Subsequently, the cuproptosis-related genes associated with DFS prognosis were further screened via the least absolute shrinkage and selection operator (LASSO) Cox regression analysis in the TCGA-THCA dataset via the R package “glmnet”. Finally, these selected genes were utilized to construct the cuproptosis-related Cox proportional-hazards model with the R package “survival” and “survminer”. The cuproptosis-related risk score (CRRS) of each patient was then computed using “predict” function from the R package “stats” based on this model. PTC patients were divided into two subgroups, high- and low-CRRS groups, based on the median value of CRRS. To delve into the disparities in biological characteristics between these two groups, GSEA software (version 4.2.3) was installed to conduct gene set enrichment analysis (GSEA) with the hallmark gene sets. These sets encapsulate and represent specific well-defined biological processes and display coherent expression obtained from the Human Molecular Signatures Database (MSigDB). Furthermore, to explore the molecular function of FDX1, a Gene Ontology (GO) enrichment analysis and Reactome Gene Sets analysis were executed for FDX1 and its strongly correlated gene (|r| > 0.1 and p  < 0.01) via the online tool Metascape ( https://metascape.org ). To evaluate the immune landscape difference in the two CRRS-based subgroups of PTC patients, we initially calculated the immune score and stromal score of each sample using the R package “estimate”. Then, the single sample Gene Set Enrichment Analysis (ssGSEA) with the signatures of 28 types of immune cells, to determine the proportion of these cells within the tumor microenvironment of each sample via the R package “GSVA”. Moreover, we compared the expression levels of 12 types of immune checkpoint-related genes between the two subgroups. To explore the functional mechanism of FDX1 in PTC, the web tool STRING ( https://cn.string-db.org/ ) was used to generate a protein-protein interaction (PPI) network with FDX1 and its associated genes (|r| > 0.1 and p  < 0.01). A combined score was calculated to evaluate the protein-protein interaction probability. The Cytoscape software was utilized to visualize the network, and the node degree was used to identify the hub genes in the network. To compare the DFS prognosis prediction accuracy of the CRRS versus clinical features, we first constructed the cuproptosis-related prognosis model. The CRRS of each patient was calculated as described above. Subsequently, a clinical prognosis model was constructed based on clinical data including age, sex, and tumor TNM stage. Time-dependent ROC analyses were performed to individually estimate the survival prediction accuracy of these features separately. Finally, we developed a combined prognosis model that constructed with the cuproptosis-related genes and the clinical features. The index of the Integrated Discrimination Improvement (IDI) was calculated using the R package “survIDINRI” to assess the improvement significance in prognostic accuracy by the CRRS. The combined prognosis model was constructed as described above and the risk features incorporated in the model were then used to calculate the risk score of each patient. Kaplan-Meier (KM) survival analysis was used to identify the prognosis difference significance in overall survival (OS) between the high and low-risk groups based on the median risk scores. Moreover, a nomogram plot featuring weighted scores, calculated using the risk features from the combined model, was devised to predict the 1–5 years OS of PTC patients. Finally, calibration curves were generated to compare the accuracy of the survival prediction made by the combined model versus an ideal model. This study was ethically approved by the Medical Ethics Committee of Xiangya Hospital, Central South University. The fixed clinical samples were embedded in paraffin. Three-micron sections were placed on glass slides, dewaxed in xylene and hydration in 100%, 95%, 85%, and 75% ethyl alcohol and washed in PBS three times, repaired in 90 ℃ sodium citrate for 10 min, and then blocked with PBS containing 0.3% Triton-X-100, 3% hydrogen peroxide and 10% normal goat serum for 30 min at room temperature. Sections were then incubated with the following primary antibodies at the indicated dilution in blocking buffer overnight at 4℃: FDX1( 12592-1-AP, Proteintech, Wuhan, China). Sections were then washed in PBS 3 times before incubating in species-specific anti-IgG secondary antibodies diluted 1:500 in blocking buffer for 1 h at room temperature. Sections were then washed in PBS 3 times, then DAB color rendering under a microscope. Total RNA was extracted using the TransZol Up Plus RNA Kit (Cat. No. ER501-01; Transgen) and cDNA was synthesized from 1 µg of total RNA using the NovoSceipt Plus All‐in‐one 1st Strand cDNA synthesis kit (Cat. No. E047; Novoprotein, Shanghai, China). The cDNA was then amplified with 2× SYBR Green qPCR Master Mix (Cat. No. B21202 ; Bimake, Texas, USA) in an FTC‐3000 real‐time PCR system (Funglyn Biotech Inc, Toronto, Canada). The relative standard curve method (2 − ΔΔCt) was used to determine the relative gene expression and Gapdh was used as a housekeeping gene for internal normalization. The PCR primers used in this study were as follows: human-FDX1: forward, 5″-TTCAACCTGTCACCTCATCTTTG‐3″, and reverse, 5″‐TGCCAGATCGAGCATGTCATT‐3″; human-CYCS: forward, 5″‐CTTTGGGCGGAAGACAGGTC‐3″, and reverse, 5″‐TTATTGGCGGCTGTGTAAGAG‐3″. PTC cell lines of TPC-1 and K1 both purchased from the American Type Culture Collection (ATCC), these cells both incubated in RPMI 1640 (VivaCell Biosciences, Shanghai, China) containing 10% fetal bovine serum (FBS; Gibco, Grand Island, USA) and 1% penicillin-streptomycin (PS; Solarbio, Beijing, China). Cells were maintained at 37 °C in a 5% CO2 humidified atmosphere. FDX1 overexpression plasmid was synthesized by Youbao (Youbio, Hunan, China). In brief, TPC-1 or K1 were transfected with FDX1 overexpression plasmid, supplemented by DNA transfection reagent (Cat. No. TF20121201, Neofect, Beijing, China), incubated at 37 °C for 24 h. A cell counting kit-8 (CCK-8) assay was used to determine cell viability. In total, 5 × 103 cells were seeded in 96-well plates overnight. After transfection with plasmid overexpression with FDX1 for 24 h,48 h, 72 h,96 h,10 µL CCK-8 dye was added to each well, and the cells were incubated for 1 h at 37 °C. Then, the absorbance was determined at 450 nm. Apoptosis was assessed in K1 and TPC-1 using FITC annexin V apoptosis detection kit with propidium iodide (BioLegend) following the supplier’s protocol. Cells were centrifuged at 1,000 r/min for 5 min and then the supernatant was discarded. Subsequently, the cells were resuspended and counted. A total of 100,000 cells were taken and resuspended with 90 µl of Annexin V-FITC binding solution, followed by addition of 5 µl of Annexin V-FITC reagent and 5 µl of propidium iodide. The cells were gently mixed with the reagent, and placed at room temperature for 15 min in the dark. Finally, the sample was detected by a flow cytometer. Wound healing assays were conducted to evaluate the impact of FDX1 overexpression on the migratory capacity of K1 and TPC-1 cells. The cells were seeded in culture-insert wells with 6-well plates and allowed to reach confluence. Upon achieving a confluent monolayer, uniform scratches were introduced to create consistent cell-free gaps. The culture medium was then replaced with fresh medium containing 2% FBS to minimize proliferation-driven interference. Cell migration was monitored at 24-hour intervals using bright-field microscopy, and images were systematically captured for quantitative analysis. Transwell chambers (8-uM pore size; Labselect) were employed to assess cell migration in vitro. Briefly, 200 µl of serum-free cell suspension (1 × 10^6 cells/ml of TPC-1 and K1) was seeded into the upper chamber (chambers for migration assays require no pre-treatment, while invasion assays necessitate Matrigel pre-treatment). The bottom well was filled with 600 µl of complete medium supplemented with 10% FBS as a chemoattractant. Following incubation, no-migratory cells on the upper surfaces were gently removed using cotton swabs. Migration cells on the lower surfaces were fixed with 4% paraformaldehyde, stained with 0.5% crystal violet, and quantified. For each membrane, at least three randomly selected fields were imaged at 10× magnification and migrated cells were evaluated. Cells were lysed in RIPA buffer to extract total proteins. Protein concentrations were determined by BCA assay, and equal amounts of protein were denatured in SDS sample buffer by boiling at 95 °C for 5 min. The denatured proteins were separated by SDS-PAGE and subsequently transferred onto PVDF membranes. The membranes were incubated with an anti-N-cadherin antibody (1:2000, Proteintech, Wuhan, China), anti-E-cadherin antibody (1:2000, Proteintech, Wuhan, China), anti-Vimentin antibody (1:3000, Proteintech, Wuhan, China), anti-BRCA1 antibody (1:2000, ABclonal, Wuhan, China), anti-RAD51 antibody (1:2000, ABclonal, Wuhan, China)at 4 °C overnight, followed by washing with TBST (Tris-Buffered Saline with Tween 20). Subsequently, the membranes were incubated with a corresponding secondary anti-rabbit horseradish peroxidase-conjugated antibody (1:5000, Solarbio, Beijing, China) for 1.5 h at room temperature. Antibody-antigen complexes were detected using a chemiluminescent ECL reagent. In this study, all statistical analyses were conducted in R software (version 4.2.1) with appropriate R packages and GraphPad Prism 9 (version 9.2.0). Wilcoxon rank-sum tests and unpaired, two-tailed Student’s t-test with a confidence level of 95% were applied to compare the mean values between two groups, and Kruskal-Wallis tests were applied for three or more groups. We confirm that all methods were performed in accordance with the relevant guidelines and regulations. The threshold of statistical significance was set at p  < 0.05. The p values were displayed as * ( p  < 0.05), ** ( p  < 0.01), *** ( p  < 0.001), **** ( p  < 0.0001), and ns for not significant.

PTC remains the most prevalent form of thyroid cancer and is currently one of the most common malignancies. However, approximately 20% of PTC patients still face a poor prognosis due to distant metastasis and regional recurrence 4 , 5 . Currently, the prognosis of PTC patients is predominantly determined by clinical and histopathologic parameters, including tissue type, primary tumor size, glandular invasion, vascular invasion, BRAF mutation, distant metastasis, and several others 14 – 16 . Numerous researchers have identified different prognostic markers, gene signatures, and prediction models specific to PTC 17 – 19 . Recently, given the crucial role that cuproptosis plays in the development of malignancy, exploring the association between curproptosis and the prognosis of PTC is of utmost importance. Cupropotosis has been linked to various diseases in recent studies. It has been reported as an apoptosis mode related to the occurrence and development of tumors, and the copper concentration is significantly increased in kinds of tumor cells or serum (e.g. breast cancer, lung cancer, prostate cancer, kidney cancer, and colon cancer) 20 , 21 . What’s more, cuproptosis also takes part in Wilson disease, Menkes disease, obesity, neurogenerative disease, and cardiovascular disease 20 . FDX1 serves as a key regulator of cuproptosis, a novel form of cell death induced by copper, and operates as an upstream regulator of protein lipoylation. FDX1 knockout results in a complete loss of protein lipoylation, shielding cells from copper toxicity 12 . Researchers found that FDX1 is responsible for reducing Cu2 + to Cu+, a more toxic form, and rapidly increased under high-level copper circumstances. FDX1 has a higher expression in normal tissues than those in disease conditions, such as Crohn’s disease, endometriosis, glioblastoma, lung fibrosis, cancer from the stomach, colorectum, breast, bile duct, kidney, lung, and thyroid 22 – 26 . However, expression of FDX1 increased in osteosarcoma, osteoarthritis, rheumatoid arthritis, polycystic ovary syndrome, temporal lobe epilepsy, liver cancer, thymoma, large-B-cell lymphoma, prostate cancer, pancreatic cancer, esophageal cancer and so on 26 – 31 . In gastric cancer, acetylated METTL16 upregulates FDX1 mRNA and protein levels via m6A-modification on FDX1 mRNA, which ultimately induces cuproptosis 32 . In glioma, FDX1 may promote its proliferation and migration related to the PI3KAKT/mTOR pathway 33 . In non-small-cell lung cancer, FDX1 is inhibited by METTL3 through copper death-associated pri-miR-21-5p maturation to promote cancer growth and metastasis 34 . Besides, researchers regarded FDX1 as a biomarker to predict the prognosis of kidney renal clear cell carcinoma, hepatocellular carcinoma, gastric cancer, and glioma 35 – 39 . However, the role of FDX1 in PTC development and clinical prognosis remains to be further explored and validated. Thus, our study innovatively investigated and analyzed the association between FDX1 gene expression and distant metastasis of PTC. In this study, we found that the expression of FDX1 in PTC was significantly reduced compared to normal tissues through bioinformatic data mining in TCGA-THCA, ENA_PRJEB11591, GSE33630 , GSE60542 , GSE50760 and GSE191288 . We also validated this finding by performing immunohistochemical analyses and an RT-qPCR analysis of patient tissue samples, revealing a consistent decrease of both FDX1 protein and mRNA expression in tumor tissues. Furthermore, our research shows that FDX1 expression correlates with tumor size, ETE, and tumor TNM classification of PTC patients and significantly reversely correlates with DFS among PTC patients. Consequently, it is well-founded that dysregulation of FDX1 could affect the prognosis of PTC patients. Besides, the statistically significant prediction accuracy of the improvement in 1- and 3-year DFS prognosis indicated that FDX1 could serve as a good prognosis predictor for PTC patients. Besides, by utilizing the LASSO regression model in the TCGA-THCA cohort, we identified six key genes closely associated with cuproptosis in PTC, including HIGD1A, SPCS2, PPIL1, CLCNKB, FHL1, and IQGAP3. Based on these genes, we constructed a prognostic model for PTC. Our model has a good prognosis prediction accuracy surpassing that of clinical features and enhances the prognosis prediction accuracy of the traditional clinicopathological staging systems in PTC patients. Among these genes, HIGD1A has been recognized as a mitochondrial protein that may participate in mitochondrial respiratory chain complex IV biogenesis, particularly under hypoxia, oxidative or metabolic stress 40 , 41 . In tumor-related studies, it has been observed that DNA damage promotes HIGD1A translocation from the mitochondria to the nucleus, ultimately promoting homologous recombination and radio/chemo-resistance in several cancer cells, like colorectal cancer, pancreatic cancer, glioma, hepatocellular, etc 42 – 46 . Additionally, SPCS has emerged as a new and novel classification system based on senescence axis regulators, revealing tumor microenvironment heterogeneity. Specifically, SPCS2 subtype indicates activated oncogenic signaling pathways and metabolic signatures that promote cancer expansion in tumor-related studies, such as clear cell renal carcinoma, and hepatocellular carcinoma 47 , 48 . PPIL1 encodes a protein that shares 46.0% identity in amino acid sequence identity with human cyclophilin A, involved in cell cycle progression and cell survival 49 , 50 . Multiple studies have demonstrated that PPIL1 promotes the growth of cancer cells, like colon cancer, breast cancer, melanoma, and so on 51 , 52 . CLCNKB, encoding CLC-Kb, a member of the CLC family of Cl- channels/transporters, was one of the candidate biomarkers in renal cell carcinoma and esophageal cancer 53 – 55 . FHL1 is identified as a tumor suppressor protein, which acts to inhibit tumor cell growth and migration and has been found markedly down-regulated in lung cancer, papillary thyroid cancer, and so on 56 – 58 . IQGAP3, belonging to the IQ motif containing the GTPase Activating Protein (IQGAP) family, plays a crucial role in cell growth and proliferation, cell migration, as well as tumor invasiveness 59 – 61 . What’s more, IQGAP3 has been linked to the poor prognosis of several cancers, leading researchers to describe it as a new diagnostic or detection index for colorectal cancer, breast cancer, and lung cancer 62 – 64 . Taken together, all of the six crucial genes regulate metabolism-related pathways and are involved in the development and progression of tumors. Though researchers found CLCNKB and FHL1 were downregulated in PTC, no evidence shows the association between FDX1/cuproptosis and these two genes in PTC 65 . The available researches show that the other four have not been reported as genes linked to papillary thyroid cancer. The results of the GO analysis showed most FDX1-related genes were significantly enriched in the mitochondrial respiratory-related biological process. Furthermore, PPI network analysis showed that cytochrome c, somatic (CYCS) occupies a key role in the regulatory network. CYCS is viewed as an electron carrier protein participating in the mitochondrial electron transport chain. Additionally, CYCS also plays a role in apoptosis. When mitochondrial membrane permeability is altered either by the suppression of anti-apoptosis members or the activation of pro-apoptosis members of the Bcl-2 family. Subsequently, it binds to Apaf-1 and triggers the activation of caspase-9, which then accelerates apoptosis by activating other caspases 66 . Our studies in clinical samples and TPC-1 and K1 cells show that CYCS expression is linked to FDX1, resulting in the proliferation and apoptosis of papillary thyroid cancer cells. Accordingly, we proposedFDX1 may be involved in PTC pathogenesis by influencing apoptosis through CYCS. However, the precise interaction between CYCS and FDX1 remains undefined and warrants further investigation.

In conclusion, we demonstrated that FDX1 as a key regulator of cuproptosis was significantly downregulated in PTC and has a significant correlation with the tumor recurrence of PTC patients. A novel prognosis indicator, CRRS, was created with six cuproptosis-related genes and displayed good prognosis prediction accuracy. However, the patient volume included in this study from public databases is still limited, and the findings remain to be further verified with additional cohorts. The detailed mechanism of cuproptosis in promoting cancer development and progression remains to be further elucidated. In summary, our findings may help improve prognosis stratification and guide clinical treatment management for PTC patients and may provide new clues for cuproptosis-targeted therapies.

Thyroid cancer (TC) currently ranks as the ninth most prevalent malignancy and tops the list of endocrine cancers, and its incidence has seen a notable surge over recent decades 1 , 2 . Papillary thyroid carcinoma (PTC), the most common thyroid neoplasm, carries the most favorable prognosis of its kind 3 . However, distant metastasis (DM) and regional recurrence, which affect approximately 20% of PTC patients, contribute to poor prognosis while DM is widely regarded as the leading cause of TC-related fatalities 4 , 5 . Despite advancements in radiation therapy, 131 I therapy, endocrine therapy, and targeted therapy, the PTC recurrence rates remain high. Consequently, further exploration of the mechanism of PTC could potentially enhance therapy outcomes and prediction of recurrence for PTC patients. The Nomenclature Committee of Cell Death has established guidelines that categorize cell death modes into accidental cell death (ACD) and regulated cell death (RCD), based on morphological, biochemical, and functional criteria 6 . ACD refers to a biologically uncontrolled process of cell death triggered by accidental injury stimuli 7 . Conversely, RCD is characterized by controlled signaling pathways that play key roles in organismal development or tissue renewal 8 . Various lethal subroutines during RCD, including apoptosis, pyroptosis, necrosis, and ferroptosis, differently influence tumor progression and therapeutic responses 9 . As a result, RCD pathways, exemplified by apoptosis, have gained increasing significance as targets for cancer medication development in recent years 10 . Nonetheless, since tumor cells exhibit evasion mechanisms against RCD, leading to treatment resistance and recurrence, a deeper comprehension of RCD mechanisms may pave the way for establishing diverse potential RCD-related therapeutic strategies. Recently, researchers have discovered a copper-dependent form of cell death distinct from other known types of RCD, termed Cuproptosis 11 . Cuproptosis takes place through the direct binding of copper to lipoylated components of the tricarboxylic acid (TCA) cycle. This process leads to the aggregation of lipoylated protein and the subsequent loss of iron-sulfur cluster proteins, ultimately resulting in proteotoxic stress and cell death 12 . Cuproptosis is modulated by mitochondrial ferredoxin 1 (FDX1)-mediated protein lipoylation. FDX1 serves as an upstream regulator of protein lipoylation and can promote the aggregation of toxic lipoylated DLAT-one of four lipoylated mitochondrial enzymes that accumulate during mitochondrial stress-which is the key mechanism of Cuproptosis. Accumulating evidence indicates that FDX1 expression is significantly lower in most cancer tissues compared to normal adjacent tissue. Additionally, FDX1 expression positively correlates with immune infiltrating cells in most tumors and plays a crucial role in tumor initiation and progression 13 . However, the specific role of cuproptosis-related genes (CRGs) in PTC development and clinical prognosis remains to be further explored and validated. In this study, we used data from the Gene Expression Omnibus (GEO), the Cancer Genome Atlas (TCGA), and the European Nucleotide Archive (ENA) databases to investigate the expression, prognostic significance, and molecular function of FDX1 in PTC. The molecular function of FDX1 was further verified in vitro experiments. Moreover, a cuproptosis-related prognosis model was created with six cuproptosis-related genes (HIGD1A, SPCS2, PPIL1, CLCNKB, FHL1, IQGAP3), a novel prognosis indicator termed CRRS was created with these risk features and displayed a good performance in predicting prognosis. In summary, our findings offer new and novel insights into prognosis and potential cuproptosis-targeted therapies for PTC patients.

Below is the link to the electronic supplementary material. Supplementary Material 1 Supplementary Material 1