Prognosis-related autophagy genes in female lung adenocarcinoma

Abstract To screen the prognosis-related autophagy genes of female lung adenocarcinoma by the transcriptome data and clinical data from The Cancer Genome Atlas (TCGA) database. In this study, screen meaningful female lung adenocarcinoma differential genes in TCGA, use univariate Cox proportional regression model to select genes related to prognosis, and establish the best risk model. In this study, Gene Ontology and Kyoto Encyclopedia of Genes and Genomes were applied for carrying out bioinformatics analysis of gene function. The gene expression and clinical data of 264 female lung adenocarcinoma patient samples were downloaded from TCGA. Twelve down-regulated genes: NRG3, DLC1, NLRC4, DAPK2, HSPB8, PPP1R15A, FOS, NRG1, PRKCQ, GRID1, MAP1LC3C, GABARAPL1. Up-regulated 15 genes: PARP1, BNIP3, P4HB, ATIC, IKBKE, ITGB4, VMP1, PTK6, EIF4EBP1, GAPDH, ATG9B, ERO1A, TMEM74, CDKN2A, BIRC5. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analysis showed that these genes were significantly associated with autophagy and mitochondria (animals). Multifactor Cox analysis of autophagy-related genes showed that ITGA6, ERO1A, FKBP1A, BAK1, CCR2, FADD, EDEM1, ATG10, ATG4A, DLC1, VAMP7, ST13 were identified as independent prognostic indicators. According to the multivariate Cox proportional hazard regression model, there was a significant difference in the survival rate observed between the high-risk group (n = 124) and the low-risk group (n = 126) during the 10-year follow-up (P < .05). Univariate Cox analysis showed that tumor stage, T, M, and N stages, and risk score were all related to the survival rate of female lung adenocarcinoma patients. Multivariate Cox analysis found that autophagy-related risk scores were independent predictors, with an area under curve (AUC) value of 0.842. At last, there is autophagy genes differentially expressed among various clinicopathological parameters: ATG4A, BAK1, CCR2, DLC1, ERO1A, FKBP1A, ITGA6. The risk score can be used as an independent prognostic indicator for female patients with lung adenocarcinoma. The autophagy genes ITGA6, ERO1A, FKBP1A, BAK1, CCR2, FADD, EDEM1, ATG10, ATG4A, DLC1, VAMP7, ST13 were identified as prognostic genes in female lung adenocarcinoma, which may be the targets of treatment in the future.


Introduction
Lung cancer is the most common malignancies with high morbidity worldwide, which represents the second leading cause of cancer-related deaths. [1] Data from the National Cancer Registry [2] show a significant increase in lung cancer incidence in women and a gradual decrease in men, while men have an overall survival (OS) advantage in the patients with non-small cell lung cancer (NSCLC). [3] Although current treatment approaches, including maximum safe resection, chemotherapy, radiotherapy, and targeted therapy, have been adopted, [4] the survival time has not improved significantly, and the 5-year survival rate is 4% to 17%. [5] Despite the progress of experimental technologies and therapeutic regimens in the field, lung cancer remains incurable. It is thus necessary to explore the novel therapeutic targets or treatment methods.
Autophagy is a genetically conserved cellular process and widely exists in eukaryotic cells, which is one of the necessary ways of organelle renewal and metabolism dynamic balance, [6] especially for maintaining homeostasis in vivo and tumor growth of lung cancer, thus it prevents the p53 activation, growth arrest, apoptosis, senescence, and immune response activation. Autophagy in NSCLC can maintain mitochondrial mass and regulate its abundance, while the genetic deficiency of essential autophagy genes in tumors will reduce tumor growth. [7] Even DNA methylation of autophagy genes plays an important role in cancer development, and the epigenetic-mediated aberrant gene silencing has been considered to be closely related to the pathogenesis of cancers, including lung cancer. [8] Generally, autophagy is thought to prevent carcinogenesis by eliminating oncogenic protein substrates, misfolded proteins, and damaged organelles. [9] Therefore, more and more studies have shown that autophagy has an important role in NSCLC. [10] 2. Materials and methods

Methods
Gene expression data, clinical characteristics, and survival information of 264 samples from female patients with lung adenocarcinoma were downloaded from The Cancer Genome Atlas database on February 27, 2021. A total of 260 autophagyrelated genes were extracted from Human Autophagy Database (http://www.autophagy.lu/index.html).

Statistical analysis
A series of gene function enrichment analysis was performed after screening the differentially expressed autophagy-related genes, to determine the main biological properties. Univariate Cox proportional hazard regression analysis was used to assess the association between OS and the selected autophagy genes. Next, multivariate Cox proportional hazard regression analysis was performed on the identified autophagy genes. Independent prognostic factors were determined by multivariate Cox proportional hazards regression analysis, and regression coefficient and risk ratio (HR) were calculated by Cox regression model for establishing an autophagy related model, which can be used as the prognosis model of female lung adenocarcinoma. The risk score of patients was divided into high-risk group and lowrisk group according to the median score. Survival curve, risk curve, receiver operating characteristic (ROC) curve of risk score, and other clinical characteristics, and beeswarm of correlation between autophagy gene, the risk score, and clinical characteristics were drawn.
All the statistical analyses including the different expression of autophagy-related genes in univariate and multivariate Cox regression models, ROC curve, and K-M survival were performed using R version 3.6.6 (The R Foundation for Statistical Computing). GO plot software (COLORSPACE, STRINGI, GGPLOT2, DOSE, CLUSTERPROFILER, ENRICH-PLOT) package was used to Gene Ontology (GO) enrichment analysis, and the same tool is also used for the enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. P < .05 was used as the significant threshold.

Identification analysis of ATGs
The expression values of 260 ATGs were integrated from the female patients with lung adenocarcinoma, and 12 downregulated genes and 15 upregulated genes were identified. Boxplots (Fig. 1A), heat maps (Fig. 1B), and volcano maps (Fig. 1C) were used to reveal the expression patterns of these differentially expressed genes in tumor and non-tumor tissues.

Functional enrichment analysis
A total of 27 differentially expressed genes were performed functional analysis to provide biological function for these genes. In order to inquire the potential signal pathways related 27 ATGs in female lung adenocarcinoma, we analyzed them with GO and KEGG. GO analysis ( Fig. 2A-C) showed that these ATGs were enriched in some basic biological processes, including the process of autophagy and apoptotic signaling pathways. The cellular components affected those of autophagy and mitochondria. Based on molecular functions, the genes were most enriched in cellular response to oxidative stress, protein phosphatase binding, and phosphatase binding. KEGG analysis (Fig. 3A,B) showed that the 27 ATGs were significantly associated with autophagy and mitochondria (animals), and the most scores of enriched pathways were greater than 0, indicating that these pathways were more likely to be enhanced.

Construction of prognostic markers in female lung adenocarcinoma, diagnostic accuracy study (LUAD)
Twenty-eight ATGs were analyzed by univariate Cox regression. Eight genes were identified as protective factors (HR < 1), while another 20 genes were identified as risk factors (HR > 1). The 28 genes were further analyzed by multivariate Cox regression analysis, and finally 12 genes (ITGA6, ERO1A, FKBP1A, BAK1, CCR2, FADD, EDEM1, ATG10, ATG4A, DLC1, VAMP7, ST13) were identified as independent prognostic indicators, which might be helpful in predicting prognosis. The 12 genes were further analyzed by multivariate Cox regression analysis, the expression coefficients of each independent risk gene were obtained and shown in Table 1, and finally the model for predicting prognosis based was developed using the following formulas: prognostic index (PI) = (-0.162 * ITGA6 expression level) + (0.466 * ERO1A expression level) + (0.872 * FKBP1A expression level) + (0.373 * BAK1 expression   The risk score of each patient was calculated on the basis of the relevant mRNA expression level and risk coefficient of each ATG. The risk score is used to predict the prognosis of female LUAD, and the median risk score is the critical value to divide patients into high-risk and low-risk groups (Fig. 5B). Heatmap was drawn to show gene expression profiles in high-risk and low-risk female LUAD groups (Fig. 5C). The genes with HR > 1 were considered to be dangerous genes, while the genes with HR < 1 were the protective genes (Fig. 4B). As shown in Figure 4, patients in the high-risk group have more possibilities to express the risk genes. In contrast, patients in the low-risk group have a disposition to express the protective genes (Fig. 4B), and the Kaplan-Meier cumulative curve showed that the survival time of patients with low-risk score was significantly longer than that of patients with high-risk score at 10-year follow-up (P < .05) ( Fig. 4A; Fig. 5A).

Autophagy acted as an independent prognostic factor
We assessed the prognostic value of autophagy-related risk score by univariate Cox regression analysis and multivariate Cox regression analysis, the autophagy-related risk score in univariate analysis was significantly correlated with OS (HR = 1.235, 95% confidence interval (CI) = 1.173-1.299, P < .001) in T, M, N stage (Fig. 6A). Multivariate analysis further showed that autophagy-related risk score was an independent prognostic indicator (HR = 1.208, 95% CI = 1.143-1.277) (Fig. 6B). The results confirmed that the autophagy-related risk score could be used as an independent prognostic factor in clinic practice.

Multiindex ROC curve of risk score and other indicators
The ROC curve of OS has applied for revealing the predictive performance of the 12 autophagy-related gene risk scores (Fig. 6C). The area under curve (AUC) value of risk score (AUC = 0.842) was significantly larger than that of other indicators including age (AUC = 0.510), tumor stage (AUC = 0.755), T stage (AUC = 0.682), N stage (AUC = 0.662), and M stage (AUC = 0.578), suggesting that the risk would better predict survival in female patients with LUAD than other individual indicators. However, due to the lack of clinical data such as treatment regimen, personal history, and tumor grade, we were unable to perform ROC analysis for other clinical factors.

Association between autophagy-related risk characteristics and clinicopathological characteristics in female patients with LUAD
These related genes were differentially expressed among various clinicopathological parameters. As shown in the Figure 7, different ATG4A expression was found in different tumor stages (Fig. 7C), and different BAK1 expression was found in different M stages (Fig. 7E). Different CCR2 expression was found in different tumor stages (Fig. 7G), different DLC1 expression was found in different tumor stages and N stages (Fig. 7B, I), different ERO1A expression was found in different tumor stages and M stages (Fig. 7D, F), different FKBP1A expression was found in different N stages (Fig. 7H), and different ITGA6 expression was found in different age stages (Fig. 7J). Different risk score expression was found in different tumor stages (Fig. 7A).

Discussion
Recent years, many studies have shown that autophagy is closely related to the prognosis of LUAD, Considering the importance of autophagy in LUAD, we can reasonably speculate that autophagy related genes have broad prospects in the prognosis evaluation of female LUAD, but autophagy plays a complex role in the lung, which can produce protective and injurious effects on the development of lung diseases. Therefore, the design of autophagy as an effective therapeutic intervention strategy needs careful consideration and further study. [11] The role of autophagy is not constant, so in this study, based on the existing genetic data of female lung adenocarcinoma, after screening autophagy-related genes and identifying 12 key  prognostic genes that may be prognostic markers or potential therapeutic targets for female LUAD, the mechanisms of these 12 genes in lung cancer are summarized as follows.
ITGA6, a member of the integrin family, plays a key role in the interaction between many cell types and participates in multiple biological processes, including cell proliferation and invasion. [12] ITGA6 has been shown to play a carcinogenic role in NSCLC. [13] Co-expression analysis of ERO1A gene in cancer patients showed that ERO1A was significantly upregulated in lung adenocarcinoma, promoting cell migration, and invasion. [14] High expression of ERO1A is associated with poor OS in lung adenocarcinoma patients. [15] BAK1 is a key effector of apoptosis, [16] especially for mitochondrial apoptosis. [17] Several studies have confirmed that it can induce proliferation inhibition and apoptosis of NSCLC cells. [18] CCR2 expression is associated with tumor stage and metastasis. [19] Experiments have confirmed that endothelial CCR2 expression is necessary for tumor cell extravasation and lung metastasis. Endothelial cells lacking CCR2 can greatly reduce lung metastasis, but the growth of primary tumors is not affected. [20] ATG10, located in chromosome 5q14, is an autophagic E2-like enzyme and participates in autophagosome formation. [21] ATG10 expression is up-regulated in lung cancer tumors, and its high expression is a risk factor for death in lung cancer patients. [22] ATG4A has a proteolytic effect and plays an effective role in autophagy. [23] Variation of intron SNP rs807185 in ATG4A was significantly associated with reduced risk of lung cancer and may be a protective factor for lung cancer. [24] DLC1 is located on chromosome 8p21 3-22, the region shows loss of heterozygosity in many human cancers, including prostate cancer, colon cancer, breast cancer, ovarian cancer, liver cancer, lung cancer and bladder cancer. [25] Loss of DLC1 leads to dysregulation of the Rho pathway, which may contribute to cancer progression. [26] Among smokers with lung adenocarcinoma, DLC1 downregulation and TP53 mutations are as usual. [27] VAMP7 is involved in the transport of autophagosomes to the periphery of cells. [28] Meanwhile, VAMP7 is a key regulator of MT1-MMP-mediated matrix degradation and invasion. MT1-MMP-mediated pericellular proteolytic activity is necessary for cancer cells to break through the basement membrane and migrate mesenchymally in the fibrous collagen matrix. [29] FADD has been observed to be highly expressed in NSCLC and is thought to be associated with increased aggressive behavior and predictive markers of prognosis in NSCLC. [30,31] It can be seen that these genes play 2 roles in the development of cancer, especially in lung cancer, which are consistent with the consensus that autophagy plays a bidirectional role in tumors. However, in this study we found that autophagy related genes which showed consistent effects might promote the development of lung cancer in some way. This investigation suggests autophagy may be a promising target for female LUAD treatment in the future.
Among these genes, fkbp1a, edem1, ST13, and VAMP7 are more attractive, which have not been found in lung cancer research before and can be regarded as potential prognostic biomarkers of female LUAD.
In addition, GO and KEGG analyses were performed to reveal abundant molecular and biological pathways. The results showed that the most abundant GO terms were highly correlated with autophagy in biological process and cell composition. Consistent with GO analysis, KEGG analysis also enriched the most important autophagy pathway. Therefore, we speculate that this autophagy pattern may act as a tumor promoter in the occurrence and development of female LUAD.

Conclusions
In summary, our study identified the 12 autophagy-related genes (ITGA6, ERO1A, FKBP1A, BAK1, CCR2, FADD, EDEM1, ATG10, ATG4A, DLC1, VAMP7, ST13) that were associated with female LUAD survival. Based on the 12 genes, a risk marker was established and proven to be an independent prognostic factor for female patients with LUAD. Our results are expected to be applied in clinical practice, which may become a potential targeted therapy for female lung adenocarcinoma, and further research will provide us with its exact mechanism. However, the limitations of this study were retrospective analysis. The prospective studies should be carried out to verify the prognostic value of autophagy related genes, and the multicenter data will confirm our findings.