Prognostic association of starvation-induced gene expression 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 in vitro. 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 found to be highly expressed relative to solid normal tissue in cancer and the survival rate was 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.


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 significant 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 significant 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% confidence 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.

Results
Effects of serum starvation on the biological activity of HNSCC cells. We first examined the process of characterizing the response of HNSCC cells to serum starvation on cell proliferation, migration, and morphology. Following serum starvation, SAS cells proliferation and migration were considerably diminished, while cell morphology did not change ( Fig. 1A-D). These findings were consistent with Ca9-22 cells proliferation, migration, and morphology mediated by serum starvation (Fig. S1A-D), and suggest that there is no much of differences between HNSCC cell lines. To determine morphological changes at the hyperfine structure level due to serum starvation, SAS cells were further investigated through TEM. Most SAS cells maintained in the presence of serum contained intact mitochondria that were distributed throughout the uniform cytoplasm. In cells cultured for 24 h in the absence of serum, autophagosomes and/or autolysosomes containing degraded mitochondria and dense structures, characteristic of autophagic cells not present in control cells, were observed (Fig. 1E). Based on the results in experiments on the biological activity of HNSCC cells under serum starvation, we decided to perform RNA-sequencing of SAS cells under serum starvation and extracted genes showing large fluctuations. We also investigated their expression in tumor tissues of TCGA-HNSCC patients and the relationship between their expression and the prognosis of patients. Fig. S2 shows the schedule of these experiments.
RNA-sequencing of serum-starved SAS cells and altered expression of genes related to autophagy, cell growth, cell death, cell migration, cell proliferation, cell cycle, and cell adhesion. SAS cells were cultured in the absence of serum for 2 and 24 h and ARG expression was examined by RNA sequencing (Table 1 and Table S1), then we performed principal components analysis (PCA) and found that the expression profile did not change significantly after 2 h of starvation, but after 24 h of starvation, the expression profile changed (Fig. S3). At first, the altered expression of ARGs due to serum starvation was examined. To be consistent with PCA, a slight change in gene expression was observed 2 h after the onset of starvation. Some genes, such as DDIT3 and ERN1, were down-regulated after 2 h of starvation but increased after 24 h. After 24-h serum starvation, more than two-fold up-regulation was observed for 12 genes (ATP6V0A2,  ATP6V1B1, ATP6V1C2, DDIT3, ERN1, NHLRC1, NUPR1, PIM2, TMEM150A, TRIB3, WIPI1, and XBP1)  (Table 1). On the other hand, down-regulation of 50% or more was observed for 13 genes (BNIP3, BNIP3L, C10orf10, DAPK2, GAPDH, HMOX1, MEFV, PLK2, RRAGD, SESN3, SRPX, S100A8, and S100A9) ( Table 1).
Expression of serum starvation-induced genes in TCGA-HNSCC patients. From the 70 genes altered in HNSCC cells by 24 h serum deficiency, the top two genes showing significant 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 finally 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 www.nature.com/scientificreports/ www.nature.com/scientificreports/ tissue. Significant 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 significant differences observed in the ATOH8 and CCL2 genes (Fig. 2B). A heat map was also created to represent the level of up-or down-regulated expression profiles of 21 genes (Fig. 2C).
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-specific 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 (ER) stress, ER stress response, and positive regulation of cell cycle arrest (Fig. 3A). KEGG analysis showed that the prognostic genes were significantly enriched in pathways of transcriptional misregulation in cancer and protein processing in the endoplasmic reticulum (Fig. 3B). In PPI network analysis, 21 genes were subdivided into 4 clusters (I-IV). In cluster I, up-regulated genes, CALR, HSPA5, DDIT3, and TRIB3, formed a close interaction network (Fig. 3C). PIWIL4 and PYCR1 in cluster IV were not associated with other up-regulated genes.
Prognostic significance of 21 SIGs in TCGA-HNSCC patients. We investigated whether SIGs that the differentially expressed SIGs in tumors and normal tissues of TCGA-HNSCC patients was associated with prognosis. 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 classified 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 expres- www.nature.com/scientificreports/ sion of CALR (Fig. 4C), FOXA2 (Fig. 4F), HSPA5 (Fig. 4G), and TRIB3 (Fig. 4K) was correlated with decreased patient survival. FOXA2 was excluded in subsequent studies because it was not significantly up-regulated in tumors compared to normal solid tissue. Conversely, high expression at AGR2 (Fig. 4A) and PIWIL4 (Fig. 4H) were correlated with significant improvement in overall patient survival. On the other hand, high expression of BST2 (Fig. 4B), CD22 (Fig. 4D), DDIT3 (Fig. 4E), PYCR1 (Fig. 4I), and SGK3 (Fig. 4J) was not associated with patient survival. The Kaplan-Meier method was also applied to down-regulated genes, but there was no association between gene expression and patient survival in TCGA-HNSCC patients (Fig. S4). When the survival curve was recalculated based on the expression of CALR, HSPA5, and TRIB3, the probability of survival in the high-and high group combinations was much lower than in the low-and low group combinations, predicting patient prognosis. It shows the high ability of group combination to do (Fig. 5).
Cox regression analysis of the association of SGIs and classical prognostic factors with survival in the TCGA-HNSCC patients. Expression of CALR, HSPA5, and TRIB3 was correlated with reduced overall survival in patients with TCGA-HNSCC, so these genes were further analyzed. Univariate and multivariate analysis (Cox proportional hazard model) was performed using the three genes and classical risk factors, such as gender, HPV, smoking, age, and TNM stage, as independent variables. In univariate analysis, CALR-High  (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 identified 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 identified 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 efficiently induce autophagy, and RNA sequencing was used to examine the expression of ARGs.    www.nature.com/scientificreports/ 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 influence cell proliferation, migration, survival, and morphogenesis 46 . Under serum starvation, SAS cells, a high-risk HPV-negative HNSCC cell line 47 , showed no significant changes in cell morphology after 24 h, but cell growth and migration capacity were suppressed. Serum starvation showed no significant effect on deforming cell morphology under microscopy. However, electron micrographs revealed the presence of autophagosomes and mitochondrial phagocytosis, being consistent with the features during autophagy of SAS cells 27 . This suggested that autophagy was induced in this serum-deficient situation. After 24-h 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), again supporting the induction of autophagy of SAS cells under serum starvation. These genes differed from the ARGs previously reported to predict the prognosis of HNSCC patients 29,33,34 . This starvation-induced approach may be beneficial in extrapolating ARGs that have not been previously identified as differentially expressed Survival curves for high CALR-high HSPA5 group, high CALR-low HSPA5 group, low CALR-high HSPA5 group, and low CALR-low HSPA5 group. (B) Survival curves for high CALR-high TRIB3 group, high CALR-low TRIB3 group, low CALR-high TRIB3 group, and low CALR-low TRIB3 group. (C) Survival curves for high HSPA5high TRIB3 group, high HSPA5-low TRIB3 group, low HSPA5-high TRIB3 group, and low HSPA5-low TRIB3 group. www.nature.com/scientificreports/ 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 (Fig. S2). Finally, 21 SIGs that showed significant 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 significantly up-regulated in cancer tissues. GO and KEGG analyses revealed the involvement of ATF6-mediated unfolded protein responses and PERKmediated 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 TRIB3 was demonstrated by PPI analysis as a cluster. 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. FOXA2 was excluded because its expression was not significantly increased in tumors compared to normal tissues in TCGA-HNSCC patients. In contrast, some in vitro up-regulated SIGs, such as AGR2 and CD22, showed no significant difference, but survival was inversely proportional to that of CALR, FOXA2, HSPA5, and TRIB3. This may be due to the fact that there was no significant difference in AGR2 and CD22 expression between tumors and normal tissues ( Fig. 2A,B). www.nature.com/scientificreports/ On the other hand, if there is a clear difference in survival, high expression of these genes may be applicable to predict a better prognosis for patients. Recalculation of the survival curve between CALR, HSPA5, and TRIB3 showed that comparing the combination of the two high groups with te combination of the low groups significantly reduced the probability of survival (Fig. 5). Furthermore, cox regression analysis confirmed that three SIGs (CALR, HSPA5, and TRIB3), sex, M-stage, and N-stage 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 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] . CALR has been shown to promote cell motility and enhance resistance to anoikis through STAT3-CTTN-AKT pathway of esophageal SCC 53 . Positive CALR staining was observed in the majority of tumor case (96%) of the oral cavity, whereas the incidence was lower in non-cancerous matching tissue cases (32%). It was also been reported that stable knockdown of CALR in oral cancer cells reduced cell proliferation 50 . 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 [54][55][56][57] . In situations where protein production is required for tumor growth, USPA5 is overactivated to process a high flux of protein passing through the ER, maintaining ER homeostasis. Expression of HSPA5 is induced by glucose starvation 58,59 . Correspondingly, HSPA5 has been reported to be up-regulated in tumors of various organs such as breast, liver, stomach, esophagus, brain, prostate, head and neck, and melanoma, and may be accompanied by aggressive tumor behavior and recurrence 60,61 . A comprehensive proteomic analysis of oral SCCs also showed up-regulation of three members of the HSP family, including HSP90, HSPA5 and HSPA8 62 .
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 deficiency 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 [63][64][65] . In the tongue SCC, both TRIB3 and AKT were highly expressed compared to adjacent non-cancerous tissues, correlating TRIB3 overexpression with tumor pathological T stage, lymph node metastasis, and tumor recurrence. However, when TRIB3 was overexpressed in tongue SCC cells using a viral vector, phosphorylated AKT protein was reduced 66 . All of these genes, CALR, HSPA5, and TRIB3, are associated with ER stress. Their up-regulation may be a promising biomarker for predicting the prognosis of HNSCC.
There are several 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 identified 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 identified. Expression of HSPA5, TRIB3, and CALR in SAS cells was up-regulated by in vitro serum starvation and upregulated 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.