Identification of an autophagy-related 12-lncRNA signature and evaluation of NFYC-AS1 as a pro-cancer factor in lung adenocarcinoma

Objective: To develop an autophagy-related lncRNA-based risk signature and corresponding nomogram to predict overall survival (OS) for LUAD patients and investigate the possible meaning of screened factors. Methods: Differentially expressed lncRNAs and autophagy genes were screened between normal and LUAD tumor samples from the TCGA LUAD dataset. Univariate and multivariate Cox regression analyses were performed to construct the lncRNA-based risk signature and nomogram incorporating clinical information. Then, the accuracy and sensitivity were confirmed by the AUC of ROC curves in both training and validation cohorts. qPCR, immunoblot, shRNA, and ectopic expression were used to verify the positive regulation of NFYC-AS1 on BIRC6. CCK-8, immunofluorescence, and flow cytometry were used to confirm the influence of NFYC-AS1 on cell proliferation, autophagy, and apoptosis via BIRC6. Results: A 12-lncRNA risk signature and a nomogram combining related clinical information were constructed. Furthermore, the abnormal increase of NFYC-AS1 may promote LUAD progression through the autophagy-related gene BIRC6. Conclusion: 12-lncRNA signature may function as a predictive marker for LUAD patients, and NFYC-AS1 along with BIRC6 may function as carcinogenic factors in a combinatorial manner.


Introduction
Lung cancer accounts for 12.7% of the newly diagnosed cancer patients and up to 22.4% of deaths caused by cancer in 2020, according to the National Cancer Institute of the United States. About 50% of lung cancer cases belong to the lung adenocarcinoma (LUAD) type. Evaluation at early stages may help identify patients with a high-risk re-occurrence rate, develop a specific treatment strategy, and improve the overall survival rate of patients. Thus, the development of meaningful prognostic markers and the construction of related risk evaluation are essential for the diagnosis and treatment of LUAD patients.
Long noncoding RNA (lncRNA) is defined as a class of transcripts longer than 200 nucleotides (bps) and not translated into proteins (Iyer et al., 2015). LncRNA has been recently found to have many biological functions such as transcriptional regulation, translation, epigenetic modification, and cell fate decision (Ulitsky and Bartel, 2013;Li and Izpisua Belmonte, 2015). LncRNA has been demonstrated to be the driver of carcinogenesis, cancer progression, and metastasis in various cancers (Hirata et al., 2015;Yue et al., 2016;Yin et al., 2018;Tamang et al., 2019). For example, lncRNA-DANCR is overexpressed in tumor initiation cells, identifying DANCR as a prognostic biomarker and therapeutic target for hepatocellular carcinoma treatment (Yuan et al., 2016). Regarding LUAD, LncRNAs have been linked to various processes such as epithelial-mesenchymal transition (EMT), angiogenesis, and metastasis (Peng et al., 2018;Guan et al., 2019;Wang et al., 2020a).
Autophagy is an evolutionarily highly conserved catabolic process in eukaryotic cells (Han et al., 2014). It maintains the intracellular homeostasis of the environment by degrading and recycling the cellular components via the formation of autophagosome vesicles to engulf dysfunctional cytoplasmic organelles and the formation of autolysosome to degrade the contents of vesicles in the end (Klionsky, 2007;Fulda and Kögel, 2015). Autophagy plays a dual-faced role in tumor progression. On the one hand, it provides essential circulating metabolic substrates for biomolecule synthesis and supports the survival of cancer cells (Fulda and Kögel, 2015). On the other hand, autophagy could interplay with apoptosis to induce the death of cancer cells, including LUAD cells (Gewirtz, 2014;Liu et al., 2017). In addition, it plays an even more complex role in drug resistance, which depends on the activation status of autophagy and cellular environments . Recently, accumulating evidence has shown that lncRNAs regulate the autophagy network via transcriptional regulations of autophagyrelated genes in LUAD cells . Given the relevance of autophagy and lncRNAs in LUAD initiation and Frontiers in Genetics frontiersin.org 02 progression, this study focuses on investigating the expression patterns of autophagy-related lncRNAs from TCGA and Human Autophagy Database (HADb) to construct a useful risk signature to predict the prognosis of LUAD patients.

Materials and methods
Data collection and processing A detailed workflow for the construction of risk signature was developed, as shown in Figure 1. lncRNA and mRNA FPKM (fragments per kilobase of transcript per million fragments mapped) (level 3) sequencing profiles and associated clinical information of patients with LUAD were obtained from the TCGA data portal (https://tcgadata.nci.nih. gov/tcga/) before 1 November 2020. The GSE31210 microarray expression data were downloaded from the Gene Expression Omnibus (GEO) database, which included data of 226 LUAD samples, 20 normal tissue samples, and their associated clinical information. All autophagy-related genes were downloaded from the Human Autophagy Database (HADb) (http://www.autophagy.lu/). The "Limma" package of R software was used to screen differential autophagy-related gene expression and lncRNAs between adjacent and tumor tissues when | log 2 fold change (FC)| ≥ 0.5, and adjusted p < 0.05 was used as a filtering threshold. "Limma" package of R software was used to analyze the correlation between each differential lncRNA and every differential autophagy-related mRNA, and autophagy-related lncRNAs were screened when corFilter ≥0.3, and p < 0.05 was used as a filtering threshold.
Identification of autophagy-related lncRNA prognostic signature After the exclusion of LUAD patients with unavailable or insufficient follow-up clinical information, 458 LUAD patients were submitted for univariate Cox regression analysis. The "Survival" package of R software was used to perform the univariate regression analysis for the above lncRNAs, and prognosis-related lncRNAs were initially selected as candidates when p < 0.05. Next, the "Survival" package of the R software was used to perform multivariate COX regression analysis, and interdependent lncRNAs were identified as the risk signature, further validated by the Kaplan-Meier survival analysis. The risk score of each patient was based on a linear combination of lncRNA expression level (Exp) multiplied by a regression coefficient (β), represented as the following: Risk score n i 1 Expi βi.
Then, based on the median risk score value, patients in the cohort were divided into high-risk and low-risk groups.

Visualization of expression correlation network of autophagy-related genes and lncRNAs
After lncRNAs were identified as the risk signature, Cytoscape of version 3.6.1 was used to generate the expression correlations network between autophagy-related genes and lncRNAs when corFilter ≥ 0.3. The relations among high-or low-risk, autophagy-related genes and lncRNAs were also illustrated as a Sankey diagram to show the possible influence of autophagy-related lncRNAs on clinical outcomes.

Risk signature assessment
After the risk signature was constructed, the Kaplan-Meier survival analysis was first used to evaluate the overall survival between high-and low-risk groups. Furthermore, to evaluate the predictive ability of the constructed risk signature, the timedependent receiver operating characteristic (ROC) analyses were performed by the "Survival ROC" package of R software to test the signature's sensitivity and specificity.
Development and validation of a nomogram incorporating the lncRNA signature with clinical factors Next, based on univariate and multivariate Cox regression results, a novel nomogram incorporating the 12-lncRNA signatures and related clinical factors was developed by the "rms" package of the R software. The predictive ability of the risk signature was further evaluated with the area under the curve (AUC) in the ROC analysis. "Caret," "foreign," and "survival" packages of R software were also used to randomly select half of the samples from the TCGA LUAD dataset as the internal cohort for validation. All the analyses were performed in the TCGA internal cohort and GSE dataset.

Gene set enrichment analysis
Based on the median risk scores, 551 LUAD samples in the TCGA cohort were divided into two groups (high-and low-risk groups). Then, gene set enrichment analysis (GSEA) was performed to identify the significantly differentially expressed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between the high-and low-risk groups using the Java-based GSEA software. A nominal (NOM) p < 0.05 was chosen as the cutoff value.
Quantitative reverse transcription polymerase reaction assays for lung cancer belong to lung adenocarcinoma patient specimens and further validation in GSE40791 datasets Ten paired LUAD tumor and adjacent normal tissues were obtained from patients of the First Affiliated Hospital of Anhui University of Science and Technology. All these patient tissues were transferred on ice and stored in liquid nitrogen. Total RNA was extracted from tissues with TRIzol reagent (Invitrogen, China) according to the manufacturer's guide manual. Reverse transcription was performed according to the manufacturer's instructions using the ReverTra Ace qPCR RT Kit (Toyobo, Japan). The SYBR Green Realtime PCR Master Mix (Toyobo) was used in the quantitative reverse transcription-polymerase chain reaction (qRT-PCR) experiment. Gene expression level was calculated with the 2-DDCt method. Primers in the qRT-PCR test are listed in Supplementary Table S1. Informed consent was signed by all the participants. The study was approved by the Ethics Committee of Chongqing Medical University. To further validate the expression patterns of these lncRNAs, we use GSE40791 datasets, which included 100 non-neoplastic (N) lung samples, and 69, 12, and 13 stages I, II, and II lung adenocarcinoma (AD) frozen tissues, respectively. The box plot is implemented by the "ggplot2" package of R software based on the Wilcoxon rank sum test.
Cell culture, plasmid construction, and transfection A549 cells were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China) and cultured in Dulbecco's MEM medium (Hyclone, Los Angeles, United States) supplemented with 10% fetal bovine serum and a 1% penicillin-streptomycin solution. All cells were incubated at 37°C in a humidified incubator containing 5% CO 2 . pEGFP-BIRC6 overexpression plasmid was constructed as previously described (Zhang et al., 2020a) by replacing it with the BIRC6 CDS region. For knockdown of NFYC-AS1, the specific shRNAs were synthesized and constructed into the same vector as previously described (Zhang et al., 2020a) by replacing targeted sequences with "CACCTGTAA TCCCAGCACTTT" (sh1) or "CACACCTGTAATCCCAGCATT" (sh2). Assayed cells were transfected with Lipofectamine ® 2000 Reagent (Invitrogen) following the supplier's instruction manual.

Cell proliferation and apoptosis assays
Cells were seeded in 96-well plates at a density of 2,000 cells per well. The growth rate of cells was evaluated

Primer
Sequence Frontiers in Genetics frontiersin.org using the CCK-8 cell proliferation kit (Dojindo Laboratories, Kumamoto, Japan), according to the manufacturers' instructions. OD detection at 450 nm was carried out by infinite 200Pro (Tecan). Cell apoptosis was analyzed using the Annexin V/7-AAD Apoptosis Detection kit (Keygen Biotech, Nanjing, China) on a CytoFLEX flow cytometer (Beckman, California, United States). Results were further analyzed by Flowjo 10.4.

Western blotting, immunofluorescence assay, and immunohistochemistry
For immunoblotting experiment, cells were lysed in RIPA buffer (50 mM Tris (pH 7.4), 150 mM NaCl, 1% Triton X-100, 1% sodium deoxycholate, 0.1% SDS). An equal amount of protein was loaded and simultaneously subjected to electrophoresis in SDS-polyacrylamide gel and transferred to 0.22 μm pore-sized PVDF membranes (Roche, Basel, Switzerland). Membranes were briefly blocked with 5% skim milk and incubated with the primary antibodies overnight at 4°C, followed by incubation with the speciesmatched secondary antibody conjugated with HRP (Cell Signaling Technology, Danvers, United States) for 1 h at room temperature prior to chemiluminescence detection. For the immunofluorescence (IF) experiment, cells with less than 50% confluency were seeded and grown on coverslips overnight. Cells were fixed with 4% paraformaldehyde and permeabilized with Triton-100. After blocking in 5% goat serum, samples were incubated with the appropriately diluted primary and secondary antibodies. Eventually, cells were observed and photographed by fluorescence microscopy (Nikon C2). For the immunohistochemistry experiment, all tissues were fixed with 4% paraformaldehyde overnight at 4°C and then embedded in paraffin. The samples were subsequently sectioned into thin slices to be mounted on slides, followed by deparaffinization in xylene and rehydration through a series of ethanol-water solutions. Antigen retrieval was carried out by immersing the sections in citrate acid buffer with heating by a microwave oven. Slides Frontiers in Genetics frontiersin.org

FIGURE 4
The Kaplan-Meier survival analysis of 12 screened lncRNAs.

FIGURE 5
The interaction network of autophagy genes and OS-associated lncRNAs in LUAD patients. (A) The co-expression network of the autophagyrelated mRNAs and lncRNAs was constructed and visualized using Cytoscape. (B) The co-expression network with risk-type information was visualized using the Sankey diagram.

Immunodeficient xenograft mouse tumor model
BALB/c-nude mice (4-5 weeks of age, 18-20g) were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (China). All experimental procedures were approved by the Institutional Animal Care and Use Committee of Chongqing Medical University. Briefly, cells were trypsinized, washed, and resuspended in phosphate-buffered saline at a density of 10 7 cells/ml. A total of fifteen mice were inoculated subcutaneously with 5×10 6 cells of A549 Vector cells, A549 NFYC-AS1 knockdown cells, and A549 NFYC-AS1 knockdown with BIRC6 overexpression cells, respectively. Mice were sacrificed 7 days after tumor cell implantation, and the tumors were removed and weighed. GraphPad Prism 8 and unpaired t-test were used for statistical analysis.

Results
Differential expression of autophagyrelated genes and lncRNAs in lung cancer belong to lung adenocarcinoma patients A detailed data processing flowchart is shown in Figure 1. Clinical and diagnostic information of LUAD patients in the TCGA and validation cohorts are shown in Table 1. A total of 232 autophagy genes were collected from the HADb database, and we identified 76 differentially expressed autophagy genes between 497 tumor tissues and 54 normal tissues from TCGA at the criteria of | log 2 (FC)| ≥ 0.5 and an adjusted p < 0.05 (Figure 2A). A total of 507 differentially expressed lncRNAs were identified, including 373 upregulated and 134 downregulated lncRNAs ( Figure 2B). Then, a total of 107 lncRNAs related to autophagy were screened using the "Limma" package of R software when the criteria were set to corFilte ≥ 0.3 and p < 0.05 (Supplementary File S1).
Identification of autophagy-related lncRNA prognostic signature for lung cancer belonging to lung adenocarcinoma Univariate Cox regression analysis identified 28 out of 107 lncRNAs as prognostic factors ( Figure 3A). Next, multivariate Cox regression was performed and 12 lncRNAs were further found to be independent factors in predicting the prognosis of LUAD patients ( Figure 3B). Results of the Kaplan-Meier survival analysis further confirmed the relevance between these 12 lncRNAs and OS (Figure 4). These 12 lncRNAs are NFYC-AS1, SNHG10, HCG18, LINC00857, LINC01116, LINC00996, CRNDE, CASC15, TMPO-AS1, LINC00654, LINC01138, and ZNF790-AS1. An expression correlation network was built on the 12 OS-associated lncRNAs and 16 corresponding associated autophagy-related genes, including FOXO1, FOXO3, MLST8, MTOR, BIRC6, RPS6KB1, CFLAR, IFNG, DAPK2, ITPR1, MBTPS2, ULK2, ATIC, GAPDH, CASP1, NLRC4, PRKCQ, ITGB4, SPHK1, (0.537*EXP_LINC01138) + (0.372*EXP_ZNF790-AS1). A total of 458 LUAD patients were divided into high-and low-risk groups based on the median risk score. Kaplan-Meier survival analysis showed a significant OS advantage of the low-risk group over the high-risk group, as shown in Figure 6A (p < 0.001). Moreover, more than 0.7 of 1-, 3-, and 5-year area under the curve (AUC) in receiver operating characteristic (ROC) analysis exhibited a satisfactory diagnostic performance ( Figure 6A). We further validated this risk signature in both TCGA internal cohort and the GSE validation dataset. Once again, Kaplan-Meier survival analysis showed an improved OS in the low-risk group compared to the high-risk group, as shown in Figures 6B,C. 0.7 of 1-, 3-, and 5-year AUC in ROC analysis also showed a fair sensitivity and specificity of this risk signature. Then, based on the multivariate Cox proportional hazards regression analyses, a nomogram including 12-lncRNA signatures was constructed by the ''rms'' package in R

FIGURE 9
Nomogram incorporating risk score and other clinical information. Top, nomograms to predict 3-or 5-year OS of patients with LUAD. Bottom, calibration curves of the nomogram for 3-or 5-year OS prediction.
Frontiers in Genetics frontiersin.org software to help predict the 3-or 5-year survival of LUAD patients ( Figure 6D).

Nomogram development and validation
Univariate and multivariate Cox proportional hazards regression analyses were performed to evaluate the predictive value of different clinical parameters. Univariate Cox regression analysis showed that stage, T status, N status, and risk score were significantly associated with OS of LUAD patients in TCGA datasets, as shown in Figure 7A (p < 0.001). On the other hand, multivariate Cox regression analysis only showed that stage and risk score are the independent predictive factor of OS for LUAD patients ( Figure 7B). In the GSE31210 dataset, univariate Cox regression analysis showed that the stage and risk score are significantly related to the OS of LUAD patients ( Figure 7C). Multivariate Cox regression showed that stage and risk score are the independent predictive factors of OS ( Figure 7D). Also, risk signature based on 12 lncRNA showed 1-, 3-, and 5-year AUC of above 0.69 in ROC analysis along with other clinical factors in both TCGA and GEO datasets (Figure 8). Then, we constructed a nomogram incorporating 12-lncRNA signatures and other clinical information by the ''rms'' package in R software to predict the 3-and 5-year survival of LUAD patients. The following calibration curve showed a satisfactory consistency between the predictive and observed results (Figure 9).

Validation of the expression of lncRNAs in lung cancer belong to lung adenocarcinoma patient samples and GSE40791
The expression of selected lncRNAs from our constructed signatures (HCG18, TMPO-AS1, NFYC-AS1, LNC00996) was investigated in 10 paired tumor and adjacent tissues from LUAD patients using qRT-PCR. The results indicated that HCG18, TMPO-AS1, and NFYC-AS1 were significantly upregulated in tumors compared with adjacent samples. In contrast, LNC00996 were significantly downregulated in tumors compared with adjacent samples (Figure 10) (p < 0.05). These results were consistent with the expression pattern from the TCGA LUAD dataset (Supplementary Figure S1). Also, results from the GSE40791 dataset, including 100 normal lung and 94 LUAD lung tissues, further validate the expression pattern of these lncRNAs (Supplementary Figure S2).

Frontiers in Genetics
frontiersin.org lncRNA NFYC-AS1 promotes lung adenocarcinoma cell proliferation through BIRC6 To investigate the possible connection and its influence between lncRNAs and related autophagy genes, we selected two potential players, NFYC-AS1 and BIRC6, based on previous literature  and their high degree of correlation ( Figure 5; Supplementary Figure S3A). Correlation analysis from GEPIA confirmed a positive correlation between the NFYC-AS1 and BIRC6 expression in TCGA LUAD dataset (p < 0.001, R = 0.35) but a relatively weak correlation in GTEx lung dataset (p < 0.05, R = 0.15) (Supplementary Figures S3B,C). Immunoblot results from lung adenocarcinoma cell line showed that knockdown of NFYC-AS1 did reduce the BIRC6 protein level (Figures 11A,C), suggesting a possible regulation of BIRC6 expression by NFYC-AS1. To determine the subsequent effect of this regulation, we explored the influence of NFYC-AS1 and BIRC6 on A549 cells. Results show that the knockdown of NFYC-AS1 clearly inhibits cell growth ( Figure 11B), whereas the restoration of BIRC6 expression reverses the inhibitory effect of NFYC-AS1 knockdown ( Figure 11D), suggesting that the pro-proliferative effect of NFYC-AS1 is via the regulation of BIRC6. More importantly, IHC of BIRC6 from 34 LUAD patients showed significantly lower expression in NFYC-AS1 low or normal lung tissues and a markedly elevated expression in tumor tissues with high NFYC-AS1 expression ( Figures 12A,B,C). In addition, mRNA expression pattern of BIRC6 from the TCGA database was consistent with the above results from patient samples ( Figure 12D), both indicating that NFYC-AS1 may promote tumor progress and act as a high-risk factor of poor prognosis in our constructed signature through the upregulation of BIRC6.
Given that BIRC6 is known to be the regulator of both autophagy and apoptosis (Ebner et al., 2018), we investigated the influence of this regulation on cell autophagy and apoptosis. Results reveal that the knockdown of NFYC-AS1 indeed alters the expression of some autophagy-related markers, including LC3B, Beclin 1, and SQSTM1/p62 in A549 cells (Figures 13A,B). Furthermore, the knockdown of NFYC-AS1 also induced evident apoptosis in A549 cells ( Figures 14A,B), whereas BIRC6 overexpression seemingly offset this effect ( Figures  14C,D). More importantly, we performed in vivo tumor xenograft experiment in nude mice and found that the knockdown of NFYC-AS1 did exhibit the inhibitory effect on  Figure S5).

Discussion
As the most common pathological subtype of lung cancer, the proportion of lung adenocarcinoma (LUAD) is up to 40%-70% of all lung cancer patients. Although certain progress has been made in the research and practice of LUAD treatment, it is still one of the most difficult types of cancers to treat due to its highly metastatic and malignant nature. The dysregulation of autophagy is thought to play an important role in lung cancer. For example, in a mouse model of lung cancer, upregulation of autophagy via caspase-3 and mTOR inhibition promotes the efficacy of radiotherapy. C-myc/miR-150/EPG5 axis mediated dysfunction of autophagy induced increased cellular ROS levels Frontiers in Genetics frontiersin.org and DNA damage response and promoted NSCLC development (Li et al., 2019a). Also, autophagy could protect LUAD cells to Src inhibitors, whereas microRNA-106a could target autophagy kinase ULK1 and compromise this protective effect (Rothschild et al., 2017). In addition, autophagy inhibition of cancer stem cells could promote the efficacy of cisplatin against NSCLC in both A549 cell lines and NOD/SCID mice . LncRNA, along with miRNA and mRNA expression, has been used to develop risk signatures to prevent and predict the development of various tumors. Increasing evidence has shown significant correlations between lncRNAs and autophagy in lung cancer. For instance, the lncRNA BLACAT1 promoted ATG7 expression through miR-17, facilitated autophagy, and promoted the chemoresistance of NSCLC cells through the miR-17/ATG7 axis (Huang et al., 2019). lncRNA CASC2 inhibited autophagy and promoted apoptosis in NSCLC cells via the miR-214/TRIM16 signaling pathway (Li et al., 2018a). LncRNA MSTO2P promotes lung cancer cell proliferation and autophagy by upregulating EZH2 . Moreover, lncRNA LCPAT1 mediates smoking or PM 2.5induced autophagy and epithelial-mesenchymal transition via RCC2, implying the possible role of lncRNA in rendering lung cancer cells into a more invasive state (Lin et al., 2018).
In this study, we identified and validated an autophagyrelated 12-lncRNA risk signature that was highly associated with the OS of patients with LUAD. The signature showed accuracy and robustness by calculating the AUC in the ROC analysis. Multivariate Cox regression analyses suggested that age, N stage, and the risk signature were independent risk factors for OS in the primary cohort. A novel nomogram incorporating clinical factors was constructed and validated to predict the prognosis for LUAD patients. The expression levels of lncRNA signatures were also validated in the clinical LUAD samples by qRT-PCR. These results suggested that both the 12-lncRNA risk signature and the novel nomogram were effective prognostic indicators in patients with LUAD.
Compared to other members in the risk signature, NFYC-AS1 received relatively little attention, except that it has been demonstrated as a prognostic biomarker in the miRNA-lncRNA-mRNA regulatory network in LUAD patients (Li et al., 2016). BIRC6 (Baculoviral IAP repeat-containing 6) is an IAP (inhibitor of apoptosis proteins) that not only plays an anti-apoptotic role through the inhibition of pro-apoptotic protein SMAC or caspases (Hauser et al., 1998;Hao et al., 2004;Qiu and Goldberg, 2005;Ren et al., 2005) but also is Frontiers in Genetics frontiersin.org associated with other cellular processes; for instance, it regulates autophagosome-lysosome fusion in autophagy (Ikeda, 2018;Jia and Bonifacino, 2019). In addition, overexpression of BIRC6 has been found in various tumors, including colorectal cancer (Bianchini et al., 2006), gastric carcinomas (Salehi et al., 2016), and lung cancer (Dong et al., 2013). In our study, the NFYC-AS1 expression shows a high positive correlation with autophagy genes, especially with BIRC6 (p < 0.001, R = 0.35), and the knockdown of NFYC-AS1 decreased the BIRC6 expression and induces both autophagy and apoptosis in A549 cells, whereas the rescue of BIRC6 reverses the inhibitory effect (Figures 13, 14). Although there is no evidence showing a direct regulation of BIRC6 by NFYC-AS1 until now, combining all the above results, it is possible that NFYC-AS1 enhances BIRC6 expression via binding to specific proteins or regulating miRNAs, which in turn leads to the inhibition of autophagy, apoptosis, and the progression of tumors. Gene set enrichment analysis (GSEA) further revealed that the constructed risk signature may regulate several KEGG pathways to influence the progression of LUAD, including the ''p53 signaling,'' "proteasome," "protein export," and "pentose phosphate" pathways. Tumor suppressor p53 acts as the guardian of the genome, which senses the stress in the cellular environment, initiates DNA repair, or leads to the cell cycle arrest. Recently, growing evidence has revealed that p53 can activate the transcription of autophagy-related genes, whereas autophagy could suppress p53 via the regulation of ROS or AMPK. In contrast, proteasome is considered responsible for proteolysis and organelle homeostasis by interconnecting with autophagy pathways. As a crucial part of autophagy, lysosomes carry out the degradation of extracellular particles from endocytosis and of intracellular components from autophagy, which contains more than 50 membrane proteins that control the specificity and timing of cargo influx and protein export. The pentose phosphate pathway (PPP) plays a key role in facilitating tumor cells with the glycolytic process and countering the damage of reactive oxidative species. G6PD is the limiting enzyme of the PPP associated with cancer progression and drug resistance, the inhibition of which indirectly leads to the dysfunction of autophagy. Although GESA analysis implicated the involvement of lncRNAs and autophagy in LUAD tumor progression, the detailed association of how autophagyrelated tumor progression and lncRNAs are included in our signatures should be further investigated.
This study has some obvious advantages. First, it included 551 samples from the LUAD TCGA dataset and 246 samples from GSE31210, providing a sufficient sample size to avoid statistical bias and make a satisfactory result. Indeed, the constructed signature and incorporated nomogram were further validated in both TCGA train and GSE datasets and showed a good specialty and sensitivity. Second, clinical samples further confirmed that the expressions of several lncRNAs are consistent with the results from the database. Then, a novel nomogram incorporating the risk signature and related clinicopathological factors is constructed and provides a helpful prediction for the prognosis of LUAD patients. Third, we further investigated the possible regulation of BIRC6 by NFYC-AS1 and related influences on lung cancer cells, proving that NFYC-AS1 may act not only as an indicator for poor prognosis but also as a potential driver of carcinogenesis through the negative modulation of autophagy and apoptosis. Although, in current experimental settings, it is difficult to distinguish which effect is dominant or whether BIRC6 is directly responsible for this inhibitory effect through autophagy, all the above results suggest the importance of NFYC-AS1 as a pro-cancer factor in LUAD patients.
In the meantime, there are some limitations to the current study. First, datasets from GEO or TCGA have various kinds of lncRNAs due to their unique sequencing methods. Thus, on the one hand, using both GEO and TCGA datasets may increase the sample size and help the validation of risk signatures. On the other hand, the overlap of two different datasets may neglect some lncRNAs that are significantly deregulated in one database but not sequenced at all in another. Second, the TCGA and GSE31210 databases did not cover important information such as medical history, treatment strategy, and family history, which could alter the prognosis of LUAD patients. Third, although our results showed a possible regulation of NFYC-AS1 on BIRC6 and BIRC6 expression indeed increased in tumor tissues (Supplementary Figure S4), the direct evidence of this regulation still needs to be further discovered. Last, in the current study, we investigated the expression of four lnRNAs in 10 patients' samples. To assure the accuracy of the expression patterns, clearly more clinical samples should be analyzed.

Conclusion
Collectively, our current study analyzed autophagyrelated lncRNAs in LUAD patients comprehensively and constructed a lncRNA-based risk signature. Our results revealed that high-risk scores are associated with a worse prognosis in LUAD patients. Our study also suggested that the 12 autophagy-related lncRNAs have significant value in predicting the prognosis of LUAD patients. Furthermore, NFYC-AS1 and BIRC6 may be potential therapeutic targets due to previous research and their significance in our study. More specific experiments and research based on clinical samples should be carried out to validate the results of our study.
Frontiers in Genetics frontiersin.org

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Ethics statement
The studies involving human participants were reviewed and approved by the Ethics Committee of Chongqing Medical University. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

Author contributions
TF and XL performed most of the research work. SX collected clinical samples and carried out some of the experiments. MZ designed the project and wrote the manuscript. All authors reviewed and approved the final manuscript for submission.

Funding
This work was supported by the Xinglin program of Chongqing TCM/TCM-integrated Key discipline.