Skip to main content

Identification and validation of a pyroptosis-related prognostic signature for thyroid cancer

Abstract

Background

Pyroptosis is a form of programmed cell death triggered by inflammasomes. However, the roles of pyroptosis-related genes in thyroid cancer (THCA) remain still unclear.

Objective

This study aimed to construct a pyroptosis-related signature that could effectively predict THCA prognosis and survival.

Methods

A LASSO Cox regression analysis was performed to build a prognostic model based on the expression profile of each pyroptosis-related gene. The predictive value of the prognostic model was validated in the internal cohort.

Results

A pyroptosis-related signature consisting of four genes was constructed to predict THCA prognosis and all patients were classified into high- and low-risk groups. Patients with a high-risk score had a poorer overall survival (OS) than those in the low-risk group. The area under the curve (AUC) of the receiver operator characteristic (ROC) curves assessed and verified the predictive performance of this signature. Multivariate analysis showed the risk score was an independent prognostic factor. Tumor immune cell infiltration and immune status were significantly higher in low-risk groups, which indicated a better response to immune checkpoint inhibitors (ICIs). Of the four pyroptosis-related genes in the prognostic signature, qRT-PCR detected three of them with significantly differential expression in THCA tissues.

Conclusion

In summary, our pyroptosis-related risk signature may have an effective predictive and prognostic capability in THCA. Our results provide a potential foundation for future studies of the relationship between pyroptosis and the immunotherapy response.

Introduction

Thyroid cancer (THCA) is the most common form of endocrine cancer worldwide and the number of THCA cases is icreasing [1]. Thyroid nodules are one of the most common clinical findings [2]. THCA can be divided into at least four subtypes, papillary thyroid carcinoma (PTC), follicular thyroid cancer (FTC), medullary thyroid cancer (MTC) and anaplastic thyroid cancer (ATC), based on the histological features of a THCA tumor. PTC is the most frequent histological subtype and accounts for more than 90% of all thyroid cancer cases [3]. The majority of PTC cases have a relatively better prognosis after surgery and 131I treatment compared to that of the other THCA subtypes [4]. However, cervical lymph node metastasis (LNM) is a potential factor that can lead to some patients suffering local recurrence and poor prognosis [5,6]. Therefore, it is important to construct novel prognostic models or find novel biomarkers which will make targeted therapies more feasible and improve the survival of patients with PTC.

Pyroptosis is a form of programmed cell death triggered by inflammasomes [7,8]. Pyroptosis has been found to be closely associated with some diseases like diabetic nephropathy and atherosclerosis. Some studies have found that pyroptosis is involved in the proliferation, invasion and metastasis of tumors. Pyroptosis results in cell swelling, plasma membrane lysis, chromatin fragmentation and the release of intracellular proinflammatory compounds. Pyroptosis is distinguished from other forms of programmed cell death morphologically although it shares certain similar characteristics with apoptosis. Generally, cells undergoing pyroptosis exhibit DNA damage and chromatin condensation during the early stage, followed by plasma membrane blebbing as well as caspase activation without losing cell membrane integrity [9]. Caspase-1 activation leads to the canonical inflammasome-induced pyroptosis pathway. Human caspase-4,5 and murine caspase-11 activation leads to the non-canonical inflammasome-induced pyroptosis pathway [10,11]. The pore-forming domain, the main executor of pyroptosis, found in these caspases is similar to the one found in the crystal structure of the human gasdermin (GSDM) superfamily (GSDMA, GSDMB, GSDMC, GSDMD, DFNA5 and DFNB59). Multiple studied have shown the abnormal expression of the GSDM family in human cancers, which implicates the potential roles in the tumorigenesis and the development. The association between pyroptosis and cancer is complicated. Pyroptosis appears to exert a dual function in cancer progression and treatment. Not only does pyroptosis result in the release of inflammatory factors which stimulate the transformation of normal cells into tumor cells, but it can promote tumor cell death. Pyroptosis plays different roles in many different types of cancer. It may be proved beneficial in preventing colorectal tumor development, and it inhibits tumor growth in hepatocellular carcinoma [12,13]. Recent studies have explored and identified novel pyroptosis-related signatures in some cancers. For example, a pyroptosis-related signature was constructed to predict patient prognosis and response to immunotherapy in gastric cancer [14]. Pyroptosis-related genes also play an important role in tumor immunity and can be used to predict the prognosis of ovarian cancer [15]. A prognostic signature for lung adenocarcinoma was built based on pyroptosis-related regulators [16]. However, the prognostic value of pyroptosis-related genes in THCA has not yet been elucidated.

Therefore, our study was aimed at developing a novel prognostic signature based on pyroptosis-related genes to systematically explore the relationship between the signature and clinicopathological features and overall survival (OS) in THCA patients. Furthermore, tumor immune microenvironment (TIME), mutation profile and the response to ICI treatment associated with the signature in THCA were further explored. The signature could predict the prognosis and immunotherapy response. In addition, this study provides a better understanding of the relationship between pyroptosis and immunotherapy response in THCA patients.

Materials and methods

Data collection

A flowchart was illustrated in Additional file 1: Figure S1 to show the research methodology. We downloaded 568 gene expression profiles (58 normal samples and 510 tumor samples) of THCA and OS clinical information from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Patients were randomly divided into a training set (n = 251) and a test set (n = 251) (Additional files 3, 4, 5). There were no significant differences in clinical variables between the two sets (Table 1). A total of thirty-three pyroptosis-related genes were obtained from prior reviews [17,18,19].

Table 1 The clinical characteristics in training, test and total sets

Differentially expressed genes (DEGs) identification

We identified DEGs in all tumor and normal samples using the “limma” package in R language (version 4.0.4). A p value of < 0.05 was set as the screening criterion. The DEGs were signed with * if p < 0.05, ** if p < 0.01 and *** if p < 0.001.

Construction of the protein–protein interaction (PPI) network

A PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://www.string-db.org/) to explore the interactions between these DEGs.

Consensus clustering of pyroptosis-related genes

THCA patients were clustered into different subgroups based on the pyroptosis-related DEGs using the “ConsensusClusterPlus” package in R [20].

Construction of pyroptosis-related prognostic signature

Patients with THCA were divided into a training set and a test set at a 1:1 ratio. The training set was used to identify prognostic pyroptosis-related genes and develop a prognostic risk signature. The predictive capability was validated in the test set and total set. A univariate Cox proportional hazard regression was employed to identify the pyroptosis-related genes with prognostic values of OS. To prevent omissions, a cut-off p value < 0.2 was set to identify prognostic variables. Subsequently, we used a least absolute shrinkage and selection operator (LASSO) penalized Cox proportional hazards regression to avoid overfitting and constructed the prognostic signature with the “glmnet” package [21]. The model was determined by penalty parameter (λ) with tenfold cross-validation following the minimum criteria (i.e. the value of λ corresponding to the lowest partial likelihood deviance). The risk scores of each THCA patient were calculated based on the gene expression level and its coefficient. The risk score was calculated as follows: risk score = sum (pyroptosis gene expression level × corresponding coefficient). Patients were classified into high- and low-risk groups according to the median risk score. Principal Component Analysis (PCA) and t-distributed Stochastic Neighbor Embedding (t-SNE) were implemented using the “stats” and “Rtsne” packages, respectively. To validate the predictive power, Kaplan–Meier survival curves were analyzed using the “survival” and “survminer” packages and the area under the curves (AUCs) were calculated with the “survivalROC” package [22].

Functional enrichment analysis

All samples were divided into high- and low-risk groups according to the prognostic signatures. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using the “clusterProfiler” package in R software according to the DEGs (|log2FC|≥ 1 and FDR < 0.05) between the high- and low-risk groups. Meanwhile, GSEA was performed in the Hallmark gene set “h.all.v7.4.symbols.gmt” to analyze the enriched biological pathways of key genes using GSEA 4.1.0. A NOM p-value of < 0.05 was considered statistically significant.

Estimation of tumor-infiltrating immune cells

The immunoscore of every patient was obtained from the ESTIMATE algorithm using the “estimate” package [23]. CIBERSORT is a deconvolution algorithm based on RNA-Seq data to estimate the composition ratio of immune cells [24]. We calculated the relative proportions of 21 types of infiltrating immune cells in all tumor samples based on THCA transcriptional profiles. A Wilcoxon rank-sum test was used to evaluate the difference in the level of immune cell infiltration in high- and low-risk groups.

Evaluation of immune status

Single-sample GSEA (ssGSEA) was used to calculate the scores of 16 infiltrating immune cells and the activity of 13 immune-related pathways in the high- and low-risk groups using the “GSVA” package of R [25]. We also compared the expression of the HLA gene between the high- and low-risk groups.

Mutation analysis

The mutation data of THCA patients were also obtained from the TCGA data portal (https://portal.gdc.cancer.gov/). The data were further analyzed using the “maftools” package [26]. We calculated the tumor mutation burden (TMB) score of every patient as follows: (total mutation ÷ total covered bases) × 10^6 [27].

Quantitative PCR

Sixty-five matched tumorous and non-tumorous tissue specimens of PTC were collected from the First Affiliated Hospital of China Medical University. The clinicopathological characteristics of 65 THCA patients from our hospital are displayed in Table 2. Total RNA was extracted from tissue samples using RNAiso (Takara, Dalian, China), then RNA was reverse transcribed into cDNA with the QuantiTect Reverse Transcription Kit (Takara, Shiga, Japan). Quantitative Real-Time PCR (qRT-PCR) analyses were performed with SYBR-Green (Takara, Shiga, Japan) to validate gene expression, and the level of GAPDH served as an internal control. The relative expression was calculated based on the comparative Ct (2−ΔΔCt) method [28]. The primers’ sequences are listed in Table 3.

Table 2 The clinicopathological features of THCA (N = 65)
Table 3 Premier sequences for qRT-PCR analysis

Results

Identification of DEGs in the TCGA cohort of THCA

Based on the P-value < 0.05, a total of 22 differentially expressed genes were identified from 33 pyroptosis-related genes. Among them, 15 were up-regulated (CASP1, CASP3, CASP5, CASP6, GSDMA, GSDMB, GSDMD, NOD1, NOD2, ELANE, NLRC4, PRKACA, GPX4, PYCARD, IL18) and 7 were down-regulated (NLRP6, IL6, TNF, PJVK, SCAF11, TIRAP and CASP9) in THCA tumors. The heatmap shows the RNA expression levels of these genes (Fig. 1A). In addition, we analyzed the correlation among pyroptosis-related genes: GPX4 was significantly negatively correlated with SCAF11 (Cor = − 0.67), while CASP1 and NOD2, and CASP4 were significantly positively correlated (Cor = 0.84) (Fig. 1B). A protein–protein interaction (PPI) network with the minimum required interaction score (highest confidence 0.9) was constructed to further explore the interactions among these pyroptosis-related genes (Fig. 1C).

Fig. 1
figure1

Expression of the pyroptosis-related genes in THCA. A The heatmap showed the expression levels of 33 pyroptosis-related genes in normal and tumor samples. *P < 0.05, **P < 0.01, ***P < 0.001. B Pearson correlation analysis of the 33 pyroptosis-related genes in THCA. C PPI network indicated the interactions of the pyroptosis-related genes

Consensus clustering of pyroptosis-related genes

In order to investigate the connections between the expression of the pyroptosis-related DEGs and THCA subtypes, we performed a consensus clustering analysis using the ConsensusCluserPlus package based on the 22 pyroptosis-related DEGs. The number of clusters was represented by the letter “k”. The k = 3 was identified with optimal clustering stability from k = 2 to 9. Finally, the THCA patients were clustered into three subtypes, namely, Cluster1 (n = 270), Cluster2 (n = 210) and Cluster3 (n = 22) (Fig. 2A). We also analyzed the gene-expression pattern between three subtypes of PCA. The results suggested that Cluster1, Cluster2 and Cluster3 could gather together (Fig. 2B). There were no statistical differences between the three clusters and OS (Fig. 2C). When the gene expression profile and clinicopathological features between the three subtypes were compared with a heatmap, no significant correlations were found between clinicopathological features in the three subtypes (Fig. 2D).

Fig. 2
figure2

Identification of consensus clusters based on the pyroptosis-related genes. A Consensus clustering matrix for k = 3. B Principal Component Analysis (PCA) of the RNA expression profile in TCGA cohort. C Heatmap and clinicopathologic features of the three clusters. D Kaplan–Meier curves of overall survival (OS) in three clusters

Construction of pyroptosis-related risk signature

We performed a univariate Cox regression analysis to identify the prognosis-related genes. The 12 genes that met a criterion with P-value < 0.2 were significantly associated with survival in the training set. Among them, 4 genes (GPX4, IL18, PRKACA and NOD1) were protective genes with HR < 1, while 8 genes (TNF, IL1B, IL6, GSDMC, NLRC4, ELANE, PJVK and GSDME) were detrimental genes with HR > 1 (Fig. 3A). To minimize overfitting, the set underwent LASSO Cox regression analysis, and 4 of the 12 genes were chosen to construct a risk signature based on the optimum λ value (Fig. 3B, C). The formula of the four-gene signature was as follows: Risk score = (− 0.059 × IL18) + (0.2427 × GSDMC) + (0.2389 × PJVK) + (− 0.023 × NOD1). Subsequently, all patients with THCA in the training set were classified into the high- and low-risk groups based on the median risk score. PCA and t-SNE indicated that patients with different risks were distributed into two directions (Fig. 4A, B). Patients with high-risk scores had a significantly poorer OS than patients in the low-risk group (p = 0.005) (Fig. 4C). The areas under the curve (AUC) of the risk signature were 0.886 at 1-year, 0.733 at 3-year and 0.844 at 5-year (Fig. 4D). In addition, we ranked the patients’ risk scores and analyzed their distributions in the training set (Fig. 4E). The survival status of THCA patients in the training set was presented in the dot plot (Fig. 4F). The heatmap showed the expression patterns of 4 prognostic genes between two risk groups (Fig. 4G).

Fig. 3
figure3

Univariate Cox regression analysis and LASSO analysis. A Forest plot showing the result of univariate Cox regression analysis of OS, 12 genes with p < 0.2. B Cross-validation for tuning parameter selection in the LASSO regression. C LASSO analysis of 12 prognostic pyroptosis-related genes

Fig. 4
figure4

Construction of risk signature in the training set. A PCA plot, B t-SNE analysis of TCGA cohort in training set. C Kaplan–Meier curves for OS of THCA patients in high- and low-risk groups. D Time-dependent ROC analysis. E The distribution of risk score, (F) survival status, and (G) the expression patterns of 4 pyroptosis-related genes in high- and low-risk groups

Furthermore, the predictive capability of risk signature was verified in the test set and total set. The risk scores of every patient were calculated and the patients were divided into high- and low-risk groups in two sets as previously described. PCA and t-SNE confirmed that the patients in the different subgroups were separated into two clusters (Figs. 5A, B, 6A, B). Kaplan–Meier survival curves indicated that the OS of high-risk patients was lower than that of the low-risk groups in the test group (p = 0.041) (Fig. 5C). The 1-year AUC was 0.578, the 3-year AUC was 0.715 and 5-year AUC was 0.767 (Fig. 5D). The distribution of the risk score, survival status and the expression of 4 pyroptosis-related genes in the test set are presented in Fig. 5E–G.

Fig. 5
figure5

Validating the 4 pyroptosis-related genes risk signature in the testing set. A PCA plot, B t-SNE analysis for THCA in testing set. C Kaplan–Meier curves for OS of THCA patients in high- and low-risk groups. D Time-dependent ROC analysis. E The distribution of risk score, F survival status, and (G) the expression patterns of 4 pyroptosis-related genes in high- and low-risk groups

Fig. 6
figure6

Validating the 4 pyroptosis-related genes risk signature in the total set. A PCA plot, (B) t-SNE analysis for THCA in total set. C Kaplan–Meier curves for OS of THCA patients in high- and low-risk groups. D Time-dependent ROC analysis. (E) The distribution of risk score, (F) survival status, and (G) the expression patterns of 4 pyroptosis-related genes in high- and low-risk groups

The results in the total set were similar to those in the training set and test set. The OS was significantly different between the two risk groups (p < 0.001) (Fig. 6C). The AUCs for 1-year, 3-year and 5-year were 0.779, 0.738 and 0.804 (Fig. 6D). The distribution of risk score, patients’ survival status and expression heatmap of 4 prognostic genes are also displayed in Fig. 6E–G.

Independent prognostic value of the risk signature

A univariate Cox regression analysis was performed to explore the relationship between clinicopathological variables and risk score on OS of THCA patients in the total set (Table 4). The risk score could serve as an independent prognostic factor for THCA in the total set referring to the results of the multivariable Cox regression analysis. Meanwhile, the ROC curve indicated that the AUC increased when combining the risk score with other clinicopathological features, which suggested that the risk score was an independent prognostic factor (Fig. 7A–C).

Table 4 Univariate and multivariate Cox regression analyses
Fig. 7
figure7

The time-dependent ROC to evaluate the prognostic power based on risk score and clinical factors in (A) 1-year, (B) 3-year, (C) 5-year

Risk signature and prognostic analysis of clinicopathological factors

We analyzed the relationships between the risk signature and clinicopathological factors. The results suggested that there were significant differences between the different ages and N stages. The risk score was significantly lower in patients with N1 stage and age below 65. Nevertheless, the risk signature was not correlated with gender, T stage, M stage and clinical stage (Fig. 8A). Additionally, survival analysis showed patients with high-risk scores were inclined to have a poorer OS in all subgroups (Additional file 2: Figure S2). The heatmap displayed the relationship between prognostic gene expression and clinical factors (Fig. 8B).

Fig. 8
figure8

The relationships between risk signature and clinic-pathologic parameters (A) (age, gender, T stage, M stage, N stage and clinical stage). B Heatmap showing the connections between clinicopathologic factors and high- and low-risk groups

Functional analyses based on the risk signature

To further elucidate the potential biological functions and pathways that are correlated with the risk score, GO enrichment and KEGG pathway analyses were performed according to the DEGs between the high- and low-risk groups. The results suggested that the DEGs were mainly enriched in immune response, cytokine–cytokine receptor interaction and chemokine signaling pathway (Fig. 9A, B).

Fig. 9
figure9

Functional enrichment analysis based on the DEGs. A Bar plot graph for GO enrichment. B Bubble graph for KEGG pathways. C Gene set enrichment analysis (GSEA) showed the significantly enriched hallmarks of tumor sets based on the risk signature in TCGA

Gene set enrichment analyses (GSEA)

The transcript messages of THCA patients classified by risk score into high- and low-risk subgroups were analyzed by GSEA. The results revealed that the majority of pyroptosis-related prognostic signature genes regulate the immune and malignant hallmarks of THCA tumors. Specifically, biological pathways such as allograft rejection, apoptosis, IL2-Stat5 signaling, IL6-Jak-Stat3 signaling, inflammatory response, P53 pathway signaling, and TNFA signaling via NF-κB were found to be enriched in the low-risk group (Fig. 9C).

Difference of the tumor-infiltrating immune cell populations between high- and low-risk groups

Immune cells played an important part in the tumor immune microenvironment (TIME) and we analyzed the difference of the tumor-infiltrating immune cell population in the high- and low-risk groups to evaluate the relationship between the prognostic signature and TIME. The CIBERSORT algorithm was utilized to calculate the relative proportion of 22 types of immune cells in each of the THCA patients. The ratios of naive B cells, activated CD4 memory T cells, resting dendritic cells, and resting mast cells were significantly higher in the low-risk group (Fig. 10A). Moreover, the correlation analyses between the risk score and degree of immunocyte infiltration indicated that activated NK cells (R = 0.23, P = 0.0025) and activated mast cells (R = 0.16, P = 0.039) were positively correlated with risk score. However, the risk score was negatively correlated with activated CD4 memory T-cells (R = − 0.25, P = 0.00086), resting dendritic cells (R = − 0.25, P = 0.001), and resting mast cells (R = − 0.19, P = 0.011) (Fig. 10B).

Fig. 10
figure10

The associations of tumor-infiltrating immune cells and risk scores and immunotherapy gene expression analysis. A The infiltrating levels of immune cells in high- and low-risk groups. B Correlation between the risk score and infiltration abundances of immune cells. C The gene expression of PD-1, PD-L1, PD-L2, CTLA4, TIGIT and TIM-3 in high- and low-risk groups

Numerous studies have shown that immunotherapy is emerging as a new hope in cancer treatment, and immune checkpoint proteins play important parts in the immune response. Therefore, we compared the expression levels of common immune checkpoint proteins in the high- and low-risk groups. The results indicated that low-risk patients had significantly higher expression levels of PD-1 (programmed cell death 1), PD-L1 (programmed cell death ligand 1), PD-L2 (programmed cell death ligand 2), CTLA-4 (cytotoxic T-lymphocyte-associated protein 4), TIGIT (T Cell Immunoreceptor with Ig and ITIM Domains) and TIM-3 (T-cell immunoglobulin and mucin-domain containing-3) (p < 0.01) (Fig. 10C). The results indicated that patients with low pyroptosis-related signature score might have a better opportunity for ICI treatment.

Comparing the immune status between the high- and low-risk groups

To further evaluate the relationship between the immune status and the risk score, we quantified the infiltrating scores of 16 immune cells and the activity of 13 immune-related pathways between the high- and low-risk groups in the TCGA cohort of THCA using single-sample gene set enrichment analysis (ssGSEA). The heatmap shows the immune status of 29 immune signature gene sets in the high- and low-risk groups (Fig. 11A). The scores of aDCs, DCs, iDCs, macrophages, mast cells, neutrophils, pDCs, Tfh, Th1 cells, Th2 cells, TIL and Treg were significantly different between the two subgroups (Fig. 11B). Except for the cytolytic activity pathway, the immune-related pathways had higher activity in the low-risk group than in the high-risk group (Fig. 11C). The low-risk group showed not only more immune activities, but significantly lower tumor purity (Fig. 11D). Additionally, we analyzed the expression levels of HLA related genes. The results suggested that the low-risk group had higher expression levels of HLA genes compared to those in the high-risk group (Fig. 11E).

Fig. 11
figure11

Comparison of the ssGSEA scores in high- and low-risk groups. A The immune status, (D) tumor purity and (E) the expression of HLA related genes in high- and low-risk groups. The boxplot showed the enrichment scores of (B) 16 immune cells and (C) 13 immune-related functions between high- and low-risk groups

A pyroptosis-related risk signature and mutation profile

Gene mutation is one of the significant factors in tumorigenesis and development. We evaluated the relationship between the signature and mutation profile in THCA patients. The top three mutated genes in THCA patients were BRAF, NRAS and HRAS. The most frequently mutated genes in the high- and low-risk groups are shown in Fig. 12A. The level of TMB was markedly higher in the high-risk group than that in the low-risk group (p = 0.0026) (Fig. 12B). Furthermore, we observed that TMB was associated with OS (p = 0.033) (Fig. 12C).

Fig. 12
figure12

The mutation profile and TMB in high- and low-risk groups. A Mutation profile of THCA patients in high- and low-risk groups. B The relationship between the risk signature and TMB. C The association between TMB and OS

The expression levels of four prognostic genes

We further validated the expression of the four prognostic genes (IL18, GSDMC, PJVK and NOD1) in 65 pairs of clinical samples from patients with PTC using qRT-PCR analysis according to the bioinformatics analysis results. The results of the qRT-PCR showed that the mRNA expression of NOD1 and IL18 were significantly higher in PTC tissues (p < 0.05). However, the expression of PJVK was decreased in PTC samples (p < 0.001) (Fig. 13), which was consistent with the results of bioinformatic analysis.

Fig. 13
figure13

The expression levels of PJVK (A), NOD1 (B) and IL18 (C) quantified using qRT-PCR analysis in 65 paired thyroid cancer tissues and no-tumorous samples

Discussion

THCA is the most common endocrine malignancy and PTC accounts for more than 85% of all thyroid cancer cases. Thyroid nodules are very common. Ultrasound guided fine needle aspiration (US-FNA) is the most commonly available way to evaluate thyroid nodules. However, rapid On-Site Evaluation (ROSE) is found to be the most helpful with small sized nodules or nodules that are more difficult to sample for less experienced radiologists [2]. Patients with PTC have a > 90% 10-year survival rate after proper surgical treatment and radiotherapy [29]. In radiotherapy, we must consider harmful effects of radiation for normal tissues surrounding tumor tissues. Accurate calculation of out-of-field dose to be critical for informing risk estimates, such as estimation of out-of-field dose variation using Markus ionization chamber detector [30]. Although the incidence of PTC has increased significantly over the past few decades, the presence of LNM can lead to locoregional recurrence and mortality. Therefore, exploring novel therapeutic targets of THCA is still a major challenging issue.

Recently, the possible beneficial effects of cancer therapies promoting pyroptosis have attracted considerable attention. Pyroptosis, an inflammatory form of programmed cell death, influences the proliferation, invasion and metastasis of tumor cells. It is a more recently identified pathway of programmed cell death that is stimulated by a range of microbial infections and non-infectious stimuli [31]. Pyroptosis is regulated via a caspase-1-dependent or caspase-1-independent mechanism [32]. Pyroptosis exerts a dual function in cancer progression and treatment mechanisms [16]. It plays a vital role in cellular lysis and release of pro-inflammatory cytokines when a host defends against infections [33]. Pyroptosis results in the release of intracellular proinflammatory contents and induces an inflammatory response leading to the death of adjacent healthy cells, which contributes to the development and progression of malignancies [9]. Additionally, pyroptosis can promote tumor cell death which makes pyrolysis a potential novel therapeutic target for cancer treatment. However, the potential role of pyroptosis-related genes in THCA remains unknown. Therefore, we aimed to discover potential diagnostic markers for targeted therapy of pyroptosis to improve the survival of patients with THCA as well as explore the prognostic and diagnostic value of pyroptosis. Our study suggested that using immunotherapy to induce pyroptosis may be an effective therapeutic direction to improve patient prognosis.

In our study, we analyzed the mRNA expression patterns of 33 pyroptosis-related genes in THCA samples and normal samples and 22 were differentially expressed. Among these genes, 15 were up-regulated and 7 were down-regulated. We identified three subgroups of THCA using consensus clustering analysis according to the expression of pyroptosis-related genes and no significant differences were found in the clinicopathological features. We then derived four prognostic risk signatures from these pyroptosis-related genes based on the univariate Cox regression analysis and LASSO Cox regression analysis. The prognostic value of the four prognostic-relevant risk signatures was evaluated in THCA patients and validated using a TCGA internal dataset. Interleukin (IL)-18, belonging to the IL-1 superfamily, is a proinflammatory and immune regulatory cytokine. IL-18 was originally identified as an interferon (IFN)-γ-inducing factor and involved in Th1 and Th2 responses in T cells, natural killer (NK) cells and macrophages [34]. IL-18 plays a dual role in cancer, as it promotes tumor development, progression and metastasis and it enhances anti-tumor immunity and reduces tumor growth in a matter depending on cancer progression [35]. IL-18 in combination with IL-12, through the activation of NK and cytotoxic T-cells, produced IFN-γ, which contributed to tumor immunity and had anti-tumor activity in different preclinical models. At present, IL-18 has been studied as a novel treatment approach and immune checkpoint therapy to significantly improved cancer treatment. Thus, IL-18 combined with immune-checkpoint therapy might be a potential treatment for early-stage tumors [36, 37]. The human gasdermin (GSDM) family (GSDMA, GSDMB, GSDMC, GSDMD, DFNA5 and DFNB59), which modulates multifunctional signal processes, can regulate cell pyroptosis [38]. The abnormal expression of the GSDM family in human cancers has been previously demonstrated, which implies their potential roles in tumorigenesis [39]. Multiple studies have shown that dysregulation of GSDMC expression is correlated with the biological processes of multiple cancers. Saeki et al. [38] found that GSDMC inhibition of tumor cell growth behaved like a potential tumor suppressor in the gastrointestinal epithelium. Watabe et al. [40] proved that GSDMC overexpression promoted tumor cells metastasis and proliferation in B16 melanoma cells. We found that high GSDMC expression was associated with poorer survival, which suggested that GSDMC might participate in tumor cell tumorigenesis and progression in thyroid cancer. PJVK, also called DFNB59, lacks the C-terminal domains. All members of the GSDM superfamily except for PJVK have a complete two-domain structure [41]. We regarded PJVK as a pyroptosis-related gene with a complete N-terminal domain and a similar pore-forming activity to other GSDMs [15]. PJVK has been found expressed in heart, brain and kidney, however, the regulatory roles of PJVK are not well understand [42]. At present, PJVK was demonstrated to be associated with hearing impairment in humans and is located on chromosome 2 [43]. We found high PJVK expression in tumor tissues from patients with a relatively poor prognosis. The nucleotide-binding oligomerization domain protein-1 (NOD1), one of the most important members of the NOD-like receptor (NLR) family, can induce pro-inflammatory responses and is involved in the apoptotic signaling pathway in some tumor cells [44]. Several studies have indicated that NOD1 plays an important role in the development and progression of gastric cancer, colon cancer, breast cancer and cervical cancer [45, 46]. NOD1 was downregulated in tumor samples in our study and its high expression predicted better survival, which suggested that the NOD1 may be a tumor suppressor gene. The risk score was calculated from the risk signatures and classified the patients into high- and low-risk groups. Kaplan–Meier survival curves showed that the OS of high-risk patients was lower than that of the low-risk groups. The AUC of the ROC curve showed the risk signature was efficient in predicting survival prognosis. Univariate and multivariate Cox regression analyses indicated that the risk score was not only an independent risk factor for prognosis, but could predict the clinical characteristics of THCA. The functional enrichment analyses indicated that immune-related pathways were significantly enriched in the low-risk groups. Therefore, we reasonably speculated that the cell pyroptosis could participate in the TIME. To elucidate the association between infiltrating immune cells and THCA, we estimated the infiltration of tumor immune cells between high- and low-risk groups of patients and found that the high-risk groups had higher proportions of naive B-cells naïve, activated CD4 memory T-cells, resting dendritic cells, and resting mast cells, while the low-risk groups’ scores were positively correlated with the proportions of activated NK cells and activated mast cells. Patients with low-risk scores had higher overall immune activity based on the ssGSEA analysis. The tumor purity was significantly enriched in the high-risk group, which suggested the lower infiltration of stromal and immune cells. Recently, cancer immunotherapy has gained wide acceptance as a potential therapeutic agent or an alternative to standard chemotherapy and has made great progress in the field of cancer therapy [47, 48]. Therefore, we explored the response of common immune checkpoints inhibitors and the expression of PD-1, PD-L1, PD-L2, CTLA4, TIGIT and TIM-3 increased significantly in low-risk patients. Our study suggested that patients with low-risk score had higher expression of common immune checkpoint molecules. According to above findings, we speculate that low-risk patients might have a better, more beneficial response from treatment with checkpoint inhibitors of PD-1, PD-L1, PD-L2, CTLA4, TIGIT and TIM-3. TMB was a significant independent predictor of the responses to immunotherapy in diverse cancers. In this study, the predictive value of the risk signature is independent of TMB. We found that the low-risk groups showed lower TMB.

However, at present, there are very few studies of pyroptosis in the thyroid cancer field, and investigation of the potential mechanisms may be meaningful in the future. The exploration of prognostic value of pyroptosis-related genes sets the stage for future mechanism research. Our study still has some limitations that need to be considered. All analyses were performed using a small TCGA cohort, a larger external validation cohort should be used to verify the predictive power of the signature in THCA, preferably validated using GEO datasets. Unfortunately, no survival information of THCA could be obtained from the GEO cohort. Furthermore, some basic experiments are also required to investigate the correlation between the model and the tumor microenvironment.

Conclusion

Taken together, we constructed a pyroptosis-related prognostic signature of genes that possessed predictive power based on a comprehensive analysis of RNA sequencing data and clinical data of THCA available in TCGA database. The signature was significantly associated with the tumor immunity. This study provided a better understanding of the relationship between pyroptosis and immunotherapy response in THCA patients. The pyroptosis-related signature could provide new possibilities to predict the prognosis and contribute to the development of new individualized therapeutic strategy for future studies of patients with THCA.

Availability of data and materials

Gene expression profiles, clinical information and mutation data of THCA in this study are available from the public database (TCGA, https://portal.gdc.cancer.gov/).

References

  1. 1.

    Cabanillas ME, McFadden DG, Durante C. Thyroid cancer. Lancet. 2016;388:2783–95. https://doi.org/10.1016/S0140-6736(16)30172-6.

    CAS  Article  PubMed  Google Scholar 

  2. 2.

    Aly AK, et al. Rapid on-site evaluation (ROSE) for fine needle aspiration of thyroid: is it helpful? SciMed J. 2021;3:1–7. https://doi.org/10.28991/SciMedJ-2021-0301-1.

    Article  Google Scholar 

  3. 3.

    Wang X, et al. Identification and validation of m(6)A RNA methylation regulators with clinical prognostic value in Papillary thyroid cancer. Cancer Cell Int. 2020;20:203. https://doi.org/10.1186/s12935-020-01283-y.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  4. 4.

    Xie Z, et al. Analysis of the expression and potential molecular mechanism of interleukin-1 receptor antagonist (IL1RN) in papillary thyroid cancer via bioinformatics methods. BMC Cancer. 2020;20:1143. https://doi.org/10.1186/s12885-020-07620-8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  5. 5.

    Kim SK, et al. Predictive factors of lymph node metastasis in follicular variant of papillary thyroid carcinoma. Ann Surg Oncol. 2017;24:2617–23. https://doi.org/10.1245/s10434-017-5912-5.

    Article  PubMed  Google Scholar 

  6. 6.

    Li P, et al. Downregulation of CDH16 in papillary thyroid cancer and its potential molecular mechanism analysed by qRT-PCR, TCGA and in silico analysis. Cancer Manag Res. 2019;11:10719–29. https://doi.org/10.2147/CMAR.S229631.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  7. 7.

    Fang Y, et al. Pyroptosis: a new frontier in cancer. Biomed Pharmacother. 2020;121: 109595. https://doi.org/10.1016/j.biopha.2019.109595.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Kovacs SB, Miao EA. Gasdermins: effectors of pyroptosis. Trends Cell Biol. 2017;27:673–84. https://doi.org/10.1016/j.tcb.2017.05.005.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  9. 9.

    Ruan J, Wang S, Wang J. Mechanism and regulation of pyroptosis-mediated in cancer cell death. Chem Biol Interact. 2020;323: 109052. https://doi.org/10.1016/j.cbi.2020.109052.

    CAS  Article  PubMed  Google Scholar 

  10. 10.

    Vande Walle L, Lamkanfi M. Pyroptosis. Curr Biol. 2016;26:R568–72. https://doi.org/10.1016/j.cub.2016.02.019.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Xu YJ, Zheng L, Hu YW, Wang Q. Pyroptosis and its relationship to atherosclerosis. Clin Chim Acta. 2018;476:28–37. https://doi.org/10.1016/j.cca.2017.11.005.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Zaki MH, Vogel P, Body-Malapel M, Lamkanfi M, Kanneganti TD. IL-18 production downstream of the Nlrp3 inflammasome confers protection against colorectal tumor formation. J Immunol. 2010;185:4912–20. https://doi.org/10.4049/jimmunol.1002046.

    CAS  Article  PubMed  Google Scholar 

  13. 13.

    Ma X, et al. Loss of AIM2 expression promotes hepatocarcinoma progression through activation of mTOR-S6K1 pathway. Oncotarget. 2016;7:36185–97. https://doi.org/10.18632/oncotarget.9154.

    Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Shao W, et al. The pyroptosis-related signature predicts prognosis and indicates immune microenvironment infiltration in gastric cancer. Front Cell Dev Biol. 2021;9: 676485. https://doi.org/10.3389/fcell.2021.676485.

    Article  PubMed  PubMed Central  Google Scholar 

  15. 15.

    Ye Y, Dai Q, Qi H. A novel defined pyroptosis-related gene signature for predicting the prognosis of ovarian cancer. Cell Death Discov. 2021;7:71. https://doi.org/10.1038/s41420-021-00451-x.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  16. 16.

    Lin W, Chen Y, Wu B, Chen Y, Li Z. Identification of the pyroptosisrelated prognostic gene signature and the associated regulation axis in lung adenocarcinoma. Cell Death Discov. 2021;7:161. https://doi.org/10.1038/s41420-021-00557-2.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  17. 17.

    Xia X, et al. The role of pyroptosis in cancer: pro-cancer or pro-“host”? Cell Death Dis. 2019;10:650. https://doi.org/10.1038/s41419-019-1883-8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  18. 18.

    Karki R, Kanneganti TD. Diverging inflammasome signals in tumorigenesis and potential targeting. Nat Rev Cancer. 2019;19:197–214. https://doi.org/10.1038/s41568-019-0123-y.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  19. 19.

    Wang B, Yin Q. AIM2 inflammasome activation and regulation: a structural perspective. J Struct Biol. 2017;200:279–82. https://doi.org/10.1016/j.jsb.2017.08.001.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  20. 20.

    Wang W, Sun B, Xia Y, Sun S, He C. RNA N6-methyladenosine-related gene contribute to clinical prognostic impact on patients with liver cancer. Front Genet. 2020;11:306. https://doi.org/10.3389/fgene.2020.00306.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Wang H, Lengerich BJ, Aragam B, Xing EP. Precision Lasso: accounting for correlations and linear dependencies in high-dimensional genomic data. Bioinformatics. 2019;35:1181–7. https://doi.org/10.1093/bioinformatics/bty750.

    CAS  Article  PubMed  Google Scholar 

  22. 22.

    Lorent M, Giral M, Foucher Y. Net time-dependent ROC curves: a solution for evaluating the accuracy of a marker to predict disease-related mortality. Stat Med. 2014;33:2379–89. https://doi.org/10.1002/sim.6079.

    Article  PubMed  Google Scholar 

  23. 23.

    Yoshihara K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612. https://doi.org/10.1038/ncomms3612.

    CAS  Article  Google Scholar 

  24. 24.

    Newman AM, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12:453–7. https://doi.org/10.1038/nmeth.3337.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  25. 25.

    Liang JY, et al. A novel ferroptosis-related gene signature for overall survival prediction in patients with hepatocellular carcinoma. Int J Biol Sci. 2020;16:2430–41. https://doi.org/10.7150/ijbs.45050.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28:1747–56. https://doi.org/10.1101/gr.239244.118.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Robinson DR, et al. Integrative clinical genomics of metastatic cancer. Nature. 2017;548:297–303. https://doi.org/10.1038/nature23306.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  28. 28.

    Higuchi R, Fockler C, Dollinger G, Watson R. Kinetic PCR analysis: real-time monitoring of DNA amplification reactions. Biotechnology (N Y). 1993;11:1026–30. https://doi.org/10.1038/nbt0993-1026.

    CAS  Article  Google Scholar 

  29. 29.

    Leboulleux S, et al. Prognostic factors for persistent or recurrent disease of papillary thyroid carcinoma with neck lymph node metastases and/or tumor extension beyond the thyroid capsule at initial diagnosis. J Clin Endocrinol Metab. 2005;90:5723–9. https://doi.org/10.1210/jc.2005-0285.

    CAS  Article  PubMed  Google Scholar 

  30. 30.

    Abdelaal AM, Attalla EM, Elshemey WM. Estimation of out-of-field dose variation using Markus ionization chamber detector. SciMedicine Journal. 2020;2:8–15. https://doi.org/10.28991/SciMedJ-2020-0201-2.

    Article  Google Scholar 

  31. 31.

    Bergsbaken T, Fink SL, Cookson BT. Pyroptosis: host cell death and inflammation. Nat Rev Microbiol. 2009;7:99–109. https://doi.org/10.1038/nrmicro2070.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  32. 32.

    Man SM, Karki R, Kanneganti TD. Molecular mechanisms and functions of pyroptosis, inflammatory caspases and inflammasomes in infectious diseases. Immunol Rev. 2017;277:61–75. https://doi.org/10.1111/imr.12534.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  33. 33.

    Guo H, Xie M, Zhou C, Zheng M. The relevance of pyroptosis in the pathogenesis of liver diseases. Life Sci. 2019;223:69–73. https://doi.org/10.1016/j.lfs.2019.02.060.

    CAS  Article  PubMed  Google Scholar 

  34. 34.

    Li Z, Yu X, Werner J, Bazhin AV, D’Haese JG. The role of interleukin-18 in pancreatitis and pancreatic cancer. Cytokine Growth Factor Rev. 2019;50:1–12. https://doi.org/10.1016/j.cytogfr.2019.11.001.

    CAS  Article  PubMed  Google Scholar 

  35. 35.

    El-Deeb MMK, El-Sheredy HG, Mohammed AF. The possible role of interleukin (IL)-18 and nitrous oxide and their relation to oxidative stress in the development and progression of breast cancer. Asian Pac J Cancer Prev. 2019;20:2659–65. https://doi.org/10.31557/APJCP.2019.20.9.2659.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Ma Z, et al. Augmentation of immune checkpoint cancer immunotherapy with IL18. Clin Cancer Res. 2016;22:2969–80. https://doi.org/10.1158/1078-0432.CCR-15-1655.

    CAS  Article  PubMed  Google Scholar 

  37. 37.

    Yasuda K, Nakanishi K, Tsutsui H. Interleukin-18 in health and disease. Int J Mol Sci. 2019;20:649. https://doi.org/10.3390/ijms20030649.

    CAS  Article  PubMed Central  Google Scholar 

  38. 38.

    Saeki N, et al. Distinctive expression and function of four GSDM family genes (GSDMA-D) in normal and malignant upper gastrointestinal epithelium. Genes Chromosomes Cancer. 2009;48:261–71. https://doi.org/10.1002/gcc.20636.

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    Lutkowska A, et al. Analysis of rs8067378 polymorphism in the risk of uterine cervical cancer from a polish population and its impact on gasdermin B expression. Mol Diagn Ther. 2017;21:199–207. https://doi.org/10.1007/s40291-017-0256-1.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Watabe K, et al. Structure, expression and chromosome mapping of MLZE, a novel gene which is preferentially expressed in metastatic melanoma cells. Jpn J Cancer Res. 2001;92:140–51. https://doi.org/10.1111/j.1349-7006.2001.tb01076.x.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Shi P, et al. Loss of conserved Gsdma3 self-regulation causes autophagy and cell death. Biochem J. 2015;468:325–36. https://doi.org/10.1042/BJ20150204.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Van Laer L, et al. Nonsyndromic hearing impairment is associated with a mutation in DFNA5. Nat Genet. 1998;20:194–7. https://doi.org/10.1038/2503.

    Article  PubMed  Google Scholar 

  43. 43.

    Liu X, et al. Gasdermins: pore-forming activities and beyond. Acta Biochim Biophys Sin. 2020;52:467–74. https://doi.org/10.1093/abbs/gmaa016.

    CAS  Article  PubMed  Google Scholar 

  44. 44.

    Liu X, et al. HPV-mediated down-regulation of NOD1 inhibits apoptosis in cervical cancer. Infect Agent Cancer. 2020;15:6. https://doi.org/10.1186/s13027-020-0272-3.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  45. 45.

    Velloso FJ, Trombetta-Lima M, Anschau V, Sogayar MC, Correa RG. NOD-like receptors: major players (and targets) in the interface between innate immunity and cancer. Biosci Rep. 2019;39. 10.1042/BSR20181709.

  46. 46.

    Jiang HY, et al. Activation of the pattern recognition receptor NOD1 augments colon cancer metastasis. Protein Cell. 2020;11:187–201. https://doi.org/10.1007/s13238-019-00687-5.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  47. 47.

    Liu J, et al. Identification of a prognostic signature of epithelial ovarian cancer based on tumor immune microenvironment exploration. Genomics. 2020;112:4827–41. https://doi.org/10.1016/j.ygeno.2020.08.027.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Charoentong P, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 2017;18:248–62. https://doi.org/10.1016/j.celrep.2016.12.019.

    CAS  Article  PubMed  Google Scholar 

Download references

Acknowledgements

We would like to thank TopEdit (www.topeditsci.com) for English language editing of this manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (81902726), China Postdoctoral Science Foundation (2018M641739).

Author information

Affiliations

Authors

Contributions

PW performed the data analysis and wrote the manuscript. JS performed the experiments. WS and HZ reviewed and revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Hao Zhang.

Ethics declarations

Ethics approval and consent to participate

The research was approved by the Institutional Research Ethics Committees of the First Affiliated Hospital of China Medical University. Informed consent for publication was obtained from all patients for collection of tissue samples prior to the surgery.

Consent for publication

Not applicable.

Competing interests

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Additional information

Publisher's Note

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

Supplementary Information

Additional file 1

: Figure S1: The flowchart of the study.

Additional file 2

: Figure S2: The stratification analysis of clinical factors based on age, gender, TNM stage and clinical stage.

Additional file 3

. The data of the total set.

Additional file 4

. The data of training set.

Additional file 5

. The data of test set.

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 http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) 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

Verify currency and authenticity via CrossMark

Cite this article

Wu, P., Shi, J., Sun, W. et al. Identification and validation of a pyroptosis-related prognostic signature for thyroid cancer. Cancer Cell Int 21, 523 (2021). https://doi.org/10.1186/s12935-021-02231-0

Download citation

Keywords

  • THCA
  • Pyroptosis
  • Prognostic signature
  • TMB
  • ssGSEA