Screening and identification of autophagy‐related biomarkers for oral squamous cell carcinoma (OSCC) via integrated bioinformatics analysis

Abstract Increasing evidences have showed that autophagy played a significant role in oral squamous cell carcinoma (OSCC). Purpose of our study was to explore the prognostic value of autophagy‐related genes (ATGs) and screen autophagy‐related biomarkers for OSCC. RNA‐seq and clinical data were downloaded from The Cancer Genome Atlas (TCGA) database following extracting ATG expression profiles. Then, differentially expressed analysis was performed in R software and a risk score model according to ATGs was established. Moreover, comprehensive bioinformatics analyses were used to screen autophagy‐related biomarkers which were later verified in OSCC tissues and cell lines. A total of 232 ATGs were extracted, and 37 genes were differentially expressed in OSCC. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis demonstrated that these genes were mainly located in autophagosome membrane and associated with autophagy. Furthermore, the risk score on basis of ATGs was identified as potential independent prognostic biomarker. Moreover, ATG12 and BID were identified as potential autophagy‐related biomarkers of OSCC. This study successfully constructed a risk model, and the risk score could predict the prognosis of OSCC patients accurately. Moreover, ATG12 and BID were identified as two potential independent prognostic autophagy‐related biomarkers and might provide new OSCC therapeutic targets.

Autophagy is characterized by sequestration of bulk cytoplasm, long-lived proteins and damaged cellular organelles encapsulated as autophagosomes and delivered for lysosomal degradation to recycle the nutrients, and it is a cellular self-consumption process. Autophagy is regulated by the conversion of the cytoplasmic microtubuleassociated protein 1 light chain 3 (LC3-I) into the membrane form of protein 2 light chain 3(LC3-II) which promotes the autophagosomal degradation. 4 In summary, autophagy is a lysosomal-mediated catabolic complex process which involves the cytoplasmic organelles and proteins to maintain metabolism and homeostasis in cells. 5 Recently, numerous studies have showed that autophagy dysregulation played significant roles in a variety of human malignancies, including colorectal cancer, 6 renal cell carcinoma, 7 non-small cell lung cancer 8 and so on. In addition, majority of evidences have showed that autophagy might play significant roles in OSCC carcinogenesis.
For example, autophagy-mediated cell apoptosis to promote tumour progression via the AKT/mTOR pathway in OSCC. 9 Dysregulation of autophagy process was relevant to tumorigenesis and prognosis in OSCC. 10 Therefore, bioinformatics analysis of ATGs may reveal its prognostic value and provide potential therapeutic targets for OSCC treatment.
The purpose of this study was to analyse the differentially expressed ATGs in OSCC and establish a cox regression model to predict the overall survival of OSCC patients. Furthermore, survival analyses combined with stratification analyses were performed to identify the accuracy of cox formula. The present study may provide a novel insight into the potential mechanisms in OSCC initiation and progression and new therapeutic targets for OSCC treatment.
Owing to its half-baked clinical data, two samples were excluded.
Finally, 317 OSCC samples and 32 normal controls were enrolled in our study. Subsequently, a total of 232 ATGs were extracted from transcriptome profiles in R software (Version 3.6.1). Then, differentially expressed analysis was performed in R software EdgeR package with the cut-off criteria |log 2 (fold change [FC])|>1.0 and FDR (adjusted P-value)<.01.

| GO and KEGG analysis
Gene Ontology (GO) analysis was performed to analyse these differentially expressed ATGs in DAVID database (https://david.ncifc rf.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 < .05 was taken as the cut-off standard.

| Construction of cox regression model
The ATG profiles were transformed and normalized in a log 2 (x + 1) manner. 11 Firstly, all ATGs were enrolled in univariate Cox regression to screen prognostic genes with P-value < 0.05. Then, stepwise regression analysis was used to construct the Cox risk model according to these prognostic genes. Finally, 9 ATGs were enrolled in risk Cox regression, and the visualizations of Cox model were performed in R software according to the risk score for each patient. In addition, the predictive power of the signature was evaluated using the receiver operating characteristic (ROC). And univariate and multivariate cox regression analyses were used to analyse the clinical characteristics including age, gender, grade, stage, T, N stage in TNM system and risk score. M classification was eliminated because of lots of missing data.

| Survival analysis based on cox model
On basis of cox model, OSCC patients in TCGA were segmented into high-and low-risk level group. Survival analysis and stratification analysis including age, gender, grade, stage, T and N classification in TNM system along with a log rank p test were applied to validate its accuracy in R software survival package.

| Western blot
RIPA lysis buffer was used to extract OSCC cell lines' protein.

| Statistical analysis
SPSS23.0 software (IBM) was used for statistical analysis. T test and chi-square test were used in our study. P-value less than .05 was considered statistically significant.

| Identification of differentially expressed ATGs
317 OSCC patients and 32 normal controls RNA-seq and corresponding clinical data were downloaded from TCGA database. 232 ATGs expression level were extracted from the transcriptome data subsequently, and differentially expressed analysis was performed in R software EdgeR package ( Figure 1A). With the cut-off criteria |log 2 FoldChange|>1 and FDR < 0.01 ( Figure 1B,C), 11 up-regulated and 26 down-regulated ATGs were sorted out. GO enrichment analysis indicated that these genes were mainly located in autophagosome membrane and associated with autophagy ( Figure S1A,C).
KEGG pathway analysis showed that most of significant ATGs were enriched in apoptosis, platinum drug resistance, ErbB signalling pathway and TNF signalling pathway (Figure S1B,D). were divided into high-risk and low-risk group according to Cox formula median. Survival analysis indicated that the overall survival rate of high-risk group was significantly lower than that of the lowrisk group ( Figure 2C). Moreover, the expression levels of protective ATGs in the low-risk group were higher than that of high-risk group.

| Establishment of cox regression model
On the contrary, the expression levels of potential oncogenes were higher in high-risk group ( Figure 2D). In addition, the risk scores combined with survival data were visualized in R software ( Figure 2E,F).

| Identification of Cox regression model
Firstly, the ROC on basis of cox model was plotted, and its area under the curve (AUC) was 0.76 which was markedly higher than other clinical characteristics ( Figure 3A). Furthermore, risk score in early stage was significantly lower than that in advanced stage ( Figure 3B) indicating that the risk score on basis of ATGs might realize early diagnosis in OSCC. Moreover, univariate and multivariate Cox regression analysis indicating that the risk score might be regarded as an independent prognostic factor ( Figure 3C,D).
Moreover, ATG12 and BID were identified as 2 independent autophagy-related biomarkers according to univariate, multivariate Cox regression analysis (Figure 2A,B) and survival analysis ( Figure 3E,F). Unfortunately, ATG12 was not the differentially expressed ATGs. However, integrated analysis indicated that ATG12 might play a carcinogenic role in OSCC. The reasons why ATG12 was not differentially expressed ATGs in OSCC might be that the cut-off criteria was too high. Therefore, statistical analysis between high and low ATG12 expression group was performed, and the result showed that ATG12 was significantly up-regulated in tumour samples ( Figure S2A).

| Survival analysis
The prognostic value of the risk score for different clinicopathological parameters including age, gender, T and N in TNM system, grade and stage was further investigated. M stage in TNM system was excluded because of numerous data missing. Survival analysis combined with stratification analysis including age, gender, T, N, grade and stage demonstrated that low-risk group had significantly higher overall survival rate than high-risk group (Figure 4).

| Identification of potential independent prognostic biomarkers
Comprehensive bioinformatics analysis indicated that ATG12 and BID might be associated with the overall survival and played a carcinogenic role in OSCC. Univariate and multivariate Cox regression analysis showed that ATG12 and BID might be selected as potential independent prognostic biomarkers in our study. Therefore, the expression levels of ATG12 and BID were validated in OSCC cell lines and tissues by qRT-PCR, Western blot and immunohistochemistry.
Our results revealed that expression of ATG12 and BID at mRNA level were up-regulation in OSCC cell lines ( Figure 5A,B) and in 50 OSCC patients than MNTs ( Figure 5C,D), which were similar with the results in TCGA database. Moreover, ATG12 and BID were also up-regulated in OSCC cell lines and tissues by Western blot assay and IHC staining, respectively ( Figure 5E, Figure 5F,G, Figure S2B,C).
However, the mRNA levels of ATG12 and BID had no significant correlation with the clinical parameters ( Table 1, Table 2).

| D ISCUSS I ON
Despite that advancements in surgery, radiation and chemotherapy for the treatment of OSCC, the exact pathogenesis of OSCC is not clear and the 5-year overall survival is still not satisfied. 13 Hence, it is crucial to explore the effective biomarkers and therapeutic targets to improve the overall survival of OSCC. Recently, increasing evidences revealed that autophagy might play a significant role in OSCC initiation and progression. 10,14,15 Therefore, our study aimed to analyse the ATGs by bioinformatics methods and then filter the potential therapeutic targets for OSCC.
A total of 232 ATGs were enrolled in this study. 37 genes were differentially expressed in OSCC which were associated with autophagy and apoptosis by GO enrichment analysis, indicating that these differentially expressed ATGs were relevant to cancer progression. 16,17 The KEGG pathway analysis revealed that these ATGs were mainly involved in 3 pathways ( Figure S1B), of which platinum drug resistance pathway might play an important role in OSCC treatment strategies 18 and ErbB signalling pathway was relevant to to head and neck squamous cell carcinoma (HNSCC) treatment, 19 indicating that these ATGs might provide potential therapeutic targets for OSCC. In fact, some In the present study, we established a risk model according to 27 prognostic ATGs to predict the overall survival of OSCC. Ultimately, 9 ATGs were selected in Cox regression model and the risk score based on ATG expression level could predict the overall survival accurately. Meanwhile, univariate, multivariate Cox regression and survival analysis combined with stratification analysis indicated that the risk score could be regarded as an independent prognostic factor and distinguish early stage from advanced OSCC, which might be helpful for OSCC early diagnosis and prediction. Furthermore, ATG12 was identified as a potential prognostic biomarker in OSCC which is the human homolog of a yeast protein involved in autophagy. 21  ATGs should be further investigated to explore the roles in OSCC in an autophagy manner.
In summary, we analysed ATGs comprehensively and identified that these differentially expressed ATGs were relevant to OSCC initiation and progression significantly. Furthermore, a 9 ATG gene signature was successfully constructed which was positively associated with overall survival of OSCC. The risk formula might provide potential therapeutic targets for OSCC early diagnosis and treatment.
However, there are some limitations in our study. First, the sample sizes in TCGA database and OSCC specimens are markedly inadequate. In addition, the mechanism of ATG12 and BID tumorigenic role in OSCC should be further investigated.

| CON CLUS ION
This study successfully constructed a risk model to predict the prognosis of patients with OSCC through comprehensively analysing ATGs, and the risk score might be an independent prognostic biomarker in OSCC. Moreover, ATG12 and BID were also identified as two potential independent prognostic biomarkers for OSCC diagnosis and treatment.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

E TH I C A L A PPROVA L
The Nanfang Hospital ethics committee (AF/SC-09/03.2) approved this study.

CO N S E NT FO R PU B LI C ATI O N
All co-authors were consent for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data for this study are available from corresponding authors if required.