A Pyroptosis-Associated Signature Plays a Role in Prognosis Prediction And Indicating Inltration in Immune Microenvironment in Clear Cell Renal Cell Carcinoma

Background: About 90% renal malignancies are RCCs (renal cell carcinomas) and the main subtype in histology is ccRCC (clear cell RCC). In recent years, pyroptosis was considered as a kind of inammation-related programmed cell death, which played a part in the invasion, metastasis and proliferation of tumor cells, thereby inuencing tumor prognosis. Nonetheless, the expression level of pyroptosis-associated genes in RCCs and its relationship with prognosis are still not clear. Results: In our research, 44 regulators for pyroptosis were identied which were expressed differentially between normal kidney and ccRCC tissues. ccRCC cases were categorized into 2 subgroups with a signicant difference between them in OS (overall survival) according to the DEGs (differentially expressed genes). The prognostic value of pyroptosis-associated genes was assessed as a signature based on the cohort of TCGA (The Cancer Genome Atlas). Following Cox regression with DEGs and LASSO (least absolute shrinkage and selection operator), a 6-gene signature had been set up and all the ccRCC cases in the cohort of TCGA were categorized into LR (low-risk) group or HR (high-risk) group (P < 0.001). In combination with clinical features, risk scores were considered as a predicting factor of OS in ccRCC cases. KEGG (Kyoto Encylopedia of Genes and Genomes) and GO (Gene ontology) analyses implied increased immunity and enrichment of genes related to immunity in the HR group. Conclusions: Our ndings indicated that the genes related to pyroptosis had an important role in tumor immunity, which were potentially used to predict the prognosis of ccRCCs. Conceptualization, Z.-Y.L. and Z.-N.X.; methodology, Z.-Y.L. and Y.-P.Y.; software, Z.-Y.L. and B.-W.W.; validation, L.-C.C., W.-G.J. and T.-D.W.; formal analysis, Z.-Y.L. and Z.-N.X.; investigation, Y.-P.Y. and W.-X.; resources, Z.-Y.L. and Z.-N.X.; data curation, Z.-Y.L.; writing—original draft preparation, Z.-Y.L and X.-Y.W.; writing—review and editing, Z.-Y.L., C.-Y.W. and C.-Z.; visualization, Z.-Y.L., P.-Z. and W.-H.Y.; supervision, Z.-Y.L., C.-Y.W. and C.-Z.; project administration, C.-Y.W. and C.-Z.; funding acquisition, C.-Y.W. and C.-Z. All authors have read and agreed to the published version of the manuscript.


Introduction
RCCs (renal cell carcinomas) are derived from renal tubular epithelial cells, accounting for approximately ninety percentage of all kinds of renal malignancies. The main histological type is clear cell RCC (ccRCC) [1]. Recently, RCCs incidence in most countries has elevated, with 400,000 new cases every year around world and over 175,000 deaths. Its mortality rate and incidence rank 3rd among urological malignancies worldwide [2,3]. For regional renal cancer in early stage, radical nephrectomy is the rst-line treatment, even though distant metastasis or tumor recurrence still remains unsolved after operation in over 20% of patients [4]. Meanwhile, RCCs are characterized by tolerance to chemotherapy and radiotherapy.
Developing a speci c prognostic strategy is of great importance to improve therapeutic effects.
Pyroptosis is a novel way of programmed cell death, which is recognized as in ammation related necrosis of cells [5], which is induced by a variety of stimuli, including heart attack, bacterial or virus infections, cancer and stroke [6]. Except for autophagy, apoptosis, and ferroptosis, this kind of cell death has attracted much attention recently.
The characteristics of pyroptosis is cell swelling, large bubbles moving from the membrane, and lysis [7].
Caspase-1 [ICE for IL (interleukin)-converting enzyme], a member in the in ammatory caspase family, was the rst kind of caspases discovered involved in facilitating pro-IL-1b to form mature IL-1b [8,9]. Caspase-1-dependent plasma membrane pores require the gradients of cellular ions, contributing to an increase in the osmotic pressure, which results in cell swelling and water in ux [10]. Caspase-1 dependence, which regulates cell lysis, is a speci c characteristic in pyroptosis but not apoptosis [11][12][13]. The gasdermin family is a major executor in pyroptosis which contains pejvakin(PJVK or DFNB59) and gasdermin-A to gasdermin-E [14]. In ammasomes participate in the activation of Caspase-1, resulting in GSDMD (gasdermin D) cleavage and the maturation and secreting pro-in ammatory factors, including IL-1B and IL-18 [15]. Except for GSDMD, cleavage of gasdermin proteins of other kinds also induce formation of pores on plasma membrane. Particularly, Caspase-3 participates in the cleavage of GSDME (Gasdermin E) to induce pyroptosis [16,17].
According to the previous studies, it was discovered that pyroptosis had a pivotal part in the development of malignancies and antitumor activities. For instance, recent research found a novel gene signature related to pyroptosis which might be used in the prognosis of Skin Cutaneous Melanoma [18].
Nonetheless, the speci c effects of pyroptosis remain to be explored in ccRCC. Thus, this systematic study was conducted to explore the expression difference of genes related to pyroptosis in ccRCC and normal tissues. Meanwhile, a new PRGs prognostic risk signature in ccRCC has been established to predict survival. Furthermore, the phenotypes related to prognosis, and the relationships between the immune microenvironment of tumor and pyroptosis were investigated.

Collection of TCGA Data
The sequencing data of transcriptome RNA and the clinical data were downloaded (website: https://portal.gdc.cancer.gov/; TCGA database) with a total of 611 ccRCC cases (72 cases of normal samples, and 539 cases of tumor samples). Samples without intact information were excluded.
Identi cation of DEGs related to pyroptosis A total of 52 genes associated with pyroptosis were extracted from prior reviews [19][20][21][22][23][24][25][26] and shown in Table S1. In order to nd DEGs associated with pyroptosis between tumor and normal tissues, the"limma"package was applied. Version 11.0 STRING (Search Tool for the Retrieval of Interacting Genes; website: https://string-db.org/) was used to construct PPI network. R programming language and pheatmap package were utilized to acquire heatmaps of DEGs."igraph" and "reshape2" R packages were applied to evaluating the relationship between the selected DEGs.

Consensus Clustering
To categorize the ccRCC according to consensus clustering, "limma" and "ConsensusClusterPlus" R packages were applied. The relationships between clinical characteristics (OS, overall survival) and clusters were evaluated via using the R package "survival" and chi-square test. The KM (Kaplan-Meier) curves and heatmaps were used to present the results by using R packages "pheatmap", "survival", and "survminer".

Development of ccRCC Prognostic Models based on PRGs
The Cox regression analysis (least absolute shrinkage and selection operator, abbreviated as LASSO) was conducted by employing the R package "glmnet" to construct a model of prognosis using candidate genes. The minimum parameters were used to determine the penalty parameter (λ). The following equations were used to calculate the RS (risk score): RS equals to e Σi(Coe ·Expi) , where Expi stands for expression of each retained gene and Coe for coe cient. PCA (principal component analysis) was performed by applying R package"Rtsne" according to the risk score. The R packages "survival" and "survminer" were applied for KM analyses. And the 1-, 3-and 5-year ROC analyses were performed by using the R package "survivalROC".
Analysis on Prognostic Value of the RS The clinical features (stage age, grade, gender, and T and M classi cation) of the patients were obtained in the TCGA cohort. The RS and these extracted variables were analyzed together in the regression model, with multivariable and univariate Cox regression models used.

Analysis of Gene Sets Enrichment
In the TCGA cohort, patients with ccRCC were assigned into 2 subgroups based on the median RS. DEGs between the HR (high-risk) group and the LR (low-risk) group were screened based on the criteria (FDR< 0.01 and |log2FC| ≥1.5). According to the DEGs, KEGG and GO analyses were conducted via using the "clusterPro ler" package. The ssGSEA was conducted by using "pheatmap package", "gsva" package to assess activities of the pathways linked to the immune responses after determination of the scores of IIC (in ltration immune cells).

Analysis of TIMER Database
Comprehensive analysis is conducted using TIMER database (website: https://cistrome.shinyapps.io/timer) and visualization of IIC in over ten thousand tumors of thirty-two types of cancers is facilitated [27]. Six subsets of tumor-in ltrating immune cells (macrophages, dendritic cells, CD8T cells, neutrophils, CD4T cells, and B cells) were contained in TIMER. A new statistical approach was used to assess abundance of the 6 types of IIC in the tumor microenvironment. The genomic, immunological, and clinical features of the tumor were comprehensively investigated by using the TIMER database. In ccRCC, for each hub gene in the risk score model, SCNA (Somatic Copy Number Alterations) module of TIMER tool was used to compare the in ltration between samples containing different SCNA, including high ampli cation, diploid/normal, deep deletion, arm-level deletion, and armlevel gain [28]. Additionally, the in ltration level in ccRCC samples was harvested using the TIMER database for determination of the relationship with the system of RS. Statistics R software (version 4.1.0) and the packages mentioned above were used for statistical analysis. The logrank test and K-M method were used to perform survival analyses. The signi cance of prognostic factors was evaluated by using multivariate and univariate Cox regression analyses. In subgroup differential analyses, Kruskal-Wallis test and Wilcoxon rank-sum test were applied. All statistical tests were twosided. P < 0.05 was regarded as statistically signi cant. To further investigate the interactions between these PRGs, PPI (protein-protein interaction) analysis was conducted via applying the STRING platform. The result was illustrated in Fig. 1B. The minimum interaction score required for PPI analysis was 0.9 (the highest con dence). According to the number of nodes, the top 30 hub genes were listed in Fig. 1C. Moreover, except for CHMP4B and HMGB1, DEGs between normal tissues and tumor tissues were shown. The network of correlation including all the genes linked to pyroptosis was demonstrated in Fig. 1D (blue color, negative; red color, positive).

Tumor categorization according to the DEGs linked to pyroptosis
In order to investigate the relationship of the ccRCC subtypes with the 44 DEGs levels, CCA (consensus clustering analysis) was performed by using the data of 539 ccRCC (Table 1) samples in the cohort of TCGA. When the clustering variable (k) was raised from two to ten, the intragroup correlation reached the peak. Moreover, the intergroup correlation was low when k = 2, indicating that 539 ccRCC patients could be well divided into two groups according to 44 DEGs, as shown in Fig. 2A. The clinical characteristics and gene expression pro le are presented in a heatmap. After comparing the cluster 1 with cluster 2, we observed that there were signi cant differences in clinical stage (stage I-IV), distant metastases (M0-1), Fuhrman grade (G1-G4), and tumor size (T1-T3), whereas the differences in the number of LNM (lymph node metastases), gender and age were not statistically signi cant (Fig. 2B). We also compared OS of the 2 clusters, which was presented in Fig. 2C. OS was signi cantly poorer in cluster 1 than cluster 2(P<0.001, Fig. 2C).

Prognostic Value of PRGs Expression Signature in ccRCC
There were 530 ccRCC samples matched with corresponding cases with intact information about survival. We used univariate Cox regression analysis to initially screen for survival-related PRGs, and 8 genes with P < 0.05 were used for analysis (Fig. 3A). To screen candidate genes for construction of the model of prognosis, Cox regression analysis (LASSO) was performed. Based on the optimal λ value, 6 genes ( Fig. 3B, C) and their coe cients (Table 2) were eventually preserved. According to the median score, 530 patients were equally assigned into high-and low-risk subgroups (Fig. 3D). As shown in Fig.  3E, PCA demonstrated that subjects with alternative risks were divided into 2 clusters. The death rate was higher and survival time shorter in the high-risk group compared with the low-risk group (as shown in Fig.   3F). As shown in Figure 3G, there was signi cant difference in ccRCC. In the high-risk subgroup, life span was shorter compared with low-risk subgroup(P<0.001). The speci city and sensitivity of the prognostic model were assessed by applying ROC (receiver operating characteristic) curve, and AUC (area under the ROC curve) was 0.706, 0.640 and 0.720 for 5-year, 3-year, and 1-year survival, respectively (Fig. 3H). To identify the prognostic value of the risk model, we conducted Kaplan-Meier analysis to con rm whether these genes involved in the construction of risk model were associated with the prognosis of ccRCC. In Fig. 4A, higher levels of ELANE, AIM2, GSDMB, IL6, NLRP1, and NOD2 were positively correlated with poor prognosis. We then analyzed if RS of the gene signature model was a factor predicting prognosis. The results suggested that RS in the cohort of TCGA was a predictive factor of prognosis (HR=3.065, 95% CI: 2.295-4.092, Fig. 4B). According to the results of multivariate analysis, after adjustment of confounding factors, RS was a predictive factor for the prognosis of patients with ccRCC in the cohort of TCGA (HR=2.251, 95% CI: 1.659-3.055, Fig. 4C). Additionally, for the cohort of TCGA, a heatmap of clinical characteristics was generated (Fig. 4D). It was discovered that there was signi cant difference in the tumor grade and stage between the subgroups. These ndings suggested that the model of prognosis based on PRGs was robust and independent in the prediction of ccRCC prognosis.

Identifying the Prognostic Model-Related BP (Biological Processes)
It was of great importance to nd out what BP were affected by the model of prognostic risk in the role of prediction. The"limma" R package was utilized to obtain the differentially expressed genes: FDR < 0.01, |log2FC |≥1.5. As a result, a total of 381 DEGs in the cohort of TCGA were selected.378 of 381 genes in the high-risk group were over-expressed, and reduced expression was indicated in the remaining 3 genes (Table S3). KEGG pathway analysis and GO enrichment analysis were conducted by using the selected DEGs (Fig. 5). Of interest, the most enriched biological processes were closely related to immune responses, in ammatory cell chemotaxis, and chemokine-mediated signalling pathways, including B cell regulated immunity, responses of humoral immune, CXCR chemokine receptor binding, and TNF/ NF−kappa B/ IL−17 pathways. These ndings suggested that the model of prognostic risk based on PRGs was associated with immune responses.

Evaluation of Immune Cell In ltration between subgroups
According to the above ndings, it was proposed that the functions of PRGs in the prediction of ccRCC prognosis might be associated with the immune microenvironment. The relationship was evaluated between the in ltration of immune cells and prognosis-related genes. The changes of in ltration were explored by using the samples with copy number alteration of ELANE, AIM2, GSDMB, IL6, NIRP1, and NOD2, respectively. The results showed that copy number of these prognosis-related PRGs was associated with the immune microenvironment in ccRCC. Mutants of AIM2, ELANE, GSDMB, NLRP1, and NOD2 inhibited the in ltration of some kinds of immune cells (Fig. 6A-F). Then, the immune estimation dataset was downloaded from the TIMER database and the correlation was analyzed between the in ltration of 6 types of immune cell (macrophages, dendritic cells, neutrophils, B cells, CD 8 +T cells, and CD 4+T cells) and the risk score of the prognostic model. The data suggested that the risk score of the PRGs-related prognosis model was positively correlated with immune in ltration (Fig. 7A-F).
The activities of thirteen pathways related to immune responses and enrichment scores of sixteen kinds of immune cells were then compared between the HR group and the LR group in the cohort of TCGA using ssGSEA (single-sample gene set enrichment analysis). The results revealed that in ltration was increased in the HR group, especially DCs, Th cells (Th2, Tfh, and Th1 cells), CD8+ T cells, Macrophages, B cells, Treg cells, pDCs, TILs (tumor-in ltration lymphocytes) and Neutrophils, compared to the LR group, as shown in Fig. 8A. In Fig. 8B, it was illustrated that all the13 immune pathways in the high-risk group were more active compared with the low-risk group.
The study design and grouping are shown in Fig. 9.

Discussion
Pyroptosis is a newly recognized kind of programmed cell death, which has a double-edged role in the progression of malignancy and mechanisms of treatment. In ammatory cytokines are secreted in pyroptosis and normal cells are stimulated, contributing to the progression of cancer [19]. Meanwhile, pyroptosis can enhance cellular death in malignancy, which makes pyrolysis a possible therapeutic and prognostic target for malignancy [29]. However, how the genes related to pyroptosis interact with each other and whether they have impact on the survival time in ccRCC remain unrevealed.
In this research, the mRNA expression of 52 currently known PRGs were rst discovered in normal and ccRCC tissues. It was found that 44 of them were differentially expressed. What makes sense is that the 2 clusters in the consensus clustering analysis according to the pyroptosis-related DEGs showed signi cant differences in clinical characteristics. This implied that pyroptosis in cancer tissues was different in ccRCC, which resulted in different OS.
To further evaluate the prognostic value of the regulators related to pyroptosis, a 6-gene risk signature was constructed. OS in cases of alternative subgroups was different. The functional analyses presented that the DEGs between the high-and low-risk groups showed signi cant differences in immune-related pathways, consistent with our expectations. Pyroptosis could contribute to the accumulation of various in ammatory factors, which was also resulted from the activation of in ammasomes [6,30]. Additionally, we also enriched Cytokine, NF−kappa B, IL-17, Ras and TNF signaling pathways which were closely associated with the development of RCC.
Meanwhile, the pathway activation and immune in ltration were compared between the HR group and the LR group, and the HR group showed increased activities of pathways related to immune responses and elevated levels of IIC in comparison with the LR group. Enhanced in ltration of immune cells was related to poor prognosis, which was consistent with previous studies [31]. Another key result in our research was that the above 6 prognostic genes related to pyroptosis had a signi cant correlation with immune in ltration, which further proved the fact that pyroptosis had an essential part in the immune microenvironment of tumor.
The current study found a signature of 6 genes linked to pyroptosis (ELANE, AIM2, GSDMB, IL6, NIRP1, and NOD2) and discovered that it could play a role in predicting OS of ccRCC patients. AIM2 (Absent in melanoma 2) was rst found in melanoma with decreased expression [21]. AIM2 is a member of the IFNinducible PYHIN (pyrin and HIN200 domain-containing) family, and works as a cytoplasmic senor for DNA which can connect with dsDNA (double-stranded DNA) [32]. AIM2 can activate CASP-1 via junctional proteins regulated by ASC to facilitate secretion and maturation of IL-18 and IL-1β, thus promoting pyroptosis [33]. Previous research demonstrated that AIM2 worked as a suppressor in multiple kinds of tumors such as prostate cancer [34], colon cancer [35], melanoma [36], melanoma [37], and breast cancer but as a promoter in NSCLC (non-small cell lung cancer) [38], OSCC (oral squamous cell carcinoma) [39], and HPV (human papillomavirus)-associated cervical cancer [40]. Therefore, AIM2 might had different effects in different tumors. In this research, the expression of AIM2 was signi cantly increased in tumor tissues in comparison with normal ones. Furthermore, an increased expression level of AIM2 was closely related to poor survival and gene mutations in AIM2 could ameliorate in ltration of immune cells. Thus, it was suggested that AIM2 was more like a pro-oncogene. The molecular mechanisms of AIM2 in the development of ccRCC currently remained unrevealed and our ndings in AIM2 might provide a new insight into further research.
GSDMB (Gasdermin B) was the gene most linked to RS in the prognostic model, prompting that GSDMB might be highly involved in ccRCC. According to previous research on human malignancies, it was found that GSDMB was upregulated in tumor tissues, including breast, uterine, gastric, and cervical cancers [41].
It was demonstrated that GSDMB was located in amplicons, the genomic regions frequently ampli ed in cancer development [42]. Thus, GSDMB might participate in cancer development and metastasis.
GSDMB can be cleaved into two fragments by caspasa-1. One cleavage form is the N-terminal of GSDMB protein with a molecular weight of 20 kDa. Cell pyroptosis can be caused by the secretion of the Nterminal domain. On the contrary, the full-length N-terminal domain and C-terminal fragment would not cause cell pyroptosis [43,44]. In general, GSDMB may be the downstream protein of pyroptosis pathway.
The key is whether some factors can trigger the upstream mechanism of GSDMB and cause pyroptosis.
At the same time, the speci c mechanism of GSDMB in ccRCC is not clear. The increased level of GSDMB in ccRCC is related to poor prognosis. This nding may provide help for the study of tumor treatment targets.
ELANE (neutrophil elastase gene) is a main serine protease produced by neutrophils, which can activate proin ammatory cytokines including IL-1β, IL-18, and TNF-α [45,46], which are regarded as promoters for pyroptosis. Kambara et al. [47] presented that ELANE could cleave and activate GSDMD and subsequently induce pyroptosis in neutrophils. The ELANE expression level in the HR group was remarkably increased compared with the LR group, but paradoxically, the neutrophil in ltration score was incredibly higher compared with the LR group. These ndings may be resulted from many complex factors driving the difference of gene expression between alternative tissues, particularly levels of genes and PRGs linked to in ammation, such as the proportion of in ltrated immune cells and the differentiation of ccRCC [48,49].
Such factors might not in uence application of PRG expression signature in diagnosing and predicting ccRCC prognosis. The relationship between the expression of pyroptosis genes and in ltration of immune cells, ccRCC differentiation status, and other factors requires further investigation, which may provide new enlightenment for predicting the diagnosis and prognosis of ccRCC.
NLRP1 (NLR family, pyrin domain containing1), a bipartite adaptor protein, is considered as apoptosisassociated speck-like protein with a ASC (caspase-recruitment domain). NLRP1 can promote the recruitment process of pro-caspase-1 to the complex of in ammasomes [50]. GSDMD can be cleaved by active caspase-1, allowing the N-terminal domain of GSDMD to form pores in the plasma membrane, and then triggers the pyroptosis mechanism [10,15,51,52]. NOD2 (nucleotide-binding oligomerization domain-containing protein 2), which can initiate NF-κB (nuclear factor-κB)-dependent and MAPK (mitogen-activated protein kinase)-dependent gene transcription. For macrophages, NOD2 could be used as its promoter to induce the activation of in ammatory bodies [53]. In our study, NOD2 was highly expressed in the HR group, which might be one of the reasons why macrophages obtained higher scores in the HR group compared with the LR group. At the same time, the activation effect of NOD2 on NF-κB pathway might make this pathway more enriched in KEGG analysis. Interestingly, as one of the genes in the risk prognostic model, activated NF-kB pathway was able to increase the amount of mRNA for IL-6 [54]. In advanced metastatic breast cancer cells, over-activated NF-kB promoted chromatin accessibility of the IL-6 promoter region and enhanced transcription of the IL-6 gene [55]. Appropriate expression of IL-6 is of great importance for human immune defense, but sustained production of IL-6 has a pivotal role in the occurrence of multiple in ammation-related diseases and cancer [56]. Whereas in our study, IL-6 was signi cantly highly expressed in the HR group, we speculated that NOD2 might act as an upstream initiator of IL-6 to promote IL-6 expression by activating the NF-kB pathway.
The purpose of our research was to categorize cases with ccRCC into different subtypes, screen DEGs, set up a model of prognosis, and connect pyroptosis with the prognosis. Though we conducted multi-angle and multi-omics validation, there were still limitations in this study. All the analyses were performed using the TCGA KIRC cohort, and GEO cohort was recommended to verify them. Additionally, in-vivo and in-vitro experiments are required to further verify our ndings. Pyroptosis, particularly the mechanism in ccRCC, is not su ciently investigated. We initially investigated the prognostic value of the 6 genes related to pyroptosis in the risk prognostic model, which provided theoretical support for future studies.

Conclusions
In conclusion, a comprehensive and systematic bioinformatics analysis was performed and we found a prognostic gene signature related to pyroptosis including six genes (ELANE, AIM2, GSDMB, IL6, NIRP1, and NOD2) for ccRCC patients. Moreover, the risk score in the prognostic model according to the 6 PRGs was an independent prognostic factor of cases with ccRCC, which was associated with the immune microenvironment. Availability of data and materials

Declarations
In the study, diferent web-based datasets were used for data analysis. The web links to all the original data sources were listed as below: The sequencing data of transcriptome RNA and the clinical data were downloaded from The Cancer Genome Atlas Program (TCGA) (https://portal.gdc.cancer.gov/) data portal. All data generated from the analysis process of this study are available from the corresponding author on reasonable request.      In ltration levels of immune cells of the six PRGs mutants linked to prognosis. A AIM2, B ELANE, C GSDMB, D IL6, E NLRP1, F NOD2. ***p < 0.001, **p < 0.01, *p <0.05. Figure 8 ssGSEA scores concerning immune pathways and cells. A Enrichment scores of sixteen types of immune cells between LR group (blue box) and HR group (red box). B Enrichment scores of thirteen pathways linked to immunity between high-(red box) and low-(blue box) risk groups. ***P < 0.001, **P< 0.01, *P < 0.05.