Comprehensive bioinformatics analysis and experimental validation: An anoikis-related gene prognostic model for targeted drug development in head and neck squamous cell carcinoma

We analyzed RNA-sequencing (RNA-seq) and clinical data from head and neck squamous cell carcinoma (HNSCC) patients in The Cancer Genome Atlas (TCGA) Genomic Data Commons (GDC) portal to investigate the prognostic value of anoikis-related genes (ARGs) in HNSCC and develop new targeted drugs. Differentially expressed ARGs were screened using bioinformatics methods; subsequently, a prognostic model including three ARGs (CDKN2A, BIRC5, and PLAU) was constructed. Our results showed that the model-based risk score was a good prognostic indicator, and the potential of the three ARGs in HNSCC prognosis was validated by the TISCH database, the model’s accuracy was validated in two independent cohorts of the Gene Expression Omnibus database. Immune correlation analysis and half-maximal inhibitory concentration were also performed to reveal the different landscapes of TIME between risk groups and to predict immuno- and chemo-therapeutic responses. Potential small-molecule drugs for HNSCC were subsequently predicted using the L1000FWD database. Finally, in vitro experiments were used to verify the database findings. The relative ARG mRNA expression levels in HNSCC and surrounding normal tissues remained consistent with the model results. BIRC5 knockdown inhibited anoikis resistance in WSU-HN6 and CAL-27 cells. Molecular docking, real-time PCR, cell counting kit-8 (CCK-8), plate clone, and flow cytometry analyses showed that small-molecule drugs predicted by the database may target the ARGs in the prognostic model, inhibit HNSCC cells survival rate, and promote anoikis in vitro. Therefore, we constructed a new ARG model for HNSCC patients that can predict prognosis and immune activity and identify a potential small-molecule drug for HNSCC, paving the way for clinically targeting anoikis in HNSCC.


Introduction
Head and neck cancer (HNC) is the sixth most common malignancy and the leading cause of death and decreased life expectancy worldwide. Oral cavity, oropharynx, and nasopharynx tumors are included in HNC. Over 500,000 fatal HNC cases occurred in 2020, and over one million new cases were diagnosed [1]. Epidemiological studies show that the incidence of HNC is rising annually, with a progressive trend toward a higher incidence in younger people. Squamous cell carcinoma accounts for over 90% of HNC [2]. Although tremendous strides have been made in head and neck squamous cell carcinoma (HNSCC) treatment, the 5-year survival rate remains below 50% since most HNSCC patients have advanced disease at diagnosis [3]. Ninety percent of cancer patients die from metastases [4]; lymph node metastases are the most common in HNSCC and are a key prognostic indicator. The 5-year survival rate is considerably poorer for HNSCC patients who experience lymph node metastases [5,6]. It is essential to develop a predictive model for HNSCC metastasis that occurs before lymph node metastasis to improve patient prognosis.
In general, most cells survive by adhering tightly to the extracellular matrix (ECM) or other cells, and cell death occurs once cells are separated from the ECM. This form of cell death was first named "anoikis" in 1994 [7]. Anoikis is a special type of apoptotic cell death with different cell morphology and biochemical markers from those of ferroptosis, necroptosis, and autophagy [8,9]. Anoikis exhibits some features of apoptosis, mainly depends on the cell-matrix interactions, and plays a critical role in cancer metastasis [7]. Tumor metastasis, by contrast, assumes that cancer cells have evolved resistance to anoikis. When cancer cells are liberated from their cell-ECM and cell-cell adhesion states, they survive, disseminate, and metastasize in the circulatory system by resisting anoikis-induced tumor cell death via autocrine or paracrine pathways [10,11].
The significance of anoikis in HNC is becoming clearer as research and genetic testing procedures develop, and anoikis resistance acquisition has been hypothesized to be a critical step in oral cancer metastasis [12]. The signaling pathways associated with anoikis resistance in HNSCC have been the subject of numerous recent investigations [13,14]. For instance, it has been proposed that the epithelial-mesenchymal transition (EMT) plays a significant role in the development of anoikis resistance in cancer cells [13]. Several scientists have investigated the changes in gene expression that occur after HNSCC cells develop anoikis resistance [14]. However, which genes are essential for HNSCC anoikis resistance and how anoikis-related genes (ARGs) affect patient prognosis remains uncertain. In addition, targeting ARGs may help develop effective drugs that can reduce metastasis in HNSCC patients.
We searched for novel prognostic indicators to explore the role of ARGs in prognosis, predicted potential effective drugs for HNSCC, and introduced new methods for the existing relatively mature comprehensive sequence therapy model. We used bioinformatics methods to construct a high-accuracy ARG-based prognostic model. Subsequently, a small-molecule drug was obtained from the L1000FWD database, and its targeting of the prognostic model genes was validated.

Data and clinical sample collection
We obtained RNA sequencing (RNA-seq) and clinical information about HNSCC from The Cancer Genome Atlas (TCGA) Genomic Data Commons (GDC) portal (https:// portal.gdc.cancer.gov/) and the Gene Expression Omnibus (GEO) database (datasets: GSE27020 n = 99 and GSE41613 n = 76), then removed data with missing survival status, unknown survival time, or survival < 30 days. The single cell RNA sequencing (scRNA-seq) of HNSCC patients was obtained from the TISCH (http://tisch.comp-genomics.org/). (dataset: HNSCC_GSE103322, cell number = 5902). GeneCards was applied to obtain (https://www.genecards.org/) a total of 419 genes (relevance score with anoikis > 0.4) (Suppl . Table S1). Three pairs of HNSCC samples were obtained for real-time PCR (RT-PCR) and were preserved in liquid nitrogen and frozen at −80°C. Our study protocol was approved by the Biomedical Ethics Committee of Peking University Stomatological Hospital. The code used in this study can be found on GitHub (https://github.com/ qiulin961028/anoikis-analysis.git), and the study flow chart is shown in Fig. 1.

Bioinformatics analysis
RNA-seq data of HNSCC patients were analyzed by log2 transformation. Differentially expressed ARGs were subsequently screened and clustered using the R software (4.1.2) LIMMA package. Univariate and multivariate COX analyses were performed to identify HNSCC prognosisrelated ARGs for subsequent prognostic model construction.
Construction and validation of the ARG prognostic model After the multivariate COX regression analysis, the risk regression coefficients and expressions of each ARG were combined to establish the risk score formula, which led to the prognostic model, calculated as follows: Coefi Â Xi Coef represents the risk regression coefficient of ARGs, and X represents the ARG expression level. The risk score of each HNSCC patient was measured, and all patients were divided into high-or low-risk groups in accordance with the median risk score. Immediately afterward, the overall survival (OS) outcomes of HNSCC patients were compared between the two risk groups by survival analysis. Receiver operating characteristic (ROC) curves were used to test the accuracy of model predictions. The role of the risk score in predicting prognosis was investigated by univariate and multivariate COX regression analyses.

Small interfering RNA transfection
The cells were cultured in a 6-well plate at a density of 2.0 × 10 5 cells/well, and two groups were established: the control group (si-NC) and the interference group (si-BIRC5). According to the recommended protocol, the working concentration of siRNA was 50 nM, and the Lipo8000 (Beyotime Biotechnology Co., Ltd., Shanghai, China) transfection reagent was used. After 24-72 h of culture, the cells were collected for follow-up experiments. The expression of ARG after knockdown was determined by RT-PCR. The siRNAs were synthesized by Tsingke Biotechnology Co., Ltd., China, with the following sequences: si-BIRC5: CCGCATCTCTACATTCAAGAA.

Molecular docking
Molecular docking was performed to analyze the binding capacity of radicicol and ARGs. The PubChem database (https://pubchem.ncbi.nlm.nih.gov/) was accessed to download the radicicol structure, and the protein structure was obtained from the Protein Data Bank (PDB; http:// www.rcsb.org/). Hydrogen atoms were added, charges were calculated, and charges and bonds on small molecules were adjusted using AutoDockTools (version 1.5.6), and the data were stored in a PDBQT file. Vina was used to calculate the binding energy, and PyMOL (version 4.6.0) was used to visualize the optimal binding model. Statistical analysis R software (version 4.1; PBC, Boston, MA, USA) and GraphPad Prism 7.0 (version 8.0; La Jolla, CA, USA) were used for the statistical analyses. All data are expressed as the mean ± standard deviation. T test and one-way ANOVA were used to compare data between two or more groups that conformed to the normal distribution, respectively. The Wilcoxon rank sum test was used for data that did not conform to a normal distribution, and p < 0.05 was considered to indicate a significant difference.

Prognostic model establishment and verification
A univariate COX analysis was performed on the 50 ARGs, and variables with p < 0.05 in the univariate COX analysis were examined (Fig. 3A) in a multivariate analysis. A total of three ARGs (CDKN2A, BIRC5, and PLAU) were significantly associated with prognosis, leading to the construction of a prognostic model (Fig. 3B). We determined risk scores for all cases with expression levels for risk regression coefficients and ARGs. Risk score = CDKN2A expression × (−0.01388) + BIRC5 expression × 0.013027 + PLAU expression × 0.001909. All patients were ranked according to their risk scores, and the differential expression profiles for three ARGs between the high-and low-risk groups were displayed in a heatmap (Fig. 3C), showing that the CDKN2A expression level was obviously decreased in the high-risk group, whereas the BIRC5 and PLAU expression levels were both significantly increased. Further analysis of patient survival status showed that patients in the high-risk group had worse survival outcomes and a higher likelihood of death than those in the low-risk group (Figs. 3D and 3E). The subsequent survival analysis results confirmed that the overall survival (OS) rate of the high-risk group was lower than that of the low-risk group (p < 0.001) (Fig. 3F). ROC curves were used to demonstrate the predictive performance of the model, with areas under the curve (AUC) of 0.714, 0.703, and 0.704 for 1-, 3-, and 5year OS rates, respectively (Fig. 3G), indicating the good predictability of the prognostic model for HNSCC patients.

Risk scores independently predicted HNSCC prognosis
The relationship of each clinical parameter to the risk score was investigated by clinicopathological analysis. The results were plotted as a heatmap (Fig. 4A), which showed that the risk scores were significantly associated with N stage (p < 0.05), T stage (p < 0.001), tumor stage (p < 0.001), and grade (p < 0.01) but were not relevant to age and sex. The relationship between the risk score and M stage could not be assessed due to the lack of M stage information. The difference in the prognosis predicted by the risk score with other clinicopathological features was assessed by the multivariate ROC curve. As shown in Fig. 4B, the AUC value of the risk score was 0.714, which significantly exceeded that of the additional clinical indicators. The combined ROC curve for T stage, N stage, and risk score was then plotted (Fig. 4C). The AUC value of the combined ROC curve was 0.768, which was significantly higher than that of any other clinical feature. These results suggest that the risk score is a more accurate prognostic predictor for HNSCC patients than other clinical features and a useful complementary means to predict HNSCC patients' prognosis according to the TNM stage. In addition, we performed univariate and multivariate COX regression analyses to verify the independent prognostic significance of the risk score for HNSCC patient OS outcomes. The univariate COX regression analysis revealed significant differences in age, sex, grade, tumor stage, and risk score (Fig. 4D), while the multivariate COX regression analysis showed that the risk score could be an independent prognostic factor for HNSCC patients (Fig. 4E) (HR = 2.134, 95% CI = 1.524-2.988). In conclusion, the superiority of the model for predicting the prognosis of HNSCC patients was demonstrated. The model-based risk score can be used as an independent prognostic factor in HNSCC patients.
Investigation of immunity status and clinical treatment response analyses in high-and low-risk groups We next evaluated the relationship between the risk model and immune infiltration. The immune cell infiltration analysis showed that the low-risk group had higher immune cell infiltration levels, such as B cells, CD8+ T cells, interdigitating dendritic cells (iDCs), macrophages, mast cells, NK cells, plasmacytoid dendritic cells (pDCs), T helper cells, follicular helper T cell (Tfh), T helper 2 cells (Th2), and tumor-infiltrating lymphocytes (TILs) (Fig. 5A). In addition, immune function was compared between the two groups, with the low-risk group demonstrating more prosperous immune function (Fig. 5B). The above results  and immunotherapeutic agents (e.g., Ribociclib) (Fig. 5C), than the high-risk group. All 41 chemicals were shown in Suppl. Fig. S1. Twelve drugs had lower IC 50 in the high-risk group (Suppl. Fig. S2), of which some drugs including Luminespib were associated with immunotherapy (Fig. 5D).  Fig. S3). These results were consistent with the results of TCGA database analysis, and again confirmed that CDKN2A, BIRC5, and PLAU have great potential as prognostic markers in HNSCC patients. Then risk scores were calculated for each patient in the GSE27020 and GSE41613 datasets using the same formula used for the external validation of the prognostic model. All patients in the two datasets were divided into high-and low-risk groups as previously described using the median risk score for each dataset. The survival curves of the two datasets were as expected, with significantly lower OS rates in the high-risk group than in the low-risk group (p < 0.01) (Figs. 6E and 6F). The ROC curves in GSE27020 had AUC values of 0.719, 0.697, and 0.671 for 1-, 3-, and 5-year OS rates, respectively (Fig. 6G). The AUC values in GSE41613 were 0.729, 0.763, and 0.732 for 1-, 3-, and 5-year OS rates, respectively (Fig. 6H). These results again confirmed the ability of the prognostic model to predict prognosis. Unfortunately, the lack of clinical information, such as sex, T stage, and N stage, prevented the comparison of the predictive ability difference between clinical features and risk scores in the dataset. Nevertheless, the results were sufficient and highlighted the outstanding ability of the model to predict prognosis.

L1000FWD analysis to identify potential target drugs
We searched for potential target drugs of HNSCC by uploading up-and down-regulated ARGs to the L1000FWD database and obtaining the top ten candidates. The basic drug information is shown in Suppl. Table S3. We selected the top three smallmolecule drugs (radicicol, dasatinib, and BRD-K85660637) for visualization. The 2D and 3D structures are shown in Fig. 7. Radicicol, the top-ranked fungicide, was selected for verification in a series of in vitro experiments.

Exploration of the expression characteristics of ARGs and their effect on anoikis in vitro
The mRNA expression level of three ARGs was detected in three pairs of matched HNSCC (T), paracancerous normal tissues (N), HOK cells, and two HNSCC cell lines. As shown in Figs. 8A and 8B, the relative mRNA expression levels of CDKN2A, BIRC5, and PLAU were higher in HNSCC tissues than in paracancerous normal tissues, while the relative mRNA expression levels of CDKN2A, BIRC5, and PLAU were also higher in both HNSCC cell lines than in HOK cells. Therefore, the expression level of ARGs was consistent with the results of our model analysis.
Subsequently, we prevented cell adhesion by cultivating HNSCC cells on nonadherent plates to mimic the environment of cancer cells to verify the correlation between the ARG and anoikis. We selected one of the ARGs (BIRC5) for experimental verification and knocked down BIRC5 in WSU-HN6 and CAL-27 cells, then cultured the cell suspension for 24 h. Flow cytometry was used to assess the proportion of cells undergoing anoikis. The proportion of apoptotic cells was increased significantly in the BIRC5knockdown group compared with that in the si-NC group (Figs. 8C and 8D). Combined with previous results, the overexpression of ARGs in HNSCC tissues and cells can significantly promote anoikis resistance in tumor cells.

Radicicol can stably bind to ARGs and inhibit their expression level
The lower the binding energy between the ligand and the receptor, the more stable it is. Therefore, a binding energy ≤ -5.0 kcal/mol was chosen as the screening condition. In this study, the binding energies of the three ARGs and radicicol were much lower than the set conditions (the binding energies of CDKN2A, BIRC5, and PLAU with radicicol were -7.0, -7.0, and -7.8 kcal/mol, respectively), indicating a good binding effect between radicicol and each ARG (Fig. 9A). We measured the expression levels of CDKN2A, BIRC5, and PLAU to dissect the relationship between radicicol and the three ARGs. RT-PCR showed that the relative mRNA expression levels of CDKN2A, BIRC5, and PLAU were decreased after radicicol treatment (Figs. 9B and 9C). Therefore, radicicol can regulate the mRNA expression level of ARGs, suggesting that radicicol may target the ARGs in the prognostic model.  ). Moreover, the number of colonies in the three radicicol groups was reduced compared with that of the control group (Fig. 10C). Radicicol can regulate the cell cycle and promote anoikis Since cell survival rate was inhibited in the radicicol groups, we hypothesized that one or more cell cycle phases might be blocked during this process. Flow cytometry analysis was performed to test this hypothesis. The radicicol groups had a significantly higher proportion of cells in the G 2 /M phase and a significantly lower proportion in the G1 and S phases than that in the control group (p < 0.001) (Figs. 11A and 11B), indicating that the cell cycle was blocked in the G 2 /M phase. As previously mentioned, HNSCC cells were cultured on nonadherent panels to investigate the effect of radicicol on anoikis. Radicicol was cocultured with suspension culture cells for 24 h, and the cells were stained using Annexin V-FITC/PI to analyze the level of anoikis through the intensity of the cell fluorescence signal. Annexin V-FITC-positive and PI-negative cells represented early anoikis, while Annexin V-FITC-positive and PI-positive cells represented late anoikis. The results showed that the percentage of anoikis cells (percentage of early + late anoikis) was significantly increased in the radicicol groups vs. the control group (Figs. 11C and 11D). A dose-dependent effect was demonstrated in all of the above experiments, and the effect of radicicol was more pronounced in CAL-27 cells. These results confirm that radicicol can block the G 2 /M phase and promote anoikis.

Discussion
The prognosis of HNSCC depends on various factors, including age, lifestyle habits, and treatment, with metastases predominating the prognosis of HNSCC [16]. Several models currently use tumor metastasis to predict the prognosis of HNSCC. TNM staging is currently an accepted model for predicting prognosis; however, N staging only examines the number and size of positive lymph nodes, which are affected by the type of neck lymphatic dissection and the total number of lymph nodes removed [17,18]. The lymph node ratio (LNR) can be a valid predictor of prognosis in HNSCC [19]; however, the LNR cannot be used to assess patient prognosis when the number of positive lymph nodes is zero. Therefore, the creation of a prognostic discriminating model for patients with HNSCC based on tumor metastasis capacity is critical for post-operative patient observation and directing clinical medication use.
Anoikis resistance is a crucial molecular mechanism for survival during the metastatic cascade of tumor cells and is one of the hallmarks of tumorigenesis in EMT and a signature trait of tumor stem cells [20]. The development of anoikis resistance leads to an increased potential for tumor cell metastasis, an expansion of cancer stem cell subpopulations, chemoresistance, and a higher likelihood of recurrence, all of which are significantly related to a poor prognosis in HNSCC [21]. ARG-based predictive models play a significant role in determining the prognosis of many malignancies [22][23][24]; thus, anoikis resistance has major therapeutic implications. In this study, we used bioinformatics analysis to build an ARG prognostic model for HNSCC to predict patient prognosis and immune activity. CDKN2A, BIRC5, and PLAU were found to be closely related to patient prognosis, consistent with the results of previous studies. CDKN2A was first identified by Kamb and can encode two proteins, namely, p16INK4a and p14ARF, thus exerting cell cycle regulation through multiple pathways [25]. CDKN2A gene expression abnormalities, primarily deletions, mutations, and abnormal hypermethylation have now been reported in a range of malignancies, including HNSCC [26]. Patients with CDKN2A deletion in HNSCC have a generally poor prognosis and are considerably more likely to experience HNSCC recurrence [27,28]. And targeting CDKN2A and/or the PI3K-AKT-mTOR pathway may be a valuable direction to develop precise therapy for HNSCC [29]. BIRC5 was first isolated and identified in a human gene bank screen by Ambrosini in 1997 and is believed to be the smallest member of the apoptosis suppressor protein family [30]. Moreover, studies have shown that abnormal expression of BIRC5 can be used as a diagnostic marker in HNSCC patients [31] and that BIRC5 is an important predictor of poor prognosis [32]. The PLAU gene is located on human chromosome 10q22.2 and encodes urokinase-type plasminogen activator [33]. PLAU has been shown to be closely related to tumor diagnosis and patient prognosis [34]. Chen demonstrated PLAU may function as an oncogene in HNSCC and regulate the EMT signaling pathway in vitro and in vivo [35]. Li demonstrated that PI3K-Akt pathway might underly the mechanism of PLAU's oncogene role in HNSCC [36]. The elevation of CDKN2A, BIRC5, and PLAU expression seen in our tissue samples was consistent with prior studies; however, more clinical samples are needed to determine the link between these three genes' expression levels and metastasis. In addition, although CDKN2A, BIRC5, and PLAU are associated with the prognosis of HNSCC, their effects on anoikis have not been experimentally verified. We selected BIRC5 to transfect siRNA in vitro to detect its effect on anoikis. According to the flow cytometry results, the inhibition of BIRC5 expression increased tumor cell apoptosis and decreased HNSCC anoikis resistance.
According to current studies, the tumor microenvironment (TME) provides a permissive environment for tumor progression and metastasis [37]. Therefore, considering that the TME can regulate anoikis resistance, we conducted ssGSEA to explore the immune status of high-and low-risk groups. The immune cells (B cells, CD8+ T cells, iDCs, macrophages, mast cells, NK cells, pDCs, T helper cells, Tfh, Th2, and TIL) and immune roles (checkpoint, cytolytic activity, HIL, inflammationpromoting, T cell co-inhibition, T cell co-stimulation, and type I (or II) IFN response) were more active in the low-risk group. These results imply that we can use the risk scores of ARGs to effectively distinguish the immune infiltration status of HNSCC and develop personalized immunotherapy. As for drug sensitivity, the low-risk group was more sensitive to multiple inhibitors of the PI3K/AKT/mTOR pathway, such as OSI-027, which not only enhanced immunotherapeutic effects [38] but also blocked the progression of HNSCC [39]. One study has found that PI3K/AKT/mTOR axis is highly activated in HNSCC, which is related to the proliferation, migration, invasion, and other biological behaviors of tumor cells [29]. As a result, our study can guide clinical treatment and provide evidence for improving future immunotherapy and seeking appropriate target populations.
In this study, differentially expressed ARGs were divided into up-regulated and down-regulated groups, with a view to using the database to develop effective drugs against anoikis. The radicicol ranks first among the obtained small molecule drugs. Radicicol, a macrolide antibiotic isolated from Monosporium bonorden by Delmotte et al. [40], is still in its infancy as a novel Hsp90 inhibitor. Since Hsp90 has regulatory effects on a variety of substrate proteins, inhibition of Hsp90 can regulate many signaling pathways, thereby inhibiting tumor proliferation, metastasis, and other processes [41], which has made clinically feasible Hsp90 small-molecule inhibitors a research focus. Interestingly, in the previous drug sensitivity analyses, we also found an Hsp90 inhibitor had a lower IC 50 in the high-risk group, indicating that Hsp90 inhibitor has great potential to become a novel drug targeting ARGs prognostic model. However, radicicol was discovered only a short time ago and has the disadvantage of low activity in vivo; thus, it has not yet entered the clinical research stage, especially for HNSCC. A recent study demonstrated that radicicol promoted anoikis in glioblastoma by causing endoplasmic reticulum stress and preventing AKT/mTOR/p3S3K phosphorylation activation, thus confirming that radicicol was closely associated with PI3K/AKT/mTOR and anoikis [42]. In addition, we demonstrated in vitro that radicicol modulated mRNA expression levels of the three ARGs included in the prognostic model, and that radicicol inhibited HNSCC cells survival rate and promoted anoikis. In summary, combined with the studies of scholars, we hypothesized that radicicol may promote the anoikis of HNSCC by inhibiting the expressions of CDKN2A, BIRC5, and PLAU, thereby blocking the PI3K/AKT/mTOR signaling pathway. However, the specific molecular mechanism still needs to be confirmed by subsequent experiments.
Our study provides a new perspective for exploring the role of ARGs in the development and metastasis of HNSCC.
We constructed a prognostic model, identified novel drug that promotes anoikis, and validated the results of the database analysis in vitro. However, it must be acknowledged that there are still some limitations in this study. The model was derived from a public database analysis; thus, its applicability and accuracy need to be further explored in clinical HNSCC patients. The effect of radicicol on anoikis in vivo and its molecular mechanisms still need to be further explored. In conclusion, the ARG prognostic model is expected to be a novel biomarker for HNSCC diagnosis and treatment decisions, and ARGs can be considered drug development targets for reducing HNSCC metastasis.
Acknowledgement: The authors would like to thank the TCGA and GEO database for the availability of the data. Availability of Data and Materials: All data used and analyzed in this study are available from the corresponding author on reasonable request.
Ethics Approval: The studies involving human participants were reviewed and approved by the Biomedical Ethics Committee of Peking University Stomatological Hospital. The ethical code is PKUSSIRB-202274058.

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