Identification and molecular typing of disulfidptosis-related biomarkers in anaplastic thyroid carcinoma

Article type: Research Article

Keywords: Cancer, Cell death

Authors: Weidong Teng ⚬, Yawen Guo, Lingling Ding, Aoni Zhou, Yehao Guo, Jiantong He, Lirong zhang, Haolan Yu, Zekai Tao, Jiafeng Wang, Jiajie Xu, Zhuo Tan, Liehao Jiang

Affiliations: https://ror.org/05gpas306grid.506977.a0000 0004 1757 7957Otolaryngology & Head and Neck Center, Cancer Center, Department of Head and Neck Surgery, Zhejiang Provincial People’s Hospital (Affiliated People’s Hospital), Hangzhou Medical College, Hangzhou, China; https://ror.org/014v1mr15grid.410595.c0000 0001 2230 9154Hangzhou Normal University, Hangzhou, China; Zhejiang Provincial Clinical Research Center for Head & Neck Cancer, Hangzhou, China; Zhejiang Key Laboratory of Precision Medicine Research on Head & Neck Cancer, Hangzhou, China; https://ror.org/05gpas306grid.506977.a0000 0004 1757 7957Hangzhou Medical College, Hangzhou, China; https://ror.org/00rd5t069grid.268099.c0000 0001 0348 3990Wenzhou Medical University, Wenzhou, China; https://ror.org/04epb4p87grid.268505.c0000 0000 8744 8924Zhejiang Chinese Medical University, Hangzhou, China

License: © The Author(s) 2026 CC BY 4.0 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Article links: DOI: 10.1038/s41420-026-03089-9  |  PubMed: 42020364  |  PMC: PMC13237300

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (8.0 MB)

Introduction

Anaplastic thyroid carcinoma (ATC) is a highly aggressive malignancy originating from the epithelial cells of thyroid follicles. It is characterized by extensive local invasion, necrosis, and hemorrhage [ref. 1]. Clinically, ATC predominantly affects elderly women over 65 years old [ref. 2], and typically presents as a rapidly enlarging neck mass with associated symptoms such as neck pain, dysphagia, and airway compression. Approximately 50% of ATC patients succumb to airway obstruction or asphyxia, and nearly 40% present with lymph node or distant metastases at diagnosis [ref. 1]. Despite the availability of surgical intervention, radiotherapy, chemotherapy, and targeted therapy, only a minority of patients qualify for surgery, and radiotherapy provides limited efficacy. Single-modality treatment is often insufficient to control disease progression or improve prognosis, leading to poor survival outcomes [ref. 3–ref. 5]. Thus, there is an urgent need to identify novel biomarkers and elucidate potential mechanisms to enhance our understanding of ATC pathogenesis and develop targeted therapeutic strategies.

Disulfidptosis is a newly characterized mode of cell death mediated by SLC7A11-induced disulfide stress, triggered by excessive intracellular cystine accumulation. In cells with high SLC7A11 expression under glucose-deficient conditions, excessive cystine accumulation leads to the formation of aberrant disulfide bonds between actin cytoskeletal proteins, resulting in cytoskeletal contraction, detachment from the cell membrane, and eventual cell death. Current research on disulfidptosis primarily focuses on its prognostic implications and therapeutic relevance in malignancies such as lung adenocarcinoma and colorectal cancer [ref. 6, ref. 7]. However, the involvement of disulfidptosis-related genes in ATC remains unclear.

This study aims to identify and characterize disulfidptosis-related biomarkers in ATC through bioinformatics analysis, thereby uncovering potential therapeutic targets involved in ATC-associated disulfidptosis.

Results

Common DEGs were obtained by intersection genes of DEGs

In the GSE9115 dataset, 1318 DEGs were detected, including 630 upregulated and 688 downregulated genes (Fig. 1A and Supplementary Fig. 1A). Similarly, in the GSE65144 dataset, a total of 2609 DEGs were identified between ATC and normal thyroid tissue samples, comprising 1370 upregulated and 1239 downregulated genes (Fig. 1B and Supplementary Fig. 1B). By intersecting the DEGs from both datasets, 113 common DEGs were identified (Table S2).

Disulfidptosis and ATC-related genes were obtained through WGCNA

GSVA analysis using DRGs as the reference gene set demonstrated a significant difference in disulfidptosis scores between ATC and normal tissue samples (Fig. 1C). To identify genes associated with both disulfidptosis and ATC, WGCNA was performed. Clustering analysis of the samples confirmed the absence of outliers (Supplementary Fig. 1C). The optimal soft threshold was determined as 6 (R² = 0.861) for network topology construction (Fig. 1D). Gene modules were generated using the dynamic tree-cutting method with a clipping height of 0.4 (Supplementary Fig. 1D, E). Among these, the brown module exhibited the strongest correlation with disulfidptosis scores and ATC (|cor | > 0.6, P < 0.05) and was designated as the key module. A total of 2990 genes from this module were identified as disulfidptosis- and ATC-related genes for further analysis (Fig. 1E).

Candidate genes were identified by intersection genes of common DEGs and disulfidptosis- and ATC-related genes

A total of 32 candidate genes were identified by intersecting common DEGs with disulfidptosis- and ATC-related genes (Supplementary Fig. 2A). Chromosomal localization analysis indicated that these genes were primarily located on chromosomes 1, 3, 7, 10, 20, and 22 (Supplementary Fig. 2B). The PPI network analysis revealed a network comprising 20 nodes and 68 edges (Supplementary Fig. 2C), with ZWINT, CENPA, and PRC1 emerging as core target genes.

Functional enrichment analysis was conducted on 32 candidate genes to determine their potential biological roles, identifying a total of 422 GO terms and 86 KEGG pathways. GO analysis revealed significant enrichment in biological processes such as nucleoside monophosphate biosynthesis, mitotic nuclear division, and microtubule cytoskeleton organization during mitosis (Supplementary Fig. 2D). KEGG pathway analysis indicated that candidate genes were primarily involved in nucleotide metabolism, axon guidance, and the p53 signaling pathway (Supplementary Fig. 2E).

Biomarkers were obtained by machine learning algorithms

To identify signature genes, candidate genes were analyzed using multiple machine learning algorithms. The LASSO regression model results are presented in Fig. 2A. The SVM-RFE algorithm identified 17 signature genes with high classification accuracy, including LGALS1, ADAM12, PFDN2, GNAI3, ELMO1, COL1A2, FOXM1, ATP1B3, MAP4, RAI14, UCK2, COL5A1, CITED2, GMPS, HPSE, TFF3, and PLAUR (Fig. 2B). The Boruta algorithm identified 31 signature genes (Fig. 2C). By intersecting the results of these machine learning approaches, five biomarkers—ATP1B3, TFF3, LGALS1, ADAM12, and COL1A2—were selected (Fig. 2D). To validate the robustness of feature selection, we performed 100 bootstrap LASSO resamplings on the candidate genes. The results showed that eight genes exceeded the 30% selection frequency threshold. Importantly, all five of our identified core biomarkers—LGALS1 (98%), COL1A2 (72%), ADAM12 (58%), TFF3 (46%), and ATP1B3 (39%)—ranked among the top five high-frequency genes, with selection frequencies significantly higher than those of other candidate genes (e.g., RAI14: 37%) (Fig. 2E). These findings fully demonstrate the importance and stability of the five core biomarkers in the model.

The expression analysis and clinical sample validation of biomarkers

The expression patterns of these biomarkers were analyzed in the GSE65144 and GSE9115 datasets, demonstrating upregulation of ATP1B3, LGALS1, ADAM12, and COL1A2 in ATC, whereas TFF3 was downregulated (Fig. 3A, B). Validation in the GSE27155 and GSE53072 datasets confirmed that biomarker expression trends were consistent with findings from the training datasets, with significant differences observed between groups (Fig. 3C, D). Further analysis using the UALCAN database indicated significant expression differences in other thyroid cancers (Fig. 3E).

After completing these experiments, we analysed the available clinical samples (The data of proteomics have been uploaded in the ProteomeXchange Consortium as reported previously (iProX database, project ID: IPX0008941000))(including 99 ATC paracancerous tissues and 115 ATC tissue samples), and analysed three genes (ATP1B3, COL1A2 and LGALS1) with certain significance. Therefore, we carried out the next step of expression measurement for the three genes (Fig. 3F), in which it is obvious that ATP1B3 and LGALS1 were significantly elevated in ATC compared to normal paracancerous thyroid tissues (NTHY), and the trend of ATP1B3 was more obvious, but the trend of COL1A2 was indeed a downward trend. The heatmap also reflects a view of the appeal (Fig. 3G).

Nomogram model had high prediction accuracy for ATC

A nomogram model incorporating ATP1B3, TFF3, LGALS1, ADAM12, and COL1A2 was constructed to predict ATC occurrence (Fig. 4A). Calibration curves indicated excellent predictive accuracy of the model (Fig. 4B). The ROC curve analysis revealed an AUC of 1.0 for the nomogram model, an AUC of 1.0 for four biomarkers (TFF3, LGALS1, ADAM12, and COL1A2), and an AUC of 0.987 for ATP1B3 (Fig. 4C). To rigorously address potential concerns of overfitting and to evaluate the model’s generalizability, we performed an internal validation using the bootstrap method with 1000 resamples on both the GSE65144 and GSE9115 datasets. The results demonstrated that in the GSE9115 dataset, the bias-corrected AUCs for TFF3, LGALS1, and ATP1B3 were all 1.000 (95% CI: 1.000–1.000), with effective sampling rates all exceeding 0.99, confirming that their perfect classification performance was not attributable to small-sample artifacts. The corrected AUCs for ADAM12 and COL1A2 also remained stable in the range of 0.947–0.949 (Fig. 4D). In the GSE65144 dataset, LGALS1 similarly maintained an AUC of 1.000, while the corrected AUCs for the other biomarkers (ranging from 0.921 to 0.962) were highly consistent with the original values (Fig. 4E). Decision curve analysis demonstrated a strong likelihood of clinical benefit for both individual biomarkers and the nomogram model, with overlapping curves for TFF3, LGALS1, ADAM12, and the overall nomogram model (Fig. 4F). The nomogram model was validated in the GSE9115 dataset (Fig. 4G–I), and the corrected C-index was 0.95, indicating strong predictive ability of the nomogram. The ROC curves of each biomarker and the nomogram model all exceeded 0.8, suggesting high diagnostic accuracy. The decision curve analysis showed that both single genes and the overall model offered substantial clinical benefits, with the curves of COL1A2, ADAM12, and the Nomogram Model overlapping. Furthermore, the external validation conducted using the GSE29265 dataset demonstrated that the AUC value derived from the ROC curve analysis indicated its excellent discriminatory ability (Fig. 4J).

Functional enrichment of candidate genes and biomarkers

To further explore the functional relevance of the biomarkers, GSEA was performed. Single-gene GO enrichment analysis demonstrated that most of the five biomarkers were significantly associated with biological processes such as collagen fibril organization, chromosome segregation, and sister chromatid segregation. KEGG analysis identified the most significantly enriched pathways as ECM-receptor interaction, cell cycle regulation, and focal adhesion. The top five GO terms and KEGG pathways are displayed in Supplementary Fig. 3A–J.

Biomarkers were associated with infiltrating cells

Using the CIBERSORT algorithm, the distribution of 22 immune cell types in normal and ATC tissue samples from the GSE65144 dataset was analyzed (Fig. 5A). Five immune cell types—follicular helper T cells, M2 macrophages, plasma cells, activated mast cells, and gamma delta T cells—exhibited significant differences in abundance between ATC and normal tissues (Fig. 5B), suggesting a potential role in ATC progression. Correlation analysis revealed strong associations between differentially abundant immune cells and biomarkers (Fig. 5C). M2 macrophages exhibited positive correlations with ATP1B3, LGALS1, ADAM12, and COL1A2 but were negatively correlated with TFF3. Conversely, follicular helper T cells demonstrated a negative correlation with ATP1B3, LGALS1, ADAM12, and COL1A2 while being positively correlated with TFF3.

Two subtypes were identified through consensus clustering

Consensus clustering analysis indicated that the CDF plot stabilized when the 12 ATC samples from the GSE65144 dataset were categorized into two distinct subtypes (Fig. 6A). PCA confirmed that the two subtypes (C1 and C2) were well separated (Fig. 6B). Differential expression analysis of biomarkers between the subtypes revealed significant differences in LGALS1, TFF3, and COL1A2 expression, whereas ATP1B3 and ADAM12 showed no significant variation (Fig. 6C).

A total of 626 DEGs were identified between the two subtypes, consisting of 365 upregulated and 261 downregulated genes. GO and KEGG enrichment analyses demonstrated that these DEGs were predominantly involved in neutrophil migration, granulocyte chemotaxis, and myeloid leukocyte migration (Fig. 6D), as well as classical signaling pathways such as NF-kappa B and PI3K-Akt (Fig. 6E). Analysis of 28 immune cell types in the C1 and C2 subtypes of the GSE65144 dataset revealed significant differences in four immune cell types: activated dendritic cells, follicular helper T cells, regulatory T cells, and myeloid-derived suppressor cells (Fig. 6F). The Wilcoxon test indicated statistically significant differences in these immune cell populations between subtypes C1 and C2 (Fig. 6G). Correlation analysis further confirmed that ATP1B3, LGALS1, ADAM12, and COL1A2 were positively associated with activated dendritic cells and regulatory T cells, whereas TFF3 was negatively correlated with these immune cells (Fig. 6H).

ATP1B3 knockdown may inhibit malignant behaviour of ATC cells by triggering disulfidptosis

As a key gene in the diagnostic model, ATP1B3 was selected for further validation through cellular and animal experiments. Two ATC cell lines, 8505C and CAL62, were used for these studies. Stable ATP1B3 knockdown was achieved via viral transfection, and knockdown efficiency was confirmed by Western blot (WB) and PCR (Fig. 7A, B). Sequences 2 and 3 demonstrated higher knockdown efficiency and were therefore selected for subsequent preliminary functional experiments.

Functional assays revealed that ATP1B3 knockdown significantly suppressed proliferation and colony formation (Fig. 7C, D). Sequence 3 demonstrated a more pronounced inhibitory effect than Sequence 2; consequently, all subsequent experiments employed Sequence 3 for further investigation. Moreover, migration and invasion capabilities were markedly reduced in the ATP1B3 knockdown group compared to the control group (Fig. 7E), suggesting that ATP1B3 plays a critical role in maintaining ATC cell viability. Given that disulfidptosis is often associated with elevated expression of SLC7A11 (XCT), SLC7A11 expression was examined in both cell lines. WB analysis confirmed an increase in SLC7A11 expression in ATP1B3 knockdown cells compared to controls (Fig. 7F). Flow cytometry demonstrated a significant increase in cell death following ATP1B3 knockdown (Fig. 7G). Furthermore, NADP + /NADPH assay results revealed significantly increased NADPH consumption in the knockdown group. For the 8505 C cell line, the average NADPH concentration in the control group was 39.640 pmol/250,000 cells, whereas the knockdown group exhibited an average concentration of 14.763 pmol/250,000 cells. In the CAL62 cell line, the average NADPH concentration in the NC group was 78.703 pmol/250,000 cells, while the average NADPH concentration in the knockdown group was 20.108 pmol/250,000 cells (Fig. 7H), consistent with previous research findings. Subsequently, we employed confocal microscopy to observe that F-actin exhibited marked contraction in ATP1B3-knockdown cells, further corroborating its role in disulfide-bond apoptosis (Fig. 7I). Finally, using a cysteine assay kit, we observed a marked reduction in cellular cysteine levels following ATP1B3 gene knockdown compared to the control group. The decrease in NADPH content prevented the reduction of cystine to cysteine within the cells (Fig. 8A).

Rescue trial

Simultaneously, we procured overexpression plasmids for the target gene. Cells exhibiting favourable post-transfection morphology were selected for plating. Following complete cell confluence, further transfection experiments were conducted using Lipofectamine 3000 as the primary transfection reagent. Cultures were maintained in 1640 medium supplemented with double FBS. Upon reaching the requisite cell proliferation levels, Western blot (Fig. 8B) assays were promptly performed to validate overexpression efficiency. Having confirmed overexpression efficiency, we proceeded with selected phenotypic assays: CCK8 proliferation assays (5000 cells per well) (Fig. 8C) and colony formation assays (10,000 cells per well, fixed and stained immediately upon reaching the required proliferation level) (Fig. 8D). These tests further corroborated the impact of ATP1B3 knockdown on cellular growth and proliferation.

Consistent in vivo tumor growth after ATP1B3 knockdown

To confirm the in vitro findings, a xenograft mouse model was established to evaluate the effect of ATP1B3 knockdown on tumor growth. Tumor growth was significantly suppressed in the ATP1B3 knockdown group compared to the control group (Fig. 8E, F). Consistent with tumor volume measurements, ATP1B3 knockdown also resulted in a significant reduction in tumor weight (Fig. 8G), highlighting its strong anti-tumorigenic effect. Immunohistochemistry analysis further confirmed decreased expression levels of HE, Ki67, and ATP1B3 in ATP1B3 knockdown tumors, reinforcing its potential role in ATC suppression. The reduction of F-actin also confirms the existence of a greater or lesser relationship with disulfide death in a sideways manner (Fig. 8H). ATP1B3 expression was significantly reduced in the ATP1B3 knockout group (P < 0.01), as was Ki67 expression (P < 0.05), further confirming the inhibitory effect of ATP1B3 knockout on tumour growth (Fig. 8H).

Discussion

ATC is a rare but highly aggressive malignancy, accounting for only 1–2% of all thyroid cancers but contributing to the majority of thyroid cancer-related deaths. The five-year survival rate is a mere 7%, with a median survival time of approximately five months. Despite extensive research efforts, no significant advancements have been made in improving the overall survival of ATC over the past three decades.

Disulfidptosis, a novel mode of cell death first described in 2023 by Prof. Gan Poyi et al., is triggered by SLC7A11-mediated disulfide stress due to intracellular cystine accumulation. This form of cell death is characterized by excessive disulfide bond formation between actin cytoskeletal proteins under glucose deprivation in cells with high SLC7A11 expression. These aberrant disulfide bonds cause actin filament contraction, detachment from the plasma membrane, cytoskeletal collapse, and ultimately cell death. By intersecting differentially expressed genes with those related to both disulfidptosis and ATC, five key genes were identified and used to construct a molecular diagnostic model that demonstrated high predictive accuracy for ATC [ref. 6].

Among the identified biomarkers, ATP1B3 was previously reported to be significantly upregulated in hepatocellular carcinoma (HCC) and associated with immune cell infiltration and immune-related cytokine expression [ref. 8]. In estrogen receptor (ER)-positive breast cancer, ATP1B3 overexpression has also been linked to increased cisplatin sensitivity, suggesting its value as both a prognostic biomarker and a potential therapeutic target in HCC and ER-positive breast cancer [ref. 9, ref. 10]. In this study, ATP1B3 knockdown in ATC cells significantly inhibited proliferation, and its high expression correlated with poorer prognosis, highlighting its potential as a biomarker for prognosis assessment and therapeutic intervention in ATC.

LGALS1 (Galectin-1), a β-galactoside-binding lectin, is overexpressed in non-small cell lung cancer (NSCLC) and promotes NSCLC cell proliferation, migration, and invasion through the NCAPG/LGALS1/SPARC axis [ref. 11, ref. 12]. High expression of LGALS1 has also been associated with adverse outcomes in ovarian and oral cancers [ref. 13, ref. 14]. Although its exact role in ATC remains undefined, its known pro-metastatic and immune-regulatory functions suggest that LGALS1 may enhance ATC cell invasiveness and facilitate immune evasion, thereby contributing to tumor progression. Further investigation is warranted to elucidate its precise mechanisms in ATC and assess its potential as a therapeutic target.

TFF3 (Trefoil factor 3), a member of the trefoil factor family, is upregulated in various malignancies, including breast, gastric, colorectal, thyroid, and cervical cancers. Its oncogenic effects are thought to be mediated through the MAPK/ERK pathway, promoting tumor proliferation and metastasis [ref. 15]. Notably, TFF3 polymorphisms have been linked to susceptibility to papillary thyroid carcinoma [ref. 16]. Consistent with our findings, Li et al. also identified TFF3 as a diagnostic biomarker for ATC with robust diagnostic performance [ref. 17]. These data suggest that TFF3 may play a critical role in ATC pathogenesis and serve as a promising diagnostic and therapeutic target.

COL1A2, which encodes type I collagen, has been implicated in tumor differentiation. Its inhibition suppresses glioma cell proliferation and invasion [ref. 18]. Interestingly, while COL1A2 is upregulated in gastric and epithelial ovarian cancers [ref. 19, ref. 20], it exhibits tumor-suppressive functions in bladder cancer [ref. 21] and head and neck cancers [ref. 22]. This dual role suggests a context-dependent function, and further investigation is needed to clarify its mechanistic role in ATC and assess its therapeutic potential.

ADAM12, a member of the ADAM (a disintegrin and metalloproteinase) family, contributes to tumor development and progression in renal clear cell carcinoma [ref. 23], breast cancer [ref. 24], and hepatocellular carcinoma [ref. 25]. It promotes renal carcinoma progression via the EGFR/ERK signaling pathway [ref. 23] and has been proposed as a diagnostic biomarker for hepatocellular carcinoma [ref. 25]. In breast cancer, ADAM12 (that is regulated by TWIST1) enhances invasion and metastasis through pseudopodia and focal adhesion formation [ref. 26]. These findings support the hypothesis that ADAM12 may play a similar role in ATC, particularly in promoting metastasis.

To further investigate the functions of these biomarkers, enrichment analysis indicated significant associations with biological processes such as nucleoside monophosphate biosynthesis during mitosis, mitotic nuclear division, and microtubule cytoskeleton organization. KEGG pathway analysis revealed enrichment in nucleotide metabolism, axon guidance, and the p53 signaling pathway, consistent with previous studies and offering insights for future exploration. To experimentally validate these findings, ATP1B3 was selected for in vitro functional studies. Western blot and PCR analyses confirmed that ATP1B3 knockdown significantly inhibited ATC cell proliferation, colony formation, migration, and invasion. Moreover, Western blot showed increased XCT expression in ATP1B3-knockdown cells. Flow cytometry using PI staining revealed a significantly elevated apoptosis rate in the knockdown group. Confocal microscopy identified morphological alterations, including reduced cell size and nuclear wrinkling, indicating that ATP1B3 may be mechanistically linked to disulfidptosis in ATC cells.

While this study is the first to identify disulfidptosis-associated biomarkers in ATC, several limitations should be acknowledged. First, the marked heterogeneity of ATC may contribute to variability in gene expression profiles, and the relatively small sample size may limit statistical power and generalizability. Second, this study focused on a single ATC subtype, and its relevance to other thyroid cancer subtypes remains to be validated. Additionally, the absence of survival and clinical outcome data limits the ability to fully assess the translational potential of these biomarkers.To overcome these limitations, we shall expand our sample size and incorporate additional subtypes in future studies. By integrating multi-omics analyses with cellular and animal model experiments, we aim to systematically validate the functional mechanisms of these biomarkers. Specifically, we plan to conduct further rescue experiments in ATP1B3 knockdown models to clarify its causal role in ATC progression and disulfidptosis. Furthermore, quantitative assessment of total disulfide bond levels via DTNB assay will furnish the most direct biochemical evidence for ATP1B3’s regulation of disulfidptosis, ultimately elucidating its molecular mechanisms in ATC. Concurrently, extending observation periods in animal models and conducting long-term clinical follow-up will facilitate evaluation of biomarker-guided therapy efficacy and safety, thereby establishing robust preclinical and translational foundations for personalised ATC treatment.

Materials and methods

Data source

The four ATC-related datasets (GSE65144, GSE9115, GSE27155, and GSE53072) were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). The GSE65144 dataset contained 13 normal and 12 ATC thyroid tissue samples, while GSE9115 included 4 normal and 5 ATC thyroid tissue samples. Similarly, GSE27155 comprised 4 normal and 4 ATC thyroid tissue samples, and GSE53072 contained 4 normal and 5 ATC thyroid tissue samples. In addition, the GSE29265 dataset was also included for external validation, which contained 20 normal thyroid tissue samples and 29 ATC tissue samples. Additionally, 16 disulfidptosis-related genes (DRGs; Table S1) were retrieved from previously published literature [ref. 27].

Differential expression analysis

Differential expression analysis was conducted to separately identify differentially expressed genes (DEGs) between ATC and normal samples in GSE65144 and GSE9115 using the ‘limma’ R package (version 3.52.4) [ref. 28], with the thresholds set at |Log₂FC | > 1 and adj. P value < 0.05. Common DEGs were then determined using the Robust Rank Aggregation (RRA) method. The expression patterns of common DEGs were visualized using volcano and heat maps generated with the ggplot2 R package (version 3.3.2) [ref. 29] and the Pheatmap R package (version 1.0.12) [ref. 30].

Identification of disulfidptosis- and ATC-related genes using Weighted Gene Co-expression Network Analysis (WGCNA)

Single-sample Gene Set Enrichment Analysis (GSVA) was performed to calculate the disulfidptosis score for each sample based on DRGs using the GSVA R package [ref. 28]. Co-expression networks were constructed and analyzed in the GSE65144 and GSE9115 datasets using the WGCNA R package. First, clustering analysis was conducted to identify potential outliers. A suitable soft threshold was selected based on the nearly scale-free topology criterion. The topological matrix was then clustered using differences between genes, and neighborhood and dissimilarity coefficients were calculated. Modules were identified using a dynamic tree-cutting method (minimum module size: 100 genes). Correlation analysis between modules and phenotypes (disulfidptosis score and ATC) was conducted, and modules with |cor | > 0.6 and P < 0.05 were designated as key modules. Genes within these modules were identified as disulfidptosis- and ATC-related genes.

Identification of candidate genes

Common DEGs overlapping with disulfidptosis- and ATC-related genes were considered candidate genes. Chromosomal localization of these candidate genes was visualized using the RCircos R package (version 1.2.2) [ref. 31]. A protein-protein interaction (PPI) network was constructed using the STRING database (https://string-db.org/). Additionally, correlation analysis of candidate genes was performed using the ‘igraph’ R package (version 14.1), with a threshold of |R | > 0.7.

Machine learning algorithms

Biomarkers were identified using three machine learning algorithms: least absolute shrinkage and selection operator (LASSO), support vector machine recursive feature elimination (SVM-RFE), and Boruta. Candidate genes were analyzed using the glmnet (Version 4.1-4) [ref. 32], e1071 (Version 1.7-13), and Boruta (Version 8.0.0) [ref. 33] packages, respectively. The feature genes identified by these algorithms were overlapped using the VennDiagram R package (Version 1.7.3) [ref. 34] to determine biomarkers.

The expression levels of biomarkers were evaluated in the training sets (GSE65144 and GSE9115) and validation sets (GSE27155 and GSE53072). Additionally, the expression levels of biomarkers in other thyroid cancers (THCA) were analyzed using the UALCAN database (http://ualcan.path.uab.edu). Finally, the prognostic value of biomarkers in THCA was assessed using survival analysis in the GEPIA database (http://gepia.cancer-pku.cn/). To evaluate the stability of the identified biomarkers, we performed bootstrap LASSO regression with 100 resamples on the candidate genes. The frequency of each gene being selected across all bootstrap samples was calculated.

Construction of a nomogram model

To predict ATC occurrence, a nomogram model was constructed based on the identified biomarkers using the rms R package (Version 6.6-0) [ref. 32]. Calibration curves were plotted to evaluate the predictive accuracy of the nomogram model. Receiver Operating Characteristic (ROC) curves were generated for both the biomarkers and the nomogram model to assess diagnostic performance. Additionally, Decision Curve Analysis (DCA) was performed to determine the clinical utility of the nomogram model. Subsequently, the nomogram model was validated in the GSE9115 dataset. To assess the robustness and potential overfitting of the nomogram model, we performed internal validation using the Bootstrap method with 1000 resamples. The bias-corrected 95% confidence intervals (CIs) for the area under the receiver operating characteristic curve (AUC) were calculated to evaluate the model’s discriminative ability. To further validate the robustness and generalizability of the nomogram model, external validation was performed using the GSE29265 dataset.

Enrichment analysis

To explore the biological functions of candidate genes and biomarkers, enrichment analysis was conducted. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the clusterProfiler R package (version 4.4.4) [ref. 35]. The adjusted P value < 0.05 was set as the cut-off criteri [ref. 36, ref. 37]. Furthermore, Gene Set Enrichment Analysis (GSEA) was conducted to investigate the functional roles of biomarkers, with thresholds set at |NES | > 1 and P < 0.05.

Immune infiltration analysis

Immune cell infiltration analysis was performed using the CIBERSORT algorithm [ref. 38] on ATC and normal tissue gene expression data from the GSE65144 dataset. Significantly different immune cell populations between the two groups were identified using the Wilcoxon test. To investigate the relationship between differential immune cells and biomarkers, Spearman correlation analysis was conducted between significantly different immune cells and biomarker expression levels.

Consensus clustering analysis

Consensus clustering is a widely applied unsupervised method for subtype classification. In this study, unsupervised clustering was performed on 12 ATC samples from the GSE65144 dataset based on the identified biomarkers using the ConsensusClusterPlus R package (Version 1.60.0) [ref. 39]. The optimal number of subtypes was determined by evaluating the cumulative distribution function (CDF) curves. To validate the classification accuracy, principal component analysis (PCA) was conducted to assess sample distribution across different subtypes. Subsequently, differential expression analysis of biomarkers between subtypes was performed, and DEGs were identified using the thresholds of P < 0.05 and |log₂FC | > 1. Functional enrichment analysis of DEGs, including GO and KEGG pathway analysis, was conducted using the ClusterProfiler package. The extent of immune cell infiltration between subtypes in the GSE65144 dataset was quantified using the GSVA algorithm, and significant differences were determined using the Wilcoxon test. Finally, correlations between biomarkers and immune cells were analyzed using Spearman’s correlation test.

Statistical analysis

All statistical analyses were conducted using R software. The Wilcoxon test was applied to compare scores between two groups, while correlation analysis was performed using Spearman’s test. A P value < 0.05 was considered statistically significant.

Clinical sample selection

The procedures used in this study comply with the ethical standards of the Ethics Committee of Zhejiang Provincial People’s Hospital. Thyroid cancer samples and adjacent thyroid tissue were collected from Zhejiang Provincial People’s Hospital, Zhejiang Cancer Hospital, and Shaoyifu Hospital (January 2010 to January 2022). Each patient signed a written informed consent form agreeing to have their samples and clinical information used for research purposes, and the study protocol was approved by the Ethics Committee of Zhejiang Provincial People’s Hospital(Approval No: 202510151505000115513).

Cell culture

ATC cell lines 8505c, and CAL62 were obtained from Procell (Wuhan, China), and the Shanghai Cell Bank of the Chinese Academy of Sciences, respectively. Cells were cultured in RPMI 1640 medium (Hyclone) supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin under standard conditions (37 °C, 5% CO₂). Once a stable state was achieved, cells were transfected with either negative control (NC) or ATP1B3 knockdown lentiviruses (Jikai Genetics, China). After 24 h of transfection, cells were transferred to 10-cm culture dishes for further incubation.

Western blotting

Following transfection, cells were lysed in protein lysis buffer (Cell lysis buffer for Western and IP; Beyotime, China) containing protease inhibitor PMSF. Equal amounts of protein were separated via SDS-PAGE and transferred to a PVDF membrane. Protein expression was visualized using an enhanced chemiluminescence (ECL) detection kit (Fude Biological, China). The primary antibodies used were ATP1B3 (67554-1-Ig, Proteintech Group, USA), β-Tubulin pAb (Ac008, Abcam, USA), SLC7A11 (26846-1-AP, Proteintech Group, USA), and HSP60 (15282-1-AP, Proteintech Group, USA).

Reverse transcription-polymerase chain reaction (RT-PCR)

Following transfection, RNA was extracted from both control and knockdown cell lines using an RNA Rapid Extraction Kit. Reverse transcription was performed using the PrimeScript RT Master Mix (Takara, Japan), and gene expression was analyzed via quantitative PCR (qPCR) using the Hieff^® qPCR SYBR Green Master Mix (Yeasen Biotechnology, China). 18S was used as an internal control, and gene expression changes were quantified using the 2^-ΔΔCt method. Primer sequences used for PCR are listed in Table 1.

CCK8

The CCK-8 assay was employed to evaluate the proliferative capacity of ATC cells. Following validation via PCR and Western blot analysis, the more efficient knockdown sequences 2 and 3 were selected for subsequent experiments. 8505c and CAL62 cells were seeded at a density of 5000 cells per well in 96-well plates. CCK-8 reagent was added every 24 h, and cell viability was assessed by measuring the optical density (OD) at 450 nm. Each experiment was performed in triplicate.

Cloning formation

Cell lines used for colony formation assays were the same as those used in previous experiments. Control and experimental groups were seeded in six-well plates at a density of 1000 cells per well. Cells were incubated until colony formation reached an optimal density. Colonies were fixed with 1 mL of paraformaldehyde for 20 min at room temperature, followed by staining with 500 µL of crystal violet for 15 min. After three washes with PBS, plates were air-dried and imaged.

Transwell migration and invasion assay

Cell migration and invasion were assessed using Transwell plates (LABSELECT, Hefei, Anhui, 6.5 mm). ATC cell lines (8505c, CAL62) were seeded into the upper chamber at a density of 2–4 × 10⁴ cells per well in serum-free medium (200 μL), while 800 μL of medium containing serum was added to the lower chamber. After 24 or 48 h of incubation, cells that had migrated or invaded were fixed with 0.1% crystal violet for 30 min. The number of migrated cells was quantified under an M7000 microscope by selecting representative fields of view. The average number of migrated cells was calculated from five random fields across three independent experiments.

Flow cytometry analysis

Apoptosis analysis was performed using flow cytometry. Control and experimental groups of 8505c and CAL62 cells were seeded in six-well plates at a density of 1 × 10⁵ cells per well. After reaching adherence, cells were detached using trypsin, centrifuged, and washed 2–3 times with PBS. Cells were stained with apoptosis detection kits using single propidium iodide (PI) staining. Samples were analyzed using a flow cytometer (Beckman Colter, Inc., Ireland), and data were processed using FlowJo software.

Confocal laser scanning microscope

For confocal imaging, 8505c (30,000 cells/well) and CAL62 (50,000 cells/well) were seeded in 35 mm glass-bottom dishes with quadruple compartments. After cell adherence, cells were washed three times with PBS (5 min per wash). Cells were then fixed with 4% paraformaldehyde for 1 h, followed by incubation with 200 µL of primary antibody dilution for 40 min. After three PBS washes, 80 µL of F-actin (ab205, 1:1000) was applied overnight at 4 °C in a humidified, dark chamber. Subsequently, cells were washed with PBS and incubated with 200 µL of Alexa Fluor® 488-conjugated AffiniPure Goat Anti-Rabbit IgG H&L (1:1000) for 1 h in the dark. DAPI staining (200 µL) was applied for 5 min. Fluorescence images were captured using a TCS SP8 confocal microscope.

NADP + /NADPH

NADP + /NADPH levels were measured in control and ATP1B3 knockdown groups. One million cell precipitates were collected from both groups for analysis. NADP + /NADPH levels were quantified using the Biyunfei NADP + /NADPH Assay Kit (WST-8 method) according to the manufacturer’s instructions. In accordance with the experimental protocol, a blank control group, a standard sample group, and an experimental group were established. A standard curve was plotted using the standard sample group. The absorbance values measured at 450 nm for the experimental group were substituted into the curve to calculate the corresponding NADPH concentrations. Data analysis and graphing were performed using Prism 10 software.

Cysteine testing

Establish four groups: (NC group, NC group supplemented with glucose-free medium, ATP1B3 knockdown group, knockdown group supplemented with glucose-free medium). Following cell processing the previous day, seed 300,000 cells into each well of a 6-well plate. The following day, after cells have adhered, replace with the corresponding medium. Continue culture for 6–8 h before harvesting the cell pellet. Cysteine was subsequently extracted using the Cysteine Content Assay Kit (boxbio, Cat. No.: AKAM004M). Absorbance was measured at 600 nm using a microplate reader, and a curve was plotted against the generated standard curve to determine cysteine content. Data were analysed and graphically represented using Prism 10 software.

Overexpression rescue assay

Following completion of the aforementioned experiments, cells exhibiting favourable condition were selected for plating at a density of 200,000 cells per well. After allowing the cells to adhere overnight, transfection with the overexpression plasmid (Shanghai Genechem Co., Ltd.) was performed. Once cells reached the desired proliferation level, they were promptly processed for Western blot analysis to validate gene overexpression efficiency. Upon confirmation of overexpression efficiency, baseline phenotypic assays were immediately conducted, including the CCK8 proliferation assay and colony formation assay.

In vivo xenograft tumor model

Female nude mice (BALB/c, 3–5 weeks old) were purchased from Hangzhou Qizhen Experimental Animal Technology (Hangzhou, China) and maintained at the laboratory animal center at Zhejiang Provincial People’s Hospital (Approval No. IACUCA20230309001) according to the requirements of feeding. All mice were randomly assigned to different groups for the test. One week later, CAL62 control and ATP1B3 knockdown cell suspensions were injected subcutaneously into the right abdominal cavity of nude mice. Moreover, the volume (0.5 × (shortest diameter)2 × (longest diameter)) and body weight were recorded every other day. Once tumors reached the required experimental volume, mice were euthanized humanely, and tumors were excised, weighed, and fixed in formalin for further analysis.

Immunohistochemistry

Paraffin-embedded tumor tissues were sectioned into 4-μm thick slices and subjected to deparaffinization and rehydration. Antigen retrieval was performed using citrate buffer (pH 6.0), followed by blocking with 3% bovine serum albumin for 30 min. Sections were incubated with primary antibodies at 25 °C for 1 h, followed by HRP-labeled secondary antibodies at room temperature for 30 min. Nuclei were counterstained with hematoxylin, and positive 3,3’-diaminobenzidine (DAB) staining was visualized as a brownish-yellow coloration. Images were captured using a fluorescence microscope (Konfoong Biotech, Ningbo, China). Concurrently, we performed further H-Score staining on immunohistochemical sections using Visiopharm software. The entire tissue area was selected for H-Score analysis of the target object. H-Score = ∑(pi × i) = (percentage of weakly stained cells × 1) + (percentage of moderately stained cells × 2) + (percentage of strongly stained cells × 3), where i denotes the positivity grade classification: negative (no staining) = 0 points; weakly positive (pale yellow) = 1 point; moderate positivity: brownish-yellow, scored 2 points; strong positivity: dark brown, scored 3 points. pi denotes the percentage of positive cells. The H-score ranges from 0 to 300, with higher values indicating greater overall positive intensity.

References

1. 1.Giuffrida DGH. Anaplastic thyroid carcinoma: current diagnosis and treatment. Ann Oncol. 2000. 10.1023/a:1008322002520.
2. V Bisof, Z Rakusic, M Despot. Treatment of patients with anaplastic thyroid cancer during the last 20 years: whether any progress has been made?. Eur Arch Otorhinolaryngol, 2015. [DOI | PubMed]
3. 3.Xu T, Zhu C, Chen J, Song F, Ren X, Wang S, et al. ISG15 and ISGylation modulates cancer stem cell-like characteristics in promoting tumor growth of anaplastic thyroid carcinoma. J Exp Clin Cancer Res. 2023;42:182. 10.1186/s13046-023-02751-9.
4. V Tiedje, M Stuschke, F Weber, H Dralle, L Moss, D Führer. Anaplastic thyroid carcinoma: review of treatment protocols. Endocr-Relat Cancer, 2018. [DOI | PubMed]
5. ND Perrier, JD Brierley, RM Tuttle. Differentiated and anaplastic thyroid carcinoma: major changes in the American Joint Committee on Cancer eighth edition cancer staging manual. CA: A Cancer J Clin, 2017. [DOI]
6. 6.Yan Y, Teng H, Hang Q, Kondiparthi L, Lei G, Horbath A, et al. SLC7A11 expression level dictates differential responses to oxidative stress in cancer cells. Nat Commun. 2023;14:3673. 10.1038/s41467-023-39401-9.
7. X Liu, L Zhuang, B Gan. Disulfidptosis: disulfide stress-induced cell death. Trends Cell Biol, 2024. [DOI | PubMed]
8. 8.Lu S, Cai S, Peng X, Cheng R, Zhang Y. Integrative transcriptomic, proteomic and functional analysis reveals ATP1B3 as a diagnostic and potential therapeutic target in hepatocellular carcinoma. Front Immunol. 2021. 10.3389/fimmu.2021.636614.
9. 9.Li X-x, Wang D-q, Sui C-g, Meng F-d, Sun S-l, Zheng J, et al. Oleandrin induces apoptosis via activating endoplasmic reticulum stress in breast cancer cells. Biomed Pharmacother. 2020. 10.1016/j.biopha.2020.109852.
10. 10.Acconcia F. Evaluation of the sensitivity of breast cancer cell lines to cardiac glycosides unveils ATP1B3 as a possible biomarker for the personalized treatment of ERα expressing breast cancers. Int J Mol Sci. 2022. 10.3390/ijms231911102.
11. 11.Sun H, Zhang H, Yan Y, Li Y, Che G, Zhou C, et al. NCAPG promotes the oncogenesis and progression of non-small cell lung cancer cells through upregulating LGALS1 expression. Mol Cancer. 2022. 10.1186/s12943-022-01533-9.
12. M Tichet, V Prod’Homme, N Fenouille, D Ambrosetti, A Mallavialle, M Cerezo. Tumour-derived SPARC drives vascular permeability and extravasation through endothelial VCAM1 signalling to promote metastasis. Nat Commun, 2015. [DOI | PubMed]
13. JM Corral, LD Puerto-Nevado, M Cedeno, A Rio-Vilariño, I Mahillo-Fernández, C Galeano. Galectin-1, a novel promising target for outcome prediction and treatment in SCLC. Biomed Pharmacother, 2022. [DOI | PubMed]
14. JM Li, CW Tseng, CC Lin, CH Law, YA Chien, WH Kuo. Upregulation of LGALS1 is associated with oral cancer metastasis. Ther Adv Med Oncol, 2018. [DOI | PubMed]
15. X Lin, H Zhang, J Dai, W Zhang, J Zhang, G Xue. TFF3 Contributes to epithelial-mesenchymal transition (EMT) in papillary thyroid carcinoma cells via the MAPK/ERK signaling pathway. J Cancer, 2018. [DOI | PubMed]
16. DXG Pei, JF Wu, DM Wang, XJ Qin. Clinical analysis between trefoil factor 3 gene polymorphsims and susceptibility to papillary thyroid carcinoma. Lin Chuang Er Bi Yan Hou Tou Jing Wai Ke Za Zhi, 2019. [DOI | PubMed]
17. C Li, X Dong, Q Yuan, G Xu, Z Di, Y Yang. Identification of novel characteristic biomarkers and immune infiltration profile for the anaplastic thyroid cancer via machine learning algorithms. J Endocrinol Invest, 2023. [DOI | PubMed]
18. Y Wang, M Sakaguchi, H Sabit, S Tamai, T Ichinose, S Tanaka. COL1A2 inhibition suppresses glioblastoma cell proliferation and invasion. J Neurosurg, 2023. [DOI | PubMed]
19. J Li, Y Ding, A Li. Identification of COL1A1 and COL1A2 as candidate prognostic factors in gastric cancer. World J Surg Oncol, 2016. [DOI | PubMed]
20. S Li, H Li, Y Xu, X Lv. Identification of candidate biomarkers for epithelial ovarian cancer metastasis using microarray data. Oncol Lett, 2017. [DOI | PubMed]
21. K Mori, H Enokida, I Kagara, K Kawakami, T Chiyomaru, S Tatarano. CpG hypermethylation of collagen type I alpha 2 contributes to proliferation and migration activity of human bladder cancer. Int J Oncol, 2009. [DOI | PubMed]
22. K Misawa, T Kanazawa, Y Misawa, A Imai, S Endo, K Hakamada. Hypermethylation of collagen alpha2 (I) gene (COL1A2) is an independent predictor of survival in head and neck cancer. Cancer Biomark, 2011. [DOI | PubMed]
23. 23.Xu J, Wang Y, Jiang J, Yin C, Shi B. ADAM12 promotes clear cell renal cell carcinoma progression and triggers EMT via EGFR/ERK signaling pathway. J Transl Med. 2023. 10.1186/s12967-023-03913-1.
24. 24.Wang R, Godet I, Yang Y, Salman S, Lu H, Lyu Y, et al. Hypoxia-inducible factor-dependent ADAM12 expression mediates breast cancer invasion and metastasis. Proc Natl Acad Sci USA. 2021. 10.1073/pnas.2020490118.
25. 25.Du S, Sun L, Wang Y, Zhu W, Gao J, Pei W, et al. ADAM12 is an independent predictor of poor prognosis in liver cancer. Sci Rep. 2022. 10.1038/s41598-022-10608-y.
26. MA Eckert, M Santiago-Medina, TM Lwin, J Kim, SA Courtneidge, J Yang. ADAM12 induction by Twist1 promotes tumor invasion and metastasis via regulation of invadopodia and focal adhesions. J Cell Sci, 2017. [DOI | PubMed]
27. 27.Liu H, Tang T. Descriptive pan-cancer genetic analysis of disulfidptosis-related gene set 1. Cancer Genet. 2023. 10.1101/2023.02.25.529997.
28. 28.Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47–e47. 10.1093/nar/gkv007.
29. 29.Coleman A, Bose A, Mitra S. Metagenomics data visualization using R. In: Metagenomic data analysis. New York, NY: Springer USA; 2023. 359–92.
30. Z Gu, D Hübschmann, C Alkan. Make interactive complex heatmaps in R. Bioinformatics, 2022. [DOI | PubMed]
31. 31.Zhang HMP, Davis S. RCircos: an R package for Circos 2D track plots. BMC Bioinforma. 2013;14:244.
32. 32.Friedman JHT, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33:1–22.
33. 33.Wang X, Jiang G, Zong J, Lv D, Lu M, Qu X, et al. Revealing the novel ferroptosis-related therapeutic targets for diabetic foot ulcer based on the machine learning. Front Genet. 2022. 10.3389/fgene.2022.944425.
34. 34.Gao C-H, Yu G, Cai P. ggVennDiagram: an intuitive, easy-to-use, and highly customizable R package to generate Venn diagram. Front Genet. 2021. 10.3389/fgene.2021.706907.
35. G Yu, L-G Wang, Y Han, Q-Y He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A J Integr Biol, 2012. [DOI]
36. S Sun, Y Shen, J Wang, J Li, J Cao, J Zhang. Identification and validation of autophagy-related genes in chronic obstructive pulmonary disease. Int J Chron Obstruct Pulmon Dis, 2021. [DOI | PubMed]
37. 37.Zhou Y, Zhang C, Zhou Z, Zhang C, Wang J. Identification of key genes and pathways associated with PIEZO1 in bone-related disease based on bioinformatics. Int J Mol Sci. 2022. 10.3390/ijms23095250.
38. 38.Sturm G, Finotello F, List M. Immunedeconv: an R package for unified access to computational methods for estimating immune cell fractions from bulk RNA-Sequencing data. In: Bioinformatics for cancer immunotherapy. New York, NY: Springer USA; 2020. p. 223–32.
39. MD Wilkerson, DN Hayes. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics, 2010. [DOI | PubMed]