Abstract

Objective

   Hepatocellular carcinoma (HCC) ranks as the third leading cause of
   cancer-related deaths, constituting 75%–85 % of all primary liver
   cancers. The objective of this study was to develop a
   necroptosis-related gene signature using single-cell and bulk RNA
   sequencing to predict HCC patient prognoses.

Methods

   A total of 25 key necroptosis regulators were identified from previous
   literature. We evaluated the necroptosis scores of different cell types
   using single-cell sequencing data from HCC and analyzed 168
   necroptosis-related genes. The Cancer Genome Atlas Liver Hepatocellular
   Carcinoma (TCGA-LIHC) dataset served as the training set for
   establishing a novel necroptosis-related gene risk signature, employing
   univariate and multivariate Cox regression analyses. Additionally, the
   study examined the model's relevance in immunity and immunotherapy, and
   predicted chemosensitivity in HCC patients based on the gene signature.
   The key genes were validated by the biological experiments.

Results

   Compared to other cell types, hepatoma cells displayed the lowest
   necroptosis scores. A new six-gene necroptosis-related signature
   (S100A11, MAGEC2, MAGEA6, CTP2C9, SOX4, AKR1B10) was developed using
   the TCGA database and validated in the ICGC database. Patients in the
   high-risk category had poorer prognoses, with the risk score serving as
   an independent prognostic indicator beyond other clinical factors.
   These high-risk patients also exhibited greater immune infiltration but
   demonstrated a weaker anti-tumor response due to elevated expression of
   immune checkpoints. Pathways involving hypoxia, glycolysis, and P53, as
   well as the frequency of P53 somatic mutations, were notably heightened
   in the high-risk group. Additionally, the six genes in the model showed
   significantly different mRNA expression in hepatoma cell lines compared
   to normal hepatocytes, with the role of MAGEA6 in liver cancer being
   elucidated through critical experiments.

Conclusions

   This study successfully developed a six-gene necroptosis-related
   signature to predict prognoses in HCC patients. It further explored the
   roles of necroptosis in hepatoma cells and the tumor microenvironment.

   Keywords: Hepatocellular carcinoma, Necroptosis, Prognostic signature,
   Single-cell

1. Instruction

   In terms of cancer incidence and mortality rates, hepatocellular
   carcinoma (HCC) ranks first, accounting for 75%–85 % of all primary
   liver cancer patients [[35]1]. In 2018, it was estimated that there
   were 9.3 new cases of liver cancer worldwide per 100,000 people, with
   an 8.5 death rate [[36]2]. Currently, researchers are still
   investigating the mechanisms of tumorigenesis and progression.
   Treatment options available for advanced-stage patients are very
   limited [[37]3]. According to the mainstream hypothesis, tumorigenic
   factors cause genetic mutations in somatic cells, which contributes to
   the derangement of normal genes and disruption of gene expression,
   resulting in the formation of tumor cells with morphology, metabolism,
   and function different from normal cells [[38]4]. The approach to
   treating cancers has now progressed to precision medicine, which
   incorporates the patient's treatment plan with the cancer's molecular
   characteristics, such as tumor mutational burden (TMB), microsatellite
   instability (MSI), and programmed cell death protein 1 (PD1) expression
   level. Patients with advanced HCC can only hope for targeted therapy
   and immunotherapy because curative treatments such as surgical
   resection, or liver transplantation are not options. However, a limited
   proportion of HCC patients can undergo clinical remission from targeted
   therapy and immunotherapy, and many acquire drug resistance
   [[39]5,[40]6]. Therefore, precise diagnosis and individualized
   treatment strategies are necessary to achieve precision medicine for
   HCC patients. Research on novel biomarkers is urgently needed to help
   in early clinical diagnosis, risk assessment, and therapy effectiveness
   evaluation.

   One of the most crucial characteristics of cancer is the evasion of
   cell death [[41]7]. To address this problem, researchers believe that
   novel therapy that circumvents the apoptotic pathway will be effective.
   Emerging research has revealed connections between anti-tumor immunity
   and non-apoptotic regulated cell death. Particularly, anti-cancer
   immune responses are synergistically induced by autophagy, necroptosis,
   pyroptosis, paraptosis, and ferroptosis, while simultaneously possibly
   being inhibited by these processes [[42]8]. Necroptosis is a
   caspase-independent necrotic cell death that is primarily mediated by
   receptor-interacting protein [RIP] kinase 1 (RIPK1), RIPK3, and mixed
   lineage kinase domain-like pseudokinase (MLKL) [[43]9]. Necroptosis, a
   distinct form of programmed cell death, can overcome apoptosis
   resistance and activate and boost anti-tumor immunity. Necroptosis in
   tumor has a complex biological mechanism. Crucial necroptosis
   regulators frequently exhibit downregulated expression in tumor cells
   [[44]10,[45]11], revealing that tumor cells could survive by avoiding
   necroptosis. However, expression levels of these essential mediators
   are increased in some forms types of cancer [[46]12]. For instance,
   more immunostimulatory cytokines were released in the tumor
   microenvironment (TME) due to RIPK3 activation-dependent derepression
   of TRIM28 in cancer cells, which helps promote effective anti-tumor
   immunity [[47]13]. However, RIPK1 was discovered to be upregulated in
   glioblastoma and was associated with a bad outcome via inhibiting p53
   activation [[48]14]. As a result, necroptosis may trigger a robust
   adaptive immune response that prevents tumor growth; The induced
   inflammatory response, however, can encourage carcinogenesis and cancer
   spread, so necroptosis might produce an immunosuppressive TME [[49]15].
   Immunotherapy is also greatly influenced by the altered immune
   microenvironment caused by necroptosis [[50]16]. Therefore, it is
   crucial to investigate how necroptosis affects HCC.

   Even though numerous studies have shown the significance of necroptosis
   in HCC, the precise mechanism is still unknown. According to the latest
   research, tumor-infiltrating lymphocytes (TILs), more specifically
   CD8^+ T cells, accumulated in tumors with increased expression of genes
   related to necroptosis. Furthermore, the expression value of RIPK1,
   RIPK3, and MLKL-p was related to improved prognosis [[51]17]. However,
   increased necroptosis in the liver is related to aging, which promotes
   chronic inflammation of the liver, leading to hepatic fibrosis and
   disease, including liver cancer [[52]18]. Additionally, Heparanase
   could promote HCC metastasis by inducing necroptosis in microvascular
   endothelial cells [[53]19]. Prognostic models based on ribonucleic acid
   (RNA) sequencing have made significant advances as high-throughput
   sequencing has expanded. The prognostic model built by
   necroptosis-related genes may potentially result in new treatment
   options in addition to improving prognosis prediction of HCC patients
   [[54]20,[55]21]. In this work, we combined bulk RNA and single-cell
   sequencing to analyze necroptosis genes and found that hepatoma cells
   had the lowest necroptosis score compared to other cell types,
   suggesting that hepatoma cells can escape necroptosis. We further
   identified genes associated with necroptosis in HCC. Researchers
   developed a prognostic signature based on necroptosis-related genes
   that could precisely forecast the outcome of HCC patients and recommend
   immunotherapeutic treatment. In addition, we investigated the
   significance of prognostic models in immune function and immunotherapy.
   The association between prognostic models and several pathways,
   including the P53 and apoptosis, was also investigated.

2. Methods and materials

2.1. Data acquisition

   We obtained bulk RNA sequencing data and associated clinical
   information from three databases: the Cancer Genome Atlas (TCGA), the
   International Cancer Genome Consortium (ICGC, LIRI-JP), and the
   [56]GEO14520 database. These databases provided 424 (comprising 374
   tumor and 50 adjacent normal tissues), 231 tumor samples, and 220 tumor
   samples, respectively. Patients lacking survival data were excluded.
   Concurrently, we accessed single-cell RNA sequencing data from
   [57]GSE149614 [[58]22], which encompassed 10 tumor and 7 normal
   tissues. Given our research's emphasis on the heterogeneity of the
   Tumor Microenvironment (TME) in tumor and peri-tumor tissues, we
   selected seven pairs of tumor and normal samples for analysis.

2.2. Single cell analysis

   All single cells were enrolled in the quality control process performed
   by the Seurat R package [[59]23]. Cells were included in the analysis
   if they met the following criteria: (a) RNA count range between 200 and
   500, and (b) mitochondrial gene expression below 5 %. This resulted in
   the analysis of 55,341 cells, comprising 26,660 cells derived from
   tumors and 28,681 from normal tissues. Data normalization was performed
   using the “LogNormalize” algorithm with a scale factor of 10,000. The
   impact of Unique Molecular Identifiers (UMIs) and mitochondrial content
   percentages were mitigated using Seurat's ScaleData function. The
   Harmony R package was employed to eliminate batch effects. For cell
   clustering and t-Distributed Stochastic Neighbor Embedding (t-SNE)
   visualization, the top 30 principal components and the 2000 most
   variable genes were chosen [[60]24]. Cell cluster-specific marker genes
   were identified using the FindAllMarkers function, and cell types were
   annotated based on canonical marker genes documented in prior studies
   (refer to [61]Table S1).

2.3. The Acquisition of necroptosis-related genes by scRNA-seq

   [62]Table S2 represents the twenty-five key necroptosis regulators that
   were extracted from earlier reviews [[63]15,[64]25]. Utilizing
   twenty-five essential necroptosis regulators, we computed the
   necroptosis-related scores for each cell using the AddModuleScore
   function in the Seurat R package. Notably, hepatoma cells exhibited the
   lowest necroptosis scores. We subsequently isolated these cells for
   further examination and observed that, apart from those in cluster 1,
   all hepatoma cells consistently showed low necroptosis scores. This led
   us to identify differentially expressed genes (DEGs) between hepatoma
   cells from cluster 1 and those from other clusters, using the
   FindAllMarkers function. These DEGs are hypothesized to be associated
   with necroptosis in HCC.

2.4. Construction of a prognostic signature associated with necroptosis

   Our analysis commenced by intersecting 168 differentially expressed
   genes (DEGs) with genes sequenced in the ICGC-LIHC and [65]GEO14520
   datasets. We then examined the expression levels of 125 overlapping
   DEGs between tumor and normal tissues in the TCGA cohort, identifying
   38 DEGs that exhibited significant differential expression (P < 0.05,
   log fold change >1). The TCGA LIHC cohort served as a training set,
   integrating the gene expression matrix of these 38 DEGs with survival
   status and duration. A univariate Cox regression analysis, conducted
   using the survival R package, determined the predictive relevance of
   these 38 DEGs. Briefly, univariable Cox regression analyses were
   conducted to screen DEGs with prognostic significance for inclusion in
   the multivariable model, using a liberal p-value threshold of <0.10.
   Then, a multivariable Cox proportional hazards model was built using a
   forward stepwise approach. At each step, the model automatically
   assessed the statistical significance (p < 0.05) of each variable and
   decided whether to retain or remove it from the model.

   The regression coefficients, hazard ratios, and 95 % confidence
   intervals were reported for each predictor variable in the final Cox
   model.

2.5. Evaluation and Validation of the risk model

   The risk score for each patient was calculated in both the training
   (TCGA) and testing (ICGC) sets, with patients divided into high- and
   low-risk groups based on the median risk score. Survival rates of these
   groups were compared using Kaplan-Meier curves. Additionally, Receiver
   Operating Characteristic (ROC) curve analysis was conducted to assess
   the predictive accuracy of the model in both sets. Further, univariate
   and multivariate Cox regression analyses were performed to evaluate
   whether the risk score independently predicted the prognosis of
   hepatocellular carcinoma (HCC) patients, considering factors such as
   tumor size, stage, age, grade, gender, and risk score. We also explored
   the association between the risk score and clinical characteristics,
   including tumor size, stage, gender, and grade.

2.6. Function enrichment analysis

   he immune function scores of each patient were calculated using the
   GSVA R package. We then assessed the average expression levels of
   immune checkpoints in both high- and low-risk groups, with a
   significance threshold set at P < 0.05. Gene Set Enrichment Analysis
   (GSEA) was employed to identify enriched gene sets within these groups.
   We sourced the 50 hallmark gene sets from MsigDB
   ([66]https://www.gsea-msigdb.org/gsea/index.jsp). Additionally, Gene
   Ontology (GO) pathway enrichment analysis and Kyoto Encyclopedia of
   Genes and Genomes (KEGG) analysis for the differentially expressed
   genes (DEGs) were conducted using the clusterProfiler R package
   [[67]26]. The functions of DEGs were further explored using Metascape
   ([68]http://metascape.org), an advanced online tool for gene function
   annotation.

2.7. Evaluation of Tumor Mutation Burden and Drug Sensitivity

   Tumor Mutation Burden (TMB) was computed utilizing somatic mutation
   data from the TCGA datasets, with comparative analyses performed
   between the low- and high-risk subsets. The "maftools" R package
   facilitated the creation of waterfall plots, illustrating the variance
   in gene mutation characteristics between these subsets. We also
   investigated the combined impact of TMB and risk score on the survival
   of hepatocellular carcinoma (HCC) patients. Drug sensitivities were
   assessed using the pRRophetic R package, which calculates the
   half-maximal inhibitory concentration (IC50) of various drugs, thereby
   determining the differential drug sensitivities in high-versus low-risk
   groups.

2.8. Cell culture

   Human hepatocellular carcinoma (HCC) cell lines, including SMMC-7721,
   HepG2, Huh-7, and normal human hepatocytes L-02, were procured from the
   Cell Bank of the Chinese Academy of Sciences, Shanghai, China. The
   hepatoma cell lines were cultured in Dulbecco's Modified Eagle Medium
   (DMEM) supplemented with 10 % fetal bovine serum (FBS) and 1 %
   penicillin-streptomycin, in a humidified incubator maintained at 37 °C
   with 5 % CO2. Conversely, the L-02 cell lines were grown in RPMI Medium
   1640.

2.9. qRT-PCR

   Cells were harvested for RNA extraction once they occupied
   approximately 80%–90 % of a 10 cm culture dish. Total RNA was extracted
   using TRIzol reagents according to the standard protocol. Subsequently,
   cDNA was synthesized using reverse transcription kits, and the mRNA
   expression levels were quantified via qPCR assays. GAPDH
   (glyceraldehyde-3-phosphate dehydrogenase) served as the internal
   reference for mRNA, and expression levels were calculated using the
   2-ΔΔCt method. The primers for SOX4, AKR1B10, MAGEC2, CYP2C9, S100A11,
   and MAGEA6, obtained from Sangon Biotech, are detailed in
   [69]Supplemental Table S2. Additionally, protein expression levels of
   the six key necroptosis-related genes were examined using the Human
   Protein Atlas (HPA) online platform.

2.10. Western blot analysis

   In 10 % polyacrylamide gels, whole cell proteins were electrophoresed
   under decreasing conditions. For 1 h, the electrophoresis was performed
   in MOPS buffer at 180 V. The membrane was blocked in 5 % nonfat dry
   milk and then incubated with an antibody (Abcam, Cambridge, UK) at a
   dilution of 1:1000 for a whole night at 4 °C. Using an ECL western
   blotting kit (Amersham Biosciences, Little Chalfont, UK) and following
   the manufacturer's instructions, all band intensities were assessed.
   Image-J software was then used to examine the data.

2.11. Colony assay

   One million five hundred cells per six-well plate were planted with
   cells. After the cells were combined, they were cultivated in culture
   media with 10 % FBS and L-15 for ten days. The following standards were
   taken into account when assessing the outcomes: Colonies were defined
   as groups of at least thirty cells.

2.12. Cell viability assays

   2000 cells per well of 96-well plates were used for cell seeding. The
   Cell Counting Kit-8 assay (CCK-8, Dojindo, Japan) was utilized to
   detect cell viability in compliance with the manufacturer's protocol.
   Using an automatic microplate reader (BioTek, Winooski, VT, USA), the
   absorbance at 450 nm was determined. For seven days in a row,
   measurements were taken every 24 h.

2.13. Cell migration assays

   In the migration assay, 1.5 × 105 cells were introduced to 200 μl
   serum-free DMEM in the upper chamber of a Transwell plate (Corning, NY,
   USA), and 700 μl DMEM with 20 % serum was added to the lower chamber.
   Following a 24-h incubation period at 37 °C, the Transwell chamber was
   taken out, cleaned once using PBS, and then fxed for 30 min using a
   methanol:acetone ratio of 1:1. It was then stained for 30 min using
   0.5 % crystal violet. After carefully wiping away the cells on the
   upper side of the filter, the chamber was washed with PBS. Using
   ImageJ, the relative cell density was calculated.

2.14. Statistical analysis

   Univariate and multivariate analyses determined the prognostic marker
   genes. Analyses between the two groups were performed using the Wilcox
   test, and a one-way ANOVA test was used to compare three or more
   groups. All statistical analyses were performed using R version 4.1.3
   (R Core Team) and GraphPad Prism version 8.3 (GraphPad Software, San
   Diego, CA, USA).The cutoff for significance was chosen at P < 0.05.

3. Results

3.1. Identified the necroptosis-related genes by using scRNA

   A total of 55341 single cells fulfilled the quality control (QC)
   described in the materials and methods, including 26660 cells derived
   from cancer tissue and 28681 cells derived from peri-tumor samples.
   These cells were organized into 43 clusters and annotated into 13 main
   cell types ([70]Fig. 1A).Then we used 25 Necroptosis genes
   ([71]Supplemental Table S1) and the function AddModuleScore by the
   Seurat package to annotate score of the necroptosis for the each single
   cell of HCC. [72]Fig. 1B showed the score of the necroptosis in the
   distribution of scRNA data. We also listed the expression of each
   necroptosis gene in each cell type([73]Fig. 1C). The necroptosis score
   in the epithelial and endothelial cells were lower than other cell
   types ([74]Fig. 1D). When apoptosis is suppressed, a regulated cell
   death called necroptosis causes inflammation and acts as a backup
   defense mechanism [[75]27]. Hence, we proposed that these hepatoma
   cells can promote tumorigenesis and progression by hindering
   necroptosis [[76]28]. Among these cell lineages, the clusters
   1,5,6,13,15,19,37,42 which annotated as hepatoma cells were remained to
   further analysis. Necroptosis scores of cluster 1 hepatoma cells were
   low ([77]Fig. 1E). To confirm the genes associated to necroptosis in
   HCC, the DEGs between hepatoma cells from cluster 1 and other clusters
   were obtained using the FindAllMarkers function. We obtained 168
   Necroptosis-related DEGs by scRNA-seq analysis in HCC ([78]Fig. 1F and
   [79]Supplemental Table S3). Then, the intersection of 168 DGEs with
   sequenced genes among multi datasets, 125 genes were all expressed in
   the ICGC-LIHC and ICGC-LIHC datasets([80]Fig. 1G). Among them, 38 DEGs
   were discovered to have significantly different expression values
   between tumor and normal samples in the TCGA LIHC cohort ([81]Fig. 1H).
   Using SCP function and GO function analysis, we found that these DEGs
   were associated with several function and KEGG pathways involved in
   necroptosis, such as fibrinolysis, response to negative regulation of
   wound healing, acute−phase response, negative regulation of blood
   coagulation, Complement and coagulation cascades, cholesterol
   metabolism and linoleic acid metabolism. etc([82]Fig. 1I).

Fig. 1.

   [83]Fig. 1
   [84]Open in a new tab

   Analysis of necroptosis-related genes in HCC using single-cell and bulk
   RNA sequencing data.

   (A) Distribution of 55,023 cells across 13 cell types in the HCC
   single-cell RNA (scRNA) dataset. (B) Assessment of necroptosis scores
   in single-cell RNA data using the AddModuleScore method. (C) Boxplot
   illustrating variations in necroptosis scores across 13 cell types. (D)
   Expression profiles of necroptosis genes in the scRNA dataset. (E)
   Necroptosis score expression in various cancer cell clusters from the
   scRNA dataset. (F) Differential gene expression between cluster 1 and
   other clusters in cancer cells. (G) Venn diagram (VeenPlot) showing the
   overlap of 168 DEGs among scRNA DEGs, ICGC-LIHC, and [85]GEO14520
   datasets. (H) Expression of 125 DEGs between normal and tumor samples
   in the TCGA-LIHC cohort. (I) Gene Ontology Biological Process (GO-BP),
   Kyoto Encyclopedia of Genes and Genomes (KEGG), and WikiPathways
   analysis for the 125 DEGs, comparing high and low necroptosis-related
   scores.

3.2. Establishment and Validation of necroptosis-related modeling signature

   We identified 12 DEGs connected to prognosis in HCC using univariate
   Cox regression analysis, with 11 unfavored and 1 favored gene
   ([86]Supplemental Table S4). These prognosis-related DEGs were
   incorporated in the subsequent multivariate Cox regression algorithm to
   improve the model's independent predicting capacity. Finally, six DEGs
   were integrated as necroptosis-related modeling genes (NRMGs) to
   generate the necroptosis-Related modeling Score (NRMG-Score). The
   following formula was then utilized to determine the risk-score:
   Riskscore = (0.0007∗ S100A11) + (0.0152∗ SOX4) + (0.0005∗
   AKR1B10) + (0.0089∗ MAGEC2) + (−0.0019∗ CYP2C9) + (0.01829∗ MAGEA6).

   In the TCGA cohort (Training dataset), patients in the high-risk subset
   possessed a poor outcome, which was further confirmed by ICGC dataset
   ([87]Fig. 2A). Additionally, we also observed high-risk patients in
   [88]GSE14520 cohort possessed a negative prognostic tendency
   (p = 0.08). Every patient in the training and test subsets had their
   risk-score determined, and the median risk-score value was applied to
   categorize patients into high- and low-risk subsets ([89]Fig. 2B and
   C).Meanwhile, the death rate for the high-risk subset was also
   significantly greater than for the low-risk subgroup ([90]Fig. 2D).
   Furthermore, our ROC analysis suggested the risk model performed with
   good predictive accuracy, with area under curve (AUC) values of 0.706
   in the one-year ROC curve, 0.668 in the two-year ROC curve, and 0.642
   in the three-year ROC curve in TCGA sets. AUC values of one, two, and
   three-year ROC curves in ICGC sets were 0.724, 0.666, and 0.672,
   respectively ([91]Fig. 2E).

Fig. 2.

   [92]Fig. 2
   [93]Open in a new tab

   Evaluation of survival, risk scores, and predictive accuracy in the
   necroptosis-related model for HCC patients.

   (A) Kaplan-Meier (K-M) curves for overall survival (OS) comparing
   high-risk and low-risk groups across TCGA, ICGC, and [94]GSE14520
   datasets. (B) Expression patterns of the six genes in the model across
   high- and low-risk groups. (C) Classification of HCC patients into
   high- and low-risk categories based on median risk score. (D)
   Distribution of survival status and risk scores of HCC patients in
   TCGA, ICGC, and [95]GSE14520 datasets. (E) Receiver Operating
   Characteristic (ROC) curves of the Necroptosis-Related Prognostic Model
   for 1-, 3-, and 5-year OS predictions.

3.3. NRMGs score was an independent clinical risk factor

   PCA plot displayed high-risk patients could be distinguished from
   low-risk patients using our NRMG-SCORE predict model ([96]Fig. 3A).
   Furthermore, univariate and multivariate Cox regression algorithms
   revealed the risk-score was a prognostic indicator independent of
   certain other clinical characteristics in both TCGA and ICGC sets
   ([97]Fig. 3B and C). The risk-score is correlated with tumor size,
   stage, and tumor grade ([98]Fig. 3D), and patients with high tumor
   grade and stage owned higher risk-score (P < 0.05).

Fig. 3.

   [99]Fig. 3
   [100]Open in a new tab

   Association of necroptosis-related model gene (NRMG) scores with
   Clinicopathological Features in HCC.

   (A) Principal Component Analysis (PCA) illustrating distinct
   distribution of patients in different risk groups within TCGA and ICGC
   datasets. (B) and (C) Univariate (upper) and multivariate (bottom)
   analyses indicating the risk score as an independent prognostic factor
   in comparison to other clinical characteristics in both TCGA (B) and
   ICGC datasets. (D) Variation in NRMG scores across different TNM stages
   and grades in the TCGA cohort.

3.4. Immune infiltration for NRMG-SCORE in bulk sequence RNA level

   We calculated the immune cell infiltrates for high- and low-risk
   subsets to examine the potential roles of NRMG-SCORE in tumorigenesis
   and development, including SOX4, AKR1B10, MAGEC2, CTP2C9, S100A11, and
   MAGEA6. Interestingly, individuals in the high-risk subset exhibited
   higher immune cell infiltrates. Treg cells, Follicular helper T cells,
   Th2 cells, and macrophages were enriched in the high-risk subset
   ([101]Fig. 4A). However, when we calculated the degree of immune
   checkpoint expression, we found that most of these genes, including
   PDCD1(programmed cell death protein 1) and CTLA4, were significantly
   overexpressed in high-risk subset ([102]Fig. 4B). We proposed although
   patients in the high-risk subset possessed higher levels of immune
   infiltration, the high expression of immune checkpoints resulted in a
   weak anti-tumor state. These patients may gain high response rates
   during immunotherapy. Our necroptosis-related risk model might offer a
   fresh way to forecast the effectiveness of immunotherapy. The most
   important component of TME, tumor immune microenvironment (TIME), is
   essential for the tumorigenesis and treatment resistance [[103]29].

Fig. 4.

   [104]Fig. 4
   [105]Open in a new tab

   Analysis of immune infiltration, Functional enrichment, and drug
   sensitivity in Relation to necroptosis-related model gene (NRMG)
   scores.

   (A) Bar plot depicting differences in immune cell infiltration between
   high-risk and low-risk groups. (B) Comparative expression levels of
   immune checkpoints in high- and low-risk groups. (C) Gene Set
   Enrichment Analysis (GSEA) results highlighting enrichment of WNT
   beta-catenin signaling, P53 pathways, glycolysis pathways, and hypoxia
   in the high-risk group. (D) Waterfall plot showing somatic mutation
   characteristics in patients with high and low risk scores. Each column
   represents an individual patient, with the upper bar plot indicating
   Tumor Mutation Burden (TMB), and the right-hand numbers showing
   mutation frequency in each gene. (E) Correlation between risk score and
   TMB. (F) Kaplan-Meier curve for overall survival in high- and low-TMB
   groups. (G) Kaplan-Meier curve evaluating patient overall survival
   based on a combination of risk score and TMB. (H) Comparison of
   estimated IC50 values for various drugs between high-risk and low-risk
   groups.

3.5. Function enrichment analysis for NRMG-score in bulk sequence RNA level

   GSEA results revealed that WNT beta-catenin signaling, P53 pathways,
   glycolysis pathways, and hypoxia were enriched in the high-risk subset
   ([106]Fig. 4C). Inadequate blood supply or hypoxia is a typical
   microenvironmental characteristic of almost all solid tumors due to
   their rapid and uncontrolled proliferation. Hypoxic regions promote
   abnormal angiogenesis [[107]30,[108]31], adhesion proliferation, and
   inflammation [[109]32]. GO analysis displayed immune-related pathways,
   like leukocyte-mediated immunity, T cell activation, leukocyte
   migration, and cytokine-mediated signaling pathway, were significantly
   enriched in high-risk subset patients ([110]Fig. S1). Moreover, KEGG
   analysis uncovered that glycolysis/gluconeogenesis, IL-17 signaling
   pathway, ECM-receptor interaction, and cytokine-cytokine receptor
   interaction were up-regulated in the low-risk subset ([111]Fig. S2).

3.6. The mutation Landscape in Necroptosis-Related Prognostic Model

   According to the GSEA results, theP53 pathway was concentrated in the
   high-risk subset. We then estimated the somatic mutation rate in the
   two risk subsets. The P53 mutation rate was 37 % in 178 high-risk
   subset patients, which was obviously greater than the 15 % P53 mutation
   rate in the low-risk subset patients ([112]Fig. 4D). TP53 is a tumor
   suppressor gene mutated in more than half of all human cancers. Tumors
   with P53 mutations progress faster, respond poorly to anti-cancer
   therapy, and have a worse outcome [[113]33]. P53 mutations are more
   common in high-risk patients, implying that NRMGs are involved in the
   P53 pathway. The TMB values were compared between the high- and
   low-risk subsets, and it was obvious that patients in the high-risk
   subset had significantly higher TMB than those in the low-risk subset
   (p = 0.0025) ([114]Fig. 4E). Patients with cancer typically have a poor
   prognosis if their TMB is high. In our research, the HCC patients’
   outcome with a high TMB was also worse ([115]Fig. 4F). We combined TMB
   and NRMGs for survival analysis due to their synergistic effect on
   prognosis. According to our findings, patients with a high-risk score
   and a high TMB owned the worst prognosis ([116]Fig. 4G). These analyses
   suggested that the combination of NRMGs and TMB could further refine
   the prognosis prediction of patients.

3.7. Analyses of drug sensitivity and Cell Localization of six modeling genes

   We calculated the concentration causing a 50 % reduction growth (IC50)
   of drugs to determine which drugs have significantly different
   sensitivities between high- and low-risk subsets using the pRRophetic R
   package to integrate prognostic models with clinical practice. The
   high-risk subset patients possessed higher sensitivity to Sorafenib,
   Thapsigargin, Tubastatin, Ruxolitinib, Rapamycin, 5-Fluorouracil, AKT
   inhibitor VIII, Doxorubicin, Imatinib, and Lisitinib ([117]Fig. 4H).

3.8. Potential necroptosis regulated model genes promote the progress of HCC

   We used single-cell sequencing to determine the expression levels of
   the modeling genes in each cell type. As seen in [118]Fig. 5A and B,
   SOX4, AKR1B10, MAGEC2, CTP2C9, and MAGEA6 were primarily expressed in
   hepatoma cells and hepatocytes, while S100A11 was detected in hepatoma
   cells, hepatocytes, immune cells, fibroblasts, and endothelial cells.
   We also detected the protein expression values of six modeling genes
   using the HPA database. HCC tumor tissue exhibited overexpression of
   SOX4, AKR1B10, MAGEC2, and CTP2C9 relative to normal tissue ([119]Fig.
   5C). Simultaneously, we compared the mRNA expression levels of
   L-02 cell lines and three HCC cell lines. MAGEA6, AKR1B10, and S100A11
   had higher mRNA expression levels in SMMC-7721, HepG2, and Huh-7, than
   in L-02. The mRNA expression of MAGEC2 was upregulated in SMMC-7721 and
   Huh-7 and decreased in HepG2. In contrast to L-02, CTP2C9 and SOX4
   levels were lower in HCC cell lines ([120]Fig. 5D).

Fig. 5.

   [121]Fig. 5
   [122]Open in a new tab

   Elucidating the Central role of MAGEA6 in HCC based on
   necroptosis-related model gene (NRMG) score.

   (A) UMAP plot showing expression levels of model genes across different
   cell types. (B) Boxplot comparing expression of six key genes between
   tumor and normal tissues in single-cell RNA (scRNA) data. (C)
   Immunohistochemical images from the Human Protein Atlas (HPA) dataset
   illustrating expression of the six model genes in tumor and adjacent
   non-tumor tissues.(D) Quantitative real-time PCR (qRT-PCR) analysis
   revealing mRNA levels of six NRMGs in hepatoma and normal hepatocyte
   cell lines. (E) Immunofluorescence assays indicating highest
   fluorescence intensity in SMMC-7721 cells. (F) qPCR results
   demonstrating significant inhibition of MAGEA6 expression in
   SMMC-7721 cells by Si-MAGEA6-1 and Si-MAGEA6-2. (G) Elevated expression
   levels of necroptosis-related proteins (p-RIPK1, p-RIPK3, and p-MLKL)
   following MAGEA6 knockdown.

   In order to explore the expression of MAGEA6 in HCC cells, we conducted
   immunofluorescence experiments, and the results showed that the
   fluorescence intensity was the highest in SMMC-7721 cells ([123]Fig.
   5E), which proved that MAGEA6 was the highest in SMMC-7721 cells, and
   we selected SMMC-7721 cells for further experiments. To gain insights
   into the function of MAGEA6, we characterized the oncogenic phenotypes
   in SMMC-7721 cells with MAGEA6 silencing (Si-MAGEA6-1 and Si-MAGEA6-2).
   The qPCR results showed that Si-MAGEA6-1 and Si-MAGEA6-2 could
   significantly inhibit MAGEA6 expression in SMMC-7721 cells ([124]Fig.
   5F), and we chosen Si-MAGEA6-1 for cell function experiments. Next,to
   figure out the effect of si MAGEA6 on cell necroptosis, Western blot
   assay was performed. As shown in [125]Fig. 5C, the expression levels of
   cell necroptosis-related proteins (p-RIPK1, p-RIPK3 and p-MLKL) were
   up-regulated after MAGEA6 knockdown.

3.9. The biological role of MAGEA6 in SMMC-7721 cell lines

   We investigated the role of the MAGEA6 in the proliferation of HCC
   cells by colony-formation, CCK8 and immunofluorescence assay. To assess
   the influence of the MAGEA6 on cell death and apoptosis, a LDH, PI and
   flow cytometry assays was used to detect. colony-formation ([126]Fig.
   6A), CCK8 ([127]Fig. 6B) and immunofluorescence ([128]Fig. 6C) assay
   revealed that knockdown of MAGEA6 inhibited the proliferation of HCC
   cells. Additionally, knockdown of MAGEA6 significantly promoted HCC
   cell death and apoptosis ([129]Fig. 6D–F). [130]Fig. 6 showed the
   transwell assays revealed that knockdown of circCCDC134 impaired the
   migration of HCC cells.

Fig. 6.

   [131]Fig. 6
   [132]Open in a new tab

   Investigating the impact of MAGEA6 on hepatocellular carcinoma (HCC)
   Cell proliferation.

   (A) Colony formation assay results. (B) Cell Counting Kit-8 (CCK8)
   assay outcomes. (C) Findings from immunofluorescence assays. (D-F)
   Results indicating that knockdown of MAGEA6 notably enhances cell death
   and apoptosis in HCC cells. (G) Transwell assays demonstrating reduced
   cell migration following the knockdown of circCCDC134 in HCC cells.

4. Discussion

   Human cells destroy themselves to preserve biological equilibrium and
   shield the body from pathogens, a process known as regulated cell death
   (RCD). Moreover, RCD is classified as apoptotic and non-apoptotic forms
   [[133]34]. For more than three decades, clinical researchers have been
   working on therapies that promote the effective elimination of cancer
   cells via apoptosis. However, tumor cells can avoid apoptosis,
   resulting in treatment resistance and cancer relapse [[134]35].
   Researchers are focusing on alternative cancer cell death processes,
   including autophagy, necroptosis, pyroptosis, paraptosis, and
   ferroptosis, and converting apoptosis to these RCDs has the potential
   to overcome apoptosis resistance. Contradictory findings from research
   on necroptosis' role in cancer suggest that necroptosis has distinct
   impacts at different stages of cancer cell growth and metastasis
   [[135]12,[136]36,[137]37]. Tumor cells and non-tumor cells interact to
   generate a complicated TME. Researchers hope to restrict tumor growth
   by inducing necroptosis in malignant cells. However, necroptosis would
   change the tumor microenvironment to a certain degree, and the impact
   on TME is uncertain. In contrast, necroptosis may prevent tumors from
   anti-cancer immune responses by promoting an immunosuppressive
   environment [[138]13,[139]38]. Contrarily, necroptosis may cause strong
   anti-tumor immunity, and its processes and ICIs may synergistically
   improve anti-tumor activity, even in cancer types resistant to ICIs
   [[140]39]. Over the last decade, targeted therapies have paved the path
   for personalized precision medicine. Targeted therapy like Sorafenib,
   Lenvatinib, and ICIs, however, may only be effective for a small
   percentage of HCC patients. Furthermore, drug-related side effects and
   resistance continue to be a concern. Therefore, exploring the
   mechanisms of necroptosis in HCC and their impact on TIME will pave the
   way for ICI-based combination therapy. Meanwhile, exploring the
   necrosis-related genes may lead to new therapeutic targets or enhance
   the prediction of HCC prognosis. In this work, we explored genes
   related to necroptosis using single-cell sequencing to see if they
   could predict HCC prognosis. After univariate and multivariate Cox
   regression algorithms, a necroptosis-related genes risk model
   comprising six genes was established using the TCGA database and
   verified in the ICGC database. Patients in the high-risk group
   possessed a poor outcome, and the risk-score was a prognostic indicator
   that could be used independently of other clinical traits.
   Additionally, the high-risk subset patients exhibited increased levels
   of immune infiltration, and the high level of immune checkpoints
   resulted in a weak anti-tumor state. Hypoxia, glycolysis, and P53
   pathway were enriched in the high-risk subset; the rate of P53 somatic
   mutations was also considerably greater. Furthermore, expression values
   of six modeling genes differed significantly between normal hepatocyte
   and hepatoma cell lines.

   We used HCC single-cell sequencing data for the first time to assess
   the necroptosis score for each cell type. Hepatoma cells exhibited the
   lowest necroptosis score when compared with other cell types, revealing
   they may have a defense mechanism against necroptosis. Several studies
   have found that inducing necroptosis in hepatoma cell lines can inhibit
   the growth of tumor cells. Xiang et al. identified shikonin, a
   necroptosis inducer, may be more effective than apoptosis inducers in
   eliminating HCC cells with high Cx32 expression [[141]40]. According to
   recent research, excessive Sorbitol dehydrogenase in HCC cells hinders
   the growth of tumor cells by promoting necroptosis [[142]41].
   Therefore, stimulating necroptosis to kill hepatoma cells is a
   promising therapeutic approach. However, the downstream mechanisms of
   necroptosis have not been fully studied. We further explored
   necroptosis-related genes using single-cell sequencing techniques to
   find more targets associated with necroptosis. Using single-cell data,
   168 DEGs associated with necroptosis were identified. They were
   enriched in response to a toxic substance, acute inflammatory response,
   and neutrophil degranulation, all of which were involved in
   necroptosis. Based on these DGEs, we successfully constructed a
   prognostic model that precisely estimates the difference in prognosis
   between HCC patients in the high and low-risk cohorts. ROC analysis
   indicated the gene signature worked well in predictive accuracy, with
   AUC values of 0.706 in the one-year ROC curve, 0.668 in the two-year
   ROC curve, and 0.642 in the three-year ROC curve in TCGA sets. The AUC
   value of the one, two, and three-year ROC curves in ICGC sets were
   0.724, 0.666, and 0.672, respectively. Besides, our risk model was a
   prognostic indicator independent of other clinical characteristics,
   such as tumor size, stage, and grade. Therefore, our CCGs prognostic
   model may be employed to evaluate the outcome and survival time of HCC
   patients.

   Six NRMGs have been studied in HCC, including S100A11, MAGEC2, MAGEA6,
   CYP2C9, SOX4, and AKR1B10. During necroptosis, a large quantities of
   damage associated molecular patterns (DAMPs), including S100 protein
   family members, are released [[143]42]. A calcium-binding protein
   called S100 is connected to numerous pathological and disease
   processes, including inflammation and carcinogenesis. Numerous members
   of the S100 protein family appear to be crucial in the emergence of
   hepatic steatosis, non-alcoholic steatohepatitis (NASH), and HCC
   [[144]43]. High-grade HCC and poor prognosis are highly correlated with
   overexpression of S100A11, which promotes inflammation/fibrosis and
   early dysregulated steatosis [[145]44]. In intrahepatic
   cholangiocarcinoma, P38/MAPK signaling pathway activation, colony
   formation, and proliferation of RBE and HCCC9810 cells were all
   dramatically reduced when S100A11 was silenced In vitro. In vivo, tumor
   growth was also prevented [[146]45]. MAGEA6 and MAGEC2 are members of
   the MAGE family, and the MAGE gene originally received attention among
   oncologists because of its strong tumor antigenic characteristics in
   melanoma. MAGE genes are found to be dysregulated in a variety of
   tumors [[147]46]. miR-448 promotes the AMPK signaling pathway through
   the downregulation of MAGEA6, thereby inhibiting stemness maintenance
   and proliferation of HCC stem cells [[148]47]. Similarly, MAGEC2 has
   been associated with metabolic reprogramming in HCC [[149]48]. The
   expression of CYP2C9, one of the most prevalent and significant
   xenobiotic metabolizing enzymes, is downregulated in HCC. CYP2C9 is
   involved in developing hepatocarcinogenesis and regulates cell
   proliferation, apoptosis, differentiation, and homeostasis [[150]49].
   SOX4 activates CXCL12 to induce endothelial cell chemotaxis,
   angiogenesis, and tumor growth [[151]50]. Meanwhile, targeting SOX4 can
   inhibit cell proliferation and migration [[152]51]. Higher levels of
   AKR1B10 expression are linked to a worse long-term prognosis, making it
   a viable biomarker for HCC. According to multiple experimental studies,
   silencing AKR1B10 causes cell death and inhibits tumor development and
   spread [[153]52]. Meanwhile, the down-regulated expression of AKR1B10
   in NK cells promotes glycolysis while inhibiting migration, invasion,
   and metastasis of hepatoma cells [[154]53]. In this research, we
   explored the mRNA value of six NRMGs in hepatoma cells and hepatocyte
   cell lines using qrt-PCR. The six modeling genes expression
   significantly differed between normal hepatocyte and hepatoma cell
   lines, indicating they might be a potential target for treatment.

   Effective anti-tumor immunity can fight cancer cells and boost the
   prognosis of HCC patients. Cancer cells suppressed the immune cell's
   function to promote immune escape by creating a hypoxic and
   hyper-lactic environment and reprogramming immune cells [[155]54].
   Therefore, stimulating anti-tumor immunity is a promising therapeutic
   approach, and ICIs are developed based on this theory. Currently,
   researchers are investigating tumor immunity and looking for fresh
   targets to activate anti-tumor immune responses. Individual in the
   high-risk subset in our necroptosis-related prognostic signature
   exhibit greater immune cell infiltration. Meanwhile, most immune
   checkpoints were over-expressed in the high-risk subset, and immune
   checkpoints activation can suppress T lymphocyte function and hence the
   auto-immune response. Immune checkpoint expression has become a
   predictor of immunotherapy effectiveness, with higher expression
   resulting in higher response rates. Hence, the different risk-scores
   based on the necroptosis-related rick model might serve as a promising
   indicator for immunotherapy strategy. In addition, TMB was noticeably
   higher in the high-risk subset, and patients in this subset owned a
   poor prognosis. TMB is a commonly used immunotherapy response biomarker
   in clinical practice, with high TMB levels indicating a high rate of
   immunotherapy response. This is consistent with our previous inference
   that patients in the high-risk subset are probably receiving benefits
   from ICIs. Finally, compared to the low-risk subset, the high-risk
   subset exhibited a much greater rate of P53 mutations, with most of
   these mutations being missense mutations. TP53 is a tumor suppressor
   gene, suggesting that NRMGs are involved in the P53-related pathway,
   which was also confirmed by GSEA analysis. Functional analysis
   suggested that NRMGs are involved in immune-related pathways including
   leukocyte-mediated immunity, T-cell activation, leukocyte migration,
   and cytokine-mediated signaling pathways. DAMPs, as well as
   inflammatory cytokines and chemokines, can be released in large
   quantities during necroptosis, resulting in an inflammatory response
   [[156]55,[157]56]. According to Zhang et al., ICI resistance in
   melanoma could be overcome by inducing ZBP1-dependent necroptosis in
   cancer-associated fibroblasts [[158]57]. Overall, necroptosis and its
   immunobiology cannot be generalized or applied to a complicated
   pathology like cancer in a “black and white” way. To effectively
   utilize necroptosis as a therapeutic target, the specific chemical
   mechanisms and signaling systems involved must be understood.

   In sum up, a six NRMGs risk model was developed for the independent
   prediction of the prognosis of HCC patients. The riskscores could be
   employed as an indicator for immunotherapy.

Funding

   This work was supported by the National Natural Science Foundation
   project of China (82372194), Tianjin Health Science and Technology
   Project of China (TJWJ2023MS012).

Data availability

   The data involved in our work are available in the TCGA
   ([159]https://gdc.cancer.gov/cancer-genome-atlas-tcga), ICGC
   ([160]https://dcc.icgc.org/), and GEO
   ([161]https://www.ncbi.nlm.nih.gov/geo/).

CRediT authorship contribution statement

   Youcheng Zhang: Writing – original draft, Methodology, Data curation.
   Dapeng Chen: Writing – original draft, Methodology. Bing Ang: Data
   curation. Xiyue Deng: Software. Bing Li: Validation. Yi Bai: Software.
   Yamin Zhang: Writing – review & editing, Methodology, Investigation.

Declaration of competing interest

   The authors declare that they have no known competing financial
   interests or personal relationships that could have appeared to
   influence the work reported in this paper.

Footnotes

   ^Appendix A

   Supplementary data to this article can be found online at
   [162]https://doi.org/10.1016/j.heliyon.2024.e37711.

Appendix A. Supplementary data

   The following are the Supplementary data to this article.
   Multimedia component 1
   [163]mmc1.docx^ (22.9MB, docx)
   Table S1

   The 25 crucial necroptosis regulators extracted from earlier reviews.
   [164]mmc2.xlsx^ (9KB, xlsx)
   Table S2

   The six modeling genes' sequences of primers. Table S3:The DEGs between
   hepatoma cells from cluster1 and other hepatoma cells. Table S4: The
   univariate Cox regression analysis results for the 38 DEGs.
   [165]mmc3.xlsx^ (8.8KB, xlsx)
   Multimedia component 4
   [166]mmc4.xlsx^ (18.6KB, xlsx)
   Multimedia component 5
   [167]mmc5.xlsx^ (9.6KB, xlsx)

Fig. 1,S2.

   Fig. 1,S2
   [168]Open in a new tab

   Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)
   analyses of genes upregulated in the high-risk group.

References