Integrated analysis of autophagy-related genes(ATGs) reveals the prognostic value of oral squamous cell carcinoma

Background: Increasing evidence demonstrated that autophagy paly a crucial role in initiation and progression of OSCC. The aim of this study was to explore the prognostic value of autophagy-related genes(ATGs) in patients with OSCC. RNA-seq and clinical data were downloaded from TCGA database following extrating ATGs expression profiles. Then, differentially expressed analysis was performed in R software EdgeR package, and the potential biological function of differentially expressed ATGs were explored by GO and KEGG enrichment analysis. Furthermore, a risk score model based on ATGs was constructed to predict the overall survival. Moreover, univariate, multivariate cox regression and survival analysis were used to select autophagy related biomarkers which were identified by RT-qPCR in OSCC cell lines, OSCC tissues and matched normal mucosal tissues. Results: Total of 232 ATGs were extrated and 37 genes were differentially expressed in OSCC. GO and KEGG analysis indicated that these differentially expressed genes were mainly located in autophagosome membrane, and associated with apoptosis, platinum drug resistance, ErbB signaling pathway and TNF signaling pathway. Furthermore, a risk score model including 9 variables was constructed and subsequently identified with univariate, multivariate cox regression, survival analysis and Receiver Operating Characteristic curve(ROC). Moreover, ATG12 and BID were identified as potential autophagy related biomakers. Conclusion: This study successfully constructed a risk model to predict the prognosis of patients with OSCC, and the risk score may be as a independent prognostic biomarker in OSCC. ATG12 and BID were identified as potential biomarkers in tumor diagnosis and treatment of OSCC.


Introduction
Oral squamous cell carcinoma (OSCC) is the most common head and neck squamous cell carcinoma and one of the most common problems in many parts of the world [1]. Smoking, drinking, betel nut and human papillomavirus(HPV) infection are most conmon risk factors [2]. Despite advances in medical equipment and treatment, the overall survival rate of OSCC is still unsatisfied. It mainly owed to the lack of effective biomarkers for early-stage diagosis, recurrence and lymph node metastases [3]. Therefore, it is crucial and pressing to identify the effective biomarkers and therapeutic targets for OSCC to improve the overall survival.
Autophagosome is formed by encapsulating part of the cytoplasm and components such as organelles and proteins. Then, autophagosome is fused with lysosomes to form autophagosome, which degrades the wrapped contents, so as to realize the metabolic needs of the cell itself and the renewal of some organelles. In summary, autophagy is a lysosomal-mediated catabolic complex process which involves the cytoplasmic organelles and proteins to maintain metabolism and homeostasis in cells [4].
Recently, increasing number of evidence showed that autophagy dysregulation played a crucial role in a variety of human malignant tumors, including colorectal cancer [5], renal cell carcinoma [6], nonsmall cell lung cancer [7] and so on. Numerous evidence showed that autophagy may play a significant role in OSCC carcinogenesis. Such as autophagy mediated cell apoptosis to promote tumor progression via the AKT/mTOR pathway in OSCC [8]. In addition, dysregulation of autophagy was relevant to tumorigenesis and prognosis in OSCC [9]. However, the role of whole subsets of ATGs in OSCC has not been explored. Therefore, comprehensive analysis of ATGs may reveal the prognostic value and potential therapeutic targets for OSCC diagnosis and treatment.
The purpose of this study was to analyze the differential expression of ATGs in OSCC and establish a cox regression model to predict the survival of OSCC patients. Furthermore, survival analysis was performed to identify this cox regression model. The present study provides novel insight for better understanding of ATGs in OSCC and useful resource for identifcation of novel biomarkers of OSCC.

Identification of differentially expressed ATGs
316 OSCC patients and 32 normal controls RNA-seq and corresponding clinical data were downloaded from TCGA database. Then, 232 autophagy related genes expression level were extracted from the transcriptome data. Subsequently, differentially expressed analysis was performed in R software EdgeR package(Figure1A,B,C).Finally, 11 upregulated and 26 downregulated autophagy related genes were sorted out With the cut-off criteria |log2FoldChange|>1 and FDR < 0.01.

Functional Analysis of differentially expressed autophagy related genes
GO enrichment analysis indicated that these genes were mainly located in autophagosome membrane, and associated with cytokine receptor binding and regulation of apoptotic signaling pathway(Figure2A,C). In addition, KEGG pathway analyses showed that most of significant autophagy related genes were enriched in apoptosis, platinum drug resistance, ErbB signaling pathway and TNF signaling pathway(Figure2B,D).

Establishment of cox regression model
Through univariate and multivariate cox regression analysis, total of 9 variables including BID, ATG12,

Survival analysis
The prognostic values of the risk score for different clinicopathological parameters including age, gender, T and N in TNM system, grade and stage were further investigated. M classification in TNM system were excluded because of numerous data missing. Survival analysis demonstrated that low risk group had significantly longer overall survival than high risk group in the stratification analysis based on age, gender, T, N, grade and stage(Figure5).

Identification of potential independent prognostic biomakers
Comprehensive bioinformatics indicated that ATG12 and BID were associated with overall survival in OSCC. And univariate and multivariate cox regression showed that ATG12 and BID might be selected as potential independent prognostic biomarkers in our study. Therefore, ATG12 and BID genes expression were validated in OSCC cell lines and tissues by qRT-PCR. Our results revealed that ATG12 and BID were upregulation in OSCC cell lines and tissues than MNTs, which were similar with the results in TCGA database (Figure6A,B). Unfortunately, the correlation between ATG12,BID and clinical parameters showed no significant difference (Table1, Table2).

Discussion
Owing to low 5-year overall survival rates and high recurrence rate [10], it is crucial to explore the effective therapeutic targets to improve the overall survival in OSCC. Recently, increasing number of studies indicated that autophagy may play a significant role in the genesis and development of OSCC [11][12][13]. However, the role of whole subsets of ATGs in OSCC keep unclear. Therefore, our study aimed to analyze the ATGs comprehensively to reveal the potential therapeutic targets of OSCC.
Total of 232 ATGs were enrolled in this study and 37 genes were differentially expressed in OSCC.
According to GO enrichment analysis, differentially expressed ATGs were accociated with apoptosis which was relevant to cancer progression [14,15]. In addition, KEGG enrichment analysis showed that these ATGs were relevant to platinum drug resistance, ErbB signaling pathway and TNF signaling pathway. Among these pathways, platinum drug resistance may play a important role in OSCC treatment strategies [16]. In addition, ErbB signaling pathway was also related to head and neck squamous cell carcinoma (HNSCC) treatment [17].
In the present study, we established a risk score model through bioinformatics analysis of 37 differentially expressed ATGs to predict overall survival of OSCC. Finally, 9 ATGs were selected in regression model and it can accurately predict overall survival. Meanwhile, univariate, multivariate cox regression and survival analysis based on stratification analysis indicated that risk score can be regarded as independent prognostic factor. And it also can distinguish early stage tumors and terminal stage tumors. Therefore, the risk score might hold potential in OSCC early diagnosis, forecast and therapy. ATG12 is the human homolog of a yeast protein involved in autophagy. Mizushima N et al demostated that ATG12 system is well conserved and may function in autophagy also in human cells [18]. Upregulation of ATG12 was correlated with advanced TNM stage in gastric cancer [19] but downregulation of ATG12 induced oncosis in tumor cells [20]. BID was found on chromosome 22q11.21 and encodes for a protein associated with apoptosis, which is heterodimerized with apoptotic activator BAX or negative apoptotic regulator BCL2. associated with poor prognosis of clear cell Renal Cell Carcinoma (ccRCC) patients in a autophagy related manner [21]. Unfortunately, few studies have reported the relationship between ATG12,BID and OSCC. In our study, bioinformatics analysis indicated that ATG12 and BID were potential independent prognostic biomarkers and subsequentially validated in OSCC cell lines and specimans.Therefore, ATG12 and BID would play a role as novel biomarkers in tumor diagnosis and treatment of OSCC. However, the correlation between ATG12,BID and clinical parameters showed no significant difference. The exact role of ATG12 and BID in OSCC need further study. In addtion, BAK1, SHPK1, LAMP1, ATIC were also associated with viarious cancers in a autophagy-dependent manner [22] [23][24] [25]. Ney PA demonstrated that BNIP3 was accociated with mitochondrial autophagy [26]. However, there are fewer studies indicated that other genes in the risk formula play important roles in the development and progression of cancer in a autophagy manner. Numerous researches reported that autophagy played a important role in OSCC [11,27]. And these genes regulated autophagy in a variety of ways ,whether these autophagy related genes were also related to OSCC development and progression in a autophagy manner remains to be further investigated.
We ananlyzed autophagy related genes comprehensively indicating that these differentially expressed genes were relevant to OSCC initiation and progression. Furthermore, a 9 autophagy related gene risk score model was successfully constructed which is positively associated with overall survival in OSCC. And this risk formula might provided potential in OSCC early diagnosis, forecast and therapy. However, there are sevearl limitations in our study. First, the sample numbers in TCGA database are markedly inadequate. Furthermore, it is necessary for further validation in our following research.

Conclusion
This study successfully constructed a risk model to predict the prognosis of patients with OSCC through conprehensively analyzing ATGs, and the risk score may be as a independent prognostic biomarker in OSCC. Furthermore, ATG12 and BID were identified as potential biomarkers in tumor diagnosis and treatment of OSCC.

Data downloading
All transcriptome profiling of OSCC was downloaded from The Cancer Genome Atlas(TCGA) database(https://portal.gdc.cancer.gov/). Owing to its half-baked clinical data, 3 samples were excluded. Finally, 316 OSCC samples and 32 normal controls were enrolled in our study.
Subsequently, total of 232 autophagy related genes were extracted from transcriptome profile in R software(Version 3.6.1). Subsequently, differentially expressed analysis was performed in R software EdgeR package with the cut-off criteral |log2 (fold change [FC])|>2.0 and FDR(adjusted P-value) <0.01.

GO and KEGG analysis
Gene Ontology (GO) was perfomed to analyze these differentially expressed ATGs in DAVID (https://david.ncifcrf.gov/). In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was employed to annotate the functions. The significance level of p < 0.05 was taken as the cut-off standard.

Construction of cox regression model
The differentially expressed ATGs were transformed and normalized in a log2(x+1) manner [28].
Firstly, univariate cox regression was performed to screen prognosis relevant genes with the cut-off criteria p-value < 0.05). Then, stepwise regression was used to construct the cox risk model according to these prognostic genes to predict the overall survival of OSCC patients. Furthermore, the risk score for each OSCC patient was visualized in R software on basis of regression model.

Identification of cox regression model
The predictive power of the signature was evaluated using the Receiver operating characteristic(ROC). In addition, univariate and multivariate cox regression were used to analyze the clinical characteristics including age, gender, grade, stage, TNM and risk score. Distant metastasis was excluded owing to numerous missing datum.

Survival analysis based on cox model
According to cox model, OSCC patients in TCGA were divided into low and high risk group. Kaplan-Meier survival curves and stratification analysis along with a logrank p test were applied to validate its accuracy in R software survival package.

Collection of OSCC specimens
50 cases of OSCC tissue and matched normal mucosal tissues (MNTs) were collected from Nanfang Hospital, Southern Medical University. MNT at least 1.5 cm from the edge of the tumor was defined as a normal control. All 50 patients with OSCC were confirmed by pathology. All patients agreed to the study and signed the consent form.

RNA extraction and quantitative real-time PCR (qRT-PCR)
Total RNAs of tissues and OSCC cells were extracted with TRIzol reagent and reversed to cDNA.