Skip to main content

ORIGINAL RESEARCH article

Front. Immunol., 15 June 2023
Sec. Cancer Immunity and Immunotherapy
This article is part of the Research Topic Roles of Pan-cell Death of Tumor Cells in Regulating Tumor Immune Microenvironment View all 14 articles

PANoptosis-based molecular subtyping and HPAN-index predicts therapeutic response and survival in hepatocellular carcinoma

Fei Song&#x;Fei Song1†Cheng-Gui Wang&#x;Cheng-Gui Wang1†Jia-Zhen Mao&#x;Jia-Zhen Mao1†Tian-Lun WangTian-Lun Wang1Xiao-Liang LiangXiao-Liang Liang1Chen-Wei HuChen-Wei Hu1Yu ZhangYu Zhang1Lu HanLu Han2Zhong Chen*Zhong Chen1*
  • 1Department of Hepatobiliary Surgery, Affiliated Hospital of Nantong University, Medical School of Nantong University, Nantong, China
  • 2Jiangsu Vocational College of Medicine, Yancheng, China

Background: Hepatocellular carcinoma (HCC) is a highly prevalent and fatal cancer. The role of PANoptosis, a novel form of programmed cell death, in HCC is yet to be fully understood. This study focuses on identifying and analyzing PANoptosis-associated differentially expressed genes in HCC (HPAN_DEGs), aiming to enhance our understanding of HCC pathogenesis and potential treatment strategies.

Methods: We analyzed HCC differentially expressed genes from TCGA and IGCG databases and mapped them to the PANoptosis gene set, identifying 69 HPAN_DEGs. These genes underwent enrichment analyses, and consensus clustering analysis was used to determine three distinct HCC subgroups based on their expression profiles. The immune characteristics and mutation landscape of these subgroups were evaluated, and drug sensitivity was predicted using the HPAN-index and relevant databases.

Results: The HPAN_DEGs were mainly enriched in pathways associated with the cell cycle, DNA damage, Drug metabolism, Cytokines, and Immune receptors. We identified three HCC subtypes (Cluster_1, SFN+PDK4-; Cluster_2, SFN-PDK4+; Cluster_3, SFN/PDK4 intermediate expression) based on the expression profiles of the 69 HPAN_DEGs. These subtypes exhibited distinct clinical outcomes, immune characteristics, and mutation landscapes. The HPAN-index, generated by machine learning using the expression levels of 69 HPAN_DEGs, was identified as an independent prognostic factor for HCC. Moreover, the high HPAN-index group exhibited a high response to immunotherapy, while the low HPAN-index group showed sensitivity to small molecule targeted drugs. Notably, we observed that the YWHAB gene plays a significant role in Sorafenib resistance.

Conclusion: This study identified 69 HPAN_DEGs crucial to tumor growth, immune infiltration, and drug resistance in HCC. Additionally, we discovered three distinct HCC subtypes and constructed an HPAN-index to predict immunotherapeutic response and drug sensitivity. Our findings underscore the role of YWHAB in Sorafenib resistance, presenting valuable insights for personalized therapeutic strategy development in HCC.

GRAPHICAL ABSTRACT
www.frontiersin.org

Graphical Abstract Graphical abstract of the study.

Introduction

Hepatocellular carcinoma (HCC), the world’s third most common solid malignant tumor, is currently facing a critical situation (1, 2). In 2020, HCC was responsible for the second-highest cancer-related deaths worldwide, following only lung cancer (3, 4). Unfortunately, this proportion has steadily increased yearly, ranked third in 2017 and fourth in 2015 (5). The specific heterogeneity of liver cancer patients, coupled with the limited number of tests and treatments available, is a significant factor contributing to this situation (6, 7). Although surgery remains the primary treatment for primary liver cancer, it is only accessible to a small percentage of patients.

Moreover, even after successful surgical treatment, some patients remain at risk of recurrence and metastasis for several years thereafter (8). Despite remarkable advancements in the treatment of various cancers, such as lung cancer and melanoma, the effectiveness of immunotherapy in treating hepatocellular carcinoma (HCC) has been unsatisfactory. Specifically, the objective response rate (ORR) of PD-1 immune checkpoint blockers (ICBs) for advanced HCC hovers around 15%, which is not sufficient (9, 10). Similarly, sorafenib, a small representative molecule targeted drug, has limited survival benefits for patients with advanced liver cancer due to tumor resistance (11, 12). As such, there is an urgent need to accurately and effectively screen patients suitable for ICBs or targeted drug sensitivity, enabling them to receive the most suitable treatment.

In recent years, scholars have defined a novel cell death pathway, PANoptosis (‘P’ for Pyroptosis; ‘A’ for Apoptosis; ‘N’ for Necroptosis), which has been defined by scholars (1315). This pathway involves the activation of a cytoplasmic multiprotein complex called PANoptosome, which can trigger multiple forms of programmed cell death, including pyroptosis, apoptosis, and necroptosis (1618). The dysregulation of PANoptosis has been associated with various human diseases, including autoinflammatory diseases, cancer, and infectious and metabolic disorders. Some biomarkers associated with PANoptosis, including NLRP3, caspase-1 for pyroptosis, ZBP1, IRF1, caspase-8 for apoptosis, and RIPK3/RIPK1 for necroptosis, have shown considerable benefits in suppressing cancer (1921). For instance, IRF1 functions in both myeloid and epithelial cells to counteract AOM/DSS-induced colorectal tumorigenesis, while RIPK3 activation in colon cancer cells leads to increased cytokine expression in the tumor microenvironment, contributing to robust cytotoxic anti-tumor immunity (19, 22). It is widely recognized that cell death resistance is a hallmark feature of hepatocellular carcinoma, and tumor cells have developed various strategies, such as the loss of TP53 tumor suppressor function, to limit apoptosis, which also plays a pivotal role in the failure of traditional cancer treatment (23, 24).

Cancer immunotherapy is a promising modality that stimulates the immune system to eliminate cancer cells with minimal side effects by modulating inherent immunosurveillance (25). Although some immune checkpoint blockade (ICB) therapies, particularly anti-PD-L1/PD-1, have shown clinical efficacy for patients with advanced stages of cancer, the objective response rate and survival benefits remain limitation (26, 27). One important reason for this is the inability of ICBs to induce programmed cell death (PCD) is essential for organismal development, host defense against pathogens, and maintaining homeostasis (28). However, resistance to PCD has been shown to promote tumor development, highlighting the need for novel PCD-based cancer therapies (29, 30). As a pivotal inflammatory PCD pathway, PANoptosis possesses critical features of pyroptosis, apoptosis, and necroptosis, which cannot be accounted for by any of these three PCD pathways alone (31). PANoptosis triggers systematic inflammation by releasing pro-inflammatory intracellular contents, making it a promising avenue for solid tumor immunotherapy (32). Thus, a deeper understanding of the mechanisms underlying PANoptosis can offer new opportunities to develop effective strategies for hepatocellular carcinoma immunotherapy.

There is a gap in current research on the role of PANoptosis in HCC. In the present study, we conducted a Consensus-Cluster-Plus analysis to identify three subgroups based on differentially expressed genes associated with PANoptosis and HCC (HPAN_DEGs). We then investigated these subgroups’ immune profiles and mutational landscape and constructed a PANoptosis risk score model (HPAN-index) for HCC. The HPAN-index can be used to grade the prognostic risk of HCC and to predict response to immunotherapy and chemotherapy drugs. Furthermore, we developed an integrated scoring nomogram to improve prognostic stratification and predictive accuracy for individual patients. Finally, we validated the drug response in different HPAN-index groups using public databases and in vitro trials, highlighting the enormous clinical potential of our findings in improving personalized decision-making for immunotherapy in HCC (Graphical abstract of the study).

Methods

Data acquisition and preprocessing

Data were obtained from the Cancer Genome Atlas (TCGA) in the training set, whereas the validation set sample data was sourced from the International Cancer Genome Consortium (ICGC) database. The test set sample data (GSE14520) was acquired from the Gene Expression Omnibus (GEO) database (3336). Additionally, the GSE109211 cohort was used as a dataset for Sorafenib resistance validation, and the GSE100797 and GSE93157 cohorts were used as datasets for immunotherapy sensitivity evaluation (3739). Please refer to Supplementary Table 1 for detailed information on the data.

To generate the PANoptosis gene list, we merged the gene lists of pyroptosis, apoptosis, and necroptosis while eliminating any redundant genes. Specifically, the pyroptosis gene list was retrieved from the Reactome pathway database, while the apoptosis gene list was integrated from three separate gene lists obtained from the AmiGO2, Reactome, and KEGG pathway databases, respectively (40, 41). Furthermore, the necroptosis gene list was sourced from the AmiGO2 database. After compiling the individual gene lists, a total of 277 non-redundant genes were identified and included in subsequent analyses (Supplementary Figure 1).

Identification of differentially expressed genes associated with HCC and PANoptosis

In our study, differential analysis was performed using the “limma” package (version 3.40.6) to identify genes differentially expressed between normal and cancer groups based on the data obtained from the TCGA and ICGC databases. To this end, we obtained the expression spectrum dataset and utilized the “lmFit” function to perform multiple linear regression. Next, the “eBays” function was utilized to calculate moderated t-statistics, F-statistics, and log-odds of differential expression via empirical Bayes moderation of standard errors directed toward an anticipated value. Subsequently, we identified the significance of variations for each gene (42). The selection criteria for identifying differentially expressed genes (DEGs) were P<0.05 and |log2FC|>1.5.

Unsupervised clustering of HCC-PANoptosis-related model genes

We performed consensus clustering analysis to identify unknown hepatocellular carcinoma (HCC) subtypes using the “Consensus-Cluster-Plus” package and model genes (43). The clustering was executed with a 1-Pearson correlation distance, and 80% resampling of the sample, and the process was repeated ten times. Empirical cumulative distribution function plots were utilized to determine the optimal number of clusters.

Functional enrichment analysis

To carry out GO and KEGG functional enrichment analyses, we employed the R packages “org.Hs.eg.db” and “clusterProfiler” (version 3.14.3). Initially, genes were annotated with GO terms using “org.Hs.eg.db” and mapped to a background set. Subsequently, the “clusterProfiler” package was utilized for GO and KEGG enrichment analyses, obtaining gene set enrichment results. In both cases, the minimum and maximum gene set sizes were set at 5 and 5000, respectively. We acquired the latest KEGG pathway gene annotations through the KEGG REST API and mapped them to a background set. Statistical significance was determined by a P value of < 0.05 and a false discovery rate (FDR) of < 0.25 for both analyses (44).

For Gene Set Enrichment Analysis (GSEA), we obtained subset collections from the Molecular Signatures Database to evaluate the relevant pathways and molecular mechanisms based on gene expression profiles and phenotype grouping (41). We performed 1000 permutations to obtain statistically significant results by P value of < 0.05 and FDR of < 0.25.

Somatic mutation analysis

To evaluate somatic mutations and assess tumor mutation burden (TMB), we utilized the “maftools” R package (45, 46). Somatic mutation data was obtained from the TCGA database and analyzed to identify non-synonymous somatic mutations. We then calculated TMB scores by dividing the number of non-synonymous somatic mutations by the total size of the genome in megabases.

Immune landscape analysis

The Tumor Immune Dysfunction and Exclusion (TIDE) framework is a computational tool that evaluates the potential for tumor immune evasion using gene expression profiles of cancer samples (47, 48). TIDE scores computed for each tumor sample serve as biomarkers to predict the response to immune checkpoint blockade, including anti-PD1 and anti-CTLA4, across different cancer types. We employed five algorithms to evaluate immune cell infiltration in the tumor microenvironment: TIMER, EPIC, xCELL, CIBERSORT, and MCPcount (49, 50). These algorithms enable a comprehensive evaluation of the immune cell landscape in the tumor microenvironment.

Chemotherapy response and small-molecule drugs

Data from the Genomics of Drug Sensitivity in Cancer (GDSC) database were analyzed to predict chemotherapy response in HCC patients (51). The half-maximal inhibitory concentration (IC50) calculated using the “pRRophetic” R package was used to indicate response to chemotherapy drugs (52). To identify potential new targets for HCC treatment, the gene expression profiles of high-risk and low-risk patient groups were compared using the Connectivity Map (CMap) reference dataset (53). Specifically, differentially expressed genes were identified and ranked based on their enrichment in the CMap dataset. A drug was considered a potential target if the enrichment score was between -1 and 0 and the adjusted p-value was less than 0.05.

Survival analysis and machine learning

We established a Lasso regression model using the “glmnet” package and utilized 10-fold cross-validation to select the optimal Lambda value, enhancing the interpretability and predictive accuracy of the model. The Lambda value of 0.0024 was optimal for minimizing the cross-validation error. We determined the coefficients of each gene using multivariate Cox analysis and generated the final regression model with the selected Lambda value. At the Lambda value of 0.0024, IRAK1, PSMD11, CHMP2A, PTRH2, SFN, YWHAB, PSMD3, TP53BP2, and PSMA4 were identified as the most important genes for predicting STATUS, with a calculated scoring formula of:.

HPANi=i=19βiEi

We initially divided patients into two groups based on the risk coefficient value using the percentile (50%) and classified them as either the high HPAN-index or low HPAN-index groups. Subsequently, we used the “survfit” function in the R software package “survival” to analyze the prognostic differences between the two groups. The log-rank test method was employed to evaluate the significance of the prognostic differences between the samples in different groups.

Cell proliferation, western blot, and invasion assays

The inhibitory effect of sorafenib on cell growth was assessed using the Cell Counting Kit-8 (CCK-8, Dojindo Kumamoto, Japan). Cells were plated at a density of 5,000 cells per well in 96-well plates. Following an initial 8-hour incubation period, the cells were treated with sorafenib at the prescribed doses or left untreated for 48 hours. Specific monitoring steps can be referred to in the instructions provided by the CCK-8 kit. Western blotting was carried out in the manner previously mentioned (54).

For the Matrigel invasion experiment, 1:8 diluted Matrigel matrix gel coating from Corning (ME) was applied to the chamber. DMEM plates without FBS were utilized to inject 2 × 106 cells per group. DMEM supplemented with 20% fetal bovine serum was added to the lower chamber, while mitomycin C was administered to the upper chamber to prevent cell proliferation. After a 48-hour incubation, the submembrane surface-invading tumor cells were fixed with 4% methanol and stained with 0.1% crystal violet. Each sample was counted across × 100 microscopic fields. All assays were performed in triplicate to ensure reliability.

Construction of sorafenib-resistant cell lines and RNA interference

The resis-PLC cells were generated from PLC cells using a protocol involving continuous exposure to increasing concentrations of sorafenib, followed by stepwise selection (55). The cells were collected every 3-4 days, passaged, and cultured in DMEM media containing progressively higher concentrations of sorafenib until they could grow steadily in its presence.

Small interfering RNA (siRNA) oligonucleotides specific to the target gene were used to knockdown expression. Cells were transfected with siRNA oligonucleotides using Lipofectamine 3000 (Invitrogen) according to the manufacturer’s protocol. The siRNA sequences for the YWHAB gene were as follows: si-1 5’-GCTGAATTGGATACGCTGAAT-3’, si-2 5’-CCAATGCTACACAACCAGAAA-3’, and si-NC 5’-UUCUCCGAACGUGUCACGUdTdT-3’. The RNA duplexes were synthesized by Genomeditech (Shanghai, China). Knockdown efficiency was assessed by western blotting.

Statistical analysis

Statistical evaluations were conducted utilizing R software (v.4.1.0). Continuous variables were displayed as mean ± standard deviation (SD), while categorical variables were shown as frequency (percentage). The Student’s t-test or Wilcoxon test examined differences between two groups concerning continuous variables contingent upon data normality assumptions. The chi-square or Fisher’s exact test was applied to categorical variables based on anticipated frequency counts. A two-sided P-value below 0.05 was deemed statistically significant across all tests. The Kaplan-Meier technique was implemented for survival assessment, and the log-rank test was adopted to compare group variations. Multivariate survival analysis utilized Cox proportional hazards regression. All analyses were conducted by a professional statistician with over 5 years of experience.

Results

Identification and functional analysis of PANoptosis-associated differentially expressed genes for HCC

The PANoptosis gene set consisted of pyroptosis (27 genes), apoptosis (259 genes), and necroptosis (15 genes) (Figures 1A, B and Supplementary Figure S1A). The differentially expressed genes for HCC comprised 5663 differentially expressed genes for hepatocellular carcinoma screened by the TCGA database and 4587 differentially expressed genes for hepatocellular carcinoma screened by the IGCG database (Figures 1C, D and Supplementary Figure S1B). We mapped the three gene sets screened for TCGA, ICGC, and PANoptosis to obtain 69 genes, defined as PANoptosis-associated differentially expressed genes for HCC (HPAN_DEGs) (Figures 1E, F).

FIGURE 1
www.frontiersin.org

Figure 1 Identification of PANoptosis-associated differential genes for HCC. (A) Concept drawing of PANoptosis (Fig-draw website, ID: YPOYA779c7); (B) The PANoptosis gene list; (C, D) Heatmap of the top 50 up- and down-regulated DEGs between HCC and normal tissue in the TCGA and ICGC databases; (E) The Vene diagram is composed of the differential genes of TCGA and ICGC respectively and the PANoptosis related dataset; (F) Protein–protein interactions among the PANoptosis-associated differential genes for HCC (HPAN_DEGs).

We performed GO, KEGG, and GSEA enrichment analyses to investigate the biological functions and related signaling pathways of HPAN_DEGs. Bioprocess (BP) analysis revealed that HPAN_DEGs are mainly enriched in signal transduction, cell communication, interleukin-1-mediated signaling pathway, and regulation of RNA stability (FDR<0.1, p value<0.05, Figure 2A). Molecular functional (MF) analysis revealed that HPAN_DEGs were mainly enriched in protein binding, enzyme binding, enzyme regulator activity, and transcription factor binding (FDR<0.1, p-value<0.05, Figure 2B). Cell composition (CC) analysis showed that HPAN_DEGs were mainly enriched in proteasome complex, endopeptidase complex, peptidase complex, and cytosol (FDR<0.1, p-value<0.05, Figure 2C). KEGG enrichment analysis suggested that HPAN_DEGs were mainly enriched in Proteasome, Apoptosis, Cell cycle, Necroptosis, p53 signaling pathway, Platinum drug resistance, and other signaling pathways (FDR<0.1, p value<0.05, Figure 2D).

FIGURE 2
www.frontiersin.org

Figure 2 GO/KEGG/GSEA enrichment analysis of the HPAN_DEGs. (A–C) GO enrichment analyses based on the HPAN_DEGs; (D) KEGG enrichment analyses based on the HPAN_DEGs; (E–G) GSEA analysis based on Hallmark, Reactome and KEGG datasets respectively.

To further clarify the biological functions undertaken by HPAN_DEGs, we conducted a GSEA analysis of HPAN_DEGs utilizing the KEGG, Hallmark, and Reactome datasets, respectively. The results showed that HPAN_DEGs were mainly enriched in the P53 signaling pathway (Hallmark and KEGG), Reactive oxygen species pathway (Hallmark), Cell cycle (Reactome and KEGG), DNA repair (Hallmark and Reactome) (NES>1, p-value <0. 05 and FDR<0.25, Figures 2E–G). Overall, our study identified 69 PANoptosis-associated differentially expressed genes for HCC (HPAN_DEGs) and revealed their enrichment in various biological functions and signaling pathways, such as the P53 signaling pathway, DNA repair, and cell cycle, indicating their potential involvement in tumor growth, metastasis, and drug resistance.

HPAN_DEG expression profiling identifies three HCC subtypes with distinct prognoses

We applied consistent clustering analysis to group the HCC cohort of TCGA based on information from the expression profiles of 69 HPAN_DEGs. When the value of K was taken as 3, the average consistency within the group was higher while ensuring that the area under the CDF curve line was as large as possible (Figures 3A, B and Supplementary Figure S1C). We named cohorts Cluster_1 (n = 131), Cluster_2 (n = 160), and Cluster_3 (n = 74) (Figure 3C). We then compared the expression levels of PAN apoptotic genes between the three Clusters and found that the PDK4 gene was up-regulated in Cluster_2 compared to Cluster_1 and Cluster_3, while SFN expression was down-regulated in the other two groups relative to Cluster_1 (Figure 3D and Supplementary Figure S1D). Additionally, survival analysis demonstrated that these three subtypes of HPAN_DEGs exhibit distinct clinical prognostic outcomes, with Cluster_1 having the poorest overall survival rate, Cluster_2 having the best, and Cluster_3 falling in between the two (Figure 3E). In summary, this study applied clustering analysis to identify three distinct subtypes of HCC based on the expression profiles of 69 HPAN_DEGs, which exhibited differential expression of PAN apoptotic genes and distinct clinical prognostic outcomes.

FIGURE 3
www.frontiersin.org

Figure 3 HPAN_DEG expression profiling identifies three HCC subtypes with distinct prognoses. (A, B) Assessment of average consistency within clusters and assessment of area under the CDF curve line when k = 2 to 10; (C) The training cohort was divided into three HCC subtypes by consensus clustering. (D) A heatmap displayed the expression of HPAN_DEGs in different HCC subtypes; (E) Kaplan-Meier survival analysis between three subtypes of HPAN_DEGs.

Distinct immunological profiles and mutational landscapes in HPAN_DEGs subgroups

Previous studies suggest that PANoptosis may influence tumor mutation and immune infiltration. To assess the immunological profile among subgroups of HPAN_DEGs, we performed an immunological landscape analysis of each of the three subgroups using several immunological algorithms, including CIBERSORT, ESTIMATE, and xCELL. The waterfall diagram of Figure 4A illustrates the distribution of the 22 immune cells in the TCGA training set. Then, we evaluated the Immune-Score, Stromal-Score, and Microenvironment-Score for HPAN_DEGs subgroups (Figures 4B, C). Our study indicates that Cluster_2 is significantly different from the other two groups in the term of Immune-Score and Stromal-Score (Cluster_2 vs. Cluster_1, 0.04 ± 0.05 vs. 0.06 ± 0.07, P = 0.03; Cluster_2 vs. Cluster_3, 0.04 ± 0.05 vs. 0.10 ± 0.15, P = 0.0025, Immune-Score) (Cluster_2 vs. Cluster_1, 0.12 ± 0.06 vs. 0.06 ± 0.04, P = 7.9E-19; Cluster_2 vs. Cluster_3, 0.12 ± 0.06 vs. 0.07 ± 0.05, P = 0.00000000073, Stroma-Score). In the assessment of the Microenvironment-Score, we found that Cluster_2 and Cluster_3 were not statistically different, while Cluster_1 was significantly different from the other two groups (Cluster_1 vs. Cluster_2, 0.12 ± 0.085 vs. 0.16 ± 0.09, P = 0.0000017; Cluster_1 vs. Cluster_3, 0.12 ± 0.08 vs. 0.17 ± 0.16, P = 0.03, Figure 4F). We also evaluated the gene expression of immune checkpoints among HPAN_DEGs subgroups, namely PD-1 (PDCD1), PD-L1 (CD274), PD-L2 (PDCD1LG2), CTLA4, LAG3, TIGIT, HAVCR2, and found that the expression of these immune checkpoints was up-regulated in Cluster_1 and Cluster_3, and downregulated in Cluster_2 (all P < 0.05, Figure 4H).

FIGURE 4
www.frontiersin.org

Figure 4 Distinct immunological profiles and mutational landscapes in HPAN_DEGs subgroups. (A) Waterfall diagram of the distribution of the 24 immune cells in the training set; (B, C) Immuno-score and Stromal-score for 3 subgroups of HPAN_DEGs. (D) Waterfall maps of the somatic mutations in different HPAN_DEGs subtypes; (E) Tumor mutation burden in different HPAN_DEGs subtypes; (F) Microenvironment score -score for 3 subgroups of HPAN_DEGs; (G) Characteristics of gene mutation in different HPAN_DEGs subtypes; (H) Expression of immune checkpoints in different HPAN_DEGs subtypes. *p < 0.05, **p < 0.01, ***p < 0.001.

Subsequently, we demonstrated the somatic mutational landscape among the HPAN_DEGs subgroups. The top 15 mutated genes in the three subgroups were TP53/CTNNB1/TTN/MUC16/ALB/PCLO/MUC4/RYR2/ABCA13/APOB/CSMD3/LRP1B/FLG/OBSCN/AXIN1. The gene with the highest mutation rate was TP53, which varied among the three clusters, with Cluster_1 (44.79%) having a higher mutation rate than Cluster_2 (23.96%) and Cluster_3 (29.17%), respectively (Figure 4D). Tumor mutational load (TMB) is an essential indicator of the number of mutations in cancer and a novel marker for evaluating the efficacy of PD-1 antibody therapy. We compared the tumor mutational burden (TMB) of the three subgroups and found that the TMB of Cluster_2 was lower than that of Cluster_1 and Cluster_3, respectively (Cluster_2 vs. Cluster_1, 1.84 ± 1.34 vs. 2.85 ± 3.99, P = 0.0042; Cluster_2 vs. Cluster_3, 1.84 ± 1.34 vs. 2.50 ± 2.08, P = 0.00034, Figure 4E). We then applied a ternary diagram showing the distribution of mutant genes among different subgroups of HPAN_DEGs (Figure 4G). In summary, our findings demonstrated substantial discrepancies in the expression of immune checkpoints and the mutational landscape among the three clusters, which could potentially have crucial ramifications in cancer immunotherapy.

Construction of a HPAN_DEGs-based PANoptosis risk score model for prognostic assessment in HCC

To further evaluate the impact of HPAN_DEGs on survival prognosis, we used LASSO, univariate and multivariate regression to screen nine gene signatures with strong prognostic associations, and finally constructed a PANoptosis risk index for hepatocellular carcinoma (HPAN-index, Figures 5A–D). HPAN-index = 0.4142* IRAK1 + 0.78337*PSMD11 - 0.33085*CHMP2A + 0.66389*PTRH2 + 0.11779*SFN + 0.82753*YWHAB - 0.95811*PSMD3 - 0.25778 *TP5 3BP2 - 0.72582*PSMA4. We classified 365 patients with complete survival information in training set into high and low HPAN-index groups utilizing the HPAN-index score (high 183 vs. low 182). Kaplan-Meier analysis demonstrated a significantly better prognosis for the low-risk group (Median Survival Time, MST = 83.8 months) than the high-risk group (MST = 29.9 months) in the training set (P < 0.001, Figure 5E). We investigated the relationship between patient prognosis, gene expression, and HPAN-index and observed a significant decrease in survival as the HPAN-index increased. As expected, CHMP2A/PSMA4/PSMD3/TP53BP2 were protective factors, whose expression was downregulated with increasing HPAN-index, while IRAK1/PSMD11/PTRH2/SFN/YWHAB were risk factors (Figure 5F).

FIGURE 5
www.frontiersin.org

Figure 5 Construction of a HPAN_DEGs-based PANoptosis risk score model for prognostic assessment in HCC (HPAN-index). (A, B) Lasso regression analysis with 10-fold cross-validation resulted in 20 genes associated with survival (Lambda = 0.0239703477847909); (C, D) Univariate and multivariate Cox analyses further screened for nine PANoptosis-associated genes associated with survival; (E) Kaplan–Meier analyses demonstrate the prognostic significance of the HPAN-index model in the training set; (F) HPAN-index distribution, survival status of each patient, and heatmaps of prognostic 9-gene signature in the training set; (G) Receiver operator characteristic (ROC) analysis of the HPAN-index model in the training set; (H) Sankey diagrams shows the interrelationship between HPAN_DEGs subtypes, the risk groups of the HPAN-index and the individual clinical characteristics; (I) A nomogram was established to predict the prognostic of HCC patients.

In addition, we evaluated the area under the curve (AUC) of the HPAN-index as a predictive model, and the results suggested that the HPAN-index was highly accurate in predicting survival at 1, 3, and 5 years (Figure 5G). We applied Sankey diagrams to visualize the relationship between the risk groups of the HPAN-index and the individual clinical characteristics, suggesting that Cluster_1 mainly converges in the high HPAN-index group. In contrast, Cluster_2 mainly converges in the low HPAN-index group (Figure 5H). Interestingly, Stage I/II in TNM staging mainly converged in the low HPAN-index group, while Stage III/IV mainly converged in the high HPAN-index group. We constructed a nomograph based on the Cox regression analysis results and found that the HPAN-index was an independent risk factor (Figure 5I and Supplementary Figure S1E). In conclusion, the HPAN_DEGs-based PANoptosis risk score model (HPAN-index) constructed by LASSO regression, univariate and multivariate regression analysis, can accurately predict the survival prognosis of hepatocellular carcinoma patients and could be a potential independent risk factor for clinical decision-making.

Validation of HPAN-index as a prognostic predictor in HCC patients across multiple cohorts

To examine the repeatability of the model HPAN-index as a predictive model, we validated the model in the ICGC_HCC cohort (Validation set) and the GSE14520 cohort (Testing set). Applying the Kaplan-Meier analysis, we can observe a significant decrease in patient survival as the HPAN-index increases (Figures 6A, C). In the Validation set, the prognosis was significantly better in the low-HPAN-index group (MST = 66.7 months) than in the high-HPAN-index group (MST = 47.3 months, P < 0.001, Figure 6B), with similar results in the Testing set (P = 0.01, Figure 6D). In conclusion, the HPAN-index model demonstrated significant predictive value for patient survival in both the Validation and Testing sets, with higher HPAN-index scores indicating poorer prognosis.

FIGURE 6
www.frontiersin.org

Figure 6 Validation of HPAN-index as a prognostic predictor in HCC patients across multiple cohorts. (A) HPAN-index distribution, survival status of each patient, and heatmaps of prognostic 9-gene signature in the validation set (ICGC, n = 243); (B) Kaplan–Meier analyses demonstrate the prognostic significance of the HPAN-index model in the validation set (ICGC, n = 243); (C) HPAN-index distribution, survival status of each patient, and heatmaps of prognostic 9-gene signature in the testing set (GSE14520, n = 242); (D) Kaplan–Meier analyses demonstrate the prognostic significance of the HPAN-index model in the validation set (GSE14520, n = 242); (E) Box plot visualizes significantly different immune cells between different HPAN-index groups; (F) Correlation of HPAN-index with immune cell infiltration evaluated using CIBERSORT in the HCC; (G) Immuno-score, Stromal-score, and ESTIMATE-score between different HPAN-index groups. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

Immune cell landscape and molecular pathways associated with HPAN-index in HCC patients

To further investigate the immune characteristics of different HPAN-index groups, we employed five distinct immune algorithms, including CIBERSORT, ESTIMATE, TIDE, TIMER, and xCell, to assess the relationship between HPAN-index and the immune microenvironment. The outcomes obtained from Figures 6E, F demonstrated a significant correlation between the HPAN-index and the expression of diverse immune cells, including B cells naive, B cells memory, T cells CD8, T cells CD4 memory activated, T cells follicular helper, T cells regulatory, T cells gamma delta, NK cells resting, NK cells activated, Monocytes, Macrophages M0, Mast cells resting, and Neutrophils, with all P values less than 0.05. Moreover, the high HPAN-index group exhibited significantly higher Stromal-Score, Immune-Score, and ESTIMATE-Score compared to the low HPAN-index group (all P<0.05, Figure 6G). Notably, the Microsatellite Instability (MSI)-score, an essential indicator reflecting tumor genome stability, significantly correlated with immune checkpoint efficacy (56). Based on our analysis, we observed that the high HPAN-index group has a significantly higher MSI-score compared to the low HPAN-index group, with similar results for other scores such as TIDE-score, IFNG-score, Merck18-score, Dysfunction-score, Exclusion-score, MDSC-score, and TAM M2-score (all P < 0.05, Figures 7A–F). Compared to the low HPAN-index group, the high HPAN-index group exhibited a significant up-regulation in the expression levels of immune checkpoints, including CD274, CTLA4, LAG3, TIGIT, HAVCR2, PDCD1, and PDCD1LG2 (all P < 0.05, Figure 7H). This observation suggests that the high HPAN-index group may be more likely to benefit from immunotherapy than the low HPAN-index group.

FIGURE 7
www.frontiersin.org

Figure 7 Immune cell landscape and molecular pathways associated with HPAN-index in HCC patients. (A–F) TIDE, immune dysfunction, immune exclusion, CD274, Merck18, CD8, and IFNγ scores in low and high HPAN-index groups; (G) Top enriched pathways for genes with specific expression in the high and low HPAN-index groups.; (H) Differential expression of immune checkpoints between the high and low HPAN-index groups. *p < 0.05, **p < 0.01, ***p < 0.001.

To further elucidate the molecular mechanisms underlying the different HPAN-index groups, we performed GSVA enrichment analysis. The results revealed that the high HPAN-index group was significantly enriched in signaling pathways such as the T cell receptor pathway, chemokine pathway, B cell receptor pathway, DNA replication, and cell cycle. In contrast, the low HPAN-index group showed significant enrichment in signaling pathways related to fatty acid metabolism, drug metabolism, glycine serine, threonine metabolism, PPAR signaling pathway, and tryptophan metabolism (Figure 7G). We further evaluated drug sensitivity in different HPAN-index groups and found that small molecule inhibitors (JAK1_8709/KRAS_G12C/Linsitinib/Nilotinib/Oxaliplatin/Niraparib/Picolinic-acid/Selumetinib/Sorafenib) showed significantly higher sensitivity in the low HPAN-index group compared to the high HPAN-index group (all P < 0.05, Supplementary Figures S2A–I). In summary, patients with high HPAN-index may be more responsive to immune checkpoint therapy, while those with low HPAN-index may be better suited for a targeted drug.

Validation of HPAN-index and identification of key molecule YWHAB in sorafenib resistance

To further validate the accuracy of HPAN-index in predicting drug resistance and explore the key molecules, we combined the correlation analysis results of gene expression data from TCGA and protein-protein interaction (PPI) topological network analysis results to identify YWHAB with higher weight in HPAN-index (Figures 8A, B). Subsequently, we utilized the DepMap database to select PLCPRF5 cells with the highest YWHAB expression levels to construct sorafenib-resistant cells (resis-PLC) (Figures 8C, D). We were thrilled to discover that knocking down YWHAB in resis-PLC not only restored sensitivity to sorafenib but actually resulted in even greater sensitivity than the non-resistant control (IC50, NC_PLC 7.265 μM, resis-PLC 11.01 μM, siYWHAB 4.288 μM, Figures 8E, F). As shown in Figures 8G, H the knockdown of YWHAB restored the ability of sorafenib to inhibit the invasion of resis-PLC. Furthermore, we further validated HPAN-index in public databases, and the results showed that the Low HPAN-index group was more sensitive to sorafenib and less responsive to immunotherapy than the High HPAN-index group (Figures 8I–L). In summary, the HPAN-index exhibits considerable advantages in predicting the response to immunotherapy and the sensitivity to targeted drugs, with YWHAB potentially playing a crucial role in the HPAN-index’s functionality.

FIGURE 8
www.frontiersin.org

Figure 8 Validation of HPAN-index and identification of key molecule YWHAB in sorafenib resistance. (A) Chord Diagram displaying the relationship between the hub HPAN_DEGs expression in the training set; (B) Protein-protein interaction network diagram between the hub HPAN_DEGs; (C) Bar graph showing the expression of YWHAB in various hepatocellular carcinoma cell lines; (D) Construction of a sorafenib-resistant PLC cell line using an incremental drug concentration method (continuous induction at low and increasing concentrations); (E) Western blot analysis of YWHAB knockdown in sorafenib-resistant cell line (resis-PLC); (F) Dose-response curve of sorafenib treatment in NC_PLC, resis-PLC, and resis-PLC with YWHAB knockdown (The IC50 values of sorafenib for NC_PLC, resis-PLC, and siYWHAB were 7.265 μM, 11.01 μM, and 4.288 μM, respectively); (G) Transwell assay for evaluating the effect of YWHAB knockdown on resis-PLC cell sensitivity recovery; (H) Quantification of invasive cells as a percentage; (I) Bar graph showing the HPAN-index of different samples from GSE109211 dataset; (J) Box plot analysis of HPAN-index for the Sorafenib-sensitive and Sorafenib-resistant groups; (K) Bar graph showing the HPAN-index of different samples from GSE100797 dataset; (L) Box plot analysis of HPAN-index for the two groups classified as the anti-PD-1 Therapy-responsive and non-responsive groups (blue) based on the GSE93157 dataset. Data are presented as mean ± standard deviation (SD) from three independent experiments, *p < 0.05, ***p < 0.001.

Discussion

As the clinical use of immunotherapeutic agents and molecularly targeted inhibitors continues to expand, liver cancer patients experience limited benefits compared to other malignancies such as melanoma, lung cancer, and kidney cancer (5759). It is primarily due to the unique characteristics of the liver, which has a remarkable regenerative capacity and serves a vital role in detoxification. Furthermore, malignant tumors originating from the liver are more heterogeneous than other tumors (60). Therefore, it is of utmost importance for scholars to address the urgent question of enabling liver cancer patients to select the most appropriate clinical drugs and achieve personalized and precise treatment for liver cancer.

In recent years, scholars have proposed the concept of PANoptosis, which highlights the complex interplay between different cell death pathways in regulating tumor development (13). It is now widely acknowledged that a single death pathway does not solely govern tumor progression but involves intricate crosstalk between various pathways (14, 19). This integration of functions often has significant implications for tumor resistance and the immune microenvironment. In this study, we developed exclusive models for hepatocellular carcinoma related to PANoptosis (HPAN-index). We further validated the predictive performance of these models in terms of prognosis, minor molecule drug sensitivity, and immunotherapy in both the validation and test sets. These findings highlight the crucial role of PANoptosis in the context of hepatocellular carcinoma and offer valuable insights for developing personalized and precise therapeutic strategies.

The team of Prof. Gao Q. identified three distinct subtypes of hepatocellular carcinoma (HCC): the Metabolic subtype, the Proliferative subtype, and the Tumor Microenvironment Dysregulation subtype (60). Our study conducted a consistency clustering analysis of 69 HPAN_DEGs in the training set and identified three subgroups that exhibit distinct characteristics regarding overall survival prognosis, mutational landscape, and immune infiltration. Cluster_1, which demonstrated the highest tumor mutational load, had the worst overall survival rate. However, Cluster_1 showed an advantage in immunotherapy with higher expression of its major immune checkpoints CD274/CTLA4/LAG3/TIGHT/HAVCR2/PDCD1/PDCD1LG2 compared to the other groups. On the other hand, Cluster_2 presented an opposite phenotype in terms of overall survival prognosis, mutational landscape, and immunomolecular profile compared to Cluster_1. We observed that Cluster_1 presented SFN+PDK4-, whereas Cluster_2 presented SFN-PDK4+. Our molecular characterization of these three subgroups revealed essential insights into the complex interplay between tumor mutational load, immune checkpoint expression, and prognosis, providing valuable information for developing personalized therapeutic strategies for hepatocellular carcinoma.

Previous studies have demonstrated that SFN is an oncogene, accelerating tumorigenesis and progression across various cancer types (61, 62). Prof. Masayuki Noguchi’s team identified that SFN specifically binds to ubiquitinated protease 8 (USP8) in lung adenocarcinoma cells, enhancing the stabilization of receptor tyrosine kinases (RTKs), including EGFR and MET, through abnormal regulation of USP8. These findings suggest that SFN may be a promising therapeutic target for lung adenocarcinoma. In line with this, our study also found that positive expression of SFN in Cluster_1 could indicate the potential for tumor proliferation in these patients. Pyruvate dehydrogenase kinase 4 (PDK4) encodes an enzyme that regulates cellular metabolism by inhibiting the phosphorylation of a key regulatory enzyme of glucose oxidation, pyruvate dehydrogenase complex (PDC) (63). High expression of PDK4 has been associated with altered metabolic pathways in tumor cells, including lactic acidification and malignant transformation (6466). In our study, the positive expression of PDK4 in Cluster_2 may suggest that the tumor type in these patients is associated with aberrant tumor cell metabolism. Overall, our study identified three subgroups based on HPAN_DEGs with distinct characteristics in terms of prognosis, mutational landscape, and immune infiltration. Furthermore, these subgroups exhibited differential expression of SFN and PDK4 genes. Our findings contribute to a better understanding of the biology of these tumor types and may provide a new basis for subgroup screening in hepatocellular carcinoma.

Using LASSO regression and univariate and multifactorial analyses, we screened nine genes strongly associated with prognosis and constructed the HPAN-index. With this model, we divide the cohort into High-index and Low-index groups, where the High-index group is mainly from Cluster_1, and the Low-index group is mainly from Cluster_2. We note that the High-index group showed a significant immune activation status. In contrast, the Low-index group showed a significant advantage in sensitivity to small molecule targeted drugs, which is consistent with Yutian Zou et al. (67). We found that the High-index group is mainly enriched in T cell receptor pathways, B cell receptor pathways, Chemokine pathways, DNA replication, and Cell cycle pathways. Professor Peter P. Lee’s research has shown that T/B receptor pathways are closely related to T cell activation and affect PD-1 expression on T cells (68). Combined with the higher expression of PANoptosis-related genes in the High-index group, we suggest that this immune activation state in the High-index group may be related to PANoptosis. Also of note, the Low-index group demonstrated significant sensitivity to some small molecule-targeted drugs, such as sorafenib, a first-line agent for the treatment of advanced primary liver cancer (69). It may be because the Low-index group is mainly enriched in Fatty acid metabolism, Drug metabolism, PPAR signaling pathway, Tryptophan metabolism, and other pathways closely related to tumor growth, metabolism, and drug resistance (70, 71). Overall, the HPAN-index may serve as an independent risk factor for predicting the prognosis of patients with hepatocellular carcinoma and as a strategy for selecting patients for immunotherapy and targeted therapeutic agents.

Apoptosis is one of the crucial mechanisms underlying tumor cell drug resistance (72). YWHAB is a gene in the human genome that encodes the 14-3-3 protein beta/alpha. The 14-3-3 protein family is a highly conserved group of molecular chaperones that participate in various cellular signaling and regulatory processes, such as metabolism, protein transport, signal transduction, apoptosis, and cell cycle (73). Silencing of YWHAB can increase the translocation of B-cell lymphoma 2 (BCL-2)-associated death promoter (BAD) from the cytoplasm to the mitochondria, thereby inducing cell apoptosis (74). The BCL-2 family is a critical group of molecules in the field of tumor drug resistance, and its mechanism of inducing drug resistance is mainly achieved by inhibiting the apoptotic pathway of tumor cells. Studies have shown that the ratio of BCL-2/Bax is higher in drug-resistant cells than in sensitive cells (75). YWHAB, as an anti-apoptotic protein, can also cause insulin resistance in cells by affecting mitochondrial polarization (76). Our research suggests that YWHAB may play an essential role in affecting cell drug sensitivity, providing insights to further study the mechanisms of drug resistance and develop new therapeutic strategies.

There are limitations to our study that should be acknowledged. We lack the necessary single-cell level sequencing data and spatial transcriptome data to comprehensively support our analysis of the immune landscape of hepatocellular carcinoma. As the immune microenvironment is a complex microscopic system, the information on the differences and interactions between cells is inevitably lost through macroscopic bulk-RNA-seq data analysis in isolation. Moreover, a large sample size of immunotherapy-related data for hepatocellular carcinoma is needed to validate the model. Therefore, we will seek to obtain such data to validate the model further. We take this opportunity to call on the scientific community to share data related to immunotherapy for liver cancer, thereby advancing the scientific understanding of this complex disease.

Conclusions

In conclusion, we screened for PANoptosis-associated differentially expressed genes (HPAN_DEGs) in hepatocellular carcinoma, which allowed us to identify three subgroups that exhibit distinct characteristics in terms of prognosis, mutational landscape, and immune infiltration. These subgroups also exhibited differential expression of SFN and PDK4, which may contribute to a better understanding of the biology underlying hepatocellular carcinoma. Additionally, we developed the HPAN-index, which is highly correlated with survival prognosis, sensitivity to small molecule-targeted drugs, and response to immunotherapy. We hope that applying this model will enable the identification of individuals more suitable for either immunotherapy or targeted therapy. This study provides a new strategy for the personalized and precise treatment of HCC and may shed light on future investigations into the mechanisms of PANoptosis in this disease.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found within the article/Supplementary Material.

Author contributions

FS performed the research and wrote the paper. FS, C-GW, and X-LL performed the data collection and normalization. C-GW and J-ZM carried out the in vitro validation. C-WH and T-LW participated in the coordination of the research. FS, YZ, and LH performed the statistical analysis. FS and ZC participated in the study design. ZC edited the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by the National Natural Science Foundation of China grants (81871927, 81070360), 2021 Changchun University of Chinese Medicine School-level Clinical Practice Teaching Reform Special Research Project (XJLCSJ202146), and the Nantong Hepatobiliary and Pancreatic Surgery Disease Research Center Construction Project (HS2015001).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1197152/full#supplementary-material

Supplementary Figure 1 | (A). Venn diagram illustrating the overlap between gene sets related to pyroptosis, apoptosis, and necroptosis; (B). Volcano plot displaying differentially expressed genes in hepatocellular carcinoma (HCC) derived from TCGA and ICGC databases; (C). Heatmap demonstrating sample clustering consistency; (D). Graphical representation of the proportion of SFN+ PDK4- and SFN-PDK4+ in three HPAN_DEGs subgroups; (E). A prognostic nomogram based on HPAN-index gene signatures was developed to predict the outcomes of HCC patients.

Supplementary Figure 2 | Efficacy of HPAN-index in predicting drug sensitivity (A-I) Box plots illustrating the comparison of IC50 values for various drugs between the high-HPAN-index group (depicted in red) and the low-HPAN-index group (depicted in green). P-values are presented in scientific notation for each comparison.

Abbreviations

HCC, Hepatocellular carcinoma; HPAN_DEGs, PANoptosis-associated differentially expressed genes for HCC; SFN, Stratifin; PDK4, Pyruvate Dehydrogenase Kinase 4; YWHAB, Tyrosine 3-Monooxygenase/Tryptophan 5-Monooxygenase Activation Protein Beta; ICBs, Immune checkpoint blockers; PCD, Programmed cell death; TCGA, The Cancer Genome Atlas; ICGC, International Cancer Genome Consortium; GEO, Gene Expression Omnibus; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; IRAK1, Interleukin 1 Receptor Associated Kinase 1; PSMD11, Proteasome 26S Subunit, Non-ATPase 11; CHMP2A, Charged Multivesicular Body Protein 2A; PTRH2, Peptidyl-TRNA Hydrolase 2; PSMD3, Proteasome 26S Subunit, Non-ATPase 3; TP53BP2, Tumor Protein P53 Binding Protein 2; PSMA4, Proteasome 20S Subunit Alpha 4; HPAN-index, PANoptosis risk index for hepatocellular carcinoma.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2021) 71:209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Llovet JM, Castet F, Heikenwalder M, Maini MK, Mazzaferro V, Pinato DJ, et al. Immunotherapies for hepatocellular carcinoma. Nat Rev Clin Oncol (2022) 19:151–72. doi: 10.1038/s41571-021-00573-2

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Sheppard S, Guedes J, Mroz A, Zavitsanou AM, Kudo H, Rothery SM, et al. The immunoreceptor NKG2D promotes tumour growth in a model of hepatocellular carcinoma. Nat Commun (2017) 8:13930. doi: 10.1038/ncomms13930

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Xia S, Pan Y, Liang Y, Xu J, Cai X. The microenvironmental and metabolic aspects of sorafenib resistance in hepatocellular carcinoma. EBioMedicine (2020) 51:102610. doi: 10.1016/j.ebiom.2019.102610

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Global Burden of Disease Cancer, C, Fitzmaurice C, Allen C, Barber RM, Barregard L, Bhutta ZA, et al Global, regional, and national cancer incidence, mortality, years of life lost, years lived with disability, and disability-adjusted life-years for 32 cancer groups, 1990 to 2015: a systematic analysis for the global burden of disease study. JAMA Oncol (2017) 3:524–48. doi: 10.1001/jamaoncol.2016.5688

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Jiang Y, Sun A, Zhao Y, Ying W, Sun H, Yang X, et al. Proteomics identifies new therapeutic targets of early-stage hepatocellular carcinoma. Nature (2019) 567:257–61. doi: 10.1038/s41586-019-0987-8

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Zender L, Spector MS, Xue W, Flemming P, Cordon-Cardo C, Silke J, et al. Identification and validation of oncogenes in liver cancer using an integrative oncogenomic approach. Cell (2006) 125:1253–67. doi: 10.1016/j.cell.2006.05.030

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Li X, Yao W, Yuan Y, Chen P, Li B, Li J, et al. Targeting of tumour-infiltrating macrophages via CCL2/CCR2 signalling as a therapeutic strategy against hepatocellular carcinoma. Gut (2017) 66:157–67. doi: 10.1136/gutjnl-2015-310514

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Finn RS, Ikeda M, Zhu AX, Sung MW, Baron AD, Kudo M, et al. Phase ib study of lenvatinib plus pembrolizumab in patients with unresectable hepatocellular carcinoma. J Clin Oncol (2020) 38:2960–70. doi: 10.1200/JCO.20.00808

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Yang JD, Heimbach JK. New advances in the diagnosis and management of hepatocellular carcinoma. BMJ (2020) 371:m3544. doi: 10.1136/bmj.m3544

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Dong ZR, Ke AW, Li T, Cai JB, Yang YF, Zhou W, et al. CircMEMO1 modulates the promoter methylation and expression of TCF21 to regulate hepatocellular carcinoma progression and sorafenib treatment sensitivity. Mol Cancer (2021) 20:75. doi: 10.1186/s12943-021-01361-3

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Wu H, Wang T, Liu Y, Li X, Xu S, Wu C, et al. Mitophagy promotes sorafenib resistance through hypoxia-inducible ATAD3A dependent axis. J Exp Clin Cancer Res (2020) 39:274. doi: 10.1186/s13046-020-01768-8

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Pan H, Pan J, Li P, Gao J. Characterization of PANoptosis patterns predicts survival and immunotherapy response in gastric cancer. Clin Immunol (2022) 238:109019. doi: 10.1016/j.clim.2022.109019

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Chen W, Gullett JM, Tweedell RE, Kanneganti TD. Innate immune inflammatory cell death: PANoptosis and PANoptosomes in host defense and disease. Eur J Immunol (2023):e2250235. doi: 10.1002/eji.202250235

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Karki R, Sharma BR, Tuladhar S, Williams EP, Zalduondo L, Samir P, et al. Synergism of TNF-alpha and IFN-gamma triggers inflammatory cell death, tissue damage, and mortality in SARS-CoV-2 infection and cytokine shock syndromes. Cell (2021) 184:149–168.e117. doi: 10.1016/j.cell.2020.11.025

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Banoth B, Tuladhar S, Karki R, Sharma BR, Briard B, Kesavardhana S, et al. ZBP1 promotes fungi-induced inflammasome activation and pyroptosis, apoptosis, and necroptosis (PANoptosis). J Biol Chem (2020) 295:18276–83. doi: 10.1074/jbc.RA120.015924

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Karki R, Lee S, Mall R, Pandian N, Wang Y, Sharma BR, et al. ZBP1-dependent inflammatory cell death, PANoptosis, and cytokine storm disrupt IFN therapeutic efficacy during coronavirus infection. Sci Immunol (2022) 7:eabo6294. doi: 10.1126/sciimmunol.abo6294

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Kuriakose T, Man SM, Malireddi RK, Karki R, Kesavardhana S, Place DE, et al. ZBP1/DAI is an innate sensor of influenza virus triggering the NLRP3 inflammasome and programmed cell death pathways. Sci Immunol (2016) 1. doi: 10.1126/sciimmunol.aag2045

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Karki R, Sharma BR, Lee E, Banoth B, Malireddi RKS, Samir P, et al. Interferon regulatory factor 1 regulates PANoptosis to prevent colorectal cancer. JCI Insight (2020) 5. doi: 10.1172/jci.insight.136720

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Lee S, Karki R, Wang Y, Nguyen LN, Kalathur RC, Kanneganti TD, et al. AIM2 forms a complex with pyrin and ZBP1 to drive PANoptosis and host defence. Nature (2021) 597:415–9. doi: 10.1038/s41586-021-03875-8

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wang Y, Pandian N, Han JH, Sundaram B, Lee S, Karki R, et al. Single cell analysis of PANoptosome cell death complexes through an expansion microscopy method. Cell Mol Life Sci (2022) 79:531. doi: 10.1007/s00018-022-04564-z

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Park HH, Kim HR, Park SY, Hwang SM, Hong SM, Park S, et al. RIPK3 activation induces TRIM28 derepression in cancer cells and enhances the anti-tumor microenvironment. Mol Cancer (2021) 20:107. doi: 10.1186/s12943-021-01399-3

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell (2011) 144:646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Muller M, Bird TG, Nault JC. The landscape of gene mutations in cirrhosis and hepatocellular carcinoma. J Hepatol (2020) 72:990–1002. doi: 10.1016/j.jhep.2020.01.019

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Morad G, Helmink BA, Sharma P, Wargo JA. Hallmarks of response, resistance, and toxicity to immune checkpoint blockade. Cell (2021) 184:5309–37. doi: 10.1016/j.cell.2021.09.020

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Blomberg OS, Spagnuolo L, Garner H, Voorwerk L, Isaeva OI, van Dyk E, et al. IL-5-producing CD4(+) T cells and eosinophils cooperate to enhance response to immune checkpoint blockade in breast cancer. Cancer Cell (2023) 41:106–123.e110. doi: 10.1016/j.ccell.2022.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Tichet M, Wullschleger S, Chryplewicz A, Fournier N, Marcone R, Kauzlaric A, et al. Bispecific PD1-IL2v and anti-PD-L1 break tumor immunity resistance by enhancing stem-like tumor-reactive CD8(+) T cells and reprogramming macrophages. Immunity (2023) 56:162–179 e166. doi: 10.1016/j.immuni.2022.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Nagata S, Tanaka M. Programmed cell death and the immune system. Nat Rev Immunol (2017) 17:333–40. doi: 10.1038/nri.2016.153

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Liu L, Li H, Hu D, Wang Y, Shao W, Zhong J, et al. Insights into N6-methyladenosine and programmed cell death in cancer. Mol Cancer (2022) 21:32. doi: 10.1186/s12943-022-01508-w

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Meng Q, Ding B, Ma P, Lin J. Interrelation between programmed cell death and immunogenic cell death: take antitumor nanodrug as an example. Small Methods (2023):e2201406. doi: 10.1002/smtd.202201406

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Karki R, Sundaram B, Sharma BR, Lee S, Malireddi RKS, Nguyen LN, et al. ADAR1 restricts ZBP1-mediated immune response and PANoptosis to promote tumorigenesis. Cell Rep (2021) 37:109858. doi: 10.1016/j.celrep.2021.109858

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Zhao P, Wang M, Chen M, Chen Z, Peng X, Zhou F, et al. Programming cell pyroptosis with biomimetic nanoparticles for solid tumor immunotherapy. Biomaterials (2020) 254:120142. doi: 10.1016/j.biomaterials.2020.120142

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Liu J, Lichtenberg T, Hoadley KA, Poisson LM, Lazar AJ, Cherniack AD, et al. An integrated TCGA pan-cancer clinical data resource to drive high-quality survival outcome analytics. Cell (2018) 173:400–416.e411. doi: 10.1016/j.cell.2018.02.052

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Consortium, I. T. P.-C. A. o. W. G Pan-cancer analysis of whole genomes. Nature (2020) 578:82–93. doi: 10.1038/s41586-020-1969-6

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Res (2013) 41:D991–995. doi: 10.1093/nar/gks1193

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Roessler S, Jia HL, Budhu A, Forgues M, Ye QH, Lee JS, et al. A unique metastasis gene signature enables prediction of tumor relapse in early-stage hepatocellular carcinoma patients. Cancer Res (2010) 70:10202–12. doi: 10.1158/0008-5472.CAN-10-2607

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Pinyol R, Montal R, Bassaganyas L, Sia D, Takayama T, Chau GY, et al. Molecular predictors of prevention of recurrence in HCC with sorafenib as adjuvant treatment and prognostic factors in the phase 3 STORM trial. Gut (2019) 68:1065–75. doi: 10.1136/gutjnl-2018-316408

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Lauss M, Donia M, Harbst K, Andersen R, Mitra S, Rosengren F, et al. Mutational and putative neoantigen load predict clinical benefit of adoptive T cell therapy in melanoma. Nat Commun (2017) 8:1738. doi: 10.1038/s41467-017-01460-0

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Prat A, Navarro A, Paré L, Reguart N, Galván P, Pascual T, et al. Immune-related gene expression profiling after PD-1 blockade in non-small cell lung carcinoma, head and neck squamous cell carcinoma, and melanoma. Cancer Res (2017) 77:3540–50. doi: 10.1158/0008-5472.CAN-16-3556

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Yu W, Ma Y, Hou W, Wang F, Cheng W, Qiu F, et al. Identification of immune-related lncRNA prognostic signature and molecular subtypes for glioblastoma. Front Immunol (2021) 12:706936. doi: 10.3389/fimmu.2021.706936

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Hein DM, Deng W, Bleile M, Kazmi SA, Rhead B, De La Vega FM, et al. Racial and ethnic differences in genomic profiling of early onset colorectal cancer. J Natl Cancer Inst (2022) 114:775–8. doi: 10.1093/jnci/djac014

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Niknafs N, Balan A, Cherry C, Hummelink K, Monkhorst K, Shao XM, et al. Persistent mutation burden drives sustained anti-tumor immune responses. Nat Med (2023) 29:440–9. doi: 10.1038/s41591-022-02163-w

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28:1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Wang T, Dai L, Shen S, Yang Y, Yang M, Yang X, et al. Comprehensive molecular analyses of a macrophage-related gene signature with regard to prognosis, immune features, and biomarkers for immunotherapy in hepatocellular carcinoma based on WGCNA and the LASSO algorithm. Front Immunol (2022) 13:843408. doi: 10.3389/fimmu.2022.843408

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Chen X, Lan H, He D, Xu R, Zhang Y, Cheng Y, et al. Multi-omics profiling identifies risk hypoxia-related signatures for ovarian cancer prognosis. Front Immunol (2021) 12:645839. doi: 10.3389/fimmu.2021.645839

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Charras G, Lenz M. A biochemical timer phases condensates in and out in cells. Nature (2022) 609:469–70. doi: 10.1038/d41586-022-01794-w

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Li CX, Wang JS, Wang WN, Xu DK, Zhou YT, Sun FZ, et al. Expression dynamics of periodic transcripts during cancer cell cycle progression and their correlation with anticancer drug sensitivity. Mil Med Res (2022) 9:71. doi: 10.1186/s40779-022-00432-w

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Wang Z, Wang Y, Yang T, Xing H, Wang Y, Gao L, et al. Machine learning revealed stemness features and a novel stemness-based classification with appealing implications in discriminating the prognosis, immunotherapy and temozolomide responses of 906 glioblastoma patients. Brief Bioinform (2021) 22. doi: 10.1093/bib/bbab032

CrossRef Full Text | Google Scholar

53. Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, et al. The connectivity map: using gene-expression signatures to connect small molecules, genes, and disease. Science (2006) 313:1929–35. doi: 10.1126/science.1132939

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Song F, Chen FY, Wu SY, Hu B, Liang XL, Yang HQ, et al. Mucin 1 promotes tumor progression through activating WNT/beta-catenin signaling pathway in intrahepatic cholangiocarcinoma. J Cancer (2021) 12:6937–47. doi: 10.7150/jca.63235

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Zhou X, Luo J, Xie H, Wei Z, Li T, Liu J. MCM2 promotes the stemness and sorafenib resistance of hepatocellular carcinoma cells via hippo signaling. Cell Death Discov (2022) 8:418. doi: 10.1038/s41420-022-01201-3

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Matsuno Y, Atsumi Y, Shimizu A, Katayama K, Fujimori H, Hyodo M, et al. Replication stress triggers microsatellite destabilization and hypermutation leading to clonal expansion in vitro. Nat Commun (2019) 10:3925. doi: 10.1038/s41467-019-11760-2

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Carlino MS, Larkin J, Long GV. Immune checkpoint inhibitors in melanoma. Lancet (2021) 398:1002–14. doi: 10.1016/S0140-6736(21)01206-X

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Liu SY, Wu YL. Ongoing clinical trials of PD-1 and PD-L1 inhibitors for lung cancer in China. J Hematol Oncol (2017) 10:136. doi: 10.1186/s13045-017-0506-z

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Awad MM, Hammerman PS. Durable responses with PD-1 inhibition in lung and kidney cancer and the ongoing search for predictive biomarkers. J Clin Oncol (2015) 33:1993–4. doi: 10.1200/JCO.2015.61.4172

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Gao Q, Zhu H, Dong L, Shi W, Chen R, Song Z, et al. Integrated proteogenomic characterization of HBV-related hepatocellular carcinoma. Cell (2019) 179:561–77. e522. doi: 10.1016/j.cell.2019.08.052

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Shiba-Ishii A, Kim Y, Shiozawa T, Iyama S, Satomi K, Kano J, et al. Stratifin accelerates progression of lung adenocarcinoma at an early stage. Mol Cancer (2015) 14:142. doi: 10.1186/s12943-015-0414-1

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Adachi H, Murakami Y, Tanaka H, Nakata S. Increase of stratifin triggered by ultraviolet irradiation is possibly related to premature aging of human skin. Exp Dermatol (2014) 23 Suppl:1, 32–36. doi: 10.1111/exd.12390

CrossRef Full Text | Google Scholar

63. Zhang Y, Ma K, Sadana P, Chowdhury F, Gaillard S, Wang F, et al. Estrogen-related receptors stimulate pyruvate dehydrogenase kinase isoform 4 gene expression. J Biol Chem (2006) 281:39897–906. doi: 10.1074/jbc.M608657200

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Sun Y, Daemen A, Hatzivassiliou G, Arnott D, Wilson C, Zhuang G, et al. Metabolic and transcriptional profiling reveals pyruvate dehydrogenase kinase 4 as a mediator of epithelial-mesenchymal transition and drug resistance in tumor cells. Cancer Metab (2014) 2:20. doi: 10.1186/2049-3002-2-20

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Blouin JM, Penot G, Collinet M, Nacfer M, Forest C, Laurent-Puig P, et al. Butyrate elicits a metabolic switch in human colon cancer cells by targeting the pyruvate dehydrogenase complex. Int J Cancer (2011) 128:2591–601. doi: 10.1002/ijc.25599

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Jha MK, Song GJ, Lee MG, Jeoung NH, Go Y, Harris RA, et al. Metabolic connection of inflammatory pain: pivotal role of a pyruvate dehydrogenase kinase-pyruvate dehydrogenase-lactic acid axis. J Neurosci (2015) 35:14353–69. doi: 10.1523/JNEUROSCI.1910-15.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Zou Y, Xie J, Zheng S, Liu W, Tang Y, Tian W, et al. Leveraging diverse cell-death patterns to predict the prognosis and drug sensitivity of triple-negative breast cancer patients after surgery. Int J Surg (2022) 107:106936. doi: 10.1016/j.ijsu.2022.106936

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Egelston CA, Avalos C, Tu TY, Simons DL, Jimenez G, Jung JY, et al. Human breast tumor-infiltrating CD8(+) T cells retain polyfunctionality despite PD-1 expression. Nat Commun (2018) 9:4297. doi: 10.1038/s41467-018-06653-9

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Bruix J, Cheng AL, Meinhardt G, Nakajima K, De Sanctis Y, Llovet J. Prognostic factors and predictors of sorafenib benefit in patients with hepatocellular carcinoma: analysis of two phase III studies. J Hepatol (2017) 67:999–1008. doi: 10.1016/j.jhep.2017.06.026

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Chen CL, Uthaya Kumar DB, Punj V, Xu J, Sher L, Tahara SM, et al. NANOG metabolically reprograms tumor-initiating stem-like cells through tumorigenic changes in oxidative phosphorylation and fatty acid metabolism. Cell Metab (2016) 23:206–19. doi: 10.1016/j.cmet.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Yu P, Wang Y, Yang WT, Li Z, Zhang XJ, Zhou L, et al. Upregulation of the PPAR signaling pathway and accumulation of lipids are related to the morphological and structural transformation of the dragon-eye goldfish eye. Sci China Life Sci (2021) 64:1031–49. doi: 10.1007/s11427-020-1814-1

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Sun M, Wang C, Lv M, Fan Z, Du J. Intracellular self-assembly of peptides to induce apoptosis against drug-resistant melanoma. J Am Chem Soc (2022) 144:7337–45. doi: 10.1021/jacs.2c00697

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Wang T, Zheng X, Li R, Liu X, Wu J, Zhong X, et al. Integrated bioinformatic analysis reveals YWHAB as a novel diagnostic biomarker for idiopathic pulmonary arterial hypertension. J Cell Physiol (2019) 234:6449–62. doi: 10.1002/jcp.27381

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Feng L, Dou C, Wang J, Jiang C, Ma X, Liu J. Upregulated 14−3−3beta aggravates restenosis by promoting cell migration following vascular injury in diabetic rats with elevated levels of free fatty acids. Int J Mol Med (2018) 42:1074–85. doi: 10.3892/ijmm.2018.3671

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Bodet L, Gomez-Bougie P, Touzeau C, Dousset C, Descamps G, Maïga S, et al. ABT-737 is highly effective against molecular subgroups of multiple myeloma. Blood (2011) 118:3901–10. doi: 10.1182/blood-2010-11-317438

PubMed Abstract | CrossRef Full Text | Google Scholar

76. da Silva Rosa SC, Martens MD, Field JT, Nguyen L, Kereliuk SM, Hai Y, et al. BNIP3L/Nix-induced mitochondrial fission, mitophagy, and impaired myocyte glucose uptake are abrogated by PRKA/PKA phosphorylation. Autophagy (2021) 17:2257–72. doi: 10.1080/15548627.2020.1821548

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: hepatocellular carcinoma (HCC), PANoptosis, immune characteristics, drug sensitivity, YWHAB, prognostic index

Citation: Song F, Wang C-G, Mao J-Z, Wang T-L, Liang X-L, Hu C-W, Zhang Y, Han L and Chen Z (2023) PANoptosis-based molecular subtyping and HPAN-index predicts therapeutic response and survival in hepatocellular carcinoma. Front. Immunol. 14:1197152. doi: 10.3389/fimmu.2023.1197152

Received: 30 March 2023; Accepted: 25 May 2023;
Published: 15 June 2023.

Edited by:

Jiajie Peng, Northwestern Polytechnical University, China

Reviewed by:

Yingcheng Wu, Fudan University, China
Shaohui Wang, Chengdu University of Traditional Chinese Medicine, China
SangJoon Lee, Ulsan National Institute of Science and Technology, Republic of Korea

Copyright © 2023 Song, Wang, Mao, Wang, Liang, Hu, Zhang, Han and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zhong Chen, chenz9806@163.com

These authors have contributed equally to this work and share first authorship

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.