An Autophagy-Related Prognostic Signature and Verification that ATIC Inhibits Autophagy Through the AMPK/mTOR Pathway in Hepatocellular Carcinoma


 Background: Autophagy regulates many cell functions related to cancer, ranging from cell proliferation and angiogenesis to metabolism. Due to the close relationship between autophagy and tumors, we investigated the predictive value of autophagy-related genes.Methods: Data from patients with hepatocellular carcinoma were obtained from the TCGA and ICGC databases. Regression analysis of differentially expressed genes was performed. Based on a prognostic model, patients were divided into a high-risk or low-risk group. Kaplan-Meier survival analyses of patients in the two groups were conducted. The immune landscapes, as shown by the single-sample gene set enrichment analysis (ssGSEA), had different patterns in two groups. The prognostic model was verified using the ICGC database and clinical data of patients collected from Zhongnan Hospital. In the results of multivariate COX regression analysis, ATIC has the largest hazard ratio, so we studied the effect of ATIC on autophagy and tumor progression through experiments in vitro and in vivo.Results: Fifty-eight autophagy-related genes were differentially expressed (FDR < 0.05, logFC > 1); 23 genes were related to the prognosis of patients. A prognostic model based on 12 genes (ATG10, ATIC, BIRC5, CAPN10, FKBP1A, GAPDH, HDAC1, PRKCD, RHEB, SPNS1, SQSTM1 and TMEM74) was constructed using the method of multivariate COX regression analysis. For the high-risk group and low-risk group distinguished by the model, there was a significant difference in survival between the two groups (P < 0.001). The model had good predictive power (AUC > 0.7). Risk-related genes were related to the terms type II IFN response, MHC class I (P < 0.001) and HLA (P < 0.05). ATIC was confirmed to inhibit autophagy and promote the proliferation, invasion and metastasis of liver cancer cells through the AMPK/mTOR signaling pathway.Conclusions: The prediction model can effectively predict the survival time of liver cancer patients. The risk score reflects patient immune cell features and immune status. ATIC inhibits autophagy and promotes the progression of liver cancer through the AMPK/mTOR signaling pathway.


Introduction
Worldwide, liver cancer is responsible for the 4th highest number of cancer-related deaths, and the number of new liver cancer cases every year is the 6th highest among all cancers [1]. The 5-year survival rate of liver cancer is only 18%, and it is the second most lethal cancer after pancreatic cancer.
Furthermore, the long-term survival rate of liver cancer varies greatly among countries and regions [2,3].
The molecular characteristics of liver cancer cells and the clinical prognosis of liver cancer are very different from those of other cancers, which leads to challenges in predicting the prognosis of liver cancer patients [4,5].
Autophagy can improve the viability of tumor cells and exacerbate tumor progression [6]. Numerous studies have reported that there is a connection between autophagy and the human immune system and immune-related functions. Recent studies have reported that autophagy can reduce cell surface inhibition of antitumor immune responses and promote immune evasion in pancreatic cancer [7]. Autophagy can regulate the balance, activation and differentiation of a variety of immune cells, as well as the secretion of a variety of cytokines [8]. In view of the close relationship between autophagy and immunity, we analyzed autophagy and immune-related functions in this study. In short, autophagy plays an important role in tumorigenesis and tumor progression, and immunotherapy is an important strategy for tumor treatment in the future [9].
Due to the close relationship between autophagy and tumors, we designed this study to explore the relationship between autophagy and the prognosis of liver cancer. Here, we identify genes that are closely related to the prognosis of liver cancer patients through bioinformatic methods such as COX regression analysis and survival analysis, and provide ideas for subsequent mechanistic studies, and establish models for predicting survival.
In the results of multivariate COX regression analysis, 5-minoimidazole-4-carboxamide ribonucleotide formyltransferase/inosine monophosphate (IMP) cyclohydrolase (ATIC) has the largest hazard ratio, suggesting that ATIC is closely related to patient prognosis. Therefore, we chose ATIC for further veri cation. ATIC is considered to be related to autophagy, but there is no experimental evidence to prove the relationship between ATIC and autophagy [10]. In this study, we determined the relationship between ATIC and tumor cell autophagy through experiments in vivo and in vitro. We found that ATIC can inhibit autophagy and promote liver cancer progression by affecting the AMPK/mTOR pathway.
This article identi ed autophagy-related genes related to the prognosis of liver cancer, analyzed the relationship of prognosis-related autophagy genes with immune cells and immune-related functions, and constructed a model to predict the survival of patients with liver cancer in the clinic uses bioinformatic methods. The results provide future research directions. Furthermore, we have proved that ATIC can inhibit autophagy and promote liver cancer progression.

Data collection
The mRNA sequencing data of 50 normal samples and 374 liver cancer samples and clinical data of 374 liver cancer patients were downloaded from The Cancer Genome Atlas (TCGA, portal.gdc.cancer.gov).
Then, mRNA sequencing data from 240 liver cancer patients were downloaded from the International Cancer Genome Consortium (ICGC, daco.icgc.org). Excluding patients with incomplete data, a total of 365 patients from the TCGA database and 231 patients from the ICGC database were included in the analysis. A list of 256 autophagy-related genes was obtained from the Human Autophagy Database (autophagy.lu). The expression data of autophagy-related genes were extracted from TCGA sequencing data.
Identi cation of differentially expressed genes and prognosis-related genes The "limma" R package was used to perform differential analysis. The screening conditions were log 2 fold-change (logFC) > 1 and false discovery rate (FDR) < 0.05. The survival R package was used to identify prognosis-related genes. The genes that were both differentially expressed and related to prognosis were determined. A protein interaction network was generated for the overlapping genes. The "igraph" and "reshape2" R packages were used to calculate the correlations between the expression levels of overlapping genes and to determine the correlation network (correlation threshold: 0.4).

Construction of prognostic models
We used both LASSO regression and multivariate Cox regression analyses to establish prognostic models based on the TCGA data. Then, we used ICGC data to verify the established models. According to the model, ICGC samples were divided into a high-risk and low-risk group. The predictive power of the models was evaluated by survival analysis, ROC curve analysis and independent prognostic analysis. The "Rtsne" R package was used for principal components analysis (PCA), and t-stochastic neighbor embedding (t-SNE) was used for dimensionality reduction and data visualization. After the two models were tested and compared, the prognostic model established by the LASSO regression method had better predictive power. Therefore, we used the prognostic model established by LASSO regression for subsequent analyses.
Identi cation of risk-related genes, gene set enrichment analysis and immune in ltration analysis The "limma" R package was used to analyze the differentially expressed genes in the high-risk and lowrisk groups. Then, GO functional enrichment analysis and KEGG pathway enrichment analysis were performed with the risk-related genes obtained from the analysis of the TCGA and ICGC data. The gene set enrichment analysis used the "clusterPro ler" and "org.Hs.eg.db" R packages, and P value < 0.05 and q value < 0.05 were used as thresholds.
Studies have shown that tumor gene expression levels are related to local immune cell activity and immune-related functions [11]. As such, we used the single-sample gene set enrichment analysis (ssGSEA) function in the "GSVA" R package to calculate the immune-related scores of each sample and the relationship between the immune-related scores and risk scores. We used the same method to calculate the relationship between immune in ltration and the expression of the representative gene ATIC in the model.

Cell culture and transfection
Huh7, HepG2 and HCCLM3 liver cancer cells were obtained from Shanghai Cell Bank (Shanghai, China). The cell line has recently been authenticated and has recently been tested for mycoplasma. Huh7, HepG2 and HCCLM3 cells were grown in DMEM (Invitrogen, Carlsbad, CA). All media were supplemented with 10% fetal bovine serum (Invitrogen). Jtsbio (Guangzhou, China) constructed and synthesized the siRNA and shRNA lentivirus vector speci c for ATIC. The sequence of shRNA is the same as siRNA 3: For overexpression of ATIC, full-length ATIC cDNA was ampli ed by PCR and cloned into the expression vector pcDNA3.1 (Invitrogen). The empty vector was used as a control. Transfection was performed using Lipofectamine 3000 reagent (Invitrogen). Forty-eight hours after transfection, cells were harvested for further analysis.

RNA extraction and RT-qPCR analysis
Total RNA was isolated from the prepared tissues and cells using TRIzol reagent (Invitrogen).
Complementary DNA was synthesized using a PrimeScript RT kit (TaKaRa, Dalian, China). Then, SYBR Green Master Mix (Applied Biosystems, Foster City, CA, USA) was used to perform qPCR with the ABI PRISM 7900 Sequence Detection System (Applied Biosystems). The data were calculated using the 2-ΔΔCt method, and GAPDH was used as an endogenous control. The primer sequences are as follows:

Transmission electron microscopy
Conventional electron microscopy was performed as follows. Cells were xed with 2.5% glutaraldehyde and then post xed with 1% osmium tetroxide. The samples were dehydrated in an ethanol series and embedded in Embed812 resin. Ultrathin sections were cut with a diamond knife. Next, the sections were mounted on copper grids and double-stained with uranyl acetate and lead citrate. The number of autophagic vacuoles was determined for a minimum of 100 cells.

Detection of cell viability and migration ability
Cell counting kit 8 (CCK8) was used to measure cell viability. An In nite M200 spectrophotometer (Tecan) was used to detect the spectrophotometric absorbance of each sample at 450 nm. Transwell assays were used to assess cell migration. BioCoat Matrigel Invasion Chambers (BD Biosciences) were used to assess cell invasion. Three random areas were assessed to determine the cell counts for the cell migration and invasion experiments. An annexin V-FITC apoptosis detection kit (Invitrogen) was used for annexin V and propidium iodide staining of cells, and the percentage of apoptotic cells was determined by ow cytometry (Beckman, USA). To assess the cell cycle, 48 hours after transfection, the cells were stained with propidium iodide and assessed by uorescence-activated cell sorting (FACS).

Colony formation assay
The treated liver cancer cells were cultured in a 6-well plate. After 14 days, the culture medium was discarded, and the cells were washed twice with PBS. Then, paraformaldehyde was used to x the cells at room temperature for 20 minutes, 0.1 mM crystal violet (500 µl) was added, the cells were stained for 10 minutes, washed twice with PBS, and the results were analyzed under a microscope.

Xenotransplantation and lung metastasis models
For these experiments, 6-week-old BALB/c nude mice were used. The treated liver cancer cells were washed and suspended, and 1×10 7 cells/0.1 ml in a single-cell suspension were inoculated into the armpits of nude mice. The long diameter (a) and short diameter (b) of subcutaneous xenograft tumors were measured every 5 days. The tumor volume was calculated according to the formula V = ab 2 /2. After 5 weeks, the nude mice were sacri ced, and the tumors were extracted. The differently treated cells in the logarithmic growth phase were selected, and a cell suspension was made and injected through the tail vein. The nude mice were sacri ced after a month and the lung tissues were collected for hematoxylineosin (H&E) staining. Generally, 1-5*10 6 cells/200 µl were used. The animal studies were approved by the Animal Research Committee of Wuhan University.

Statistical analysis
All statistical analyses were performed using GraphPad Prism 8.0 software (GraphPad Software, San Diego, California, USA) and SPSS 24.0 software (SPSS Inc., Chicago, USA). Each experiment included at least three independent experiments. The data are expressed as the mean ± standard deviation (SD). The differences between the experimental groups were assessed by Student's t-test or one-way ANOVA. A twosided P value was calculated, and a probability level of 0.05 was considered statistically signi cant.

Results
Identi cation of differentially expressed genes related to autophagy and prognosis We rst identi ed 58 genes that were abnormally expressed in tumors compared with normal tissues, of which 4 genes had abnormally low expression and 54 genes had abnormally high expression ( Fig. 1A and B). A total of 43 prognosis-related genes were screened using univariate Cox regression analysis (Fig. 1C). The intersection of the two gene sets resulted in 23 overlapping genes ( Fig. 1D and Table. S1).
The overlapping genes were all highly expressed in tumor tissues (Fig. 1E), and they were all unfavorable factors for prognosis (Fig. 1F). Among them, BIRC5, PEA15 and RAB24 were the most signi cant genes. CAPN10 and SPNS1 were the genes with the greatest impact on prognosis. To explore the interactions and coexpression relationships between the overlapping genes, we generated a protein interaction network. MAPK3, SQSTM1, CASP8 and HSP90AB1 were the core genes in the network and showed research potential (Fig. 1G). In the coexpression network, the lines between genes represent coexpression relationships, and red represents a positive correlation. IKBKE and PRKCD had the strongest coexpression relationship (Fig. 1H).
We use data from ICGC database to verify the model established based on TCGA data. The patients in the TCGA group and the ICGC group were divided into a high-risk group and a low-risk group according to risk score. The results show that the risk score is related to the patient's disease stage and grade (Table. 1). Furthermore, the difference in the prognosis of the two groups was calculated. The results showed that the survival difference between the high-risk group and the low-risk group was more obvious for the model established by LASSO regression analysis than for the model established by Cox regression analysis ( Fig. 2D-G).
ROC risk curves were drawn and the area under the curve (AUC) values were calculated to evaluate the predictive performance of the model. For the LASSO regression model, the 1-, 2-, and 3-year AUC values for TCGA patients were 0.768, 0.714 and 0.696, respectively (Fig. 2H). The ICGC group was used to test the e cacy of the model, and the AUC values were 0.745, 0.761 and 0.739, respectively (Fig. 2I). The AUC values of the model established by multivariate Cox regression analysis were slightly lower than those of the model established by multivariate COX regression ( Fig. 2J and K). Therefore, the model established by LASSO regression analysis had better predictive performance, and the model was further veri ed and analyzed. The distribution of risk scores among patients in the TCGA cohort was symmetrical (Fig. 2L). The distribution of patient survival statuses showed that the survival time of patients in the high-risk group was shorter than that of patients in the low-risk group (P < 0.01, Fig. 2M). When ICGC patient data were used for veri cation, the same distribution was observed ( Fig. 2N and O). To show the clustering of risk scores more clearly, we conducted PCA and t-SNE analysis. The results showed that in both the TCGA and ICGC groups, patients in the high-risk and low-risk groups were clearly distributed in different clusters (Fig. 2P-S).
The risk score is an independent prognostic factor and is related to immune cells and immune function We used univariate and multivariate Cox regression analyses to determine whether the risk score was an independent prognostic indicator. The results showed that in the TCGA and ICGC groups, disease stage and risk score were independent factors affecting prognosis (P < 0.001, Fig. 3A-D). We further performed an enrichment analysis of genes with differential expression (logFC > 1, adjusted P < 0.05). The results of the GO analysis showed that the genes were mainly involved in cell cycle-related functions, including nuclear division and chromatic segregation ( Fig. 3E and F). The results of the KEGG pathway enrichment analysis showed that the autophagy-related genes were closely related to the cell cycle ( Fig. 3G and H). In addition, we also found that the autophagy-related genes were also related to immune function in the GO analysis and KEGG analysis. Related GO terms included humoral immune response, chemokine production, IL-17 signaling pathway and cytokine activity.
The differences in immune cells and functions were quanti ed and analyzed. Between the high-risk group and the low-risk group in the TCGA cohort, there were signi cant differences in the activation degree of immune cells related to cellular immunity, including macrophages, natural killer (NK) cells, T helper 1 (Th1) cells, Th2 cells and regulatory T (Treg) cells. In terms of immune-related functions, CCR and the type II IFN response, which are also related to cellular immunity, were signi cantly differentially activated between the high-risk and low-risk group. In addition, immune terms related to antigen presentation, such as APC costimulation and MHC class I, were also signi cantly differentially enriched between the two groups ( Fig. 3I and J). In the ICGC cohort, the degree of activation of macrophages and Th2 cells was obviously different between the high-risk and low-risk group, which was the same as the result in the TCGA group. The changes in immune function in the ICGC group were the same as those in the TCGA group, and signi cant enrichment of MHC class I and the type II IFN response was found in both cohorts ( Fig. 3K and L).
For the key gene ATIC in the model, we divided patients into a high expression group and a low expression group according to the expression level of ATIC. In the TCGA cohort, the expression of ATIC was related to the in ltration of CD8 + T cells, macrophages and mast cells (Fig. 3M). ATIC was related to immune checkpoints and type I and type II IFN response immune functions (Fig. 3N). In the ICGC cohort, ATIC was still associated with CD8 + T cells and NK cells but not with macrophages (Fig. 3O). In the ICGC cohort, ATIC was related to immune checkpoints and IFN response, which is consistent with the results in TCGA (Fig. 3P).
Based on the differential genes between the high-risk group and the low-risk group, we generated a connectivity map to identify small molecules with therapeutic potential using online tool cMAP (Table  S2). Among them, meclofenamic acid, UNC-0321, fananserin and icariin showed great therapeutic potential. Doxorubicin and erlotinib are used in the treatment of liver cancer. For high-risk patients, these two drugs may have a better effect ( Fig. 3Q and R).

ATIC is highly expressed in liver cancer tissues and associated with a poor prognosis
The results of the multivariate Cox regression analysis showed that ATG10, ATIC, HDAC1, RHEB, SPNS1 and TMEM74 had a signi cant effect on prognosis. Among them, the abnormally high expression of ATIC and HDAC1 was the most obvious harmful factor affecting prognosis (P < 0.05, hazard ratio > 1). The immunohistochemical results in The Human Protein Atlas database (https://www.proteinatlas.org) veri ed that the expression of ATIC, HDAC1, RHEB, SPNS1 and TMEM74 in liver cancer tissues was signi cantly higher than that in normal liver tissues ( Fig. 4A and B).
Since the hazard ratio of ATIC was greater than 1 and it had a large coe cient in the established model, we decided to functionally verify the effects of ATIC. According to follow-up and gene expression data from TCGA, high expression of ATIC was closely related to poor prognosis (P < 0.001, Fig. 4C). Compared with the normal liver cell line LO2, the HepG2 and Huh7 cell lines showed signi cantly higher expression of ATIC (Fig. 4D). We collected tumor samples and paracancerous samples of 52 liver cancer patients from the the Biological Repositories, Zhongnan Hospital of Wuhan University and veri ed that ATIC was signi cantly high expressed in liver cancer tissue samples (P < 0.001, Fig. 4E).
We analyzed the association between ATIC and the clinicopathological characteristics of liver cancer patients. As shown in Table 2, it was observed that ATIC expression was associated with lymph node invasion (P = 0.048) and prognosis (Fig. S1). However, the expression of ATIC was independent of age (P = 0.393), gender (P = 0.760) and portal vein tumor thrombosis (PVTT, P = 0.405) ( Table 2). We also analyzed the relative risk of ATIC in the prognosis of liver cancer. The results of Cox regression analysis showed that, compared with low-expressing liver cancer patients, high ATIC expression was closely related to the poor prognosis of liver cancer patients (P = 0.002). Multivariate Cox regression analysis further con rmed that the expression of ATIC was signi cantly correlated with the prognosis of liver cancer (P = 0.006) ( Table 3).   We used HepG2 and Huh7 cell lines for the next experiment. After transfection, the mRNA and protein expression levels of ATIC were signi cantly decreased (Fig. 4F). The results of westernblot are consistent with the results of PCR, showing that ATIC is highly expressed in liver cancer tissues (Fig. 4G) and liver cancer cell lines (Fig. 4H). siRNA not only reduced the mRNA level of ATIC, but also reduced the protein level (Fig. 4I).

ATIC inhibits autophagy and apoptosis of tumor cells and promotes malignant tumor behavior
siRNA negative control (siRNA NC) and siRNA 3 were used to transfect the liver cancer cell lines Huh7 and HepG2, and CCK8 was used to detect cell proliferation. After ATIC knockdown, the proliferation ability of the cells decreased signi cantly (P < 0.05, Fig. 5A). The results of the colony formation experiment showed that after 2 weeks of culture, the number of colonies formed by the cells with ATIC knockdown was signi cantly lower than that formed by control cells (P < 0.05, Fig. 5B). The transwell and wound healing experiments showed that knockdown of ATIC expression reduced the invasion and migration ability of cells ( Fig. 5C and D). Consistent with the enrichment analysis results, knockdown of ATIC affected the cell cycle of tumor cells (Fig. 5E). The levels of LC3 I and LC3 II re ect cell autophagy. An increase in the proportion of LC3 II represents an increase in the level of autophagy. Apoptosis and autophagy are closely related. Knockdown of ATIC induced cell apoptosis (Fig. 5F). After transfection of pcDNA/ATIC, the level of ATIC protein increased signi cantly (Fig. 5G). Correspondingly, after ATIC knockdown, the levels of cleaved caspase 3 in cells increased, and at the same time, the proportion of the autophagy-related protein LC3 II increased (Fig. 5H). The levels of cleaved caspase 7 and 9 did not change, indicating that ATIC can not only affect the level of autophagy but also affect the caspase 3related apoptosis pathway.
ATIC inhibits autophagy through the AMPK/mTOR pathway After knocking down ATIC in Huh7 cells, intracellular autophagosomes (shown by the red arrow) were observed under a transmission electron microscope. The number of intracellular autophagosomes after siRNA treatment was signi cantly increased, which was consistent with the changes in LC3 levels (Fig. 6A). In the transplanted tumor mouse model, the growth rate of tumor cells with ATIC knockdown was signi cantly slower than that of control cells, and the tumor volume was signi cantly smaller than that of the control group (Fig. 6B). The results of immunohistochemistry showed that the expression of ATIC protein was lower in tumors treated with shRNA than in control tumors (Fig. 6C). The lung metastasis ability of tumor cells was signi cantly reduced after ATIC knockdown (Fig. 6D). After transfection of pcDNA/ATIC, the proliferation and migration of tumor cells were signi cantly enhanced ( Fig. 6E and F). The mTOR pathway is closely related to autophagy, so we tested the AMPK/mTOR pathway after knocking down ATIC. After knocking down ATIC, the phosphorylation of AMPKα increased, while the phosphorylation of mTOR decreased (Fig. 6G). The AMPK/mTOR pathway inhibitor compound C (C-C) was used for rescue experiments. After adding C-C, the increase in p-AMPKα levels and the decrease in p-mTOR levels caused by siRNA treatment were reversed. Similarly, the use of AMPK blockers reversed the changes in LC3 protein levels caused by ATIC knockdown (Fig. 6H). Furthermore, C-C also prevented the decrease in cell proliferation and migration caused by knocking down ATIC (Fig. 6I-K).
When siRNA and C-C were used to treat cells at the same time, the number of autophagosomes observed by transmission electron microscopy was signi cantly reduced compared with that seen for siRNAtreated cells (Fig. 6L). The apoptosis induced by ATIC knockdown was also reversed by C-C (Fig. 6M).

Discussion
Given the heterogeneous nature of HCC, it is di cult to carry out satisfactory risk assessment and clinical management for patients with liver cancer [12,13]. Therefore, we want to establish a prognostic model to evaluate the risk and survival of patients to provide guidance for clinical treatment.
Autophagy is regulated by autophagy-related genes, and the abnormal expression of autophagy-related genes affects the occurrence and development of tumors by modulating autophagy [14]. Due to the relationship between autophagy-related gene expression and tumor progression, establishing a liver cancer prognostic model based on autophagy-related genes is reasonable. The prognostic model constructed in this study can predict whether a patient is at high-risk of poor prognosis and the long-term survival rate of the patient. In the GO analysis and KEGG analysis results, autophagy-related genes were closely related to mitosis and the cell cycle. Consistent with this nding, ATIC is known to affect the cell cycle and promote cell proliferation.
It has been reported that autophagy-related genes can inhibit tumor growth by affecting immune cell in ltration, but there are still few studies on the relationship between autophagy and immune in ltration [15]. The immune in ltration analysis results suggest that the autophagy-related genes in the model are related to the tumor immune microenvironment, the activity of immune cells, and immune function, but the mechanisms underlying these relationships are not yet fully clear. HDAC1 has been reported to affect antigen presentation and immune activation. Therefore, the relationship between autophagy-related genes such as ATIC, PRKCD and SPNS1 and immune function deserves further study [16]. We grouped patients according to the expression of ATIC and analyzed the difference in immune cell in ltration between the two groups. It was found that CD8 + T cells, macrophages, immune checkpoints and IFN immune responses were different between the two groups. The main immune cell related to immunotherapy is CD8 + T cells [17,18]; thus, we speculate that ATIC is related to the tumor immune microenvironment and immunotherapy. There are many mechanisms by which autophagy acts on immunity, including abnormal autophagy interfering with the survival and activity of T cells or autophagic activation promoting and inhibiting the secretion of cytokines [19,20]. These ndings suggest that the induction or suppression of autophagy combined with immunotherapy may be a prospective treatment strategy [21].
We screened small molecule drugs with therapeutic potential based on changes in the gene levels of patients in the high-risk group, but in vivo and in vitro experimental veri cation was not performed. The changes in the gene pro le of cells treated with these small molecule drugs were exactly the opposite of the gene pro le of patients in the high-risk group. This suggested that these small molecule drugs have a potential therapeutic effect on liver cancer. Among them, doxorubicin had been widely used in the clinical treatment of liver cancer [22]. Icariin was reported to inhibit the proliferation of tumors, but it had not been studied in liver cancer [23,24]. As for the other molecules, we will conduct research on drug treatment effects and drug resistance in the future.
ATIC is related to tumor cell proliferation, rheumatoid arthritis and the e cacy of radiotherapy [25,26]. There has been no previous report on the correlation between ATIC and autophagy and tumor immune microenvironment. In this study, we have proved through in vivo and in vitro experiments that ATIC can promote the progression of liver cancer through the AMPK/mTOR pathway. And, the study found that our model is related to the patient's prognosis and immune cell in ltration, and it has been veri ed by external databases, which further proves the reliable prognostic ability of our model. Previous studies have shown that autophagy is related to immune cell in ltration and tumor immune tolerance [27,28]. We found that ATIC, as an autophagy-related gene, is related to immune cell in ltration through bioinformatics analysis, but the speci c mechanism is still unclear. We selected ATIC for veri cation in this study, but in subsequent studies, we will also verify the relationship between several other genes and the immune microenvironment and conduct in-depth research on the mechanism of ATIC affecting immune responses.
In the analysis of correlation between ATIC expression and clinicopathologic characteristics of liver cancer patients, we found that the expression of ATIC is related to lymph node invasion and patient prognosis. The TNM stage of the ATIC high expression group was higher, but there was no statistical difference. We think it may be because our sample size is too small, and we will include more patients in the future.
The prognostic model we have established can effectively predict the prognostic risk of patients and provide guidance for clinical decision-making. For example, after a patient has undergone liver cancer resection, are adjuvant treatments such as transcatheter arterial chemoembolization (TACE) necessary [29]? If the patient is a high-risk patient based on the model score, other adjuvant treatments are also recommended even if the tumor is a well-differentiated liver cancer.
At the same time, we have screened out many genes with research value through bioinformatics analysis.
For example, SPNS1 is closely related to the prognosis of liver cancer patients according to the result of multivariate COX regression analysis, but there has not been any report on its mechanism to promote tumor occurrence and tumor development. Therefore, we hope that our analysis results can provide ideas for future research on the role of autophagy in tumors.
The constructed prognostic model can be used to predict the survival rate of patients, and the predictive performance has been veri ed in an external ICGC dataset to support its validity. That the key gene ATIC affects tumor progression by modulating autophagy has been veri ed in vitro and in vivo, and the effect of ATIC on autophagy in liver cancer has been clari ed here for the rst time.

Competing interests
The authors declare that the research was conducted in the absence of any commercial or nancial relationships that could be construed as a potential con ict of interest.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.