Skip to main content

Immunological characteristics of immunogenic cell death genes and malignant progression driving roles of TLR4 in anaplastic thyroid carcinoma


Anaplastic thyroid carcinoma (ATC) was a rare malignancy featured with the weak immunotherapeutic response. So far, disorders of immunogenic cell death genes (ICDGs) were identified as the driving factors in cancer progression, while their roles in ATC remained poorly clear. Datasets analysis identified that most ICDGs were high expressed in ATC, while DE-ICDGs were located in module c1_112, which was mainly enriched in Toll-like receptor signalings. Subsequently, the ICD score was established to classify ATC samples into the high and low ICD score groups, and function analysis indicated that high ICD score was associated with the immune characteristics. The high ICD score group had higher proportions of specific immune and stromal cells, as well as increased expression of immune checkpoints. Additionally, TLR4, ENTPD1, LY96, CASP1 and PDIA3 were identified as the dynamic signature in the malignant progression of ATC. Notably, TLR4 was significantly upregulated in ATC tissues, associated with poor prognosis. Silence of TLR4 inhibited the proliferation, metastasis and clone formation of ATC cells. Eventually, silence of TLR4 synergistically enhanced paclitaxel-induced proliferation inhibition, apoptosis, CALR exposure and release of ATP. Our findings highlighted that the aberrant expression of TLR4 drove the malignant progression of ATC, which contributed to our understanding of the roles of ICDGs in ATC.

Peer Review reports


Anaplastic thyroid carcinoma (ATC) was a rare malignancy that accounted for 50% of thyroid cancer deaths [1, 2]. ATC progressed rapidly and was characterized by local and distant metastases. Currently, surgical resection was the main clinical treatment, supplemented by radiotherapy and chemotherapy [3,4,5]. However, the prognosis of ATC remained very poor, with the median survival of patients often less than 6 months [6]. Hence, it was urgent to explore the crucial molecules that driving the malignant progression of ATC and develop new clinical treatment strategies.

In recent years, the increasing popularity of immunotherapy had greatly encouraged researchers to investigate its potential for treating ATC patients [7]. Clinical trials for patients with locally advanced or metastatic ATC have shown that sparazumab treatment resulted in an objective response rate of 19%. However, the median progression-free survival and overall survival were only 1.7 months and 5.9 months, respectively indicating limited efficacy of single-agent immunotherapy for ATC [8]. And a current clinical study of CSF-1R inhibitor combined with PD-1 inhibitor showed remission in only 37% of patients with advanced thyroid cancer [9]. In another study of immunotherapy combined with chemoradiotherapy for ATC patients, a combination of palizumab, docetaxel, and doxorubicin with volume-modulated intensity radiotherapy was performed, but regrettably all patients died within 6 months [10]. These indicated an urgent need for highly effective strategies to increase sensitivity to ATC immunotherapy.

Immunogenic cell death (ICD) was a form of cell death that was triggered by different cell stressors, such as chemotherapy drugs, oncolytic viruses, radiotherapy and photodynamic therapy [11]. This type of cell death was distinguished by the release of damage-associated molecular patterns (DAMPs), including HMGB1, calreticulin (CALR), ATP, and heat shock proteins. These DAMPs played crucial roles in anti-tumor immunity [12]. Emerging evidences suggested that ICD induction of tumor cells and transforming into anti-tumor vaccines could significantly inhibit tumor growth. Kwong et al. showed that paclitaxel (PTX) induced ICD-related DAMPs in ovarian cancer and significant antitumor responses in vivo [13]. And the combination of Trifluridine and Oxaliplatin eliminated type 2 tumor-associated macrophages, leading to higher cytotoxic CD8 + T cell infiltration and activation, thereby inducing ICD to improve curative effect in patients with colorectal cancer [14]. Oleander was found to trigger endoplasmic reticulum stress and induce ICD-mediated immune destruction of breast cancer cells [15]. Therefore, induction of ICD might act crucial roles in improving the efficacy of ATC immunotherapy. However, there was no studies associated with ICD reported in ATC.

Here, we had integrated the datasets from multiple laboratories or array platforms to conduct the largest ATC cohort gene expression studies. To address the challenges of ATC samples collection and the systematic bias between different datasets, efforts could be made to identify the unique and specifically altered ICDGs signature in ATC. We demonstrated significant ICDGs and constructed the gene co-expression network in ATC. Then we established an ICD score evaluation model to score ATC samples and divided them into two groups. The differential expressed genes (DEGs) of high and low ICD score subgroups were identified and analyzed for functional enrichment. Nonetheless, the differences in immune landscapes and immune checkpoint expression between the two subgroups were compared. Further, we identified the signature of 5 ICDGs to reveal crucial signals for driving the malignant progression of ATC. We found that the upregulation of TLR4 significantly promoted the malignant characteristics of ATC, and inhibition of TLR4 induced ICD in combination with paclitaxel. To our knowledge, this was the first time that functions of ICDGs had been studied in ATC. Moreover, our results demonstrated that the disorder of immunogenic cell death genes were crucial driving factors leading to the malignancy of ATC, and the intervention of TLR4 might be beneficial to the immunotherapy of ATC.

Materials and methods

Microarray information

The microarray datasets (GSE65144, GSE33630, GSE29265 and GSE76039) in the GEO database ( were downloaded [16], in which researchers deposit their microarray data related to various diseases. The dataset of GSE65144 included 12 ATCs and 13 normal tissues [17]; GSE33630 contained 11 ATCs, 49 papillary thyroid carcinomas (PTC) and 45 normal tissues [18]; GSE29265 dataset contained 9 ATC, 20 normal tissues and 20 PTCs; GSE76039 had 20 ATCs and 17 poorly-differentiated thyroid carcinomas (PDTCs) [19, 20]. The above datasets were based on GPL570 (Affymetrix HT HG-U133 plus PM Array). The four datasets were corrected and its dimensionality were reduced by the Remove Batch Effect function of the limma package and principal component analysis in the R4.1.3 environment [21, 22]. And analysis workflow of this research was showed in Fig. 1.

Fig. 1
figure 1

Analysis workflow of this research

Identification and MEGENA analysis of immunogenic cell death genes

The list of ICDGs was from the meta-analysis of ICD-related literature conducted by Abhishek et al. [23]. The differentially expressed analysis of ATC vs. normal groups were also performed using limma package, and the volcano plot and heatmap plot were made by ggpubr and pheatmap packages respectively. The protein-protein interaction (PPI) network analysis were conducted by the STRING database and optimized by the Cytoscape 3.6.1 software [24]. The ICDGs were analyzed by the Multiscale Embedded Gene co-Expression Network analysis (MEGENA) package for calculating a multi-scale network including different possible variations in gene interactions in R4.1.3 [25]. Then the network enrichment analysis of all node genes of c1_112 was implemented by the ClueGO v2.5.2 plug-in unit of Cytoscape 3.6.1 software.

The establishment of ICD score model and enrichment analysis

ATC samples in the dataset were scored according to ICDGs expression and ssGSEA algorithm in GSVA package. The mean score of all ATC sample was calculated, and the samples were divided into high ICD score group and low ICD score group based on this. Then the difference genes between the subgroups were identified by the limma package, whose thresholds were log2 (fold change) > 1.5 or log2 (fold change) < -1. The DEGs was mapped into a pointplot by the ggplot2 package. The GO function annotation, KEGG pathway enrichment and GSEA analysis were performed to further elucidate the biological roles of DEGs between the subgroups by the clusterProfiler package in the R4.1.3 environment.

Immune characteristics analysis

ESTIMATE, an algorithm that detect the proportion of stromal and immune cells based on gene expression characteristics in tumor samples, was utilized for the calculation of stromal score, immune score, ESTIMATE score and tumor purity of high or low ICD score groups [26]. CIBERSORT that was a deconvolution algorithm was used to calculate the fraction of 22 immune cells in each samples in two subgroups [27]. According to the gene expression of 28 published immune cell genesets combined with ssGSEA algorithm, the differences in the extent of infiltration of 28 immune cell types between the two groups were compared [28].

The establishment of Lasso-logistic model and prognostic analysis

The Lasso Logistic model by the glmnet package mainly adjusted the λ parameter to complete the variable selection, so that the estimated value of irrelevant variables was 0, and further learned and weighted to predict the ICDGs most related to ATC [29]. The prognostic analysis of ICDGs in thyroid cancers was applied by Kaplan-Meier plotter (

WGCNA analysis

The WGCNA package was used to analyze gene expression patterns in multiple samples [30]. WGCNA was utilized for co-expression analysis of DEGs between high and low TLR4 groups. The mean linkage and Pearson’s correlation methods were used to plot a sample dendrogram with stromal score, immune score, ESTIMATE score and tumor purity. The DEGs were assigned to 6 modules, and 26 key genes of red module were obtained by calculating module membership (MM) and gene significance (GS), followed by GO Enrichment analysis.

Cell culture and transfection

The 8505 C (Anaplastic thyroid carcinoma cell line) and Nthy-ori 3 − 1 (Normal human thyroid cell line) were purchased from Fenghui Biotechnology. The CAL62 (Anaplastic thyroid carcinoma cell line) was purchased from Procell Life Science&Technology. Cell Line authentication: all cell lines had been performed authenticated using Short Tandem Repeat (STR) analysis on 2021 in Procell Life Science&Technology. The CAL62, 8505 C and Nthy-ori 3 − 1 were cultured in DMEM or RPMI-1640 medium supplemented with 10% FBS (Gbico, USA) and in culture flask (Jet Biofil, TCF012050). 8505 C cell was seeded in 6-well plated with 40–50% confluence, and the siRNAs-TLR4 were transfected with jetPRIME (Polyplus, USA) according to the instruction. The sense sequences of TLR4 siRNA#1 was CAGCAAUCUCACUUCUGUA and antisense sequences was UACAGAAGUG- AGAUUGCUG. The sense sequences of TLR4 siRNA#2 was GUGUGUUUCCAU- GUCUCAU and antisense sequences was AUGAGACAUGGAAACACAC.

Western blot and qRT-PCR

The operation processes of western blot and qRT-PCR were based on the previous literature published by our team [31]. The primary antibody TLR4 (ABclonal Technology, western blot 1:1000, China) and GAPDH (Proteintech, 1:1000, China) were applied for detection of proteins by the electrophoresis apparatus ( Bio-rad, USA). qRT-PCR primers: The forward primer of TLR4 was 5’-GTGCCTCCATTTCAGCTCTG-3’ and the reverse primer was 5’-CAAAGATAC- ACCAGCGGCTC-3’. The forward primer of β-ACTIN was 5’-ACCTTCTACAATG- AGCTG- CG-3’ and the reverse primer was 5’-CCTGGATAGCAACGTACATGG-3’.


The operation processes of immunohistochemistry was based on the previous literature published by our team [31]. The primary antibody TLR4 were applied for immunohistochemistry (TLR4, immunohistochemistry 1:100). The tissue sections ATC1, ATC2 and ATC3 were acquired from the tumor tissue sample bank of Zhejiang Provincial People’s Hospital, which were clinicopathologically diagnosed as anaplastic thyroid carcinoma. The process was approved by the Ethics Committee of Zhejiang Provincial People’s Hospital.

Cell proliferation, migration, invasion and clone formation

For cell proliferation assay, the treated ATC cells were seeded in the 96-well plate with 30–40% confluence for 48 h growth, and then the viability was detected by the Cell Counting Kit-8 (Fdbio, China) on the microplate reader (Bioteck, USA). For migration assay, the treated ATC cells were seeded in the 12-well plate with 100% confluence replaced with serum-free medium. The cells were scratched at the 0 h and imaged, then observing their migration capacity at 48 h and imaged. For cell invasion assay, the transfected ATC cells were seeded in the upper chamber pre-coated with 5% Matrigel (BD Biosciences, USA). Medium containing 0% and 10% serum were added to the upper and lower chambers for incubating 48 h, respectively. The invasive cells were then fixed with methanol and stained with 0.1% crystal violet (Applygen, China). The infiltrating cells were photographed after the upper chamber cells were gently wiped. For clone formation assay, the treated ATC cells were seeded with 500/12-well to incubate for 1 week, and then cells were fixed with methanol, stained with 0.1% crystal violet (Applygen, China) and photographed.

Flow cytometry and ATP release detection

For cell apoptosis assay, the treated ATC cells were centrifuged and resuspended with 1X binding buffer (Beyotime, China). Annexin V (10 µl) and PI (5 µl) were added and incubated for 30 min at room temperature, protected from light, before detection with the Flow cytometer (Beckman, USA). For detection of exposed CALR on the cell membrane surface, treated cells were harvested and blocked with 1.5% bovine serum albumin for 30 min. Then, samples were incubated with primary antibody CALR (Proteintech, 1:200, USA) at room temperature for 40 min followed by secondary antibody PE anti-Rabbit (Abcam, 1:200, USA) for 30 min, and then detect on the flow cytometer. For ATP release detection assay, cell supernatants were collected and then the content of released ATP was detected using the ATP Assay Kit (Applygen, China) according to the manufacturer’s instructions.

Zebrafish xenograft models

Zebrafish were fed in a 28.5 °C recirculating water system with a standard 12-hour light and dark cycle and fed pellets twice a day. For tumor growth models, zebrafish embryos were mechanically delaminated 2 days after fertilization and anesthetized with 0.15 mg/mL tricaine. 8505 C-NC or siTLR4 cells were stained with 5 µg/mL DiI (Beyotime, C1036) at 37℃ for 30 min. Then, the cells were cleaned twice and resuspended with PBS, microinjected into the perivitelline space (about 500 cells/piece) with the glass capillary needle and grown at 34 °C. At the day 3, the zebrafish were anesthetized and then imaged under the fluorescence microscope. Image pro plus software was used to analyze and calculate the optical density of fluorescence signal.

Data analysis

All data were represented as the mean ± SD of three independent tests. Each experiment was independently repeated three times. The unpaired Student’s t-tests with GraphPad Prism 7 was utilized to analyze statistical differences between the two groups. *indicated the statistical significance (*P < 0.05, **P < 0.01, ***P < 0.001).


Identification and enrichment analysis of DE-ICDGs in anaplastic thyroid carcinoma

As stated in our previous study, GSE76039, GSE33630, GSE65144 and GSE29265 were captured and integrated to correct the batch effect while preserving the difference. The combined dataset was analyzed for finding the DE-ICDGs in ATC. As shown in Fig. 2A, a total of 1994 differentially expressed genes were identified from the above dataset. Subsequently, we investigated the expression of ICDGs in normal tissues and ATC (Fig. 2B-C). The analysis of the dataset revealed that a majority of the ICDGs were up-regulated in ATC, while only a few were down-regulated. Notably, LY96, CASP1, ENTPD1, TLR4, and PDIA3 displayed the most significant changes. Furthermore, we constructed a protein–protein interaction (PPI) network with the ICDGs and found that IL1B, IFNG, CD4, IL10 and TLR4 play crucial roles (Fig. 2D). The above results implied that ICDGs may be highly correlated with the progression of ATC.

Fig. 2
figure 2

Identification of differential expressed immunogenic cell death genes and network modules analysis in ATC. (A) Volcano plot of DEGs identified in the combined dataset according to the fold change and adjusted P value. (B) The expression of ICDGs in volcano plot. (C) The heatmap of the expression profiles of ICDGs in thyroid cancer samples. (D) Protein–protein interactions network of ICDGs. (E) The network modules of 1994 DEGs. (F) The submodule where DE-ICDGs was located. (G) Network analysis diagram of DE-ICDGs and their functions

The network modules and enrichment analysis of immunogenic cell death genes

The MEGENA analysis were used to construct a gene co-expression network of 1994 DEGs, and further clarify the biomolecular differences associated with the malignant features of ATC. By identifying the interactions of these DEGs, we agglutinated 83 closely linked gene co-expression network modules (Fig. 1E). Then, the co-expression network child module c1_112 which DE-ICDGs were located in was analyzed (Fig. 1F). This module was composed of two gene clusters centered on LAIR1, C1QA, C3AR1, FPR3 and MS4A6A respectively. To investigate the signaling of ATC malignant phenotypes mediated by DE-ICDGs, network enrichment analysis was conducted on all node genes of c1_112 to annotate modules with representative biological function categories (Fig. 1G). Our analysis results demonstrated that these genes can be summarized into 7 specific clusters, and these specific clusters are mainly enriched in 4 biological processes, including Toll-like receptor signaling pathway, legionellosis, pertussis and systemic lupus erthematosus. The results showed that the DE-ICDGs in ATC were highly correlated with changes in these biological processes.

The immunogenic cell death related score of ATC and underlying signal pathways in different subgroups

Considering the heterogeneity and complexity of individual immunogenic cell death genes expression and subsequent identification of key genes, we performed ssGSEA to quantify the expression of ICDGs in ATC and defined the result as the immunogenic cell death related score (ICD score). Next, ICD score was evaluated for all ATC samples. After calculating the mean value, all patients in the integrated cohort were assigned to the high ICD score and low ICD score groups (p < 0.001) (Fig. 3A-B). To explore the differences in the biological behaviors between ICD subgroups, we identified the key DEGs and signal pathways in each subgroups for comprehending the molecular mechanism in regulation of ATC progression (Fig. 3C). Here, we identified 35 upregulated genes and 5 downregulated genes. To assess the functional annotation of these genes in ATC, we conducted GO and KEGG analysis. Our findings revealed that these genes were significantly enriched in activities associated with immunity, such as humoral immune response, immunoglobulin production, production of molecular mediator of immune response and leukocyte migration in the GO analysis (Fig. 3D), as well as cytokine − cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, chemokine signaling pathway and Toll − like receptor signaling pathway in KEGG analysis (Fig. 3E). To further identify specifically activated signaling pathways in the high ICD score group, we performed a GSEA comparison between the two groups. As shown in Fig. 3F, enriched gene sets were associated with immune pathways such as defense response, humoral immune response, positive regulation of immune system process and regulation of immune system process. These results indicated that high ICD score was associated with the immune characteristics in ATC.

Fig. 3
figure 3

Identification of key differential expressed genes and underlying signal pathways between immunogenic cell death score subgroups. (A) The heatmap of ICD scores of ATC samples. (B) The boxplot of the different ICD score subgroups. (C) The pointplot of DEGs expression between the subgroups. (D) GO and (E) KEGG analysis of the DEGs. (F) The KEGG pathway enrichment of the DEGs by GSEA

Comparison of immune characteristics between the ICD score subgroups

In view of the high relationship of ICD score and immune characteristics, ESTIMATE, CIBERSORT, and ssGSEA were performed to comprehend the differences in immunological functions of the subgroups in ATC. ESTIMATE analysis showed that the high ICD score group had a significantly higher immune score, estimate score and stromal score than the low ICD score group, while the tumor purity was lower (p < 0.001) (Fig. 4A–D). These results demonstrated that the high ICD score group had higher proportions of immune and stromal cells, and lower proportion of tumor cells. Furthermore, CIBERSORT analysis demonstrated that ICD score has positive correlations with the proportions of T cells CD4 memory activated, T cells gamma delta, and negative correlations with the proportions of macrophages M0 among the 22 tumor-infiltrating lymphocyte types (p < 0.05) (Fig. 4E). Simultaneously, ssGSEA analysis showed 28 immune cell subtypes (especially activated B cells, activated CD8 T cells, effector memory CD4 T cell, effector memory CD8 T cell, immature B cell, MDSC, T follicular helper cell and type 1 T helper cell) to be highly expressed in the high ICD score group (p < 0.05) (Fig. 4F). We further compared the expression of immune checkpoint molecules in the two groups (p < 0.05) (Fig. 4G). According to the data, the mRNA expression of several immune checkpoint molecules revealed that PDCD1LG2, CD274, HAVCR2, CTLA4 and LAG3 were upregulated in the high ICD score group, but there was no significant difference in PDCD1. These results suggested that ICD score had inalienable internal relationship with the immune characteristics of ATC, and high ICD score indicated extensive alterations in the expression of immune checkpoint molecules.

Fig. 4
figure 4

Comparison of immune characteristics between the immunogenic cell death score subgroups. Comparison of the (A) stromal score, (B) immune score, (C) ESTIMATE score, (D) tumor purity based on ESTIMATE, (E) proportion of immune cells based on CIBERSORT and (F) expression of immune cells based on ssGSEA between the high and low ICD score groups. (G) Comparison of the expression of immune checkpoints between the ICD score subgroups

Identification and validation of key immunogenic cell death genes

We analyzed the DE-ICDGs in ATC, so that the crucial ICD genes driving the malignant progression of ATC could be identified. The coefficients of 5 genes (TLR4, ENTPD1, LY96, CASP1 and PDIA3) were found to be non-zero by the Lasso-logistic regression model (Fig. 5A). TLR4 with the largest coefficient, was found to be significantly upregulated in ATC and was associated with poor prognosis in patients with thyroid cancer (Fig. 5B-C). In addition, the expression of TLR4 was positively correlated with ICD score in ATC samples (Fig. 5D). For validation, we first detected the expression of TLR4 in ATC and normal thyroid cell lines. As shown in results in Fig. 5E-F, the mRNA and protein levels of TLR4 in ATC was distinctly higher than that in normal (p < 0.05). Simultaneously, similar results were obtained by the immunohistochemical staining after determining the expression of TLR4 in ATC and normal thyroid tissues (Fig. 5G). The above results suggested that immunogenic cell death gene TLR4 may act crucial roles in the malignant progression of ATC.

Fig. 5
figure 5

TLR4 regulated the malignant characteristics of ATC. (A) The cross-validation applied to adjust parameter selection in the LASSO model. (B) The expression of TLR4 between ATC and normal tissues. (C) The overall survival of patients with thyroid carcinoma based on TLR4 expression by Kaplan-Meier plotter. (D) Correlation between TLR4 expression and ICD score in ATC samples. (E) The mRNA and (F) protein level of TLR4 in different ATC cell lines (8505 C and CAL62) and normal thyroid cell lines (Nthy-ori 3 − 1, NTHY). (G) IHC staining of TLR4 in ATC samples and normal tissue. (H) Silencing effect of siRNA-TLR4. The effect of siTLR4 on (I) the proliferation for 48 h measured by CCK-8, (J) migration ability for 48 h measured by wound-healing assay, (K) invasive ability for 48 h measured by transwell assay and (L) colony formation ability for two weeks in 8505 C cells

Downregulation of TLR4 inhibited the malignant phenotypes of ATC

Since TLR4 was upregulated in ATC, whether TLR4 alteration would influence the malignant phenotypes of ATC needs to be further elucidated. Then, siRNA-knockdown ATC model of TLR4 was constructed for subsequent experiments (Fig. 5H). And our data showed that reducing the expression of TLR4 effectively impeded the proliferation of ATC cells (p < 0.001) (Fig. 5I). Scratch repair experiments showed that silencing TLR4 reduced the migration ability of ATC cells (Fig. 5J). Next, the invasion ability of ATC cells was significantly weakened after the inhibition of TLR4 (Fig. 5K). Moreover, low levels of TLR4 inhibited the colony formation of ATC cells by clonal formation experiments (Fig. 5L). These findings demonstrated that TLR4 could modulate the malignant progression of ATC. However, the molecular mechanisms which TLR4 involved in remained obscure.

WGCNA analysis of hub genes related with TLR4 and functional enrichment

To construct a weighted co-expression network, TLR4-related modules and genes needed to be screened. The ATC samples were divided into two groups based on the TLR4 expression, using the TLR4 expression mean as the reference limit. Weighted gene co-expression network analysis was performed to screen the DEGs. After a series of adjustments for WGCNA parameters, the DEGs were divided into 6 modules by average linkage hierarchical clustering (Fig. 6A-B). We performed the correlation analysis between modules and traits to identify a module associated with TLR4 and immune score (Fig. 6C). The red module contained 326 genes and had the highest correlation with immune score (R = 0.82, P = 3e-05). 26 genes in the red module were selected as hub genes (module absolute membership > 0.8, gene significance > 0.5) (Fig. 6D). To determine the bio-functions of the 26 hub genes, we performed the GO enrichment analysis (Fig. 6E). Myeloid leukocyte activation, activation of immune response, leukocyte mediated immunity, T cell differentiation and microglial cell activation were the most frequently noted pathways in the functional enrichment analysis.

Fig. 6
figure 6

Identification of module genes associated with TLR4 by WGCNA. (A) Cluster dendrogram and module colors of ATC. (B) Clusters of DEGs based on the dissimilarity measure. (C) The correlation heatmap between module eigengenes and ESTIMATE results. (D) The scatter plot of red module eigengenes. (E) The GO analysis of 26 red module eigengenes

The inhibition of TLR4 enhanced paclitaxel-induced immunogenic cell death

Moreover, with considering that TLR4 was an important ICD gene and the reported roles in PTX-induced ICD, we first silenced TLR by siRNA transfection for 24 h and then treated with 1 nM PTX for 48 h. The results indicated that the repressed expression of TLR4 could effectively enhance PTX-induced proliferation inhibition and apoptosis of ATC cells (p < 0.01) (Fig. 7A-B). Simultaneously, a correlation analysis revealed a negative correlation between the expression of TLR4 and CALR whose exposure was one of the DAMPs in ICD (Fig. 7C). For verification, we observed more CALR exposure on the ATC cell membrane by flow cytometry (Fig. 7D). And the release of ATP was significantly increased in the TLR4 inhibition and PTX synergistic groups (p < 0.05) (Fig. 7E). In addition, we also found that silencing TLR4 could significantly reduce the level of CALR and inhibit the tumors growth in zebrafish xenograft models (p < 0.05) (Fig. 7F-G). The above contents suggested that silencing TLR4 could synergize with PTX to induce ICD in ATC.

Fig. 7
figure 7

Inhibition of TLR4 enhanced PTX-induced immunogenic cell death. (A) The effect of siTLR4 synergistic with PTX (1 nM) on the proliferation for 48 h measured by CCK-8 in 8505 C. (B) The apoptosis percentages of 8505 C after the synergistic treatment of siTLR4 and PTX (1 nM) for 48 h measured by flow cytometry. (C) Correlation between CALR expression and TLR4 in ATC samples. (D) The CALR exposure and (E) ATP release of 8505 C after the synergistic treatment of siTLR4 and PTX (1 nM) for 48 h. (F) The tumorigenesis of zebrafish in NC and siTLR4 groups. (G) The CALR exposure of NC and siTLR4 groups in zebrafish


Cancer immunotherapy was to stimulate the immune system to elicit the body’s immune response to the tumor, but tumor immune escape made its clinical efficacy limited [32, 33]. Accumulating evidence clearly underscored that ICD activated the innate immune system and actively modulates cancer immunotherapy [34, 35]. Understanding the molecular signature of ICDGs had the potential to identify straregies to overcome the immune escape “cold tumor”. This knowledge opened up new avenues for enhancing the effectiveness cancer immunotherapy. Despite the considerable efforts to characterize ICD-related genes or ICD-inducers in other cancers, the dynamic profiling of ICDGs in ATC remained far from being established [36, 23, 37,38,39]. In this study, large cohort data of ATC samples were integrated by us to analyze the signature of ICDGs in ATC. We found that most ICDGs were upregulated in ATC and had close interaction. Next, MEGENA analysis was performed to identify the submodules of ATC gene network that ICDGs belonged to. Preliminary, it was found that they mainly mediated Toll-like receptor signaling pathway, and then affected the malignant characteristics of ATC.

Further, the ICD score model was established for classifying ATC samples, with significant differences between subgroups.

Obviously, these genes were found to be strongly linked to immune response, immune system regulation, and cytokine signaling, as verified through functional enrichment analysis. Following this, a comparison of the immune characteristics was conducted between two subgroups. The tumor purity was lower in the high ICD score group because of the massive infiltration of immune cells and stromal cells in its tumor microenvironment. These overabundant immune cells had the largest difference in the proportion of intermediate-class T cell subtypes, including activated CD8 T cells, effector memory CD4 T cell, effector memory CD8 T cell, immature B cell, T follicular helper cell and type 1 T helper cell. It was implied that ICDGs was associated the dysfunction of T cell in ATC. T cell dysfunction was one of the pivotal mechanisms of tumor immune escape [40]. During tumor development, T cell dysfunction/exhaustion occurred and various checkpoint inhibitory receptors were upregulated, thus limiting the survival and functions of T cells [41, 42]. Emerging evidences had shown that the increased immune checkpoint molecules were described as the hallmarks of T cell exhaustion. Indeed, in our study, the expression of PDCD1, PDCD1LG2, CD274, HAVCR2, CTLA4 and LAG3 were upregulated in the high ICD score group. These results suggested that high ICD score represented extensive alterations in immune characteristics of ATC.

Creatively, the Lasso-logistic model was performed to analyze ICDGs, and we discovered a signature of 5 ICDGs in ATC. This findings helped us gain a deeper our understanding of the polymorphism of ICD in cancer. These ICDGs (TLR4, ENTPD1, LY96, CASP1 and PDIA3) were rarely reported in ATC. Notably, the expression of ENTPD1 in ATC was higher than that in normal tissues, suggesting its significance in regulating metabolic pathway and acting important roles in thyroid cancer [43, 44]; LY96 was highly associated with increased intratumoral lymphocytic infiltration in oncocytic PDTC [45]; CASP1 was identified as the most promising biomarker for the diagnosis or treatment of PTC and involved in molecular changes in brain metastatic PTC [46, 47]; patients with low PDIA3 expression had lower etiological specific survival than those with high PDIA3 expression in PTC [48]. These researches further enhanced the authenticity and credibility of our analytical strategies. Even if we did not explore these genes in depth, their pivotal roles in the progression of ATC could not be ignored.

Under physiological conditions, TLR4 was a critical regulator of inflammatory responses and tissue regeneration, thus maintaining tissue homeostasis [49]. Notably, accumulated evidence showed abnormal expression of TLR4 signaling in several tumors. Indeed, functional expression of TLRs in tumor cells was actively involved in tumor survival and progression, immune escape and apoptotic resistance [50]. Nevertheless, there have been inconsistences in impact of TLR4 activation on cancer cells across various cancer models. TLR4-induced secretion of IL-1β shaped the immunosuppressive microenvironment of pancreatic cancer [51]. MAPK/ERK signal-dependent TLR4 expression promoted the survival signaling and stemness, thus facilitating the PTC progression with BRAF-V600E mutations [52, 53]. Conversely, the research on pituitary epithelial neoplasms had found that activation of TLR4 impeded tumor cell growth [54]. Likewise, we found that the up-regulation of TLR4 was positively correlated with ICD score in ATC, and significantly correlated with the poor prognosis of thyroid cancer patients. It was curious about the ambiguous role of TLR4 in ATC progression.

In our study, we verified that both mRNA and protein levels of TLR4 were upregulated in ATC cell lines. After the silence of TLR4 expression, the ability of proliferation, metastasis, invasion and clone formation were significantly inhibited.

These results demonstrated that TLR4 was a crucial driver of ATC malignant phenotypes. Additionally, the WGCNA analysis revealed that the red module associated with TLR4 exhibited the strongest correlation with immune score, providing further insights into the potential molecular mechanisms of TLR4 in ATC. And enrichment analysis had found the genes related TLR4 involved in myeloid leukocyte activation, activation of immune response, leukocyte mediated immunity, T cell differentiation and microglial cell activation pathways. Nonetheless, more work was urgent to explain how TLR4 regulates these immune-related pathways in ATC. Intriguingly, we found that reducing TLR4 expression significantly enhanced PTX-induced proliferation inhibition and apoptosis, as well as increased CALR exposure and ATP release. These release of DAMPs were ICD-specific, although this was in contrast to the effect of TLR4 reported in the literature [13]. In addition, whether TLR4 regulated ICD based on immune or other signals required further validation.

In conclusion, we highlighted the first dynamic panorama of ICDGs and the immunological characteristics in ATC. These results might be facilitated for ICD-based immunotherapy interventions in ATC patients. We also identified and verified TLR4, the pivotal ICD-related gene driving the progression of ATC malignancy, which held great prospect and value for inducing ICD and improving ATC treatment.

Data Availability

The microarray datasets (GSE65144, GSE33630, GSE29265 and GSE76039) were gotten from the GEO database ( In addition, all data for this study are available from the corresponding author upon reasonable request.


  1. Lim H, Devesa SS, Sosa JA, Check D. and C. M. Kitahara. Trends in thyroid Cancer incidence and mortality in the United States, 1974–2013. JAMA. 2017;317(13):1338–48.

  2. Molinaro E, Romei C, Biagini A, Sabini E, Agate L, et al. Anaplastic thyroid carcinoma: from clinicopathology to genetics and advanced therapies. Nat Rev Endocrinol. 2017;13(11):644–60.

    Article  CAS  PubMed  Google Scholar 

  3. Akaishi J, Sugino K, Kitagawa W, Nagahama M, Kameyama K, et al. Prognostic factors and treatment outcomes of 100 cases of anaplastic thyroid carcinoma. Thyroid. 2011;21(11):1183–9.

    Article  PubMed  Google Scholar 

  4. Glaser SM, Mandish SF, Gill BS, Balasubramani GK, Clump DA, et al. Anaplastic thyroid cancer: prognostic factors, patterns of care, and overall survival. Head Neck. 2016;38(Suppl 1). (E2083-90.

  5. Maniakas A, Dadu R, Busaidy NL, Wang JR, Ferrarotto R, et al. Evaluation of overall survival in patients with anaplastic thyroid carcinoma, 2000–2019. JAMA Oncol. 2020;6(9):1397–404.

    Article  PubMed  Google Scholar 

  6. Jannin A, Escande A, Al Ghuzlan A, Blanchard P, Hartl D, et al. Anaplastic thyroid carcinoma: an update. Cancers (Basel). 2022;14(4).

  7. Tuccilli C, Baldini E, Sorrenti S, Catania A, Antonelli A, et al. CTLA-4 and PD-1 ligand gene expression in epithelial thyroid cancers. Int J Endocrinol. 2018;2018:1742951.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Capdevila J, Wirth LJ, Ernst T, Ponce Aix S, Lin CC, et al. PD-1 blockade in anaplastic thyroid carcinoma. J Clin Oncol. 2020;38(23):2620–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Lu X, Bao L, Pan Z, Ge M. Immunotherapy for anaplastic thyroid carcinoma: the present and future. Zhejiang Da Xue Xue Bao Yi Xue Ban. 2021;50(6):675–84.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Chintakuntlawar AV, Yin J, Foote RL, Kasperbauer JL, Rivera M, et al. A phase 2 study of Pembrolizumab Combined with Chemoradiotherapy as initial treatment for anaplastic thyroid Cancer. Thyroid. 2019;29(11):1615–22.

    Article  CAS  PubMed  Google Scholar 

  11. Kroemer G, Galassi C, Zitvogel L, Galluzzi L. Immunogenic cell stress and death. Nat Immunol. 2022;23(4):487–500.

    Article  CAS  PubMed  Google Scholar 

  12. Zhu M, Yang M, Zhang J, Yin Y, Fan X, et al. Immunogenic cell death induction by Ionizing Radiation. Front Immunol. 2021;12(705361).

  13. Lau TS, Chan LKY, Man GCW, Wong CH, Lee JHS, et al. Paclitaxel induces immunogenic cell death in Ovarian Cancer via TLR4/IKK2/SNARE-Dependent exocytosis. Cancer Immunol Res. 2020;8(8):1099–111.

    Article  CAS  PubMed  Google Scholar 

  14. Limagne E, Thibaudin M, Nuttin L, Spill A, Derangère V, et al. Trifluridine/Tipiracil plus Oxaliplatin improves PD-1 blockade in Colorectal Cancer by Inducing Immunogenic Cell Death and Depleting macrophages. Cancer Immunol Res. 2019;7(12):1958–69.

    Article  CAS  PubMed  Google Scholar 

  15. Li X, Zheng J, Chen S, Meng FD, Ning J, et al. Oleandrin, a cardiac glycoside, induces immunogenic cell death via the PERK/elF2α/ATF4/CHOP pathway in Breast cancer. Cell Death Dis. 2021;12(4):314.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Edgar R, Domrachev M, Lash AE. Gene expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 2002;30(1):207–10.

    Article  CAS  PubMed Central  Google Scholar 

  17. von Roemeling CA, Marlow LA, Pinkerton AB, Crist A, Miller J, et al. Aberrant lipid metabolism in anaplastic thyroid carcinoma reveals stearoyl CoA desaturase 1 as a novel therapeutic target. J Clin Endocrinol Metab. 2015;100(5):E697–709.

    Article  CAS  Google Scholar 

  18. Dom G, Tarabichi M, Unger K, Thomas G, Oczko-Wojciechowska M, 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Landa I, Ibrahimpasic T, Boucai L, Sinha R, Knauf JA, et al. Genomic and transcriptomic hallmarks of poorly differentiated and anaplastic thyroid cancers. J Clin Invest. 2016;126(3):1052–66.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Wen S, Qu N, Ma B, Wang X, Luo Y, et al. Cancer-Associated fibroblasts positively correlate with dedifferentiation and aggressiveness of thyroid Cancer. Onco Targets Ther. 2021;14:1205–17.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Cangelosi R, Goriely A. Component retention in principal component analysis with application to cDNA microarray data. Biol Direct. 2007;2.

  22. Li Z, Safo SE, Long Q. Incorporating biological information in sparse principal component analysis with application to genomic data. BMC Bioinformatics. 2017;18(1):332.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Garg AD, De Ruysscher D, Agostinis P. Immunological metagene signatures derived from immunogenic cancer cell death associate with improved survival of patients with lung, breast or ovarian malignancies: a large-scale meta-analysis. Oncoimmunology. 2016;5(2):e1069938.

    Article  PubMed  Google Scholar 

  24. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25(8):1091–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Song WM, Zhang B. Multiscale Embedded Gene Co-expression Network Analysis. PLoS Comput Biol. 2015;11(11):e1004574.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.

    Article  CAS  PubMed  Google Scholar 

  27. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013;39(4):782–95.

    Article  CAS  PubMed  Google Scholar 

  29. Wu TT, Chen YF, Hastie T, Sobel E, Lange K. Genome-wide association analysis by lasso penalized logistic regression. Bioinformatics. 2009;25(6):714–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Xu T, Jin T, Lu X, Pan Z, Tan Z, et al. A signature of circadian rhythm genes in driving anaplastic thyroid carcinoma malignant progression. Cell Signal. 2022;95(110332).

  32. Hegde PS, Chen DS. Top 10 challenges in Cancer Immunotherapy. Immunity. 2020;52(1):17–35.

    Article  CAS  PubMed  Google Scholar 

  33. Kennedy LB. S. Salama. A review of cancer immunotherapy toxicity. CA Cancer J Clin. 2020;70(2):86–104.

    Article  PubMed  Google Scholar 

  34. Ahmed A, Tait SWG. Targeting immunogenic cell death in cancer. Mol Oncol. 2020;14(12):2994–3006.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Fucikova J, Kepp O, Kasikova L, Petroni G, Yamazaki T, et al. Detection of immunogenic cell death and its relevance for cancer therapy. Cell Death Dis. 2020;11(11):1013.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Calvillo-Rodríguez KM, Mendoza-Reveles R, Gómez-Morales L, Uscanga-Palomeque AC, Karoyan P, et al. PKHB1, a thrombospondin-1 peptide mimic, induces anti-tumor effect through immunogenic cell death induction in Breast cancer cells. Oncoimmunology. 2022;11(1):2054305.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Hayashi K, Nikolos F, Lee YC, Jain A, Tsouko E, et al. Tipping the immunostimulatory and inhibitory DAMP balance to harness immunogenic cell death. Nat Commun. 2020;11(1):6299.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Li Y, Zhang H, Li Q, Zou P, Huang X, et al. CDK12/13 inhibition induces immunogenic cell death and enhances anti-PD-1 anticancer activity in Breast cancer. Cancer Lett. 2020;495:12–21.

    Article  CAS  PubMed  Google Scholar 

  39. Wang X, Wu S, Liu F, Ke D, Wang X, et al. An immunogenic cell death-related classification predicts prognosis and response to Immunotherapy in Head and Neck squamous cell carcinoma. Front Immunol. 2021;12:781466.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Chauvin JM, Zarour HM. TIGIT in cancer immunotherapy. J Immunother Cancer. 2020;8(2).

  41. Crespo J, Sun H, Welling TH, Tian Z, Zou W. T cell anergy, exhaustion, senescence, and stemness in the Tumor microenvironment. Curr Opin Immunol. 2013;25(2):214–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. McLane LM, Abdel-Hakeem MS, Wherry EJ. CD8 T cell exhaustion during chronic viral Infection and Cancer. Annu Rev Immunol. 2019;37:457–95.

    Article  CAS  PubMed  Google Scholar 

  43. Fan R, Dong L, Li P, Wang X, Chen X. Integrated bioinformatics analysis and screening of hub genes in papillary thyroid carcinoma. PLoS ONE. 2021;16(6):e0251962.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Luo JH, Zhu YH, Xiang C. Favorable function of Ectonucleoside triphosphate diphosphohydrolase 1 high expression in thyroid carcinoma. Hereditas. 2021;158(1):33.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Metovic J, Vignale C, Annaratone L, Osella-Abate S, Maletta F, et al. The Oncocytic variant of poorly differentiated thyroid carcinoma shows a specific Immune-related gene expression Profile. J Clin Endocrinol Metab. 2020;105(12).

  46. Chen Y, Wang K, Shang M, Zhao S, Zhang Z, et al. Exploration of DNA methylation-driven genes in papillary thyroid carcinoma based on the Cancer Genome Atlas. J Comput Biol. 2021;28(1):99–114.

    Article  CAS  PubMed  Google Scholar 

  47. Schulten HJ, Hussein D, Al-Adwani F, Karim S, Al-Maghrabi J, et al. Microarray expression profiling identifies genes, including cytokines, and biofunctions, as diapedesis, associated with a brain Metastasis from a papillary thyroid carcinoma. Am J Cancer Res. 2016;6(10):2140–61.

    CAS  PubMed  PubMed Central  Google Scholar 

  48. Kure S, Chiba T, Ebina A, Toda K, Jikuzono T, et al. Correlation between low expression of protein disulfide isomerase A3 and lymph node Metastasis in papillary thyroid carcinoma and poor prognosis: a clinicopathological study of 1,139 cases with long-term follow-up. Endocr J. 2022;69(3):273–81.

    Article  CAS  PubMed  Google Scholar 

  49. Rakoff-Nahoum S, Medzhitov R. Toll-like receptors and cancer. Nat Rev Cancer. 2009;9(1):57–63.

    Article  CAS  PubMed  Google Scholar 

  50. Ridnour LA, Cheng RY, Switzer CH, Heinecke JL, Ambs S, et al. Molecular pathways: toll-like receptors in the Tumor microenvironment–poor prognosis or new therapeutic opportunity. Clin Cancer Res. 2013;19(6):1340–6.

    Article  CAS  PubMed  Google Scholar 

  51. Das S, Shapiro B, Vucic EA, Vogt S, Bar-Sagi D. Tumor cell-derived IL1β promotes Desmoplasia and Immune suppression in Pancreatic Cancer. Cancer Res. 2020;80(5):1088–101.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Gao Y, Wang F, Zhang L, Kang M, Zhu L, et al. LINC00311 promotes cancer stem-like properties by targeting miR-330-5p/TLR4 pathway in human papillary thyroid cancer. Cancer Med. 2020;9(4):1515–28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Peyret V, Nazar M, Martín M, Quintar AA, Fernandez EA, et al. Functional toll-like receptor 4 overexpression in papillary thyroid Cancer by MAPK/ERK-Induced ETS1 transcriptional activity. Mol Cancer Res. 2018;16(5):833–45.

    Article  CAS  PubMed  Google Scholar 

  54. Tichomirowa M, Theodoropoulou M, Lohrer P, Schaaf L, Losa M, et al. Bacterial endotoxin (lipopolysaccharide) stimulates interleukin-6 production and inhibits growth of pituitary tumour cells expressing the toll-like receptor 4. J Neuroendocrinol. 2005;17(3):152–60.

    Article  CAS  PubMed  Google Scholar 

Download references


Not applicable.


This research was funded by National Natural Science Foundation of China under Grant No. 82203858 and 82173157; Medical and Health Science and Technology Project of Zhejiang (No. 2021KY491); Clinical Research Fund Project of Zhejiang Medical Association (No. 2020ZYC-A104); Basic Scientific Research Project of Basic Scientific Research Funds of Hangzhou Medical College (No. KYQN202126); Zhejiang Provincial Program for the Cultivation of High-level Health Talents (to Ping Huang and Zongfu Pan).

Author information

Authors and Affiliations



Planed and designed the research: Tong Xu and Ping Huang. Conducted experiments and data analysis: Chaozhuang Zhu and Zongfu Pan. Wrote, polish and revised the manuscript: Feifeng Song, Wanli Zhang and Mengnan Yuan.

Corresponding author

Correspondence to Ping Huang.

Ethics declarations

Ethics approval and consent to participate

The study was performed in accordance with relevant guidelines and regulations and approved by the Ethics Committee of Zhejiang Provincial People’s Hospital. Informed consent was obtained from all subjects and their legal guardian.

Consent for publication

Not applicable.

Competing interests

All authors contributed to the research and approved the submitted version, and declared no known competing interests or personal relationships in this paper that could have appeared to influence the work reported.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

: The row bands of western blot

Rights and permissions

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

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Xu, T., Zhu, C., Song, F. et al. Immunological characteristics of immunogenic cell death genes and malignant progression driving roles of TLR4 in anaplastic thyroid carcinoma. BMC Cancer 23, 1131 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: