Prognostic association of starvation-induced genes in head and neck cancer


 Autophagy-related genes (ARGs) have been implicated in the initiation and progression of malignant tumor promotion. To investigate the dynamics of expression of genes, including ARGs, head and neck squamous cell carcinoma (HNSCC) cells were placed under serum-free conditions to induce growth retardation and autophagy, and these starved cells were subjected to transcriptome analysis. Among the 21 starvation-induced genes (SIGs) located in the autophagy, cell proliferation, and survival signaling pathways, we identified SIGs that showed prominent up-regulation or down-regulation. These included AGR2, BST2, CALR, CD22, DDIT3, FOXA2, HSPA5, PIWIL4, PYCR1, SGK3, and TRIB3. The Cancer Genome Atlas (TCGA) database of HNSCC patients was used to examine the expression of up-regulated genes, and CALR, HSPA5, and TRIB3 were highly expressed relative to solid normal tissue in cancer and the survival rate was found to be reduced in patients with high expression. Protein-protein interaction analysis demonstrated the formation of a dense network of these genes. Cox regression analysis revealed that high expression of CALR, HSPA5, and TRIB3 was associated with poor prognosis in patients with TCGA-HNSCC. Therefore, these SIGs up-regulated under serum starvation may be molecular prognostic markers in HNSCC patients.


Introduction
Head and neck cancer is the sixth most common malignancy in the world, 90-95% of which is squamous cell carcinoma (SCC). Over 60% of patients already have advanced cancer at the time of their rst visit, with an estimated 5-year survival rate of 40-50% [1][2][3][4] . Surgery, radiation, chemotherapy, and targeted therapies are used to treat head and neck squamous cell carcinoma (HNSCC) 5,6 . In recent years, immunotherapy with antibodies that target the immune checkpoint pathway has been introduced and has shown long-term effects on cisplatin-resistant cancer, distant metastases, and recurrence of poor prognosis 7,8 . However, valid cases are limited to 18-25% of advanced HNSCCs 9 . Long-term immunological side effects are also a problem 10 . Effective indications for these therapies need to be developed and new therapies need to be developed.
A major advance in recent HNSCC research is the aggregation of extensive genetic analysis results of HNSCC [11][12][13][14][15] . A typical HNSCC database is the Cancer Genome Atlas (TCGA), published in 2015. Recent technological advances have enabled TCGA and other large-scale genomics studies to determine the broader landscape and frequency of chromosomal alterations, mutations, and expressed genes that contribute to HNSCC pathogenesis, prognosis, and resistance to therapy [11][12][13][14][15] . The TCGA-HNSCC database may be used to screen for differentially expressed genes (DEGs) in cancer and normal tissue transcriptome studies in HNSCC patients 16,17 . Furthermore, using the TCGA database, many studies on the deviation of genes and signaling pathways involved in carcinogenesis and prognosis are being conducted. This includes studies on hypoxia-immune signature 18 , cancer-associated alternative splicing event-related genes 19 , the miRNA-30 family [20][21][22] , and the KEAP1-NRF2-CUL3 axis 23 . Consequently, promising biomarker genes for the prognosis of HNSCC patients have been proposed.

Cells
The human HNSCC cell line SAS cells were obtained from the Japanese Collection of Research Bioresources (Tokyo, Japan). Cells were cultured in Dulbecco's modi ed Eagle's medium supplemented with 10% fetal bovine serum (FBS) at 37°C in a humidi ed atmosphere with 5% CO 2 .

Transmission electron microscopy (TEM)
TEM was performed to observe SAS cells in serum-starved condition for 24 hours. Serum-starved cells were washed with PBS, xed in 2.5% glutaraldehyde in phosphate buffer, post-xed in 2% osmium tetroxide, dehydrated in graded ethanol, and then embedded in epoxy resin. Ultrathin sections were stained with 2% uranyl acetate and observed using a JEM-1200 EX microscope (JEOL, Tokyo, Japan).

Cell proliferation assay and migration assay
For MTT assay, starved or control SAS cells were incubated with 3-(4,5-Dimethyl-2-thiazolyl)-2,5-diphenyl-2H-tetrazolium bromide reagent (DOJINDO, Osaka, Japan) for 2 hours at 37°C. At the end of each experiment, the medium was removed and 100 μL solution of 4% HCl 1N in isopropanol was added to immediately dissolve the formazan crystals, and absorbance at 570 nm was recorded. For the migration assay, SAS cells were cultured in DMEM with 10% FCS until confluent. The cell layers were scratched using a plastic tip, as previously described 40 . The cells were further incubated in DMEM with/without FBS for 6 hours. The closure rate of each scratched area was measured using ImageJ software, as previously described 40  and all subsequent analyses were executed using this software. We set the lower limit as replacing positive numbers less than 10 with 10, and 0 counts with 8. Then, we calculated the log2 ratio against the geometric mean of the two control samples. We ltered out genes if their counts were always less than 15, or if their log2 ratios were between -0.5 and 0.5 in all samples; a total of 6,363 genes remained after ltering. We extracted candidate DEGs by a 2-fold criterion.

Analyzing TCGA-HNSC RNA-Seq data
We obtained and analyzed the RNA-Seq count data of TCGA-HNSC from the GDC Data Portal 42 with the Subio platform. The work ow of TCGA RNA-Seq was the same as that applied to our RNA-seq data except for the thresholds. The lower limit for positive counts was 50, for 0 counts was 32, and the lter on counts was 50, and that on log2 ratios was between -1 and 1. In addition, the log ratios were taken against the average of solid normal tissue samples. For each of the 21 selected genes, we divided the primary tumor samples into two groups, those with count values higher or lower than the median, to compare the survival time with the Kaplan-Meier method.

Pathways analysis and Protein and Protein interaction
The molecular pathways of the 21 selected genes were analyzed for gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) server. GO enrichment was carried out over three primary levels: cellular components (CC), biological processes (BP), and molecular functions (MF). Based on the STRING online database (https://string-db.org/), we used these genes to establish a protein-protein interaction (PPI) network. Then, the most signi cant modules in the PPI networks were visualized.

Statistical analyses
Statistical analyses were performed using the Student's t-test with Microsoft Excel (Microsoft, Redmond, WA, USA). Results were expressed as the mean ± SD. Differences were considered signi cant at P<0.05.
For the survival analysis shown in Table 3, the hazard ratio (HR) relative to the indicated reference (ref) value, its 95% con dence interval (CI), and P-value (those of < 0.05 are indicated in bold) for the Cox hazard model are shown. The HR and its 95% CI were calculated by Cox regression analysis after proper evaluation of the assumptions of the Cox regression models with the use of the survival package.

Effects of serum starvation on the biological activity of head and neck cancer SAS cells
We first examined the effect of serum starvation on cell morphology and proliferation of SAS cells. The morphology of SAS cells was not affected in the absence of serum for 24 hours (Fig. S1). When examined by MTT assay, cell proliferation decreased in a time-dependent manner in the absence of serum (Fig. 1A). As for migration assays, SAS cell migration was considerably diminished in serum-free conditions ( Fig. 1B, C). Morphological changes at the hyper ne structure level due to serum starvation were further investigated through TEM. Most SAS cells contained intact mitochondria that were distributed throughout the uniform cytoplasm in the presence of serum. When the cells were cultured in the absence of serum for 24 hours, the intracellular space was occupied by large vascular structures. There were autophagosomes and/or autolysosomes containing degraded mitochondria and dense structures (Fig. 1D).
3.3. RNA-sequencing of serum-starved SAS cells to detect the expression of genes related to cell growth, cell death, cell migration, cell proliferation, cell cycle, and cell adhesion Since 24-hour serum starvation reduced the ability of SAS cells to proliferate and migrate, RNA sequencing data obtained after starvation were further referred to as cell growth, cell death, cell migration, cell proliferation, cell cycle, and cell adhesion. Analyzed using 6 keywords ( Table 2). After 24 hours of serum starvation, more than two-fold up-regulation was observed for 425 genes. The top 5 genes were determined for each keyword . This included HSPA1A, OSGIN1, UCN, BST2, and SGK3 for cell growth,  LOC728739, UCN, NPAS2, AGR2, and PYCR1 for cell death, BST2, ADGRA2, CALR, SGK3, and HSPA5 for  cell migration, SGK3, IKZF3, SPTA1, MIR17HG, and CD22 for cell proliferation, ERN1, DDIT3, BEX2, CALR,   and PIWIL4 for cell cycle, and AMIGO1, TNXB, TNC, FOXA2, and CD22 for cell adhesion (Table 2).

Expression of serum starvation-induced genes in TCGA-HNSCC patients
From the 70 genes altered in HNSCC cells by 24-hour serum de ciency, the top two genes showing signi cant expression changes were selected for each of the 7 keywords, including autophagy, cell proliferation, cell death, cell migration, cell proliferation, cell cycle, and cell adhesion. Of the 28 genes selected, 6 were associated with replication. Moreover, the microRNA MIR17HG was excluded. Therefore, we nally focused on 21 genes. Of these, 11 were up-regulated genes and 10 were down-regulated genes.
When the expression of these SIGs was examined in TCGA-HNSCC patients, 9 of the 11 up-regulated genes were also up-regulated in the primary tumor compared to solid normal tissue. Signi cant expression differences were observed in BST2, CALR, DDIT3, HSPA5, and TRIB3 ( Fig. 2A). On the other hand, 6 out of the 10 down-regulated genes had reduced expression in tumors compared to solid normal tissue, with signi cant differences observed in the ATOH8 and CCL2 genes (Fig. 2B).

Function and PPI analysis of SIGs
GO and KEGG enrichment pathway analyses were performed to investigate the biological properties and potential signaling pathways of the 21 selected genes. Using GO enrichment analysis, enriched terms were ATF6-mediated unfolded protein response, PERK-mediated unfolded protein response, negative regulation of sequence-speci c DNA-binding transcription factor activity, and negative regulation of transcription. These GO terms were associated with several important biological processes including DNA-templated gene expression response to endoplasmic reticulum stress, endoplasmic reticulum stress response, and positive regulation of cell cycle arrest (Fig. 3A). KEGG analysis showed that the prognostic genes were signi cantly enriched in pathways of transcriptional misregulation in cancer and protein processing in the endoplasmic reticulum (Fig. 3B). In PPI network analysis, CALR, HSPA5, DDIT3, and TRIB3 formed a close interaction network among the up-regulated genes (Fig. 3C). The association with AGR2, FOXA2, CD22, and SGK3 was weak, and PIWIL4 and PYCR1 were not associated with other upregulated genes.
3.6. Prognostic signi cance of SIGs in TCGA-HNSCC patients TCGA-HNSCC patients were divided into two groups based on the expression of SIGs. The expression levels of patients in the high expression group were higher than the median, and the remaining patients were classi ed in the low expression group 43,44 . The difference in survival time determined by the Kaplan-Meier method was examined using the generalized Wilcoxon test and the long rank test (Fig. 4). Among the up-regulated SIGs, high expression of CALR ( (Table 3). Multivariate analysis showed that the combination of two genes (CALR-High and HSPA5-High) (P = 0.022) and three genes (P = 0.027) did not make a clear difference in correlation.

Discussion
Autophagy has been suggested to be a biological marker for estimating the prognosis of cancer patients.
In a previous HNSCC bioinformatics study, Li et al. 29 identi ed a novel autophagy-related signature consisting of three hub genes, MAP1LC3B, FADD, and LAMP1, that may provide promising biomarker genes for the treatment and prognosis of HNSCC. Similarly, Jin et al. 33 determined 35 genes for HNSCC and identi ed ITGA3, CDKN2A, FADD, NKX2-3, BAK1, CXCR4, and HSPB8 as prognostic ARGs. Ren et al. 34 also reported 13 ARGs as genes that predict prognosis. In the present study, HNSCC cells were cultured under serum starvation, which can e ciently induce autophagy, and RNA sequencing was used to examine the expression of ARGs.
FBS is commonly used as a supplement to animal cell culture medium 45 . Additionally, FBS consists of several compositions such as macromolecules, carrier proteins for lipoid substances and trace elements, attachment and spreading factors, low molecular weight nutrients, hormones, and growth factors 45 . Among them, growth factors were reported to in uence cell proliferation, migration, survival, and morphogenesis 46 . Under serum starvation, SAS cells, a high-risk HPV-negative HNSCC cell line 47 , showed no signi cant changes in cell morphology after 24 hours, but cell growth and migration capacity was suppressed. Electron micrographs revealed the presence of autophagosomes and mitochondrial phagocytosis, suggesting that autophagy was induced in this situation. After 24-hour starvation, mRNA sequencing of SAS cells detected 12 up-regulated ARGs (ATP6V0A2, ATP6V1B1, ATP6V1C2, DDIT3, ERN1, NHLRC1, NUPR1, PIM2, TMEM150A, TRIB3, WIPI1, and XBP1) and 13 down-regulated ARGs (BNIP3, BNIP3L, C10orf10, DAPK2, GAPDH, HMOX1, MEFV, PLK2, RRAGD, SESN3, SRPX, S100A8, and S100A9). These genes differed from the ARGs previously reported to predict the prognosis of HNSCC patients 29,33,34 . This starvation-induced approach may be bene cial in extrapolating ARGs that have not been previously identi ed as differentially expressed genes. In addition, as with ARGs, we also found aberrant expression of genes related to cell growth, cell death, cell migration, cell proliferation, cell cycle, and cell migration. Finally, 21 SIGs that showed signi cant up-regulation or down-regulation were selected. Comparing how these genes were expressed in normal and cancer tissues in TCGA-HNSCC patients, we found 11 genes that were more strongly expressed in cancer cells and 10 genes that were down-regulated in cancer tissues. Among them, BST2, CALR, DDIT3, HSPA5, and TRIB3 were signi cantly up-regulated in cancer tissues.
GO and KEGG analyses revealed the involvement of ATF6-mediated unfolded protein responses and PERK-mediated unfolded protein responses mainly in the nucleus, and the ability of SIGs to bind glycoproteins and ubiquitin protein ligases. In addition, networking between CALR, HSPA5, DDIT3, and TRITB3 was demonstrated by PPI analysis. Consistent with the PPI analysis results, when TCGA-HNSCC patients were divided into high-expression and low-expression groups, and then analyzed by the Kaplan-Meier method, CALR, FOXA2, HSPA5, and TRIB3 were found to be correlated with reduced survival. Cox regression analysis further con rmed that three SIGs (CALR, HSPA5, and TRIB3), sex, M-stage, and Nstage were associated with survival in HNSCC patients. This suggests that CALR, HSPA5, and TRIB3 are predictors of poor prognosis. Since the combination of two genes (CALR-High and HSPA5-High) and three genes did not make a clear difference in correlation (Table 3), patients will have a poor prognosis, especially when both CALR and HSPA5 are highly expressed.
Calreticulin, CALR, is a soluble multifunctional protein found in the endoplasmic reticulum (ER) lumen and is involved in calcium homeostasis, transcriptional regulation, immune response, and cellular function 48, 49 . It is expressed at higher levels in many cancerous tissues than in normal tissues. High CALR expression is correlated with both advanced clinical stage and lymph node metastasis [50][51][52] . The translocation of CALR from the ER to the plasma membrane occurs during immunogenic cell death 53 due to oncolytic virus infection of HNSCC cells, leading to enhanced tumor immunity 54 . The unfolding protein response (UPR) is a cellular stress response related with ER stress. One of the proteins involved in this UPR is Heat shock 70 kDa protein 5/glucose-regulated protein (HSPA5/GRP78). HSPA5 is the master regulator of UPR and is associated with tumor progression, tumor size, and poor prognosis [55][56][57][58] . In situations where protein production is required for tumor growth, USPA5 is overactivated to process a high ux of protein passing through the ER, maintaining ER homeostasis. Expression of HSPA5 is induced by glucose starvation 59,60 . It has also been shown that HSPA5 is up-regulated during viral infection, translocated to the cell membrane, and recognized by SARS-CoV-2 spikes 61 .
Tribbles homologue 3, TRIB3, is a member of the mammalian pseudokinase tribble family and is involved in multiple biological processes including the cellular response to glucose de ciency stress and ER stress. Several studies have shown that TRIB3 is elevated in multiple cancer cell lines and primary tumors including colorectal cancer, breast cancer, and lung cancer. In renal cancer, TRIB3 is overexpressed compared to normal tissue and is associated with tumor progression and poor prognosis [62][63][64] . On the other hand, TRIB3 binds to AKT, inhibits its phosphorylation and prevents its activation 65 . In the case of HNSCC, patients with high TRIB3 expression had a better prognosis, and overexpression of TRIB3 reduced AKT phosphorylation and inhibited AKT activation 66 . All of these genes are associated with endoplasmic reticulum stress. Their up-regulation may be a promising biomarker for predicting the prognosis of HNSCC.
There are many past studies where in vitro events and RNA-sequencing data were linked to informatics analysis of HNSCC patients 67,68 . You et al. 67 established radiation-resistant cells by repeated irradiation in vitro and identi ed radioresistant genes using the TCGA-HNSCC database. In the present study, by analyzing genes induced by serum starvation of HNSCC cells, we detected genes that could not be obtained by previous TCGA database analysis and show their usefulness in predicting the prognosis of HNSCC patients. This approach may help to understand the genetic response of cancer cells to ER stress under therapeutic processes such as radiation therapy and chemotherapy.

Conclusions
Up-regulated and down-regulated genes associated with serum starvation using HNSCC cells were identi ed. Expression of HSPA5, TRIB3, and CALR in SAS cells was up-regulated by in vitro serum starvation and up-regulated in TCGA-HNSCC tissue tumors. High expression of these genes was closely associated with reduced survival in patients with TCGA-HNSCC. These SIGs have the potential to be molecular prognostic markers in HNSCC patients.   Expression of serum starvation-induced genes (SIGs) in TCGA-HNSCC patients The expression of 21 genes that showed signi cant changes in expression by serum starvation and their relative expression levels in TCGA-HNSCC patients were determined in primary tumors and solid normal tissues. (A) 11 genes   (AGR2, BST2, CALR, CD22, DDIT3, FOXA2, HSPA5, PIWIL4, PYCR1, SGK3, and TRIB3) have been upregulated more than two-fold. (B) Ten genes (ATOH8, AXIN2, BNIP3, CCL2, CDKN2C, CNTN1, EGLN3, ID2, SERPINB3, and S100A8) were down-regulated by more than 50%.

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