A Novel Pyroptosis-Related Gene Signature for Prognostic Prediction of Head and Neck Squamous Cell Carcinoma

Introduction

HNSCC marks the 6th most commonly seen malignant cancer across the globe, which originates from the nasal cavity, larynx, hypopharynx, oropharynx, nasopharynx, oral cavity, and salivary gland.^1,2 Owing to the absence of effective methods in early diagnosis, most of the HNSCC patients are already at advanced stages when diagnosed, and will suffer from recurrence after treatment, with a poor survival rate of 35%.³ The multimodal strategies, such as surgery, radiotherapy and chemotherapy, are the main treatments for HNSCC. In spite of the recent progress in therapies, the improvement of survival rate is still slow.⁴ Therefore, novel therapeutic targets are urgently required to bring better results for HNSCC patients, and it is necessary to develop dependable new prognostic models to bring more feasibility for targeted therapies.

Pyroptosis is a cell death triggered by caspase-1/4/5/11 and results in cell swelling, chromatin fragmentation, plasma membrane cleavage, and intracellular proinflammatory release.⁵ More and more studies reveal that pyroptosis is vital for cancer progression. Some previous studies have reported that some key components of pyroptosis, such as inflammatory vesicles, proinflammatory cytokines, and gasdermin proteins are closely related to tumourigenesis, invasion, and metastasis.^6,7 A few latest studies have confirmed that pyroptosis plays a critical role in tumour immune microenvironment and antitumour processes as well.^8,9 However, it remains unclear whether pyroptosis has an effect on the development of HNSCC.

In this research, we conducted a systematic research to identify the pyroptosis-associated gene expression in healthy tissues and HNSCC tissues, and the DEGs were adopted to establish a model of prognosis so as to reveal the prognostic role of those genes. The predictive ability of the model was validated by GEO data, and the expression of mRNA of the pyroptosis-associated genes in our model was validated by RT-qPCR analysis. Eventually, our research investigated the potential mechanisms by functional enrichment method.

Materials and Methods

Data Collection

The clinical information and the mRNA expression of 528 HNSCC sufferers were obtained from the TCGA database (https://portal.gdc.cancer.gov/), and another 96 HNSCC sufferers from the GEO cohort (GSE31056) were used as the verification set. Our research was conducted without ethical examination, given that the entire data were accessible openly. Our team chose an overall 33 pyroptosis-associated genes as per past researches.^7,10

Construction and Validation of a Prognostic Pyroptosis-Related Gene Signature

Our team compared the levels of pyroptosis-related gene expression of cancer and neighboring healthy samples to determine the DEGs via the “limma” R package.¹¹ The DEGs standard registered a false discovery rate (FDR) below 0.05. The DEGs are as described below: *P < 0.05, **P < 0.01, and ***P < 0.001. A PPI network regarding the DEGs was established via STRING 11.0 (https://string-db.org/). The Cox study facilitated the capability examination of pyroptosis-associated genes to forecast the OS. Benjamini-Hochberg modulated p-scores were adopted to diminish FDR. The “glmnet” R package was applied to implement a LASSO Cox regression in order to determine whether DEGs could forecast the OS and HNSCC suffer conditions.¹² The best score of the penalty parameter (λ), correlating with the lowest partial possible abnormality, was determined by cross-verification at 10 folds. Our team computed the risk score by the formula below: risk score = esum (the standardized expression level of every gene× the regression coefficient). Afterwards, our team applied the mid-value risk score as the criterion for categorizing sufferers into the two above mentioned groups. Out team described the gene distribution in these groups via implementing PCA by the “stats” R package.¹³ The “surv_cutpoint” function of the “survminer” R package was adopted to analyze the survival in order to determine the ideal cut-off values. Subsequently, our team leveraged the “survivalROC” R package to carry out time-dependent ROC study to assess whether the genetic feature could forecast or not.¹⁴ Afterwards, our research applied the univariable and multivariable Cox regression method to identify whether the risk score enabled the independent forecast of the prognosis on OS. For the validation, a HNSCC cohort from the GEO (GSE31056) was included in the analysis. The risk score was computed via the identical formula in the TCGA cohort. The sufferers in the GSE31056 cohort were also classified into two groups likewise on account of the mid-value risk score of TCGA, and these groups were utilized for verifying the gene signature.

Functional Enrichment Analysis

Our team finished the GO enrichment and the KEGG pathway analyses for the DEGs in these groups with the “clusterProfiler” R package and |log2 (fold-change)| ≥ 1 and FDR < 0.05 were deemed as the standards of these DEGs.¹⁵ The “gsva” R package was adopted on the ssGSEA (single-sample gene set enrichment analysis) to compute the clustering value of immune cells and pathways associated with immunity.¹⁶

Patients and Specimens

Our team chose 15 HNSCC samples and 10 healthy samples from the First Hospital of Jiaxing between 2020 and 2021. This research was approved by the Ethics Committee of the First Hospital of Jiaxing and was conducted in accordance with the Declaration of Helsinki, with all individuals offering informed consents. Each sample was put in liquid nitrogen and reserved at −80°C directly.

Quantitative Real-Time PCR

After trizol extracted the entire RNA from the samples, it was turned into cDNA by reverse transcription. PCR was implemented via TB Green^® Premix Ex Taq™II kit (Takara, China). The reactive activity is demonstrated as follows: half a minute at 94°C, half a minute at 58°C, one minute at 72°C, forty times. GAPDH served as the interior control, and the comparative level of expression was identify via 2−ΔΔ CT.

Statistical Analyses

Student’s t-test was applied to contrast the gene expressions in cancer samples and neighboring healthy samples. Our team implemented the Mann–Whitney test to contrast the enrichment values of immune cells and the pathways associated with immunity in these groups. The Log rank test and Kaplan–Meier approach were utilized to describe the survival curves. Afterwards, our team carried out Cox regression study to reveal the OS-predicting factors. Our research used all statistical observation in R (Version 4.0.3). A two-tailed P value below 0.05 was important on statistics. The overall flow diagram of the study is shown in Figure 1 and graphical abstract.

Figure 1 Workflow diagram of data analysis.

Results

Identification of DEGs Between Normal and Tumor Tissues

We compared the expression levels of 33 pyroptosis-associated genes between 44 healthy and 502 cancer samples from the TCGA cohort, and we found 18 differentially expressed genes (DEGs) (P < 0.05). Among them, 17 genes (AIM2, CASP1, CASP5, CASP6, CASP8, GSDMB, GSDMD, GSDME, IL1B, NLRC4, NLRP1, NLRP6, NLRP7, NOD1, PLCG1, PYCARD, TNF) were upregulated while only 1 gene (ELANE) was downregulated in the tumor group. The heatmap presented the expression levels of these genes (Figure 2A). To further investigate the interactions of these genes, we used a protein–protein interaction (PPI) method, and the outcomes are presented in Figure 2B. Figure 2C shows the associations between these pyroptosis-associated genes (red: positive association; blue: negative association).

Figure 2 Expressions of 33 pyroptosis-related genes and interactions between them. (A) Heatmap of pyroptosis-related genes between tumor and normal tissues. P values were showed as: *P < 0.05; **P< 0.01; ***P < 0.001. (B) The PPI network showing the interactions of the pyroptosis-related genes. (C) Correlation network of pyroptosis-related genes.

Construction of a Prognostic Gene Signature in the TCGA Cohort

Our team conducted the univariate Cox regression observation using the DEGs to screen out the survival-related genes, and finally, 4 genes (GSDME, NLRP1, NLRP6, IL1B) were tightly related to OS (P < 0.05) (Figure 3A). The aforementioned 4 genes were used for the purpose of establishing a model of prognosis via LASSO Cox regression method. A 4-gene signature was then established in the light of the best score of λ (Figure 3B and C). Our research adopted the formula below to compute the risk score: e(0.8938 * the NLRP1 expression level + 1. 0402 * the GSDME expression level + 0.4383* the level of expression of NLRP6 +1. 0041 * the IL1B expression level). Our study applied the mid-value risk score as the criterion for separating cases into low-risk and high-risk groups (Figure 3D). The PCA study presented that there were two diverse distribution of these HNSCC cases (Figure 3E). People in the high-risk group displayed a greater death ratio versus the other group (Figure 3F). Likewise, the Kaplan–Meier curve revealed that the low-risk group exhibited a notably greater OS versus the high-risk group (Figure 3G, P < 0.001). Figure 3H unveils the OS forecast role of the risk score, with the region below the curve (AUC) being 0.645, 0.612, and 0.623 of 1, 2, 3 years, separately.

Figure 3 Construction of risk signature in TCGA cohort. (A) Screen out of OS-related genes with univariate cox regression analysis. (B) LASSO regression of 4 OS-related genes. (C) Cross-validation for tuning parameter selection in LASSO regression. (D) Distribution of HNSCC patients based on risk score. (E) PCA plot for HNSCCs based on risk score. (F) Survival status for each patient. (G) Kaplan–Meier curves for OS of patients in high- and low-risk groups. (H) The ROC curves demonstrated the predictive efficiency of the risk score.

Validation of the Risk Signature in the GEO Cohort

A total of 96 HNSCC sufferers from the Gene Expression Omnibus (GEO) cohort (GSE31056) were used as the verification set. Our team divided the patients into two groups as mentioned above on the foundation of the mid-value risk score of the TCGA cohort (Figure 4A). The PCA displayed acceptable separation between the two groups (Figure 4B). The individuals in the low-risk group displayed a lower mortality rate versus the other group (Figure 4C). Moreover, the Kaplan–Meier curve also demonstrated that the group with low risk exhibited a remarkably higher OS compared to the group with high risk (Figure 4D, P < 0.001). The time-dependent ROC analysis of the GEO cohort suggests that the forecast ability of the model is satisfactory, with the AUC being 0.744, 0.741, and 0.796 of 1, 2, 3 years, separately (Figure 4E).

Figure 4 Validation of risk model in GEO cohort. (A) Distribution of HNSCC patients in GEO cohort based on risk score in TCGA cohort. (B) PCA plot for HNSCCs. (C) Survival status for each patient. (D) Kaplan–Meier curves for OS of patients in high- and low-risk groups. (E) The ROC curves demonstrated the predictive efficiency of the risk score.

Independent Prognostic Value of the Risk Signature

Cox regression study was employed to verify whether the risk score could independently forecast the OS. As discovered via the univariable Cox regression study, the risk score presented a tight relevance with the OS in both TCGA and GEO cohorts (HR = 1.6629, 95% CI = 1.2390–2.2319, P = 0.0007 and HR = 5.7277, 95% CI = 1.8784–17.4645, P = 0.0022) (Figure 5A and C). Likewise, the multivariable Cox regression study revealed a remarkable correlation between the OS and risk score in both cohorts (HR = 1.6631, 95% CI = 1.2316–2.2459, P < 0.001 and HR = 4.2197, 95% CI = 1.3920–12.7916, P = 0.0109) (Figure 5B and D). Moreover, our team created a heatmap of risk score and the clinical characteristics in the TCGA cohort, and the results indicated that the survival status of patients was distinctly diverse between the two groups (Figure 5E).

Figure 5 Univariate and multivariate Cox regression analyses for risk score in TCGA and GEO cohorts. (A) Univariate analysis for TCGA cohort. (B) Multivariate analysis for TCGA cohort. (C) Univariate analysis for GEO cohort. (D) Multivariate analysis for GEO cohort. (E) The heatmap for the connections between clinicopathologic features and the risk groups (*P < 0.05).

Functional Analyses Based on the Risk Model

In order to reveal the biological roles and pathways associated with the risk score, we firstly used the “limma”R package to find out the DEGs between the two groups. The DEGs standard was a FDR below 0.05 and |log2FC | ≥ 1. A total of 1025 DEGs were detected between the two groups in the TCGA cohort. Subsequently, for these DEGs, our team carried out the GO enrichment and KEGG pathway assay. Figure 6A shows the top 10 MF, CC, and BP terms. The prime enriched Go terms were gated channel activity, postsynaptic membrane, integral and intrinsic components of synaptic membrane, and synaptic membrane. Figure 6B shows the main KEGG pathways, including Nicotine addiction, Primary immunodeficiency, Glutamatergic synapse, and cAMP signaling pathway.

Figure 6 Functional analysis based on the DEGs between the two-risk group in the TCGA cohort. (A) Barplot graph for GO enrichment. (B) Bubble graph for KEGG pathways.

Comparison of the Immune Activity

The immunity condition was calculated via ssGSEA by virtue of the enrichment score and its correlation with the risk score was studied. By comparison, the enrichment scores of TIL, Tfh, T helper cells, pDCs, Neutrophils, CD8+ T cells, B cells, and aDCs were remarkably diverse (adjusted P < 0.05, Figure 7A). These groups demonstrated a notable diversity on the values of T cell co−activation, T cell co−suppression, Inflammation promotion, HLA, Cytolytic activity, and Check point (adjusted P < 0.05, Figure 7B).

Figure 7 The comparison of the ssGSEA scores for the immune cells and immune pathways in the TCGA cohort. (A) The comparison of the enrichment scores of 16 types of immune cells between the high- and low-risk groups. (B) The comparison of the enrichment scores of 13 immune-related pathways between the high- and low-risk groups. Adjusted P values were showed as: *P < 0.05; **P < 0.01; ***P < 0.001. Abbreviation: ns, not significant.

Validation of the 4 Pyroptosis-Related Genes in HNSCC Tissues

Our team verified the levels of mRNA expression as to the 4 pyroptosis-associated genes in the HNSCC samples and healthy samples via RT-qPCR, and the outcomes revealed that GSDME (Figure 8A), NLRP1 (Figure 8B), NLRP6 (Figure 8C), and IL1B (Figure 8D) were all distinctly greater in the HNSCC samples versus the healthy samples (P < 0.001).

Figure 8 RT-qPCR detecting the mRNA expression levels of GSDME (A), NLRP1 (B), NLRP6 (C), IL1B (D) in the HNSCC and normal tissues.

Discussion

Pyroptosis is defined as a novel type of programmed cell death characterized by inflammatory response and is reported to play a dual role in cancer progress.¹⁷ Signaling pathways and inflammatory mediators in the process of pyroptosis are closely associated with the invasion, proliferation, and metastasis of cancer.^18,19 Nevertheless, the induction of tumor cell pyroptosis could be a novel treatment.²⁰ However, whether pyroptosis-related genes are related to the outcomes of HNSCC patients and how they act in HNSCC remain unknown.

In our research, we systematically analyzed the levels of mRNA expression of the 33 pyroptosis-associated genes in the HNSCC and healthy samples. The outcomes showed that more than half of the pyroptosis-related genes (55%, 18/33) were DEGs by comparing the cancer and neighboring healthy samples and four of them were associated with the OS, which hinted that pyroptosis participated in HNSCC and exhibited a potential forecast role of relevant genes. In order to explore more on the forecast power of those genes, we built a 4-gene signature through LASSO Cox regression analysis, and then the signature was validated to perform satisfactorily in a GEO cohort. We further validated the mRNA levels of these 4 genes in the HNSCC and healthy tissues via RT-qPCR analysis, and the outcomes coincided with those in the TCGA cohort. The functional analyses indicated a different immune status between the two groups.

Four pyroptosis-related genes (GSDME, NLRP1, NLRP6, IL1B), which functioned via different mechanisms, were used to establish a forecast model.^21,22 GSDME, also known as DFNA5, is from the gasdermin family which is the executor of pyroptosis. The activated caspase-3 could realize the cleavage of GSDME to release a N-terminal fragment and then promote pyroptosis.²³ GSDME was recognized as a tumor suppressor gene due to its inactivation and the methylation of DNA promoter region was found in many tumors, including breast cancer,²⁴ hepatocellular carcinoma²⁵ and so on, however, it was reported to be overexpressed in nonsmall-cell lung cancer²⁶ and esophageal squamous cell cancer.²⁷ Therefore, GSDME may play a specific role in different tumors. Intriguingly, GSDME seemed to be an oncogenic gene, as it was upregulated in tumor tissues. Due to the limited data and the conflicting results in different tumors, our results about GSDME provide a new direction for future research. NLRP1 is one of the NALP family proteins consisting of a pyrin domain, a CARD domain and an NBD-LRR domain which can induce pyroptosis.²⁸ Previous studies²⁹ reported that NLRP1 was regulated downward in colorectal cancer samples relative to healthy samples. Nevertheless, in our research, the expression of NLRP1 was distinctly higher in the HNSCC specimens versus the healthy ones, and was related to a shorter survival time. The relationship between NLRP1-mediated pyroptosis and HNSCC is still unknown and needs further study. NLRP6 is a new type of inflammasome which play critical roles in the development of cancer. The knockdown of NLRP6 encourages gastric tumor cells to proliferate, migrate, and invade.³⁰ However, our results demonstrated that higher NLRP6 expression was identified in the tumor samples versus the normal ones. Similarly, IL1B is another type of inflammasome which is released by caspase-1 and can induce pyroptosis. A recent study has presented that the IL 1B expression is lower in hepatocellular carcinoma samples.³¹ On the contrary, cervical cancer cells release more IL 1B than normal cervical cells.³² We found that IL 1B was upregulated in the tumor tissues and its overexpression could predict poor outcomes, which revealed that it was a cancer-promoting gene in our research. Further research should highlight how IL 1B acts in the process of pyroptosis and tumor development. In summary, all 4 genes (GSDME, NLRP1, NLRP6, IL1B) were regulated upward in the HNSCC sufferers and were associated with unfavourable clinical results. Nevertheless, whether those genes influenced the HNSCC prognosis by virtue of modulating pyroptosis needs more studies.

We analysed the DEGs between the two groups divided by risk scores and discovered that the DEGs primarily emerged in immune and inflammatory responses. The outcomes of GO and KEGG methods displayed that pyroptosis might play an important role in tumor immune microenvironment. The enrichment scores of the crucial anti-tumor immune cells and immune-associated pathways were lower in the high-risk group in contrast to the other one, suggesting that the unfavourable clinical outcomes in the high-risk group might be induced by the damaged anti-cancer immune activity.

At present, the underlying mechanism of pyroptosis is less studied, especially in HNSCC. Our study identified 4 genes which acted as the regulators of pyroptosis, and then our team carried out a preliminary research on the prognostic value of these pyroptosis-associated genes and offered a new direction for future investigation. Nevertheless, whether these regulators play the same role in pyroptosis-pathway in HNSCC is still elusive, and deserves further in-depth research. In addition, there were limitations in our model of prognosis due to the openly available data. It is pivotal to conduct further in vitro and in vivo researches to validate the clinical significance of this model.

Conclusion

This research found that pyroptosis was tightly associated with HNSCC, as more than half of the pyroptosis-associated genes were differently expressed by comparing the tumor and healthy samples. Furthermore, we constructed a novel forecast model by 4 pyroptosis-associated genes. The risk score according to the signature was able to independently forecast the OS of HNSCC in both TCGA and GEO cohorts. The potential mechanism regarding the correlation of between the pyroptosis-associated genes and HNSCC is still vague and deserves more in-depth researches.