Construction of an anaplastic thyroid cancer stratification signature to guide immune therapy selection and validation of the pivotal gene HLF through in vitro experiments

Article type: Research Article

Keywords: anaplastic thyroid cancer (ATC), T cell immunity, machine learning, prediction, model

Authors: Li Pengping, Yin Kexin, Xie Yuwei, Sun Ke, Li Rongguo, Wang Zhenyu, Jin Haigang, Wang Shaowen, Huang Yuqing

Affiliations: Department of Thyroid & Breast Surgery, The First People’s Hospital of Xiaoshan District, Xiaoshan Affiliated Hospital of Wenzhou Medical University, Hangzhou, Zhejiang, China; The First Affiliated Hospital of Anhui Medical University, Hefei, Anhui, China; Neuromedicine Center, The University of Hong Kong-Shenzhen Hospital, Shenzhen, Guangdong, China

License: Copyright © 2025 Pengping, Kexin, Yuwei, Ke, Rongguo, Zhenyu, Haigang, Shaowen and Yuqing CC BY 4.0 This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Article links: DOI: 10.3389/fimmu.2024.1478904  |  PubMed: 39872516  |  PMC: PMC11769983

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (18.4 MB)

Introduction

Although the prognosis for most thyroid cancer patients is favorable, anaplastic thyroid carcinoma (ATC) remains a notably aggressive form with a median survival time of only five months (ref. 1, ref. 2). ATC comprises only 2% of all types of thyroid cancers but accounts for a disproportionate 14% to 39% of deaths associated with this cancer (ref. 3). For patients who outlive the median survival time, conventional therapies such as chemotherapy and radiotherapy offer limited benefits (ref. 4, ref. 5). Besides, early diagnosis of ATC presents another significant challenge. Thus, detecting ATC in its initial stages and finding other effective therapies are both crucial for improving patient outcomes.

It has been reported that resistance to chemotherapy and radiotherapy is the main reason for poor survival. For some patients carrying the BRAF and MEK gene mutation, the combination of dabrafenib and trametinib has shown promise in extending the effective treatment period (ref. 6). For those patients without these mutations, immunotherapy (anti-PD-1 and anti-PD-L1) has manifested a 1-year survival rate of approximately 40%, making it the most promising known therapy (ref. 7). However, more research is still urgently needed to identify the subgroup of patients who can benefit from the immunotherapy.

Our study endeavors to identify ATC patients who could potentially benefit from immunotherapy. Leveraging a set of five uniquely expressed genes across various types of thyroid cancer, our research has developed a machine-learning model capable of distinguishing each type within the GEO dataset of thyroid cancer patients (GSE60542, GSE76039, GSE33630, GSE53157, GSE65144, GSE29265, GSE82208, GSE27155, GSE58545, GSE54958, and GSE32662). Utilizing these genes, we have stratified ATC into three distinct groups, each demonstrating significantly different responses to anti-PD1 therapy. Additionally, we employed weighted gene co-expression network analysis (WGCNA) to identify 12 differentially expressed genes intimately associated with both the grouping and immunotherapy outcomes. This led to the creation of a refined signature that could more accurately predict ATC’s immune responsiveness to anti-PD1 therapy, which was further corroborated using thyroid cancer cohorts and 9 melanoma cohorts from the clinical trial. Among the 12 genes analyzed, HLF emerged as significantly associated with various cancer hallmarks. Our study elucidated the mechanism by which HLF impeded anaplastic thyroid carcinoma (ATC) progression. Specifically, HLF down-regulated the epithelial-to-mesenchymal transition (EMT) pathway, reduced T cell exhaustion, and increased sensitivity to sorafenib, as demonstrated by our in-vitro experiments.

Methods

Data collection

Gene expression profiles and clinical characteristics of thyroid cancer were collected from the cancer genome atlas (TCGA, https://portal.gdc.cancer.gov) and gene expression omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). A total of 506 samples with follow-up data were collected from TCGA, and 715 samples were collected from GEO datasets GSE27155, GSE58545, GSE76039, GSE82208, GSE60542, GSE29265, GSE33630, GSE65144, GSE53157, GSE54958, and GSE32662. Among these, 114 samples were thyroid non-cancerous tissue (TNC), 78 samples were anaplastic thyroid carcinoma (ATC), 225 samples were papillary thyroid carcinoma (PTC), 40 samples were follicular thyroid carcinoma (FTC), and 54 samples were medullary thyroid carcinoma (MTC). Verified cohorts of clinical trials (anti-PD1) were collected from GSE115821 (N=37), GSE78220 (N=28), GSE91061 (N=109), Nathanson_2017 (N=24), phs000452 (N=153), and PRJEB23709 (N=91).

Bioinformatic analysis

Subgrouping signature construction

GEO cohorts of thyroid cancer (TC) were used to identify differentially expressed genes (DEGs) in ATC using R4.2.0 (package: DESeq2). The ConsensusClusterPlus package in R4.2.0 was employed to perform consensus clustering analysis based on the 5 DEGs (BCL2, BHLHE40, MICAL2, TGM2, TPO) (parameters: maxK=10, reps=50). The consensus cumulative density function (CDF) and delta area indicated that a 3-subgroup division was the optimal outcome.

Subtype identification model construction

Decision curve analysis (DCA) was used to evaluate the value of the 5-DEGs in identifying ATC or other subgroups of TC using the ggDCA package in R4.2.0. The DALEX, randomForest, kernlab, xgboost, caret, and pROC packages were applied to construct identification models.

Tumor immune index calculation

Infiltration immune cell fractions were predicted using CIBERSORT, ssGSEA (single sample GSEA), and Pompimol Charoentong’s algorithm in R4.2.0. The immune score was predicted using the estimate package in R4.0. Tumor Immune Dysfunction and Exclusion (TIDE) and anti-PD1 response were predicted using the online tool (http://tide.dfci.harvard.edu). To assess MeTIL characteristics, the individual methylation values of MeTIL markers were converted to MeTIL scores using principal component analysis (PCA). The data were converted to a unit-free Z-Score by applying the formula (x-μ)/σ. According to the median value of PDCD1, the samples were divided into high and low-expression groups. Wilcoxon Rank Sum Tests were used to compare the MeTIL scores between the two groups (ref. 8).

TIP (Tracking Tumor Immunophenotype) was a meta-server that systematically integrates two existing third-party methods, “ssGSEA” and “CIBERSORT”, for tracking, analyzing, and visualizing the anti-cancer immune state and the proportion of immune-infiltrating cells in the seven steps of the cancer immune cycle using RNA-seq or microarray data. Spearman correlation between genes and TIP scores, as well as the autocorrelation between TIP scores, were calculated, and the linkET package was used for visualization (ref. 8). The red and green lines represented positive and negative correlations, respectively, while the gray lines indicated no significance. The thickness of the lines represented the absolute value of the correlation coefficient. The correlation in the triangular region was represented by the color depth and size of the square: red/blue indicated a positive/negative correlation, with darker colors signifying more significant P-values, and larger squares representing greater absolute values of the correlation coefficient. Easier was a tool for predicting biomarker-based immunotherapies (Cytolytic score, CYT; Tertiary lymphoid structures, TLS; IFNy, T cell-inflamed, Chemokines) based on cancer-specific immune response models, aiming to predict anti-tumor immune responses from RNA-seq data. The CYT level of TCGA-BLCA was calculated using the Easier package to evaluate CYT characteristics. According to the median value of PDCD1, the samples were divided into high and low-expression groups (ref. 9). Wilcoxon Rank Sum Tests were used to compare the CYT scores between the two groups.

Cell signaling score calculation

The CancerSEA database collated 14 different functional states of tumor cells (ref. 10). The Z-score algorithm, proposed by Lee et al. (ref. 11), integrated characteristic gene expression to reflect the activity of a given pathway. Fourteen functional state gene sets were calculated using the Z-score algorithm in the R-package GSVA. The values of each gene set were enumerated separately as Z-scores. Pearson correlations between genes and the Z-scores of each gene set were then calculated.

Pan-cancer analysis of HLF expression

The expression level of HLF in pan-cancer was calculated by R4.2.0 (TCGA cohorts).

The 12-gene signature construction

Firstly, WGCNA (Weighted Gene Co-Expression Network Analysis) was applied to screen out hub genes using the WGCNA package in R4.2.0. The most correlated gene sets (both negative and positive) were collected for subsequent machine learning. AI modeling for ATC subgroup stratification was developed using six AI functions: extreme gradient boosting (XGboost, xgboost package in R4.2.0), support vector machine (SVM, e1071 packages in R4.2.0), multi-logistic (nnet packages in R4.2.0), random forest (RF, randomForest package in R4.2.0), deep learning (DL, h2o package in R4.2.0), and K-Nearest Neighbor (KNN, kknn package in R4.2.0). During model construction, 75% of the data was randomly selected as the training cohort, and 25% was randomly selected as the testing cohort. Gene expression values were standardized to range from 0 to 1 using the preProcess function (caret and tidyverse packages).

Biological experiments

Clinical sample collection

Twenty-two thyroid cancer samples were collected from 2021-12-01 to 2022-08-01. All experiments were approved by the Medical Ethics Committee of The First People’s Hospital of Xiaoshan District, Xiaoshan Affiliated Hospital of Wenzhou Medical University. All patients with ATC were confirmed by at least two pathologists.

Multiple immune fluorescence staining

The procedures for paraffin embedding, tissue sectioning, and immunohistochemistry for HLF, CD8, and PD1 expression levels were performed as previously described (PMID: 23200678 and 20571492). The working concentrations of antibodies against HLF (Proteintech, Wuhan, China), CD8 (Abcam, Shanghai), and PD1 (Proteintech, Wuhan, China) were 1:150. The protein expression levels were assessed by Mean of Integrated Option Density (IOD) with Image-Pro Plus. Briefly, the area of interest (AOI) was detected to gain the Mean of IOD (IOD/AOI, MI).

Reagents

Sorafenib was purchased from CSNpharm (A316727) and dissolved in PBS. Antibodies against beta-actin (AF7018, Affinity), CD8 (GB15068, Servicebio), HLF (DF7892, Affinity), N-cadherin (AF5239, Affinity), E-cadherin (BF0219, Affinity), Vimentin (BF8008, Affinity), Twist1 (AF4009, Affinity), Snail1 (AF6032, Affinity), PD-L1 (BF8035, Affinity), phosphorylated-JAK3 (p-JAK3, AF8160, Affinity), JAK3 (AF0008, Affinity), p-STAT3 (AF3293, Affinity), and STAT3 (BF6294, Affinity) were used for western blot.

Cell culture

ATC cell lines (CAL62, TCO1) were obtained from the cell bank of the Chinese Academy of Science in 2022 with STR matching analysis. The culture media for both ATC cell lines were DMEM with 10% fetal calf serum and 100 units/mL penicillin and streptomycin.

Small interfering RNA experiments

5 × 10^5 ATC cells were transplanted into 6-well plates for 24 hours and then transfected with three different sequences of HLF siRNA (GenePharma, Shanghai, China) for 48, 72, and 96 hours using Lipofectamine 3000 reagent (Invitrogen, USA) and Opti-MEM (Life Technologies, USA), according to the manufacturer’s instructions for optimal transfection efficiency. The three siRNA sequences for HLF were as follows:

• Sequence-1
• Forward (5′-3′): TGCAAAATGTTCAAAATTGAA
• Reverse (5′-3′): CAATTTTGACATTTTGCTAA

• Sequence-2
• Forward (5′-3′): ATTAAAAAAAAACTTTTCGGGTC
• Reverse (5′-3′): CGAAAGTTTTTTTTTAATAT

• Sequence-3
• Forward (5′-3′): AAATGTTGCTGAGCTTTTCCT
• Reverse (5′-3′): GAAAGCTCAGCAACTTTTA

Western blot

Total protein extraction: Cells were harvested using a cytology brush and lysed with RIPA lysis buffer (Sigma, USA) supplemented with phosphorylase and protease inhibitor mixture (Thermo, USA), and quantified by the BCA assay. Cytoplasmic and nuclear protein extraction: Cells were harvested using trypsin (Invitrogen), then cytoplasmic and nuclear proteins were extracted using the Cytoplasmic and Nuclear Protein Extraction Kit (Thermo Scientific, USA) according to the protocol, and quantified by the BCA assay. Protein samples were separated by SDS-PAGE (EpiZyme, China, PG113) and transferred to a 0.45 μm or 0.20 μm pore-sized PVDF membrane (Millipore, USA). The membranes were blocked with Tris-buffered saline containing 5% skim milk powder (Biosharp, BS102-500 g, China) at room temperature for 1 hour, followed by incubation with primary antibodies at 4°C overnight. The next day, after three washes with TBST solution, the membranes were incubated with secondary antibodies at room temperature for 60 minutes. Finally, immunoreactive bands were detected using an enhanced chemiluminescence kit (Biosharp, BL520B, China). The conjugation yield was calculated via gel band quantification using Image J software (ref. 12).

Migration ability assays

For trans-well assays, 50,000 cells, with or without special treatments, were transplanted into trans-well plates (24-well, 8.0μm, Corning Incorporated, Corning, NY, USA) with a 10% gradient of fetal calf serum for 48 hours. After 24 hours of incubation, the cells that had migrated to the lower surface of the filter were fixed with 4% paraformaldehyde and stained with hematoxylin and eosin. The stained cells were then observed and photographed using a light microscope. Quantification of the passed cell area was performed using Image-Pro Plus (ref. 13).

Live & dead cell staining

Live and dead cell staining was carried out using Calcein AM/PI staining. After being seeded in a 24-well plate and cultured for 24 hours, ATC cells were treated with DMSO or 5μM GEM for another 48 hours. Then, all cells were co-cultured with Calcein AM and PI and observed at 480 nm and 525 nm, respectively.

PBMCs extraction

Simply, PBMCs were isolated via Ficoll-Paque density gradient centrifugation: 5 mL of peripheral blood was collected from healthy female volunteers, diluted with PBS at a 1:1 ratio, followed by gentle mixing. Add 10 mL of the diluted blood to 2 mL of Ficoll liquid (density 1.077). The clear stratification of blood and Ficoll liquid confirmed success. Carefully transferred the sample to the centrifuge and spin at 500 g for 15 minutes. Removed the centrifuge tube with care, aspirate the white thin film layer in the middle, representing individual nucleated cells. Wash the isolated nucleated cells with 10 mL of PBS, centrifuge at 250 g for 10 minutes, and discarded the supernatant. Repeat the washing step once and the suspended cells were frozen in vials at 100 million cells/mL in HI FBS with 5% DMSO after washing. Stored in liquid nitrogen, they were revived gradually and washed in pre-warmed RPMI with FBS and pen/strep. Following a 4-5 hour incubation at 37°C, viability was assessed using Trypan blue (0.1%).

Flow cytometry

The co-cultured PBMC were stained with Fixable Viability Stain (Thermo, L34965) and Fc receptor blocking reagent [Ultra-LEAF™ Purified anti-mouse CD16/32 (101320, BioLegend)]. Next, they were stained with CD-3 (BD 557943), PD-1 (BD 561273), and CD8 antibody (thermo, A15448). The prepared single-cell suspensions were filtered through 40-μm nylon meshes (352340, Corning). Results were then acquired using BD Calibur, BD Fortessa, or Miltenyi MACSQuant systems. Data were analyzed with FlowJo_V10 software (TreeStar).

Statistical analysis

All data analyses were performed in R4.2.0. Pearson’s test was used to calculate the correlation between different genes. Wilcox rank sum test and Kruskal-Wallis rank sum test were used to assess differences in continuous variables. Univariate Cox regression was performed to calculate the hazard ratio (HR), and the log-rank test was used to compare survival differences. Heatmaps were generated using the pheatmap package in R4.2.0. Receiver operating characteristic (ROC) curves and AUC values were generated using the pROC package in R4.2.0. GO and KEGG analyses were performed using the clusterProfiler package in R4.2.0. P<0.05 was considered to indicate a statistically significant difference.

Results

Five distinct genes were expressed significantly between any pair of thyroid cancer subtypes

The main available gene expression data of 537 thyroid cancer samples from GEO database (GSE60542, GSE76039, GSE33630, GSE53157, GSE65144, GSE29265, GSE82208, GSE27155, GSE58545, GSE54958, and GSE32662) were extracted for analysis of differentially expressed genes between different subtypes of thyroid cancer (f1). Multiple comparisons between each group identified five genes (BCL2, BHLHE40, MICAL2, TGM2, TPO) that are significantly differentially expressed (P value > 0.05 and a fold change (|FC|) > 2) in any pair of all the subtypes (thyroid carcinoma (TCA), differentiated thyroid carcinoma (DTC), anaplastic thyroid carcinoma (ATC), papillary thyroid carcinoma (PTC), medullary thyroid carcinoma (MTC), follicular thyroid carcinoma (FTC), thyroid non-cancerous tissues (TNC) (f1).

The model based on the 5 genes by machine learning can distinguish each subtype well

BCL2, BHLHE40, MICAL2, TGM2, and TPO were further employed to make a model to differentially diagnose TCA from TNC (f2), ATC from DTC (f2), MTC from DTC (f2), and FTC from DTC (f2). Different machine learning algorithms, for instance, Random Forest (RF), Support Vector Machine (SVM), eXtreme Gradient Boosting (XGB), Generalized Linear Model (GLM), Gradient Boosting Machine (GBM), Kernel k-Nearest Neighbors (KKNN), Neural Network (NNET), Least Absolute Shrinkage and Selection Operator (LASSO) were used to select the best model. With six of eight AUC values more than 0.9 (f2), the model based on RF was validated as the best one in both TCGA and GEO databases.

The immune cell infiltration varied significantly among the ATC subgroups C1, C2, and C3

The immune cell infiltration between TNC and ATC, calculated by three methods (Pompimol Charoentong’s algorithm, ssGSEA and CIBERSORT), was presented in a heatmap. A majority of immune cell types could infiltrate the tumor microenvironment in ATC (f3). Consensus clustering based on the previous 5 genes was harnessed to categorize ATC into three distinct groups: C1, C2, and C3 (f3). The immune cell infiltration in the three ATC subgroups was further calculated using three well-recognized methods. Approximately 85% of the immune cell types, including activated CD4 T cells, central memory CD8 T cells, effector memory CD8 T cells, dendritic cells, and macrophages, exhibited differential infiltration among C1-C3. Significant differences in activated CD8 T cell infiltration were observed in the results from the first two methods (f3).

The predicted immunotherapy response differed among the ATC subgroups C1, C2, and C3

The expression of key immune checkpoints was compared among the three groups (f4). CTLA4, CD80, and CD86 (CTLA4 system) were most highly expressed in the C2 group and least expressed in the C1 group (f4). LAG3 and FGL1 exhibited a similar expression pattern to the CTLA4 system (f4). Inhibitory markers of CD8 T cells, including T cell exhaustion (TEX) and regulatory T cells (Treg), were more highly expressed in the C2 and C3 groups, while markers of T cells in a stress response state (T sr) were relatively higher in the C1 and C2 groups (f4). Activating markers of CD8 T cells, such as the cGAS-STING score and CD8 effector T cells (Teff), showed a similar expression pattern to TEX and Treg (f4). Other immune markers, including IFN-gamma, CD80, and dysfunction score, were also enriched in the C2 and C3 groups (f4). According to the Tumor Immune Dysfunction and Exclusion (TIDE) analysis, patients in the C1 group may benefit from anti-PD1 therapy, whereas those in the C2 and C3 groups may be more suitable for cytotoxic T lymphocyte (CTL) therapy (f4).

12 innovative genes identified by WGCNA effectively distinguished C1-C3 groups using machine learning

Based on the weighted gene co-expression network analysis (WGCNA), the purple gene module was most negatively correlated with the grouping, TEX, CTLA-4, LAG-3, PD-L1, and immune dysfunction score, while the yellow module was most positively correlated with these markers (f5). Further analysis, overlapping the differentially expressed genes between each pair of C1-C3 groups (P value > 0.05 and a fold change (|FC|) > 2) with the combined gene set of both the yellow and purple modules, identified 12 genes (HLF, BCL2, HHEX, LRP2, FOXE1, FAM189A2, TSHR, EPB41L4B, OCLN, NEBL, ATP8A1, and TMEM30B) that may play an important role in the immune therapy of ATC (f5). The expression of these 12 genes showed consistent patterns within the C1-C3 groups (high in C1, intermediate in C2, and low in C3) (f5). Moreover, the individual expression of these 12 genes was positively related to TEX (f5). The subgrouping model by SVM, based on these 12 genes, effectively distinguished the C1-C3 groups with an AUC over 0.9 (f5).

Grouping based on the 12 genes across 9 melanoma cohorts receiving anti-PD1 therapy revealed significantly different response rates

The grouping model based on the 12 innovative genes was validated in available clinical trial cohorts. Except in the melanoma-PRJEB23709 cohort, the actual response rates to anti-PD1 therapy in the C2 and C3 groups were much lower than in the C1 group, consistent with our results (f4; f6). The correlation of individual gene expression with the predicted response rate to anti-PD1 therapy, as determined by TIDE, was further analyzed (f7). HLF, ATP8A1, and NEBL stood out due to their correlation coefficients over 0.4 and AUC values over 0.75 (f7). All components of the 12-gene signature were down-regulated in the TCGA thyroid cancer samples compared to the non-cancer samples (f7). Of the three genes, only HLF was of prognostic significance in the disease-free interval (DFI) of TCGA thyroid carcinoma (THCA) patients (f7). More importantly, its pro-survival role in prognosis was further corroborated in overall survival (OS) from three head and neck squamous cell carcinoma (HNSC) cohorts (GSE41613, GSE65858, TCGA). Its role was also validated in disease-specific survival (DSS) and progression-free interval (PFI) in the TCGA-HNSC cohort (f7).

HLF was a tumor suppressor in pan-cancer and could promote T-cell infiltration in ATC

Comparative analysis between cancerous (CA) and non-cancerous (NC) tissues revealed that HLF expression was generally lower in CA, observed in 17 out of 20 types (f8). The expression of HLF was negatively related to markers of apoptosis, cell cycle regulation, differentiation, DNA damage and repair, epithelial-mesenchymal transition (EMT), hypoxia, inflammation, invasion, metastasis, proliferation, and quiescence in TCGA-THCA patients (f8). Significant differences were also observed in the methylation status of genes in tumor-infiltrating lymphocytes (MeTIL), cytolytic activity (CYT), tertiary lymphoid structures (TLS), human leukocyte antigen (HLA), immunoinhibitor family, immunostimulator family, immune cell recruitment, and other immune markers between HLF high-expression and low-expression groups (f8).

Knockdown of HLF in ATC cell line CAL62 could greatly increase the migration of this cell line. A similar trend was also observed in the treatment of sorafenib in the lower layer (4μM) (f9). Dead/live cell staining using PI and calcein AM, proliferation staining by EDU, along with apoptosis detection demonstrated that HLF-knockdown CAL62 cells exhibited a lower proportion of dead cells, a higher proportion of live cells, an increased proliferation rate, and enhanced resistance to sorafenib (f9). The EMT pathway was upregulated in HLF-knockdown CAL62 cells (f10). Co-culturing CAL62 (in the lower layer) with peripheral blood mononuclear cells (PBMCs) (in the upper layer) using transwells indicated that HLF knockdown in CAL62 induced more PD1+ CD8 T cells (f11). Further co-culturing of CAL62/TCO1 ATC cells (in the lower layer) with Jurkat T cells (in the upper layer) showed that T cell recruitment decreased following HLF knockdown in the ATC cell line (f11). Finally, in 22 ATC clinical samples from our hospital, CD8 T cells and PD-L1/CD274 were detected, revealing an increase in PD-L1 and a decrease in CD8 T cells in the HLF low-expression samples (f12).

Discussion

To uncover the most distinct genes in different thyroid cancer subtypes, we screened BCL2, BHLHE40, MICAL2, TGM2, and TPO by comparing each pair of existing subtypes (f1). BCL-2 (B-cell lymphoma 2) is a protein that plays a crucial role in regulating apoptosis (ref. 14, ref. 15). Some researchers have found that BCL-2 might help thyroid cancer evade apoptosis, making it a suitable target for therapy (ref. 16–ref. 18). In our study, BCL-2 levels were relatively lower in ATC but higher in TNC, suggesting it may not be a promising target for ATC. BHLHE40 has been reported to be involved in the aggressiveness of various cancers, including colorectal, pancreatic, and endometrial cancers (ref. 19–ref. 22). One study focusing on ATC revealed that the lncRNA/H19-miR-454-3p/BHLHE40 axis could potentiate the progression of ATC (ref. 23). Consistent with this, BHLHE40 expression was significantly higher in TCA compared to TNC, and its expression in ATC was the second highest among all subtypes. MICAL2 is known as an oncogene in many cancers, such as pancreatic, ovarian, and gastric cancers (ref. 24–ref. 26), while there is no related research on thyroid cancer. Our results indicated that MICAL2 may also play a role in ATC. A similar phenomenon was observed with TGM2, which has been proven to be an oncogene with little research on ATC (ref. 27, ref. 28). TPO has been reported to predict higher metastasis and recurrence in PTC (ref. 29). It was nearly unexpressed in ATC in our results. Our findings revealed that BCL2, BHLHE40, MICAL2, TGM2, and TPO were expressed differently in ATC, FTC, MTC, PTC, and TNC. This inspired us to construct a signature based on these genes to help differentiate the subtypes. The model using the RF machine learning method worked precisely in differentiating TCA from TNC, ATC from DTC, MTC from DTC, and FTC from DTC (f2). However, the mechanisms by which these genes differentiate the subtypes and their roles in ATC, especially MICAL2, still require further research.

Since the five genes were so distinct, we continued to divide the most aggressive type, ATC, into more subtypes by consensus clustering, which may provide clues for specific therapy selection. The great variation of immune cell infiltration in the three subgroups indicated that immune therapy selection may differ among them. Further prediction by TIDE implied that the C1 group could more likely benefit from anti-PD1 therapy, while the C2 and C3 groups may be more suitable for CTL therapy (f3, f4). Until now, immunotherapy (anti-PD-1 and anti-PD-L1) has shown the most promising 1-year survival rate of approximately 40% for ATC patients without BRAF and MEK gene mutations (ref. 7). There is no CTL therapy for ATC patients in clinical trials. Our findings imply that grouping ATC patients may maximize the efficacy of immune therapy. More real-world clinical data are still needed to test this hypothesis.

To further characterize the gene expression patterns among C1, C2, and C3 groups in ATC, WGCNA was used to identify the most associated gene modules. Twelve genes (HLF, BCL2, HHEX, LRP2, FOXE1, FAM189A2, TSHR, EPB41L4B, OCLN, NEBL, ATP8A1, and TMEM30B) were selected from the most related modules due to their distinct expression among the three subgroups. A more precise model based on these 12 genes was created to replace the previous 5-gene signature. More importantly, when the new model was retrospectively applied to the available real-world melanoma clinical cohort, it validated that the C1 group was more suitable for anti-PD1 therapy (f5, f6). However, we must acknowledge that one result from the nine cohorts did not comply with our conclusion. In the future, we hope to test the model using data from ATC patients receiving anti-PD1 therapy or CTL therapy.

HLF was chosen from the 12 genes for final experimental validation because it performed well in both immune therapy response prediction and prognosis prediction in external cohorts (f7). As a transcriptional activator, HLF has been validated as a tumor suppressor in triple-negative breast cancer and ovarian cancer (ref. 30, ref. 31). However, one study indicated that it could promote the development of hepatocellular carcinoma and resistance to sorafenib (ref. 32). No research on HLF has been found in ATC. Our study indicated that the knockdown of HLF could promote the migration, proliferation, survival and sorafenib resistance of ATC cell lines (f9). The up-regulation of EMT pathway in the HLF-knockdown group may explain its role as a tumor suppressor in ATC (f10). The increased T-cell exhaustion, indicated by up-regulated PD-1 after 48 hours, and the dampened T-cell recruitment were also observed following the decrease in HLF expression in ATC cell lines (f11, f12). More interestingly, PD-L1/CD274 was also up-regulated in the HLF-low-expression group (f10, f12). These results suggest that HLF may be necessary for immune cells to function normally in the ATC tumor environment. Although further research is needed to explore HLF’s role in the ATC microenvironment, to our knowledge, this study is the first to investigate HLF in ATC. We will continue to explore this in our future studies.

Conclusion

In summary, our study identified five genes with distinct expression patterns across all subtypes of thyroid cancer. A signature based on these five genes can precisely distinguish between the subtypes. Additionally, our group developed a 12-gene signature in ATC that can predict the response to anti-PD1 therapy to some extent. The tumor suppressor role of HLF was validated in ATC cell lines through in vitro experiments.

References

1. A Jannin, A Escande, A Al Ghuzlan, P Blanchard, D Hartl, B Chevalier. Anaplastic thyroid carcinoma: an update.. Cancers (Basel). (, 2022. [DOI]
2. M de Ridder, E Nieveen van Dijkum, A Engelsman, E Kapiteijn, HJ Klümpen, CRN Rasch. Anaplastic thyroid carcinoma: a nationwide cohort study on incidence, treatment and survival in the Netherlands over 3 decades.. Eur J Endocrinol. (, 2020. [DOI]
3. K Sehgal, A Serritella, M Liu, A ON, C Nangia, T Pappa. A phase I/II trial of sapanisertib in advanced anaplastic and radioiodine refractory differentiated thyroid carcinoma.. J Clin Endocrinol Metab. (, 2024. [DOI]
4. KA Araque, S Gubbi, J Klubo-Gwiezdzinska. Updates on the management of thyroid cancer.. Horm Metab Res. (, 2020. [DOI]
5. V Tiedje, M Stuschke, F Weber, H Dralle, L Moss, D Führer. Anaplastic thyroid carcinoma: review of treatment protocols.. Endocr Relat Cancer. (, 2018. [DOI]
6. A Ocanto, L Torres, F Couñago. Current status of anaplastic thyroid carcinoma.. World J Clin Oncol. (, 2024. [DOI]
7. ET Pavlidis, IN Galanis, TE Pavlidis. Update on current diagnosis and management of anaplastic thyroid carcinoma.. World J Clin Oncol. (, 2023. [DOI]
8. J Jeschke, M Bizet, C Desmedt, E Calonne, S Dedeurwaerder, S Garaud. DNA methylation-based immune response signature improves patient diagnosis in multiple cancers.. J Clin Invest. (, 2017. [DOI]
9. Ó Lapuente-Santana, M van Genderen, PAJ Hilbers, F Finotello, F Eduati. Interpretable systems biomarkers predict response to immune-checkpoint inhibitors.. Patterns (N Y). (, 2021. [DOI | PubMed]
10. H Yuan, M Yan, G Zhang, W Liu, C Deng, G Liao. CancerSEA: a cancer single-cell state atlas.. Nucleic Acids Res. (, 2019. [DOI]
11. E Lee, HY Chuang, JW Kim, T Ideker, D Lee. Inferring pathway activity toward precise disease classification.. PloS Comput Biol. (, 2008. [DOI | PubMed]
12. G Cao, P Li, X He, M Jin, M Li, S Chen. FHL3 contributes to EMT and chemotherapy resistance through up-regulation of slug and activation of TGFβ/smad-independent pathways in gastric cancer.. Front Oncol. (, 2021. [DOI | PubMed]
13. P Li, G Cao, Y Zhang, J Shi, K Cai, L Zhen. FHL3 promotes pancreatic cancer invasion and metastasis through preventing the ubiquitination degradation of EMT associated transcription factors.. Aging (Albany NY). (, 2020. [DOI | PubMed]
14. C Zhang, C Bo, L Guo, P Yu, S Miao, X Gu. BCL2 and hsa-miR-181a-5p are potential biomarkers associated with papillary thyroid cancer based on bioinformatics analysis.. World J Surg Oncol. (, 2019. [DOI | PubMed]
15. V Gunda, KA Sarosiek, E Brauner, YS Kim, S Amin, Z Zhou. Inhibition of MAPKinase pathway sensitizes thyroid cancer cells to ABT-737 induced apoptosis.. Cancer Lett. (, 2017. [DOI | PubMed]
16. D Champa, MA Russo, XH Liao, S Refetoff, RA Ghossein, A Di Cristofano. Obatoclax overcomes resistance to cell death in aggressive thyroid carcinomas by countering Bcl2a1 and Mcl1 overexpression.. Endocr Relat Cancer. (, 2014. [DOI]
17. D Kaloni, ST Diepstraten, A Strasser, GL Kelly. BCL-2 protein family: attractive targets for cancer therapy.. Apoptosis. (, 2023. [DOI | PubMed]
18. CS Mitsiades, P Hayden, V Kotoula, DW McMillin, C McMullan, J Negri. Bcl-2 overexpression in thyroid carcinoma cells increases sensitivity to Bcl-2 homology 3 domain inhibition.. J Clin Endocrinol Metab. (, 2007. [DOI]
19. K Asanoma, H Yagi, I Onoyama, L Cui, E Hori, M Kawakami. The BHLHE40−PPM1F−AMPK pathway regulates energy metabolism and is associated with the aggressiveness of endometrial cancer.. J Biol Chem. (, 2024. [DOI | PubMed]
20. A Ji, H Li, X Fu, Y Zhang, Y Liu. Long non-coding RNA NEAT1 induced by BHLHE40 activates Wnt/β-catenin signaling and potentiates colorectal cancer progression.. Cell Div. (, 2024. [DOI | PubMed]
21. S Yang, D Zhang, Q Sun, H Nie, Y Zhang, X Wang. Single-cell and spatial transcriptome profiling identifies the transcription factor BHLHE40 as a driver of EMT in metastatic colorectal cancer.. Cancer Res. (, 2024. [DOI]
22. Y Cao, X Wang, Y Liu, P Liu, J Qin, Y Zhu. BHLHE40 inhibits ferroptosis in pancreatic cancer cells via upregulating SREBF1.. Adv Sci (Weinh). (, 2024. [DOI | PubMed]
23. Y Wu, J Yang, H Zhang, J Cheng, P Lei, J Huang. LncRNA H19 Influences Cellular Activities via the miR-454-3p/BHLHE40 Axis in Anaplastic Thyroid Carcinoma.. Horm Metab Res. (, 2024. [DOI]
24. B Garg, S Khan, DS Babu, E Mose, K Gulay, S Sharma. MICAL2 is a super enhancer associated gene that promotes pancreatic cancer growth and metastasis.. bioRxiv. (, 2024. [DOI]
25. Z Liu, B Sun, A Xu, J Tang, H Zhang, J Gao. MICAL2 implies immunosuppressive features and acts as an independent and adverse prognostic biomarker in pancreatic cancer.. Sci Rep. (, 2024. [DOI | PubMed]
26. T Xia, F Ye, W Zhao, P Min, C Qi, Q Wang. Comprehensive analysis of MICALL2 reveals its potential roles in EGFR stabilization and ovarian cancer cell invasion.. Int J Mol Sci. (, 2023. [DOI]
27. W Chang, W Gao, D Liu, B Luo, H Li, L Zhong. The upregulation of TGM2 is associated with poor prognosis and the shaping of the inflammatory tumor microenvironment in lung squamous cell carcinoma.. Am J Cancer Res. (, 2024. [DOI]
28. N Zhang, S Gao, H Peng, J Wu, H Li, C Gibson. Chemical proteomic profiling of protein dopaminylation in colorectal cancer cells.. J Proteome Res. (, 2024. [DOI]
29. X Li, R Cheng. TPO as an indicator of lymph node metastasis and recurrence in papillary thyroid carcinoma.. Sci Rep. (, 2023. [DOI | PubMed]
30. T Han, T Chen, L Chen, K Li, D Xiang, L Dou. HLF promotes ovarian cancer progression and chemoresistance via regulating Hippo signaling pathway.. Cell Death Dis. (, 2023. [DOI | PubMed]
31. H Li, P Yang, J Wang, J Zhang, Q Ma, Y Jiang. HLF regulates ferroptosis, development and chemoresistance of triple-negative breast cancer by activating tumor cell-macrophage crosstalk.. J Hematol Oncol. (, 2022. [DOI | PubMed]
32. DM Xiang, W Sun, T Zhou, C Zhang, Z Cheng, SC Li. Oncofetal HLF transactivates c-Jun to promote hepatocellular carcinoma development and sorafenib resistance.. Gut. (, 2019. [DOI]