Identification of a pyroptosis‐based model for predicting clinical outcomes from immunotherapy in patients with metastatic melanoma

Abstract Immunotherapy has greatly improved outcomes for patients with advanced melanoma, but good predictive biomarkers remain lacking in clinical practice. Although increasing evidence has revealed a vital role of pyroptosis in the tumor microenvironment (TME), it remains unclear for pyroptosis as a predictive biomarker for immunotherapy in melanoma. RNA sequencing data and annotated clinical information of melanoma patients were obtained from four published immunotherapy datasets. LASSO regression analysis was conducted to develop a pyroptosis‐based model for quantifying a pyroptosis score in each tumor. Based on four clinical cohorts, we evaluated the predictive capability of the model using multiple immunotherapeutic outcomes, including clinical benefits, overall survival (OS), and progression‐free survival (PFS). Furthermore, we depicted the distinctive TME features associated with pyroptosis. Compared with the group with low pyroptosis scores, the group with high pyroptosis scores consistently achieved better durable clinical benefits in four independent cohorts and the meta‐cohort. ROC analysis validated that the pyroptosis‐based model was a reliable biomarker for predicting clinical benefits from immunotherapy in melanoma. Survival analyses showed that the group with high pyroptosis scores harbored more favorable OS and PFS than those with low pyroptosis scores. Molecular analysis revealed that tumors with high pyroptosis scores displayed a typical immune‐inflamed phenotype in TME, including enrichment of immunostimulatory pathways, increased level of tumor‐infiltrating lymphocytes, upregulation of immune effectors, and activation of the antitumor immune response. Our findings suggested that the pyroptosis‐related model associated with multiple immune‐inflamed characteristics might be a reliable tool for predicting clinical benefit and survival outcomes from immunotherapy in melanoma.


| INTRODUCTION
Melanoma is the most malignant cutaneous neoplasm which arises from melanocytes. 1 According to GLOBOCAN 2020, cutaneous melanoma was responsible for approximately 354,635 new cases and 57,043 deaths, accounting for 1.7% of tumor incidence and 0.6% of tumor mortality. 2 Furthermore, it causes the majority (75%) of deaths associated with cutaneous neoplasm due to high invasion and early metastasis. 3,4 In the past decade, the advent of immunotherapy has dramatically shifted the treatment pattern in melanoma. 5 The development of multimodal immunotherapy has brought several promising treatment options for metastatic or unresectable advanced melanoma, such as anti-PD-1/PD-L1 inhibitors, engineered TCR therapy, adoptive T cell therapy, vaccines, and combination strategies. However, due to the nonnegligible heterogeneity among different populations, 7 only a limited number of melanoma patients could achieve durable clinical benefits from immunotherapy in clinical practice. 8 Therefore, developing a predictive biomarker for screening patients who might respond to immunotherapy is of great clinical significance.
Pyroptosis is a novel form of programmed necrotic cell death mediated by the gasdermin (GSDM) protein family, 9,10 mainly characterized by inflammasome activation, caspase activation, and pore formation in the cell membrane. 11 Pyroptosis can cause permeable cell swelling and rupture accompanied by the production of inflammatory cytokines and chemokines, resulting in a potent inflammatory response. 12 Initially, pyroptosis was considered an inflammatory form of the innate immune response to defend the host infection from pathogens. 13 Recently, pyroptosis has been considered a promising opportunity to reduce immunosuppression and enhance the antitumor immune response. [14][15][16] Single-cell RNA sequencing analysis suggested that inducing pyroptosis could promote the recruitment of tumor-infiltrating immunostimulatory cells and inhibit the activation of immunosuppressive cells in breast cancer. 17 Furthermore, a molecular classification study found that subtypes with high pyroptosis-related gene expression belonged to the immunological "hot" phenotype with activated immune response in bladder cancer. 18 Bioinformatic analyses also identified a series of pyroptosis-based signatures associated with immune status and prognosis in melanoma. 19,20 However, the relevance of pyroptosis with the efficacy of tumor immunotherapy remains unclear in melanoma.
This study developed a pyroptosis-based model to predict clinical benefits from immunotherapy in melanoma using four published immunotherapy cohorts. [21][22][23][24] In addition, we explored the prognostic value of the pyroptosisbased model in melanoma patients with or without immunotherapy. Furthermore, we depicted the unique tumor microenvironment (TME) characteristics associated with the pyroptosis-based model to explain the predictive value for immunotherapy.

| Datasets collection
We performed a literature search in the Medline, Web of Science, and EMBASE databases to retrieve the public available immunotherapy cohorts. The keywords included melanoma, immunotherapy, and RNA sequencing. Studies were included when met the following criteria: (1) tumor samples were pathologically confirmed as melanoma; (2) patients received at least one kind of immunotherapy; (3) original data on response, clinical benefit, or survival outcomes were provided; (4) studies were published in English; (5) RNA-seq data were available. The exclusion criteria were as follows: (1) studies without detailed patient information; (2) duplicated studies; (3) review, letter, editorial, comment, abstract, and case report. As a result, this study enrolled four eligible cohorts with RNA-seq and related clinical data, including the Gide, 21 Lauss, 23 Liu,22 and Nathanson 24 cohorts. The Gide cohort contained 73 melanoma patients receiving anti-PD-1 monotherapy or the combination between anti-PD-1 and anti-CTLA-4 therapy. The Lauss cohort consisted of 25 melanoma patients receiving adoptive T cell therapy. The Liu cohort included 119 melanoma patients receiving anti-PD-1 monotherapy. The Nathanson cohort contained 21 melanoma patients receiving anti-CTLA-4 monotherapy. In addition, the clinicopathological features, including age, gender, M stage, primary tumor site, therapy type, therapeutic response, clinical benefit, and survival outcomes, were recorded in Table 1. This study aimed to develop an effective biomarker for immunotherapy in melanoma. The Gide cohort was selected as the training cohort, and the others were chosen as validation cohorts. Besides, the RNA sequencing with annotated clinical data of the TCGA-SKCM cohort was downloaded from the cBioPortal database. After removing the cases with incomplete RNA sequencing and clinical data, a total of 402 melanoma patients were enrolled in the TCGA-SKCM cohort.

| Construction of the pyroptosisbased model for predicting response to immunotherapy
First, the geneset containing 33 pyroptosis-related genes was collected from previous studies, [25][26][27] as shown in Table S1. The gene expression matrix with annotated clinical information was extracted from the Gide cohort. Second, LASSO regression analysis was conducted by the R package "glmnet" to identify the pyroptosis-related genes closely associated with clinical benefits from immunotherapy. The gene number was determined according to the ideal penalty (Lambda). Third, the gene coefficient (coef) was obtained from the LASSO regression model. Subsequently, a pyroptosis-based model, named pyroptosis score, was constructed according to the combination between the coefficient (coef) and gene expression level as follows: pyroptosis score = ∑ (coef i × gene i ). Following that, the pyroptosis score of each tumor sample in the Gide, Lauss, Liu, and Nathanson    cohorts was determined using this formula. The distribution of the pyroptosis score in each cohort was visualized using the R package "waterfalls". In addition, four machine learning methods, including random forests (RF), support vector machines (SVM), artificial neural networks (ANN), and K-nearest neighbor (KNN), were performed by the R package "Caret" to build models for predicting clinical benefits from immunotherapy in melanoma.

| ROC analysis and survival analysis
The predictive efficiency for clinical benefits from immunotherapy was tested using receiver operating characteristic (ROC) analyses in the Gide, Lauss, Liu, and Nathanson cohorts. The area under the ROC curve (AUC) was determined using the R package "pROC". In addition, Kaplan-Meier analysis was conducted to evaluate the prognostic value of the pyroptosis model in melanoma patients with immunotherapy. The median value of the pyroptosis score was adopted as the grouping threshold.

| Logistic regression analysis, Cox regression analysis, and meta-analysis
The Odds ratios (OR) and 95% confidence interval (CI) were determined using univariate logistic regression analysis to estimate the difference in DCB between the groups with different pyroptosis scores. The hazard ratio (HR) and 95% CI were determined using univariate Cox regression analysis to assess the association of the pyroptosis score with OS and PFS. Based on these results, a meta-analysis was conducted using the R package "meta" to evaluate the combined effect and the interstudy heterogeneity. The interstudy heterogeneity was defined according to the value of I 2 as follows: low (I 2 < 25%), moderate (I 2 = 25-75%), or high (I 2 > 75%). 28 A fixed-effect model was applied in the meta-analysis when I 2 was less than 50%. Otherwise, a random-effect model was applied in the meta-analysis.

| Function annotation and pathway enrichment analysis
First, considering the limited number of samples in a single cohort and the low heterogeneity among the clinical cohorts, we merged the Gide, Lauss, Liu, and Nathanson cohorts as an overall cohort for molecular analyses. The batch effect from different cohorts was corrected using the R package "sva". Then, the differentially expressed genes (DEGs) between the groups with different pyroptosis scores were screened using the R package "limma". |Fold change (FC)| > 1.5 with the adjusted p < 0.05 was adopted as the thresholds. Based on these DEGs, Gene Ontology (GO), Kyoto Gene and Genomic Encyclopedia (KEGG), and gene set enrichment analysis (GSEA) analyses were performed using the R packages "GO.db", "org.Hs.eg. db", "topGO", "DOSE", "GSEABase", "clusterProfiler", and "stringr". For GSEA analysis, the REACTOME genesets were downloaded from the Molecular Signatures Database (MsigDB). The number of random sample permutations was adopted as 1000.

| Evaluation of tumor-infiltrating immune cells, TME scores, immune effectors, and transcriptional biomarkers for immunotherapy
CIBERSORT was a deconvolution algorithm widely applied to estimate tumor-infiltrating immune cells, which output the relative infiltration level of immune cells.
In this study, we also evaluated the infiltration level of 22 immune cell types in each tumor sample using the CIBERSORT algorithm. The R package "ESTIMATE" was applied to determine the TME scores, including ESTIMATE score, stromal score, and immune score, which reflected tumor purity, stromal content, and immune cell infiltration level. In addition, immune-related genes were obtained from published literature, [29][30][31][32][33][34][35] including activated DC markers, antigen-processing related genes, CD8 + T effectors, genes downstream of IFN-γ, NK markers, cytolysis molecules, and immune checkpoints. Furthermore, we explored the correlation between the pyroptosis score and several established transcriptional biomarkers for immunotherapy, including cytolytic activity (CYT) 34 and the T cell-inflamed gene expression profile (GEP). 32 As previously reported, CYT was determined by the mean expression of GZMA and PRF. GEP was determined by the mean expression of a T cell-inflamed gene expression profile, including CD276, CD8A, HLA-E, CXCL9, CXCR6, CCL5, PDCD1LG2, HLA-DQA1, HLA-DRB1, PSMB10, CD27, IDO1, LAG3, NKG7, CMKLR1, CD274, STAT1, and TIGIT.

| Statistical analysis
All data analyses were performed in SPSS 21.0 and R software 3.6.1. The correlation between the pyroptosis score and the pyroptosis level, the immune effectors, and the transcriptional biomarkers for immunotherapy (CYT and GEP) were determined using Spearman's correlation analysis. Mann-Whitney test was applied to compare the difference in immune cell proportions and TME scores between the low and high pyroptosis groups. The comparison of clinical benefit rate was dealt with chi-squared (χ 2 ) test or Fisher's exact test. The log-rank test was used to compare Kaplan-Meier curves. A P value <0.05 was defined as statistically significant.

| Baseline characteristics of melanoma patients with immunotherapy
We collected the baseline characteristics of 238 melanoma patients with immunotherapy from four previous studies, including the Gide, 21 Lauss, 23

| Construction of a pyroptosisbased model for predicting DCB from immunotherapy in melanoma
The gene expression matrix was matched with the geneset of 33 pyroptosis-related genes and extracted from the Gide cohort. Then, LASSO regression analysis was performed to screen candidate genes significantly associated with DCB from immunotherapy. As a result, four ideal genes were enrolled in the predictive model ( Figure 1A,B). Based on the gene expression level with variable coefficient, a pyroptosis-related model was developed according to the following formula: pyroptosis score = (0.139 × CASP5) According to the predictive model, each tumor sample in four immunotherapy cohorts was assigned with a pyroptosis score. Following that, we investigate the association between pyroptosis score and 33 pyroptosis-related genes in the immunotherapy cohorts. It was found that tumor samples with higher pyroptosis scores harbored more abundant expression levels of pyroptosis-related genes in the Gide, Lauss, Liu, and Nathanson cohorts ( Figure 1C). Subsequently, we quantified the pyroptosis level in each tumor sample using the mean expression level of 33 pyroptosis-related genes. Correlation analysis showed that the pyroptosis scores were all significantly correlated with the pyroptosis level in the Gide (R = 0.77, p < 0.001), Lauss (R = 0.48, p < 0.001), Liu (R = 0.53, p < 0.001), and Nathanson (R = 0.89 p < 0.001) cohorts ( Figure 1D).

| Prediction of durable clinical benefit from immunotherapy by the pyroptosis-based model
The pyroptosis-based model was applied to predict the clinical benefits of immunotherapy in melanoma patients. First, we investigated the association of the pyroptosis score with DCB in the clinical cohorts. As shown in Figure 2A, a decreasing trend in the DCB proportion was observed along with the decrease in pyroptosis score in each cohort. We then compared the proportion of DCB between the groups with different pyroptosis scores. Expectedly, the group with high scores harbored a higher rate of DCB than the group with low scores in the Gide (78.4% vs 47.2%, p = 0.006, Figure 2B), Lauss (69.2% vs 8.3%, p = 0.004, Figure 2B), and Liu (68.3% vs 33.9%, p < 0.001, Figure 2B) cohorts. A similar trend was also observed in the Nathanson cohort, although it was not significant (63.6% vs 20.0%, p = 0.081, Figure 2B). Subsequently, ROC analysis suggested that the pyroptosis score was an effective and robust biomarker for predicting DCB from immunotherapy across various cohorts  Figure S1). Furthermore, meta-analysis validated that patients with high pyroptosis scores had a higher DCB rate than those with low pyroptosis scores (pooled OR = 4.80, 95% CI = [2.72, 8.48], p < 0.001, Figure 2D). The heterogeneity was mild among these cohorts (I 2 = 0%, tau 2 = 0, p = 0.53, Figure 2D).
To compare the predictive value of different machine learning models, we further applied four machine learning methods, including random forests (RF), support vector machines (SVM), artificial neural networks (ANN), and K-nearest neighbor (KNN), to build models for   Figure S2D). However, compared with these models, the LASSO regression model used in this study harbored a more robust predictive performance (all AUC >0.7, Figure 2C).

| Association of the pyroptosisbased model with OS and PFS in the immunotherapy cohorts
We performed Kaplan-Meier analysis to evaluate the prognostic value of pyroptosis-based model on OS and PFS in the immunotherapy cohorts. First, we observed that pyroptosis scores had a significantly protective value on OS in the immunotherapy cohorts  Figure 4D). Additionally, it remains unknown that the prognostic value of the pyroptosis-based model for melanoma patients who did not receive immunotherapy. Therefore, we further explored the association of the pyroptosis score with OS and PFS in TCGA-SKCM cohort. It was found that the differences of OS and PFS between the groups with low and high pyroptosis scores were both not significant in melanoma patients without any immunotherapy ( Figure S3).

| Predictive and prognostic value of the pyroptosis-based model in patients with other cancer types receiving immunotherapy
To evaluate the applicability of the pyroptosis-based model in other carcinomas, we collected two immunotherapy datasets from published cancer-related literature, including metastatic gastric cancer (GC) 36 and advanced clear cell renal cell carcinoma (ccRCC). 37 The Kim cohort consisted of 45 GC patients who received pembrolizumab therapy. 36 The Braun cohort contained 281 ccRCC patients treated with nivolumab. 37 Figure S4D). Taken together, we considered that the predictive value of the pyroptosis-based

| Function annotation and pathway enrichment of the pyroptosis-based model
The pyroptosis-based model harbored a reliable and robust predictive value for melanoma patients treated with immunotherapy. Therefore, it is valuable to investigate the molecular basis underlying the model's predictive value. First, we screened 673 DEGs between tumor samples with different pyroptosis scores using the R package "limma". Subsequently, we conducted GO function annotation and KEGG pathway enrichment analyses to investigate the biological functions of DEGs. Interestingly, DEGs were annotated with immune-related biological functions ( Figure 5A). At the biological process level, DEGs were mainly involved in the T cell activation, regulation of leukocyte activation, and regulation of lymphocyte activation. At the level of the cellular component, DEGs were mostly annotated with the side of the membrane, external side of plasma membrane, and plasma membrane receptor complex. At the molecular function level, DEGs presented a series of immune-related functions, such as MHC protein complex binding, cytokine receptor activity, and MHC protein binding. KEGG pathway enrichment analysis showed that immune-related pathways account for 70% of the top10 signaling pathways, including cytokine-cytokine receptor interaction, hematopoietic cell lineage, Th17 cell differentiation, Th1 and Th2 cell differentiation, intestinal immune network for IgA production, allograft rejection, and graft-versus-host disease ( Figure 5B). Furthermore, GSEA analysis showed that pyroptosis score was positively correlated with several immune-related pathways, including antigen processing cross presentation, cytokine signaling in the immune system, downstream TCR signaling, immunoregulatory interactions between a lymphoid and a non-lymphoid cell, interleukin-1 family signaling, interleukin-4 and interleukin-13 signaling, MHC class II antigen presentation, programmed cell death, TCR signaling, and TNFR2 non-canonical NF-kβ pathway (All adjusted p < 0.05, Figure 5C).

| Association between the pyroptosisbased model and tumor immune microenvironment in melanoma
Considering the significant association between DEGs and many immune-related signaling pathways, we further explore the pyroptosis-based model's role in melanoma's tumor immune microenvironment. First, we performed CIBERSORT algorithm to estimate the relative proportion of tumor-infiltrating immune cells. Then we compared the abundance of 22 immune cell types between the groups with different pyroptosis scores. We observed that tumors with high pyroptosis scores harbored significantly elevated infiltration level of CD8 + T cells (p < 0.001), activated memory CD4 + T cells (p < 0.001), polarized Macrophages M1 (p < 0.001) and M2 (p < 0.001), and plasma cells (p < 0.001) than those with low pyroptosis scores ( Figure 6A). On the contrary, tumors with low pyroptosis scores contained more resting immune cells, including naïve CD4 + T cells (p < 0.05) and macrophage M0 (p < 0.001). Second, the ESTIMATE algorithm was applied to evaluate TME scores in tumor samples. We observed that tumors with high pyroptosis scores had significantly higher ESTIMATE (p < 0.001), stromal (p < 0.001), and immune scores (p < 0.001) than those with low pyroptosis scores ( Figure 6B), suggesting that tumors with high pyroptosis scores had more abundant TME composition, more stromal content, and elevated infiltration level of the immune cells. In addition, we found that tumors with high pyroptosis scores harbored more immune effectors, including the activated DC cell markers, antigen-processing related genes, CD8 + T cell effectors, genes downstream of IFN-γ, NK cell markers, and immune checkpoints (All p < 0.05, Figure 6C and Table S2). Besides, previous studies identified CYT and GEP as the predictive biomarkers at the transcriptional level to predict response to immunotherapy in various malignancies, which reflected T cell cytolytic activity and inflammation level in the tumor microenvironment. 32,34 This study further investigated the correlation between the pyroptosisbased model and these biomarkers. There was a significantly positive correlation of the pyroptosis score with  CYT (R = 0.64, p < 0.001, Figure 6D) and GEP (R = 0.66, p < 0.001, Figure 6D). In summary, these findings support that tumors with high pyroptosis scores tended to develop a "hot" immune phenotype, whereas those with low pyroptosis scores showed a "cold" immune phenotype. Therefore, there was an increased likelihood of clinical benefits of immunotherapy for melanoma patients with high pyroptosis scores (Figure 7).

| DISCUSSION
Pyroptosis has a dual role in tumorigenesis and antitumor immunity across diverse cancers. On the one hand, the release of pyroptosis-producing cytokines could drive the accumulation of immunosuppressive cells (such as MDSC, T reg , and M2 cells) and regulate the crosstalk between tumor cells and TME components, eventually leading to tumorigenesis and progression. 14 On the other hand, pyroptosis could enhance the antitumor effect via recruiting adaptive immune cells, activating cytolytic lymphocytes, and promoting the phagocytosis of macrophages. 16 Despite these advances, most studies focus on a single infiltrating immune cell or pyroptosis-related gene. Therefore, an integrated analysis for investigating the role of pyroptosis in TME in melanoma is urgently needed. Furthermore, it remains unclear whether pyroptosis-related genes could be applied in clinical practice to predict clinical outcomes from immunotherapy. Therefore, constructing a pyroptosis-based model for immunotherapy will deepen our understanding of the clinical application of pyroptosis in tumor immunotherapy and assist oncologists in identifying the advantage population with clinical benefits from immunotherapy.
In this study, we identified a predictive model using four pyroptosis-related genes for melanoma patients with immunotherapy. The predictive efficiency of the model for clinical benefit and survival outcomes from immunotherapy was explored and validated in four independent cohorts. Although the immunostimulatory role of pyroptosis has been revealed in melanoma, 38,39 this study highlights the clinical application of the pyroptosis-based model in tumor immunotherapy. We observed that the pyroptosis-based model achieved    CD8A  CD8B  CD3D  CD3E  CD2  CD27  IDO1  CXCL10  CXCL9  HLA-DRA  IFNG  CMKLR1  PSMB10  STAT1  HLA-DQA1  HLA-DRB1  NKG7  HLA-E  GZMA  GZMB  GZMK  PRF1  CD274  PDCD1  PDCD1LG2 LAG3 CTLA4 an effective and robust prediction for clinical benefit, OS, and PFS in the Gide, Lauss, Liu, Nathanson, and meta cohorts. To our knowledge, it is the first time to investigate the predictive value of the pyroptosis-based model for immunotherapy efficacy in melanoma based on multiple immunotherapy cohorts. Furthermore, it is worth noting that the pyroptosis-based model only depends on the gene panel's expression level, which is low-cost and convenient in clinical application. Besides, we also evaluated the association between the pyroptosis-based model and survival outcomes for melanoma patients without immunotherapy using the TCGA-SKCM cohort. Interestingly, the prognostic value of the model was not significant ( Figure S3), suggesting that the pyroptosis-based model was a predictive biomarker for immunotherapy instead of a prognostic biomarker in melanoma. Therefore, on the other hand, it is recommended for patients with high pyroptosis scores rather than those with low scores to receive immunotherapy in clinical practice; on the other hand, increasing the pyroptosis level might be a promising strategy to improve sensitivity to immunotherapy for melanoma patients. Converting tumor immunophenotype from "cold" to "hot" by triggering pyroptosis may be a potentially effective strategy for tumors with immunedesert characteristics. Recent studies have developed novel concepts and methods to integrate pyroptosis-based therapy with immunotherapy to cure tumors. Zhao et al. designed a pyroptosis-associated biomaterial nanoparticle to improve the immunotherapy efficacy for tumors with low immunogenicity. 15 Zhang et al. developed a covalent organic framework for inducing pyroptosis and triggering durable antitumor immunity to boost tumor immunotherapy. 40 Yang et al. summarized a series of compounds that could activate pyroptosis in cancers, 41 combined with immunotherapy to exert an antitumor effect. Therefore, combining pyroptosis-inducing therapy with immunotherapy might be a promising strategy for patients with low pyroptosis scores. The development of the pyroptosis-based model has become a hot spot in melanoma research. Cao et al. constructed a prognostic model of genes associated with pyroptosis in melanoma. 19 Similarly, Zhang et al. also developed a pyroptosis-related gene signature to predict the prognosis of melanoma patients. 38 Wang et al. developed a pyroptosis model for predicting the characteristics of the immune microenvironment, the sensitivity to immunotherapy, and the prognosis of melanoma patients. 42 Meng et al. characterized the pyroptosis-related immune infiltration landscape in melanoma. 39 Different from Zhang et al.'s and Wang et al.'s models, our pyroptosis-based model focused on the predictive efficacy of clinical benefit from immunotherapy in melanoma rather than the prognosis of melanoma patients. More importantly, previous models were constructed based on the clinical cohorts where most patients did not receive immunotherapy, such as TCGA or GEO datasets. On the contrary, our study developed and validated the pyroptosis-based model based on the immunotherapy cohorts. It is the first time to systematically evaluate the predictive accuracy of the pyroptosis-based model using four independent immunotherapy cohorts. Summarily, we considered that the predictive value of our model was more reliable and robust than previous models.

Pyroptosis score DCB
Previous studies have revealed a complex molecular mechanism underlying the antitumor effect of pyroptosis. First, the production of proinflammatory cytokines and immunostimulatory molecules, such as IL-1β, 43 IL-18, 44 HMGB1, 45 and ATP, 46 could activate cytotoxic CD8 + T cells and Th1 CD4 + T cells but inhibit the differentiation of regulatory T cells, eventually leading to the antitumor immune response. Second, GSDM protein could activate the tumor-infiltrating lymphocytes in TME and upregulate the expression of immune effectors, including IFN-γ, GZMB, TNF-α, and PRF. 47,48 Besides, pyroptotic macrophages could trigger the cytolytic activity of NK cells against cancer cells. 49 These findings indicated that tumors with high pyroptosis levels harbored a "hot" immunophenotype, with great potential for immunotherapy response. We also explored the association between the pyroptosis-based model and various aspects of immunological characteristics in TME to explore the molecular mechanism underlying the predictive value for immunotherapy. It was found that a number of immune-related signaling pathways were significantly enriched in tumors with high pyroptosis scores, such as antigen processing cross presentation, TCR signaling, and cytokine signaling in immune system. Furthermore, tumors with high pyroptosis scores harbored more abundant infiltration levels

Likelihood for clinical benefit from immunotherapy
Low High of cytotoxic CD8 + T cells and activated CD4 + T cells. In contrast, tumors with low pyroptosis scores contained more resting immune cells, such as naïve CD4 + T cells and macrophage M0. At the level of molecules, we observed positive correlations between pyroptosis score and various immune effectors, including activated DC cell markers, antigen-processing related genes, CD8 + T cell effectors, genes downstream of IFN-γ, NK cell markers, and immune checkpoints. It could be explained that tumors with high pyroptosis scores, featured by the immune-inflamed phenotype, harbored a higher likelihood of clinical benefit from immunotherapy. Four pyroptosis-related genes were enrolled in the predictive model, including CASP5, NLRP6, NLRP7, and PYCARD. The biological functions of these genes have been reported in previous studies. CASP5, encoding a member of the caspase family, plays a key mediator in programmed cell death. 50 CASP5 initiates the cleavage of GSDMD and the production of the N-terminal GSDMD moiety that binds to membranes and forms pores, eventually leading to cell expansion and rupture. 51 Hu et al. found that patients with high expression of CASP5 had a worse survival outcome than those with low expression in kidney carcinoma. 52 Ulybina et al. found that coding polymorphisms in CASP5 might be a key regulator in predisposition to lung cancer. 53 NLRP6 is known as a sensor element in the NLRP6 inflammasome. 54 It promotes the polymerization of the inflammasome to participate in innate immunity and inflammation. 55 Patients with low expression of NLRP6 showed a worse OS in multiple cancer types, such as liver cancer, 56 breast cancer, 57 and gastric cancer. 58 NLRP7, another member of the NLR family, is involved in activating proinflammatory caspases. 59 Li et al. demonstrated that NLRP7 promotes CRC progression by polarizing macrophages M2. 60 Reynaud et al. found that NLRP7 promoted tumorigenesis and progression of choriocarcinoma by constructing an immunosuppressive TME. 61 Interestingly, in this study, NLRP7 took the largest coefficient in the pyroptosis-based model, which was positively correlated with an immune-inflamed TME in melanoma. PYCARD, encoding an adaptor protein, functions as a key regulator in inflammation and programmed cell death. 62 A prior study showed that elevated expression of PYCARD contributed to the tumor cell growth and poor prognosis of pancreatic cancer. 63 Furthermore, Strickler et al. identified the diagnostic value of PYCARD in melanoma using immunohistochemical analysis. 64 Nevertheless, the biological functions of these pyroptosis-based genes involved in tumor immunity and immunotherapy have rarely been investigated in melanoma, deserving further investigation in vitro/ in vivo experiments.
This study had several limitations. Firstly, the best cutoff value of the pyroptosis score remains a mystery for clinical application of the pyroptosis-based model. Secondly, there are several potential confounding factors, including the heterogeneity of different studies, the selection bias of cohorts, and the batch effect of RNA-seq data from different platforms. Thirdly, although our research has preliminarily established the association between the model and immunotherapy, it remains unclear which immunotherapy strategies are best suited for this model. In addition, it is valuable to explore the molecular mechanism underly the association between these pyroptosis-related genes and tumor immunotherapy in the future. Besides, although the prediction capability of the pyroptosis-based model has been explored in four independent immunotherapy cohorts, it would be ideal if our findings could be validated in prospective, large-scale, and multi-center studies. Furthermore, limited by the incomplete clinical information, the pyroptosis-based model established in this study was only based on gene expression. It would be better to develop a predicting model based on gene and clinical data, which contained more comprehensive features associated with immunotherapy. In the future, a combined model might better assist oncologists in identifying the advantage population with clinical benefits from immunotherapy.
In conclusion, our study developed a reliable and robust model based on four pyroptosis-related genes to predict clinical outcomes from immunotherapy in melanoma. Comprehensive analysis of four independent cohorts consistently supported that patients with high pyroptosis scores had favorable clinical outcomes from immunotherapy than those with low pyroptosis scores. Furthermore, we comprehensively analyzed the TME characteristics associated with pyroptosis score to explain the predictive and prognostic value of the model for immunotherapy. These findings may assist oncologists in predicting immunotherapeutic efficacy and stratifying melanoma patients, which deserve further investigations in the future.

AUTHOR CONTRIBUTIONS
Haiyong Wang designed the study. Guanghao Wu, Biying Chen, and Junjie Jiang analyzed the data and wrote the manuscript. Guanghao Wu, Junjie Jiang, Biying Chen, and Yiran Chen drew figures. Guanghao Wu, Biying Chen, and Yanyan performed a literature search. Junjie Jiang, Yanyan Chen, and Haiyong Wang reviewed the manuscript. All authors have read and approved the final version of the manuscript.