Identification and validation of a mitochondrial dysfunction-associate

Introduction

Dedifferentiated thyroid cancer (DDTC) is a relatively rare yet highly aggressive thyroid malignancy characterized by the transformation of well-differentiated thyroid carcinoma cells into more malignant phenotypes. Although 85% of differentiated thyroid cancer (DTC) achieve favorable outcomes following active surveillance, surgery, post-operative radioactive iodine therapy, and thyroid-stimulating hormone (TSH) suppression, 5–15% of patients experience locoregional recurrence or distant metastasis—most commonly to the lungs and bones. Among these, one-third lose iodine-avidity either at presentation or during treatment, evolving into dedifferentiated disease, after which 10-year survival drops to <10%.^1–3 The management of DDTC therefore remains formidable: surgery and currently available targeted therapies confer limited efficacy, and robust prognostic tools are lacking. Clinically, several prognostic markers and assessment tools have already been validated for DDTC; nevertheless, their utility remains limited. For instance, BRAF/RAS mutations or TERT promoter mutations cannot independently predict the risk of recurrence or death, the detection rate of some biomarkers is low, and the ATA risk stratification relies solely on the initial surgical pathology, precluding dynamic re-evaluation.⁴ Elucidating the underlying molecular mechanisms and identifying novel biomarkers and therapeutic strategies are thus of paramount importance.

Mitochondrial dysfunction drives dedifferentiation and aggressive progression of DDTC through multiple interconnected mechanisms. Metabolic reprogramming—characterized by suppressed oxidative phosphorylation and enhanced glycolysis—abolishes thyroid-specific markers such as the iodide transporter NIS, rendering the cells radioiodine-refractory and enabling therapeutic escape. Concurrently, accumulation of mitochondrial reactive oxygen species (ROS) stabilizes HIF-1α and activates the epithelial-to-mesenchymal transition (EMT) program, thereby promoting a metastatic phenotype.^5,6 This is further reinforced by apoptotic resistance (eg, up-regulation of BCL-2 family proteins and p53 inactivation) and dysregulated mitochondrial dynamics, exemplified by increased expression of the fission mediator Drp1, which collectively enhance cellular survival and invasive potential.⁷ Ultimately, these alterations propel the tumor toward a poorly differentiated state and facilitate distant metastasis. Given the pivotal role of mitochondrial dysfunction in the onset and progression of DDTC, exploring biomarkers based on mitochondrial dysfunction-related genes (MDRGs) may facilitate the development of practical diagnostic and preventive strategies. Furthermore, the screening of drugs targeting biomarkers associated with mitochondrial dysfunction may provide novel therapeutic insights for the management of DDTC.

Yue Wu et al elucidated the prognostic significance of oxidative-stress-related mitochondrial genes in clear-cell renal cell carcinoma. A six-gene mitochondrial function signature was constructed that predicts patient survival and immunotherapeutic response.⁸ However, analogous studies focusing on mitochondrial-function-related genes in DDTC remain scarce. Building upon the above research, our study incorporates Gene Set Enrichment Analysis (GSEA), prognostic nomogram, competing endogenous RNA (ceRNA) network, drug sensitivity, molecular docking analyses and so on to comprehensively explore the role of MDRGs in the prognosis of DDTC. Despite the lack of functional experimental verification, the study still reinforced the potential of these prognostic genes as targets for future research and therapeutic interventions in DDTC. The study flowchart was illustrated in Figure 1.

Figure 1 The flowchart of this study.

Materials and Methods

Data Sources

The TCGA-THCA cohort was extracted from the TCGA database (https://portal.gdc.cancer.gov). This cohort comprises 505 THCA samples and 59 control samples. The THCA-related datasets (GSE33630 and GSE76039) were obtained from the Gene Expression Omnibus (GEO) database (https://ww.ncbinlm.nih.gov/). Among them, 252 DDTC samples and 59 normal samples were included. The GSE33630 included 11 ATC samples and 45 normal thyroid tissue samples,^9,10 while GSE76039 encompasses 17 PDTC samples and 20 ATC samples.^11,12 Detailed clinical information (including age distribution, gender ratio, and disease stage) and sample information of these cohorts were provided in Supplementary Table 1 and Table 1, respectively. Furthermore, 1674 MDRGs were identified through Genecards Database (http://www.genecards.org) with a relevance score ≥ 5. This database was selected for its comprehensive integration of authoritative sources (eg, NCBI, GO, KEGG) and systematic coverage of mitochondrial pathways. The threshold (≥5) ensured the inclusion of genes with robust associations to mitochondrial dysfunction, referencing a similar strategy.¹³

Table 1 Sample Size Distribution Across Different Groups and Datasets

The Screening of Differentially Expressed Genes (DEGs) Between Dedifferentiated Thyroid Cancer (DDTC) and Control Groups

Based on previously published literature,¹⁴ the expression profiles of 16 genes associated with thyroid metabolism and function, namely DIO1, DIO2, DUOX1, DUOX2, FOXE1, GLIS3, NKXX2-1, PAX8, SLC26A4, SLC5A5, SLC5A8, TG, THRA, THRB, TPO, and TSHR, were obtained and used as indicators to determine the degree of thyroid-specific differentiation. A lower differentiation score indicated a more poorly differentiated thyroid cancer phenotype.¹⁵ The TCGA-THCA samples were subsequently stratified into two groups according to the median differentiation score: high-scoring samples (representing highly differentiated tumors) and low-scoring samples (representing dedifferentiated tumors). The DEGs between dedifferentiated group and the control group were identified via “edgeR” R package (version 3.36.0) with adj P < 0.05 and |log2FC| > 1, which were noted as DEGsl. Afterwards, batch effects between GEO datasets (GSE33630 and GSE76039) from the same platform were corrected using the ComBat algorithm implemented in the “sva” R package (version 3.44.0). The “limma” R package (version 3.48.3) was utilized to screen the DEGs between DDTC and control groups with |log2FC| >1 and adj P < 0.05, which were noted as DEGs2.¹⁶ Finally, the DEGs1 were intersected with DEGs2 to obtain DEGs between DDTC and control groups.

Identification of THCA Related Genes (THCA-RGs)

Using the traits of each sample (THCA and control) from the TCGA-THCA cohort, the “WGCNA” R package (1.70.3) was employed to construct a co-expression network, aiming to identify the gene module most closely associated with DDTC.¹⁷ The samples were clustered to remove outliers. Then, the scale-free network evaluation coefficient R² was set to 0.85. The optimal soft thresholding was obtained when the value of R² first exceeded the critical value of 0.85 and the average connectivity of the co-expression network was close to 0. According to the screened soft threshold, a scale-free network was constructed. The minimum number of genes in each module was set to 100, and the MEDissThres was set to 0.5, so that genes could be effectively grouped into multiple modules using the cutreeDynamic function. To obtain the appropriate soft threshold (β), the network topology analysis was performed to achieve a scale-free topology. Thereafter, the similarity between genes was calculated according to the adjacency, and the gene dendrogram and module color were established using the degree of dissimilarity. The relevance between each gene module and THCA was computed via Pearson correlation analysis (cor > 0.3, P < 0.05). The module with the highest correlation were defined as THCA related module, and genes within this module were defined as THCA-RGs.

Identification and Enrichment Analysis of Differentially Expressed MDRGs (DE-MDRGs)

The DE-MDRGs were gained through taking the intersection of DEGs, THCA-RGs, and MDRGs. To elucidate the biological functions and processes associated with DE-MDRGs in DDTC, enrichment analysis was performed via “clusterProfiler” R package (version 4.0.2), based on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (adj P < 0.05 and count ≥ 1).¹⁸

Creation of the Risk Model

In the TCGA-THCA cohort, after excluding samples with missing survival information, the samples with complete survival information from the dedifferentiation group were randomly classified into training and validation sets according to a 5:5 ratio. Subsequently, univariate Cox analysis was carried out using the “survival” R package (version 3.2.13) to evaluate whether the DE-MDRGs were associated with patients’ survival outcomes (hazard ratio (HR)≠1, P < 0.05). The least absolute shrinkage and selection operator (LASSO) of DE-MDRGs was developed using the “glmnet” R package (version 4.1.4). After performing 10-fold cross-validation, genes with non-zero coefficients at the lambda.min value were selected as prognostic genes. Risk scores for each patient were calculated based on the expression levels of the selected prognostic genes, utilizing the formula: . Where β represents the risk regression coefficient of each prognostic gene, and expi represents expression level of each prognostic gene.

Therefore, each DDTC patient could obtain a risk score by this formula. Patients in both training set and validation set were classified into high-risk and low-risk groups on the basis of median risk score. The risk curves of the high and low risk groups were drawn with the help of the “ggplot2” R package (version 3.3.5), and the heatmap of the gene expression of the model was also plotted. The difference in survival outcomes between two risk groups was analyzed via Kaplan–Meier (K–M) curves by “survminer” R package (version 0.4.9). The predictive capability of risk model was evaluated via receiver operating characteristic (ROC) curves using the “survivalROC” R package (version 1.0.3). Furthermore, samples from the training set and validation set were merged into the entire TCGA-THCA cohort, and the same method was applied to validate the risk model in this entire TCGA-THCA cohort to assess its generalizability.

Creation of the Nomogram

To appraise the significance of the risk model’s impact on prognosis for predicting overall survival (OS) in THCA patients, firstly, clinicopathological factors (including age, gender, T stage, N stage, M stage, and overall stage) as well as risk scores were included in univariate Cox regression analyses (p<0.05, HR≠1) and proportional hazards (PH) assumption tests (p>0.05). Subsequently, factors meeting the thresholds were subjected to multivariate Cox regression analyses (p<0.05, HR≠1). All the above analyses were performed using the “survival” R package (version 3.2.13). Finally, the factors derived therefrom were identified as independent prognostic factors. A nomogram containing independent prognostic factors was created to predict the survival of DDTC patients at 1-, 3- and 5-year using the “rms” R package (version 3.2.13). Furthermore, a corresponding calibration curve, decision curve analysis (DCA) curve, and receiver operating characteristic (ROC) curve were plotted to assess the reliability and accuracy of the nomogram’s predictions.

Gene Set Enrichment Analysis (GSEA)

The differential analysis of TCGA-THCA was carried out between two risk groups. To further explore the differential relevant pathways and molecular mechanisms between the risk groups, the GSEA was performed via “clusterProfiler” R package (version 4.0.2) and “org.Hs.eg.db” R package (version 3.13.0).¹⁸ The “GO:c5.go.v7.5.entrez.gmt” and “KEGG:c2.cp.kegg.v7.5.1.symbols.gmt” acquired from Molecular Signatures Database (MsigDB) (https://www.gsea-msigdb.org/gsea/msigdb) was utilized as background sets. The threshold values were |NES| > 1, NOM P < 0.05, and q < 0.25.

Immune Analysis

The infiltrating abundance of 22 immune cells of each sample was assessed in training set via CIBERSORT algorithm. The differences of infiltrating abundance of immune cells between two risk groups were compared via Wilcox Test (P < 0.05). And the relevance between signature genes and immune cells was evaluated via Spearman algorithm.

Chemotherapy Drug Sensitivity Analysis

The anti-cancer drugs were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC, www.cancerRxgene.org). To assess the sensitivity of anti-cancer drugs in the two risk groups of the training set, the IC50 of anti-cancer drugs was calculated via the “pRRophetic” R package (version 0.2). The differences of IC50 between the two risk groups were compared by Wilcoxon text (P < 0.05).

Drug Prediction of the Prognostic Genes

Each characteristic gene was used as key word to search drugs interacting with these genes in the Drug-Gene Interaction Database (DGIdb, https://dgidb.genome.wustl). Subsequently, molecular docking was performed between the characteristic genes and the drugs with the highest interaction scores. The PDB format files of protein structures of the characteristic genes were downloaded from the RCSB Protein Data Bank (https://www1.rcsb.org/). The SDF format files of the structures of the drugs were obtained from the Drugbank database (https://go.drugbank.com/). The molecular structures were pre-treated via the “AutoDock” Tools. AutoDock vina software (version 4.2.6) molecularly docked the prognostic genes with drugs.¹⁹ PyMol software (version 3.1.1) visualized the docking results.²⁰

Construction of Competing Endogenous RNA (ceRNA) Network

To elucidate the underlying mechanisms of prognostic genes in DDTC, we constructed a ceRNA network. The miRNAs regulating the prognostic genes were predicted via miRWalk Database (http://mirwalk.umm.uni-heidelberg.de/) and miRDB database (http://www.mirdb.org). The predicted miRNAs of the two databases were intersected to obtain the final miRNAs. The miRTarBase database (http://mirtarbase.mbc.nctu.edu.tw) was utilized to predict the lncRNAs with regulatory interactions with miRNAs. Finally, the Cytoscape software (version 3.7.2) was utilized to visualize the ceRNA network.²¹

Cell Culture

The Nthy-ori 3–1 (Nthy) normal human thyroid cell line, as well as the THJ16T, 8505C, Cal-62 ATC cell lines, were purchased from the Institute of Cell Science, Shanghai Academy of Biological Sciences, China. The Nthy normal human thyroid cells were maintained in F12K medium (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% FBS, 100 M/mL streptomycin, 100 U/mL penicillin. The THJ16T, 8505C, Cal-62 cells were cultured in DMEM medium (Gibco, USA) containing 10% fetal bovine serum (FBS, Gibco, USA), 1% penicillin (100 U/mL)/streptomycin (100 mg/mL; Gibco, USA) at 37°C in a humidified incubator of 5% CO2.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

In the combined dataset of GSE33630 and GSE76039, Wilcoxon test was used to analyze the expression differences of the prognostic genes between the DDTC group and the control group (p < 0.05). Then, total RNA was extracted from Nthy, THJ16T, 8505C and Cal-62 cells using TRIzol (Ambion, Austin, USA) in accordance with the manufacturer’s protocol. NanoPhotometer N50 was employed to detect the concentration of RNA. The SureScript-First-strand-cDNA-synthesis-kit (Servicebio, Wuhan, China) was utilized to perform reverse transcription of total RNA to cDNA the based on the manufacturer’s instructions. qRT-PCR was performed utilizing the 2xUniversal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China). The primer sequences for qRT-PCR were shown in Supplementary Table 2. The internal reference gene was GAPDH. The 2−ΔΔCt method was utilized to calculate the expression of prognostic genes.²² The t-test of GraphPad Prism 10 software was employed to analyse qRT-PCR data (p < 0.05).

Statistical Analysis

Statistical analysis was carried out through R software (https://www.r-project.org/). Differences were analyzed via the Wilcoxon test. p < 0.05 represented a significant difference.

Results

Between the DDTC and Control Groups, 796 DEGs Were Obtained

The 505 THCA samples in TCGA-THCA cohort were classified into 253 highly differentiated samples and 252 dedifferentiated samples on the basis of the median value of differentiation score. In total, 3669 DEGs1 were identified between the dedifferentiated group and the control group (Figure 2A and B). As depicted in Figure 2C and D, the data exhibited overall stability following the combined batch processing of GSE33630 and GSE76039, thereby validating the efficacy of the batch processing. Subsequently, after batch correction, 48 DDTC samples and 45 normal samples were obtained for the screening of DEGs between DDTC and control groups. A total of 3395 DEGs2 were detected (Figure 2E and F). Ultimately, the intersection of DEGs1 and DEGs2 yielded 796 DEGs between the DDTC and control groups, comprising 390 upregulated genes and 406 downregulated genes (Figure 2G and H).

Figure 2 Differential expression analysis in the TCGA-THCA, GSE33630, and GSE76039 datasets. (A) Volcano map of DEGs1 between the poorly differentiated group and the control group in the TCGA-THCA dataset. Red indicates up-regulated genes, blue indicates down-regulated genes, and black indicates non-significant genes. (B) Heat map of DEGs1 between the poorly differentiated group and the control group in the TCGA-THCA dataset. (C) Distribution of data from the two datasets before batch processing. (D) Distribution of data from the two datasets after batch processing. (E and F) The volcano map (E) Volcano map of DEGs2 in the merging dataset. Red indicates up-regulated genes, blue indicates down-regulated genes, and black indicates non-significant genes. (F) Heat map of DEGs2 in the merging dataset. (G) Venn diagram of up-regulated genes in TCGA and merged datasets. (H) Venn diagram of down-regulated genes in TCGA and merged datasets.

Identification of DE-MDRGs and Functional Enrichment Analysis

The results of WGCNA showed that there were no outlier samples (Figure 3A). When the optimal soft threshold was 12, the average connectivity of the co-expression network was close to 0 (Figure 3B). A total of 12 co-expression modules were identified (Figure 3C). The brown module was most significantly relevant to THCA (r = −0.82, P = 2e-75) (Figure 3D). A total of 5949 THCA-RGs were identified. Subsequently, the intersection of DEGs with THCA-RGs and MDRGs yielded 41 DE-MDRGs (Figure 3E). Enrichment analysis demonstrated that DE-MDRGs enriched in 477 GO items (Supplementary Table 3). For instance, the DE-MDRGs were involved in “extrinsic apoptotic signaling pathway”, “mitochondrial outer membrane”, “ubiquitin protein ligase binding”, etc. (Figure 3F). Of the KEGG result, 32 pathways were significantly enriched (Supplementary Table 4), including “apoptosis”, “p53 signaling pathway”, “AGE-RAGE signaling pathway in diabetic complications”, etc., all of which were associated with mitochondrial dysfunction (Figure 3G).

Figure 3 Identification of differentially expressed mitochondrial dysfunction related genes (DE-MDRGs) and functional enrichment analysis. (A) Sample clustering diagram. Branches represent samples, and the ordinate represents the height of hierarchical clustering. (B) Selection of the optimal soft-thresholding (power). Left panel: The horizontal axis represents the soft threshold (power), and the vertical axis represents the scale – free topology model fit (R2). A red cutoff line is drawn at R2 = 0.8, which is used to determine the appropriate soft threshold for achieving a scale – free network topology. The red numbers (from 1 to 25) correspond to different soft threshold (power) values. As the soft threshold (power) increases, the scale – free topology model fit rises. When it reaches the red cutoff line, it indicates that the corresponding soft threshold can make the network approximately exhibit a scale – free topology. Right panel: The horizontal axis represents the soft threshold, and the vertical axis represents the mean connectivity; the red numbers correspond to the labels under different soft thresholds, intuitively showing the change in mean connectivity values for each soft threshold. (C) Hierarchical clustering of genes and module identification. (D) Heat map of the relationships between gene modules and clinical traits (THCA). (E) Venn diagram of brown module genes, DEGs, and mitochondrial dysfunction – related genes. Blue represents brown module genes, pink represents DEGs, green represents mitochondrial dysfunction genes, and the overlapping parts represent DE – MDRGs. (F) GO enrichment analysis of DE-MDRGs. (F) KEGG enrichment analysis of DE-MDRGs. Abbreviations: BP, biological process; CC, cellular component; MF, molecular function.

The Prediction Effectiveness of Risk Models Built with DE-MDRGs Was Excellent

Both the training and validation sets comprised 126 dedifferentiated samples with complete survival data, distributed according to a 5:5 ratio. The survival information of patients in the training cohort, validation cohort, high-risk group, and low-risk group is detailed in Supplementary Figure 1A–D. In total, 11 DE-MDRGs were relevant to prognosis of DDTC patients via univariate COX analysis, specifically MTCH1, SLC25A37, SLC26A4, GJB6, SFXN3, DPP4, PMAIP1, KCNQ1, LRP2, TG, and NOX4 (Figure 4A). Subsequently, six prognostic genes (SLC26A4, SLC25A37, KCNQ1, PMAIP1, DPP4 and NOX4) were selected for the construction of a risk prediction model through LASSO regression analysis (lambda.min = 0.0031) (Figure 4B). Following that erected risk model was: risk score = h0(t)*exp((−0.51) * SLC26A4+ 1.59 * SLC25A37 + (−2.79) * KCNQ1 + (1.78) * PMAIP1 + (−0.22)* DPP4 + 0.60 * NOX4). After that, risk model was estimated in training set. Patients in the training set were classified into the high-risk and low-risk groups, each comprising 63 individuals (Figure 4C). The median follow-up durations for the high-risk and low-risk cohorts were 1056 days and 987 days, respectively. Kaplan–Meier survival analysis revealed that patients in the high-risk group exhibited a significantly lower survival rate compared to those in the low-risk group (p = 0.0464) (Figure 4D). Moreover, AUC values of ROC curves for 1-year, 2-year, 3-year, 4-year, and 5-year were 0.998, 0.992, 0.992, 0.992, and 1.000, respectively (Figure 4E). The validation set and entire TCGA-THCA cohort were utilized to verify the model, and the results were concordant with those of the training cohort. The median follow-up durations for the high-risk and low-risk cohorts in validation cohort were 813 days and 957 days, respectively. The Kaplan–Meier curves illustrated that survival probability of high risk patients was shorter than that of low risk patients (p = 0.0495). AUC values of ROC curves for 1-year, 2-year, 3-year, 4-year, and 5-year were 0.654, 0.654, 0.608, 0.700, and 0.717 respectively (Figure 5A–C). Furthermore, in the entire TCGA-THCA cohort, the Kaplan–Meier curve revealed that the survival probability of patients in the high-risk group was notably lower than that in the low-risk group (p=0.031) (Supplementary Figure 2A). ROC curve analysis demonstrated that the AUC values for 1-year, 2-year, 3-year, 4-year, and 5-year were 0.82, 0.888, 0.739, 0.782, and 0.84, respectively (Supplementary Figure 2B). These results demonstrated that the risk model had a better performance of the prognostic prediction of DDTC.

Figure 4 Screening of prognostic genes and construction of risk model. (A) Forest plot of univariate COX regression analysis. (B) LASSO regression analysis. Left: Partial likelihood deviance plot for LASSO regression. Right: Coefficient profile plot for LASSO regression. (C) Distribution plots of risk scores, survival status, and expression of prognostic genes in the training set. Red represents the high-risk group/death cases, and green represents the low-risk group/survival cases. (D) Kaplan–Meier (K–M) survival curves of high- and low-risk groups in the training set. (E) Receiver operating characteristic (ROC) curves of the risk model in the training set. Abbreviation: AUC, area under the curve.

Figure 5 Validation of the risk model. (A) Distribution plots of risk scores, survival status, and expression of prognostic genes in the validation set. Red represents the high-risk group/death cases, and green represents the low-risk group/survival cases. (B) KM survival curves of two risk groups in the validation set. (C) ROC curves of risk model in the validation set.

Robust Nomogram Built Using Independent Prognostic Factors

In total, two independent prognostic factors were identified: the risk score and pathologic T stage (Figure 6A and B, Supplementary Table 5). As shown in Figure 6C, a nomogram incorporating both the risk score and pathologic T stage was created. The concordance index (C-index) of the calibration curve was 0.927, indicating excellent agreement between predicted and observed outcomes (Figure 6D). Additionally, the AUC values of the nomogram at 1-, 3- and 5-year were 0.88, 0.83 and 0.93, respectively (Supplementary Figure 3A–C). Decision curve analysis (DCA) revealed that the net benefit of the nomogram exceeded that of the individual factors, thereby further highlighting the robust predictive capacity of the nomogram (Supplementary Figure 4). These results collectively suggested that the nomogram exhibited a high degree of accuracy in prognostic prediction for patients with DDTC.

Figure 6 Independent prognostic analysis. (A) Forest plot of univariate Cox regression analysis (B) Forest plot of multivariate Cox regression analysis. (C) Construction of the nomogram based on the independent prognostic factors. (D) Calibration curve of the nomogram.

Exploration of the Potential Functions of High and Low Risk Groups

The top10 GO items and KEGG pathways with significant enrichment were illustrated in Figure 7A and B. Of the GO results, a total of 697 GO items were markedly enriched (Supplementary Table 6). Specifically, terms such as “extracellular matrix structural constituent”, “antimicrobial humoral response”, “collagen containing extracellular matrix”, “keratinization”, “cytokine activity” were prominently enriched in the high-risk group (Figure 7A).

Figure 7 Gene Set Enrichment Analysis (GSEA) of high and low risk groups. (A) GO enrichment analysis. (B) KEGG enrichment analysis.

Regarding the KEGG results, a total of 33 KEGG pathways were significantly enriched (Supplementary Table 7). Notable pathways enriched in the high-risk group included “cytokine receptor interaction”, “ECM receptor interaction”, “hypertrophic cardiomyopathy HCM” (Figure 7B). Conversely, in the low-risk group, pathways such as “oxidative phosphorylation” and “ribosome” were markedly enriched (Figure 7B).

Disturbances of Immune Microenvironment Might Be Associated with the Development of DDTC

Four differential immune cells were obtained between the high-risk and low-risk groups, including activated dendritic cells (DC), activated CD4 memory T cells, plasma cells, and regulatory T cells (Tregs) (Figure 8A and B) Correlation analysis revealed that SLC26A4, SLC25A37, KCNQ1, PMAIP1, and NOX4 were significantly associated with the abundance of ADC (Figure 8C). Notably, SLC26A4 and KCNQ1 demonstrated a marked negative correlation with ADC, while the remaining three genes exhibited an opposing trend, showing positive correlation with activated dendritic cells. Additionally, PMAIP1 was markedly positively related to activated CD4 memory T cells, and negatively related to plasma cells. In addition, both PMAIP1 and DPP4 were found to be significantly positively relevant to Tregs.

Figure 8 Immune infiltration analysis. (A) Stacked bar plot of 22 immune cell proportions. (B) Comparison of immune cells score in samples from high and low risk groups. Red represents the high – risk group, and blue represents the low – risk group. ns, not significant; *p<0.05; **p<0.01; ***p<0.001. (C) Heatmap of correlations between immune cells and prognostic genes. Red indicates a positive correlation, and blue indicates a negative correlation. *p<0.05; **p<0.01.

Predicted Drugs Might Contribute to the Treatment of DDTC

In total, the IC50 of 67 anti-cancer drugs exhibited significantly different between two risk groups. Specifically, 40 drugs demonstrated increased sensitivity in the high-risk group, while 27 drugs showed greater sensitivity in the low-risk group (Supplementary Tables 8 and 9).

A total of 41 drugs were predicted to interact with the prognostic genes (Figure 9A). 1, 1, 13, 1, 24, and 1 drugs were predicted for SLC26A4, SLC25A37, KCNQ1, PMAIP1, DPP4 and NOX4 respectively. The interaction scores of the following pairs were the highest: KCNQ1-indomethacin, PMAIP1-bortezomib, DPP4-sitagliptin, SLC25A37-sodium citrate, SLC26A4-iodine, and NOX4-dextromethorphan. These pairs were selected for molecular docking analysis (Figure 9B–G). Typically, a docking affinity score of less than −4.25 kcal/mol indicates a standard binding capacity, while a score less than −5 kcal/mol suggests a stable binding, and a score less than −7 kcal/mol indicates a strong binding capacity. The docking affinity score of KCNQ1-indomethacin, PMAIP1-bortezomib, DPP4-sitagliptin, SLC25A37-sodium citrate, SLC26A4-iodine, and NOX4-dextromethorphan were −7.1 kcal/mol, −5.5 kcal/mol, −7.3 kcal/mol, −4.4 kcal/mol, −1.8 kcal/mol, and −9.1 kcal/mol, respectively (Table 2). These results indicated that the combinations of KCNQ1 with indomethacin, PMAIP1 with bortezomib, DPP4 with sitagliptin, SLC25A37 with sodium citrate, and NOX4 with dextromethorphan exhibit stable binding interactions.

Table 2 Molecular Docking Results

Figure 9 Drug prediction and molecular docking. (A) Construction of the gene-drug network. (B–G) The results of molecular docking based on prognostic genes and their targeted drugs.

Investigation of Potential Regulatory Mechanisms of Characteristic Genes

Through miRWalk and miRDB databases, a total of 99 miRNAs were predicted for prognostic genes. Additionally, 52 lncRNAs were identified as potential regulators of these prognostic genes. Finally, the lncRNA-miRNA-mRNA network was constructed, comprising 6 mRNAs, 99 miRNAs, and 53 lncRNAs (Figure 10). For example, IGF2-AS targeting hsa-miR-197-3p, which in turn regulates KCNQ1 and CYB561D2 targeting hsa-miR-26b-5p, which regulates SLC26A4, among others.

Figure 10 LncRNA-miRNA-mRNA regulatory network. Blue represents lncRNA, green represents miRNA, and pink represents mRNA.

Expression Validation of Prognostic Genes

In the combined dataset of GSE33630 and GSE76039, the expression of SLC26A4, SLC25A37, and KCNQ1 was significantly lower in DDTC cells than that in the cells of the control group (Figure 11A). To further validate the expression levels of prognostic genes, we performed quantitative real-time polymerase chain reaction (qRT-PCR). The results showed that the expression of SLC26A4, SLC25A37, and KCNQ1 was significantly downregulated in DDTC cells (THJ16T, 8505C, Cal-62) compared to that in the normal thyroid cell line (Nthy) (Figure 11B–D). Conversely, the expression of PMAIP1, DPP4, and NOX4 was significantly upregulated in DDTC cell lines (Figure 11E–G). Moreover, the expression trend of SLC26A4, KCNQ1, PMAIP1, DPP4, and NOX4 was consistent with the data obtained from public databases.

Figure 11 Differential expression of six prognostic genes in the GEO dataset (A) Validation of prognostic gene expression levels in the merged dataset from GEO datasets. ***p<0.001, NS represents not significant. (B–G). qRT-PCR validation. The horizontal axis represents cell lines, and the vertical axis represents the relative expression levels of prognostic genes.*p<0.05; **p<0.01; ***p<0.001; ****p<0.0001.

Discussion

Mitochondria serve as the primary controllers of the intrinsic apoptotic pathway and play a crucial role in cell survival and the production of cellular energy.²³ Accumulating evidence suggests that mitochondrial dysfunction is implicated in the development and progression of DDTC. However, the role of MDRGs in DDTC remains underexplored. Zhu et al²⁴ previously reported that the expression of mitochondrial ROS regulators was associated with the clinical features and progression-free survival (PFS) in patients with THCA. While our study specifically focused on DDTC and comprehensively investigated genes related to mitochondrial dysfunction.

We utilized datasets from TCGA and GEO to explore the potential of MDRGs as prognostic biomarkers for DDTC and to identify potential therapeutic targets for this disease. Our research identified six MDRGs that were most significantly associated with the prognosis of DDTC, specifically SLC26A4, SLC25A37, KCNQ1, PMAIP1, DPP4 and NOX4. Furthermore, qRT-PCR confirmed that the expression patterns of these genes were consistent with those observed in public databases.

SLC26A4, an iodine transporter located on the apical membrane of thyroid cells, is primarily involved in iodine transport and maintenance of thyroid function.²⁵ Studies have shown that SLC26A4 mutation can cause iodine transport disorders, thereby disrupting thyroid hormone synthesis.²⁶ Moreover, SLC26A4 is also involved in regulating the differentiation and proliferation of thyroid follicular cells.²⁷ In metastatic papillary thyroid cancer (PTC), SLC26A4 expression is reduced compared to normal thyroid tissue, positioning it as a key risk factor for PTC development through its impact on iodine metabolism.²⁸ Furthermore, SLC26A4 expression has been found to be significantly downregulated in DDTC compared to DTC, suggesting that the limited uptake and rapid excretion of radioactive iodine in DDTC may be attributed to the absence of functional apical transporters.²⁹

The protein encoded by the SLC25A37 gene is an integral component of the mitochondrial inner membrane, primarily functioning as an iron ion transporter involved in mitochondrial iron metabolism.³⁰ Thyroid hormone synthesis necessitates iron as a cofactor, facilitating the iodination process catalyzed by thyroid peroxidase (TPO).³¹ Consequently, SLC25A37 holds a pivotal role in regulating the synthesis of thyroid hormone precursors. Notably, our research has demonstrated an upregulation in the expression level of SLC25A37 in thyroid cancer cells, potentially attributed to the escalated demand for iron ions within these neoplastic cells.³²

Abnormal expression of potassium channels has been implicated in a range of human cancers.^33,34 Our study identified KCNQ1 as a prognostic biomarker for DDTC. The KCNQ1 gene encodes the Kv7.1 channel protein, which regulates voltage-dependent potassium ion flux in diverse cell types. In thyroid cells, the expression of KCNQ1 is essential for maintaining cell membrane potential and thyroid hormone release. Studies have found that KCNQ1 is related to thyroid tumorigenesis, and alterations in its expression level may affect the proliferation and differentiation of thyroid cells.^35,36 Therefore, KCNQ1 may function as a regulator during thyroid differentiation.

PMAIP1, a member of the pro-apoptotic subfamily of the BCL-2 protein family and a recognized p53 target, plays an important role in cell stress response and apoptosis signaling pathways. It has been identified as a potential oncogene in pancreatic cancer cell lines. Conversely, in lung cancer cells, PMAIP1 induces apoptosis.^37,38 Additionally, PMAIP1 expression has been reported to be downregulated in gastric cancer with bone metastasis.³⁹ In the context of thyroid differentiation, PMAIP1 affects the survival and function of thyroid cells by regulating cell cycle and apoptosis. Studies have shown that the expression of PMAIP1 in thyroid cancer cells is closely related to the aggressiveness and prognosis of cancer.⁴⁰ These findings suggest that PMAIP1 may serve as a key regulator in the differentiation and malignant transformation of thyroid cells.

NOX4 is an enzyme that produces ROS, which plays a crucial role in cell signaling and regulation of cell differentiation.⁴¹ In thyroid differentiation, NOX4 affects multiple signaling pathways by regulating ROS levels, including the MAPK/ERK and PI3K/Akt pathways. These pathways are critical for thyroid cell proliferation and differentiation.^42,43 The increased activity of NOX4 leads to increased ROS levels, which disrupt the REDOX balance of thyroid cells and consequently alter their differentiation state.⁴⁴ Consistent with our findings, several studies have reported increased NOX4 expression in thyroid cancer. It may be related to the enhanced cancer cell proliferation and survival.⁴⁵

Subsequently, we developed a prognostic nomogram based on these 6 MDRGs and clinicopathological features. Some researchers established and validated prognostic nomograms of THCA by analyzing the Surveillance, Epidemiology, and End Results (SEER) databases. However, these analyses did not include genomic data.^46,47 Liu et al⁴⁸ developed a nomogram on the basis of genes related to insulin-like growth factor signaling pathway in THCA, though this model lacked detailed pathological classification. In contrast, we constructed a nomogram based on mitochondrial dysfunction-related signature in conjunction with other clinicopathological features. This nomogram demonstrated robust predictive performance for the prognosis of DDTC patients. However, several critical issues warrant attention regarding survival analysis and predictive performance assessment based on the risk model. Firstly, in the low-risk group, all subjects remained event-free throughout the follow-up period, with no events occurring after 4000 days. This observation may be attributable to inadequate follow-up duration, which could introduce bias and limit the robustness of the findings. To enhance the scientific rigor and reliability of our study, we plan to continuously gather additional independent thyroid cancer datasets from diverse regions and research institutions, beyond the existing datasets, in future endeavors. This approach will facilitate the construction of a larger external validation cohort, thereby enabling a comprehensive evaluation of the model’s predictive performance across varied populations. Secondly, the AUC value of the 5-year ROC curve was 1 in the training set but 0.717 in the validation set. This discrepancy likely stems from differences in clinical characteristics, sample size, follow-up duration, and the number of positive events between the two cohorts, which collectively impede the model’s generalizability. To ensure the reliability of our results, the risk model was meticulously constructed using prognostic features derived from univariate Cox regression analysis and LASSO regression analysis. During the LASSO regression analysis, we employed 10-fold cross-validation to bolster the model’s stability and reliability. The integration of LASSO analysis with 10-fold cross-validation is instrumental in mitigating overfitting of the AUC value, thereby enhancing the model’s predictive accuracy and generalizability. Furthermore, we intend to conduct further validation of the model’s performance in future external validation cohorts.⁴⁹

The mitochondrion is a critical organelle responsible for energy production and involvement in numerous biological processes. Dysfunction of mitochondria can result in the disruption of tricarboxylic acid (TCA) cycle enzymes, leakage from the electron transport chain, and subsequent oxidative stress, which in turn alter cellular metabolism and signal transduction pathways, ultimately leading to apoptosis and therapeutic resistance.⁵⁰ In order to investigate relevant pathways and molecular mechanisms of mitochondrial dysfunction in DDTC, we conducted GSEA enrichment analysis for high-risk and low-risk groups. The results showed that there were significant differences in cytokine receptor, chemokine signaling pathway, oxidative phosphorylation and other functions between the two groups. These findings provided a foundation for further experimental research.

Previous studies have illuminated the significant role of immune cell infiltration in the tumor microenvironment.⁵¹ Fei et al⁵² pioneered the construction of a novel immune-related prognostic risk model for THCA patients and validated it in the TCGA cohort. The risk model had general applicability and high stability to predict the prognosis of THCA patients. Zhang et al⁵³ investigated Estrogen-related genes associated with immune infiltration in thyroid cancer, revealing significant upregulation of macrophages in the high-risk group. Thus, a deeper understanding of the immune function changes involved in the development of THCA may enhance prognostic prediction and facilitate precise, personalized therapies. The CIBERSORT algorithm which quantifies the abundance of 22 immune cell types within our training set, has provided a detailed characterization of immune cell infiltration in differentiated thyroid carcinoma (DDTC). This analysis revealed substantial differences in the infiltration of ADCs, activated memory CD4+ T cells, plasma cells, and regulatory T cells (Tregs) between different risk groups. These findings suggest that these particular immune cells may play pivotal roles in the progression of DDTC and could potentially serve as biomarkers for disease prognosis or therapeutic targets. For instance, significant correlations identified between prognostic genes (SLC26A4, SLC25A37, KCNQ1, PMAIP1, and NOX4) and ADCs underscore the intricate interplay between the genetic landscape and the immune milieu. The negative correlation of SLC26A4 and KCNQ1 with ADCs, in contrast to the positive association of the other genes, suggests a complex regulatory mechanism governing immune responses within the tumor environment. This could have implications for understanding the mechanisms of immune evasion and the development of novel immunotherapies.

Furthermore, we have constructed a lncRNA-miRNA-mRNA network associated with mitochondrial dysfunction in DDTC, laying the groundwork for future mechanistic investigations. Studies have elucidated that tumor suppressor microRNAs, including miR-204-5p, miR-202-3p, and miR-26b-5p, can mitigate cell invasion and migration in THCA.^44,54 Insulin-like growth factor 2 antisense RNA (IGF2-AS) is a long non-coding RNA that plays a crucial role in gene regulation and has been implicated in various cancers.⁵⁵ In the context of DDTC, IGF2-AS may act as a lncRNA, sequestering miRNAs such as hsa-miR-197-3p, thereby modulating the expression of oncogenes or tumor suppressor genes like KCNQ1. The interaction between IGF2-AS and miR-197-3p could be pivotal in understanding the molecular mechanisms underlying DDTC progression.

The exploration of potential therapeutic agents for DDTC has yielded promising insights, particularly regarding the interactions between specific drugs and prognostic genes identified in this study. Among the drugs analyzed, dextromethorphan, which exhibited a docking score of −9.1 kcal/mol with NOX4, is primarily known as a cough suppressant. However, emerging evidence suggests that it may also exert anti-cancer effects through the modulation of oxidative stress and inflammation.⁵⁶ The strong binding affinity with NOX4, a gene implicated in reactive oxygen species production, indicates that dextromethorphan could play a role in mitigating oxidative stress in DDTC, thereby influencing tumor progression and response to therapy. The findings suggest that existing medications may be repurposed to target specific molecular pathways involved in DDTC, paving the way for novel therapeutic strategies that could improve patient outcomes.

The present study has several limitations that warrant consideration. Firstly, although we utilized multiple datasets to enhance the robustness of our results, the relatively small sample size and the potential for batch effects and confounding variables remains a concern. Secondly, the prognostic model developed in this study, while demonstrating strong predictive capabilities, requires external validation in independent cohorts to confirm its generalizability and clinical applicability. Furthermore, the study primarily focused on the expression levels of MDRGs without functional experiment validation for these bioinformatics results raised concerns regarding the reproducibility and biological relevance of the findings. Future research should aim to provide both external cohort validation and functional experiments to confirm its prognostic robustness and mechanistic underpinnings.

Conclusion

Clinical prediction models based on these 6 MDRGs may improve prognostic precision, enable the stratification of high-risk patients, and provide clinicians with actionable guidance for personalized therapeutic planning. The exploration of immune microenvironment disturbances and the identification of potential therapeutic agents targeting MDRGs provide promising avenues for future research aimed at improving the management of DDTC. Collectively, this study identified a gene signature associated with mitochondrial dysfunction, furnishing potential biomarkers for prognostic assessment in DDTC. Despite inherent limitations, our findings delineate a roadmap for future investigations. We anticipate that external cohort validation and functional experimentation will corroborate these results and elucidate their clinical translational potential.

Ethics Approval

The studies involving human participants were approved by the Ethics Committees of the Affiliated Hospital of Qingdao University (QYFYWZLL27846).

Disclosure

The authors report no conflicts of interest in this work.

References

1. Toraih EA, Fawzy MS, Ning B, et al. A miRNA-based prognostic model to trace thyroid cancer recurrence. Cancers. 2022;14(17):4128. doi:10.3390/cancers14174128

2. Ruchong P, Haiping T, Xiang W. A five-gene prognostic nomogram predicting disease-free survival of differentiated thyroid cancer. Dis Markers. 2021;2021:5510780. doi:10.1155/2021/5510780

3. Ma B, Jiang H, Wen D, et al. Transcriptome analyses identify a metabolic gene signature indicative of dedifferentiation of papillary thyroid cancer. J Clin Endocrinol Metab. 2019;104(9):3713–3725. doi:10.1210/jc.2018-02686

4. Wang Y, Yang J, Chen S, Wang W, Teng L. Identification and validation of a prognostic signature for thyroid cancer based on ferroptosis-related genes. Genes. 2022;13. doi:10.3390/genes13060997

5. Wang SF, Tseng LM, Lee HC. Role of mitochondrial alterations in human cancer progression and cancer immunity. J Biomed Sci. 2023;30(1):61. doi:10.1186/s12929-023-00956-w

6. Foo BJ, Eu JQ, Hirpara JL, Pervaiz S. Interplay between mitochondrial metabolism and cellular redox state dictates cancer cell survival. Oxid Med Cell Longev. 2021;2021(1):1341604. doi:10.1155/2021/1341604

7. Dabravolski SA, Nikiforov NG, Zhuravlev AD, et al. The role of altered mitochondrial metabolism in thyroid cancer development and mitochondria-targeted thyroid cancer treatment. Int J Mol Sci. 2021;23(1):460. doi:10.3390/ijms23010460

8. Yi HS, Chang JY, Kim KS, Shong M. Oncogenes, mitochondrial metabolism, and quality control in differentiated thyroid cancer. Korean J Intern Med. 2017;32(5):780–789. doi:10.3904/kjim.2016.420

9. Tomás G, Tarabichi M, Gacquer D, et al. A general method to derive robust organ-specific gene expression-based differentiation indices: application to thyroid cancer diagnostic. Oncogene. 2012;31(41):4490–4498. doi:10.1038/onc.2011.626

10. Dom G, Tarabichi M, Unger K, et al. A gene expression signature distinguishes normal tissues of sporadic and radiation-induced papillary thyroid carcinomas. Br J Cancer. 2012;107(6):994–1000. doi:10.1038/bjc.2012.302

11. Landa I, Ibrahimpasic T, Boucai L, et al. Genomic and transcriptomic hallmarks of poorly differentiated and anaplastic thyroid cancers. J Clin Invest. 2016;126(3):1052–1066. doi:10.1172/jci85271

12. Wen S, Qu N, Ma B, et al. Cancer-associated fibroblasts positively correlate with dedifferentiation and aggressiveness of thyroid cancer. Onco Targets Ther. 2021;14:1205–1217. doi:10.2147/ott.S294725

13. Chen Y, Bi K, Zhang C, et al. Identification of endoplasmic reticulum stress and mitochondrial dysfunction related biomarkers in osteoporosis. Hereditas. 2025;162(1):21. doi:10.1186/s41065-025-00387-7

14. Xu W, Li C, Ma B, et al. Identification of key functional gene signatures indicative of dedifferentiation in papillary thyroid cancer. Front Oncol. 2021;11:641851. doi:10.3389/fonc.2021.641851

15. Liu T, Wang J, Xiu Y, Wu Y, Xu D. DNA methylation age drift is associated with poor outcomes and de-differentiation in papillary and follicular thyroid carcinomas. Cancers. 2021;13(19):4827. doi:10.3390/cancers13194827

16. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. doi:10.1093/nar/gkv007

17. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. 2008;9(1):559. doi:10.1186/1471-2105-9-559

18. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. 2021;2(3):100141. doi:10.1016/j.xinn.2021.100141

19. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem. 2010;31(2):455–461. doi:10.1002/jcc.21334

20. Seeliger D, de Groot BL. Ligand docking and binding site analysis with PyMOL and Autodock/Vina. J Comput Aided Mol Des. 2010;24(5):417–422. doi:10.1007/s10822-010-9352-6

21. Otasek D, Morris JH, Bouças J, Pico AR, Demchak B. Cytoscape Automation: empowering workflow-based network analysis. Genome Biol. 2019;20(1):185. doi:10.1186/s13059-019-1758-4

22. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25(4):402–408. doi:10.1006/meth.2001.1262

23. Wang W, Sun Q, Wu Z, et al. Mitochondrial dysfunction-related genes in hepatocellular carcinoma. Front Biosci. 2013;18(3):1141–1149. doi:10.2741/4169

24. Zhu J, Zheng X, Lu D, Zheng Y, Liu J. Clinical relevance and tumor growth suppression of mitochondrial ROS regulators along NADH: Ubiquinone oxidoreductase subunit B3 in thyroid cancer. Oxid Med Cell Longev. 2022;2022(1):8038857. doi:10.1155/2022/8038857

25. Everett LA, Glaser B, Beck JC, et al. Pendred syndrome is caused by mutations in a putative sulphate transporter gene (PDS). Nat Genet. 1997;17(4):411–422. doi:10.1038/ng1297-411

26. Kopp P. Pendred’s syndrome and genetic defects in thyroid hormone synthesis. Rev Endocr Metab Disord. 2000;1(1–2):109–121. doi:10.1023/a:1010024722595

27. Pesce L, Bizhanova A, Caraballo JC, et al. TSH regulates pendrin membrane abundance and enhances iodide efflux in thyroid cells. Endocrinology. 2012;153(1):512–521. doi:10.1210/en.2011-1548

28. Zhang Q, Li J, Shen H, Bai X, Zhang T, Liu P. Screening and validation of lymph node metastasis risk-factor genes in papillary thyroid carcinoma. Front Endocrinol. 2022;13:991906. doi:10.3389/fendo.2022.991906

29. Colombo C, Minna E, Gargiuli C, et al. The molecular and gene/miRNA expression profiles of radioiodine resistant papillary thyroid cancer. J Exp Clin Cancer Res. 2020;39(1):245. doi:10.1186/s13046-020-01757-x

30. Shaw GC, Cope JJ, Li L, et al. Mitoferrin is essential for erythroid iron assimilation. Nature. 2006;440(7080):96–100. doi:10.1038/nature04512

31. Anderson GJ. Mechanisms of iron loading and toxicity. Am J Hematol. 2007;82(12 Suppl):1128–1131. doi:10.1002/ajh.21075

32. Mancias JD, Wang X, Gygi SP, et al. Quantitative proteomics identifies NCOA4 as the cargo receptor mediating ferritinophagy. Nature. 2014;509(7498):105–109. doi:10.1038/nature13148

33. Peroz D, Rodriguez N, Choveau F, et al. Kv7.1 (KCNQ1) properties and channelopathies. J Physiol. 2008;586(7):1785–1789. doi:10.1113/jphysiol.2007.148254

34. Zhang P, Yang X, Yin Q, et al. Inhibition of SK4 potassium channels suppresses cell proliferation, migration and the epithelial-mesenchymal transition in triple-negative breast cancer cells. PLoS One. 2016;11(4):e0154471. doi:10.1371/journal.pone.0154471

35. Nakajo K. Gating modulation of the KCNQ1 channel by KCNE proteins studied by voltage-clamp fluorometry. Biophys Physicobiol. 2019;16:121–126. doi:10.2142/biophysico.16.0_121

36. Morokuma J, Blackiston D, Adams DS, et al. Modulation of potassium channel function confers a hyperproliferative invasive phenotype on embryonic stem cells. Proc Natl Acad Sci U S A. 2008;105(43):16608–16613. doi:10.1073/pnas.0808328105

37. Innamaa A, Jackson L, Asher V, et al. Expression and effects of modulation of the K2P potassium channels TREK-1 (KCNK2) and TREK-2 (KCNK10) in the normal human ovary and epithelial ovarian cancer. Clin Transl Oncol. 2013;15(11):910–918. doi:10.1007/s12094-013-1022-4

38. Than BL, Goos JA, Sarver AL, et al. The role of KCNQ1 in mouse and human gastrointestinal cancers. Oncogene. 2014;33(29):3861–3868. doi:10.1038/onc.2013.350

39. Ishida M, Sunamura M, Furukawa T, et al. The PMAIP1 gene on chromosome 18 is a candidate tumor suppressor gene in human pancreatic cancer. Dig Dis Sci. 2008;53(9):2576–2582. doi:10.1007/s10620-007-0154-1

40. Oda E, Ohki R, Murasawa H, et al. Noxa, a BH3-only member of the Bcl-2 family and candidate mediator of p53-induced apoptosis. Science. 2000;288(5468):1053–1058. doi:10.1126/science.288.5468.1053

41. Bedard K, Krause KH. The NOX family of ROS-generating NADPH oxidases: physiology and pathophysiology. Physiol Rev. 2007;87(1):245–313. doi:10.1152/physrev.00044.2005

42. Saito Y, Maillet M, watt AJ, et al. NOX4 regulates G2/M phase progression during endothelial cell cycle. Circ Res. 2006;98(6):837–839. doi:10.1161/01.RES.0000215985.18538.c4

43. Block K, Gorin Y. Aiding and abetting roles of NOX oxidases in cellular transformation. Nat Rev Cancer. 2012;12(9):627–637. doi:10.1038/nrc3339

44. Van Branteghem C, Augenlicht A, Demetter P, Craciun L, Maenhaut C. Unraveling the roles of miR-204-5p and HMGA2 in papillary thyroid cancer tumorigenesis. Int J Mol Sci. 2023;24(13):10764. doi:10.3390/ijms241310764

45. Lyle AN, Deshpande NN, Taniyama Y, et al. Poldip2, a novel regulator of Nox4 and cytoskeletal integrity in vascular smooth muscle cells. Circ Res. 2009;105(3):249–259. doi:10.1161/CIRCRESAHA.109.193722

46. Yu JW, Pang R, Liu B, Zhang L, Kong LY. Development and validation of a nomogram model for predicting distant metastasis of aged ≥50 patients with thyroid carcinoma: a SEER database analysis. Eur Rev Med Pharmacol Sci. 2024;28(6):2351–2362. doi:10.26355/eurrev_202403_35742

47. Li Y, Zhang H, Cao Y, et al. Establishment and verification of the first prognostic nomograms in locally advanced thyroid cancer based on the analysis of clinical and follow-up information on 2396 patients. Heliyon. 2024;10(3):e24798. doi:10.1016/j.heliyon.2024.e24798

48. Liu J, Miao X, Yao J, et al. Investigating the clinical role and prognostic value of genes related to insulin-like growth factor signaling pathway in thyroid cancer. Aging. 2024;16(3):2934–2952. doi:10.18632/aging.205524

49. Sun X, Wu C, Kang J, Lv H, Liu X. Development and validation of a risk prediction model for short-term progression of carotid atherosclerosis among early middle age adults. Atherosclerosis. 2024;397:118557. doi:10.1016/j.atherosclerosis.2024.118557

50. Luo Y, Ma J, Lu W. The significance of mitochondrial dysfunction in cancer. Int J Mol Sci. 2020;21(16):5598. doi:10.3390/ijms21165598

51. Fridman WH, Pagès F, Sautès-Fridman C, Galon J. The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer. 2012;12(4):298–306. doi:10.1038/nrc3245

52. Fei H, Han X, Wang Y, Li S. Mining prognostic biomarkers of thyroid cancer patients based on the immune-related genes and development of a reliable prognostic risk model. Mediators Inflamm. 2023;2023:6503476. doi:10.1155/2023/6503476

53. Zhang L, Zhou M, Gao X, et al. Estrogen-related genes for thyroid cancer prognosis, immune infiltration, staging, and drug sensitivity. BMC Cancer. 2023;23(1):1048. doi:10.1186/s12885-023-11556-0

54. Chen J, Yin J, Liu J, et al. MiR-202-3p functions as a tumor suppressor and reduces cell migration and invasion in papillary thyroid carcinoma. Eur Rev Med Pharmacol Sci. 2019;23(3):1145–1150. doi:10.26355/eurrev_201902_17005

55. Tian Y, Han W, Fu L, Zhang J, Zhou X. IGF2 is upregulated by its antisense RNA to potentiate pancreatic cancer progression. Funct Integr Genomics. 2023;23(4):348. doi:10.1007/s10142-023-01277-9

56. Ommati MM, Jamshidzadeh A, Saeed M, Rezaei M, Heidari R. Dextromethorphan improves locomotor activity and decreases brain oxidative stress and inflammation in an animal model of acute liver failure. Clin Exp Hepatol. 2022;8(3):178–187. doi:10.5114/ceh.2022.118299