RNA-seq data and corresponding clinical data were obtained from The Cancer Genome Atlas (TCGA; https://tcga-data.ncinih.gov/tcga/). The expression values of all genes were downloaded in fragments per kilobase of transcript per million mapped reads format. The RNA-seq data of normal skin were downloaded from the Genotype-Tissue Expression (GTEx) database (https://gtexportal.org/home/). The transcriptome samples included 471 tumor samples and 1 normal sample. Among them, 458 clinical samples were included for prognosis analysis because 13 patients had missing survival time and survival state. The clinical characteristics of the 458 patients are displayed in Supplementary Table 2. The merged RNA expression profile in TCGA-melanoma and GTEx-skin, including 471 CM samples and 234 normal samples, was normalized. The batch effects between TCGA and GTEX data were removed using the R software package “limma.” TCGA and Gene Expression Omnibus (GEO) data were normalized using R software packages “limma” and “sva” to remove biases from different platforms. We identified 40 pyroptosis-related genes and 675 inflammatory response-related genes (Supplementary Table 2) using the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/msigdb/).

To identify differentially expressed GRIPs (DEGs), we analyzed the gene transcription data of the TCGA and Genotype-Tissue Expression (GTEx) datasets using the R package limma (42) with a false discovery rate < 0.05 and log2FC > 1. Heatmaps were conducted by using the R package “pheatmap.” Analyses of the Gene Ontology (GO) database (which includes biological process, cellular component, and molecular function information) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database were used to evaluate the pathways associated with DEGs.

To obtain a prognostic GRIP signature, we first performed univariable Cox regression to determine the association between DEGs and overall survival (OS). DEGs found to be correlated with OS (p < 0.05) in the univariable Cox regression were then selected for Lasso regression to further identify the key GRIPs according to dynamic coefficient profiling and lambda. Finally, we performed multivariable Cox regression analysis, including these key Lasso genes, to construct an optimal risk score. We used a heatmap to evaluate the relevance of the risk score and clinicopathological features.

The risk score was calculated as follows: risk score = ∑i=1nβi * (expression of GRIPs), where n is the number of key GRIPs and β is the regression coefficient. The risk score for each CM patient was calculated, and patients were classified as either low-risk or high-risk using the median risk score as the cutoff. All TCGA samples were randomly divided into a training dataset and a testing dataset at a 1:1 ratio. The GEO dataset GSE65904 was used as an external validation cohort. The validity of the risk score was verified in the training, testing, and entire cohort, as well as the external cohort using receiver operating characteristic (ROC) curves, risk plots, and survival analysis.

To determine whether the risk score was an independent predictor of OS, we performed univariable and multivariable Cox regression analyses, which included the risk score, patient age, patient sex, melanoma stage, and TNM stage (T: tumor size; N: lymph node involvement; M: metastasis), in the entire cohort. Hazard ratios (HRs) and 95% confidence intervals (CIs) for each predictor were calculated.

The R packages “rgplot,” “survival,” and “rms” were used to construct nomograms (43) predictive of the 1-, 3-, and 5-year OS rates of patients with CM. The calibration curves of the nomograms in predicting 1-, 3-, and 5-year OS rates were constructed using the R package rms. The areas under the ROC curves (AUCs) for 3-year OS were predicted with independent predictors, and the nomogram was calculated to evaluate the predictive values of the GRIP signature and the nomogram.

We performed GSEA to identify the biological processes and pathways that were significantly alerted between the high-risk and low-risk subgroups based on the KEGG gene set “c2.cp.kegg.v7.4.symbols.gmt,” the Hallmark gene set “h.all.v7.4.symbols.gmt,” and other inflammation- and pyroptosis-related gene set. The TCGA gene expression was used as a phenotype label (44). The normalized enrichment score (NES) was calculated for each gene set. A nominal p-value of < 0.05 and a NES of >1.5 were used to identify significantly enriched pathways.

We used the CIBERSORT (45, 46), ESTIMATE (47, 48), MCPcounter (49, 50), single-sample GSEA (51), and TIMER algorithms (52) to compare the cellular components and immune responses between the high- and low-risk groups. Heatmaps were used to detect the differences in the immune response under different algorithms. In addition, ssGSEA was used to assess immune cell infiltration and function. Immune-related survival analysis was performed based on ssGSEA.

We systematically identified 36 immune checkpoint-related genes that might be correlated with immune checkpoint inhibitors' (ICIs) response from previous studies. To compare the effect of ICI between the high- and low-risk groups, we assessed the programmed death 1 (PD-1) and cytotoxic T-lymphocyte antigen 4 (CTLA-4) using data from The Cancer Immunome Atlas (https://tcia.at/). TIDE algorithm was used to calculate the TIDE and Dysfunction scores to assess responsiveness to immunotherapy and level of immune dysfunction (http://tide.dfci.harvard.edu/) (53).

We used the “pRRophetic” package in R to evaluate the half-maximal inhibitory concentration (IC50) of common chemotherapy and molecular drugs to estimate the predictive role of the GRIP signature for CM treatment, such as rapamycin, docetaxel, and imatinib (54, 55).

The protein expression levels of the GRIPs in normal and tumor tissues were verified using semiquantitative analysis of the immunohistochemical data extracted from The Human Protein Atlas (https://www.proteinatlas.org/), a database that includes immunohistochemistry-based expression data for ~20 of the most common types of cancers (56).

Human melanoma cell line: A375 was purchased from American Type Culture Collection (ATCC, America); SK-MEL-28 cell line was purchased from National Infrastructure of Cell Line Resource. Human immortalized keratinocytes cell line (Hacat) and human skin melanocyte cell line (PIG1) were purchased from Shanghai Guandao Biological Engineering Company with STR certifications. All cells were cultured in Dulbecco's modified Eagle's medium (Gibco) + 10% fetal bovine serum 10% fetal bovine serum under a humidified atmosphere of 37°C and 5% CO₂.

The real-time quantitative PCR (qRT-PCR) was performed to determine the relative expression levels of the genes. Total RNA was extracted using TRIzol reagent from two melanoma cells (A375, SK-MEL-28) and two control cells (Hacat, PIG1). The concentration of RNA was determined by ultraviolet spectrophotometry. Then, RNAs were reverse transcribed into cDNAs using the PrimerScript RT Master Mix. Finally, SYBR (Takara, Japan) was used to evaluate the mRNA expression levels. The mRNA expression levels of TLR1, IFNGR2, CCL25, IL15, CCL8, NLRP6, EMP3, and RTP4 were normalized by glyceraldehyde-3-phosphate dehydrogenase (GAPDH) that was used as the internal reference. The primer sequences of GAPDH and the eight genes are listed in Supplementary Table 3.

All statistical analyses were performed using R software (version 4.0.2). Cox regression and Lasso regression analyses were used to assess the predictive value of the GRIP signature. The Kaplan-Meier method was used to assess the survival of CM patients. The AUC was calculated from the ROC curve to assess prediction accuracy. Normally and nonnormally distributed variables were analyzed using the unpaired Student's t-test and the Wilcoxon test. All statistical tests were two-sided with p ≤ 0.05 being statistically significant.

The co-expression network first identified 206 inflammatory response-related genes that were highly correlated with pyroptosis-related genes (Figure 1A). Using PCA pyroptosis scores, we found that the pyroptosis level of CM tissue was significantly higher than that of normal tissue (Figure 1B). Patients with CM were divided into high and low pyroptosis score groups (Figure 1C). The OS rate of patients with high pyroptosis scores was significantly higher than patients with low pyroptosis scores (Figure 1D). Similarly, GSVA showed that CM tissue had a higher pyroptosis level than normal tissue (Figure 1E) and patients with CM and high GSVA pyroptosis scores had higher OS rates than those with low GSVA pyroptosis scores (Figure 1F). WGCNA of the 675 inflammatory response-related genes identified 473 genes in the ME blue and ME turquoise modules that were highly related to PCA and GSVA pyroptosis scores (Figures 1G–J). The 206 genes obtained in the co-expression network and the 473 module genes obtained by WGCNA were intersected to identify 185 GRIPs (Figure 1K).

Of the 185 GRIPs, 134 were differentially expressed between tumor and normal tissue samples. In the tumor samples, 106 of these genes were upregulated compared with those in the normal tissue samples and 28 were downregulated (Figures 2A,B). GO analysis of the 134 DEGs showed that positive regulation of cytokine production, T-cell activation, leukocyte cell-cell adhesion, and regulation of T-cell activation was enriched in tumor samples (Figure 2C). KEGG analysis showed that cytokine-cytokine receptor interaction, chemokine signaling, Th17 cell differentiation, Toll-like receptor (TLR) signaling, PD-L1 expression, and PD-1 checkpoint pathways were significantly enriched in tumor samples (Figure 2D).

Univariable Cox regression revealed that 83 of the 134 DEGs were correlated with OS (Figure 3A). Lasso regression identified nine prognostic GRIPs based on the optimal lambda value (Figures 3B,C). Then, the multivariate Cox regression finally identified eight GRIPs (EMP3, TLR1, IFNGR2, IL15, × CCL8, × NLRP6, CCL25, RTP4) for construction of the prognostic signature (Figure 3D). The risk score was calculated as Risk score = 0.003 × EMP3 exp-0.065 × TLR1 exp-0.012 × IFNGR2 exp-0.288 × IL15 exp-0.057 × CCL8 exp-0.633 × NLRP6 exp-0.329 × CCL25 exp-0.024 × RTP4 exp. The Sankey diagram showed the prognostic co-expression relationship between the pyroptosis genes and the GRIP signature (Figure 3E). A heatmap of the association between the risk signature and clinicopathological characteristics demonstrated that low immune scores were associated with higher risk scores (Figure 3F).

The GRIP signature had an AUC of 0.776, 0.776, and 0.713 for predicting 1-, 2-, 3-year OS in the training cohort (Figure 4A). The risk plots showed that the OS rates of patients gradually decreased as their risk scores gradually increased (Figures 4B,C). The patients in the training cohort were divided into high- and low-risk groups using the median risk score as a cutoff, and PCA showed the two groups to be well distinguished (Figure 4D). In the training cohort, the OS rate of the high-risk group was significantly lower than that of the low-risk group (Figure 4E). The prognostic GRIP signature was validated in the testing cohort and entire cohort (Figures 4F–O). External verification of the GRIP signature using the GEO dataset showed that the model had a fair ability to predict patient prognosis with an AUC of 0.638, 0.660, and 0.659 in predicting 1-, 2-, and 3-year ROC curve (Figures 4P–T).

Univariable and multivariable Cox regressions showed that the risk score, age, T, and n were independent predictors for the OS of patients with CM (Figures 5A,B). The nomogram (Figure 5C) incorporating clinicopathological characteristics and the risk score was accurate in predicting the 1-, 3-, and 5-year OS rates of patients with CM (Figure 5D). A decision curve analysis proved that the signature and nomogram had a greater predictive ability than other traditional clinicopathological features for patients with CM net benefits (Figure 5E). The nomogram had an AUC of 0.810, showing a good predictive ability for the outcomes of patients with CM, which might be applied in the clinical management of patients with CM (Figures 5F,G). Compared with existing signatures for CM, the risk signature showed its performance better (Figure 5H).

The KEGG analysis indicated that RNA polymerase, aminoacyl tRNA biosynthesis, base excision repair, lysine degradation, pyrimidine metabolism, and oxidative phosphorylation were significantly activated in the high-risk group (Figure 6A), whereas the chemokine signaling pathway, apoptosis, JAK-STAT signaling pathway, natural killer cell, Toll-like receptor signaling pathway, cytokine-cytokine receptor interaction, T-cell receptor signaling pathway, and antigen processing and presentation were significantly activated in the low-risk group (Figure 6B). The Hallmark pathway analysis showed that DNA repair, E2F targets, MYC targets v1 and v2, and oxidative phosphorylation were upregulated in the high-risk group (Figure 6C), whereas IL2-STAT5 signaling, IL6-JAK-STAT3 signaling, inflammatory response, complement signaling, apoptosis, interferon α response, and TNFA signaling were upregulated in the low-risk group (Figure 6D). In addition, the inflammation pathway, T-cell receptor pathway, immune response, immune system process, inflammasomes, PD-1 signaling, and pyroptosis were upregulated in the low-risk group (Figure 6E). These results suggest that a pyroptosis-mediated inflammatory response upregulated the immune microenvironment and activated an inflammatory response in the low-risk group.

The analysis of immunocyte infiltration indicated that immune cell subpopulations were highly enriched in the low-risk group compared to the high-risk group (Figure 7A). The ssGSEA of immune cells and functions showed that, compared with the high-risk group, the low-risk group had higher levels of immune checkpoints, cytolysis, human leukocyte antigens, inflammation, T-cell co-stimulation, T-cell co-inhibition, CD8+ T cells, T helper cells, and tumor-infiltrating lymphocytes (TILs) (Figures 7B,C). Differences in these factors, especially immune checkpoints, CD8+ T cells, B cells, inflammation, TILs, NK cells, T helper cells, and immune cells functions, had important impacts on the survival of patients (Figures 7D–O).

A total of 36 immune checkpoints, including PD-L1, PD-1, and CTLA-4, were found to be differentially expressed between the low- and high-risk groups (Figure 8A). Compared with high-risk patients, low-risk patients had better responses to anti-PD-1 and anti-CTLA-4 immunotherapy (Figures 8B–E). The TIDE and dysfunction analyses confirmed that low-risk patients were more sensitive to immunotherapy (Figures 8F,G). The prognostic performance in immunotherapy cohorts showed that the patients with CM and a low-risk and high immunotherapy score had significant elevations in prognosis (Figures 8H–M).

The estimated IC50 values for 17 common targeted drugs, such as afatinib, sorafenib, and refametinib, differed significantly between the two risk groups (Figure 9A). In addition, the low- and high-risk groups had significantly different sensitivity to 12 common chemotherapeutic drugs, such as docetaxel, rapamycin, cisplatin, and DMOG, and differed significantly between the two risk groups (Figure 9B). These results suggest that the GRIP signature can regulate sensitivity to chemotherapeutic and targeted drugs and that the risk model might be used to identify the potential biomarkers for chemotherapy and targeted therapy sensitivity.

Immunohistochemical images of IFNGR2, CCL25, IL15, RTP4, and NLRP6 were obtained from The Human Protein Atlas database, and semiquantitative analysis confirmed the expression of these proteins in CM tissues (Figure 10). Then, the expression levels of the eight genes were validated in CM cell lines by qRT-PCR (Figure 11). The results showed that the expressions of TLR1, CCL8, EMP3, IFNGR2, CCL25, and IL15 were upregulated in A375 and SK-MEL-28 cell lines compared with Hacat and PIG1 cell lines. Although the expression of RTP4 and NLRP6 showed no statistical difference between PIG1 cells, A375 and SK-MEL-28 cell lines, the two genes were upregulated in A375 and SK-MEL-28 cell lines compared with Hacat cell.

Our findings demonstrate that pyroptosis has antitumor activity and improves patient survival. Recent studies have shown that the antitumor effects of pyroptosis mainly depend on a strong inflammatory response and that pyroptosis changes the surrounding TME by triggering the release of inflammatory mediators (57, 58). Thus, we constructed a co-expression network and performed WGCNA to identify GRIPs and investigate their role in CM. Overall, our analyses identified 134 differentially expressed GRIPs. Our GO and KEGG analyses further showed that differentially expressed GRIPs were enriched in T-cell activation, cytokine and chemokine activity, and inflammatory response-related signaling pathways, such as the cytokine-cytokine receptor interaction, TLR, TNF, JAK-STAT, nuclear factor (NF)-κB, and NOD-like receptor (NLR) signaling pathways. After univariable Cox and Lasso regression analysis, we also constructed a GRIP risk signature, which we verified in the TCGA and GEO cohorts, and a nomogram scoring system for predicting the prognosis of patients with CM. Next, In addition, our GSEA showed that pyroptosis and inflammatory response have potential antitumor effects in patients with CM. Furthermore, our results also provide evidence of a complicated relationship among CM, the pyroptosis-mediated inflammatory response, and immunocyte infiltration in the TME. Drug sensitivity analysis indicated that GRIP signature might be used to identify potential biomarkers for chemotherapy and targeted therapy sensitivity. The immunohistochemical images of IFNGR2, CCL25, IL15, RTP4, and NLRP6 in both normal skin tissue and melanoma tissue confirmed the expression of the GRIPs in CM. Finally, the qRT-PCR experiments validated the expressions of our risk genes in Hacat, PIG1, A375, and SK-MEL-28 cell lines.

Many kinds of pyroptosis or inflammation signatures have shown value in CM (32–36); however, there is almost no record about the genes both related to inflammation and pyroptosis models in CM prognosis and anticancer therapies. Our signature comprising eight GRIPs (TLR1, IFNGR2, CCL25, IL15, CCL8, NLRP6, EMP3, and RTP4) was superior to Lin‘s inflammation signature, an inflammatory response-related signature predicted the 3-year OS of patients with hepatocellular carcinoma with an AUC of 0.705 (59) and other researches' pyroptosis signature (10, 35) at 0.662 to 0.70. Our nomogram scoring system with an ROC curve of AUCs as 0.810 was higher than a recent study about signature nomogram in CM (60–65). In addition, all eight GRIPs we identified directly or indirectly play vital regulatory roles in molecular mechanisms. For example, Hu et al. have confirmed that the expression of TLR4 was related to prognosis and positively correlated with the infiltration of B cells, CD4 and CD8 T cells, neutrophils, macrophages, and dendritic cells in CM (66). Goto et al. showed that TLR2, TLR3, and TLR4 are highly expressed in human melanoma cells and that downstream signal transduction factors, such as NF-κB, and inflammatory response-related factors are activated in these cells (67). In another study, Sun et al. showed that NLR signaling pathways also play an important role in the inflammatory activity of the TLR4/NF-κB pathway and are associated with pyroptosis and the expression of downstream proinflammatory cytokines such as caspase-1, IL-1β, and IL-18 (68). The TLR4 we identified in GRIP signature was similar to those reported previously in patients with CM. In addition, IL15, CCL25, and CCL8 are associated with the inflammasome, which is a multiprotein complex whose assembly and activation are responsible for the recruitment and activation of caspases-1 and−5 that play an essential role in pyroptosis (68–70). EMP3 belongs to the peripheral myelin protein 22-kDa (PMP22) gene family as novel therapeutic targets in human cancer, which is the most important indicator of progression-free and metastasis-free survival for patients with urothelial carcinoma of the upper urinary tract (71). Therefore, our results suggested that the eight GRIPs might regulate pyroptosis by producing inflammation response. The functional enrichments we identified in CM were similar to those reported previously.

In the TME, pyroptosis and the inflammatory response may enhance immune surveillance and antitumor immunity by recruiting TILs such as CD4+ and CD8+ T cells (72–74), and high levels of immune cell infiltration are associated with favorable outcomes (75). Our results further provide evidence of a complicated relationship among CM, the pyroptosis-mediated inflammatory response, and immunocyte infiltration in the TME. For instance, our signature mainly had a significantly positive correlation with the infiltration of neutrophils, CD8+ T cells, CD4+ T cells, and dendritic cells, which also had high expression levels of pyroptosis-related genes. Some studies have indicated that the level of immune infiltration in the TME (particularly that of immune cell subpopulations such as CD8+ T cells) is a crucial factor in the assessment of melanoma prognosis (76). Ahmadzadeh et al. reported that CD8+ T cells stimulate granulocytes and “produce perforin” or “produce granulocyte colony-stimulating factor and perforin” to kill tumor cells and that melanoma with high CD8+ T cell infiltration is more likely to respond to anti-PD-1/PD-L1 therapy (77). Similarly, we found that the GRIP signature was also highly correlated with PD-1 and CTLA-4 expression and with enhanced ICI efficacy and chemotherapeutic and targeted drug sensitivity in patients with CM. Thus, regulating the expression of GRIPs may increase the efficacy of ICB in melanoma, and our results also revealed the potential value of the GRIP signature to enhance the effect of immunotherapeutic, chemotherapeutic, and targeted drugs.

Although this study presents encouraging results, there are still several limitations in our study. First, the model was constructed and validated in open data source (TCGA and GEO). It would be better if its prognostic value was tested in another independent patient cohort. Second, our study lacked validation of clinical samples. Future studies are required to investigate the underlying mechanisms of GRIPs in mediating CM progression and immune microenvironment. We will continue to work on this in further studies.