The pancreatic adenocarcinoma RNA-seq data and corresponding clinical information were downloaded from the TCGA dataset¹. The cohort contains 178 tumor tissues and four normal pancreatic tissues. And, 177 PC patients with complete clinical information were extracted for further analysis. Perl language was performed to convert gene names from Ensemble IDs to a profile of gene symbols with the Ensemble database².

Autophagy-related genes (ARGs) were obtained from the Human Autophagy Database³. We extracted the lncRNA list from mRNA expression data of the GNECODE project⁴. Then, the Pearson correlation was applied to analyze the correlation between the lncRNAs and ARGs. The lncRNAs with correlation coefficient |R²| > 0.5 and p < 0.01 were considered as autophagy-related lncRNAs (ARlncRNAs).

To increase the reliability of our study, we randomly divided the entire dataset into a training set (accounting for 60%) and a validation set (accounting for 40%) by the “caret” R package (version 6.0-84)⁵ (Deist et al., 2018). And the whole dataset was considered as an entire set (n = 177). At first, we adopted the univariate cox regression analysis to identify the significant ARlncRNAs in the training set with a p < 0.01 by the “survival” R package. Then, the least absolute shrinkage and selection operator (LASSO) penalized cox regression analysis was performed to further reduce the dimension and the multivariate cox regression analysis was utilized to calculate the risk coefficients of the prognostic signature. The risk score formula is shown as follows: Risk score = ΣCoef ARlncRNAs × Exp ARlncRNAs. The Coef ARlncRNAs represents the coefficient of each ARlncRNAs and Exp ARlncRNAs is the expression of each ARlncRNAs. Based on the median risk score of the signature, the patients were divided into low-risk and high-risk groups. The survival analysis for the different groups was performed using the Kaplan-Meier (K-M) survival curve analysis and log-rank test analysis with the “survminer” R package. Moreover, we constructed the receiver operating characteristic (ROC) curve by using the “survivalROC” R package to evaluate the specificity and sensitivity of the prognostic signature.

In order to simplify the predictive model, we created a nomogram based on independent clinical prognostic factors with the “rms” R package (Iasonos et al., 2008). We plotted the calibration curve of the nomogram to value the predictive power of the prognostic signature.

Grouped samples and expression patterns were analyzed using the principal component analysis (PCA). Gene set enrichment analysis (GSEA) was performed to evaluate different functional phenotypes between low- and high-risk groups (Subramanian et al., 2005). Moreover, we analyzed the correlation between different risk groups and clinical characteristics with the chi-square test and the results were presented in a heat map.

To better investigate the relationship between the signature and immune cell infiltration, we calculated the infiltration expression of 22 immune cells in PC by using the “CIBERSORT” R package. Then, the immune-related Pearson correlation coefficients were tested for relevance in the R program.

Moreover, to explore the clinical utility of the signature to predict therapeutic effect, we used the Pearson’s analysis to calculate the correlation between risk scores with molecules of targeted therapy. The therapy targets are as follows: programmed cell death 1 (PD-1, also known as PCDC1), programmed cell death ligand 1 (PD-L1, also known as CD274), epidermal growth factor receptor (EGFR), vascular Endothelial Growth Factor Receptor 3 (VEGFR3, also known as FLT4), KIT proto-oncogene (KIT), Fms-like tyrosine kinase 3 (FLT3), MET proto-oncogene (MET), vascular Endothelial Growth Factor Receptor (VEGFR1, also known as FLT1), mammalian target of rapamycin (mTOR), platelet-derived growth factor receptor alpha (PDGFRA), and platelet-derived growth factor receptor beta (PDGFRB).

The DIANA online tools⁶ were employed to explore the miRNAs binding to lncRNA. We employed three miRNA databases, including miRDB⁷, miRTarBase⁸, and TargetScan⁹, to predict the target genes of miRNAs. To predict the expression correlation between lncRNAs and miRNAs, the threshold was set at 0.9. Subsequently, the lncRNA-miRNA-mRNA regulatory network was mapped by the Cytoscape (version 3.7.0)¹⁰ to better understand the connections.

To further explore the functional annotation and pathway analysis of the target mRNAs, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed by using the “clusterProfiler” R package with a p < 0.05, and FDR < 0.05.

Connectivity Map (CMap)¹¹ is a collection of genome-wide transcriptional expression data from cultured human cells treated with bioactive small molecules and analyzed by corresponding matching algorithms to investigate the relationship between drug and gene expression changes and phenotypes (Lamb et al., 2006). We uploaded up- and down-regulated target genes from the lncRNA-miRNA-mRNA network to CMap. A connectivity score ranging from -1 to 1 was used to reflect the degree of closeness between the expression spectrums. The drugs with negative scores were potential therapeutic molecules. Moreover, these candidate drugs were investigated in the Pubchem database¹².

Gene Expression Profiling Interactive Analysis (GEPIA)¹³ is a website for large-scale expression analysis and interactive analysis that has been used to compare the expression of signature lncRNAs (Tang Z. et al., 2019).

Blood samples of PC patients were collected from the Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology (Wuhan, China) between May 2019 and July 2019. Patients were eligible if they didn’t receive any preoperative radiation and chemotherapy and the postoperative pathology was officially diagnosed as pancreatic adenocarcinoma. Exclusion criteria are as follows: (1) patients with history of previous cancer; (2) patients with multiple tumors or PC is not a primary lesion; and (3) patients with co-morbidities of the blood system. Finally, 30 of the 37 blood samples were eligible for further study. And, 10 blood samples of healthy donors were collected as a control group.

The serum specimen was separated at 3,000 rpm for 10 min from the venous blood. All the serum samples were stored at −80°C. Ethical approval for the use of human samples was obtained from the Tongji Hospital Research Ethical Committee.

Total RNA was extracted from serum samples and cells by the TRIzol reagent (Life Technologies, Thermo Fisher Scientific, United States). The complementary DNA was synthesized with the PrimeScript^TM RT Master Mix (Takara Bio Inc, Dalian, China) according to the manufacturer’s instructions. Quantitative real-time PCR (qRT-PCR) was performed using a SYBR Green PCR kit (Thermo Fisher Scientific) following the standard protocol. And GAPDH served as the internal control. The forward primer for LINC01559 was 5′-GTCCTGCAGAACTCCCTCTT-3′, the reverse primer for LINC01559 was 5′-AGTCCTGGAGCTGCAGAAAT-3′. The forward primer for AC245041.2 was 5′-TTGCCCCCATCTTTGCCATTCC-3′, the reverse primer for AC245041.2 was 5′- TTGACCCATCTTTCCTCCCCAC-3′. The forward primer for AC005332.6 was 5′-AAGACAGCACG GTGTTAAAAAG-3′, the reverse primer for AC005332.6 was 5′-TTGAATCCAGGAGGCGGAAG-3′. The forward primer for GAPDH was 5′-GGAGCGAGATCCCTCCAAAAT-3′, the reverse primer for GAPDH was 5′-GGCTGTTGTCATACTT CTCATGG-3′. The relative expression was calculated by the 2^–ΔΔCt method.

Human PC cell lines PANC-1 and SW1990 were obtained from the American Type Culture Collection (Manassas, VA, United States). These cells were maintained in the Dulbecco’s modified Eagle medium (Gibco, Carlsbad, CA, United States) supplemented with 10% fetal bovine serum (FBS; Gibco), 100 units/mL of penicillin, and 100 mg/mL of streptomycin (Sigma-Aldrich, St. Louis, MO, United States). All cell lines were authenticated, mycoplasma-free and cultured at 37°C in a humidified incubator containing 5% CO₂.

LINC01559 si-1, LINC01559 si-2 and si-NC were purchased from the DesignGene Biotechnology (Shanghai, China) and transfected into PC cells using the Lipotransfectamine 3000 (Thermo Fisher Scientific). The lentiviral vector containing the tandem-labeled GFP-mRFP-LC3 reporter were also constructed by the DesignGene Biotechnology (Shanghai, China), and transfection was carried out according to the manufacturer’s specification.

Western blotting assay was performed to detect the expression of LC3B, p62, and GAPDH as previously described (Tian et al., 2018). The antibodies of LC3B (#2775, 1:200), p62/SQSTM1 (#5114, 1:500), cleaved caspase-3 (#9661, 1:1,000), cleaved PARP (#9544, 1:1,000), and GAPDH (#5174, 1:1,000) were purchased from the Cell Signaling Technology (CST, Danvers, MA, United States). And the intensity of bands was estimated by the Image J2X (National Institute of Mental Health, Bethesda, MD, United States). All experiments were repeated three times.

For migration assay, 1 × 10⁵/mL PC cells were suspended into the upper transwell chamber of 24-well transwell plates (8 μm pore size; Corning) containing 200 μL serum-free medium, while the lower chambers were supplied with a 500 μL complete culture medium. After 48 h co-culture, the cells on the lower surface of membrane were fixed in 4% paraformaldehyde and stained with crystal violet solution. The stained cells were then counted under a Nikon light microscope (Nikon, Japan). For invasion assay, the upper transwell chambers were coated with 60 μL Matrigel matrix gel (BD Biosciences, United States). The other operations were the same as the transwell migration experiment. All experiments were repeated three times.

Indicated PC cells (2 × 10⁵ cells/well) were seeded in 6-well plates to grow to 90% confluence, and then we scratched the wound with a 200 μL pipette tip across the center of the well. After washing three times with PBS, the cells were incubated in a serum-free medium at 37°C with 5% CO₂. Wound healing was observed with the optical microscope (Nikon, Japan) at 0 and 24 h, respectively. All experiments were repeated three times.

A Cell Counting Kit-8 (Dojindo Laboratories Co. Ltd, Kumamoto, Japan) assay was used to evaluate cell viability. Briefly, PANC-1 and SW1990 cells were seeded in 96-well plates at a density of 2 × 10³ cells per well. Each group had triplicates (n = 3). Cells were treated with gemcitabine at indicated concentration after 48 h. At the indicated time point, 10 μL of CCK-8 solution was added into each well and cells were incubated for 2 h at 37°C with 5% CO₂. Then, the absorbance was measured at 450 nm with a plate reader (Bio-Tek Elx 800, United States). To assess cell proliferation, PC cells (1 × 10³ cells/well) were placed into 96-well plates and measured at 0, 24, 48, 72, 96. All experiments were repeated three times.

The percentage of apoptotic cells was analyzed by the PE Annexin V Apoptosis Detection kit (BD Pharmingen). Briefly, 1 × 10⁵/mL cells indicated PC cells were seeded into six-well culture plates. After 48 h of relevant treatment, the cells were harvested using trypsin without EDTA and washed twice with binding buffer. Finally, the cells were resuspended in 100 μL of binding buffer containing 5 μL Annexin V-PE and 5 μL 7-ADD for 15 min in the dark at room temperature, and the apoptosis was detected by flow cytometry (BD Biosciences, Franklin Lakes, NJ, United States). All experiments were repeated three times.

Indicated PC cells (1 × 10³ cells/well) were seeded into 6-well plates and cultured in completed medium at 37°C with 5% CO₂ for 2 weeks. Then, the cells were cleaned, fixed and dyed, and finally photographed. All experiments were repeated three times.

Cells transfected with GFP-mRFP-LC3B were grown on glass coverslips. Following the indicated treatments, cells were fixed with 4% formaldehyde for 30 min and photographed under a confocal microscope (Carl Zeiss, Germany, LSM710).

All statistical analysis and plotting were performed in the R language (Version 3.6.2). Univariate and multivariate Cox regression analyses were used to identify independent prognostic factors for PC. And p < 0.05 was considered statistically significant.

A total of 13,482 lncRNAs was extracted from the TCGA dataset, 825 of which were identified as ARlncRNAs by the Pearson correlation analysis (|R| > 0.5, p < 0.01).

Firstly, all patients were divided into training cohort (n = 107) and validation cohort (n = 70) (Table 1). Then, we employed the LASSO penalized cox regression analysis by the training cohort and found 10 more representative ARlncRNAs: AC245041.2, AL354892.2, FLVCR1.DT, AC125494.2, AL162274.2, LINC01559, AC090114.2, SH3PXD2A.AS1, AC005332.6, and AC092171.2 (Figure 1A). Moreover, the stepwise multivariate Cox regression was utilized to establish a predictive signature for PC patients in the training cohort with a risk score = (0.319702425 × expression level of AC245041.2) + (−0.934877496 × expression level of AC125494.2) + (0.038664123 × expression level of LINC01559) + (−0.594425726 × expression level of AC090114.2) + (−0.110425977 × expression level of AC005332.6) + (−0.184537572 × expression level of AC092171.2) (Figure 1A).

Next, the patients in the training cohort were divided into high-risk group and low-risk group based on the median risk score. Figure 1B showed the distribution of prognostic signature, survival outcomes of PC patients in different groups, and the expression profiles of the selected lncRNAs. Notably, the Kaplan-Meier survival analysis in the training cohort revealed that the survival time of PC patients was significantly longer in the low-risk group than the high-risk group (Figure 1C). As shown in Figure 1D, the area under the ROC (AUC) for 1, 3, and 5 years of the survival were 0.938, 0.890, and 0.804 in the training cohort, suggesting that the signature exerted a certain potential property for prognostic prediction in PC patients.

To verify the accuracy of the signature, we analyzed its prognostic value in the validation cohort and entire cohort. LncRNAs expression profiles, risk distribution, and survival rate in the validation cohort and entire cohort were shown in Figures 2A,B. Similar to the results in the training cohort, the Kaplan-Meier (K-M) survival analysis in both the validation cohort and entire cohort indicated that the survival outcome of PC patients was better in low-risk group than in high-risk group (Figures 2C,D). The AUC at 1, 3, and 5 years were and 0.848, 0.677, 0.737, and 0.921, 0.808, 0.774 in the validation cohort and entire cohort, respectively (Figures 2E,F).

Furthermore, the univariate and multivariate Cox regression analysis were employed to confirm the independent prognostic role of the signature for PC patients in the entire cohort (Figures 3A,B). Another independent prognostic factor was age. Moreover, a nomogram based on the signature risk score and clinical features was constructed and the calibration curve for 1, 3, and 5 years of the nomogram showed a great predictive power of the prognostic signature (Figures 3C,D).

PCA was employed to demonstrate the significant distribution difference between low- and high-risk groups based on the risk scores (Figure 4A). Then, GSEA was implemented to explore the significant enriched pathways between the two groups. As shown in Figure 4B, the top five up-regulated and down-regulated KEGG pathways were the “mTOR signaling pathway,” “calcium signaling pathway,” “regulation of autophagy,” “RNA polymerase,” “lysine degradation,” and “pentose phosphate pathway,” “starch and sucrose metabolism,” “glycolysis gluconeogenesis,” “glycosphingolipid biosynthesis lacto and neolacto-series,” and “pentose and glucuronate interconversions,” respectively. These results indicated that the high-risk score was significantly associated with autophagy regulation and several signaling pathways may participate in the process.

To investigate the clinical utility of the signature, we explored the relationship of the signature with clinical features. We found that the high-risk score was significantly correlated with tumor grade, AJCC stage, N stage, T stage, and survival status (Figure 5A).

To investigate the relationship between the prognostic signature and immune cell infiltration. Pearson correlation analysis showed that the signature score was significantly correlated with the infiltration of activated dendritic cells (cor = 0.152, p = 0.043), plasma cells (cor = −0.155, p = 0.040), CD8 T cells (cor = −0.193, p = 0.010), M1 macrophages (cor = 0.200, p = 0.008), and neutrophils (cor = 0.152, p = 0.043) (Figure 5B).

The tumor mutation burden (TMB) has been shown to be related to the clinical efficacy of immunotherapy (Samstein et al., 2019). To explore the value of our signature for predicting the efficacy of immunotherapy in PC, we assessed the TMB of PC patients in the high-risk and low-risk groups. We found that TMB of PC patients in the high-risk group was higher than that in the low-risk group, which implied that immunotherapy may be a potentially effective treatment to those PC patients with high-risk scores (Figure 5C).

Furthermore, we analyzed the correlation of the signature score with the therapy-related targets. Pearson’s correlation analysis showed that the risk score was significantly associated with PD-L1 (cor = 0.151, p = 0.044), VEGFR3 (cor = −0.194, p = 0.010), EGFR (cor = 0.177, p = 0.019), FLT3 (cor = −0.165, p = 0.028), KIT (cor = −0.164, p = 0.029), and MET (cor = 0.358, p = 1.026e–06) (Figure 5D).

LncRNAs could interact with miRNAs to modulate mRNA expression, thereby modulating the biological characteristics of malignant tumors. To explore the regulation of these selected LncRNAs, we constructed a regulatory network consisting of six lncRNAs, 107 miRNAs, and 209 mRNAs (Figure 6).

To better understand the function of the regulatory network, the “clusterProfiler” R package was employed to conduct a KEGG pathway and GO enrichment analysis. As shown in Figure 7, these genes in the regulatory network are enriched in many cellular components (CC) and molecular functions (MF). The most significantly enriched molecular functions included “protein serine/threonine kinase activity,” “Rab GTPase binding,” and “protein serine/threonine/tyrosine kinase activity,” In terms of KEGG pathway, the main significant pathways included “autophagy,” “small cell lung cancer,” “Toll-like receptor signaling pathway,” “pancreatic cancer,” “ErbB signaling pathway,” “colorectal cancer,” “endocrine resistance,” and “hedgehog signaling pathway.” These results indicated that the regulatory network may contribute to therapeutic resistance of PC through multiple signal pathways.

To screen small molecule drugs, 209 selected mRNA were further analyzed in the Connectivity Map (CMap). The top six most significant potential small molecule drugs were listed in Figure 8A, including vorinostat (C₁₄H₂₀N₂O₃), trichostatin A (C₁₇H₂₂N₂O₃), sirolimus (C₅₁H₇₉NO₁₃), phthalylsulfathiazole (C₁₇H₁₃N₃O₅S₂), GW-8510 (C₂₁H₁₅N₅O₃S₂), and daunorubicin (C₂₇H₂₉NO₁₀). And the 2D chemical structures of these potential agents were shown in Figure 8B.

Firstly, we evaluated the expression profiles and prognostic performance of these six lncRNAs. The expression level of AC245041.2, LINC01559, and AC005332.6 was significantly upregulated in PC than in normal tissues (Figure 9A). Moreover, the Kaplan-Meier survival analysis demonstrated the prognostic power of these six lncRNAs (Figure 9B). Then, qRT-PCR was applied to value the expression level of these lncRNAs in serum. Notably, only the expression of LINC01559 was markedly increased in the serum of PC patients, indicating that LINC01559 could serve as a diagnostic biomarker (Figure 9C).

We selected LINC01559 for further analysis. As shown in Figure 10A, we successfully silenced the expression level of LINC01559 in PC cells by si-LINC01559 transfection. Next, we explored the effect of silencing LINC01559 on the PC cells proliferation, migration, and invasion. CCK8 assay showed that the inhibition of LINC01559 led to a reduced viability in the PANC-1 and SW1990 cells (Figure 10B). Also, transwell assay was performed to demonstrate that the invasion and migration ability of PC cells were suppressed under LINC01559 depletion (Figures 10C,D). Furthermore, it was proved by wound healing assays that silencing LINC01559 obviously hindered the migration ability of PANC-1 and SW1990 cells (Figure 10E). These results suggested that knockdown of LINC01559 suppressed PC cell proliferation, migration, and invasion in vitro.

The role of LINC01559 in autophagy and chemotherapeutic resistance was further explored. WB analysis were employed to show that PC cells transfected with si-LINC01559 exhibited decreased the expression of LC3I/LC3II but increased the p62 expression, indicating that autophagy was inhibited after LINC01559 depletion (Figure 11A). However, gemcitabine (10 μM) treatment induced autophagy in PC cells (Figure 11A). This observation was further confirmed by the tandem LC3B-RFP-GFP fluorescence microscopy assay. As shown in Figure 11B, gemcitabine increased the number of red-only LC3 puncta in PC cells, implying an increase of autophagic flux. Inhibition of LINC01559 reduced the number of red-only LC3 puncta in GFP-mRFP-LC3-transfected PC cells compared with the cells treated with gemcitabine. Besides cell viability assay, colony formation assay, cell apoptosis assay, and WB analysis of apoptotic markers were performed. CCK-8 results showed that the IC50 value for gemcitabine was significantly increased in LINC01559-silenced PC cells (Figure 11C). In contrast, knockdown of LINC01559 significantly induced the colony-forming capacity of PANC-1 and SW1990 cells (Figure 11D) and increased the gemcitabine-induced apoptosis rates (Figure 11E). And the protein level of cleaved caspase3 and PARP were increased in LINC01559-downregulated cells with or without gemcitabine (10 μM) treatment (Figure 11F). These results suggested that the inhibition of LINC01559 could suppress autophagy and stimulate apoptosis, which would ultimately lead to sensitize PC cells to gemcitabine.