Development of an autophagy-related gene prognostic signature in lung adenocarcinoma and lung squamous cell carcinoma

Introduction

Lung cancer is a fatal malignancy worldwide and is one of the leading causes of death caused by malignant tumors. In 2019, the mortality of lung cancer still ranks first among those of all kinds of cancers in the United States (Siegel, Miller & Jemal, 2019). For worldwide, it is also the leading cause of cancer death among men and the second leading cause of cancer death among women (Torre, Siegel & Jemal, 2016). More than half (57%) of lung cancer patients are diagnosed at the time of the distant stage (Torre, Siegel & Jemal, 2016). Even patients who underwent surgical resection, chemotherapy, radiotherapy and targeted therapy did not have significantly improved survival times. The five-year survival varies from 4–17%, leads to a need to explore new therapeutic targets and treatments (Gray et al., 2019; Hirsch et al., 2017). According to the histological classification, lung cancer is divided into small-cell lung cancer (SCLC) and non-small-cell lung cancer (NSCLC), the latter of which accounts for approximately 85% of all cases (Neal, Hamilton & Rogers, 2014). Squamous cell carcinoma and adenocarcinoma account for approximately 90% of the total NSCLC cases, which make them the most common types of lung cancer (Neal, Hamilton & Rogers, 2014). The poor therapeutic effect of NSCLC is mainly due to the lack of effective indicators for detecting the development of tumors at the early stage. Therapeutic progress of NSCLC is approached by the advances in the molecular field and the development of new drugs that target molecular abnormalities. But the existing treatment targets are prone to inducing resistance. New treatment markers and targets are needed to achieve better prognosis. The identification of aberrant genes has been a hot topic, in which the research on autophagy has a great prospect. Autophagy is the phagocytotic process of engulfing cytoplasmic proteins, complexes or organelles. The autophagosome, a cytoplasmic double-membrane structure, can be transported into lysosome and fusion with lysosome to generate the autolysosome (Galluzzi et al., 2015; Levine & Kroemer, 2008; Maiuri et al., 2007b). The degradation products can be transported back and recycled for general cell metabolism. Generally speaking, autophagy has a dual function in the apoptosis, which means it has positive and negative effects. It is suppressed by carcinogenic proteins to prevent excessive protein degradation in stressed tumor cells. Meanwhile, persistent autophagy activation leads to apoptosis (White, 2015; Ye et al., 2012). According to the present understanding, autophagy is involved in the innate and adaptive immune responses and can be induced by immune receptors such as Toll-like receptors and NLRs (nucleotide oligomerization domain-like receptors) (Cadwell, 2016). It takes part in the process of antigen presentation and the development of lymphocytes (Zhong, Sanchez-Lopez & Karin, 2016), which makes autophagy a possible target for improving immunotherapy in NSCLC. The relationship between autophagy and NSCLC has not been fully revealed, and some studies have suggested a role of autophagy in the targeted drug resistance. For example, in patients with NSCLC, EGFR tyrosine kinase inhibitors (EGFR-TKIs) and anaplastic lymphoma kinase inhibitors (ALK) can be used as effective treatments. EGFR-TKIs can induce autophagy, and high levels of autophagy after treatment with EGFR-TKIs may also lead to autophagic death of the cells (Lee et al., 2015). Hence, the combination of EGFR-TKIs with autophagy inducers may be beneficial. A similar synergistic effect can be observed with ALK inhibitor resistance (Ji et al., 2014). In this study, we revealed an autophagy-related risk signature involving six genes. This signature can be used as an independent prognostic marker for LUAD and LUSC patients. Our study indicates that autophagy may be a promising target for the treatment of NSCLC.

Material and Methods

Autophagy-related gene datasets and patient samples

The gene expression datasets and clinical information of LUAD and LUSC patients were downloaded from TCGA database on September 9, 2019. The supplementary clinical information of corresponding patients was obtained from cBioPortal (http://www.cbioportal.org). An independent microarray NSCLC cohort was extracted from the GEO database (accession number: GSE3141). Overall, the expression data from 1,102 samples (103 normal samples and 999 tumor samples) were obtained with the TCGA dataset. A total of 111 samples from the GEO dataset GSE3141 were used in the verification group. A total of 232 genes from the HADb (Human Autophagy Database) were identified as autophagy-related genes.

Results

Functional enrichment analysis of the differentially expressed genes

Functional enrichment analysis of the 41 differentially expressed genes offered a biological understanding of these genes. The GO term functional enrichment and the KEGG pathway enrichment analyses of these genes are summarized in Figs. 3 and 4.

The top enriched GO terms for biological processes were autophagy, processes utilizing autophagic mechanisms, and the intrinsic apoptotic signaling pathway. Cellular components included the autophagosome membrane, the autophagosome, and integral components of the mitochondrial outer membrane. On the basis of molecular function, genes were mostly enriched in terms of protein phosphatase binding, phosphatase binding, and p53 binding. In the KEGG pathway enrichment analysis, these genes were shown to be notably associated with pathways in the autophagy (animal), apoptosis, and bladder cancer pathways. Most of the Z-scores of enriched pathways were more than zero, indicating that most of the pathways were more likely to be enhanced.

Identification of an autophagy-related risk signature for the prognosis of LUAD and LUSC

By revealing the distinct expression patterns found in LUAD and LUSC patients, we considered that identifying an autophagy-related risk signature might be useful for predicting prognosis. A univariate Cox regression analysis was performed to establish a candidate pool of autophagy-related genes (Fig. 5A). Ultimately, five genes (HDAC1, ATG4D, TP73, EIF4EBP1 and TP63) were identified as protective factors (HR < 1), while another five genes (BNIP3, DAPK1, ATIC, ERO1A and FADD) were identified as risk factors (HR > 1). Subsequently, a multivariate Cox analysis was conducted. As a result, EIF4EBP1, TP63, BNIP3, ATIC, ERO1A and FADD were identified as independent prognostic indicators for DFS and selected for development of the prognostic signature.

According to the multivariate Cox proportional hazards regression model, we obtained the expression coefficient of each independent risk gene. Our prognostic model for predicting prognosis based on the six genes was formed using the following formula: prognosis index (PI) = (−0.170 * expression level of EIF4EBP1) + (−0.057 * expression level of TP63) + (0.117 * expression level of BNIP3) + (0.170 * expression level of ATIC) + (0.214 * expression level of ERO1A) + (0.268 * expression level of FADD). We then calculated the risk score of each patient and used the median risk value as a cutoff point for classifying patients into high-risk group (n = 732) or low-risk group (n = 370), as shown in Figs. 5D and 5E. The heatmap of these six signature-related genes and the Kaplan–Meier curve depending on risk score are also displayed (Figs. 5B and 5C). A significant difference in survival between the high-risk group and the low-risk group was observed. Patients in the high-risk group had a shorter OS than patients in the low-risk group (five-year survival rate = 36.7% vs. 44.9%, p = 0.0017). Similar results could also be seen with the DFS (median time = 0.344 years vs. 0.512 years, p <0.001). The results of Kaplan–Meier analysis also showed a prognostic ability of each single gene. The downregulation of EIF4EBP1 was strongly correlated with inferior DFS in LUAD and LUSC patients (P < 0.05; Fig. 6A). Similarly, low expression of TP63 led to inferior DFS (P < 0.05; Fig. 6B). In contrast, the upregulation of ATIC, ERO1A and FADD indicated a decreased DFS (P < 0.05; Figs. 6D, 6E and 6F). However, we did not observe a significant difference in DFS with regard to BNIP3 expression (P < 0.05; Fig. 5C). Considering that each gene had a different value in the prognostic model, a statistical difference may not occur in each survival analysis of a single gene. ROC curves of OS and DFS were used to reveal the predictive performance of the six-gene risk signature (Fig. 7). The AUC values of the signature for OS and DFS were 0.656 and 0.671, which were obviously higher than those associated with age (AUC = 0.547 and 0.478, respectively), sex (AUC = 0.551 and 0.502, respectively), tumor stage (AUC = 0.634 and 0.641, respectively), tumor T stage (AUC = 0.629 and 0.648, respectively), tumor N stage (AUC = 0.578 and 0.633, respectively) and tumor M stage (AUC = 0.501 and 0.489, respectively). These results indicated that the risk signature had a better ability to predict the survival of LUAD and LUSC patients than did clinical factors.

Associations between the autophagy-related risk signature and clinicopathologic features in LUAD and LUSC patients

An analysis was performed to explore the associations between clinical parameters and the risk signature (Fig. 8). The results showed that the signature was significantly associated with tumor stage (p = 0.006), M stage (p = 0.004), and survival outcome (p < 0.001). Additionally, Student’s t-test analysis also indicated that these signature-related genes were differentially expressed across various clinicopathological parameters. As shown in Fig. 9, differential ATIC expression was found across different tumor stages, M stages and survival outcomes. Differential expression of BNIP3 was observed across different tumor stages and M stages. EIF4EBP1 showed different expression across different sexes and survival outcomes. ERO1A, showed differential expression across sexes, tumor stages and T stages. A difference in the expression of FADD was observed across ages and sexes. The differential expression of TP63 was related to survival outcome, sex and M stage.

Discussion

Autophagy is a process carried out by cells to address nutritional deficiency and other cellular stresses. It is involved and regulated by a series of proteins and is closely correlated with a variety of cell processes and reactions. In recent years, a large number of studies have shown that autophagy is a “double-edged sword” in the occurrence and treatment of tumors. On the one hand, as a survival mechanism of tumor cells facing pressure, it plays a role in protecting cells. On the other hand, autophagy can result in killing tumor cells under certain conditions. Therefore, autophagy is considered to be a possible regulatory point for improving the therapeutic effects of tumor-targeted drugs and reducing drug resistance.

Autophagy has promise for improving the survival of NSCLC, but most studies usually focus on the role of a particular gene related to autophagy. The large-scale databases, such as TCGA and GEO, provide us with effective measures to explore gene signatures, thus providing a better understanding of the relationship between autophagy and tumors. In this study, based on the existing gene data of patients with NSCLC, we screened autophagy-related genes and identified six key prognostic genes, all of which may be possible molecular biomarkers of prognosis and potential therapeutic targets. We verified the autophagy-related genes in multiple datasets, which proved that the signature had very good prognostic ability across data from multiple centers.

The GO and KEGG analyses were also conducted to show the molecular and biological pathways enriched. The results suggested that the top enriched GO terms in terms of biological processes and cellular components were highly correlated with autophagy. On the basis of molecular function, p53 binding is closely related with the autophagy-related gene TP63 which will be discussed in detail later.

In addition, in the KEGG analysis, the most significant pathway was also enriched in autophagy processes. Because of this result, we speculated that specific autophagy pattern

The protein encoded by EIF4EBP1 binds to eukaryotic translation initiation factor 4e (EIF4E) and suppresses the EIF4E complex, thus affecting the mTOR (mammalian target of rapamycin) signaling pathway, which has been shown to promote tumorigenesis (Karlsson et al., 2013; Karlsson et al., 2011). EIF4EBP1 may promote or inhibit the development of tumors as a bi-functional factor (Armengol et al., 2007; Cai, Ye & She, 2014; Martineau et al., 2013). In general, phosphorylated EIF4EBP1 is considered to be an indicator of tumor activity, indicating a worse prognosis. Unphosphorylated EIF4EBP1 (Cai, Ye & She, 2014), on the other hand, is thought to inhibit tumor activity (Martineau et al., 2013). The TP63 gene is located on chromosome 3q27/29 and is closely related to human head, neck, esophagus, lung and skin squamous cancers (Massion et al., 2003; Senoo et al., 2001).

It encodes many subtypes of p63 transcription factors, which are members of the p53 protein family, an important hub in the transcriptional and signaling networks of the epithelial cells. Therefore, the dysregulation of TP63 is closely related to the occurrence of squamous cell carcinoma (Romano, Solomon & Sinha, 2012). The genome-wide analysis suggests that the genetic variant of TP63 may affect susceptibility to LUAD (Hu et al., 2014; Hu et al., 2011; Miki et al., 2010). However, there is still a lack of research on the specific mechanism.

BNIP3 was shown to be one of the most important players involved in autophagy. It encodes proteins belonging to the Bcl-2 family, which can regulate programmed cell death under some circumstances and may impart some pro-death activity (Vande Velde et al., 2000). It is linked to autophagy through three potential mechanisms. First, it can induce mitochondrial dysfunction to activate autophagy (Scherz-Shouval & Elazar, 2011). At the same time, through competitive binding with the BCL2 complex, BNIP3 can liberate Beclin-1 to induce autophagy (Bellot et al., 2009; Maiuri et al., 2007a). In addition, by inhibiting Rheb, an upstream activator of mammalian target of rapamycin (mTOR), BNIP3 may activate autophagy by repressing mTOR (Li et al., 2007). A study on early, operable NSCLC showed that the high expression of BNIP3 was an independent predictor of poor OS (Giatromanolaki et al., 2004). In vitro experiments also confirmed that BNIP3 participated in lung cancer cell migration by interacting with aryl hydrocarbon receptor (AhR) (Tsai et al., 2017). ATIC is a protein enzyme that acts on the last two steps of the de novo purine biosynthetic pathway (Martinez-Outschoorn et al., 2017; Yamaoka et al., 1997). Little has been revealed in regard to its function in lung cancer. Recent studies have shown that inhibiting the activity of ATIC metabolites may be important for the anti-tumorigenic effects of the drug pemetrexed, which is used against NSCLC (Racanelli et al., 2009).

ERO1A is a main regulatory factor of protein disulfide isomerase (PDI), which is one of the most abundant proteins in the endoplasmic reticulum. Endoplasmic reticulum stress is reported to be associated with tumorigenesis in a variety of cancers, including NSCLC (Kim et al., 2012; Cancer Genome Atlas Research, 2008). As an important component of the endoplasmic reticulum, PDI is believed to be a marker of poor prognosis in patients with tumors, such as glioblastoma, breast cancer and hepatocellular carcinoma (Shai et al., 2003; Thongwatchara et al., 2011; Xia et al., 2017). It is not surprising that ERO1A, as the main regulator of PDI, is also associated with poor prognosis in NSCLC (Hsu et al., 2016).

FADD was originally described as an adapter molecule for apoptosis and is the key to transmitting death signals from cell surface receptors (Mouasni & Tourneur, 2018). It is closely related to autophagic cell death and tumor development. Similarly, high expression of FADD was observed in NSCLC, and it was considered to be associated with the increased invasive behavior of the tumor and a marker for predicting prognosis (Chen et al., 2005; Luo et al., 2018).

We summarized the relevant information of the six autophagy-related genes (Table 1). It can be seen that half of them have a two-side role in the development of cancer (in addition to BNIP3, FADD and ERO1A, a consistent tumor-promoting effect is present). These results are in agreement with the current consensus that autophagy plays a bi-functional role in tumors. As we mentioned previously, most of the current studies only target an individual gene. However, in view of the complex effect of autophagy, it may not be appropriate and beneficial to simply inhibit or induce some of the autophagy-related genes based on these findings. Our study suggests that autophagy-related genes may affect the cancer development through special pattern, and by which autophagy may show a consistent effect. Revealing these specific patterns can help the clinicians identify the high-risk types and use them as new therapeutic targets. Our signature based on autophagy-related genes also confirms this conjecture.

In summary, the molecular mechanisms play an important role in the relationship between autophagy and NSCLC. Our results are expected to be applied to clinical practice, which means it may suggest potential targeted autophagy therapies for NSCLC. Further investigations will provide more information of internal mechanisms. Our study first reveals that autophagy-related pattern may affect the prognosis of patients with LUAD and LUSC. And a signature is presented to help distinguish the high-risk patients. However, a limitation of this study is its retrospective nature. More prospective studies should be conducted to validate the prognostic function of autophagy-related signatures. We also encourage multi-center data to confirm our findings. More laboratory data based on the thought of autophagy pattern can further develop our study and provide the internal mechanisms of autophagy-related network.