Analysis of multiple programmed cell death-related prognostic genes and functional validations of necroptosis-associated genes in oesophageal squamous cell carcinoma

Abstract

Background: Oesophageal squamous cell carcinoma (ESCC) is a lethal malignancy. Immune checkpoint inhibitors (ICIs) showed great clinical benefits for patients with ESCC. We aimed to construct a model predicting prognosis and response to ICIs by integrating diverse programmed cell death (PCD) forms.

Methods: Genes related to 14 PCDs were collected to generate multi-gene signatures, including apoptosis, necroptosis, pyroptosis, ferroptosis, and cuproptosis. Bulk and single-cell RNA transcriptome datasets were used to develop and validate the model. We assessed the functions of two necroptosis-related genes in ESCC cells by Western blot, co-immunoprecipitation (Co-IP), LDH release assay, CCK-8, and migration assay, followed by immunohistochemistry (IHC) staining on samples of patients with ESCC (n = 67).

Findings: We built and validated a 16-gene prognostic combined cell death index (CCDI) by combining immunogenic cell death (ICD) and necroptosis signatures. The CCDI could also predict response to ICIs in cancer, as shown by Tumour Immune Dysfunction and Exclusion (TIDE) analysis, confirmed in four independent ICI clinical trials. Trajectory analysis revealed that HOOK1 and CUL4A might affect ESCC cell fate. We found that HOOK1 induced necroptosis and inhibited the proliferation and migration of ESCC cells, while CUL4A exhibited the opposite effects. Co-IP assay confirmed that HOOK1 and CUL4A promoted and reduced necrosome formation in ESCC cells. Data from patients with ESCC further supported that HOOK1 and CUL4A might be a tumour suppressor and oncogene, respectively.

Interpretation: We constructed a CCDI model with potential in predicting prognosis and response to ICIs in cancer. HOOK1 and CUL4A in the CCDI model are crucial prognostic biomarkers in ESCC.

Funding: The Natural Science Foundation of China [82172786], The National Cancer Center Climbing Fund of China [NCC201908B06], The Natural Science Foundation of Heilongjiang Province [LH2021H077].

Keywords: ESCC; Immune checkpoint inhibitor; Programmed cell death; Single-cell RNA-Seq.

Fig. 1

Establishment of 14 prognostic models based on different cell death forms. a. The flowchart of the study. We constructed 14 prognostic models based on genes associated with various cell death forms. b–o. Expression (Normal vs Tumour tissues) (n = 358) and results of multivariable Cox regression analysis for genes included in each prognostic model are shown for different types of cell death, including intrinsic apoptosis (b), extrinsic apoptosis (c), entotic cell death (d), necroptosis (e), cuproptosis (f), oxeiptosis (g), ferroptosis (h), alkaliptosis (i), pyroptosis (j), parthanatos (k), autophagy-dependent cell death (l), lysosome-dependent cell death (m), netotic cell death (n), and immunogenic cell death (ICD) (o) (n = 179). Hazard ratio (x axis) is presented on the logarithmic scale in the forest plots. P values for differences in gene expression between tumours and normal tissues were determined by the student t-test. CCDI, combined cell death index; ESCC, oesophageal squamous cell carcinoma; TIDE, Tumor Immune Dysfunction and Exclusion; WB, Western blot; LDH, Co-IP, co-immunoprecipitation; CCK-8, cell counting kit-8; lactate dehydrogenase; IHC, Immunohistochemistry.

Fig. 2

Evaluation of the 14 prognostic models of cell death forms in the training set (GSE53625, n = 179). a. The area under the time-dependent receiver operating characteristic (ROC) curve for prognosis-predicting models of 14 cell death forms at indicated time points. b. ROC plots of necroptosis, ICD, and combined cell death index (CCDI) models. c. Time-dependent ROC curves for the CCDI. d. The Kaplan–Meier curves for high- and low-risk patient subgroups are presented for the three models (The dotted curves depict the 95% confidence interval for the survival curve, n = 179, log-rank test, P < 0.0001). e. Uniform Manifold Approximation and Projection (UMAP) analysis and visualisation estimate the discriminating ability of necroptosis, ICD, and CCDI models for the low- and high-risk groups. f. The precision–recall curves evaluate the prognostic accuracy of indicated models (Recall is the sensitivity while precision means positive predictive value). g. Correlation between necroptosis and ICD risk scores. CCDI increases while necroptosis and ICD risk scores increase. h. The network of genes and the labelled prognostic models. Line weights measure the correlation between individual gene expression levels and the risk score of corresponding risk models, as well as the correlation among the three models concerning the risk score. Pearson correlation analysis was used to assess the correlations (n = 179). AUC, area under the curve; CI, confidence interval; OS, overall survival.

Fig. 3

Comparison of the Immunotherapy efficacy between high and low-risk groups defined by the CCDI. a. Heat map of correlation between the CCDI and immune checkpoints (Blue square represents negative correlation, while red square represents positive correlation) (GSE53625, n = 179). b. Tumor Immune Dysfunction and Exclusion (TIDE) analysis indicated the distribution of responders and non-responders to immune checkpoint inhibitors (ICIs) in the training set (GSE53625, n = 179). c. Correlation of TIDE scores of patients with CCDI (GSE53625, n = 179). d. Percentages of ICI responders in the low-risk group compared to the high-risk group (GSE53625, n = 179). e. Kaplan–Meier OS curves for the CCDI-defined high- and low-risk groups in the IMvigor210 (externally validated, log-rank test, P < 0.0001) dataset (n = 220). f, h, j, l. Comparison of the response to anti-PD-L1 treatment in the high-risk and low-risk groups in the IMvigor210 study (Bladder urothelial carcinoma, n = 137) (f), GSE78220 (Melanomas, n = 28) (h), GSE67501 (Human renal cell carcinoma, n = 11) (j) and GSE165252 (resectable oesophageal adenocarcinoma, n = 40) (l) cohorts. g, i, k, m. The differences in CCDI between the high- and low-risk groups of the IMvigor210 (g), GSE78220 (i), GSE67501 (k) and GSE165252 (m) cohorts. Pearson correlation analysis was performed to assess the correlations. Differences in immunotherapy response between high and low-risk groups were compared using the chi-square test. Differences in CCDI between the two groups were analysed using a student t-test, and one-way ANOVA was used for the comparisons among three or more groups.

Fig. 4

Clustering and annotation of single-cell RNA sequencing data in GSE196756 (n = 6). a. The t-distributed stochastic neighbour embedding (t-SNE) plot of patients, tissue types, and the eight identified main cell types in ESCC tumours and normal oesophageal tissues. b. Bubble plot showing 36 marker genes expressed in the eight-cell types. The size of the bubbles represents the proportion of cells expressing marker genes, and the spectrum of colours indicates the average expression level of marker genes (log1p transformation). c. Histograms display the percentages of different cells across normal and tumour tissues of the oesophagus. d-e. The number of interactions (d) and interaction weights/strength (e) of cell and cells in the tumour microenvironment of ESCC. The colour and width of the lines represent the number of interacting pairs between cell types. f. The bubble plots indicate the key outgoing and incoming signalling patterns of the cell, respectively. The size of the dots is proportional to the calculated contribution fraction in cellular communication, with higher contribution fractions representing signalling pathways that are more abundant in cellular interactions.

Fig. 5

Pseudo-time analysis of each cell type and dynamic changes of CCDI- gene expression in trajectories in the training set (GSE196756, n = 6). a–g. Cell subpopulations (upper panels), single cell motility trajectories (middle panels), and dynamic genes during the cell differentiation (lower panels) were visualised for dendritic cells (DC) (a), CD4+ T cells (b), B cells (c), monocytes (d), fibroblasts (e), endothelial cells (f), and ESCC cells (g). UMAP plots show the main subclusters of cells. The “Monocle2” algorithm was used to perform trajectory plots for cells and determine dynamic genes. Each dot indicates a single cell coloured by its cluster, and the solid black line shows the LOESS fit in the lower panels.

Fig. 6

Validation of the implications of CCDI-gene in cell state transition with GSE188900 (n = 7). a–g. T cells (a), B cells (b), NK cells (c), Monocytes (d), fibroblast (e), endothelial cells (f), and ESCC cells (g) were divided into subgroups (upper panels), shown in single cell motility trajectories (middle panels). The dynamic genes during the cell differentiation (lower panels) were also determined for these cell populations. Each dot indicates a single cell coloured by its cluster. The solid black line represents the LOESS fit in the lower panels.

Fig. 7

Impacts of HOOK1 and CUL4A on necroptosis, proliferation, and migration of ESCC cells. a. Differential expression of HOOK1 and CUL4A in normal oesophageal epithelial cells and ESCC cells by Western blot analysis. b. The overexpression efficiency of HOOK1 and CUL4A was examined using qRT-PCR (upper panel) and Western blot analysis (lower panel). c. KYSE410 cells were treated with increasing concentration of cisplatin for 24 h and then were subjected to Western blot analysis for necroptosis biomarkers. d-f. ESCC cells (KYSE150) with HOOK1 overexpression were treated with the indicated concentrations of different necroptosis inhibitors for 24 h, including necrostatin-1 (Nec-1), necrosulfonamide (NSA), and RIPA-56. Western blot analysis was conducted to examine necroptosis biomarkers. g, Co-IP showed that HOOK1 facilitated the precipitation of RIPK3 and MLKL with RIPK1. h, Co-IP indicated that HOOK1 reduced the RIPK3 and MLKL interacting with RIPK1. i–l. The LDH release assay was used to quantify the effect of Nec-1 (15 μM) (i), NSA (10 μM) (j), and RIPA-56 (1.5 μM) (k) on HOOK1 overexpression-mediated cell damage and impacts of CUL4A overexpression on cisplatin (15 μM)-caused necroptosis (l). m-n. The effects of overexpression of HOOK1 (m) and CUL4A (n) on the proliferation of ESCC cells by CCK-8 assay. o-r. The effects of overexpression of HOOK1 (o, q) and CUL4A (p, r) on cell migration ability by transwell assay (scale bar = 5 μm) and wound healing assay (scale bar = 20 μm). The data are presented as the means ± standard deviations (SD) of three independent experiments. The P value was determined using the student t-test.

Fig. 8

Validation of expression and prognostic potential of HOOK1 and CUL4A in Clinical specimens. a. Western blot results showed the differential expression of HOOK1 and CUL4A in tumours and matched peri-tumoral tissues from 12 patients. b-c. Representative images of immunohistochemistry (IHC) staining showing high and low expression levels of HOOK1 (b) and CUL4A (c) in peri-tumoral and tumour tissues of indicated patients with ESCC (scale bars = 10 μm and 20 μm for 100x and 50x, respectively; n = 67). d-g. Relationship between the expression level of HOOK1 and tissue type (d), tumour differentiation degree (e), clinical stage (f), and overall survival (g) in patients with ESCC. h-k. Relationship between the expression level of CUL4A and tissue type (h), degree of tumour differentiation (i), clinical stage (j), and overall survival (k) in patients with ESCC. l. Kaplan–Meier curves for patients with ESCC with different expression patterns of HOOK and CUL4A. m. Results of univariate and multivariable Cox analyses of characteristic clinical factors, HOOK1 and CUL4A in ESCC. The x-axis shows the HR (logarithmic scale) in the forest plots. Differences in gene expression between two groups were analysed using the student t-test; one-way ANOVA was adopted, while three or more groups were compared. The Kaplan–Meier curves were analysed using the log-rank test. HR, Hazard ratio.