The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/repository) was used to gather RNA sequencing data and associated clinical information for 374 patients with HCC and 50 normal samples. We excluded 4 patients who did not have complete follow-up information. In addition, transcriptomic counts and FPKM data with clinical features from 231 patients with HCC were obtained from the ICGC database (https://dcc.icgc.org/projects/LIRI-JP). Standardization for gene expression was performed using the “limma” tool.

Mutation data corresponding to 370 patients with HCC were from TCGA, and four incomplete information patients were excluded. The remaining patients were divided into two risk groups and the mutational spectrum was analyzed using the MAftools package, with measurement of tumor mutational burden (TMB) [17].

Data for CNV corresponding to 370 HCC patients also were from TCGA. The lolipop plot was generated to show the gain and loss function of genes in the HCC database.

The STRING database (http://string-db.org) provides the integration of protein interactions, including direct and indirect associations. The STRING database was used to show the correlation between the PRGs.

There were 52 PRGs identified in the earlier literature. 42 PRGs showed differential expression between tumors and normal tissues. The “limma” R package was used to analyze the differentially expressed genes (DEGs) in the TCGA cohort and a p value of 0.05 was used to construct a heatmap or volcano plot paragraph. In addition, a network of protein-protein interactions was also created by the tool GENEMANIA (http://genemania.org/) for the PRGs.

TCGA-LIHC were divided into two groups of the 42 DEGs by using the “ConsensusClusterPlus” R program [18].

A Univariate Cox analysis was performed to eliminate nine survival related PRGs for further investigation. The prognostic model was constructed using the LASSO-Cox method and multivariate Cox analysis to reduce overfitting. Three genes were ultimately retained with their coefficients, and the penalty parameter was calculated using the minimum criteria. The Equation 1: risk score $=(0.311\times BAK1\text{exp}.)+(1.143\times GSDMC\text{exp}.)+(-0.402\times NLRP6\text{exp}.)$ was used to calculate the risk score.

The prognostic model was based on dividing patients with HCC into two risk clusters according to their median risk score. We then used the “limma” and “scatterplot3d” packages in R to generate a principal component analysis (PCA) for both risk groups on the basis of gene expressions. The R package “survminer” was used for survival analysis comparing the prognosis between the two risk groups, and the R packages “survival” and “timeROC” were used to construct the ROC curve. A signature of three genes was constructed using the LASSO-Cox method and the multivariate Cox regression analysis, and then to assess efficacy in the patients from the TCGA test or ICGC external validation cohort. The pyroptosis-related gene (PRG) mRNA level and the risk score was normalized using the “scale” function. Using the risk score, patients in the TCGA test or ICGC cohorts could be divided into two risk groups.

The GSEA 4.0.1 (Broad institute, Boston, MA, USA) was utilized to carry out GSEA analysis between two risk groups to determine the differential and enrichment studies performed in the “clusterProfiler” R package [19].

To examine the correlation between immune cell infiltration and PRG scores, 28 infiltrating immune cells in each case of HCC were calculated using the ssGSEA algorithm [20]. In addition, the “GSVA” R package was also used to perform a ssGSEA using gene-specific markers from the TCGA cohort in order to determine immune function between the PRG subsets. The Spearman correlation test was conducted to determine whether the PRG score was significantly correlated with immunosuppressive molecules or immune scores.

The tissue samples were approved by the Zhongshan Hospital. The patients signed an informed consent form before undergoing surgical resection. Protein was extracted from cancer and adjacent tissues from 12 patients with HCC for validation by western blotting.

Huh7 and Hep3B cells were cultured in complete DMEM medium (D5796, sigma) containing 10% FBS (16140071, Gibco) and 1% penicillin/streptomycin (3810-74, Sigma) and used in the present study. IFN-$\gamma$ (CST, 13684S) was used for PD-L1stimulation. All the cell lines used in our study is negative for mycoplasma (detected by Betotime, C0296). The Hep3B and Huh7 cell lines were purchased from American type culture collection (ATCC, Manzas, Virginia,US), which was already authenticated by Short Tandem Repeat (STR) profiling for excluding cell cross contamination. And these cell lines have not been reported to be contaminated or misidentified in ExPASy Cellosaurus databases.

Cells were lysed for extract total protein extraction and pelleted at 95 °C for 8 min. Electrophoretic separation was performed by 10% sodium dodecyl sulfate - polyacrylamide gel electrophoresis (SDS-PAGE) at a 100 V constant voltage, and at 330 mA, protein was transferred to polyvinylidene difluoride (PVDF) membrane and then incubated with 1:2000 diluted GSDMC (ab225635, Abcam, Cambridge, UK), PD-L1 (CST, 13684S, Boston, MA, USA), caspase-8 (ab119809, Abcam, Cambridge, UK) and Actin (AA128, Beyotime, Shanghai, China) antibodies for 12 h at 4 °C. The PVDF membrane was incubated with a 1:10000 diluted secondary antibody. Finally, chemiluminescence was performed.

The cDNA or shRNA (Genepharma, Shanghai, China) targeting CD274 or GSDMC was recombined into lentiviral vectors to overexpress or knockdown GSDMC. The recombinant plasmid was transfected into 293T cells, and the mature infectious lentivirus was collected after three days. Stable CD274 or GSDMC-overexpressing Hep3B and CD274 or GSDMC-knockdown Huh7 cell lines were constructed and verified by western blotting.

Cells were seeded in 96-well plates (10${}^3$ cells/well). After 8 h, CCK-8 (C0037, Beyotime, Shanghai, China) was added to the wells on day 1. Plates were placed in an incubator at 37 °C for 30 min, and cell absorbance was measured using a microplate readout reader, measurements were performed for 6 consecutive days.

The 6-well plates were used to seed the HCC cells (800 cells/well), and then were placed at 37 °C incubator for 2 weeks. After the culture medium was discarded, the cells were fixed with 4% paraformaldehyde, then the cells were stained using crystal violet (C0121, Beyotime, Shanghai, China), washed with phosphate buffered saline (PBS) (P1020, Solarbio, Beijing, China), and dried.

After 10 h, HCC cells were seeded in 6-well plates to form a dense monolayer. Lines were drawn with the tip of a pipette size 20 µL pipette perpendicular to the cell layer to form straight cell wounds. The width of the cell wound width was then washed with PBS and recorded under a microscope (Olympus, Japan) by taking photograph. Plates were placed at 37 °C incubator for 48 h, then the wound width was recorded again.

About 100 µL BD Matrigel mixture (diluted 1:5 with DMEM) was pre-coated in a chamber (3513, Corning, Corning, NY, USA) at 37 °C incubator for 5 h. DMEM without serum was used to dilute the HCC cells and 4 $\times$ 10${}^4$ cells were plated in the upper chamber. Then, 600 µL complete DMEM was added to the bottom chamber. Cells from the lower chamber were then incubated for 48 h at 37 °C, fixed, stained with crystal violet, washed with PBS, and photographed.

The DEGs between HCC and normal liver tissue were analyzed by a one-way ANOVA. The Mann-Whitney test was used to calculate the ssGSEA scores using the adjusted p-value from the BH method. The whole statistical investigations were performed using R 4.0.1 (Auckland, New Zealand). A p 0.05 was used to define statistical significance.

Fig. 1 shows the detailed protocol. A total of 52 PRGs selected from the KEGG database were included in this study. To find gene alterations of PRGs in patients with HCC, we downloaded data from the TCGA-LIHC databases. At the genetic level, 157 of the 364 samples (43.13%) showed mutations in regulatory genes linked to pyroptosis. Of these, TP53 had the highest mutation frequency (30%) followed by NLRP2 (2%). No mutations were found in other PRGs such as BAK1, CASP1, CASP5, CHMP2A, CHMP3, CHMP4A, CHMP4B, CYCS, GSDME, IL1A, CASP6, CASP9, GPX4, GSDMA, PJVK as well as SCAF11 (Fig. 2a). Our results indicated that most PRGs, except TP53, have no significant changes in gene variation.

Fig. 1.

Flow chart of the whole study. LIHC, liver hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; GSEA, gene set enrichment analysis; GSDMC, Gasdermin C.

Fig. 2.

Pyroptosis-related genes (PRGs) Screening of patients with hepatocellular carcinoma (HCC). (a) The mutation frequency of PRGs in 370 HCC patients. Each column represented a patient. The proportion of each variant type was showed by right-hand bar plots, while stacked bar plots show the conversion fractions in each sample. (b) TCGA-LIHC cohort’s CNV variation frequency for pyroptosis-related genes. A column’s height represents the frequency of alteration. Green dot: deletion frequency; red dot: amplification frequency. (c) Based on the LIHC cohort, the Copy Number Variation (CNV) alteration location of PRGs is shown. (d) The differential expression of PRGs were showed by heatmap. (e) The differential PRGs expressed showed by box plot. (f) Protein-protein interaction network. (g) Gene association network (Blueline: negative correlations; Redline: positive correlations. The depth of color reflects the relevance’s strength).

We also found copy number variation (CNV) alteration in 7 of the 52 PRGs, which showed that CASP9, HMGB1, and TP53 had a significant frequency of CNV deletions (Fig. 2b), whereas GSDMC, AIM2, GSDMD, and CHMP6 had a significant frequency of CNV amplifications (Fig. 2b). Fig. 2c showed the CNV alterations location on the chromosomes.

The expression levels of 52 PRGs were compared based on 50 normal samples and 374 HCC samples in TCGA data. Forty-two DEGs were identified (all p 0.05, Table 1), among which 39 genes were enriched in the tumor cluster, while three genes (IL-6, IL-1B and NLRP3) were down-regulated. Fig. 2d showed the heatmaps of RNA levels (blue: low expression; red: high expression) and the histograms in Fig. 2e.

We then utilized the STRING platform to analyze the network of protein-protein interactions of 52 PRGs’ regulators (Fig. 2f). The PRGs correlation network is shown in Fig. 2g (blue: negative correlations; red: positive correlations).

Our next step was to use consensus cluster analysis to explore the association of all 370 HCC subtypes in the TCGA cohort with these prognostic genes (Supplementary Fig. 1a). When κ = 2, the intragroup correlation was the smallest (Supplementary Fig. 1b,c), indicating that it was appropriate to split 370 patients into two clusters on the basis of 42 DEGs obtained previously (Cluster 1 = 160, Cluster 2 = 210). Then we found that patients in cluster 1 survived for a shorter period than those in cluster 2 (p 0.05, Supplementary Fig. 1d), indicating that these DEGs may indicate a prognosis. The baseline clinical features and gene expression profile of the two clusters are shown as a heatmap. However, no notable baseline differences were observed between the two groups (Supplementary Fig. 1e).

Through univariate Cox regression analysis, we identified nine genes (BAK1, BAX, CHMP4B, CASP8, GSDMC, NLRP6, NOD1, NOD2, and PLCG1). The curves in Fig. 3b represent the trajectories of each independent variable coefficient. The two dashed lines in Fig. 3c indicate two special $\lambda$ values: lambda. min and lambda.1se. Based on LASSON Cox analysis (Fig. 3b,c), lambda. min was chosen to construct a nine genes PRGs model. Finally, we obtained three key prognosis-related PRGs (BAK1, GSDMC, and NLRP6) by multivariate regression analysis with a cut-off of 0.1, and the other six genes were narrowed with a p-value $>$ 0.1. A total of 370 HCC specimens were randomly separated into two cohorts (training cohort, n = 186; test cohort, n = 184). In the training set, we found that three genes (BAK1, GSDMC, and NLRP6) were closely associated with HCC occurrence after performing a multivariate Cox regression analysis of PRGs (Table 2). We classified TCGA training cohort patients into two clusters (high- or low-risk clusters) according to the median risk score calculated by Equation 1 (Fig. 4a). The PCA demonstrated that patients in the different risk groups were well separated into two clusters (Fig. 4b). Patients in the high-risk clusters had a worse prognosis than those in low-risk clusters (Fig. 4c). The Kaplan–Meier curve (KM) curve also showed a notable difference in OS between the two risk groups (p = 0.002, Fig. 4d). Receiver operating characteristic (ROC) analysis was performed to evaluate the specificity and sensitivity of the prognostic model. We demonstrated that the 3-gene risk model could correctly indicate the clinical outcomes of HCC patients (AUC 0.58 for 1-year, 0.658 for 3-year, and 0.656 for 5-year survival) (Fig. 4e).

Fig. 3.

TCGA-LIHC risk signature construction. (a) Univariate analysis showed nine DEGs (p 0.05). (b) LASSO analyses of all nine genes. (c) LASSO regression parameter selection by cross-validation. (d) Three DEGs with p 0.1 based on multivariate analysis.

Fig. 4.

Validation in TCGA training cohort. (a–c) Patients were classified as baseline on the basis of their risk scores in the training cohort. (d) The Kaplan–Meier curve (KM) survival curve. (e) Analysis of time-dependent Receiver operating characteristic (ROC) scores.

We then assessed the predictive power of the risk signature in both the TCGA test cohort (n = 184) and the external cohort (ICGC, n = 231). In the TCGA test cohort, the high-risk (n = 83) and low-risk (n = 101) subgroups were well differentiated into two clusters using PCA (Supplementary Fig. 2a–c). A notably different OS was demonstrated by the KM survival curves of the two groups (p = 0.022) (Supplementary Fig. 2d). Based on the risk score, the ICGC cohort was classified into high-risk (n = 115) and low-risk (n = 116) subgroups using PCA (Supplementary Fig. 3a–c). There was a marked difference in the OS between the two groups, as represented by the KM survival curves (p 0.01) (Supplementary Fig. 3d). The AUC values for survival at 1-year, 3-year, and 5-year survival were 0.725, 0.745, and 0.580, respectively (Supplementary Fig. 3e). Taken together, these data revealed that this risk model consisting of three genes could predict the prognosis of HCC patients well, whether in the TCGA training cohort, the test cohort, or the ICGC cohort.

Cox regression analysis was used to assess whether the risk score could independently predict OS. The findings suggested that the risk score may be a promising hazard factor both in the univariate Cox regression analysis (p 0.001, HR = 2.280, 95 percent CI: 1.572–3.307, Fig. 5a) and in the multivariate analysis (p 0.001, HR = 2.171, 95 percent CI: 1.389–3.396, Fig. 5b).

Fig. 5.

Independent prognostic value assessment of risk scores. (a,b) Univariate and multivariate analysis in training cohorts. (c,d) Test cohorts by univariate and multivariate analysis. (e) The correlation of gene heatmap with between the different risk groups.

To assess the predictive power of the risk score in HCC, we performed univariate and multivariate Cox regression analyses in the TCGA test cohort of 184 patients (Fig. 5c,d). Similar to the results in previous cohorts, the risk score was proved to be a strong prognostic predictor (p 0.001, HR = 3.892, 95 percent CI: 2.097–7.226 for univariate analysis, Fig. 5c; p 0.001, HR = 3.496, 95 percent CI:1.878–6.508 for multivariate, Fig. 5d). Individuals in the high-risk group tended to have higher levels of BAK1 and GSDMC but lower level of NLRP6, as shown in the heatmap (Fig. 5e).

GSEA was used to explore the role of PRGs in HCC. Most of the genes were enriched in the apoptosis, IL-6-JAK-STAT3, P53, and TGF-$\beta$ signaling pathways, etc. (Supplementary Fig. 4a), suggesting that genes related to pyroptosis may participate in the oncogenesis and development of HCC through these pathways.

On the basis of the functional analyses, GSEA analysis was performed to compare 16 enriched immune cell types between the two risk clusters in the TCGA cohorts. In the TCGA cohort, CD56 bright cells, CD4 T cells, dendritic cells (DC), NK cells, Plasmacyloid DC, CD56 dim cells, and Memory CD4 T cells showed a positive correlation, while eosinophils, neutrophils, T helper (Th)1 cells, macrophages, mast cells, memory B cells, activated CD8 T cells, Th17, activated B cells and monocytes were negatively correlated with the risk signature (Supplementary Fig. 4b).

We explored changes in the protein levels of three PRGs using the Human Protein Atlas (HPA) database. The results indicated that BAK1 and GSDMC were upregulated in tumor tissues, whereas NLRP6 expression was lower level in HCC than in normal tissues (Supplementary Fig. 5).

PD-L1 is an immune inhibitor, that can upregulates GSDMC expression to transform apoptosis to pyroptosis in breast cancer [21]. Using the TCGA-LIHC database, we found that PD-L1 and GSDMC were positively correlated (Cor = 0.3, p 0.001) at the mRNA level in HCC (Fig. 6a). Since PD-L1 protein is seldom expressed in Huh7 or Hep3B cells without IFN-$\gamma$ stimulation, we detected GSDMC expression in IFN-$\gamma$-treated hepatoma cells. Immunoblotting assays showed that PD-L1 knockdown decreased GSDMC levels in Huh7 cells (Fig. 6b, left panel). In contrast, PD-L1 overexpression enhanced GSDMC levels in Hep3B cells (Fig. 6b, right panel).

Fig. 6.

GSDMC is upregulated by CD274 (PD-L1). (a) CD274 showed positive correlation with GSDMC (Cor = 0.3, p 0.01). (b) The full length of GSDMC, GSDMC-N and caspase-8 expression levels in PD-L1- knockdown Huh7 (sh1 and sh2) or PD-L1-overexpressing Hep3B cell (OE) lines, representative results of 3 experiments. To induce PD-L1 expression, we used IFN-$\gamma$ as a stimulator.

A previous study has reported that GSDMC was increased in lung adenocarcinoma, colon adenocarcinoma and breast cancer and correlated with reduced overall survival [22, 23]. However, the role of GSDMC in HCC has not yet been investigated. We wondered whether GSDMC is upregulated in HCC and contributes to the malignant biological behaviors of tumors. Twelve pairs of HCC samples and adjacent tissues were collected for western blot analysis. Compared with adjacent liver tissues, GSDMC was significantly upregulated in HCC tissues (Fig. 7a). Then stable GSDMC-overexpressing Hep3B (Fig. 7b) cell line was successfully constructed. The proliferation ability was enhanced in GSDMC-overexpressing Hep3B cells, as shown by the CCK-8 proliferation assay and colony formation assay (Fig. 7c,d). Wound healing and invasion experiments demonstrated that the cell motility and invasion ability were substantially enhanced in GSDMC-overexpressing Hep3B cells (Fig. 7e,f).

Fig. 7.

Overexpression of GSDMC significantly promotes the proliferation and invasion of HCC cell. (a) The full length GSDMC, GSDMC-N, and caspase-8 expression levels in HCC and adjacent hepatic tissues. n = 12. (b) Validation of GSDMC-overexpressing Hep3B cell lines, n = 3. (c) Effects of overexpression of GSDMC on the proliferation of HCC cells, n = 3. (d) Representative plot and quantitation of colony formation assay in HCC cells overexpressing GSDMC, n = 3. (e,f) Effects of GSDMC overexpression on HCC cell motility (e) and invasion (f), n = 3. **p 0.01, ***p 0.001, unpaired two-tailed t test or one way ANOVA.

Next, we constructed a GSDMC knockdown Huh7 cell line (Fig. 8a). The results showed that the proliferative capacity was decreased (Fig. 8b,c), and cell motility and invasion ability were impaired in GSDMC-knockdown Huh7 cells (Fig. 8d,e).

Fig. 8.

Knockdown of GSDMC significantly inhibits the proliferative and invasive ability of HCC cell. (a) Validation of GSDMC-knockdown Huh7 cell lines (sh-GSDMC), n = 3. (b) Effects of knockdown of GSDMC on the proliferation of HCC cells, n = 3. (c) Representative plot and quantification of colony formation assay in HCC cells knockdown by GSDMC, n = 3. (d,e) Effects of GSDMC knockdown on HCC cell motility (d) and invasion (e), n = 3. **p 0.01, ***p 0.001, unpaired two-tailed t test or one-way ANOVA.

Taken together, these findings revealed a novel tumorigenic function of GSDMC in HCC.

Since GSEA analysis revealed an enrichment of immune cell-related pathways, we hypothesized that infiltration of immune cells plays a critical role in the tumorigenic effect of GSDMC in HCC. GSDMC expression showed a positive correlation with immune cells and stromal cells (Fig. 9a). Furthermore, we investigated the association between GSDMC expression and 24 types of tumor-infiltrating immune cells by using the TIMER database. It is demonstrated that GSDMC was strongly correlated with macrophages and Th2 cells (Cor $>$0.3, p 0.05) (Fig. 9b).

Fig. 9.

The relationship between GSDMC expression and tumor microenvironment (TME) in HCC. (a) GSDMC was positively correlated with both immune cells and stromal cells. (b) High GSDMC expression had the strongest correlation with macrophages and Th2 cells (cor $>$ 0.3, p 0.05). (c) GSDMC was correlated with various immune checkpoints in the TISIDB database. (d) GSDMC was mainly expressed in keratinocytes and epithelial cells. (e) GSDMC did not express in immune cells in HCC single-cell sequencing datasets (gSE140228 and GSE98638).

Next, we explored the correlation between GSDMC and immunosuppressive checkpoint. As shown in the heatmap, GSDMC is correlated with various immune inhibitor proteins such as CD274 (PD-L1), IL-10, and IDO1 in the TISIDB database (Fig. 9c). We further examined the expression pattern of GSDMC. It is mainly expressed in keratinocytes and epithelial cells (Fig. 9d), according to the HPA dataset. GSDMC was not expressed in immune cells in HCC single-cell sequencing datasets (GSE140228 and GSE98638, Fig. 9e) using TISCH database.