Identification of autophagy-associated MicroRNAs and their prognostic significance in patients with laryngeal squamous cell carcinoma

Background: Autophagy is a biological process that cells engulf their cytoplasmic proteins or organelles to achieve the needs of metabolic and renewal. microRNAs can affect the development of cancer by regulating cell autophagy. We aim to identify the autophagy-associated MicroRNAs in laryngeal cancer and further explore their roles in the development of cancer. Results: we finally identified 6 key autophagy-associated microRNAs (AAMRs, hsa-miR-100-5p, hsa-miR-143-3p, hsa-miR-3607-5p, hsa-miR-454-3p, hsa-miR-455-5p, hsa-miR-99a-5p) were significantly correlated with prognosis of laryngeal squamous cell carcinoma (LSCC). Pathway enrichment analysis revealed that all key AAMRs were significantly enriched in the mTOR signaling pathway. The Risk index of each sample was calculated according to the result of multivariate COX analysis. Patients with lower risk index have better survival results. The area under the ROC curve for 1-year, 3-year, and 5-year survival rates were 0.7, 0.776 and 0.78 respectively. Conclusion: We identified 6 key AAMRs, which may act as new potential therapeutic targets. A new risk index model based on AAMRs can predict the prognosis of laryngeal squamous cancer.


Background
Laryngeal squamous cell carcinoma (LSCC) is common tumors of head and neck, accounting for 13.9% of all head and neck carcinomas. The incidence and mortality of laryngeal cancer were 0.8% and 0.6% respectively [1]. Although the recent report showed that the incidence of LSCC has dropped by 24% in the past 10 years, the five-year survival rate did not change significantly [1]. The chances of survival for patients with LSCC depend significantly on the initial stage of cancer. Compare with 80%-90% survival rate of stage I and II LSCC patients, the chances for survival fallen to 40% in LSCC patients with stage IV [2]. For the reason above, early detection and accurate diagnosis is a crucial part for patients to receive appropriate treatment.
Autophagy is a process that cells engulf their cytoplasmic proteins or organelles and coats them into vesicles. These vesicles can fuse with lysosomes to form autophagy lysosomes, which degrades the contents they encapsulate, to achieve the needs of metabolic and renewal. Autophagy can occur in various pathological processes, however, whether its role is positive or negative has not been fully elucidated, especially in the study of tumors. In the past decade, many studies have tried to inhibit the tumor growth process by targeting the autophagy pathway [3][4][5]. Besides, recent studies have demonstrated that autophagy inhibitor can inhibit the tumor growth and therapeutic resistance of laryngeal cancer cells [6]. Therefore, the exploration of potential autophagy targets is important for LSCC.
MicroRNA is a kind of non-coding RNA that can regulate the expression of mRNA at the posttranscriptional level [7], and an increased number of studies have revealed the dysregulation of microRNA remarkably affects various biological processes of cancer including autophagy [8,9]. Indeed, a great number of studies have revealed that microRNA can regulate many key proteins to affect various steps in the autophagy process [9]. For example, Recent research shows that miR-101-3p may inhibit autophagy through targeting EZH2 in LSCC [10]. Our study aims to establish a co-expression network of autophagyassociated microRNAs (AAMRs) by using bioinformatics methods, which may provide some theoretical basis for the targeted therapy of LSCC.
All in all, considering that AAMRs may play a significant role in LSCC and could be served as potential therapeutic targets. We constructed correlation networks to screen AAMRs and its prognostic value was also be systematically analyzed in this study.

identification of AAMRs
A total of 117 laryngeal cancer tissue samples and 12 normal samples' microRNA-seq data were downloaded from TCGA. 32 laryngeal cancer tissue samples and 5 normal samples' microRNA-seq data were downloaded from GSE124679. The number of differentially expressed microRNAs identified from TCGA and GSE124679 was 120 and 91 respectively (Fig. 2). The ARGs' expression data of corresponding LSCC tissue were also obtained from two databases. Two autophagy-microRNA co-expression networks were constructed respectively to screen candidate AAMRs (correlation coefficient > 0.3 and p value < 0.05).
After intersecting candidate AAMRs from two databases, 13 AAMRs were identified.

Identification of key AAMRs and their functional annotation
As described above, we have identified 13 AAMRs. Univariate Cox regression was then performed on these 13 microRNAs and the result showed that 6 of these AAMRs (hsa-miR-100-5p, hsa-miR-143-3p, hsa-miR-3607-5p, hsa-miR-454-3p, hsa-miR-455-5p, hsa-miR-99a-5p) were significantly correlated with the prognosis of LSCC (Fig. 3). These 6 AAMRs were defined as key AAMRs. Functional annotation was conducted based on key AAMRs by using the DIANA-mirPath tool. The result of pathway enrichment analysis showed that all key AAMRs were involved in the pathway of the 'mTOR' signaling pathway, which plays a significant role in various neoplastic processes including autophagy (Table 1).

calculation of risk index
To construct a reliable prognosis model, all six key AAMRs were included in the multivariate Cox analysis. The result showed that there are three AAMRs (hsa-miR-100-5p, hsa-miR-3607-5p and hsa-miR-99a-5p) could be independent factors significantly related to prognosis (Fig. 3B). We used Cytoscape software to construct a network of prognostic microRNAs with their co-expressed autophagy-related genes (Fig. 4). Based on the three prognostic AAMRs and the result of multivariate Cox regression. We calculated the level of the Risk index of each LSCC samples. The formula is as follows: Risk index = expression level of hsa-miR-100-5p * 0.672 + expression level of hsa-miR-3607-5p * -0.3248 + expression level of hsa-miR-99a-5p * -0.4086.

prognostic value of Risk index in LSCC patients
We used the median level of Risk index to divided LSCC samples into high Risk-index group and low Risk-index group. Figure 5A shows us the distribution of risk index in all included LSCC samples. Figure 5B shows the survival status of LSCC patients in two groups. We also use the heatmap to show the expression situation of prognostic AAMRs in different groups, it is clear that hsa-miR-100-5p expressed higher in high Risk-index group while hsa-miR-3607-5p and hsa-miR-99a-5p expressed higher in low Risk-index group to play a significant role in oncogenesis and tumor progression by acting either as cancer suppressor or proto-oncogenes. Autophagy is a really important biological phenomenon. In some models, cancer may stimulate the occurrence of autophagy to maintain the function of mitochondrial and energy balance [11]. Therefore, inhibiting the process of autophagy could be beneficial for cancer treatment. A great number of studies have shown that microRNAs can inhibit the development of various cancer by regulating the process of autophagy [12][13][14]. However, up to now, there has been no systematic method to identify Autophagy-associated microRNAs in LSCC. Therefore, it is necessary to find a way to screen AAMRs and explore their roles in the occurrence and development of LSCC.
In our study, we intersect the candidate AAMRs screened from two datasets (TCGA and Specifically, all 6 key AAMRs were enriched in the mTOR signaling pathway, which plays a vital role in the regulation of cell growth, survival and autophagy. mTOR protein could function as a sensor of energy and nutrient levels, thus further negatively regulated the occurrence of autophagy [21]. The research showed that inhibition of mTOR can significantly affect the proliferation and apoptosis in laryngeal cancer cell lines [22]. Our result indicates that AAMRs may affect the development of LSCC through autophagy and mTOR signaling pathways.
Our study dose has some limitations. Firstly, the data source for our analysis was based on TCGA and GEO database, thus making it impossible to obtain all patient information.

Conclusion
In conclusion, through constructing the autophagy-microRNA correlation network and using Cox regression analysis, 6 key autophagy-associated microRNAs, which have prognostic value for LSCC, were identified. These 6 key AAMRs may act as new potential therapeutic targets. Besides, we also constructed a new model based on the risk index to better predict the prognosis of LSCC patients. LSCC samples in high-risk group and lowrisk group may show different autophagy states.

screen of microRNAs and autophagy-related genes
The profiles of microRNAs and autophagy-related genes were downloaded from the GEO dataset (GSE124679) and the TCGA dataset. The list of autophagy-related genes (ARGs) (n = 223) was obtained from the Human Autophagy Database (http://autophagy.lu/clustering/index.html). Differential analyses were performed on microRNAs data between tumor tissue and normal tissue. The Pearson correlation network was constructed to explore the relationship between ARGs and differentially expressed microRNAs (DE-miRNAs). We defined the correlation coefficient > 0.3 and p value < 0.05 as the inclusion criteria of candidate autophagy-associated microRNAs (AAMRs). Finally, we intersected the candidate AAMRs from two datasets as AAMRs. A detailed flow chart for screening AAMRs was shown in Fig. 1. The differentially expressed microRNAs in two datasets. Figure 2A Volcano plot showed the DE-miRNAs in TCGA. Figure 2B Volcano plot showed the DE-miRNAs in GSE124679.

Figure 3
The result of Cox regression. Figure 3A, The result of univariate Cox regression for key AAMRs. Figure 3B, The result of multivariate Cox regression for key AAMRs.

Figure 4
Network of prognostic microRNAs with their co-expressed autophagy-related genes Figure 5 prognostic value of Risk index. Figure 5A Distribution of risk index. Figure 5B Survival status of LSCC patients. Figure 5C expression situation of prognostic AAMRs in different groups. Figure 5D Kaplan-Meier curve of LSCC patients based on the level of risk index. forest plot. Figure 6A Forest plot of univariate Cox regression for LSCC patients.