Pyroptosis-Related Signature and Tumor Microenvironment Infiltration Characterization in Head and Neck Squamous Cell Carcinoma

Background: Head and neck squamous cell carcinoma (HNSCC) is the sixth most widespread and deadly cancer. Until now, very few studies have systematically evaluated the role of pyroptosis-related genes (PRGs) and lncRNAs in HNSCC patients. Methods: We integrated the genomic data to comprehensively assess the role of pyroptosis with the tumor microenvironment cell-infiltrating characteristics in HNSCC. In addition, we also constructed a set of the scoring system to calculate the pyroptosis dysfunction in each patient. Results: The analysis of the CNV alteration frequency displayed that CNV changes were common in 33 PRGs, and the frequency of copy number gain and loss was similar. CASP8 demonstrated the highest mutation frequency. Considering the individual heterogeneity, a scoring system to quantify the pyroptosis pattern in each patient was constructed based on these phenotypic-related genes, which we named as the PyroptosisScore. The results indicated that the low PyroptosisScore group experienced increased extensive TMB than the high group, with the most significant mutated genes being TP53 and TTN. Finally, we tried to find some useful pyroptosis-related lncRNAs, and 14 differentially expressed lncRNAs were selected as independent prognosis factors of HNSCC patients based on the multivariate Cox analysis. Conclusion: This work suggests the pyroptosis features and the potential mechanisms of the tumor microenvironment. The exploration may assist in identifying novel biomarkers and help patients predict prognosis, clinical diagnosis, and management.


BACKGROUND
Head and neck squamous cell carcinoma (HNSCC) arises in the stratified mucosa that lines the upper aero-digestive tract such as the tongue, lips, tonsil, nasopharynx, and larynx (Peyrade et al., 2021). There are over 655,000 new HNSCC cases per year, with over half a million deaths and ranking as the sixth most widespread and deadly cancer worldwide (Ferlay et al., 2010). The important danger elements involve tobacco smoking, alcohol, and high-risk human papillomavirus infection (Geng et al., 2021). Epstein-Barr virus infection is closely linked to nasopharyngeal carcinoma (Hitt et al., 1989). Although many forms of treatment, such as surgery, radiation, and chemotherapy, the asymptomatic nature and the absence of timely detection result in a 5-year survival rate of less than 50% in cancer patients (Bishop et al., 2021). Moreover, 50-60% of patients with over-stage III disease relapse locoregionally (Jacobs et al., 1992). Hence, an improved understanding of preclinical models is critical for determining prognosis between tumors and their host microenvironment and thus developing an appropriate therapeutic method.
Pyroptosis, also known as cell inflammatory necrosis, is a kind of programmed cell death, which is manifested by the continuous expansion of cells until the cell membrane ruptures, which causes the release of cell contents and activates a strong inflammatory response (Rühl et al., 2018). The body can recruit and activate caspase-1 protein through the classical inflammasome pathways (canonical inflammasome pathways) or activate the caspase 4/5/ 11 protein through the non-canonical inflammasome pathway to directly cleave and activate gasdermin protein D (GSDMD), eventually resulting in the formation of membrane pores and cell death (Jorgensen et al., 2017), while the gasdermin protein is the sole executor of cell pyroptosis. The GSDM protein family includes six members (GSDMD, GSDMA, GSDMB, GSDMC, GSDME (DFNA5), and DFNB59), which can be cleaved by cysteine protease caspase through a variety of ways, exposing its active N-terminus, and then, pores are formed on the cell membrane surface, causing cell swelling, thereby inducing cell pyroptosis and ultimately participating in the occurrence and development of inflammation and malignant tumors (Qiu et al., 2017). Lots of evidence shows pyroptosis in tumors is on the ascendancy. The study indicated that granzyme B can straightly cleave GSDME and initiate pyroptosis, promote activating the antitumor immune response, and suppress tumor growth . Exogenous calcium treatment can result in a significant increase in the caspase-1 activity, and the N-terminal effector domain of GSDMD was exceedingly increased in the CD38-overexpressing HNSCC . The study also found that the suppression of caspase-3mediated GSDME-derived apoptosis contributed to noncancerous tissue defense in patients with oral squamous cell carcinoma among cisplatin chemotherapy (Huang et al., 2020). Interestingly, the study indicated knocking down LncRNA-XIST increased ROS levels and activated NLRP3 inflammasome in lung cells (Liu et al., 2019), and the generation of ROS was reported to induce pyroptosis (Cui et al., 2019). However, until now very few studies have systematically evaluated the role of pyroptosis-related genes (PRGs) and lncRNAs in HNSCC patients. In the current study, we integrated the genomic data to comprehensively assess the role of pyroptosis with the tumor microenvironment (TME) cell-infiltrating characteristics in HNSCC. In addition, we also constructed a set of the scoring system to calculate the pyroptosis dysfunction in each patient.

Data Collection
RNA-sequence (487 tumors and 42 normal, FPKM value) data based on 511 patients and the somatic mutation data were collected from The Cancer Genome Atlas (TCGA-HNSCC) database. Meanwhile, the GSE65858 cohorts were downloaded from Gene Expression Omnibus (GEO) in this study for further analysis. The clinical information of the patients is demonstrated in Table 1. The FPKM value is converted to transcripts per kilobase million (TPM) value. The batch effect from the nonbiotech bias is corrected through the "combat" algorithm based on the SVA R package. The copy number variation (CNV) and long non-coding RNA (lncRNA) data were acquired from the TCGA database. We extracted 33 PRGs from prior reviews (Karki and Kanneganti, 2019;Xia et al., 2019;Wang and Yin, 2017;Man and Kanneganti, 2015), which are shown in Supplementary  Table S1. GTF files were got from Ensembl for annotation to identify the lncRNAs and mRNAs. The connection was regarded as significant if the correlation coefficient |R (Ferlay et al., 2010)| >0.3 and p < 0.001 between mRNAs and lncRNAs. The differential expression of lncRNAs or mRNAs was decided at FDR<0.05 and |log 2 FC|≥1 using R software, limma package.

Unsupervised Clustering for 33 PRGs
Based on the expressions of 33 PRGs, unsupervised cluster analysis was performed to identify different pyroptosis modification patterns, and patients were classified for further analysis. The number and stability of clustering are determined by using a consensus clustering algorithm (Hartigan, 1979). ConsensusClusterPlus R package was used to process the aforementioned analysis, and 1000 repeats were performed to ensure the stability of the classification (Wilkerson and Hayes, 2010). In order to determine the differences in biological processes between pyroptosis modification models, we performed gene set variation analysis (GSVA) enrichment analysis based on the "GSVA" R package. The gene set of "c2.cp.kegg. V6.4. symbols" was obtained from the MSIGDB online tool for GSVA analysis. GSVA is a nonparametric, unsupervised method for estimating genome-enrichment variation from samples of expressed datasets.

Generation of the PRG Signature
To quantify the pyroptosis modification pattern in single cancer, a scoring system was constructed to evaluate the pyroptosis gene marker in a single HNSCC patient, which we named PyroptosisScore. The steps to construct the signature of the pyroptosis gene signature are as follows : first, the differentially expressed genes (DEGs) determined from individual PyroptosisCluster were standardized in all HNSCC samples, and the overlapping genes were collected. By using unsupervised cluster analysis to overlap the DEGs, the patients were divided into several categories for subsequent analysis. The number and stability of gene clusters are determined by the consensus clustering algorithm. Second, univariate Cox regression analysis determined the prognostic gene in the signature. Finally, principal component analysis (PCA) was conducted to establish the PRG signature and choose the main components 1 and 2 to serve as signature scores. Thus, the PyroptosisScore was acquired based on the formula (Zeng et al., 2019): PyroptosisScore = (PC1i + PC2i). Here, i means the expression value of phenotype PRGs.

Estimation of TME Cell Infiltration
In order to quantify the subsets of tumor-infiltrating immune cells among the two groups and to assess their immune function, a single-sample gene-set enrichment analysis (ssGSEA) was applied. Enrichment scores obtained by ssGSEA analysis mirror the relative abundance of individual TME-infiltrated cells in every subject. Spearman correlation analysis was applied to analyze the connection between the risk score values and the TIICs (CIBERSORT, TIMER, XCELL, QUANTISEQ, EPIC, and MCPcounter algorithms). The correlation coefficients of the outcomes were demonstrated in a lollipop diagram by using R ggplot2 packages.

Construction of the Pyroptosis-Related lncRNA Signature Based on TCGA
The prognostic value of pyroptosis-related lncRNAs determined by univariate and multivariate Cox regression analyses and the significant lncRNA were selected as characteristic lncRNAs. The prognostic signature as risk score = e sum (lncRNA's expression×coefficient) . The risk score of individual patients was also counted. Based on the median risk score value, the set was divided into two (high-risk and low-risk) groups, and ROC curves were used to predict the accuracy of prognostic signatures. The pyroptosis-related lncRNA signature on patient survival was analyzed by applying the Kaplan-Meier curve.

STATISTICAL ANALYSIS
The normally distributed variables and the non-normally distributed variables were performed by the unpaired student's t-test and the Wilcoxon test, respectively. The signature and clinicopathological connections were assessed based on the chi-squared test. Logistic regression analyses to determine if the signature can be used as an independent clinical prognostic predictor. The "surv-cutpoint" function repeatedly tested all possible cut points to find the maximum rank statistic, which was classified as PyroptosisScore and then divided the patients into the high-PyroptosisScore group and the low-PyroptosisScore group based on the maximum selected logrank statistic to reduce the batch effect of the calculation. The CNV landscape of 33 PRGs in 23 pairs of chromosomes used the R package of RCircos. The analysis of mutations in patients with high-and low-PyroptosisScore subtypes were based on the waterfall function of the maftools R package in the TCGA-HNSCC cohort. For each analysis, the statistically significant setting was p < 0.05.

Landscape of the Genetic Variation of 33 PRGs
The summarized occurrence of CNV and somatic mutations of 33 PRGs in TCGA-HNSCC is shown in Figure 1A. Among 506 samples, 127 experienced mutations with a frequency of 25.1% and CASP8 demonstrated the highest mutation frequency, and most samples experienced nonsense mutation. The analysis of the CNV alteration frequency displayed that CNV changes were common in 33 PRGs, and the frequency of copy number gain and loss was similar, while GSDMC and GSDMD had a general frequency of CNV gain, and IL18, ELANE, and GPX4 had a widespread frequency of CNV loss ( Figure 1B). Meanwhile, the location of the CNV alteration of these genes on chromosomes is shown in Figure 1C. The expression of these genes between normal and tumor cells in TCGA-HNSCC is shown in Figure 1D. In general, PRGs with amplified CNV were significantly higher expressed in tumor tissues and vice versa. The aforementioned analysis revealed a high level of heterogeneity in the genetic and expression patterns of the pyroptosis in HNSCC, indicating that pyroptosis plays a key role in the tumor occurrence and development.

Pyroptosis Patterns in HNSCC
Tumor cells avoid host immunosurveillance by upregulating tumor mutational burden (TMB), while TMB in turn as an indicator of immunological reaction and tumor behavior. We then merged the transcriptome data of TCGA-HNSCC and GSE65858 and obtained the mRNA expression data of 33 PRGs. A univariate Cox regression analysis determined the Frontiers in Cell and Developmental Biology | www.frontiersin.org May 2022 | Volume 10 | Article 702224 3 prognostic values of these genes in HNSCC (Supplementary Figure S1). We noticed that the expression of most genes in the CASP family affects the prognosis of patients. Thus, as CASP8 demonstrated the highest mutation frequency, we next explored the correlation between the expression of CASP8 and gene mutant and wild types. Compared with wild-type tumors, CASP1 was significantly upregulated, and IL18 was significantly downregulated in CASP8 mutant tumors (Supplementary Figure S2). Meanwhile, we also demonstrated the expressions of four (CASP1, CASP4, CASP8, and CASP9) CASP family genes based on the Human Protein Atlas database (Supplementary Figure S3). The comprehensive landscape of 33 PRG gene interactions, regulator relation, and their prognosis for HNSCC patients were portrayed within the regulator network Figure 2A. The aforementioned results suggested that pyroptosis may act as an important role in TME cell infiltrating characterization among single tumors. Based on the expressions of 33 PRGs, the ConsensusClusterPlus R package was applied to classify patients. Finally, three different clusters were identified by the unsupervised clustering method, involving 172 cases of model A, 354 cases of model B, and 227 cases of model C. The cluster B model indicated a particularly outstanding survival benefit ( Figure 2B). In addition, the heatmap for the association between PRGs and clinicopathological manifestations is also analyzed Figure 2C.

Cluster and Biological Function
GSVA enrichment analysis investigates the biological function of these different clusters. PyroptosisCluster-A was mainly enriched in pathways related to immune activation such as the cytokine-cytokine receptor interaction, chemokine signaling pathway, Toll-like receptor signaling pathways, and cell adhesion molecules cams ( Figure 2D and Supplementary  Table S2). PyroptosisCluster-B mostly included enrichment pathways linked to cancer degradation repairs such as DNA/ RNA degradation, base excision repair, and mismatch repair  Figure 2E and Supplementary Table S3). However, PyroptosisCluster-C was prominently correlated with metabolism pathways such as linoleic acid metabolism, arachidonic acid metabolism, and ether lipid metabolism ( Figure 2F and Supplementary Table S4). We then applied the ssGSEA approach to contrast the element differences of immune cells between the three clusters Figure 3A. The result suggested that there were significant differences in the compositions of immune cell types among the three clusters, which indicated that pyroptosis influences infiltrating cell types of tumors through a different biological process. We next identified the DEGs in each cluster and used the Venn diagram to identify the co-expressed genes. Finally, 2284 co-expressed genes in three clusters were identified using the limma package ( Figure 3B). KEGG analysis showed that these 2284 genes mainly enriched in pathways related to carcinogenic activation pathways such as the TNF signaling pathway, Rap1 signaling pathway, Fc epsilon RI signaling pathway, transcriptional misregulation in cancer, and MAPK signaling pathway ( Figure 3C and Supplementary Table  S5). These genes displayed enrichment of biological process significance associated with carcinogenic pathways, which committed once more that pyroptosis played a key role in tumor occurrence.

Construction of the Pyroptosis Gene Signatures and Functional Annotation
To further explore the possible regulatory mechanism, unsupervised cluster analysis was conducted based on the acquired 2284 pyroptosis phenotypic-related genes to classify patients into distinct genomic subtypes. Three distinct pyroptosis genomic phenotypes were determined and named as pyroptosis gene clusters A/B/C. The analysis showed that three different gene clusters had different characteristic genes ( Figure 3D). Gene cluster C demonstrated the outcome of a poorer prognosis Figure 3E. In the three pyroptosis gene clusters, the main differences in the expression of PRGs were noticed, which was inconsistent with the prospective outcomes of PyroptosisCluster patterns Figure 3F. Considering the individual heterogeneity, a scoring system to quantify the pyroptosis pattern in each patient was constructed based on these phenotype-related genes, which we named as PyroptosisScore. The ggalluvial diagram was used to visualize the relationship between PyroptosisScore, PyroptosisCluster, and GeneCluster Figure 4A. The correlation between the infiltrating cell types and the PyroptosisScore is presented in Figure 4B. Then, the patients were divided into low-or high-PyroptosisScore groups based on the median value. The value of PyroptosisScore in predicting patients' prognosis was also explored, and patients with high PyroptosisScore indicated a bad survival outcome Figure 4C. The Kruskal-Wallis test determined a significant difference in PyroptosisScore among pyroptosis GeneClusters and PyroptosisCluster ( Figure 4D,E). PyroptosisCluster A and B significantly increased when compared to PyroptosisCluster A which indicated that cancer-related pathways indicated high scores Figure 2D-F. Studies reported that the tumor mutational burden (TMB) plays an important role in tumor prognosis. Thus, we next explored the connection between the high/low score group and TMB, and the results have shown that the low score group demonstrated high TMB ( Figure 4F). Interestingly, while the high TMB group suggested a poorer prognosis, when high TMB was combined with a low PyroptosisScore, there was a beneficial outcome indicated which may result in heterogeneity and needs more in-depth research studies ( Figure 4G,H).

Characteristics of the PyroptosisScore
We then used the maftools package to analyze the differences in the distribution of somatic mutations among low and high PyroptosisScores in the TCGA-HNSCC cohort. The results indicated that the low-PyroptosisScore group experienced increased extensive TMB than the high group, with the most significant mutated genes being TP53 and TTN Figure 5A,B. A detailed Spearman correlation analysis was also performed using different algorithms, a lollipop shape to exhibit the result, as shown in Figure 5C. From the results, we can find that most immune cells have a positive correlation with the PyroptosisScore. The relationship between the clinical characteristics and PyroptosisScore was also explored, and we found that a high score may be related to an advanced tumor stage Figure 6A-C. In addition, a high score group means a bad prognosis in stage III-IV patients Figure 6D. We finally used the Frontiers in Cell and Developmental Biology | www.frontiersin.org May 2022 | Volume 10 | Article 702224 6 immunophenoscore (IPS) function based on The Cancer Immunome Atlas database to predict the responsiveness to CTLA-4 and PD-1. The results showed that the high-and low-score groups had a good response to PD-1 and CTL single positive or double negative ( Figure 6E), which can provide potentially useful insights for further analysis. Meanwhile, we also found a statistically significant difference in OS between immunotherapy and risk scores based on the IMvigor210 cohort (Mariathasan Turley Nickles, 2018) (Supplementary Figure S4A). In addition, the AUC value of the PyroptosisScore signature was 0.735, indicating that the prediction accuracy is moderate (Supplementary Figure S4B).

DISCUSSION
Currently, the combined application of multiple immunotherapies and the combined application of immunotherapy, radiotherapy, chemotherapy, and targeted therapy may improve the prognosis of cancer patients by inhibiting tumor progression through the synergistic action of multiple mechanisms. In cancer, pyroptosis is considered to be an autonomous tumor suppressor mechanism, which has a profound impact on the inhibition of tumor progression. Thus, inducing pyroptosis cell death or activating the related signaling pathways is suggested to be a novel way of cancer therapy. This study determined a novel and effective pyroptosisrelated prognostic characterization and the TME cell infiltration dysfunction in the HNSCC patients. The analysis may help us understand the pyroptosis status and TME antitumor immune response and thus provide potential biomarkers for clinical therapeutic intervention.
As a newly discovered method of cell death, pyroptosis has two-sided effects on tumors. On the one hand, the key inflammasomes in the process of pyroptosis can promote tumor cell death and inhibit tumor cell proliferation and metastasis; on the other hand, the accumulation of inflammasomes is conducive to the formation of a microenvironment suitable for tumor cell growth, proliferation, and metastasis. Translocation of GSDMD protein is a necessary procedure for the formation of membrane pores within pyroptosis. A recent study reported that DRD2 can regulate the tumor microenvironment and facilitate M1 polarization of macrophages and activated GSDME-executed pyroptosis in breast cancer (Tan et al., 2021a). Silencing GSDME alleviates cisplatin-or doxorubicin-induced HK-2 cell pyroptosis by raising cell viability and reducing LDH release (Shen et al., 2021). Our study indicated that GSDMC and GSDMD had a general frequency of CNV amplification, and CASP8 demonstrated the highest mutation frequency. CASP8 suggested the molecular switch that controls apoptosis, pyroptosis, and necroptosis, and protects against tissue damage (Fritsch et al., 2019). Caspase-8 can directly cleave caspase-3 to activate apoptosis or GSDMD to trigger pyroptosis, as well as also cleaving RIPK1 and RIPK3 to inhibit necroptosis (Orning and Lien, 2021). Moreover, membrane cholesterol is essential for the activation and recruitment of caspase-8 and its non-apoptotic functions in cancer cells (Kumar et al., 2018). Until now, there is no related study about the role of CASP8 in HNSCC, and our analysis can provide useful information for the in-depth exploration.
Pyroptosis plays an antitumor role by activating the immune response. Due to the defect of the cell membrane, the pyroptotic cells liberate a large amount of cellular substance to cause a strong inflammatory response and a large infiltration of lymphocytes. Lymphocyte infiltration was significantly increased, which further induced caspase-3-independent and -dependent tumor cell pyroptosis, thus forming positive feedback to enhance the anti-tumor function . Colocalization of GSDMD with granzyme B was shown in the proximity of immune synapses, and GSDMD deficiency lessens the cytolytic capacity of CD8 + T cells (Xi et al., 2019). The study also presented that promoting tumor-localized necroptosis and pyroptosis may finally increase responses to immune checkpoint inhibitors (ICI) and tumor-infiltrating lymphocytes (Rosenbaum et al., 2021). Pyroptosis even appears to reprogram the tumor microenvironment to an immunostimulatory state (Tsuchiya, 2021). In our study, biological function analysis revealed that pyroptosis PyroptosisCluster was mainly enriched in pathways related to immunity and cancer, and ssGSEA also suggested that there were significant differences in the compositions of immune cell types among the clusters.
Long non-coding RNA (lncRNA) participates in numerous biological regulatory courses and thus is closely linked to the incidence, progress, and metastasis of cancers. Recent studies have indicated that lncRNAs participated in the malignant phenotypes of tumors not only by genomic or transcriptomic alterations but also through changing the immune microenvironment. A recent study suggested that silencing of lncRNA HOTTIP results in the inhibition of cell proliferation and NLRP1 inflammasome-mediated pyroptosis in ovarian cancer (Tan et al., 2021b). LncRNA ADAMTS9-AS2 served as a tumor suppressor and improved cisplatin sensitivity in gastric cancer cells through activating NLRP3-mediated pyroptotic cell death by sponging miR-223-3p (Ren et al., 2020). In our study, we identified 14 differentially expressed lncRNAs (LINC01963, AL132989.1, AC060780.1, AL451085.2, AC098487.1, PTOV1-AS2, LINC00460, AC012467.2, LINC00958, CTBP1-DT, AC116914.2, AC103702.2, AL022328.2, and AL450992.2), which were chosen as independent prognosis factors of HNSCC patients. Although pyroptosis plays a critical role in tumors, the pyroptosis-related lncRNAs that affect the responsible genes have not been determined. This study integrated some lncRNA biomarkers to assess their effect on patients' prognosis; this can contribute to the determination of novel biomarkers and accurate medical targets in HNSCC. In addition, this study can assist with prognosis prediction, diagnosis, and management strategies of patients with HNSCC. However, further independent studies are required to prove these outcomes and the predictive pyroptosis-related lncRNA. This study also provides some novel insights into cancer immunotherapies that target the genes involved in apoptosis and further reverse the undesirable TME cell infiltration characteristics, which is the transformation of "cold tumor" to "hot tumor," and may contribute to the development of new drug combination strategies or new immunotherapeutic agents. Finally, there were limitations to this study. First, the outcomes were not validated in clinical samples. Second, the results do not provide accurate clinical data due to the relatively small number of patients. Finally, it is unclear whether some genes (such as CASP9 and NLRP7) are pyroptosis-related genes, so we should be cautious when interpreting our results. Although this research explored the probability of establishing a prognostic prediction model, improvement is required. Frontiers in Cell and Developmental Biology | www.frontiersin.org May 2022 | Volume 10 | Article 702224 9