Autophagy-related biomarkers in hepatocellular carcinoma and their relationship with immune infiltration

Background Autophagy regulation plays vital roles in many cancers. We aimed to investigate the expression, prognostic value, and immune infiltration of autophagy-related genes in hepatocellular carcinoma (HCC) by bioinformatics analysis. Method Human autophagy-related differentially expressed genes (DEGs) between adjacent and HCC tissues were identified. We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. We also evaluated immune infiltration and the response to tumor-sensitive drugs. Finally, we verified the expression of these proteins in clinical samples by immunohistochemistry (IHC), RNA isolation and real-time reverse transcription polymerase chain reaction (RT‒PCR). Results A total of 57 autophagy-related DEGs were identified. The HUB genes (BIRC5, CDKN2A, SPP1, and IGF1) were related to the diagnosis and prognosis of HCC. The HUB genes were significantly enriched in immune-related pathways. Furthermore, correlation analysis revealed that HUB gene expression was associated with immune infiltration. We identified 35 tumor-sensitive drugs targeting the HUB genes. Finally, by IHC, we discovered that the protein of CDKN2A, BIRC5, and SPP1 were upregulated in HCC tissues, while IGF1 was downregulated in HCC tissues compared with the levels in paracarcinoma tissues; by RT‒PCR, we discovered that the mRNA of CDKN2A, BIRC5, and SPP1 were upregulated in HCC tissues, while the mRNA of IGF1 was downregulated in HCC tissues compared with the levels in paracarcinoma tissues. Conclusion We screened and validated four autophagy-related genes associated with immune infiltration and prognosis in patients with HCC.


Introduction
Hepatocellular carcinoma (HCC) is a prevalent malignancy and ranks third in cancer-related mortality worldwide [1].Chronic infection with hepatitis C virus (HCV) or hepatitis B virus (HBV) is the leading cause of HCC [2].Among primary liver cancers, HCC is the most common type [3].Considering the harmfulness and complexity of HCC, identifying potential molecular biomarkers and therapeutic targets is vital for the early diagnosis of HCC.
Autophagy is a lysosome-dependent pathway that maintains cellular homeostasis and controls cellular components by promoting the clearance of misfolded proteins and damaged organelles [4].Autophagy plays dual roles in HCC, as it not only inhibits the formation of malignant tumors but also contributes to the persistence of cancer [5].Evidence suggests that autophagy is critical for tumor growth, metastasis, and therapeutic resistance [6].For example, autophagy

Functional enrichment analysis of autophagy-related genes
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed with the R package clusterProfiler (4.4.4).A p < 0.05 was considered to indicate significant enrichment.The visual bubble plots for the GO-BP and KEGG enrichment analyses were generated with the R package ggplot2 (3.3.6).

Protein-protein interaction (PPI) network construction and hub gene identification
The String database (https:// string-db.org/) was used to construct a PPI network for autophagy-related genes with a minimum required score of 0.4.Hub genes were further characterized by cytoHubba.

Human Protein Atlas (HPA)
The HPA database is an open database that allows researchers to explore the human proteome.In this study, we used the HPA database (http:// www.prote inatl as.org/) to verify the protein expression of four hub genes selected from adjacent noncancerous tissues and hepatocellular carcinoma tissues by immunohistochemistry.

Verification of the mRNA expression of IGF1, CDKN2A, BIRC5 and SPP1
The differences in the mRNA expression levels of IGF1, CDKN2A, BIRC5, and SPP1 between tumor and nontumor tissues were verified from the GSE84402 dataset, in which 28 samples were included, consisting of 14 hepatocellular carcinoma samples and 14 nontumor samples.
The TCGA database was also used to validate the expressions of IGF1, CDKN2A, BIRC5, and SPP1 in HCC.There were 424 samples, including 374 tumor tissues and 50 tumor-adjacent tissues.Statistical significance was indicated by p < 0.05.

Analysis of the sensitivity to potential therapeutic drugs
Using the CADSP database (https:// smuon co.shiny apps.io/), we screened 288 drugs and identified antitumor drugs that were relatively sensitive to the DEGs.

Sample collection
Tumor tissue and adjacent healthy tissue samples were collected from 22 patients with hepatocellular carcinoma at the Fifth Affiliated Hospital of Zhengzhou University from July 2023 to August 2023.

RNA isolation and real-time reverse transcription polymerase chain reaction (RT-PCR)
Total RNA was extracted from the paraffin-embedded tumor and adjacent carcinoma tissues using an RNA kit (G3013, Wuhan Servicebio Technology Co., Ltd.).cDNA was amplified using primers specific for the reverse transcription kit (G3337, Wuhan Servicebio Technology Co., Ltd.).RT-PCR was performed using a LightCycler 480 (Roche) under the

Identification of DEGs related to autophagy
Our workflow is shown in Fig. 1.The GSE112790 dataset included 183 HCC tumors and 15 adjacent liver samples.The differentially expressed genes (DEGs) were screened with the following criteria: |logFC|> 1 and adj P < 0.05.As shown in Fig. 2A, these autophagy-related genes could distinguish between paracarcinoma tissue samples and HCC tissue samples, suggesting that autophagy-related genes play essential roles in the development of HCC.A total of 1408 differentially expressed mRNAs (DEMs) were identified between cancer tissue samples and paracarcinoma tissue samples in the dataset (Fig. 2B).Venn diagrams were drawn using an online Venn diagrams tool.There were 57 autophagy-related DEGs (Fig. 2C), including 25 upregulated DEGs (Fig. 2D, Table 1) and 32 downregulated DEGs (Fig. 2E, Table 1).

Enrichment analysis of autophagy-related genes
In this study, we performed gene set enrichment analysis to observe the functional enrichment of autophagy-related DEGs in HCC (Table 1).
Gene Ontology (GO) enrichment analyses were conducted for the upregulated (Table 2) and downregulated DEGs related to autophagy (Table 3).
The top 10 pathways identified by the GO analyses are displayed in a bubble diagram.The coupregulated DEGs were mostly enriched in protein catabolic process, macroautophagy, peptidyl serine phosphorylation, peptidyl-serine modification and regulation of cellular protein catabolism process (Fig. 3A), and the codownregulated DEGs were mostly enriched in regeneration, response to extracellular stimuli, epithelial cell proliferation, response to nutrient levels, and positive regulation of cellular catabolic processes and regeneration (Fig. 3C).
Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted for the upregulated (Table 2) and downregulated DEGs related to autophagy (Table 3).
The top 10 pathways identified by the KEGG analyses are displayed in a bubble diagram.The coupregulated DEGs were mostly enriched in fluid shear stress and atherosclerosis, focal adhesion and the PI3K-Akt signaling pathway (Fig. 3B), and the codownregulated DEGs were mostly enriched in autophagy-animal and cytokinecytokine receptor interactions (Fig. 3D).

Construction of a protein-protein interaction (PPI) network via STRING and identification of HUB genes via Cytoscape software
STRING was used to construct a network for the DEGs associated with autophagy (Fig. 4A).After further analysis with Cytoscape, IGF1, CDKN2A, BIRC5, and SPP1 were selected as the HUB genes and visualized (Fig. 4B).IGF1, CDKN2A, BIRC5, and SPP1 were selected as potential biomarkers (Table 4) and further validated in subsequent studies.

Overall survival analysis for patients with HCC stratified by IGF1, CDKN2A, BIRC5 and SPP1 levels
The overall survival of patients with HCC stratified by IGF1, CDKN2A, BIRC5, and SPP1 levels was analyzed by Kaplan-Meier plotter to further investigate whether IGF1, CDKN2A, BIRC5, or SPP1 had an effect on overall survival.Our findings suggested that high levels of CDKN2A, BIRC5, and SPP1 were associated with poor overall survival in patients with HCC, suggesting that BIRC5, CDKN2A, and SPP1 are associated with HCC progression and could be used as tumor biomarkers in patients with HCC.IGF1 was at lower levels in tumor tissues than in adjacent tissues, possibly because of its association with lower survival (Fig. 7E-H).

Analysis of the associations of tumor node metastasis (TNM) stages with IGF1, CDKN2A, BIRC5 and SPP1 levels
We analyzed the associations of TNM stage with IGF1, CDKN2A, BIRC5, and SPP1.We found significant differences in different T stages compared to regular level, suggesting that HUB genes may be related to the progression of HCC (Fig. 7I-L).

GSEA revealed pathways associated with IGF1, CDKN2A, BIRC5 and SPP1 expression in HCC
GSEA indicated that IGF1-related genes were mainly enriched in immune-related pathways, such as disease of the immune system, immunoregulatory interactions between lymphoid and nonlymphoid cells, drug metabolism, other enzymes, and drug ADME (Fig. 8A).CDKN2A-related genes were mainly enriched in the mitotic G1 phase and G1 transition, integrated cancer pathway, regulation of TP53 activity, and disease of programmed cell death (Fig. 8B).BIRC5-related genes were mainly enriched in the mitotic G1 phase and G1 transition, immunoregulatory interaction, lymphoid and nonlymphoid cells, retinoblastoma gene in cancer, and FceRI-mediated NFkb activation (Fig. 8C).SPP1-related genes were mainly enriched in pathways related to cancer, immunoregulatory interactions between lymphoid and nonlymphoid cells, photodynamic therapy-induced NFkb survival signaling, focal adhesion PI3K, AKT, and the mTOR signaling pathway (Fig. 8D).

Correlation analysis of IGF1, CDKN2A, BIRC5 and SPP1 with the TME in HCC
As demonstrated in Fig. 9A, patients with lower level of IGF1 had lower numbers of memory dendritic cells (DCs), aDCs, cytotoxic cells, mast cells, neutrophils, Tgd, Th1 cells, Th17 cells, and Tregs.Additionally, as depicted in Fig. 9B-D, patients in the groups with high CDKN2A, BIRC5, and SPP1 had high Th2 cells and low Th 17 cells and eosinophils.According to Fig. 9E, patients with high IGF1 had negative correlations with NK and CD56bright cells but positive correlations with neutrophils, DCs, cytotoxic cells, Tregs, Th1 cells, Th17 cells, Mast cells, Tgd and NK CD56dim cells.As demonstrated in Fig. 9F-H, patients with high CDKN2A, BIRC5 and SPP1 had negative correlations with Th17 cells and eosinophils but positive correlations with Th2 cells.Neutrophils, Th17 cells, and eosinophils may be the most prevalent differential immune cells.Th2 cells may play a role in promoting cancer.In this research, we further assessed the correlations between immune infiltration levels and the expression levels of IGF1, CDKN2A, BIRC5, and SPP1 in HCC via the TIMER database.As shown in Fig. 10, IGF1 was negatively correlated with the infiltration of macrophages (p = 1.13 × 10 − 2); moreover, CDKN2A, BIRC5, and SPP1 levels were positively correlated with the infiltration of macrophages, B cells, CD4 + T cells, CD8 + T cells, dendritic cells, and neutrophils.Figure 11A-D shows that there was variability in the stromal, immune, and estimated scores of patients with HCC, indicating some differences in the TME according to the level of IGF1, CDKN2A, BIRC5, and SPP1.Various immune checkpoints were associated with IGF1, CDKN2A, BIRC5, and SPP1 (Fig. 11E).

Possible sensitivity to therapeutic drugs
Using the CADSP database, we screened 288 drugs and identified antitumor drugs that were relatively sensitive to DEGs.Patients with high level of CDKN2A and BIRC5 had the lowest half maximal inhibitory concentrations (IC50) for Epothilone B, while SPP1 was inhibited by bortezomib.Shikonin (SHK) had the lowest IC50 in patients in the group with high IGF1 level.We identified 35 tumor-sensitive drugs targeting the HUB genes (Fig. 12).

External validation in patients with HCC from the hospital
We validated the protein levels of IGF1, CDKN2A, BIRC5, and SPP1 in HCC by IHC (Fig. 13).CDKN2A, BIRC5, and SPP1 were upregulated in HCC tissues by IHC.However, IGF1 was downregulated in HCC tissues compared with adjacent tissues.
We also validated the expressions of IGF1, CDKN2A, BIRC5, and SPP1 in HCC by RT-PCR (Fig. 14).CDKN2A, SPP1, and BIRC5 were upregulated in HCC tissues according to RT-PCR; IGF1 was downregulated in HCC tissues compared with adjacent tissues.

Discussion
In this study, we investigated autophagy-related genes and evaluated immune cell infiltration in HCC via bioinformatics methods, which could provide new perspectives for immunotherapy.We first collected the GSE112790 dataset, which included 183 HCC tumor samples and 15 adjacent liver samples.Then, 57 autophagy-related genes were identified between the tumor and adjacent tissue samples.Subsequently, we constructed a PPI network and identified four HUB genes, namely, IGF1, CDKN2A, BIRC5, and SPP1.The results of these analyses were confirmed by IHC and RT-PCR.The protein of CDKN2A, BIRC5, and SPP1 were upregulated in HCC tissues by IHC.However, IGF1 was downregulated in HCC tissues compared with adjacent tissues.The mRNA of CDKN2A, SPP1, and BIRC5 were upregulated in HCC tissues according to RT-PCR; the mRNA of IGF1 was downregulated in HCC tissues compared with adjacent tissues.
IGF1 is a RAS pathway-related protein that is present in plasma in all tissues [8].The protein encoded by this gene is similar to insulin in function and structure and is a member of a family of proteins involved in mediating growth and the development of multiple diseases [9].Our study revealed that IGF1, an autophagy-related protein, is expressed at lower levels in HCC tissues than in para-adjacent tissues and may be involved in the prognostic regulation of HCC.
The IGF1 level was greater in adjacent nonneoplastic liver tissues than in tumor tissues, which was correlated with significantly worse survival after HCC resection [10].Multivariate analysis revealed that higher IGF-1 in adjacent nonneoplastic liver tissue remained an independent predictor of poor outcome [10].Activation of the IGF1 signaling cascade leads to upregulation of downstream mitogens, including MAPK and PI3K/Akt.During liver carcinogenesis, the secretion of IGF1 by adjacent hepatocytes may lead to paracrine stimulation of HCC and more aggressive tumor behavior [11].Increased levels of insulin-like growth factor 1 are associated with cancer development, and circulating IGF1 is positively associated with breast cancer risk [12].The serum IGF1 concentration is greater in patients with melanoma than in healthy individuals [13].Notably, elevated IGF1 is closely associated with drug resistance.After gefitinib treatment, liver cells exhibit increased phosphorylation of IGF1R and AKT in HCC, indicating that activation of IGF1R signaling may contribute to gefitinib resistance [14].Suppressor of IGF1 signaling plays essential roles in cancer therapy.Therapeutic strategies, including IGF1 and IGF1R antibodies and tyrosine kinase inhibitors, have been developed, and both show good responses [15].Targeting IGF1-associated lncRNAs has potential in cancer therapy.For example, silencing lncNEAT1 inhibited chemoresistance in gastric cancer [16].These findings indicate that IGF1 plays a specific role in the development of cancer.Most importantly, elevated IGF1 is strongly associated with resistance to various chemotherapies, which may be modulated by IGF1-associated lncRNAs [17].Finally, targeting IGF1-associated lncRNAs is also expected to prevent and reverse chemoresistance.More precise in vivo and in vitro experiments are needed to verify this mechanism in the future.
In the present study, GSEA indicated that IGF1-related genes were mainly enriched in immune-related pathways, such as disease of the immune system, immunoregulatory interactions between lymphoid and nonlymphoid cells, drug metabolism, other enzymes, and drug ADME.Therefore, we hypothesized that regulating IGF1, CDKN2A, BIRC5, and SPP1 levels may affect the tumor immune microenvironment in HCC.Previous reports have demonstrated a close relationship between immunity and autophagy [18].Recent studies have shown that autophagy can direct the immune response The results showed that IGF1, CDKN2A, BIRC5, and SPP1 levels correlated with infiltrating macrophage levels (Fig. 10).
In the tumor-immune microenvironment, tumor-associated macrophages (TAMs) are recruited and activated by different chemokines and differentiate into proinflammatory and antitumor-active M1-TAMs and into anti-inflammatory and tumor-promoting M2-TAMs.The two forms of polarization are interconverted in the presence of certain stimuli [20].TAMs can induce autophagy in HCC cells and attenuate the toxic effects of oxaliplatin.This autophagy-mediated drug resistance mechanism provides a novel therapeutic strategy [21].Regulating autophagy in TAMs may be a promising strategy for inhibiting HCC progression.Modulated by genetic engineering techniques, macrophages can target specific antigens, such as CD19, CD22, and Her2, to identify tumor cells [22].Chimeric antigen receptor macrophage (CAR-M) cells can phagocytose tumor cells, secrete proinflammatory cytokines to alter the tumor immune microenvironment, and present tumor antigens to T cells to activate the immune response [23].CAR-M technology is constantly evolving; however, due to the heterogeneity of HCC, it remains challenging to explore specific targets in liver tumor cells and to engineer macrophages in the liver TME.We also found significant associations of CD276, CD40, and TNFRSF4 with IGF1, CDKN2A, BIRC5 and SPP1.Targeting these immune checkpoints, such as inhibitors that target checkpoint molecules, may enhance the body's immune function.Combining ICB with other treatments may improve the immunological conditions of the TME, thus enhancing the antitumor immune response.Immune checkpoint blockade (ICB) has been successful in treating more immunogenic tumors [24].However, ICB is still largely ineffective in patients with tumors that are not infiltrated by immune cells (cold tumors).ICB can be combined with compounds capable of converting noninflammatory tumors into inflammatory tumors to further increase the number of appropriate tumor types.This may, in turn, increase the sensitivity of these tumors to ICB therapy.Finally, by using the CADSP database, we screened 288 drugs and identified antitumor drugs that are relatively sensitive to IGF1, CDKN2A, BIRC5, and SPP1.Epothilone B had the lowest IC50 for CDKN2A and BIRC5, while bortezomib had the lowest IC50 for SPP1.The patients in the group with high levels of IGF1 had the lowest IC50 for SHK; thus, these patients are more sensitive to SHK.Although the effect of chemotherapy on liver cancer is limited, it can improve the prognosis of specific patients and can extend their lives to a certain extent.SHK is a natural molecule isolated from peramina that shows great potential in anticancer therapy.SHK can also modulate the immunosuppressive TME by inhibiting tumor cell glycolysis and repolarizing tumor-associated macrophages [25].M2-TAMs exert anti-inflammatory and tumorigenic effects.SHK can reduce the production of lactate acid, an essential driver of TAM2 polarization, thus repolarizing M2-TAMs [26].PD-L1 is the primary ligand for PD-1.PD-L1 expression is an immune evasion mechanism in various malignancies [27].SHK promotes PD-L1 degradation and suppresses the immune escape of pancreatic cancer cells by inhibiting the NF-κB/STAT3 and NF-κB/CSN5 signaling pathways [28].Autophagy is a form of programmed cell death; therefore, the regulation of autophagy can be an effective intervention strategy for cancer therapy.Enhanced autophagy is accompanied by RIPK1-and RIPK3-dependent necrotizing apoptosis induced by SHK stimulation [29].However, SHK is a hydrophobic natural molecule with unsatisfactory solubility and rapid intestinal absorption, resulting in low oral bioavailability [30].Given the shortcomings of SHK, significant efforts have been made to develop various nanodrug delivery systems.A mannosylated lactoferrin nanosystem (Man-LF NPs) was prepared to deliver SHK and JQ1, which target colon cancer cells and TAMs [31].When combined with other immunodrugs, JQ1 effectively reduced PD-L1 expression in tumor cells [32].
In this study, we screened and validated four autophagy-related genes associated with immune infiltration and prognosis in patients with HCC.It provided new insights into the HCC immune microenvironment and provide new perspectives on autophagy gene-targeted immunotherapy.
A limitation of this study is that our data are mainly derived from public data, which is available in limited amounts; moreover, the clinical validation samples were limited.Trials with sufficiently large sample sizes are needed for validation in the future.

Conclusions
We identified and validated autophagy genes associated with HCC and revealed that HUB genes were correlated with the TME in HCC.These findings can provide new perspectives on autophagy gene-targeted immunotherapy.

Fig. 1
Fig. 1 Flow chart of the study

Fig. 2
Fig. 2 Identification of autophagy-related genes using GEO and TCGA datasets.A PCA array.B DEGs between the 183 HCC tumors and 15 adjacent liver samples.C Heatmap showing the expression levels of 57 autophagy-related DEGs in normal and tumor tissues.D Twentyfive overlapping genes identified as autophagy-related upregulated DEGs.E Thirty-two overlapping genes identified as autophagy-related downregulated DEGs

Fig. 3
Fig.3Enrichment analysis of autophagy-related genes by TCGA database.A GO-BP enrichment analysis of autophagy-related upregulated genes.B KEGG enrichment analysis of autophagy-related upregulated genes.C GO-BP enrichment analysis of autophagy-related downregulated genes.D KEGG enrichment analysis of autophagy-related downregulated genes

Fig. 4
Fig.4 Identification of HUB genes by the STRING database.A PPI network of overlapping DEGs.B The most significant 15 node degree genes calculated by the cytoHubba app in Cytoscape.IGF1, CDKN2A, BIRC5, and SPP1 were selected as the HUB genes.The node color intensities representing different genes correlate with the degree of expression values

Fig. 5
Fig. 5 Expression levels of IGF1, CDKN2A, BIRC 5, and SPP1 in HCC from the GEO and TCGA databases.A-H IGF1, CDKN2A, BIRC5, and SPP1 mRNA expression levels are based on the TCGA database, including paired (A-D) and nonpaired samples (E-H).I-L The mRNA expression levels of IGF1, CDKN2A, BIRC5, and SPP1 mRNA are based on GSE84402

Fig. 6
Fig. 6 Levels of IGF1, CDKN2A, BIRC5, and SPP1 in normal individuals and patients with HCC from the HPA database

Fig. 13 Fig. 14
Fig. 13 Levels of IGF1, CDKN2A, BIRC5, and SPP1 in patient tissues by IHC.Levels of CDKN2A, BIRC5, and SPP1 in HCC were increased compared with tumor-adjacent tissues.However, IGF1 expression was downregulated in HCC tissues compared with adjacent tissues

Table 1
Lists of DEGs upregulated/downregulated in dataset GSE112790

Table 3
GO-BP enrichment and KEGG pathway analyses of downregulated DEGs associated with autophagy Ontology