Feng Cao
Yong Qi
Wenyong Wu
Xutong Li
Chuang Yang- 1Department of General, Visceral and Transplantation Surgery, University Hospital RWTH Aachen, Aachen, Germany
- 2Department of General Surgery, The First Hospital of Anhui Medical University, Hefei, China
- 3Department of Infectious Diseases, The First Affiliated Hospital of Anhui Medical University, Hefei, China
- 4Department of Visceral, Transplant, Thoracic and Vascular Surgery, University Hospital of Leipzig, Leipzig, Germany
Background: Cuproptosis, as a recently discovered type of programmed cell death, occupies a very important role in hepatocellular carcinoma (HCC) and provides new methods for immunotherapy; however, the functions of cuproptosis in HCC are still unclear.
Methods: We first analyzed the transcriptome data and clinical information of 526 HCC patients using multiple algorithms in R language and extensively described the copy number variation, prognostic and immune infiltration characteristics of cuproptosis related genes (CRGs). Then, the hub CRG related genes associated with prognosis through LASSO and Cox regression analyses and constructed a prognostic prediction model including multiple molecular markers and clinicopathological parameters through training cohorts, then this model was verified by test cohorts. On the basis of the model, the clinicopathological indicators, immune infiltration and tumor microenvironment characteristics of HCC patients were further explored via bioinformation analysis. Then, We further explored the key gene biological function by single-cell analysis, cell viability and transwell experiments. Meantime, we also explored the molecular docking of the hub genes.
Results: We have screened 5 hub genes associated with HCC prognosis and constructed a prognosis prediction scoring model. And the model results showed that patients in the high-risk group had poor prognosis and the expression levels of multiple immune markers, including PD-L1, CD276 and CTLA4, were higher than those patients in the low-risk group. We found a significant correlation between risk score and M0 macrophages and memory CD4+ T cells. And the single-cell analysis and molecular experiments showed that BEX1 were higher expressed in HCC tissues and deletion inhibited the proliferation, invasion and migration and EMT pathway of HCC cells. Finally, it was observed that BEX1 could bind to sorafenib to form a stable conformation.
Conclusion: The study not only revealed the multiomics characteristics of CRGs in HCC but also constructed a new high-accuracy prognostic prediction model. Meanwhile, BEX1 were also identified as hub genes that can mediate the cuproptosis of hepatocytes as potential therapeutic targets for HCC.
1 Introduction
The latest global cancer statistics indicate that there are more than 900,000 new cases of liver cancer and 800,000 deaths per year, and its morbidity and mortality rank seventh and third among all cancers, respectively (Sung et al., 2021), which pose a serious threat to human health and safety (Llovet et al., 2021). Among them, nearly 90% of the cases are hepatocellular carcinoma (HCC) and it is also the most common type of pathology (Llovet et al., 2016). HCC progresses rapidly and stealthily, with a high degree of malignancy, the prognosis of HCC is very poor and the 5-year survival rate is below 20% (Craig et al., 2020). Although alpha-fetoprotein (AFP) is currently widely used in early HCC screening and prognosis assessment, it is affected not only by many non-hepatic carcinoma-related factors, but its expression level is also significantly increased in other diseases, such as acute viral hepatitis (AHA), resulting in low sensitivity and specificity (Huo et al., 2018; Cai et al., 2019). For lack of effective biomarkers, most patients are in the advanced stages when they are diagnosed with HCC. Although a variety of modalities including surgical resection, radiofrequency ablation, trans-arterial chemoembolization, systemic chemotherapy and liver transplantation have been significantly developed, the clinical effect is still very limited. In addition, the complex pathogenesis of HCC is also an important factor in its high mortality, which involves multiple molecular mechanisms, such as cell death regulation, genetic mutation, tumor heterogeneity, immune regulation, and dysregulation of the tumor microenvironment (Kurebayashi et al., 2018; Demirtas and Gunduz, 2021; Torrens et al., 2021). Thus, it is vital to identify new HCC-specific biomarkers, extensively study the pathogenesis of HCC and explore new precise therapeutic targets.
In recent years, with the in-depth understanding of cell death, a growing number of studies have confirmed that a serious of programmed cell death (PCD) programs, such as autophagy, ferroptosis and pyroptosis, play an indispensable role in tumorigenesis and progression (Mou et al., 2019; Li et al., 2020; Hou et al., 2020). Cuproptosis, as a new PCD program, has the characteristics of copper dependence and copper regulation. The mechanism is mainly through copper binding directly to lipoylated components in the tricarboxylic acid (TCA) cycle to mediate lipoylated protein aggregation and iron-sulfur cluster protein loss, which ultimately leads to proteotoxic stress and cell death (Tsvetkov et al., 2022). Copper has the functions as a metabolic cofactor at the active site as well as a dynamic signaling metal and metalloallosteric regulator, which is connected with various clinical diseases, especially cancer, because the growth and metastasis of tumors have higher requirements for this metal nutrient (Ge et al., 2022). As the Mortada study found, the copper levels in the plasma and bladder tissue were significantly higher in patients with bladder cancer than in those non-bladder cancer patients (Mortada et al., 2020). Interestingly, Atakul et al. (2020) found that the serum copper level of patients with endometrial cancer was significantly lower than that of normal people, and it was negatively correlated with the degree of tumor invasion. Relevant studies have shown that unbalanced copper homeostasis can induce many forms of cell death including apoptosis, autophagy and ferroptosis, through various mechanisms, such as reactive oxygen species accumulation, proteasome inhibition and mitochondrial dysfunction (Jiang et al., 2022). However, the specific mechanism of copper ion-induced PCD is still unclear. In addition, with the role of copper in tumor proliferation, invasion and metastasis, its antitumor potential has also been highlighted. Li et al. (2020) showed that the combination of disulfiram with copper can greatly improve its antitumor efficacy. Similarly, Mariani et al. (2021) also demonstrated that the combined use of copper complexes and cisplatin enhanced the antitumor effect against melanoma, lung cancer and breast cancer. Therefore, this could be a new strategy for cancer treatment by using copper ion metal carriers to eliminate cancer cells. However, the role of copper-induced cuoproptosis in HCC has not been reported. Therefore, this study aimed to explore the possible molecular markers and drug targets of copper death in HCC, comprehensively analyze the multiomics characteristics of cuproptosis related genes (CRGs), including genomics, transcriptomics and tumor microenvironment (TME), and extensively investigate the latent function of cuproptosis in the TME, clinical characteristics and prognosis of HCC to provide a new strategy and prediction model for clinical diagnosis, treatment and prognosis evaluation.
2 Methods and materials
2.1 Data sources
The workflow chart of this study is shown in Supplementary Figure S1. We downloaded 424 HCC patients with ribonucleotide (RNA) sequences, clinical information and GSE76427, GSE52018, GSE149614 datasets from the Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) databases and Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). All raw files were normalized and annotated by using “limma” package in R, and the RNA sequences of fragments per kilobase million (FPKM) in TCGA were converted to transcripts per kilobase million (TPM) sequences. Then, batch effects of the three datasets were eliminated by the “Combat” algorithm. After integrating all datasets and excluding patients lacking overall survival (OS) data, the clinical data of 526 HCC patients were saved for succeeding analysis.
2.2 Difference analysis and cluster analysis of CRGs
We obtained 10 CRGs from previous publications [13]. According to the expression levels of CRGs in the HCC genomic dataset, the R language “reshape2” and “ggpubr” packages were used to analyze and plot the expression differences of CRGs in tumor and normal samples. Then, Kaplan‒Meier survival curves and interaction networks were analyzed and plotted using the HCC clinical dataset. Finally, according to the gene expression levels of CRGs, consensus unsupervised clustering (CUC) analysis was performed using the k-means method in the “ConensusClusterPlus” package to classify all HCC patients into different molecular subtypes. The principal component analysis (PCA) graph was drawn by analyzing the expression levels of CRGs and typing results by the “ggplot2” software package.
2.3 Difference analysis between CRG-related subtypes
We compared the prognosis and clinical characteristics of CRGS-related subtypes according to the clinical information files and classification results of HCC patients and drew Kaplan‒Meier survival curves and clinical characteristic heatmaps using R language for visualization. Then, we used the R language “GSVA” package to complete the gene set variation analysis (GSVA) enrichment analysis. At the same time, the ssGSEA algorithm in the “GSVA” package was used to quantitatively analyze the immune cells to compare the immune infiltration fraction between different subtypes.
2.4 Difference analysis between DEG-related subtypes and functional annotation
At first, we used R language “limma” package with parameters set to fold-change of 1.5 and adjusted p-value of <0.05 to extract CRG-related genes. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes enrichment (KEGG) functional enrichment analyses were performed. Then, univariate cox regression analysis was performed on CRG-related genes to screen for prognosis-related differentially expressed genes (DEGs). The prognostic and clinical features between DEG-related subtypes were compared using the R “survival” package and the “PheATmap” package and visualized using Kaplan‒Meier survival curves and heatmaps. In addition, the differences in the expression of CRGs among the related types of CRG-related genes were compared again using the R packages “Reshape2” and “GGPubR.”
2.5 Prognosis model construction and validation
First, 526 HCC patients were randomly divided into training and test cohorts with 263 cases each according to the 1:1 ratio. The R package “Glmnet” was used to include prognosis-related DEGs in the Least Absolute shrinkage and Selection Operator (LASSO) and multivariate Cox regression analyses. Then, risk scores (RSs) were calculated for genes with nonzero regression coefficients. The RS was calculated as follows: RS =
2.6 Establishment of a nomogram scoring system
Based on the clinical data and grouping results, the “pheatmap” package was used to visualize the correlation of RS with survival status and the difference in expression of prognosis-related DEGs. Then, the “rms” package was used to construct a nomogram scoring system for predicting prognosis, and the “calibrate” function was used to draw a calibration curve for evaluating the accuracy of the prognostic prediction model.
2.7 Analysis of immune cell infiltration, TME and CSCs
The 22 immune cell infiltration degrees of each sample were evaluated by the “preprocessCore” and “e1071” packages, and then the “CiberSort” algorithm was used to calculate the correlation between immune cell infiltration degree and RS- and prognosis-related DEGs. Next, we used the “estimate” package to calculate the stromal score, immune score and estimation score and further evaluated tumor purity. In addition, we evaluated the relationship between CSCs and RSs using cancer stem cell (CSC) score files.
2.8 Somatic mutation and drug susceptibility analysis
Differences in somatic mutations of HCC patients were analyzed using the “maftools” package and displayed in a waterfall diagram. We continued to use the “pRophetic” software package and the genomics of drug sensitivity in cancer (GDSC) database (https://cancerrxgene.org) to calculate the semi-inhibitory concentration (IC50) value of commonly used chemotherapy drugs in tumors, and the Wilcoxon signed rank test was used to compare the difference in IC50 values between the risk group and the low-risk group. The filter condition was set to p < 0.001.
2.9 Western blotting
We randomly collected tumor tissues and paracancerous tissues from 3 HCC patients treated at the First Affiliated Hospital of Anhui Medical University from March to April 2022. Tissue and normal LO2 hepatocytes and Hep3b, Huh7 and LM3 HCC cells were extracted with RIPA (Beyotime, Shanghai, China) after grinding and filtration. Then, the tissue was lysed on ice and centrifuged. We collected the supernatant and then the protein concentration was determined by a BCA kit (Beyotime). After gel preparation, electrophoresis, and transfer, the membrane was incubated with primary (BEX1, 1:5000, 12390-1-AP, Proteintech), (G6PC, 1:5000, 66860-1-Ig, Proteintech), (NEIL3, 1:5000, 11621-1-AP, Proteintech), (GCLM, 1:5000, 66808-1-Ig, Proteintech), (NT5DC2, 1:5000, YLS-K0806, Yilisa) and secondary antibodies (Anti-rabbit IgG, HRP-linked Antibody, 7074S, CST). For Western blotting analysis, the proteins underwent separation by SDS‒PAGE, nitrocellulose membrane transfer, quick blocker kit (Beyotime) blocking. After that, the bands were detected by ECL Plus (EMD Millipore, Billerica, MA, United States).
2.10 Single-cell analysis
GSE149614 raw data were processed in R using the “Seurat” package, and cells within the tissue were filtered and visualized with parameters of 500–6,000 expressed genes and mitochondrial ratio >10%. After performing quality control (QC) and selecting cells using the “CreateSeuratObject” algorithm. Then, we used the Seurat function of FindIntegrateionAnchors to merge sample files with common anchors among variables (dims = 1:20). We performed PCA on the scaled data using the “JackStraw” algorithm and ranked the PCA using the “ElbowPlot” function. The FindClusters function was used to perform T-distributed stochastic neighbor embedding (tSNE) dimensionality reduction clustering (resolution = 0.2) on the first 20 PC data. The FindMarkers function of Seurat and the CIBERSORT algorithm were used to determine the expression of marker genes in each cluster.
2.11 Cell viability, migration and invasion assays
LM3 cells were purchased from. Then cells were cultured with RPMI-1640 supplemented with 10% FBS (Beyotime, Shanghai, China). Subsequently, two siRNAs were used to knock out the expression of BEX1. Specifically, the sequences of siRNAs for BEX1 were listed. Then, LM3 cells were cultured until the density reached about 60%–70%, and plasmids, shRNA and Lipofectamine™ 2000 diluents were prepared by serum-free Opti-MEM medium to prepare transfection complexes. The transfection complex was added to LM3 cells and cultured for 48 h. Then 5,000 cells with BEX1 siRNAs were cultured in 96-well plate and incubated for 48 h at 37°C. Then the cell viability was measured using CCK-8 kit. In addition, 24-well inserts were used to perform cell migration and invasion assays. Briefly, LM3 cells were infected with BEX1 siRNAs for 48 h. Then 5,000 cells were seeded in the upper chamber with (migration assay) or without Matrigel (invasion assay). Meanwhile, 500 μl fresh medium were transferred to the low chamber and incubated for 24 h. After that, cells were fixed by 3.7% paraformaldehyde and stained with crystal violet. Cell images were obtained by fluorescence microscope and counted using ImageJ.
2.12 Molecular docking
We used AutoDock Vina software for molecular docking (Trott and Olson, 2010). Sorafenib was used as a ligand, and the key genes BEX1, NEIL3, GCLM, G6PC and NT5DC2 were used as receptors. The PDB format files were downloaded from the RSCB PDB database (http://www.rcsb.org/). Convert the Suolafeini PDB format to MOL2 format using Chem3D. Then, AutoDockTools 1.5.6 (https://autodock.scripps.edu/) was used to process receptor proteins and small molecule ligands and saved as PDBQT format files. During the molecular docking process, the Lamarckian algorithm was used to identify the most binding mode. The search space volume was > 27,000 A^3, the exhaustiveness was set to 8 and the maximum number of conformations output was set to 15.
2.13 Statistical analyses
The study was carried out under the R version 4.1.3, Strawberry Perl version 5.32.1.1, GraphPad Prism 7 and ImageJ. *, **, *** indicated p < 0.05, p < 0.01, p < 0.001, respectively. A p < 0.05 was considered statistically significant.
3 Results
3.1 Differential expression of CRGs and identification of cuproptosis-related subtypes
A total of 10 CRGs were obtained in this study for subsequent analysis. First, the results of copy number variation (CNV) analysis showed that there were clear CNVs in all genes except MIF1. CDKN2A had the highest gene mutation frequency and was mainly amplified on chromosome 9. The remaining CRGs also had different degrees of deletion and amplification variation, and their positions on the chromosome were also shown (Figures 1A, B). In addition, we obtained results showed that the remaining CRGs, except FDX1 and MIF1, were significantly highly expressed in tumor tissues. Interestingly, CDKN2A, which had the highest frequency of CNVs, also had the most significant differences in expression (Figure 1C). This suggests that the expression levels of genes in tumor tissues may be regulated by their CNVs. The results of the interaction network showed that CDKN2A and FDX1 were negatively correlated and were not found to be associated with other CRGs, while other genes interacted closely. Among them, FDX1, DLD and LIAS may be anticancer factors (Figure 1D). Further survival analysis found that CDKN2A, DLAT, FDX1, PDHA1, GLS and LIPT1 were associated with HCC prognosis (Figures 1E–J). Patients with highly-expressed FDX1 had a better prognosis, which further verified the anticancer effect of FDX1. The remaining five genes are the opposite. In addition, HCC patients were classified by CUC analysis of CRGs. The obtained results showed that k = 2 was the best choice, and the correlation between samples in subtypes A and B was the highest (Supplementary Figure S2). Besides, HCC samples can be separated into subtype A and subtype B, which was also verified by PCA (Supplementary Figure S3).
FIGURE 1. Genetic and transcriptional alterations of CRGs in HCC. (A) CNV of CRGs in 424 samples in TCGA. (B) Chromosomal localization of CRGs with CNV. (C) Expression of CRGs. (D) CRGs interaction network. (E–J) Kaplan-Meier survival analysis of 6 CRGs associated with HCC prognosis.
3.2 Differences in clinical features and functional annotations of CRG-related genes across subtypes
Kaplan‒Meier survival curves between different subtypes showed that patients with subtype A had a longer OS time than patients with subtype B (Figure 2A). Compared with subtype B, patients with subtype A had lower TNM stage and clinical stage than those with subtype B, the expression of FDX1 was relatively high, while the other genes were highly expressed in subtype B (Figure 2B). Moreover, GSVA suggested that subtype A was primarily enriched in fatty acid, bile acid and amino acid metabolism pathways, while cell cycle regulation, nucleic acid repair, nucleic acid synthesis and metabolism pathways were more likely to be found in subtype B (Figure 2C). The dysregulation of these pathways can cause gene mutation which can lead to abnormal cell metabolism or cell death. In the differences of immune cell infiltration. We found that activated B cell, activated CD8 T cell, eosinophilna, macrophagena, etc., degree was significantly higher in subtype A than in subtype B, whereas those of activated CD4 T cells and type 2 T helper cells were lower than those of subtype B (Figure 2D). To explore the underlying molecular behaviors of cuproptosis pattern, we identified 614 CRG-related DEGs and performed functional annotation analysis. These CRG-related DEGs were mainly enriched in biological processes such as cell cycle regulation and DNA replication (Supplementary Figure S4A). KEGG analysis showed that these genes were mainly enriched in cell cycling and DNA replication signaling pathways (Supplementary Figure S4B).
FIGURE 2. Difference analysis in distinct subtypes and functional annotations of CRGs related genes. (A) Survival curves between CRGs subtypes. (B) Clinical characteristics of CRGs subtypes. (C) GSVA enrichment analysis. (D) Differences in immune cell infiltration.
3.3 Identification of gene subtypes and construction of a prognostic model
First, we obtained 382 DEGs that were associated with prognosis. We divided HCC patients into two subtypes, A and B, by CUC analysis (Supplementary Figure S5). Kaplan‒Meier survival curves showed that patients of subtype A survived longer than those of subtype B (Figure 3A). Combined with clinical characteristics, a positive correlation was found between B subtype pattern and advanced TNM staging, BCLC staging, and high expression of prognosis-related genes (Figure 3B). In addition, FDX1, DLD, CDKN2A, DLAT, PDHB, PDHA1, GLS and LIPT1 also showed significant differences in expression between the two subtypes. Among them, only FDX1 was highly expressed in subtype A, which further verified the role of FDX1 as an anticancer factor (Figure 3C). We then constructed a predictive model based on CRG-related genes. LASSO regression analysis was performed on CRG-related DEGs with prognostic value, and the risk coefficient of each CRG-related DEGs was evaluated. Thirteen prognostic CRG-related DEGs were retained according to the minimum partial likelihood deviance (Figures 3D, E). Then, we performed multivariate Cox regression analysis on 13 genes related to prognosis and finally found that BEX1, G6PC, GCLM, NEIL3 and NT5DC2 CRG-related genes are independent influencing factors of HCC. The Sankey diagram shows the distribution and correlation of patients with different prognoses in distinct subtypes and RS subgroups (Figure 3F). Among the different gene subtypes, the RS of subtype B was higher than that of subtype A (Figure 3G). Interestingly, the difference in CRG expression between the two groups was consistent with the difference between gene subtypes. Only FDX1 was highly expressed in the low-risk group, while DLD, CDKN2A, DLAT, PDHB, PDHA1, GLS and LIPT1 were all highly expressed in the high-risk group (Figure 3H). This further demonstrates the anticancer effect of FDX1. In addition, there was also a clear correlation between RS and immune checkpoints, including programmed cell death protein 1 (PD-1), programmed cell death 1 ligand 1 (PD-L1), cytotoxic T-lymphocyte-associated protein 4 (CTLA4) and cluster of differentiation 44 (CD44), which were highly expressed in the high-risk group (Figure 3I).
FIGURE 3. Identification of hub genes and construction of the prognostic model. (A) Survival curves between the hub genes subtypes. (B) Clinical characteristics of DEGs. (C) Expression of CRGs between DEGs subtypes. (D,E) LASSO regression analysis and partial likelihood deviance. (F) Alluvial diagram of subtype and RS distributions. (G) Differences of RS in DEGs clusters. (H,I) Expression of CRGs and immune checkpoints.
3.4 Validation of the prognostic model and development of a nomogram
Based on the median RS, the patients in the training and test groups were sorted according to the RS (Figures 4A, B). Survival status showed that patients in high-risk group of training group had a worse outcome than those in the low-risk group. We also found the similar results in the test group (Figures 4C, D). The heatmaps showed that NEIL3, GCLM and NT5DC2 had a higher expression level in the high-risk group, and G6PC and BEX1 were expressed at low levels (Figures 4E, F). Furthermore, in both the training and test groups, Kaplan‒Meier survival curves showed that the survival rate in the high-risk group was significantly lower than that in the low-risk group (Figures 5A, B). This may be related to the high expression of BEX1 and G6PC in the low-risk group, indicating that BEX1 and G6PC may have a synergistic effect with FDX1. In addition, to assess the efficacy of the prognostic model, the areas under the curve (AUCs) of the training group at 1, 3, and 5 years were 0.792, 0.756, and 0.740, respectively. Likewise, the test group also had the similar performance, with AUCs of 0.720, 0.688, and 0.641 at 1, 3, and 5 years, respectively (Figures 5C, D). Therefore, this study created a nomogram that integrated multiple molecular markers and clinicopathological parameters to predict the prognosis (Figure 5E). Meanwhile, the calibration diagram also demonstrated that the model had a perfect reliability (Figure 5F).
FIGURE 4. Prognostic value of the CRGs signature. (A,B) RS distribution. (C,D) Survival status. (E,F) Expression of the 5 hub genes.
FIGURE 5. Prediction model and nomogram. (A,B) Kaplan-Meier curves in the training and test cohorts. (C,D) ROC curves estimate prognosis value. (E) Nomogram for predicting the OS of HCC patients. (F) Calibration curves of the nomogram.
3.5 Correlation of RS with immune cells, TME score and somatic mutation and drug susceptibility analysis
The protein interaction network between CRG-related DEGs showed that five key prognosis-related genes in the prognostic model were negatively correlated with FDX1. Among them, NEIL3 and NT5DC2 were positively correlated with CDKN2A, and the correlation was the highest (Figure 6A). In addition, we found RS was positively correlated with M0 macrophages and negatively correlated with resting memory CD4+ T cells (Figures 6B, C). Among the five key prognosis-related genes, G6PC, GCLM, NEIL3 and NT5DC2 were significantly associated with a variety of immune cells, among which G6PC was positively associated with M1 macrophages, and NEIL3 was positively correlated with activated memory CD4+ T cells and negatively correlated with resting memory CD4+ T cells (Figure 6D). TME difference analysis found that the stromal score in the low-risk group was significantly higher than that in the high-risk group (Figure 6E). Our study also assessed the correlation between CSC and RS and concluded that RS was positively correlated with CSC (Figure 6F). Furthermore, we analyzed the somatic mutation profile. The obtained results showed that the top 10 mutated genes in the two groups were TP53, CTNNB1, TTN, MUC16, ALB, PCLO, APOB, RYR2, MUC4 and FLG. Among them, the PCLO mutation frequency was higher in the low-risk group, while other genes were lower than those in the high-risk group (Figures 6G, H). Finally, we evaluated the association between RS and drug sensitivity. We found that the IC50 values of various chemotherapeutics were significantly different between the high- and low-risk groups, including axitinib, gefitinib and erlotinib, sorafenib, vinorelbine, gemcitabine, nilotinib and tipifarnib (Supplementary Figures S6A–H).
FIGURE 6. TME characteristics and drug susceptibility. (A)The interaction network of CRGs and hub genes. (B,C) Correlation between immune cells and RS. (D) Correlations between the immune cells and 5 hub genes. (E) Differences in the StromalScore and ImmuneScore. (F) Correlation of RS with CSCs. (G,H) Mutation of genes in distinct RS group.
3.6 BEX1 are differentially expressed in hepatocytes
Detecting protein expression in HCC cells found that these five key genes were abnormally expressed to varying degrees in different HCC cells, especially BEX1 (Figure 7A). Compared with normal tissues, the expression levels of BEX1 in tumor tissues were significantly different, while the GCLM, G6PC, NEIL3 and NT5DC2 protein bands were unclear, which may be related to their low expression levels (Figure 7B). This was confirmed by single-cell analysis, which showed that there was no significant expression in hepatocytes (Figures 7K, L). To further explore the cellular expression of BEX1 in HCC tumors, 18 HCC samples in the GSE149614 dataset were analyzed. First, we controlled the effects of low-quality cells, mitochondrial genes, ribosomal genes, and hemoglobin. The correlations between the total number of unique molecular identifiers (UMIs) in each cell and the mitochondrial ratio, the total number of genes, and the hemoglobin ratio were 0.11, 0.91, and −0.01, respectively (Supplementary Figures S7A, B). No significant separation trend of HCC cells was observed when PCA was used to reduce the dimension, and we finally selected the top 20 PCs for further analysis based on the elbow plots (Supplementary Figures S5C, D). tSNE analysis classified HCC cells into 30 clusters (Figure 7C). The cells were annotated and divided into immune cell clusters and non-immune cell clusters (Figure 7D). We then proceeded to annotate eight subclusters using single-cell markers, including B cells, endothelial cells, hepatocytes, macrophages, monocytes, NK cells, smooth muscle cells and T cells (Figure 7E). In addition, normal and tumor tissues were similarly clustered and annotated by tSNE (Figures 7F, G). Finally, we tested the expression patterns of target genes in HCC cell clusters. The obtained results showed that the expression levels of BEX1 were the same as Western blotting, and they were only expressed in hepatocytes (Figures 7H–J). This finding indicates that BEX1 may be the key genes mediating cuproptosis in HCC cells.
FIGURE 7. Western blotting and single cell analysis of the 5 hub genes. (A,B) Western blotting in HCC cell and tumor tissue. (C–G) Cell clusters and annotates for GSE149614 of HCC patients. (H–L) Expression pattern of 5 hub genes at the single-cell level in normal and tumor cell clusters through t-SNE analysis.
3.7 Knockdown of BEX1 inhibited HCC cell proliferation, invasion and migration
We found that BEX1 had the most significant differences in the five hub genes after detection of liver cancer cells, tissues and single-cell analysis. In order to further explore the biological function of BEX1, we used two methods to knockdown BEX1, both of which could effectively silence the expression of BEX1 (Figure 8C). The results of cell viability experiment showed that the viabilities of liver cancer cells decreased significantly after 48 h of BEX1 knockdown compared with the control group (Figure 8A). In addition, transwell results showed that the invasion and migration of LM3 cells with low BEX1 were significantly reduced (Figure 8B). Finally, we also detected the expression of epithelial mesenchymal transition (EMT) markers including E-cadherin, N-cadherin and Vimentin. The results indicated that BEX1 knockdown promoted E-cadherin expression and decreased the expression of N-cadherin and Vimentin. Which demonstrated that inhibition of BEX1 suppressed liver cancer cells EMT pathway (Figure 8C). Those results indicated that BEX1 deletion inhibited the proliferation, invasion and migration and EMT pathway of HCC cells.
FIGURE 8. BEX1 biological function and molecular docking. (A) Cell viability experiment of BEX1 knockdown. (B,C) Transwell results and cadherin proteins expression in HCC cell with sh-BEX1. (D) Molecular docking of sorafenib and the 5 hub genes.
3.8 Sorafenib specifically binds to key genes
To further explore whether the target gene may be a potential target of HCC immunotherapy, we used sorafenib as the ligand and the hub gene as the receptor for molecular docking analysis. The obtained results showed that the binding energies of sorafenib to G6PC, GCLM, BEX1, NT5DC2 and NEIL3 were −9.6, −8.5, −6.7, −8.3, and −8.0 kcal/mol, respectively. These hub genes can bind to sorafenib and form a stable conformation (Figure 8D). This finding suggests that the high expression of BEX1 in hepatocytes may be a potential therapeutic target to mediate cuproptosis in HCC.
4 Discussion
Numerous studies have indicated that PCD plays a key role in tumorigenesis and progression of HCC and antitumor response, with the powerful immunotherapeutic potential (Demirtas and Gunduz, 2021). However, the specific function in anti-tumor remains unknown, and research on cuproptosis is even scarcer. Therefore, this study extensively explored the copy number variation, clinical feature relationship, TME and immune infiltration of CRGs in HCC at the transcriptome level, single-cell level and network pharmacology level and revealed that GCLM and BEX1 may be potential therapeutic targets mediating cuproptosis in HCC.
CNV is a characteristic change in neoplastic diseases that has been gradually recognized in the postgenomic era. The human genome contains numerous repetitive sequences of varying frequencies and intensities, and CNV is defined as a type of alteration involving deletion, insertion, replication and multilocus variation of gene segments ranging from 1 Kb to 3 Mb (Kinross et al., 2012; Shi et al., 2021). Current studies suggest that CNV is not only the basis of individual genetic differences but also plays an essential role in tumorigenesis, invasion and metastasis. CNV-related indicators may become ideal tumor diagnostic markers (Behroozi et al., 2020; Pariyar et al., 2021). Huang et al. (2021) found that estrogen-related receptor alpha (ESRRA) CNV was significantly correlated with the histological grade of ovarian cancer (OC). The results indicated that CNV has the function in affecting biological phenotype and heterogeneity of tumors and promoting tumor progression. However, the above mentioned study evaluated only the CNV of genes alone and lacked a comprehensive study combining CNV and gene expression. In this study, we analyzed them together, and the characteristics of cuproptosis genes in HCC were more accurately identified. We found clear CNVs in all CRGs except MIF1. Among them, CDKN2A had the highest mutation frequency and expression difference, while FDX1 had no difference in expression between normal tissue and tumor tissue, although there was a higher amplification mutation. Similarly, previous studies by Ghaffari et al. (2016) showed that the CNV of the baculoviral inhibitor of apoptosis repeat-containing 5 gene is highly increased in tumor tissue and may become a marker for early cancer detection and prognosis. Therefore, this study suggests that CNV detection of these highly expressed CRGs may contribute to the diagnosis of HCC. Of note, CNVs can also activate proto-oncogenes and reduce the activity of tumor suppressor genes, thereby mediating the pathogenesis and prognostic mechanisms of various tumors, including liver cancer (Diskin et al., 2009; Nik-Zainal et al., 2016). In the analysis of prognosis and clinical features in this study, it was found that CNV genes were closely related to OS, disease progression and immune invasion in HCC patients. This finding indicates that CNV genes, such as CDKN2A, plays a key role in the escape from immune surveillance.
Furthermore, this study also screened CRG-related genes to construct an HCC risk model for predicting prognosis, searching for immunotherapy strategies and drug targets. The obtained results indicated that five CRG-related genes were independent prognostic factors of HCC. We further constructed a new prognostic model based on these genes. Compared with traditional TNM staging, this model integrates multiple molecular markers, clinicopathological parameters and other multilevel prognostic indicators, which can not only identify HCC patients with different risks but also more accurately and dynamically monitor tumor progression and prognosis. Similar studies have been conducted in previous studies, such as by Song et al. (2021), who screened pyroptosis-related genes and constructed a colorectal cancer (CRC) risk model. Based on these model results, Liao et al. (2022) found that FOXP2 promotes CRC pyroptosis by interacting with caspase-1. In addition, several studies have reported the role of cuproptosis in HCC. Zhao et al. (2022) explored cuproptosis related genes (CRGs) related to HCC survival and clinical features. And they found 20 CRGs were correlated to HCC outcomes and might be used as a prognostic biomarker for HCC. Besides, Zhang et al. (2022) investigated the association between FDX1 expression and cancer stages and outcomes in HCC. They constructed the model based on FDX1 and its related genes. And the AUC values of cuproptosis-related risk score (CRRS) in predicting the OS were 0.72 and 0.68 at 1 and 3 years respectively. Similarly, Yan et al. (2022) developed a predictive model (GCSH, LIPT1 and CDKN2A) based on CRGs in HCC. The AUC values of ROC analysis for 1 year OS were from 0.614 to 0.683. And they found LIPT1 might be a target in the treatment of HCC. In our study, a predictive model was constructed consists of five genes (BEX1, NEIL3, GCLM, G6PC and NT5DC2) through comprehensive bioinformatics analysis. The model was validated by an internal dataset with an AUC was 0.720–0.792 for 1 year, 0.688–0.756 for 3 years and 0.641–0.740 for 5 years, respectively, which showed a robust performance than previous models. Besides, molecular docking and experiments indicated that BEX1 may mediate the cuproptosis of hepatocytes as potential therapeutic targets for HCC. Therefore, the five key genes and risk models screened in this study play an extremely vital role in the diagnosis and treatment of HCC. On the other hand, the accuracy of quantitative analysis of these genes is not less than that of whole transcriptome sequencing, and it is more economical and clinically feasible.
With the great breakthrough in immunotherapy, an increasing number of researchers have realized that tumor cells do not exist in isolation, and the TME in which they are located plays an indispensable role in tumor progression. The TME is mainly composed of immune cells including lymphocytes, macrophages and granulocytes in the center and fibroblasts, inflammatory cells and various signaling molecules in the surroundings (Turley et al., 2015; Seager et al., 2017). Numerous evidence have shown that the cellular components in the TME are strongly linked to the progression, metastasis and efficacy of HCC (Chew et al., 2017; Satilmis et al., 2021). For example, Mano et al. (2019) found that bone morphogenetic protein 4 (BNP-4) could enhance the aggressiveness of HCC by activating fibroblasts (CAFs) to secrete cytokines in the TME. Besides, tumor-associated macrophages (TAMs) also exhibited different activation states in the TME depending on the different stimulus. Previous studies have shown that TAMs in HCC tumor stroma produce various proinflammatory cytokines, including TNF-α, IL-β, IL-6 and IL-23, which induce the expansion of CD4+ Th17 cells, according to the overexpression of PD-L1, CTLA4 to inhibit antitumor immunity (Kuang et al., 2010). Previous studies have mainly concentrated on the innate T cell immune response. However, in the ongoing exploration of advanced HCC and other cancers, immunotherapy methods increasingly focus on immune checkpoint inhibitors (ICIs), such as CTLA-4, PD-1 and PD-L1, which lock the immune checkpoint inhibition pathways (Kudo, 2019a). These results indicate that ICIs based on TME changes will be a new therapeutic strategy for HCC.
However, in this research, we found that multiple immune markers, such as PD-L1, CD276, CD80 and CTLA4, were significantly higher in the high-risk group with worse prognosis than in the low-risk group, and RS was significantly correlated with M0 macrophages and memory CD4+ T cells. This suggests that HCC hepatocytes may achieve immune resistance or immune escape by overexpressing PD-L1 and binding to PD-1 on the surface of specific cytotoxic T cells. Of course, this mechanism is more complicated in practice. For example, Kudo (2019b) found that activated CD8+ T cells released interferon gamma (IFN-γ) during the process of initial cellular immunity, which not only attacked tumor cells but also engaged with receptors on the surface of cancer cells to upregulate PD-L1 and inhibit antitumor effects. Currently, research on PD-1 and PD-L1 inhibitors is in full swing and has become an important part of the systemic treatment of HCC in clinical practice. Numerous studies demonstrated that nivolumab and pembrolizumab have shown efficacy in the treatment of HCC especially as an alternative strategy after sorafenib failure or unacceptable toxicity (El-Khoueiry et al., 2017; Zhu et al., 2018). In addition, as a transmembrane receptor on the surface of activated T cells, the expression of CTLA4 is strictly regulated. In resting or naive T cells, CTLA4 is mainly located in intracellular vesicles (Valk et al., 2008). However, on activated T cells, CTLA4 directly competes with CD28 for B7 ligands to mediate tumor immune responses (Hathcock et al., 1993). In addition, CTLA4 activation also supported the transformation of CD4+ T cells into regulatory T cells by increasing transforming growth factor-β (TGFβ) secretion and forkhead box protein 3 (FOXP3) expression (Zheng et al., 2006). In fact, many mechanisms are involved in the tumor immune response in the TME, but the specific effects are still controversial.
In addition, we detected the protein level of five important genes and found that the expression of BEX1 in tumor cells and tissues was higher than that in normal liver cells and tissues. To further explore the biological function of BEX1, cell viability experiment and transwell results showed that BEX1 deletion inhibited the proliferation, invasion and migration of HCC cells. These results suggest that BEX1 may play a key role in HCC tumorigenesis and development. However, only few studies described the role of BEX1 in HCC. BEX1 (brain-expressed X-linked protein 1), attached to the BEX family and consists of five proteins with unclear functions (Kazi et al., 2015). BEX1 was initially thought to be associated with retinoic acid differentiation in teratomas (Faria et al., 1998). Later, multiple research teams successively reported its expression changes in many cancers (Quentmeier et al., 2005; Foltz et al., 2006; Kazi et al., 2015). Current research suggests that in addition to being involved in the regeneration of neuronal axons and regulating the cell cycle, BEX1 is also involved in the proliferation and invasion (Vilar et al., 2006; Khazaei et al., 2010; Doi et al., 2020; Lee et al., 2021) showed that BEX1 and BEX4 can improve the tumor formation and radio resistance of glioblastoma multiforme cells. In breast cancer, overexpression of BEX1 and BEX2 can inhibit the apoptosis of tumor cells (Naderi et al., 2007). Deficiency of Bex1 expression led to the decrease of cell proliferation, colony and tumor formation, and the increase of cell apoptosis in acute myeloid leukemia (Lindblad et al., 2015). These studies all suggest that BEX1 may be an oncogene, but there is still a lack of studies on the gene, especially in HCC. Sagawa et al. (2015) showed that BEX1 was upregulated in Cx32ΔTg rat liver, and knockdown of BEX1 could significantly inhibit the growth of rat hepatoma cell lines. In human HCC studies, Wang et al. (2021) discovered that BEX1 is an oncofetal protein that interacts with RUNX family transcription factor 3 (RUNX3) in hepatoblastoma (HB) and CSC-HCC to block β-catenin transcription and activate the Wnt/β-catenin signaling pathway, thereby regulating the self-renewal of hepatic CSCs. These results indicate that BEX1 is a promising therapeutic target forHB and CSC-HCC, and targeting the BeX1-mediated Wnt/β-catenin signaling pathway may help to resolve the heterogeneity and high recurrence rate of HCC.
Previous studies also explored the potential targets of cuproptosis in HCC. For example, FDX1, LIPT1, DLAT, CDKN2A and GLS were reported to have potential value as cuproptosis targets in HCC. Zhang et al. (2022) found FDX1 was downregulated in HCC, and higher expression of FDX1 was related to a better outcome. Beside, FDX1 inhibited the proliferation and colony formation of tumor cells in the presence of copper ions and this inhibitory effect was diminished by using tetrathiomolybdate (TTM) to chelate copper ions, which indicated that FDX1 exerted anti-tumor effects through cuproptosis (Li et al., 2022). And Yan et al. (2022) revealed that LIPT1 was higher expressed in HCC and patients with low LIPT1 expression had longer OS than those high LIPT1 expression. Cell proliferation assay showed that LIPT1 depletion inhibited cancer viability in HCC cells. Besides, LIPT1 knock down significantly suppressed cell migration and invasion capacity in HCC cells. Similar results were found in DLAT (Yang et al., 2023), CDKN2A and GLS (Ma et al., 2023). In our study, we found BEX1 deletion inhibited the proliferation, invasion and migration and EMT pathway of HCC cells. Besides, BEX1 can bind to sorafenib and form a stable conformation which further suggests that the BEX1 in hepatocytes may be a potential therapeutic target to mediate cuproptosis in HCC.
Furthermore, we detected the protein level of five important genes and found that the expression of the other genes except BEX1 in tumor cells and tissues was higher than that in normal liver cells and tissues. To further explore the specific cellular expression patterns of these genes in HCC, we conducted scRNA-seq data to identify the distribution of G6PC, GCLM NEIL3 and NT5DC2 in different cell subclusters in HCC tissues including HCC tumor cells and non-tumor cells. And the results clearly revealed that GCLM and G6PC were only highly expressed in HCC cells, while lowly expressed in non-tumor cells. The expression of NEIL3 and NT5DC2 was not obvious in hepatocellular carcinoma cells. These results suggest that GCLM and G6PC regulate HCC development directly on the target of hepatocytes rather than other cells, which provides potential therapeutic targets.
Glucose-6-phosphatase (G6PC) is a component critical for catalyzing glycogenolysis, and the downregulated of G6PC enhances glucose storage in premalignant cells (Franco et al., 2005; Resaz et al., 2014), and glycogen accumulation is a key carcinogenic event in the malignant transformation of the liver (Liu et al., 2021). Bioinformatics analysis demonstrated that low expression of G6PC was related to poor outcomes in HCC (Tian and Liao, 2022). However, there still is an empty in studies related to the regulation of copper by BEX1 and GP6C.
GCLM is a glutamate cysteine ligase modifier subunit that is the main component of glutathione (GSH) synthetase and participates in the synthesis and metabolism of GSH on multiple levels, and increasing evidence has shown that GSH metabolic dysregulation is involved in the pathophysiological mechanisms of various diseases, including diabetes, liver fibrosis, alcoholic liver disease, and malignant tumors (Liu and Gaston, 2010; Lv et al., 2019; Shen and Wang, 2021). Previous studies have shown that HCC patients showed higher levels of oxidative stress markers and low levels of GSH and GSH-related antioxidant enzymes in plasma compared with nonalcoholic steatohepatitis patients (Shimomura et al., 2017). The imbalance between high oxidative stress and low antioxidant capacity may be an important reason for the occurrence and development of HCC. Therefore, regulating the expression of GCLM may indirectly play an essential function in the development of HCC. In addition, studies have shown that GCLM expression is mainly regulated by transcription factors, including activator protein-1 (AP-1), nuclear factor kappa B (NFκB) and nuclear factor erythroid 2 related factor 2 (Nrf2) (Jaiswal, 2004; Yang et al., 2005; Roos et al., 2020) confirmed that treatment of human hepatocellular carcinoma (HepG2) cells with the receptor tyrosine kinase inhibitor lapatinib activates the Kelch-like ECH-associated protein 1 (Keap1)-Nrf2 signaling pathway, thereby upregulating GCLM levels and inducing GSH synthesis. Most interestingly, Thai et al. (2021) found that copper nanoparticles (Cu NPs) have the most significant effect on NrF2-mediated cytotoxicity, and upregulated GCLM can be used as a biomarker for Cu NP exposure in HCC cells. This is similar to our study, indicating that GCLM is more sensitive to copper damage and may be a key target in mediating cuproptosis in hepatocytes.
Finally, we verified the potential role of these hub genes as therapeutic targets in HCC. We performed drug-target network and molecular docking analyses with sorafenib as the ligand and the hub gene as the receptor, and the obtained results showed that BEX1, G6PC and GCLM could bind to sorafenib and form stable conformations. Sorafenib, as a multiple-target tyrosine kinase inhibitor (TKI), has the functions of anti-angiogenesis and anti-proliferation, besides it can also prolong the overall median survival time of patients with advanced HCC (Llovet et al., 2008). Two important clinical trials, the Asia-Pacific and Sorafenib HCC Assessment Randomized Protocol (SHARP), also showed that sorafenib has a powerful function in improving the prognosis of patients with HCC (Cheng et al., 2009; Vogel and Saborowski, 2020). In addition, Sorafenib is an effective first-line treatment in patients with advanced HCC (Xing et al., 2021). In vitro experiments revealed that sorafenib inhibits tumor cell viability, promotes cell apoptosis in HCC cells (Xie et al., 2018). Therefore, we used sorafenib as the ligand and the hub gene as the receptor for molecular docking analysis.
Although current ICIs have opened up a new strategy in treating malignant tumors, sorafenib is still the first-line drug for HCC chemotherapy (Benson et al., 2021). In recent years, studies on the discovery of new efficacy and potential therapeutic targets for known drugs are not uncommon; Pushpakom et al. (2019) found that Raloxifene was initially used to treat osteoporosis and was later approved for the treatment of breast cancer. Additionally, a study by Pang et al. (2019) is similar to ours; they confirmed that C3 and ANXN1 can stably bind to Vorinostat, which can be used as potential therapeutic targets for papillary renal cell carcinoma through genetic screening and molecular docking. Therefore, it is reasonable to speculate that BEX1, G6PC and GCLM can be used as potential therapeutic targets of HCC to mediate the mechanism of cuproptosis in hepatocytes.
Nevertheless, there are still several limitations in this study. First, this is a prospective study, outcomes can be affected by the introduction of unknown variables during follow-up, or changes in known variables including treatment modality disease and surgery-related complications. Second, our prognostic model was obtained based on TCGA and GEO data but lacked validation with external datasets, limited the breadth of its use. Third, we investigated the cellular molecular function of BEX1 in HCC but did not explore the interaction between BEX1 and copper ions in hepatic tumor cells, therefore, we could not clarify whether the BEX1-regulated cellular activity was related to copper ions. Finally, the specific mechanism by which BEX1 regulates HCC progression remains unclear. In the future studies, external datasets or one’s own data to assess the stability of the model is needed. Furthermore, in vivo and in vitro experiments are needed to investigate the exact mechanism by which BEX1 regulates cuproptosis especially in the presence of copper ions or when activity is inhibited.
5 Conclusion
This study extensively explored the multi-omics features of CRGs in HCC, including CNV, clinicopathological indicators, prognosis, immune infiltration, and TME, and constructed a prognostic prediction model integrating multiple molecular markers and clinicopathological parameters, which offers a new method for clinical diagnosis and prognosis evaluation. Meanwhile, GCLM and BEX1 were identified as hub genes, which are potential therapeutic targets to mediate the cuproptosis program in HCC cells.
Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.
Ethics statement
The study was carried out under the permission of the Ethics Committee of the First Affiliated Hospital of Anhui Medical University. Patients or their relatives signed written informed consent. All the studies were performed in accordance with the Declaration of Helsinki.
Author contributions
FC: Data curation, Writing–Original draft preparation. YQ: Software, Writing–Reviewing and Editing. WW: Visualization, Investigation. XL and CY: Methodology, Software and Revision. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the China Scholarship Council under No. 202008340087 and No. 202108340033. We acknowledge support from the German Research Foundation (DFG) and Leipzig University within the program of Open Access Publishing.
Acknowledgments
We thank the investigators who participated and provided data unselfishly in TCGA and GEO databases.
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/fcell.2023.1240390/full#supplementary-material
References
Atakul, T., Altinkaya, S. O., Abas, B. I., and Yenisey, C. (2020). Serum copper and zinc levels in patients with endometrial cancer. Biol. TRACE Elem. Res. 195 (1), 46–54. doi:10.1007/s12011-019-01844-x
Behroozi, J., Shahbazi, S., Bakhtiarizadeh, M. R., and Mahmoodzadeh, H. (2020). ADAR expression and copy number variation in patients with advanced gastric cancer. BMC Gastroenterol. 20 (1), 152. doi:10.1186/s12876-020-01299-8
Benson, A. B., D'Angelica, M. I., Abbott, D. E., Anaya, D. A., Anders, R., Are, C., et al. (2021). Hepatobiliary cancers, version 2.2021, NCCN clinical practice guidelines in oncology. J. Natl. Compr. CANC NE 19 (5), 541–565. doi:10.6004/jnccn.2021.0022
Cai, J., Chen, L., Zhang, Z., Zhang, X., Lu, X., Liu, W., et al. (2019). Genome-wide mapping of 5-hydroxymethylcytosines in circulating cell-free DNA as a non-invasive approach for early detection of hepatocellular carcinoma. Multicent. Study 68 (12), 2195–2205. doi:10.1136/gutjnl-2019-318882
Cheng, A. L., Kang, Y. K., Chen, Z., Tsao, C. J., Qin, S., Kim, J. S., et al. (2009). Efficacy and safety of sorafenib in patients in the asia-pacific region with advanced hepatocellular carcinoma: A phase III randomised, double-blind, placebo-controlled trial. LANCET Oncol. 10 (1), 25–34. doi:10.1016/S1470-2045(08)70285-7
Chew, V., Lai, L., Pan, L., Lim, C. J., Li, J., Ong, R., et al. (2017). Delineation of an immunosuppressive gradient in hepatocellular carcinoma using high-dimensional proteomic and transcriptomic analyses. P Natl. Acad. Sci. U. S. A. 114 (29), E5900-E5909–9. doi:10.1073/pnas.1706559114
Craig, A. J., von Felden, J., Garcia-Lezana, T., Sarcognato, S., and Villanueva, A. (2020). Tumour evolution in hepatocellular carcinoma. Nat. Rev. GASTRO Hepat. 17 (3), 139–152. doi:10.1038/s41575-019-0229-4
Demirtas, C. O., and Gunduz, F. (2021). Programmed cell death 1 and hepatocellular carcinoma: an epochal story. J. Gastrointest. CANC 52 (4), 1217–1222. doi:10.1007/s12029-021-00758-z
Diskin, S. J., Hou, C., Glessner, J. T., Attiyeh, E. F., Laudenslager, M., Bosse, K., et al. (2009). Copy number variation at 1q21.1 associated with neuroblastoma. NATURE 459 (7249), 987–991. doi:10.1038/nature08035
Doi, T., Ogawa, H., Tanaka, Y., Hayashi, Y., and Maniwa, Y. (2020). Bex1 significantly contributes to the proliferation and invasiveness of malignant tumor cells. Oncol. Lett. 20 (6), 362. doi:10.3892/ol.2020.12226
El-Khoueiry, A. B., Sangro, B., Yau, T., Crocenzi, T. S., Kudo, M., Hsu, C., et al. (2017). Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. LANCET 389 (10088), 2492–2502. doi:10.1016/S0140-6736(17)31046-2
Faria, T. N., LaRosa, G. J., Wilen, E., Liao, J., and Gudas, L. J. (1998). Characterization of genes which exhibit reduced expression during the retinoic acid-induced differentiation of F9 teratocarcinoma cells: involvement of cyclin D3 in RA-mediated growth arrest. Mol. Cell Endocrinol. 143 (1-2), 155–166. doi:10.1016/s0303-7207(98)00127-0
Foltz, G., Ryu, G. Y., Yoon, J. G., Nelson, T., Fahey, J., Frakes, A., et al. (2006). Genome-wide analysis of epigenetic silencing identifies BEX1 and BEX2 as candidate tumor suppressor genes in malignant glioma. CANCER Res. 66 (13), 6665–6674. doi:10.1158/0008-5472.CAN-05-4453
Franco, L. M., Krishnamurthy, V., Bali, D., Weinstein, D. A., Arn, P., Clary, B., et al. (2005). Hepatocellular carcinoma in glycogen storage disease type ia: A case series. J. Inherit. Metab. Dis. 28 (2), 153–162. doi:10.1007/s10545-005-7500-2
Ge, E. J., Bush, A. I., Casini, A., Cobine, P. A., Cross, J. R., DeNicola, G. M., et al. (2022). Connecting copper and cancer: from transition metal signalling to metalloplasia. Nat. Rev. CANCER 22 (2), 102–113. doi:10.1038/s41568-021-00417-2
Ghaffari, K., Hashemi, M., Ebrahimi, E., and Shirkoohi, R. (2016). BIRC5 genomic copy number variation in early-onset breast cancer. Iran. Biomed. J. 20 (4), 241–245. doi:10.7508/ibj.2016.04.009
Hathcock, K. S., Laszlo, G., Dickler, H. B., Bradshaw, J., Linsley, P., and Hodes, R. J. (1993). Identification of an alternative CTLA-4 ligand costimulatory for T cell activation. SCIENCE 262 (5135), 905–907. doi:10.1126/science.7694361
Hou, J., Zhao, R., Xia, W., Chang, C. W., You, Y., Hsu, J. M., et al. (2020). PD-L1-mediated gasdermin C expression switches apoptosis to pyroptosis in cancer cells and facilitates tumour necrosis. Nat. Cell Biol. 22 (10), 1264–1275. doi:10.1038/s41556-020-0575-z
Huang, X., Ruan, G., and Sun, P. (2021). Estrogen-related receptor alpha copy number variation is associated with ovarian cancer histological grade. J. Obstet. Gynaecol. Re. 47 (5), 1878–1883. doi:10.1111/jog.14741
Huo, T. I., Hsu, C. Y., and Liu, P. H. (2018). Magic mirror on the wall: which is the best biomarker for hepatocellular carcinoma? HEPATOLOGY 67 (6), 2482–2483. doi:10.1002/hep.29869
Jaiswal, A. K. (2004). Nrf2 signaling in coordinated activation of antioxidant gene expression. FREE Radic. BIO Med. 36 (10), 1199–1207. doi:10.1016/j.freeradbiomed.2004.02.074
Jiang, Y., Huo, Z., Qi, X., Zuo, T., and Wu, Z. (2022). Copper-induced tumor cell death mechanisms and antitumor theragnostic applications of copper complexes. NANOMEDICINE-UK 17 (5), 303–324. doi:10.2217/nnm-2021-0374
Kazi, J. U., Kabir, N. N., and Ronnstrand, L. (2015). Brain-Expressed X-linked (BEX) proteins in human cancers. Biochim. Biophys. Acta. 1856 (2), 226–233. doi:10.1016/j.bbcan.2015.09.001
Khazaei, M. R., Halfter, H., Karimzadeh, F., Koo, J. H., Margolis, F. L., and Young, P. (2010). Bex1 is involved in the regeneration of axons after injury. J. Neurochem. 115 (4), 910–920. doi:10.1111/j.1471-4159.2010.06960.x
Kinross, K. M., Montgomery, K. G., Kleinschmidt, M., Waring, P., Ivetac, I., Tikoo, A., et al. (2012). An activating Pik3ca mutation coupled with Pten loss is sufficient to initiate ovarian tumorigenesis in mice. J. Clin. Invest. 122 (2), 553–557. doi:10.1172/JCI59309
Kuang, D. M., Peng, C., Zhao, Q., Wu, Y., Chen, M. S., and Zheng, L. (2010). Activated monocytes in peritumoral stroma of hepatocellular carcinoma promote expansion of memory T helper 17 cells. HEPATOLOGY 51 (1), 154–164. doi:10.1002/hep.23291
Kudo, M. (2019b). Combination cancer immunotherapy with molecular targeted agents/anti-CTLA-4 antibody for hepatocellular carcinoma. LIVER CANCER 8 (1), 1–11. doi:10.1159/000496277
Kudo, M. (2019a). Immuno-oncology therapy for hepatocellular carcinoma: current status and ongoing trials. LIVER CANCER 8 (4), 221–238. doi:10.1159/000501501
Kurebayashi, Y., Ojima, H., Tsujikawa, H., Kubota, N., Maehara, J., Abe, Y., et al. (2018). Landscape of immune microenvironment in hepatocellular carcinoma and its additional impact on histological and molecular classification. HEPATOLOGY 68 (3), 1025–1041. doi:10.1002/hep.29904
Lee, S., Kang, H., Shin, E., Jeon, J., Youn, H., and Youn, B. (2021). BEX1 and BEX4 induce GBM progression through regulation of actin polymerization and activation of YAP/TAZ signaling. Int. J. Mol. Sci. 22 (18), 9845. doi:10.3390/ijms22189845
Li, H., Wang, J., Wu, C., Wang, L., Chen, Z. S., and Cui, W. (2020b). The combination of disulfiram and copper for cancer treatment. DRUG Discov. TODAY 25 (6), 1099–1108. doi:10.1016/j.drudis.2020.04.003
Li, X., Dai, Z., Liu, J., Sun, Z., Li, N., Jiao, G., et al. (2022). Characterization of the functional effects of ferredoxin 1 as a cuproptosis biomarker in cancer. Front. Genet. 13, 969856. doi:10.3389/fgene.2022.969856
Li, X., He, S., and Ma, B. (2020a). Autophagy and autophagy-related proteins in cancer. Mol. CANCER 19 (1), 12. doi:10.1186/s12943-020-1138-4
Liao, P., Huang, W. H., Cao, L., Wang, T., and Chen, L. M. (2022). Low expression of FOXP2 predicts poor survival and targets caspase-1 to inhibit cell pyroptosis in colorectal cancer. J. CANCER 13 (4), 1181–1192. doi:10.7150/jca.62433
Lindblad, O., Li, T., Su, X., Sun, J., Kabir, N. N., Levander, F., et al. (2015). BEX1 acts as a tumor suppressor in acute myeloid leukemia. Oncotarget 6 (25), 21395–21405. doi:10.18632/oncotarget.4095
Liu, Q., Li, J., Zhang, W., Xiao, C., Zhang, S., Nian, C., et al. (2021). Glycogen accumulation and phase separation drives liver tumor initiation. Cell 184 (22), 5559–5576.e19. doi:10.1016/j.cell.2021.10.001
Liu, R. M., and Gaston, P. K. (2010). Oxidative stress and glutathione in TGF-beta-mediated fibrogenesis. FREE Radic. BIO Med. 48 (1), 1–15. doi:10.1016/j.freeradbiomed.2009.09.026
Llovet, J. M., Kelley, R. K., Villanueva, A., Singal, A. G., Pikarsky, E., Roayaie, S., et al. (2021). Hepatocellular carcinoma. Nat. Rev. Dis. Prim. 7 (1), 6. doi:10.1038/s41572-020-00240-3
Llovet, J. M., Ricci, S., Mazzaferro, V., Hilgard, P., Gane, E., Blanc, J. F., et al. (2008). Sorafenib in advanced hepatocellular carcinoma. NEW Engl. J. Med. 359 (4), 378–390. doi:10.1056/NEJMoa0708857
Llovet, J. M., Zucman-Rossi, J., Pikarsky, E., Sangro, B., Schwartz, M., Sherman, M., et al. (2016). Hepatocellular carcinoma. Nat. Rev. Dis. Prim. 2, 16018. doi:10.1038/nrdp.2016.18
Lv, H., Zhen, C., Liu, J., Yang, P., Hu, L., and Shang, P. (2019). Unraveling the potential role of glutathione in multiple forms of cell death in cancer therapy. Oxid. Med. Cell Longev. 2019, 3150145. doi:10.1155/2019/3150145
Ma, Y. L., Yang, Y. F., Wang, H. C., Yang, C. C., Yan, L. J., Ding, Z. N., et al. (2023). A novel prognostic scoring model based on copper homeostasis and cuproptosis which indicates changes in tumor microenvironment and affects treatment response. Front. Pharmacol. 14, 1101749. doi:10.3389/fphar.2023.1101749
Mano, Y., Yoshio, S., Shoji, H., Tomonari, S., Aoki, Y., Aoyanagi, N., et al. (2019). Bone morphogenetic protein 4 provides cancer-supportive phenotypes to liver fibroblasts in patients with hepatocellular carcinoma. J. Gastroenterol. 54 (11), 1007–1018. doi:10.1007/s00535-019-01579-5
Mariani, D., Ghasemishahrestani, Z., Freitas, W., Pezzuto, P., Costa-da-Silva, A. C., Tanuri, A., et al. (2021). Antitumoral synergism between a copper(II) complex and cisplatin improves in vitro and in vivo anticancer activity against melanoma, lung and breast cancer cells. BBA-GEN Subj. 1865 (10), 129963. doi:10.1016/j.bbagen.2021.129963
Mortada, W. I., Awadalla, A., Khater, S., Ahmed, A., Hamam, E. T., El-Zayat, M., et al. (2020). Copper and zinc levels in plasma and cancerous tissues and their relation with expression of VEGF and HIF-1 in the pathogenesis of muscle invasive urothelial bladder cancer: A case-controlled clinical study. Environ. Sci. Pollut. R. 27 (13), 15835–15841. doi:10.1007/s11356-020-08113-8
Mou, Y., Wang, J., Wu, J., He, D., Zhang, C., Duan, C., et al. (2019). Ferroptosis, a new form of cell death: opportunities and challenges in cancer. J. Hematol. Oncol. 12 (1), 34. doi:10.1186/s13045-019-0720-y
Naderi, A., Teschendorff, A. E., Beigel, J., Cariati, M., Ellis, I. O., Brenton, J. D., et al. (2007). BEX2 is overexpressed in a subset of primary breast cancers and mediates nerve growth factor/nuclear factor-kappaB inhibition of apoptosis in breast cancer cell lines. CANCER Res. 67 (14), 6725–6736. doi:10.1158/0008-5472.CAN-06-4394
Nik-Zainal, S., Davies, H., Staaf, J., Ramakrishna, M., Glodzik, D., Zou, X., et al. (2016). Landscape of somatic mutations in 560 breast cancer whole-genome sequences. NATURE 534 (7605), 47–54. doi:10.1038/nature17676
Pang, J. S., Li, Z. K., Lin, P., Wang, X. D., Chen, G., Yan, H. B., et al. (2019). The underlying molecular mechanism and potential drugs for treatment in papillary renal cell carcinoma: A study based on TCGA and cmap datasets. Oncol. Rep. 41 (4), 2089–2102. doi:10.3892/or.2019.7014
Pariyar, M., Johns, A., Thorne, R. F., Scott, R. J., and Avery-Kiejda, K. A. (2021). Copy number variation in triple negative breast cancer samples associated with lymph node metastasis. NEOPLASIA 23 (8), 743–753. doi:10.1016/j.neo.2021.05.016
Pushpakom, S., Iorio, F., Eyers, P. A., Escott, K. J., Hopper, S., Wells, A., et al. (2019). Drug repurposing: progress, challenges and recommendations. Nat. Rev. DRUG Discov. 18 (1), 41–58. doi:10.1038/nrd.2018.168
Quentmeier, H., Tonelli, R., Geffers, R., Pession, A., Uphoff, C. C., and Drexler, H. G. (2005). Expression of BEX1 in acute myeloid leukemia with MLL rearrangements. LEUKEMIA 19 (8), 1488–1489. doi:10.1038/sj.leu.2403820
Resaz, R., Vanni, C., Segalerba, D., Sementa, A. R., Mastracci, L., Grillo, F., et al. (2014). Development of hepatocellular adenomas and carcinomas in mice with liver-specific G6Pase-alpha deficiency. Dis. Model Mech. 7 (9), 1083–1091. doi:10.1242/dmm.014878
Roos, N. J., Aliu, D., Bouitbir, J., and Krahenbuhl, S. (2020). Lapatinib activates the kelch-like ECH-associated protein 1-nuclear factor erythroid 2-related factor 2 pathway in HepG2 cells. Front. Pharmacol. 11, 944. doi:10.3389/fphar.2020.00944
Sagawa, H., Naiki-Ito, A., Kato, H., Naiki, T., Yamashita, Y., Suzuki, S., et al. (2015). Connexin 32 and luteolin play protective roles in non-alcoholic steatohepatitis development and its related hepatocarcinogenesis in rats. CARCINOGENESIS 36 (12), 1539–1549. doi:10.1093/carcin/bgv143
Satilmis, B., Sahin, T. T., Cicek, E., Akbulut, S., and Yilmaz, S. (2021). Hepatocellular carcinoma tumor microenvironment and its implications in terms of anti-tumor immunity: future perspectives for new therapeutics. J. Gastrointest. CANC 52 (4), 1198–1205. doi:10.1007/s12029-021-00725-8
Seager, R. J., Hajal, C., Spill, F., Kamm, R. D., and Zaman, M. H. (2017). Dynamic interplay between tumour, stroma and immune system can drive or prevent tumour progression. Converg. Sci. Phys. Oncol. 3, 034002. doi:10.1088/2057-1739/aa7e86
Shen, H., and Wang, W. (2021). Effect of glutathione liposomes on diabetic nephropathy based on oxidative stress and polyol pathway mechanism. J. LIPOSOME Res. 31 (4), 317–325. doi:10.1080/08982104.2020.1780607
Shi, X., Qu, M., Jin, X., Liu, L., Meng, F., and Shen, H. (2021). Relationship between TSHR, BRAF and PIK3CA gene copy number variations and thyroid nodules. ENDOCRINE 73 (1), 116–124. doi:10.1007/s12020-020-02587-9
Shimomura, Y., Takaki, A., Wada, N., Yasunaka, T., Ikeda, F., Maruyama, T., et al. (2017). The serum oxidative/anti-oxidative stress balance becomes dysregulated in patients with non-alcoholic steatohepatitis associated with hepatocellular carcinoma. Intern. Med. 56 (3), 243–251. doi:10.2169/internalmedicine.56.7002
Song, W., Ren, J., Xiang, R., Kong, C., and Fu, T. (2021). Identification of pyroptosis-related subtypes, the development of a prognosis model, and characterization of tumor microenvironment infiltration in colorectal cancer. ONCOIMMUNOLOGY 10 (1), 1987636. doi:10.1080/2162402X.2021.1987636
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA-CANCER J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660
Thai, S. F., Jones, C. P., Robinette, B. L., Ren, H., Vallant, B., Fisher, A., et al. (2021). Effects of copper nanoparticles on mRNA and small RNA expression in human hepatocellular carcinoma (HepG2) cells. J. Nanosci. Nanotechnol. 21 (10), 5083–5098. doi:10.1166/jnn.2021.19328
Tian, L., and Liao, Y. (2022). Identification of G6PC as a potential prognostic biomarker in hepatocellular carcinoma based on bioinformatics analysis. MEDICINE 101 (33), e29548. doi:10.1097/MD.0000000000029548
Torrens, L., Montironi, C., Puigvehi, M., Mesropian, A., Leslie, J., Haber, P. K., et al. (2021). Immunomodulatory effects of lenvatinib Plus anti-programmed cell death protein 1 in mice and rationale for patient enrichment in hepatocellular carcinoma. HEPATOLOGY 74 (5), 2652–2669. doi:10.1002/hep.32023
Trott, O., and Olson, A. J. (2010). AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 31 (2), 455–461. doi:10.1002/jcc.21334
Tsvetkov, P., Coy, S., Petrova, B., Dreishpoon, M., Verma, A., Abdusamad, M., et al. (2022). Copper induces cell death by targeting lipoylated TCA cycle proteins. SCIENCE 375 (6586), 1254–1261. doi:10.1126/science.abf0529
Turley, S. J., Cremasco, V., and Astarita, J. L. (2015). Immunological hallmarks of stromal cells in the tumour microenvironment. Nat. Rev. Immunol. 15 (11), 669–682. doi:10.1038/nri3902
Valk, E., Rudd, C. E., and Schneider, H. (2008). CTLA-4 trafficking and surface expression. TRENDS Immunol. 29 (6), 272–279. doi:10.1016/j.it.2008.02.011
Vilar, M., Murillo-Carretero, M., Mira, H., Magnusson, K., Besset, V., and Ibanez, C. F. (2006). Bex1, a novel interactor of the p75 neurotrophin receptor, links neurotrophin signaling to the cell cycle. EMBO J. 25 (6), 1219–1230. doi:10.1038/sj.emboj.7601017
Vogel, A., and Saborowski, A. (2020). Current strategies for the treatment of intermediate and advanced hepatocellular carcinoma. CANCER Treat. Rev. 82, 101946. doi:10.1016/j.ctrv.2019.101946
Wang, Q., Liang, N., Yang, T., Li, Y., Li, J., Huang, Q., et al. (2021). DNMT1-mediated methylation of BEX1 regulates stemness and tumorigenicity in liver cancer. J. Hepatol. 75 (5), 1142–1153. doi:10.1016/j.jhep.2021.06.025
Xie, L., Zeng, Y., Dai, Z., He, W., Ke, H., Lin, Q., et al. (2018). Chemical and genetic inhibition of STAT3 sensitizes hepatocellular carcinoma cells to sorafenib induced cell death. Int. J. Biol. Sci. 14 (5), 577–585. doi:10.7150/ijbs.22220
Xing, R., Gao, J., Cui, Q., and Wang, Q. (2021). Strategies to improve the antitumor effect of immunotherapy for hepatocellular carcinoma. Front. Immunol. 12, 783236. doi:10.3389/fimmu.2021.783236
Yan, C., Niu, Y., Ma, L., Tian, L., and Ma, J. (2022). System analysis based on the cuproptosis-related genes identifies LIPT1 as a novel therapy target for liver hepatocellular carcinoma. J. Transl. Med. 20 (1), 452. doi:10.1186/s12967-022-03630-1
Yang, H., Magilnick, N., Ou, X., and Lu, S. C. (2005). Tumour necrosis factor alpha induces co-ordinated activation of rat GSH synthetic enzymes via nuclear factor kappaB and activator protein-1. Biochem. J. 391 (2), 399–408. doi:10.1042/BJ20050795
Yang, Q., Zeng, S., and Liu, W. (2023). Roles of cuproptosis-related gene DLAT in various cancers: A bioinformatic analysis and preliminary verification on pro-survival autophagy. PEERJ 11, e15019. doi:10.7717/peerj.15019
Zhang, Z., Zeng, X., Wu, Y., Liu, Y., Zhang, X., and Song, Z. (2022). Cuproptosis-related risk score predicts prognosis and characterizes the tumor microenvironment in hepatocellular carcinoma. Front. Immunol. 13, 925618. doi:10.3389/fimmu.2022.925618
Zhao, X., Chen, J., Yin, S., Shi, J., Zheng, M., He, C., et al. (2022). The expression of cuproptosis-related genes in hepatocellular carcinoma and their relationships with prognosis. Front. Oncol. 12, 992468. doi:10.3389/fonc.2022.992468
Zheng, S. G., Wang, J. H., Stohl, W., Kim, K. S., Gray, J. D., and Horwitz, D. A. (2006). TGF-beta requires CTLA-4 early after T cell activation to induce FoxP3 and generate adaptive CD4+CD25+ regulatory cells. J. Immunol. 176 (6), 3321–3329. doi:10.4049/jimmunol.176.6.3321
Zhu, A. X., Finn, R. S., Edeline, J., Cattan, S., Ogasawara, S., Palmer, D., et al. (2018). Pembrolizumab in patients with advanced hepatocellular carcinoma previously treated with sorafenib (KEYNOTE-224): A non-randomised, open-label phase 2 trial. LANCET Oncol. 19 (7), 940–952. doi:10.1016/S1470-2045(18)30351-6
Keywords: hepatocellular carcinoma, cuproptosis, programmed cell death, immune microenvironment, single-cell RNA-sequencing, molecular docking
Citation: Cao F, Qi Y, Wu W, Li X and Yang C (2023) Single-cell and genetic multi-omics analysis combined with experiments confirmed the signature and potential targets of cuproptosis in hepatocellular carcinoma. Front. Cell Dev. Biol. 11:1240390. doi: 10.3389/fcell.2023.1240390
Received: 14 June 2023; Accepted: 24 August 2023;
Published: 08 September 2023.
Edited by:
Lin-Lin Bu, Wuhan University, ChinaReviewed by:
Kaichao Song, Chinese Academy of Medical Sciences and Peking Union Medical College, ChinaPeng Chen, Zhejiang University, China
Copyright © 2023 Cao, Qi, Wu, Li and Yang. 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: Xutong Li, MTg2NjEzOTA1NTdAMTYzLmNvbQ==; Chuang Yang, Y2h1YW5nLnlhbmdAbWVkaXppbi51bmktbGVpcHppZy5kZQ==
†These authors have contributed equally to this work and share first authorship