Thyroid differentiation score-related genes and prognostic model for thyroid cancer

Introduction

Thyroid cancer (THCA) is the most common endocrine malignancy, with an increasing incidence worldwide (1). Despite the generally favorable prognosis associated with differentiated thyroid cancer (DTC), recurrent or progressive disease remains the major challenge in the management of anaplastic thyroid cancer (ATC) (2,3). Current conventional treatments, including radiotherapy and chemotherapy, provide limited survival benefits (4). Therefore, there is a pressing need for more precise and personalized prognostic tools that can guide clinical decision-making and improve patient outcomes.

One promising approach in the prognostication of THCA involves the use of molecular and genetic markers. Recent advances in high-throughput sequencing and bioinformatics have enabled the identification of numerous molecular alterations associated with THCA, including mutations in genes such as BRAF, RAS, and TERT promoter (5-7), as well as chromosomal rearrangements involving RET/PTC and PAX8/PPARγ (8). These genetic alterations have been linked to specific clinicopathological features and outcomes (9). However, the clinical utility of these markers as independent prognostic tools remains limited, as their predictive power often depends on the context of other clinical and pathological factors (10).

The thyroid differentiation score (TDS), proposed by members of The Cancer Genome Atlas Research Network (11), is an emerging molecular tool that quantifies the extent of thyroid-specific gene expression, reflecting the degree of differentiation of THCA (12). Tumor differentiation is a critical determinant of prognosis in THCA, with well-differentiated tumors generally associated with better outcomes. The TDS integrates the expression levels of key thyroid-specific genes such as thyroglobulin, thyroid peroxidase, and sodium-iodide symporter (13). A study has shown that a lower TDS, indicating dedifferentiation, is associated with more aggressive tumor behavior, resistance to radioactive iodine therapy, and poorer clinical outcomes (12). Given its potential prognostic value, the TDS can serve as a foundation for constructing a prognostic model that integrates both molecular and clinical factors to predict outcomes in THCA more accurately. Moreover, identifying prognostic biomarkers based on the TDS may provide new therapeutic targets, particularly for patients with poorly differentiated or ATC (14).

By leveraging large-scale genomic data and advanced statistical methods, this study aimed to construct a prognostic model for THCA by incorporating the TDS and identifying novel prognostic biomarkers associated with thyroid differentiation. This approach has the potential to significantly improve patient management and outcomes, particularly for those with high-risk disease. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-460/rc).

Methods

Data acquisition

We downloaded the RNA sequencing (RNA-seq) data with clinical information for THCA patients from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/). A total of 568 samples were enrolled in this study, including 510 disease samples and 58 normal samples, and the messenger RNA (mRNA) expression profiles from the TCGA-THCA cohort were extracted for further analysis. Besides, gene expression data in the GSE33630 dataset (containing 60 disease samples and 45 normal samples) was acquired from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) and served as a validation set. Single-cell sequencing data were acquired from the GSE191288, including 3 THCA samples and 1 normal sample. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Acquisition of TDS

We calculated each sample’s TDS using the mRNA expression levels of 16 thyroid function-related genes (DIO1, DIO2, DUOX1, DUOX2, FOXE1, GLIS3, PAX8, NKX2-1, SLC26A4, SLC5A5, SLC5A8, TG, TPO, TSHR, THRA, and THRB) (13) from the TCGA-THCA dataset. The TDS is a metric to evaluate the association between thyroid differentiation and genetic or epigenetic alterations. We first centered the log2-normalized RNA-seq by expectation maximization (RSEM) values at the median across all samples to calculate the TDS, resulting in the log2 fold change (FC). This value was then averaged across the 16 genes for each sample in TCGA-THCA cohort:

TDS=MeanofLog2FCacrossthese16genes[1]

Based on the median value of TDS, patients in the TCGA-THCA were divided into high TDS group or low TDS group.

Identification of differentially expressed genes (DEGs)

DEGs between the normal and THCA samples or DEGs between the high TDS and low TDS groups in the TCGA-THCA cohort were identified using the Limma package in R (P<0.05, FC =1.5). The DEGs were visualized through a volcano plot using the R package ggplot2. Then, the overlapping DEGs related to THCA and TDS were identified through a Venn diagram and used for further analyses. The Venn diagram was drawn using the R package VennDiagram.

Functional enrichment analysis

To reveal the potential function of overlapping DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using R packages “clusterProfiler” and “org.Hs.eg.db”, and results were visualized using “Goplot” package.

Construction and evaluation of a prognostic risk model

Firstly, the least absolute shrinkage and selection operator (Lasso) analysis was performed on the overlapping DEGs to screen candidate genes related to prognosis. Then Cox regression analysis was conducted, and a prognostic risk model was established based on multivariate Cox regression analysis. The risk score for each patient in the TCGA-THCA cohort was calculated using expression levels and regression coefficients of risk genes in the multivariate Cox regression analysis:

Riskscore=∑1iβi×expi[2]

Samples were divided into high- and low-risk groups based on the optimal truncation value of the risk score. The optimal truncation value was calculated using “surv_cutpoint” in R package survival. Kaplan-Meier (K-M) survival curves were generated using the survival package to visualize the overall survival (OS) between the two risk groups. “pROC” package was employed to create receiver operating characteristic (ROC) curves, and the area under the curve (AUC) values for survival predictions were calculated. The “ggDCA” package was used to perform decision curve analysis (DCA). The GSE33630 dataset served as the validation set to assess the model’s predictive performance. Additionally, the difference in clinical factors such as age, gender, and pathological stage between the two risk groups was analyzed. Then, we conducted univariate and multivariable Cox regression analyses using the R package “survival” within TCGA cohorts to determine whether the risk score and clinical factors function as independent prognostic indicators.

Clinical analysis of the model genes

First, K-M curves were utilized to compare OS between groups with low and high expression levels of these genes. The expression levels of the model genes in tumor versus normal samples were analyzed in the TCGA-THCA cohort and subsequently validated using the GSE33630 dataset.

Gene set enrichment analysis (GSEA) and drug sensitivity analysis for model genes

GSEA was conducted to screen model gene-related pathways, and the results were visualized using ggplot2. To identify model gene-related drugs, the half maximal inhibitory concentration (IC50) of potential drugs for THCA was obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database. The oncoPredict package was used to calculate the drug sensitivity score for each patient. Then, Pearson correlation analysis was performed to determine the relationship between each gene and the drugs, with thresholds set at “correlation >0.4 and P<0.05”. The intersection of these results yielded 10 drugs, from which four commonly used drugs were selected for visualization.

Immune cell infiltration analysis

Firstly, the immune scores of 64 types of immune cells were calculated for patients in the TCGA-THCA cohort using the “IOBR” package and the xCell method. Based on the median expression value of the model genes, all patients were divided into high-expression and low-expression groups. Immune cells that showed significant differences in the high- or low-expression group of model genes were identified. After taking the intersection, 23 immune cell types were found to exhibit differences across high- and low-expression groups of all model genes. The correlation between these immune cells and the model genes was then calculated and displayed using a lollipop chart.

Single-cell RNA sequencing (scRNA-seq) analysis

Preprocessing and filtering of scRNA-seq were conducted using the Seurat package. Quality control criteria were set as nFeature_RNA >500, 1,000< nCount_RNA <20,000, and percent.mt <20. After uniformizing using the Scale Data function, the data was processed to find significant principal components through principal component analysis (PCA). Then, t-distributed stochastic neighbor embedding (tSNE) analysis was performed to screen out cell clusters. These cell clusters were annotated using SingleR 2.0.0 in R. The distribution of different cell types in tumor and normal samples was explored through Fisher’s test. The Wilcoxon-Mann-Whitney test was used to calculate the expression difference of each gene in the model in different samples.

Statistical analysis

All statistical analyses were performed using R software version 4.1.2. For survival analysis, survival differences between groups were evaluated by the log-rank test. Differences between the two groups were assessed using the Wilcoxon-Mann-Whitney test. P<0.05 was considered statistically significant.

Results

Functions of TDS-related DEGs in THCA

Based on the TCGA-THCA cohort, 2,448 DEGs in tumor and normal (Figure 1A) and 1,445 DEGs in high TDS and low TDS were identified (Figure 1B). Venn diagram revealed 990 overlapping genes between THCA-DEGs and TDS-DEGs (Figure 1C). Subsequently, GO and KEGG enrichment analyses were performed to determine the functions of these 990 overlapping genes. These 990 genes may function in extracellular region and participate in regulation of response to stimulus, regulation of multicellular organismal process, regulation of developmental process, response to external stimulus, and biological adhesion (Figure 1D). KEGG enrichment analysis showed that these genes were mainly enriched in pathways such as cell adhesion, apoptosis, and thyroid hormone synthesis (Figure 1E).

Figure 1 Functions of TDS-related DEGs in THCA. (A) DEGs between tumor and normal samples in TCGA-THCA cohort. (B) DEGs between high TDS and low TDS groups in TCGA-THCA cohort. (C) Venn diagram revealed overlapping genes between THCA-DEGs and TDS-DEGs. (D) GO enrichment analysis of 990 overlapping genes. (E) KEGG enrichment analysis of 990 overlapping genes. DEGs, differentially expressed genes; ECM, extracellular matrix; FDR, false discovery rate; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; TCGA, The Cancer Genome Atlas; TDS, thyroid differentiation score; THCA, thyroid cancer.

Construction of TDS-related prognostic risk model for THCA

To construct a risk model for predicting THCA prognosis, Lasso regression analysis was performed on the overlapping DEGs. The model was most accurate when the variable number was nine (Figure 2). The nine genes were then subjected to univariate Cox regression, and seven genes significantly related to prognosis were identified (Table 1). Then, multivariate Cox regression analysis was performed on the seven genes, and genes with P<0.05 were selected for model construction (Table 2). Finally, a four-gene prognostic model, including ATPase secretory pathway Ca2⁺ transporting 2 (ATP2C2), mast cell expressed membrane protein 1 (MCEMP1), FAM111 trypsin-like peptidase B (FAM111B), and uronyl 2-sulfotransferase (UST), was established. The risk score for this model was calculated as follows:

Riskscore=1.656∗ATP2C2+0.965∗MCEMP1−1.763∗FAM111B−1.146∗UST[3]

Figure 2 Lasso regression analysis based on 990 overlapping DEGs. DEGs, differentially expressed genes; Lasso, least absolute shrinkage and selection operator.

Evaluation of risk model’s performance and clinical relevance

Furthermore, TCGA-THCA patients were divided into high-risk or low-risk groups based on the optimal truncation value of risk scores. K-M curve revealed that patients in the high-risk group had lower OS than those in the low-risk group (Figure 3A). The AUC value of the risk model for predicting OS was 0.895 (Figure 3B). In an external dataset, the AUC value was 0.748 (Figure 3C). Additionally, DCA demonstrated good performance when the high-risk threshold exceeded 0.6 (Figure 3D). These results suggest that the risk model has good predictive performance for THCA prognosis.

Figure 3 Evaluation of risk model’s performance and clinical relevance. (A) Kaplan-Meier curve for TCGA-THCA patients. (B) ROC analysis of the risk model to predict the overall survival of THCA patients based on the TCGA cohort. (C) ROC analysis of the risk model to predict the overall survival of THCA patients based on the GSE33630 dataset. (D) Decision curve analysis for the model. (E) Kaplan-Meier curve for TCGA-THCA patients after adjusting for clinical features. (F) Risk of THCA patients with different ages, genders, and pathological stages. ns, no significance; **, P<0.01; ****, P<0.0001. AUC, area under the curve; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; THCA, thyroid cancer; yr, year.

Furthermore, as shown in the K-M curve, patients in the high-risk group still had significantly poorer prognoses after adjusting for the clinical features (Figure 3E). Therefore, we compared risk differences among THCA patients with different clinical characteristics. As shown in Figure 3F, individuals over 55 years old or in advanced stages exhibited a significantly higher risk of death (P<0.01), while no significant difference in risk was observed between genders. Univariate logistic regression analysis further revealed that age, pathological stage, and risk score were correlated with THCA prognosis (P<0.05, Table 3). Among these factors, age and the risk score were independent predictors of OS in THCA patients (P<0.001, Table 3).

Prognostic significance of four model genes in THCA

Additionally, a significant difference in OS was observed between the two patient groups, differentiated by the four genes’ low or high expression levels. As illustrated in Figure 4A-4D, patients with high expression of ATP2C2 and MCEMP1 had shorter OS, and those with high expression of FAM111B and UST had longer OS (P<0.01). To further investigate the relationship between these risk genes and THCA, we analyzed their expression in tumor versus normal samples. In the TCGA-THCA cohort, expression levels of ATP2C2 and UST were reduced whereas expression levels of MCEMP1 and FAM111B were elevated in the tumor group compared to the normal group (P<0.0001, Figure 4E). The four risk genes’ expression levels were then validated in the GSE33630 dataset and showed consistent expression with those in the TCGA-THCA cohort (P<0.05, Figure 4E). We then explored the association of hub genes and clinical features with OS of TCGA-THCA patients. As shown in Table 4, age, stage, ATP2C2, MCEMP1, FAM111B, and UST were all significantly associated with OS of THCA patients, however, only age, ATP2C2, MCEMP1, FAM111B, and UST could independently predict the OS. The finding indicated that these four risk genes could be biomarkers for THCA prognosis.

Figure 4 Prognostic significance of four risk genes in THCA. (A-D) Different overall survival between low and high ATP2C2, MCEMP1, FAM111B, and UST expression was represented in the Kaplan-Meier curves. (E) Expression levels of four risk genes between normal and tumor groups based on the TCGA-THCA and GSE33630 datasets. *, P<0.05; ****, P<0.0001. ATP2C2, ATPase secretory pathway Ca²⁺ transporting 2; FAM111B, FAM111 trypsin-like peptidase B; MCEMP1, mast cell expressed membrane protein 1; N, normal; T, tumor; TCGA, The Cancer Genome Atlas; THCA, thyroid cancer; UST, uronyl 2-sulfotransferase.

Potential pathways and drugs related to four risk genes in THCA

To reveal significant functional enrichment within four risk genes, GSEA was performed. According to KEGG pathway analysis, ATP2C2 was negatively related to apoptosis, focal adhesion, leukocyte transendothelial migration, and p53 signaling pathway (Figure 5A). MCEMP1 was positively associated with apoptosis, autoimmune thyroid disease, p53 signaling pathway, and Toll-like receptor signaling pathway (Figure 5B). Likewise, FAM111B was positively correlated to apoptosis, autoimmune thyroid disease, leukocyte transendothelial migration, and T cell receptor signaling pathway (Figure 5C). In addition, UST was negatively linked to apoptosis, autoimmune thyroid disease, focal adhesion, and p53 signaling pathway (Figure 5D). Notably, apoptosis and p53 signaling pathways were associated with more than three risk genes, suggesting that apoptosis and p53 signaling pathways may play critical roles in how these risk genes influence THCA progression.

Figure 5 Pathways related to four risk genes in THCA. (A-D) Gene set enrichment analysis revealed the pathways related to ATP2C2, MCEMP1, FAM111B, and UST, respectively. ATP2C2, ATPase secretory pathway Ca²⁺ transporting 2; KEGG, Kyoto Encyclopedia of Genes and Genomes; MCEMP1, mast cell expressed membrane protein 1; FAM111B, FAM111 trypsin-like peptidase B; UST, uronyl 2-sulfotransferase; THCA, thyroid cancer.

Moreover, using the GDSC database, we conducted a drug sensitivity analysis to identify potential drugs targeting the four risk genes in THCA. By applying criteria of correlation >0.4 and P<0.05, we identified 10 common drugs linked to all these genes, including dasatinib, Z-LLNle-CHO, EphB4, HG-5-88-01, AZD5438, PARP, avagacestat, paclitaxel, SU11274, and amuvatinib. The relationships between the risk genes and four commonly used anti-tumor drugs—dasatinib, avagacestat, paclitaxel, and amuvatinib—are illustrated in Figure 6. Notably, the drug sensitivity to dasatinib, avagacestat, and paclitaxel was significantly higher in the ATP2C2 high-expression group compared to the low-expression group (P<0.01, Figure 6A). A similar increase in sensitivity to dasatinib, avagacestat, and paclitaxel was observed in the UST high-expression group (P<0.05, Figure 6D). Amuvatinib’s sensitivity was elevated in both the MCEMP1 and FAM111B high-expression groups (P<0.001, Figure 6B,6C).

Figure 6 Drug sensitivity analysis. (A-D) Drug sensitivity of dasatinib, avagacestat, paclitaxel, and amuvatinib in low and high ATP2C2, MCEMP1, FAM111B, and UST expression. *, P<0.05; **, P<0.01; ****, P<0.0001. ATP2C2, ATPase secretory pathway Ca²⁺ transporting 2; FAM111B, FAM111 trypsin-like peptidase B; MCEMP1, mast cell expressed membrane protein 1; UST, uronyl 2-sulfotransferase.

Relationship of four risk genes with immune cells in THCA

A previous study demonstrated that there is a correlation between clinical outcomes in THCA patients and the density and type of tumor-infiltrating immune cells (15). Therefore, we further investigated immune cells associated with four risk genes. The results are illustrated in Figure 7. Some immune cells were significantly related to four risk genes, such as macrophages, macrophages M1, regulatory T cells (Treg), monocytes, class-switched memory B cells, natural killer T cells (NKT), dendritic cells (DC), and immature DC (iDC).

Figure 7 Relationship of four risk genes with immune cells in THCA. (A-D) Correlation between immune cells and four risk genes ATP2C2, MCEMP1, FAM111B, and UST, respectively. ATP2C2, ATPase secretory pathway Ca²⁺ transporting 2; Cor, correlation; DC, dendritic cells; FAM111B, FAM111 trypsin-like peptidase B; HSC, hematopoietic stem cell; iDC, immature dendritic cells; MCEMP1, mast cell expressed membrane protein 1; MPP, multipotent progenitor; NKT, natural killer T cells; THCA, thyroid cancer; tSNE, t-distributed stochastic neighbor embedding; UST, uronyl 2-sulfotransferase.

Four risk genes’ distribution in immune cells of THCA patients

Subsequently, single-cell sequencing analysis was conducted to explore the specific cells where the risk genes might function. In the GSE191288 dataset, 19 cell clusters were identified and annotated into eight cell types, including Treg, natural killer (NK) cells, B cells, endothelial cells, epithelial cells, CD4⁺ memory cells, monocytes, and fibroblasts (Figure 8A,8B). It was observed that Treg, NK cells, and CD4⁺ memory cells were mainly enriched in THCA tumor tissues compared to non-tumor (NT) tissues (Figure 8C). We then analyzed the distribution of four risk genes in these eight cell types. As shown in Figure 9A, UST was found in most of these cells. The bubble plot further revealed that UST was mostly expressed in the CD4⁺ memory cells, followed by epithelial cells, NK cells, and fibroblasts (Figure 9B). These results indicated that UST may play critical functions in THCA, therefore, it was selected for further analysis.

Figure 8 Identification of cell clusters based on single-cell RNA sequencing data. (A) All cells were clustered into 19 clusters. (B) A total of 8 cell types were identified. (C) Cell distribution between NT tissues and THCA tumor tissues. NK, natural killer; NT, non-tumor; THCA, thyroid cancer; tSNE, t-distributed stochastic neighbor embedding.

Figure 9 Association of four risk genes with 8 cell types. (A) Distribution of four risk genes in 8 cell types. (B) Expression of four genes in 8 cell types. ATP2C2, ATPase secretory pathway Ca²⁺ transporting 2; FAM111B, FAM111 trypsin-like peptidase B; MCEMP1, mast cell expressed membrane protein 1; NK, natural killer; tSNE, t-distributed stochastic neighbor embedding; UST, uronyl 2-sulfotransferase.

UST was associated with the p53 signaling pathway in CD4⁺ memory cells

GSEA showed the enrichment of apoptosis and p53 signaling pathways, and UST was negatively linked to these two pathways. Figure 10A exhibited the pathway score of apoptosis and the p53 signaling pathway in eight cell types. The pathway scores of these two pathways were especially high in T cells, NK cells, and CD4⁺ memory cells. Subsequently, the association among UST, two pathways, and eight cell types was explored. As shown in Figure 10B, UST was significantly related to the p53 signaling pathway in CD4⁺ memory cells. These findings suggested that UST may play a significant role in modulating the p53 signaling pathway, particularly within CD4⁺ memory cells, indicating its potential influence on THCA.

Figure 10 UST was associated with the p53 signaling pathway in CD4⁺ memory cells. (A) Pathway scores of apoptosis and p53 signaling pathways in cells. (B) Association of UST with apoptosis and p53 signaling pathway in 8 cell types. NK, natural killer; tSNE, t-distributed stochastic neighbor embedding; UST, uronyl 2-sulfotransferase.

Discussion

The current study presented a comprehensive analysis of TDS-related DEGs in THCA, leading to the construction of a robust prognostic risk model (comprising ATP2C2, MCEMP1, FAM111B, and UST). The analysis not only highlighted the potential functional roles of these DEGs but also delved into their clinical significance, including their impact on prognosis and potential as therapeutic targets. Notably, among the identified genes, UST emerged as a critical player in the modulation of the p53 signaling pathway, particularly in CD4⁺ memory cells, which might have profound implications for THCA progression and treatment.

Research suggests that many ATC result from the dedifferentiation of DTC, and some genes may significantly influence the occurrence and progression of dedifferentiated DTCs through multiple pathways (14,16). We identified 990 overlapping DEGs between THCA-related genes and TDS-related genes. GO and KEGG enrichment analyses revealed that these DEGs were involved in key processes such as cell adhesion, apoptosis, and thyroid hormone synthesis. Impaired cell adhesion has been associated with enhanced metastatic potential in THCA (17). The enrichment of these DEGs in pathways related to apoptosis and thyroid hormone synthesis also aligns with the known pathophysiology of THCA, where dysregulated apoptosis contributes to uncontrolled cell proliferation (18,19) and alterations in thyroid hormone pathways are closely linked to THCA (20). The functional analysis of these overlapping DEGs thus reinforces their potential role in driving THCA progression and highlights the importance of targeting these pathways in therapeutic strategies.

Furthermore, based on the overlapping DEGs, a TDS-related prognostic risk model was constructed and resulted in the identification of four key genes—ATP2C2, MCEMP1, FAM111B, and UST—that were significantly associated with THCA prognosis. A recent study demonstrated that TDS quantified the survival outcomes of THCA and could serve as a potential prognostic biomarker in THCA (12). Likewise, the risk model in our study demonstrated an AUC value of 0.895. Additionally, the inclusion of ATP2C2 and MCEMP1 as risk factors and FAM111B and UST as protective factors is particularly noteworthy. ATP2C2, involved in calcium ion transport, has been implicated in cancer progression through its role in calcium signaling, which is essential for various cellular processes, including proliferation, migration, and invasion (21,22). ATP2C2 has been identified as an independent prognostic factor in THCA (23). Similarly, MCEMP1, which is associated with immune response modulation, has been linked to poor prognosis in several cancers (24-26). In contrast, FAM111B, a gene involved in DNA replication and repair, has been identified as an oncogene in bladder cancer and hepatocellular carcinoma (27,28). However, in THCA, FAM111B has been reported to inhibit migration, invasion, and glycolysis both in vitro and in vivo (29). In this study, the prognostic significance of these genes, particularly their varying expression levels between tumor and normal tissues, suggests that they could serve as valuable biomarkers for THCA prognosis.

Among the identified genes, UST stands out due to its association with the p53 signaling pathway and its expression in CD4⁺ memory cells. Previous studies have not extensively explored UST’s role in THCA or its involvement in the p53 signaling pathway. However, UST (also named 2OST), an enzyme involved in the modification of proteoglycans, has been implicated in other cancers, where its dysregulation has been associated with altered cell signaling and tumor growth (30,31). In our study, THCA patients with high expression of UST had a better prognosis. Therefore, we speculate that UST may be a tumor suppressor gene in THCA. The p53 signaling pathway is frequently dysregulated in various cancers (32). In THCA, mutations or alterations in the p53 pathway have been linked to poor prognosis and resistance to treatment (33,34). Our study revealed that UST was negatively associated with the p53 signaling pathway in THCA, particularly within CD4⁺ memory cells. CD4⁺ memory cells are a subset of T cells that play a crucial role in the immune response by maintaining long-term immunity and facilitating rapid responses upon re-exposure to antigens (35). In the patients with undifferentiated THCA, it is reported that CD4⁺ memory cells were significantly decreased compared with the DTC patients (36). Another study demonstrated that aggressive THCA at presentation and during follow-up is linked to the reduction of T cells, especially the CD4⁺ T cells (37). The expression of UST in CD4⁺ memory cells suggests that UST may influence the function of CD4⁺ memory cells, thereby impacting the immune response against THCA. We propose that in THCA, UST may inhibit the p53 signaling pathway, thus activating CD4⁺ memory cells and enhancing their ability to carry out immune surveillance. The p53 pathway is known to regulate various aspects of the immune response, including function of T cell (38). Watanabe et al. reported that p53 downregulation by T cell receptor signaling is critical for antigen-specific CD4⁺ T cell responses (39). Another study demonstrated that interferon regulatory factor eight maintained tumor cell sensitivity to cytotoxic T cells by inhibiting p53 expression (40). These studies further support our speculation.

There are several limitations in this study. Firstly, the reliance on public datasets may introduce selection bias, limiting the generalizability of the results. Secondly, the findings are primarily based on bioinformatic analyses, and experimental validation of the identified genes and pathways is necessary to confirm their biological relevance and therapeutic potential.

Conclusions

In conclusion, this study provides valuable insights into the roles of TDS-related DEGs in THCA and highlights the potential of the prognostic model in predicting patient outcomes. Additionally, our study identifies UST as a critical regulator of the p53 signaling pathway in CD4⁺ memory cells, with significant implications for THCA progression and prognosis. These findings offer promising avenues for developing personalized therapies and improving the clinical management of THCA.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-460/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-460/prf

Funding: This work was supported by the 2022 Wenzhou Basic Research Project (No. Y20220866).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-460/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The Ethics Committee of The Second Affiliated Hospital and Yuying Children’s Hospital of Wenzhou Medical University deemed that this research is based on open-source data, so the need for ethics approval was waived.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Ju SH, Song M, Lim JY, et al. Metabolic Reprogramming in Thyroid Cancer. Endocrinol Metab (Seoul) 2024;39:425-44. [Crossref] [PubMed]
2. Zhang L, Feng Q, Wang J, et al. Molecular basis and targeted therapy in thyroid cancer: Progress and opportunities. Biochim Biophys Acta Rev Cancer 2023;1878:188928. [Crossref] [PubMed]
3. Boucai L, Zafereo M, Cabanillas ME. Thyroid Cancer: A Review. JAMA 2024;331:425-35. [Crossref] [PubMed]
4. Han X, Liao J, Zhou Y, et al. Thyroid cancer prognostic biomarker ARL4A and its relationship with immune infiltration. Int J Clin Exp Pathol 2024;17:108-20. [Crossref] [PubMed]
5. Sun D, Zhang X, Jin X, et al. BRAF(V600E) mutation is associated with better prognoses in radioactive iodine refractory thyroid cancer patients treated with multi-kinase inhibitors: a retrospective analysis of registered clinical trials. Thyroid Res 2025;18:5. [Crossref] [PubMed]
6. Yang J, Chen Y, Zhang S, et al. Clinical significance of RETN gene expression and rs3219175 G > a polymorphism in cancer. Nucleosides Nucleotides Nucleic Acids 2024; Epub ahead of print. [Crossref]
7. Wei M, Wang R, Zhang W, et al. Landscape of gene mutation in Chinese thyroid cancer patients: Construction and validation of lymph node metastasis prediction model based on clinical features and gene mutation marker. Cancer Med 2023;12:12929-42. [Crossref] [PubMed]
8. Leeman-Neill RJ, Brenner AV, Little MP, et al. RET/PTC and PAX8/PPARγ chromosomal rearrangements in post-Chernobyl thyroid cancer and their association with iodine-131 radiation dose and other characteristics. Cancer 2013;119:1792-9. [Crossref] [PubMed]
9. Zhang W, Lin S, Wang Z, et al. Coexisting RET/PTC and TERT Promoter Mutation Predict Poor Prognosis but Effective RET and MEK Targeting in Thyroid Cancer. J Clin Endocrinol Metab 2024;109:3166-75. [Crossref] [PubMed]
10. Ruchong P, Haiping T, Xiang W. A Five-Gene Prognostic Nomogram Predicting Disease-Free Survival of Differentiated Thyroid Cancer. Dis Markers 2021;2021:5510780. [Crossref] [PubMed]
11. Cancer Genome Atlas Research Network. Integrated genomic characterization of papillary thyroid carcinoma. Cell 2014;159:676-90. [Crossref] [PubMed]
12. Wang JR, Zafereo ME, Cabanillas ME, et al. The Association Between Thyroid Differentiation Score and Survival Outcomes in Papillary Thyroid Carcinoma. J Clin Endocrinol Metab 2025;110:356-63. [Crossref] [PubMed]
13. Boucai L, Seshan V, Williams M, et al. Characterization of Subtypes of BRAF-Mutant Papillary Thyroid Cancer Defined by Their Thyroid Differentiation Score. J Clin Endocrinol Metab 2022;107:1030-9. [Crossref] [PubMed]
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. [Crossref] [PubMed]
15. Zhang N, Zhang H, Wu W, et al. Machine learning-based identification of tumor-infiltrating immune cell-associated lncRNAs for improving outcomes and immunotherapy responses in patients with low-grade glioma. Theranostics 2022;12:5931-48. [Crossref] [PubMed]
16. Zhu X, Hu C, Zhang Z, et al. PD-L1 and B7-H3 are Effective Prognostic Factors and Potential Therapeutic Targets for High-Risk Thyroid Cancer. Endocr Pathol 2024;35:230-44. [Crossref] [PubMed]
17. Pinto D, Parameswaran R. Role of Truncated O-GalNAc Glycans in Cancer Progression and Metastasis in Endocrine Cancers. Cancers (Basel) 2023;15:3266. [Crossref] [PubMed]
18. Li C, Tang Y, Li Q, et al. The prognostic and immune significance of C15orf48 in pan-cancer and its relationship with proliferation and apoptosis of thyroid carcinoma. Front Immunol 2023;14:1131870. [Crossref] [PubMed]
19. Lu K, Wei W, Hu J, et al. Apoptosis Activation in Thyroid Cancer Cells by Jatrorrhizine-Platinum(II) Complex via Downregulation of PI3K/AKT/Mammalian Target of Rapamycin (mTOR) Pathway. Med Sci Monit 2020;26:e922518. [Crossref] [PubMed]
20. Lukasiewicz M, Zwara A, Kowalski J, et al. The Role of Lipid Metabolism Disorders in the Development of Thyroid Cancer. Int J Mol Sci 2024;25:7129. [Crossref] [PubMed]
21. Dang DK, Makena MR, Llongueras JP, et al. A Ca(2+)-ATPase Regulates E-cadherin Biogenesis and Epithelial-Mesenchymal Transition in Breast Cancer Cells. Mol Cancer Res 2019;17:1735-47. [Crossref] [PubMed]
22. Zhao M, Zhang Q, Song Z, et al. ATP2C2 as a novel immune-related marker that defines the tumor microenvironment in triple-negative breast cancer. Transl Cancer Res 2023;12:1802-15. [Crossref] [PubMed]
23. Zhao H, Zhang S, Shao S, et al. Identification of a Prognostic 3-Gene Risk Prediction Model for Thyroid Cancer. Front Endocrinol (Lausanne) 2020;11:510. [Crossref] [PubMed]
24. Wang D, Gu Y, Huo C, et al. MCEMP1 is a potential therapeutic biomarker associated with immune infiltration in advanced gastric cancer microenvironment. Gene 2022;840:146760. [Crossref] [PubMed]
25. Ling L, Hu T, Zhou C, et al. Low Expression MCEMP1 Promotes Lung Adenocarcinoma Progression and its Clinical Value. Curr Cancer Drug Targets 2025;25:281-93. [Crossref] [PubMed]
26. Huang P, Liu Y, Jia B. The Expression, Prognostic Value, and Immunological Correlation of MCEMP1 and its Potential Role in Gastric Cancer. J Oncol 2022;2022:8167496. [Crossref] [PubMed]
27. Huang N, Peng L, Yang J, et al. FAM111B Acts as an Oncogene in Bladder Cancer. Cancers (Basel) 2023;15:5122. [Crossref] [PubMed]
28. Li F, He HY, Fan ZH, et al. Silencing of FAM111B inhibited proliferation, migration and invasion of hepatoma cells through activating p53 pathway. Dig Liver Dis 2023;55:1679-89. [Crossref] [PubMed]
29. Zhu X, Xue C, Kang X, et al. DNMT3B-mediated FAM111B methylation promotes papillary thyroid tumor glycolysis, growth and metastasis. Int J Biol Sci 2022;18:4372-87. [Crossref] [PubMed]
30. Liu C, Sheng J, Krahn JM, et al. Molecular mechanism of substrate specificity for heparan sulfate 2-O-sulfotransferase. J Biol Chem 2014;289:13407-18. [Crossref] [PubMed]
31. Ferguson BW, Datta S. Role of heparan sulfate 2-o-sulfotransferase in prostate cancer cell proliferation, invasion, and growth factor signaling. Prostate Cancer 2011;2011:893208. [Crossref] [PubMed]
32. Hernández Borrero LJ, El-Deiry WS. Tumor suppressor p53: Biology, signaling pathways, and therapeutic targeting. Biochim Biophys Acta Rev Cancer 2021;1876:188556. [Crossref] [PubMed]
33. Nappi A, Miro C, Pezone A, et al. Loss of p53 activates thyroid hormone via type 2 deiodinase and enhances DNA damage. Nat Commun 2023;14:1244. [Crossref] [PubMed]
34. Liu JF, Zou B, Xiang C, et al. Comprehensive bioinformatics analysis unveils THEMIS2 as a carcinogenic indicator related to immune infiltration and prognosis of thyroid cancer. Sci Rep 2024;14:8156. [Crossref] [PubMed]
35. Künzli M, Masopust D. CD4(+) T cell memory. Nat Immunol 2023;24:903-14. [Crossref] [PubMed]
36. Lee SE, Koo BS, Sun P, et al. Neutrophil diversity is associated with T-cell immunity and clinical relevance in patients with thyroid cancer. Cell Death Discov 2024;10:222. [Crossref] [PubMed]
37. Kotwal A, Gustafson MP, Bornschlegl S, et al. Circulating immunophenotypes are potentially prognostic in follicular cell-derived thyroid cancer. Front Immunol 2023;14:1325343. [Crossref] [PubMed]
38. Brummer T, Zeiser R. The role of the MDM2/p53 axis in antitumor immune responses. Blood 2024;143:2701-9. [Crossref] [PubMed]
39. Watanabe M, Moon KD, Vacchio MS, et al. Downmodulation of tumor suppressor p53 by T cell receptor signaling is critical for antigen-specific CD4(+) T cell responses. Immunity 2014;40:681-91. [Crossref] [PubMed]
40. Poschel DB, Kehinde-Ige M, Klement JD, et al. IRF8 Regulates Intrinsic Ferroptosis through Repressing p53 Expression to Maintain Tumor Cell Sensitivity to Cytotoxic T Lymphocytes. Cells 2023;12:310. [Crossref] [PubMed]