Development of a prognostic model for anoikis and identifies hub genes in hepatocellular carcinoma

Considering the high fatality of hepatocellular carcinoma (HCC), current prognostic systems are insufficient to accurately forecast HCC patients' outcomes. In our study, nine anoikis‑related genes (PTRH2, ITGAV, ANXA5, BIRC5, BDNF, BSG, DAP3, SKP2, and EGF) were determined to establish a risk scoring model using LASSO regression, which could be validated in ICGC dataset. Kaplan–Meier curves and time-dependent receiver operating characteristic (ROC) curve analysis confirmed the risk score possessed an accurate predictive value for the prognosis of HCC patients. The high-risk group showed a higher infiltration of aDCs, macrophages, T-follicular helper cells, and Th2 cells. Besides, PD-L1 was significantly higher in the high-risk group compared to the low-risk group. Several anoikis‑related genes, such as ANX5, ITGAV, BDNF and SKP2, were associated with drug sensitivity in HCC. Finally, we identified BIRC5 and SKP2 as hub genes among the nine model genes using WGCNA analysis. BIRC5 and SKP2 were over-expressed in HCC tissues, and their over-expression was associated with poor prognosis, no matter in our cohort by immunohistochemical staining or in the TCGA cohort by mRNA-Seq. In our cohort, BIRC5 expression was highly associated with the T stage, pathologic stage, histologic grade and AFP of HCC patients. In general, our anoikis-related risk model can enhance the ability to predict the survival outcomes of HCC patients and provide a feasible therapeutic strategy for immunotherapy and drug resistance in HCC. BIRC5 and SKP2 are hub genes of anoikis‑related genes in HCC.

under physiological or pathological conditions.Besides, anoikis can limit cancer progression by preventing the dissemination of cancer cells to distant organs 8 .Studies have shown that the induction of the anoikis occurs through the interaction of the intrinsic pathway (mitochondria approach) and extrinsic pathway (death receptor approach).It has been reported cancer cells develop anoikis resistance by escaping the death signals and activating the pro-survival signals 9,10 .Meanwhile, acquisition of anoikis resistance is a hallmark of cancer cell invasiveness, metastasis, treatment resistance, and relapse 11 .In addition, researchers have identified that different players, including cell adhesion molecules, growth factors, hypoxia, stemness, autophagy, long non-coding RNAs, and signaling pathways, are associated with anoikis resistance, allowing cell survival in the cellular environment [12][13][14][15] .However, the precise mechanisms involved in anoikis resistance are not fully understood.Therefore, the early investigation into the mechanisms of anoikis resistance may greatly help to improve the prognostic assessment of cancer and further efficacious treatments for cancer.
Currently, enormous efforts have shown that many anoikis-related genes are associated with the proliferation of several tumors, including gastric cancer, colorectal cancer, lung cancer, and pancreatic cancer [16][17][18] .For example, GDH1 provides signals of anoikis resistance and tumor proliferation in LKB1-deficient lung cancer.Also, GDH1 knockdown can not only improve cell sensitivity to anoikis induction but also reduce tumor metastasis 19 .In high-grade serous ovarian carcinoma, CBX2 has been verified to drive anoikis resistance and tumor metastasis 20 .TMPRSS4, a cell surface-anchored serine protease, is reported to promote anoikis resistance, enhancing survival of circulating tumor cells and leading to early metastasis in prostate cancer 21 .Therefore, anoikis-related genes are considered as promising diagnostic markers and potential therapeutic targets.However, the construction of a prognostic risk model based on anoikis-related genes remains unclear in HCC.Hence, we exploit the relationship between anoikis-related genes and HCC clinical information to forecast the outcomes of HCC patients.
Given the limited effects of standard treatments, increasing attention is focused on immunotherapy in recent years 22,23 .Studies highlighted immune cells producing cytokines and components of the ECM as response to liver injury, which promotes the survival of cancer stem cells 24 .Currently, immunotherapy is considered as an additional treatment and witness substantial progress 25,26 .Thus, exploring the underlying mechanism between the immune system and cancer to enhance the immunotherapy of HCC is urgently required.
Herein, multiple mRNA expression databases were utilized to develop and validate a risk scoring model related to anoikis for HCC.Moreover, we explored the potential relevance between the prognostic risk score and tumor stemness, immunity as well as drug sensitivity.Finally, we identified BIRC5 and SKP2 as hub genes in HCC.They are over-expressed in HCC tissues, and were associated with poor prognosis, no matter in the TCGA cohort or our TMA cohort.All of these might not only promote the capability to forecasting the survival outcomes of HCC patients, but also provide prominent reference for molecular mechanisms of immunotherapy and drug resistance in HCC.

Methods
Data collection process.The mRNA expression matrix of 374 liver hepatocellular carcinoma samples and 50 adjacent normal samples, and the corresponding clinicopathological information were downloaded from The Cancer Genome Atlas (TCGA) database (https:// portal.gdc.cancer.gov/).Besides, we also downloaded the mRNA expression data and the clinical information of another 231 HCC patients from the International Cancer Genome Consortium (ICGC) database (https:// dcc.icgc.org/).From the GeneCards (https:// www.genec ards.org/), 82 anoikis-related genes were extracted by screening with a relevance score > 1.8.

Identification of prognostic anoikis-related DEGs and functional analysis.
To obtain differentially expressed genes (DEGs), we used the "limma" package of R language to conduct differential analysis on the expression of the above 82 anoikis-related genes in normal and tumor samples from the TCGA database.Thirtyeight anoikis-related DEGs were obtained, and the screening criteria were as follows: |log (fold change)| > 1 and adjusted p < 0.05 27 .Through univariate Cox regression analysis, 26 anoikis-related DEGs were associated with the prognosis of HCC (prognostic anoikis-related DEGs, p < 0.05).
STRING Database (Version 11.0, https:// STRIN GDb.org/) was used to map the protein-protein interaction (PPI) network of these 26 prognostic anoikis-related DEGs.The Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses were performed based on the aforementioned prognostic anoikis-related DEGs 28 .

Establishment and verification of anoikis risk scoring model. The least absolute contraction and
selection operator (LASSO) method was used to punish the factor data.The risk score was calculated according to the standardized expression level of each gene and its corresponding regression coefficient.Then an anoikis risk scoring model was established with nine genes.According to the median risk score, all cases were divided into the high-risk scoring group and the low-risk scoring group.The "stats" and "Rtsne" packets of R language were used for principal component analysis (PCA) and t-distribution random neighbor embedding (t-SNE) to detect the internal correlation between the two groups.Kaplan-Meier (K-M) curve was used to evaluate the survival difference between the two groups.The reliability of the risk model was analyzed by using the "timedependent receiver operating characteristic (ROC)" R package.The statistical significance of all analyses was set as p < 0.05.The prognostic accuracy of the anoikis risk score was validated in the ICGC cohort.Univariate and multivariate Cox analyses were used to analyze the independent prognostic role of age, gender, grade, stage, and the anoikis risk score in the TCGA cohort and the ICGC cohort.
Patients and tissue samples.This study was approved by the Human Ethics Committee of the Second Hospital of Dalian Medical University to collect tissue samples.We collected paraffin-embedded liver cancer tissues and adjacent non-tumor tissues from 88 HCC patients who underwent surgical treatment at the Second Hospital of Dalian Medical University from 2017 to 2019.These HCC patients ranged from 35 to 77 years old, including 19 females and 69 males.All patients signed an informed consent form as required by the ethics committee.Complete clinical data and follow-up information were available for all patients.
Immunohistochemistry. We first produced tissue microarrays (TMAs) of these 88 HCC tissues and matched adjacent non-cancerous tissues.Three μm sections were obtained from paraffin-embedded tissues for BIRC5 and SKP2 immunohistochemical studies according to standard procedures.Briefly, a brief protein hydrolysis digestion was performed using IHC enzyme antigen recovery agent and peroxidase blocking was performed using 3% H 2 O 2 .Then sections were incubated overnight at 4 °C with BIRC5 (Survivin) polyclonal antibody (1:200, Bioss, China) or SKP2 polyclonal antibody (1:100, Bioss, China).After PBS wash to remove the unbound primary antibody, the secondary antibody (enzyme-labeled goat anti-rabbit IgG, dilution, 1:1000) was coupled with peroxidase at incubate for 30 min at 37 °C.Sections were stained with 3,3′diaminobenzidine (DAB) substrate and washed with water.Finally, sections were counterstained with hematoxylin for 20 s, then washed with hydrochloric acid and ammonia, dehydrated, cleared, and cover-stained.

Evaluation of staining.
Stained sections were scanned using an Aperio GT450 section scanner (Leica Biosystems Imaging, USA), and markers were pathologically evaluated in the Aperio ImageScope image viewer according to the percentage of positive cells and staining intensity.Protein expression was scored by two independent and experienced pathologists.The degree of staining was defined as the percentage of hepatocytes or paracancerous tissue in the positively stained area as a percentage of the entire tissue area, scored on a scale of 1-4, 1 (0-25%); 2 (26-50); 3 (51-75%); 4 (76-100%).Staining intensity was defined as a score of 0-3, 0 (negative); 1 (weakly positive); 2 (positive); 3 (strongly positive).The percentages and intensity scores were summed to give each patient a final immunoreactivity score (IRS).Low expression of BIRC5 was judged when the IRS was 1 and high expression was judged when the score was 2-3.Low expression of SKP2 was judged when the IRS was 1-3 and high expression was judged when the IRS was 4-6.
Statistical analysis.SPSS 26.0 software was used for statistical analysis, and the chi-square test was used to analyze the relationship between protein/mRNA expression and clinicopathological characteristics of HCC patients.Univariate and multivariate Cox regression analyses were used to analyze the factors influencing survival in HCC patients.The Kaplan-Meier method was used for the detection of survival analysis and the logrank test was used to assess the differences.The relative risk of death was expressed as an adjusted hazard ratio (HR) with a 95% confidence interval (CI).If p < 0.05 all tests were considered statistically significant.
Institutional review board statement.The study was conducted in accordance with the Declaration of Helsinki, and The Ethics Committee of the Second Affiliated Hospital of Dalian Medical University approved this research.

Screening and identification of prognostic anoikis-related DEGs in HCC. The general strategy
for the present study is shown in Supplementary Fig. S1.We screened anoikis-related genes whose relevance score > 1.8 from the GeneCards and obtained 82 genes (Supplementary Table S1).Subsequently, a total of 40 genes related to HCC prognosis were obtained by univariate Cox regression analysis in the TCGA cohort (Fig. 1A).Among the above 82 genes, 38 genes were DEGs between tumor tissues and normal tissues.We intersected them and finally got 26 prognostic anoikis-related genes (Fig. 1B).Then, we performed the heatmap analysis on the expression of 26 overlapping genes in tumor tissues and normal tissues.In general, the 26 anoikis-related genes were overexpressed in tumor tissues compared with normal tissues (Fig. 1C).Next, we used the PPI network to find the intrinsic relationship between the 26 genes.In Fig. 1D Construction of the prognostic model based on anoikis-related genes in the TCGA cohort.Subsequently, we developed the risk score model based on 26 prognostic anoikis-related genes (Supplementary Fig. S3).The Lasso coefficient map of 26 candidate gene expression levels was shown in Fig. 1E.The 26 candidate genes were subjected to Lasso regression analysis for dimensionality reduction, and cross-validation was used to obtain the optimal penalty value (λ) in Lasso regression (Fig. 1F).An anoikis prognostic risk score model was established according to the regression coefficient, and the risk score was calculated using the following formula: risk score = 0.093040 × PTRH2 + 0.014600 × ITGAV + 0.071630 × ANXA5 + 0.099028 × BIRC5 + 0.145664 × BD NF + 0.176856 × BSG + 0.156909 × DAP3 + 0.130661 × SKP2 + 0.041248 × EGF.Based on the median risk score calculated by the risk score formula, HCC patients in the TCGA cohort were divided into high-risk or low-risk groups (Fig. 2A).The results of PCA and t-SNE analysis showed that HCC patients in different risk groups were well-separated (Fig. 2B,C).In addition, we further investigated the relationship between survival status, survival time and the risk score.As shown in Fig. 2D, HCC patients presented higher mortality and shorter survival time in the high-risk group than that in the low-risk group.The prognostic significance of this risk model in HCC was further identified utilizing K-M analysis.In Fig. 2E, the OS of the high-risk group was greatly lower than that of the low-risk group (p = 2.338e−05).We also used the time-ROC analysis curve to evaluate the sensitivity and specificity of the prognostic risk score model.The area under the receiver operating characteristic curve (AUC) was 0.778 at 1 year, 0.713 at 2 years, and 0.697 at 3 years (Fig. 2F).

External validation of anoikis prognostic risk model in the ICGC cohort.
To test the reliability of the prognostic risk model established in the TCGA cohort, we selected HCC patients from the ICGC cohort and also divided them into the low-risk group or high-risk group by the median risk score calculated in the TCGA cohort (Fig. 3A).Similarly, PCA and t-SNE analysis exhibited excellent separation between the two groups (Fig. 3B,C).HCC patients in the high-risk group also have a shorter survival time compared with the low-risk group (Fig. 3D).Besides, the K-M survival curve showed the OS of the high-risk group was vastly shorter compared with that in the low-risk group (p = 3.712e−08, Fig. 3E).At the same time, the time-ROC curve revealed the strong predictive power of our model in the ICGC cohort, and the AUC predictive values for 1-year, 2-years, and 3-years survival rates were 0.733, 0.700, and 0.722 respectively (Fig. 3F).These results confirmed that our prognostic risk model was reliable and exportable.
Independent prognostic value of the risk model.To determine whether the risk score could be an independent prognostic indicator, we used the univariate and multivariate Cox regression analysis to evaluate it and other clinical characteristics (age, gender, pathological TNM stage and grade).According to the results of univariate regression analysis, the risk score and pathological TNM stage were significantly correlated with OS of HCC patients in the TCGA cohort (all p values < 0.001, Fig. 4A).The risk score, pathological TNM stage, and gender were significantly associated with OS of HCC patients in the ICGC cohort (p = 0.031 in gender, and other p values < 0.001, Fig. 4B).Furthermore, we conducted the multivariate Cox regression analysis to correct the influencing factors, and the risk score was an independent predictor of OS (TCGA cohort: HR = 3.300, 95% CI = 2.154-5.056,p < 0.001; ICGC cohort: HR = 3.461, 95% CI = 1.572-7.622,p = 0.002, Fig. 4C,D).In particular, the HR of our risk score was higher than that of the pathological TNM stage.
In addition, we compared the risk score of HCC patients in different clinical characteristic groups.The score was equivalently insignificant between the different age and gender groups (TCGA cohort: p = 0.150 and 0.071, Fig. 4E; ICGC cohort: p = 0.87 and 0.63, Fig. 4F).In the TCGA cohort, the risk score of grade 3-4 patients was higher than that of grade 1-2 (p = 4.2e−09, Fig. 4G).The risk score of pathological TNM III-IV was also found to be higher than that of stage I-II (p = 0.0027, Fig. 4H).In the ICGC cohort, the risk score of stage III-IV was also found to be higher than that of stage I-II (p = 0.00031, Fig. 4I).These results suggested the high risk score of our model was connected with high-level (grade 3-4) and advanced stage HCC patients.

Relationship between immune infiltration and anoikis risk score of HCC.
To explore the potential clinical relevance between immune status and anoikis, it is worth conducting ssGSEA to calculate the enrichment scores of various immune cells, related functions, or pathways in the TCGA and ICGC cohorts.As shown in Fig. 5A and B, aDCs, macrophages, T-follicular helper cells, and Th2 cells all showed high infiltration in the high-risk groups in the TCGA and ICGC cohorts (p < 0.05).Regarding immune-related pathways, MHC class I was significantly upregulated in the high-risk group, while type I IFN response and type II IFN response were downregulated in the TCGA cohort and the ICGC cohort (all adjusted p < 0.05, Fig. 5C,D) and the.In addition,   www.nature.com/scientificreports/expression-based mRNAs represent stem cell expression of the transcriptome.The risk score had no relationship with DNAss, but positively correlated with RNAss (Supplementary Fig. S4A,B).Besides, immune and stromal scores were used to evaluate tumor microenvironment (TME).There was a negative correlation between the risk score and stromal score, but no relationship between the risk score and immune score (Supplementary Fig. S4C,D).The expression level of PD-L1 was significantly higher in the high-risk group compared with the low-risk group (p = 0.00048, Supplementary Fig. S4E), and it was positively correlated with the risk score (R = 0.18, p = 0.00047, Supplementary Fig. S4F).
Relationship between the expression of anoikis-related model genes and drug sensitivity.Anoikis-related genes play crucial roles in regulating drug resistance.We used the CellMiner database to explore the relationship between nine model genes and drug sensitivity.The top 16 gene-drug pairs ranked by Pearson correlation coefficient were shown in Fig. 6, which includes Dexrazoxane, LDK-378, Ibrutinib, Erlotinib, Oxaliplatin, Dasatinib, Lomusti, Tamoxifen, Tyrothricin, Nelarabine, Bleomycin, Ifosfamide, Palbociclib and Lenvatinib, respectively.We found that a total of 4 genes were associated with drug sensitivity, namely ANX5, ITGAV, BDNF, and SKP2.SKP2 was positively correlated with the sensitivity of chemotherapy, while ANXA5 was negatively correlated with the sensitivity of targeted drugs.Interestingly, ANXA5 was negatively associated with seven drug sensitivity tests, including chemotherapy drugs and targeted drugs.
WGCNA identified hub genes in our anoikis model.To find the most critical genes among the 9 model genes for further study, we performed WGCNA analysis on the TCGA database.Based on the scale-free topological model and average connectivity of the WGCNA analysis, the soft threshold was set to 2 (Supplementary Fig. S5A).The clustered gene dendrogram showed HCC patient co-expression modules, with each leaf representing a gene and branches on the tree representing co-expression modules.(Supplementary Fig. S5B).Eleven correlations were determined by WGCNA analysis between modules and normal samples and HCC, with the pink module having the highest positive correlation in HCC (cor = 0.4, p = 3e−25, Fig. 7A).We identified two hub genes (BIRC5 and SKP2) by intersecting the pink module genes with our model genes through Venn diagrams (Fig. 7B).
Based on the TCGA database, we explored the expression differences of these 2 hub genes between normal tissues and HCC tissues.The results revealed that both BIRC5 and SKP2 were significantly highly expressed in HCC tissues compared to normal tissues (Fig. 7C,D).K-M curves showed that both BIRC5 and SKP2 high www.nature.com/scientificreports/expression groups were associated with poor prognosis in patients with HCC compared to low expression groups, both in terms of OS (Fig. 7E,F) and progression-free interval (PFI, Fig. 7G,H).
Correlation between the expression levels of hub genes and clinicopathological characteristics in the TCGA cohort.The correlation between BIRC5 and SKP2 mRNA expression and clinicopathological characteristics of HCC was assessed based on the TCGA database.Supplementary Table S2 showed that BIRC5 mRNA expression levels were significantly correlated with T-stage, pathologic stage, histologic grade, and Alpha-fetoprotein (AFP).In addition, Supplementary Table S3 displayed that SKP2 mRNA expression levels were significantly correlated with age, histologic grade, and AFP.
Our study further investigated the effects of BIRC5 and SKP2 mRNA expression on OS and PFI in HCC by univariate Cox and multivariate Cox regression analysis.Univariate Cox regression showed that the T stage, M stage, pathologic stage, BIRC5 mRNA expression and SKP2 mRNA expression were correlated with OS in HCC (Supplementary Table S4A).Multivariate Cox regression showed that BIRC5 mRNA expression was an independent factor (Supplementary Table S4B).As for PFI, univariate Cox regression suggested that the T stage, M stage, pathologic stage, vascular invasion, BIRC5 mRNA expression and SKP2 mRNA expression were correlated with it in HCC (Supplementary Table S5A).BIRC5 mRNA expression was an independent factor influencing PFI in HCC (Supplementary Table S5B).
We also analyzed the expression of these two factors in several other cancers, including lung adenocarcinoma (LUAD), breast invasive carcinoma (BRCA), as well as most other digestive tumors, including cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), pancreatic adenocarcinoma (PAAD) and rectum adenocarcinoma (READ), to verify the specificity of their prognostic value in HCC.The median expression levels in these cancers themselves were used as cutoff.K-M curves showed that only BIRC5 high expression group were associated with poor OS in PAAD, and BIRC5 has no prognostic value in other cancers.SKP2 has no prognostic value in all of the above cancers (Supplementary Fig. S6).So overall, the value of these two factors for poor prognosis in HCC is somewhat specific.

Immunohistochemical staining of BIRC5 and SKP2 and correlation with clinicopathological characteristics in our TMA cohort.
To further experimentally verify the expression levels of BIRC5 and SKP2, 88 HCC and adjacent non-cancerous tissue samples were collected in this study, and the expression levels of BIRC5 and SKP2 were detected by immunohistochemical staining (Fig. 8A,F).BIRC5 was mainly expressed in the nucleus of HCC cells.The IRS of BIRC5 protein was significantly higher in HCC than in adjacent noncancerous tissues (Fig. 8B-E).In addition, SKP2 was mainly expressed in the cytoplasm of HCC cells.SKP2 protein levels were also significantly higher in HCC tissues than in adjacent non-cancerous tissues (Fig. 8G-K).www.nature.com/scientificreports/ We collected detailed information on the clinical characteristics of the TMA cohort to further validate the correlation between BIRC5 or SKP2 expression and clinicopathological characteristics of HCC.Table 1 showed that BIRC5 expression levels were significantly correlated with T-stage, pathologic stage, histologic grade, and AFP.For SKP2, its expression levels were only significantly correlated with age and smoking history (Supplementary Table S6).
In order to further investigate the potential differential expressions of BIRC5 or SKP2 in liver tissues affected by chronic liver diseases such as hepatitis and cirrhosis, as compared to normal liver tissues, we re-classified our original liver control tissues based on the presence or absence of hepatitis and cirrhosis.We first conducted an analysis using Pearson's Chi-squared test and Fisher's exact test.The results revealed that the expression level of BIRC5 was found to be independent of the history of cirrhosis and hepatitis among HCC patients (Table 2A).On the other hand, the expression level of SKP2 showed a significant correlation with the cirrhosis history of HCC patients (p = 0.042), but not with their hepatitis history (Table 2B).However, when examining the IRS scores, no significant difference was observed in the expression of BIRC5 or SKP2 between tissues with or without cirrhosis, as well as between tissues with or without hepatitis (Supplementary Fig. S7).These findings suggest that although there is a higher prevalence of SKP2 positivity among patients with cirrhosis, there are still individuals without cirrhosis who exhibit high SKP2 expression.Consequently, there is no significant difference in the mean values between these two groups.Therefore, further investigation is required to elucidate the role of SKP2 in cirrhosis.
Prognostic value of BIRC5 and SKP2 in our TMA cohort.The K-M survival curves of the TMA cohort revealed that high expression of BIRC5 was associated with poor prognosis in HCC (Fig. 8L, PFI, p < 0.05; Fig. 8M, OS, p < 0.05).Similarly, PFI (p < 0.05, Fig. 8N) and OS (p < 0.01, Fig. 8O) were significantly lower in SKP2-high expressing HCC compared to SKP2-low expressing HCC.These findings are consistent with the results of our bioinformatics analysis.
Univariate Cox regression showed that T-stage, pathologic stage, histologic grade, BIRC5 expression, and SKP2 expression affected OS in the HCC TMA cohort (Table 3A).Multivariate Cox regression showed that pathologic stage and histologic grade were independent factors affecting OS (the above clinicopathological characteristics and BIRC5 or SKP2 expression, Table 3B).As for PFI, T-stage, pathologic stage, histologic grade, BIRC5 expression, and SKP2 expression were correlated to PFI by univariate Cox regression (Supplementary Table S7A).Multivariate Cox regression showed that the pathologic stage was an independent factor influencing PFI (Supplementary Table S7B).

Discussion
Given the high heterogeneity observed in HCC patients, it is significant to develop a new prognostic system to predict the outcomes of the patients more accurately.It has been proved that anoikis can block cancer development and resistance in anoikis leads to cancer progression and metastasis.So we obstructed and verified a prognostic risk model based on anoikis-related genes in our study.In addition, we identified BIRC5 and SKP2 as hub genes and preliminarily validated them in the TCGA cohort and our own TMA cohort.
The TNM system, widely applied in the majority of substantial tumors, fails to account for the degree of liver dysfunction and performance status of patients 29 .Although the BCLC system may provide more prognostic information, it still lacks universal applicability in clinical practices 6 .Under the situation that prognostic systems for HCC patients have certain shortages, several signatures based on hub genes have been established in recent research 30,31 .Here, our model has higher predictive power in contrast to previous studies, for our model was successfully tested in both internal and external databases.Particularly, the time-ROC curve shows that AUC predictive values for the 1-year survival rate are 0.778 in internal validation and 0.733 in external validation, respectively.We noted the AUCs of the prognostic model based on DNA methylation-driven genes for 1 year, 2 years, and 3 years were 0.650, 0.696, and 0.664 32 .Hence, our model is probably more robust and effective.We found that the high-risk score was correlated with adverse clinical characteristics, such as low tumor differentiation and advanced stage, which contribute to poor survival outcomes of high-risk patients.Besides, univariate and multivariate Cox regression analyses revealed that our prognostic signature is capable of serving as an independent prognostic factor for patients with HCC.Furthermore, We not only investigated the role of anoikis in the immune microenvironment but also explored the relationship between anoikis-related DEGs expression and drug sensitivity.
In the present study, we identified a risk model consisting of the nine anoikis-related genes (PTRH2, ITGAV, ANXA5, BIRC5, BDNF, BSG, DAP3, SKP2, and EGF) for predicting the prognosis of HCC.Recently, Chen et al. constructed a prognostic anoikis model in HCC using part of the TCGA database as a training dataset 33 .Since we used higher criteria for screening the genes of anoikis (relevance score > 1.8) and the whole TCGA dataset as a training set, only BSG is the same gene among our 9 model genes and the above model.Therefore, our model and their model are completely different.More importantly, we have identified hub genes, BIRC5 and SKP2 in HCC anoikis model using the WGCNA method, which is an innovative finding.Furthermore, we validated the overexpression and prognostic value of these two hub genes in liver cancer through immunohistochemical staining in our cohort, providing further confirmation.
Multiple studies have reported that most of the nine anoikis-related genes are related to cancer.IGTAV, involved in cell-cell and ECM adhesions, serves as an independent prognostic factor in esophageal adenocarcinoma (EAC) without neoadjuvant therapy 34 .In colon cancer cells, activation of TP53 suppresses ITGAV expression, resulting in cell survival 35 .Weiler et al. confirmed that overexpression of IGTAV is associated with poor clinical outcomes in patients with HCC 36 .Cilengitide is an ITGAV antagonist, and Mas-Moruno et al. identified that cilengitide reduces angiogenesis and inhibits breast cancer bone metastasis in vivo preliminary results 37,38 .Annexin A5 (ANXA5) was reported to participate in HCC pathogenesis via integrin and MEK-ERK pathways 39 .Kang et al. suggested that the administration of ANXA5 alleviates the immunosuppressive www.nature.com/scientificreports/properties of TME generated by chemotherapy and improves the efficacy of anti-tumor therapy.So the ANXA5 administration following chemotherapy could be considered as a promising immune checkpoint inhibitor for cancer treatment 40 .Basigin (BSG), also known as cluster of differentiation 147 (CD147), is overexpressed and considered to be a prognostic marker in various human tumors [41][42][43][44][45] .It has been observed that CD147 increases cathepsin B expression by activating β-catenin signaling, which mediates CD147-induced invasive phenotype in HCC 46 .Studies have shown that BSG-directed monoclonal antibodies excite clinical success in HCC treatment 47 .Additionally, the overexpression of brain-derived neurotrophic factor (BDNF) and tropomyosinrelated kinase B (TrkB) receptors is associated with lower OS in HCC 48 .It was reported that the BDNF/TrkB axis can promote several important biological processes in several types of cancer, such as cell proliferation, tumor immunosuppression, drug resistance, resistance to anoikis, activation of PI3K/AKT and JAK/STAT3 pathways, EMT, and angiogenesis 49 .Several studies have suggested that the downregulation of BDNF suppresses cell viability and metastasis in multiple solid tumors 50,51 .Therefore, the BDNF/TrkB axis can be considered as a potential therapeutic target for several types of cancer.DAP3 was reported to function as a positive mediator of anoikis.Han et al. found that the expression of DAP3 is significantly associated with the malignant properties of HCC cells.DAP3 suppresses adenosine-to-inosine (A-to-I) RNA epitome to further promote cancer progression 52 .It was observed that EGF is highly expressed in HCC and positively relates to higher tumor grade.In addition, EGF promotes the migration and aggressiveness of HCC cells mainly via the induction of fibronectin (FN) in vitro 53 .
To identify pivotal genes, we utilized the WGCNA algorithm to screen key genes and acquired the overlapping genes among these key genes and our model genes.We obtained the two hub genes: BIRC5 and SKP2.It is necessary to further investigate the potential mechanisms of the two hub genes.
Survivin (BIRC5), a member of the inhibitor of the apoptosis family, can promote cell division and tumor proliferation 54 .Hori et al. found that BIRC can predominantly suppress anoikis 55 .Xu et al. suggested that BIRC5 was upregulated in HCC and associated with a poor prognosis 56 .Co-suppression of OCT4 and BIRC5 can inhibit cancer proliferation by inducing cancer cell apoptosis and cell cycle arrest in HCC 57 .Also, the downregulation of BIRC5 leads to the inhibition of migration and invasion in penile cancer 58 .Thereby, BIRC5 should be regarded as a potential target for anti-cancer drugs and prognosis prediction 59 .
S-phase kinase-associated protein 2 (SKP2), also known as p45, is involved in cancer progression, migration, and invasion 60 .High expression of SKP2 was widely observed in various human tumors, indicating the close correlation between tumor progression and SKP2 [61][62][63] .It is reported that the knockdown of the endogenous SKP2 expression by RNA interference reduces cell proliferation, blocks the cell cycle at the G1 phase, and increases apoptosis in HCC 64 .Delogu et al. reported that SKP2 cooperates with N-Ras and AKT proto-oncogenes to promote carcinogenesis in HCC 65 .Currently, SKP2 is not only an independent prognostic factor for HCC but also a new target for anti-tumor drugs and gene therapy 52,66 .
We further preliminarily verified the two hub genes based on the TCGA database and HCC tissues we collected.We not only analyzed the mRNA expression of the two hub genes in the TCGA cohort, but also performed immunohistochemical staining to explore the protein expression in our TMA cohort.In the TCGA cohort, we observed the upregulation of the two hub genes in HCC patients, accompanied by poor survival outcomes.The expression level of BIRC5 mRNA was associated with the T stage, pathologic stage, histologic grade and AFP.SKP2 mRNA expression level was closely related to age, histologic grade and AFP.In our TMA cohort, all of the HCC patients have indications for surgery, so the percentage of early-stage patients is more than the TCGA cohort (pathologic stage III & IV is 25.7% in the TCGA cohort, and it is 15.9% in our cohort).Therefore, the OS of HCC patients in our TMA cohort seems relatively longer than that in the TCGA cohort.Besides, HCC tissues and adjacent non-cancerous tissues in our cohort are paired, from the same patients.Even though there are some differences between our patients and the TCGA cohort patients, BIRC5 and SKP2 were still both overexpressed in HCC tissues compared with tumor-adjacent tissues in our TMA cohort.Moreover, high expression of BIRC5 or SKP2 protein also indicated poor prognosis in our TMA cohort.As previously mentioned, Fang et al. suggested that silencing of AFP may reduce intracellular BIRC5 mRNA expression to promote apoptosis in HepG2 cells 67 .In lung adenocarcinoma, the high level of BIRC5 was related to advanced AJCC stage, T and N stages 68 .Finally, Taken together, the results above indicated that the expression of the hub genes has a strong correlation with the prognosis and clinical features in HCC, no matter mRNA or protein, which further validated the stability and accuracy of our risk model.
We further employed functional analysis and found that the prognostic model was enriched in anoikis regulation and tumor-related pathways, including PI3K-Akt, ErbB, and P53 pathways.Additionally, our model was significantly associated with virus infection, such as HPV, Human immunodeficiency virus (HIV), Hepatitis B virus (HBV), and Hepatitis C virus (HCV).As revealed by Kakavandi and his colleagues, viral infection promotes resistance to anoikis and provides the opportunity for cancer cells to spread to distant organs, which is in line with our research results [69][70][71] .
Currently, immunotherapy opens a new chapter for cancer treatment.Investigation of the potential mechanism between immune microenvironment and tumor progression is needed.Studies revealed that macrophages have two extreme functional phenotypes, including M1 and M2 macrophages.Tumor-associated macrophages (TAMs) usually acquire an M2-like phenotype and play vital roles in the TME 72 .TAMs produce several cytokines to suppress other immune cells.Besides, TAMs could actively enhance the remodeling of the tumor stroma and facilitate the metastatic process of cancer cells by producing several proteolytic enzymes [73][74][75] .In addition, Li et al. found that the high level of Th2 cells is significantly correlated with poor prognosis in HCC 76 .In our study, the high-risk group has higher infiltration of macrophages and Th2 cells, which means macrophages and Th2 cells are risky immune cells for HCC.We also discovered that PD-L1 was dramatically upregulated in the high-risk group, pointing out that PD-L1 inhibitors could be more effective in HCC patients in the high-risk group.Besides, it is important to consider the contributions of BIRC5 and SKP2 to the tumor microenvironment.BIRC5 is a potential biomarker and inducer of myeloid-derived suppressor cell (MDSC) infiltration in HCC, leading to T cell Vol:.(1234567890 www.nature.com/scientificreports/rejection or dysfunction of the tumor immune microenvironment, ultimately reducing the response to immune checkpoint inhibitors (ICIs) 77 .Cytotoxic T lymphocytes (CTLs) are the primary lymphocytes responsible for killing cancer cells in liver cancer 78 .They primarily restrict tumor growth by inducing G1 phase cell cycle arrest, a mechanism that may involve downregulation of SKP2 79 .Additionally, SKP2 knockdown has been found to significantly increase the expression of immunoinfiltration-related genes in osteosarcoma (OS) mouse models, suggesting that SKP2 may mediate immune rejection in the tumor microenvironment 80 .
Resistance of cancer cells to chemotherapy remains a big challenge in cancer treatment.Resistance to cisplatin treatment is a major clinical issue in HCC 81 .Studies have demonstrated that suppressing BIRC5 can enhance the cytotoxicity of chemotherapy agents like cisplatin 82 .Asechi et al. discovered that BIRC5 expression, mediated by PI3K, enables rat hepatoma cell lines to resist cisplatin-induced apoptosis 83 .Sorafenib, the firstline chemotherapy agent for advanced HCC, sensitizes resistant HCC cells to radiation-induced apoptosis by inhibiting STAT3 phosphorylation and reducing BIRC5 expression 84 .Furthermore, the combination of Sorafenib (SF) and BIRC5 shRNA has shown a synergistic effect in reversing multidrug resistance (MDR).Regarding SKP2, studies have reported that both lonafarnib and troglitazone induce G1-S phase block by down-regulating SKP2.Conversely, ectopic overexpression of SKP2 leads to resistance to troglitazone, resulting a decrease in the number of G1-phase blocked cells 85,86 .Hepatitis B virus X has been found to downregulate SHIP, promoting HCC metastasis and chemical resistance by inducing SKP2 expression 87 .Oxaliplatin (OXA) is one of the most common chemotherapy drugs for patients with HCC 88 .In HCC, cyanidin significantly increases OXA sensitivity and inhibits the EMT induced by OXA via PI3K/Akt signaling 89 .In addition, PKI-587, the dual PI3K/mTOR inhibitor, increases the chemosensitivity of OXA in HCC 90 .In our study, the high level of ANXA5 and ITGAV is negatively related to the low sensitivity to OXA.The underlying mechanism between low sensitivity to OXA and ANXA5 and ITGAV needs further exploration.Currently, Lenvatinib is considered a new first-line targeted therapy for patients with advanced HCC.In our study, high expression of ITGAV is positively related to the high sensitivity of lenvatinib, which indicates that HCC patients with high expression of ITGAV could benefit from lenvatinib.In summary, the combination of specific target inhibitors and chemotherapy drugs appears to be a promising option for the treatment of HCC.
However, a few limitations of this study should be taken into consideration.First, the construction and validation of our prognostic model are based on retrospective databases, so our results should be further confirmed by prospective clinical research.Second, further experimental studies of the nine anoikis-related genes in HCC progression are needed.

Conclusions
Our study established a novel risk signature based on nine anoikis-related genes, which perform a good capability in predicting the prognosis of HCC patients.In addition, our study reflected that risk score was closely related to the immune microenvironment and analyzed the relationship between anoikis-related genes and drug resistance.Finally, we identified BIRC5 and SKP2 as hub genes among the nine model genes, and verified them at the mRNA and protein levels.The results confirmed that BIRC5 and SKP2 should be considered as potential prognostic predictors in HCC.Altogether, our study can not only provide potential therapeutic targets for HCC but also bring new insights into drug resistance and ameliorate the effect of HCC treatment.

Figure 1 .
Figure 1.Identify prognostic anoikis-related DEGs and construct the prognostic model in HCC in the TCGA cohort.(A) P-value and hazard ratio (HR) of 40 anoikis genes related to HCC prognosis by univariate Cox regression analysis.(B) Venn diagram showed anoikis-related DEGs and prognostic anoikis-related genes.(C) Heatmap showed the expression of 26 overlapping genes between tumor tissues (red) and normal tissues (blue).(D) Correlation network of 26 prognostic anoikis-related DEGs.The red lines indicated positive correlations.(E) LASSO coefficient profiles of the expression of anoikis-related DEGs.(F) Penalty value obtained through LASSO cross-validation.An optimal log λ value was indicated by the vertical black line in the figure.

Figure 2 .
Figure 2. Validation of anoikis prognostic risk model in the TCGA cohort.(A-D) Risk score (A) distribution, (B) PCA plot, (C) t-SNE, and (D) OS status analysis of the TCGA cohort.(E) Kaplan-Meier survival curves were used to compare OS between the low-risk and the high-risk groups.(F) ROC curves verified the prognostic performance of the risk model.

Figure 3 .
Figure 3. External Validation of anoikis prognostic risk model in the ICGC cohort.(A-D) Risk score (A) distribution, (B) PCA plot, (C) t-SNE, and (D) OS status analysis of the ICGC cohort.(E) Kaplan-Meier survival curves were used to compare OS between the low-risk and the high-risk groups.(F) ROC curves verified the prognostic performance of the risk model.

Figure 4 .
Figure 4. Univariate and multivariate Cox regression analyses for the anoikis risk score and their levels between different clinicopathologic characteristics.(A,B) Univariate Cox regression analysis of OS in (A) the TCGA cohort and (B) the ICGC cohort.(C,D) Multivariate Cox regression analysis of OS in (C) the TCGA cohort and the (D) ICGC cohort.(E-I) Anoikis risk score levels of age (≤ 65 and > 65), gender (female and male), groups (group 1-2 and group 3-4), and pathological TNM stages (stage I-II and stage III-IV) in the TCGA cohort (E, G, and H) and the ICGC cohort (F and I).

Figure 5 .
Figure 5. Relationship between immune infiltration and anoikis risk score of HCC.(A-D) Box plots presented 16 remarkable immune cells and 13 immune-related functions between two groups in the TCGA cohort and the ICGC cohort.The adjusted p-value was shown as ns, which is not significant; *p < 0.05; **p < 0.01; and ***P < 0.001.(E) Comparison of risk scores in four immune cell subtypes.

Figure 6 .
Figure 6.Scatter plots for the correlation between nine model genes expression and drug sensitivity (top 16 pairs ranked by p value).

Figure 7 .
Figure 7. WGCNA identified hub genes in our anoikis model.(A) Heatmap of the correlation of gene modules with normal tissues and HCC tissues.(B) Venn diagram analysis showed that the overlap of WGCNA analysis and LASSO model led to two hub genes: BIRC5 and SKP2.(C,D) Differential expression of (C) BIRC5 and (D) SKP2 in normal tissues and HCC tissues based on TCGA database.(E,F) K-M curves showing the difference of OS between the (E) BIRC5/(F) SKP2 high and low expression groups.(G,H) K-M curves showing the difference of PFI between the (G) BIRC5/(H) SKP2 high and low expression groups.

Figure 8 .
Figure 8. Experimental validation of immunohistochemical staining of BIRC5 and SKP2 in oure TMA cohort and their prognostic value.(A) Immunohistochemical staining of BIRC5 in the HCC TMA cohort.(B) The IRS of BIRC5 in HCC and adjacent non-cancerous tissues.****p < 0.0001; Data are presented as the mean ± standard error of mean.(C-E) The intensities of BIRC5 immunostaining were (C) positive, (D) weakly positive and (E) negative, no strongly staining was observed in HCC tissues.Left image, magnification ×100; right image, magnification ×400; the red squares indicate the area shown in higher magnification.(F) Immunohistochemical staining of SKP2 in the HCC TMA cohort.(G) The IRS of SKP2 in HCC and adjacent non-cancerous tissues.****p < 0.0001; Data are presented as the mean ± standard error of mean.(H-K) The intensities of SKP2 immunostaining were (H) strongly positive, (I) positive, (J) weakly positive and (K) negative.Left image, magnification ×100; right image, magnification ×400; the red squares indicate the area shown in higher magnification.(L,M) Kaplan-Meier survival curves to compare (L) OS and (M) PFI between the high-BIRC5 expression and low-BIRC5 expression patients.(N,O) Kaplan-Meier survival curves to compare (N) OS and (O) PFI between the high-SKP2 expression and low-SKP2 expression patients. https://doi.org/10.1038/s41598-023-41139-9

Table 1 .
Correlation of BIRC5 expression levels with clinicopathological characteristics of HCC in the TMA cohort.

Table 2 .
Correlation of BIRC5 and SKP2 expression levels of liver control tissues with cirrhosis or hepatitis in the TMA cohort. A.

Table 3 .
Univariate Cox regression analysis and multivariate regression analyses of OS in the TMA cohort. A.