This article is mentioned in:

Introduction

Pancreatic adenocarcinoma (PAAD) is one of the most lethal malignancies worldwide, with an increasing incidence, a lack of specific symptoms, high invasiveness and a 5-year survival rate of ~10% (1,2). Approximately 60,430 new PAAD cases and 48,220 deaths of PAAD occurred in the USA in 2021 (2). PAAD has become one of the most difficult cancer types to treat and is projected to become the second leading cause of cancer-related death by 2030 worldwide (3). The efficacy of clinical PAAD treatments, including surgery, chemotherapy, radiotherapy and immunotherapy have improved; however, the 5-year overall survival (OS) rate remains poor (4,5). Therefore, it is necessary to identify novel prognostic biomarkers and therapeutic targets for PAAD to guide clinical practice.

Ferroptosis is an iron-dependent form of regulated cell death that is characterized by overwhelming lipid peroxide accumulation and redox imbalance (6,7). Recent studies have reported that ferroptosis is correlated with the progression and therapeutic response of cancers, and that activation of the ferroptosis process is a novel strategy for cancer treatment, especially for malignances resistant to traditional therapies (8,9). It has been reported that inhibition of ferroptosis confers sorafenib resistance and is associated with reduced OS in hepatocellular carcinoma (10,11). Moreover, ferroptosis may serve an important role in the development of PAAD. Using pancreatic glutathione peroxidase 4 (GPX4) conditional-knockout mice or high-iron diet models, Dai et al (12) reported that ferroptotic damage promotes KRAS-driven PAAD through the activation of the TMEM173/STING-dependent DNA sensor pathway in mice. Badgley et al (13) reported that the import of oxidized cysteine is a critical dependency of PAAD and that induction of ferroptosis by cysteine depletion may be a promising target for PAAD treatment. Another study reported that inhibition of aspartate aminotransaminase facilitates pancreatic cancer cell death by inducing ferroptosis (14). Therefore, it was hypothesized that novel ferroptosis-related biomarkers could have great potential for predicting the prognosis of patients with PAAD and identifying therapeutic targets for PAAD.

Long non-coding RNAs (lncRNAs) are RNA transcripts >200 nucleotides in length that lack a protein-coding capacity. As a newly reported type of gene regulator, lncRNAs regulate various important biological processes and participate in the development and progression of numerous diseases, including PAAD (15–17). For example, the lncRNA ENST00000480739 inhibits tumor cell invasion by the regulation of OS-9 and hypoxia inducible factor-1α in PAAD (18). In another study, upregulation of LINC01232 was reported to be significantly related to the poor prognosis of patients with PAAD, and LINC01232 exerted oncogenic roles in PAAD through the regulation of transmembrane 9 superfamily member 2 (19). In addition, lncRNA glutaminase-antisense (GLS-AS) is downregulated in PAAD, and loss of GLS-AS could impair GLS-mediated metabolism and suppress cancer progression (20). Several other lncRNAs, including PVT1 (21), AGAP2-AS1 (22) and PLACT1 (23), have been reported as potential therapeutic targets or prognostic biomarkers of PAAD. Therefore, lncRNAs are implicated in the development and progression of PAAD, and are potential targets for PAAD treatments and biomarkers for prognosis.

Ferroptosis is an important regulated form of cell death for cancers that is associated with lncRNAs (24). Numerous lncRNAs participate in the regulation of ferroptosis and serve important roles in tumorigenesis and progression through modulation of the transcription and post-transcriptional modification of ferroptosis-related genes (25,26). A recent study reported that the lncRNA NEAT1 improved ferroptosis sensitivity in non-small cell lung cancer by regulating the expression of acyl-CoA synthetase long-chain family member 4 (27). Mao et al (28) reported that the lncRNA P53RRA is downregulated in cancers and promotes ferroptosis through the nuclear sequestration of p53. Furthermore, it was reported that LINC00336 regulated the expression of cystathionine-β-synthase to induce ferroptosis in lung cancer by functioning as a competing endogenous sponge of microRNA-6852 (29). As these ferroptosis-related lncRNAs serve important roles in the pathophysiology of malignant diseases, they may be useful for prognosis prediction and may serve as therapeutic targets for patients with PAAD; however, the general role of these lncRNAs in the prognosis prediction remains unknown. Ferroptosis is rapidly gaining attention in cancer treatment and cancer genomics databases have been implemented to identify promising prognostic factors and to explore molecular mechanisms for multiple types of cancers in recent years. The present study used lncRNA expression profiles to identify ferroptosis-related lncRNA signatures (FRLSs), which may serve as novel biomarkers to predict the prognosis of patients with PAAD and may be utilized as potential therapeutic targets based on ferroptosis.

Materials and methods

Cell lines and tissue samples

The human pancreatic ductal epithelial cell line (GENNIO BIO, Guangzhou, China; Cat.no. JNO-089) and human pancreatic cancer cell lines AsPC-1, BxPC-3, PANC-1 and SW1990 (ATCC, USA) were cultured in complete growth medium with 10% fetal bovine serum (Gibco, USA), as recommended by the manufacturer. The human pancreatic ductal epithelial cell line and pancreatic cancer cell lines at purchase were first passage. Cells were incubated at 37°C in a humidified incubator with 5% CO2. LINC01133-specific small interfering (si)RNAs (LINC01133-siRNA-1, sense 5′-AAUGGAUCCAUUCCCUGCAACUG-3′; and LINC01133-siRNA-2, sense 5′-GAAGUGGAAGCAAAGUUCUCCAAAG-3′) and a negative control siRNA (5′-TTCTCCGAACGTGTCACGT-3′) were purchased from Shanghai GeneChem Co., Ltd. and transiently transfected using Lipofectamine® 3000 (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol.

A total of 30 paired freshly frozen PAAD tissues and adjacent nontumor (ANT) tissues were obtained from the Pancreatic Tumor Bank, Institute of Hepatopancreatobiliary Surgery, Southwest Hospital, Third Military Medical University (Chongqing, China) between May and November 2020. The clinicopathological characteristics of the patients are presented in Table SI. The patient inclusion criteria were as follows: i) Pathologically diagnosed with pancreatic adenocarcinoma; ii) no preoperative anti-tumor treatment; and iii) had complete clinicopathological data. Exclusion criteria were as follows: i) Previous chemotherapy or other anti-tumor treatment before surgery; and ii) lacked complete clinical information. The diagnoses were independently reviewed by two senior pathologists and classified according to World Health Organization criteria. PAAD and ANT tissues were stored for total RNA isolation (snap-frozen in liquid nitrogen and then stored at −80°C). Written informed consent was obtained from all patients prior to the study. The present study was approved by the ethical committee of Southwest Hospital (approval no. 20210312).

Cell Counting Kit-8 (CCK-8) assays

AsPC-1 and BxPC-3 cells transfected with si-control and si-LINC01133 were treated with 10 µM erastin for 24 h. Then, cell viability was assessed using a Cell Counting Kit-8 assay (Dojindo, Shanghai China Co., Ltd.) according to the manufacturer's protocols. The data were assessed in at least three independent experiments. Erastin, a ferroptosis-inducing agent, was purchased from APExBIO Technology LLC.

Reverse transcription-quantitative PCR (RT-qPCR)

Total RNA was extracted from the tumor cells or tissues using an Ultrapure RNA kit (CoWin Biosciences). The RNA concentration was assessed using a NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific, Inc.). RT-qPCR was performed using a PrimeScript RT reagent kit and SYBR® Premix Ex Taq™ kit (Takara Biotechnology Co., Ltd.) according to the manufacturers' protocol. The following thermocycling conditions were used for qPCR: 95°C for 30 sec, followed by 40 cycles at 95°C for 5 sec and 60°C for 30 sec. Relative mRNA expression levels were calculated using the 2−ΔΔCq method (30) and normalized to GAPDH mRNA expression levels. The primer sequences used were as follows: LINC01133 forward (F), 5′-TGGATCCATTCCCTGCAACT-3′ and reverse (R), 5′-GTGCTGGGCTCTGGATTCTT-3′; CASC8 F, 5′-GCAGTGAGCCAAGGAGCAAT-3′ and R, 5′-AACCGCAACACTGGTTGTGT-3′; SLC7A11 F, 5′-TCCCTCTATTCGGACCCATTTA-3′ and R, 5′-TTCTTCTGGTACAACTTCCAGT-3′; GPX4 F, 5′-CGATACGCTGAGTGTGGTTTGC-3′ and R, 5′-CATTTCCCAGGATGCCCTTG-3′; and GAPDH F, 5′-CAGGAGGCATTGCTGATGAT-3′ and R, 5′-GAAGGCTGGGGCTCATTT-3′.

Collection of data on patients with PAAD

The Cancer Genome Atlas (TCGA)-PAAD dataset and clinical information were identified and downloaded from the TCGA database (portal.gdc.cancer.gov, accessed on 20 September 2021) using the key words ‘pancreas’, ‘transcriptome profiling’, ‘HTSeq-FPKM’ and ‘TCGA-PAAD’. Data from 181 patients with PAAD were used as a training cohort. Data on another 90 PAAD samples from the Pancreatic Cancer-AU (PACA-AU) dataset were downloaded from the International Cancer Gene Consortium (ICGC) database (https://dcc.icgc.org, accessed on 24 September 2021) and were used as the validation cohort. The inclusion criteria were samples with complete follow-up information (survival time >30 days) and complete clinicopathological data. The exclusion criterion were samples without survival information.

Gene Expression Profiling Interactive Analysis (GEPIA)

GEPIA (http://gepia.cancer-pku.cn/index.html) is a web server for large-scale expression analysis and interactive analysis (31). This database is a highly cited resource for the comparison of the expression levels of signature lncRNAs (32–35).

Construction and evaluation of the prognostic FRLS

The list of ferroptosis-related genes was obtained from the FerrDb Database (zhounan.org/ferrdb/, accessed on 20 September 2021). The expression matrixes of lncRNAs and mRNAs were identified using classification and annotation analyses. Potential ferroptosis-related lncRNAs were identified using Pearson's correlation analysis (r>0.4; P<0.001). The limma R package was used to screen the differentially expressed ferroptosis-related lncRNAs between normal and tumor tissues with a false discovery rate (FDR) <0.05 and |log2 fold change (FC)|>1. Subsequently, univariate Cox regression analysis was performed to screen prognostic ferroptosis-related lncRNAs (threshold, P<0.01). Next, iterative least absolute shrinkage and selection operator (LASSO) Cox regression was used with the glmnet R package to identify the optimal prognostic signature. Consensus genes were assessed with a frequency >100 within 1,000 LASSO Cox regressions. Consensus genes were then incorporated into the multivariate Cox regression analysis to construct the risk score. Using the median risk score as the cutoff point, patients with PAAD were divided into a high-risk group (over the median risk score) and a low-risk group (no more than the median risk score). Kaplan-Meier survival analysis was performed to evaluate the OS rates between the two groups with the survival and survminer R packages. Where two survival curves intersected, the two-stage procedure (36) was used, which had good applicability and robustness in different crossing situations. Receiver operating characteristic (ROC) curve analysis was performed using the survivalROC R package; P<0.05 was used as the cutoff criterion. Furthermore, principal component analysis (PCA) was performed using the Rtsne and ggplot2 packages to assess the classification ability of the risk signature. Then, the rms package was used to generate the nomogram for the prediction of the 1, 3 and 5-year OS rates of patients with PAAD. A calibration curve was constructed to assess the consistency between the actual OS rates and the nomogram-predicted OS rates. Finally, decision curve analysis (DCA) was performed to evaluate the net benefits with the different predictors using the rmda package.

Gene set enrichment analysis (GSEA)

GSEA (http://www.gsea–-msigdb.org/gsea/downloads.jsp, version 4.1.0) was performed to evaluate the high and low-risk groups of the prognostic FRLS and explore the possible cellular pathways. Gene Ontology (GO) term enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were annotated using the c5.go.v7.4.symbols.gmt and c2.cp.kegg.v7.4.symbols.gmt gene sets from the GSEA Molecular Signatures Database. Random assortment times were set to 1,000. Biological processes GO terms and KEGG pathways were selected on the sorted samples using default weighted enrichment statistics, and their enrichment scores and P-values were ranked. Gene clusters with adjusted P<0.05 were considered as significantly enriched genes.

Assessment of the correlation of the risk score with immune infiltration

Tumor IMmune Estimation Resource (TIMER) (37), CIBERSORT (38), CIBERSORT-ABS (38), QUANTISEQ (39), Microenvironment Cell Populations-counter (MCP-counter) (40), XCELL and Estimating the Proportion of Immune and Cancer cells (EPIC) (41) algorithms were used to evaluate and compare the abundance of immune cells between the high- and low-risk groups based on the FRLS. The immune cell infiltration and functions were evaluated between the high and the low-risk groups using the two-sample Wilcoxon test. P<0.05 was used as the cutoff criterion. Single-sample GSEA (ssGSEA) was used to assess the immune functions and cells between the high and low-risk groups using the GSVA package. The correlation of the risk score and immune checkpoints was evaluated using the difference in gene expression levels between the two groups.

Statistical analysis

All computational and statistical analyses were performed using R software (https://www.r–project.org/, version 4.1.0). Fisher's exact test or the χ2 test were used for analyses of clinical data. Unpaired or paired two-tailed Student's t-test was used to compare two groups; one-way ANOVA followed by Tukey's post hoc test was used for multiple group comparison. The correlation between ferroptosis-related genes and lncRNAs was assessed using Pearson's correlation analysis. The ferroptosis-related lncRNA and gene co-expression network were constructed and visualized using Cytoscape software (https://cytoscape.org/, version 3.8.2), and the ggalluvial package was used to generate the Sankey diagram. Univariate and multivariate Cox regression analyses were performed to assess whether the prognostic model was independent of other traditional clinicopathological factors (including age, Sex, tumor grade and TNM stage) in the prediction of OS of patients with PAAD. The hazard ratio (HR) and 95% confidence interval (CI) were evaluated using the survival R package. The Cox model was presented as the hazard function denoted by h(t). Briefly, the hazard function could be interpreted as the risk of dying at time t (42). The Cox regression was fitted using the covariates as follows: Age, Sex, grade, stage and risk score. Univariate Cox analyses for all these variables were first calculated; then multivariate cox analyses were fitted to describe how these factors jointly impacted on survival. For each analysis, statistical significance was set at P<0.05.

Results

Identification of prognostic ferroptosis-related lncRNAs in patients with PAAD

A total of 253 ferroptosis-related genes were identified using the FerrDb database (43), and expression data were available for 246 of these genes in the TCGA-PAAD dataset (Table SII). A total of 2,749 ferroptosis-related lncRNAs were identified using Pearson's correlation analysis (r>0.4; P<0.001). Subsequently, 215 differentially expressed lncRNAs were identified (log2FC>1; P<0.05; Table SIII). Furthermore, univariate Cox proportional hazards analysis was performed and 26 prognostic ferroptosis-related lncRNAs were identified (P<0.01; Fig. 1). Among them, 14 lncRNAs were identified as poor prognostic factors (HR >1), and 12 were identified as favorable prognostic factors (HR <1).

Figure 1. Identification of prognostic ferroptosis-related lncRNAs in pancreatic adenocarcinoma. Forest plot for the P-value and hazard ratio of selected lncRNAs evaluated using univariate Cox regression analysis (all P<0.01). Twelve lncRNAs were significantly correlated with a good prognosis (blue) and 14 lncRNAs were significantly correlated with a poor prognosis (red). lncRNA, long non-coding RNA.

Construction and validation of the prognostic FRLS

To construct a risk signature, the aforementioned 26 ferroptosis-related lncRNAs were screened using the LASSO regression algorithm in the TCGA database (Fig. 2A and B); 17 lncRNAs were obtained(LINC01133, CASC8, AL356740.3, LINC02535, LINC01091, HNF1A.AS1, TRPC7.AS1, UNC5B.AS1, AC015660.1, LINC01843, AC130456.2, AC068580.2, LINC02004, LINP1, AC092171.3, DCST1.AS1 and LINC02178). Of these, nine lncRNAs were further identified by multivariate Cox proportional hazards regression analysis (Table I), including LINC01133, CASC8, AL356740.3, LINC02535, LINC01091, AC068580.2, LINC02004, AC092171.3 and AC015660.1, which were used to construct a prognostic FRLS. Correlations between the expression levels of these nine lncRNAs and ferroptosis-related genes are presented (Fig. 2C; Table SIV). Among these candidates, LINC01091, AC092171.3 and AL356740.3 were predicted to be protective lncRNAs, whereas the remaining lncRNAs were predicted to be risk lncRNAs (Fig. 2D). The risk score of each sample was calculated by taking the sum of the expression level of each lncRNA and multiplying it by the corresponding coefficients for each sample.

Figure 2. Construction of the ferroptosis-related lncRNA signature in The Cancer Genome Atlas cohort. (A) Profiles of LASSO coefficients. (B) LASSO penalized Cox regression analysis. The vertical dashed line is at the optimal log(λ) value. (C) Construction of ferroptosis-related lncRNA and mRNA regulatory networks using Cytoscape. Red hexagons represent prognostic lncRNAs and green ellipses represent ferroptosis-related mRNAs. (D) Sankey diagram presenting the relationships among ferroptosis-related genes, lncRNAs and risk type. LASSO, least absolute shrinkage and selection operator; lncRNA, long non-coding RNA.

I. Multivariate Cox regression analysis for the risk model based on nine ferroptosis-related lncRNAs.

I.

Multivariate Cox regression analysis for the risk model based on nine ferroptosis-related lncRNAs.

lncRNA	Coefficient	HR	HR.95L	HR.95H	P-value
LINC01133	0.010370471	1.01042443	1.001045639	1.019891091	0.029285882
CASC8	0.108503871	1.114609223	0.982588369	1.264368436	0.091625469
AL356740.3	-0.454385767	0.634837786	0.364897336	1.104472342	0.107778175
LINC02535	0.620347204	1.859573581	1.247702734	2.771504628	0.002311884
LINC01091	-0.926940642	0.39576264	0.116636711	1.342871085	0.137008934
AC068580.2	0.121716427	1.129433779	0.984743052	1.295384273	0.081831429
LINC02004	0.885501211	2.424199123	1.138107086	5.163610229	0.021716368
AC092171.3	-0.666412005	0.513547883	0.316408144	0.833516562	0.006998521
AC015660.1	0.176288966	1.192782683	1.003061085	1.418388721	0.046092594

[i] lncRNA, long non-coding RNA; HR, hazard ratio; 95L and 95H, 95% confidence interval.

Patients with PAAD were subsequently stratified into high and low-risk groups using the median risk score as the cutoff value. The classification ability of the prognostic FRLS was assessed using PCA in both the TCGA training and ICGC validation cohorts (Fig. 3A and B, respectively). The prognostic effectiveness of this model for OS status of patients with PAAD was assessed using Kaplan-Meier curve analysis. The results from the TCGA cohort demonstrated that the OS rate of patients in the high-risk group was significantly lower (P<0.001; Fig. 3C), which suggested that the newly developed risk signature may effectively predict survival. Furthermore, a risk score distribution dot plot demonstrated that more patients survived in the low-risk group compared with the high-risk group (Fig. 3E and G). The expression levels of the nine ferroptosis-related lncRNAs in the risk signature is presented as a heatmap in Fig. 3I. Furthermore, the predictive accuracy of the model was evaluated using time-dependent ROC analysis, with area under the curve (AUC) values for 1, 2 and 3-year OS of 0.748, 0.748 and 0.785, respectively (Fig. 3K).

Figure 3. Training and validation of the prognostic ferroptosis-related lncRNA signature in the TCGA cohort and ICGC cohort. Principal component analysis between the high and low-risk groups in (A) the training (TCGA) and (B) validation (ICGC) cohorts. Kaplan-Meier curves of overall survival in the (C) TCGA cohort (log-rank test) and (D) ICGC cohort (two-stage test). Association between the prognostic signature and survival outcome of patients with pancreatic adenocarcinoma, distribution of the (E and F) risk scores, (G and H) survival time and (I and J) ferroptosis-related lncRNA levels. Receiver operating characteristic curve for the survival prediction model in the (K) TCGA and (L) ICGC cohorts. AUC, area under the curve; ICGC, International Cancer Gene Consortium; lncRNA, long non-coding RNA; TCGA, The Cancer Genome Atlas.

To assess whether the prognostic significance of FRLS remained in other populations, the same analyses were performed in the ICGC cohort. In accordance with the results in the TCGA training cohort, patients in the high-risk group in the ICGC validation cohort demonstrated a significantly shorter OS compared with patients in the low-risk group (Fig. 3D). LncRNAs risk distribution, survival rate and expression profiles in the validation cohort were shown in Fig. 3F, H and J. The AUCs at 1, 2 and 3 years OS reached 0.733, 0.711 and 0.705, respectively (Fig. 3L). All these results suggested that the prognostic FRLS may accurately and stably predict the clinical outcome of patients with PAAD.

Correlation of the prognostic FRLS with clinicopathological factors

The independent value of the prognostic FRLS was assessed using univariate and multivariate Cox regression analyses. Univariate Cox regression analysis demonstrated that the risk score, age and tumor grade were associated with the OS of patients with PAAD (P<0.001; Fig. 4A). Multivariate Cox regression analysis demonstrated that the risk score was an independent prognostic factor for the prediction of the OS of patients with PAAD (P<0.001; Fig. 4B).

Figure 4. Association between the risk score and different clinicopathological factors. Forest plots for (A) univariate and (B) multivariate Cox regression analyses of prognostic value in The Cancer Genome Atlas cohort. The hazard ratios of stratified clinical factors and risk score were assessed using the Cox model, including age (>60 vs. ≤60), sex (male vs. female), tumor grade (G3 + G4 vs. G1 + G2), tumor stage (S3 + S4 vs. S1 + S2) and risk score (high vs. low). (C) The heatmap demonstrated the distribution of the nine ferroptosis-related long non-coding RNAs and clinicopathological factors in the high-risk and low-risk patients with pancreatic carcinoma. There were significant differences in tumor grade between high and low-risk groups (*P<0.05). Kaplan-Meier curves of the overall survival compared between the high and low-risk groups stratified by clinical variables: (D and E) Age (≤65 and >65 years); (F and G) sex (female and male); (H and I) tumor grade (G1-2 and G3-4); and (J and K) N stage (N0 and N1-3). T, tumor; N, node; M, metastasis.

To further evaluate the roles of the FRLS in the development of PAAD, the correlations of the risk score with clinicopathological features were evaluated. The associations between the prognostic model based on the nine-gene FRLS and the clinicopathological factors are presented as a heatmap (Fig. 4C). Stratification analysis was performed according to the following clinical factors: Age (≤65 and >65), sex (female and male), tumor grade (G1-2 and G3-4) and N stage (N0 and N1-3). The risk signature demonstrated significant differences in prognosis, with high-risk patients having significantly poorer OS compared with low-risk patients in all subgroups (Fig. 4D-K). In summary, these results suggested that the FRLS was closely related to PAAD progression.

Construction of a predictive nomogram

To assess the potential clinical practicality of the prognostic model based on the nine-gene FRLS, a nomogram was constructed using the risk score and clinicopathological features to estimate the 1, 3 and 5-year OS rates of patients with PAAD (Fig. 5A). The calibration curve of the nomogram demonstrated that the predicted OS rate was close to the actual OS rate at 1, 3 and 5 years (Fig. 5B). Furthermore, the DCA results demonstrated that this nomogram based on the nine-gene FRLS had better clinical practicality for the prognosis prediction of patients with PAAD (Fig. 5C). Finally, a 5-year OS time-dependent ROC curve was generated. The AUC value of the clinical prognostic nomogram was 0.748, which was markedly higher compared with the AUC values of age, sex, grade and stage (Fig. 5D), which further suggested the discriminative ability of the nine-gene FRLS combined with tumor pathological characteristics to predict the survival time of patients with PAAD.

Figure 5. Clinical prognostic nomogram for survival prediction. (A) Clinical prognostic nomogram developed to predict 1, 3 and 5-year OS. (B) Calibration curves presenting nomogram predictions for 1-year, 3-year and 5-year OS. (C) Decision curve analysis of the clinical practicality of the nomogram for prognosis prediction. (D) Time-dependent receiver operating characteristic curve analyses for predicting OS at 5 years by risk score, age, sex, grade and stage. AUC, area under the curve; G, grade; N, node; OS, overall survival; T, tumor.

Identification of FRLS-associated biological pathways

To evaluate the significant changes in functional phenotypes between the high- and low-risk groups of the prognostic FRLS, GSEA between the two risk groups was used. The GO terms ‘antigen processing and presentation of peptide antigen via MHC class 1’ [normalized enrichment score (NES)=1.73; P=0.008], ‘antigen processing and presentation of exogenous antigen via MHC class 1’ (NES=1.66; P=0.02), ‘interleukin-1-mediated signaling pathway’ (NES=1.59; P=0.002) and ‘innate immune response in mucosa’ (NES=1.62; P=0.002) were significantly enriched in patients with PAAD with high risk scores. However, ‘voltage-gated cation channel activity’ (NES=−1.95; P=0.000), ‘voltage-gated potassium channel activity’ (NES=−1.89; P=0.000) and ‘mast cell activation involved in the immune response’ (NES=−1.60; P=0.004) were significantly enriched in patients with PAAD with low risk scores (Fig. 6A-E; Table SV). Notably, several cell death-related biological processes, such as ‘positive regulation of cell aging’, ‘regulation of necrotic cell death’ and ‘cell redox homeostasis’, were enriched in PAAD samples with high risk scores (Table SV).

Figure 6. Functional analysis of the low-risk and high-risk groups. In the GO enrichment analyses, three GO items, namely, (A) antigen processing and presentation of peptide antigen via MHC class 1, (B) cadherin binding involved in cell adhesion and (C) innate immune response in mucosa, were significantly enriched in the high-risk group. GO items (D) voltage-gated cation channel activity and (E) voltage-gated ion channel activity showed significantly differential enrichment in the low-risk group. In the KEGG enrichment analyses, three KEGG items, namely, (F) p53 signaling pathway, (G) proteasome and (H) glycolysis/gluconeogenesis, were significantly enriched in the high-risk group. KEGG items (I) calcium signaling pathway and (J) glycine, serine and threonine metabolism showed significantly differential enrichment in the low-risk group. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MHC, major histocompatibility complex.

KEGG pathway analysis demonstrated that the lncRNAs in the newly developed signature were enriched in classical tumor pathways and metabolism (Fig. 6F-J; Table SVI). The ‘p53 signaling pathway’ (NES=1.72; P=0.009, ‘proteasome’ (NES=1.75; P=0.013) and ‘glycolysis/gluconeogenesis’ (NES=1.64; P=0.008) were significantly enriched in the high-risk group. The ‘calcium signaling pathway’ (NES=−1.68; P=0.002), ‘glycine, serine and threonine metabolism’ (NES=−1.66; P=0.010) and ‘tryptophan metabolism’ (NES=−1.56; P=0.037) were significantly enriched in the low-risk group. These results indicated that the gene sets were connected with the tumor immune microenvironment, as well as carcinogenesis and tumor progression pathways.

Correlation of the FRLS with immune infiltration

As the GSEA results suggested that the FRLS was associated with the tumor immune microenvironment, the correlation of the nine-gene FRLS with immune cell infiltration was assessed. Heatmaps of immune infiltration were generated using the TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCP-counter, XCELL and EPIC algorithms (Fig. 7A). Comparative analysis of immune functions and immune cells was performed to assess the differences in cytolytic activity, major histocompatibility complex (MHC) class I, T cell co-inhibition, T cell co-stimulation, type I IFN response, type II IFN response, B cells, CD8+ T cells, mast cells, neutrophils, plasmacytoid dendritic cells (pDCs), T helper (Th) cells, follicular helper T cells (Tfhs) and tumor-infiltrating lymphocytes (TILs) between the two risk groups (P<0.05; Fig. 7B and C). Furthermore, violin plots were generated based on the CIBERSORT algorithm and demonstrated the relationship between the risk score and immune cell infiltration (Fig. 7D). Patients with PAAD in the low-risk group demonstrated significantly higher ratios of CD8+ T cells (P<0.001), activated memory CD4+ T cells (P=0.003) and monocytes (P=0.01) compared with the high-risk group. However, M0 macrophages (P=0.002) were significantly positively associated with a high-risk score. In addition, correlation scatter plots demonstrated the same trends between the prognostic signature and immune cell infiltration in CD8+ T cells (r=−0.28; P=0.002), activated memory CD4+ T cells (r=−0.23; P=0.011), monocytes (r=−0.18; P=0.045) and M0 macrophages (r=0.25; P=0.005) (Fig. 7E-H).

Figure 7. Correlations of the risk score with the immune microenvironment and immune checkpoints compared between the high and low-risk groups. (A) Heatmaps for immune responses using the TIMER, CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCP-COUNTER, XCELL and EPIC algorithms in the two groups. P<0.05 was considered to indicate a statistically significant difference. Boxplots presenting the single-sample gene set enrichment analysis (B) immune function and (C) immune cell scores compared between the two groups. (D) Violin plot presenting the immune-infiltrating lymphocytes in the low- and high-risk groups. (E-H) Scatter plots presenting the correlation between the risk score and immune cell infiltration. (I) Expression of immune checkpoints compared between the low- and high-risk groups. The Kruskal-Wallis test was used to assess differences between the groups. *P<0.05, **P<0.01 and ***P<0.001. APC, antigen presenting cell; CCR, chemokine receptor; DCs, dendritic cells; HLA, human leukocyte antigen; iDCs, inflammatory dendritic cells; Tfh, T follicular helper cell; MHC, major histocompatibility complex; NK, natural killer; ns, not significant; pDCs, plasmacytoid dendritic cells; TIL, tumor-infiltrating lymphocytes; Treg, regulatory T cells.

Given the importance of immune checkpoint blockade immunotherapy, the relationship between the expression of immune checkpoints and the two risk groups was evaluated (Fig. 7I). The low-risk group tended to demonstrate higher mRNA expression levels of 22 immune checkpoint genes, including TIGIT, inducible T-cell co-stimulator, adenosine A2a, CD40 ligand, CD27, CTLA4, transmembrane and immunoglobulin domain containing 2, B and T-lymphocyte attenuator, lymphocyte activating 3, BTNL2, programmed cell death 1, leukocyte-associated immunoglobulin-like receptor 1, CD200 receptor 1, CD244, CD28, CD200, KIR3DL1, IDO2, CD48, CD160, CD44 and TNFRS, whereas the other three immune checkpoint genes (TNFSF9, HHLA2 and CD276) were highly expressed in high-risk patients (all P<0.05). These results suggested that the FRLS may be a potential biomarker for immune checkpoint inhibitor therapy.

Validation of ferroptosis-related lncRNA expression in cells and tissues

LINC01133 and CASC8 were selected for validation, as their LASSO coefficients were the top two among the risk lncRNAs. RT-qPCR was performed to evaluate the expression levels of LINC01133 and CASC8 in pancreatic cancer tissues, and the results demonstrated that LINC01133 and CASC8 were significantly increased compared with ANT tissues (Fig. 8A and B). Furthermore, the mRNA expression levels of SLC7A11 and GPX4, two essential ferroptosis-related genes (43), were significantly increased in pancreatic cancer tissues compared with ANT tissues (Fig. 8C and D). The SLC7A11/GPX4 axis serves as the canonical defense against ferroptosis by facilitating intracellular glutathione biosynthesis and alleviating lipid peroxidation (44). Next, the GEPIA (31) database was used to verify the expression levels of ferroptosis-related lncRNAs in pancreatic cancer tissues, ANT tissues and normal pancreatic tissues. The RNA expression levels of LINC01133 and CASC8 were significantly upregulated in pancreatic cancer tissues compared with normal tissues (Fig. 8E and F). Moreover, Kaplan-Meier survival analysis demonstrated that upregulation of these two lncRNAs were associated with shorter overall survival in the GEPIA database (Fig. 8G and H). The expression levels of LINC01133 and CASC8 were assessed in pancreatic cancer cells and pancreatic ductal epithelial cell lines. The results demonstrated that LINC01133 and CASC8 expression levels were significantly upregulated in pancreatic cancer cells compared with pancreatic ductal epithelial cells (Fig. 8I and J). To evaluate the role of LINC01133 in the ferroptosis process of pancreatic cancer, LINC01133 expression was then knocked down in AsPC-1 and BxPC-3 cells, both of which have high levels of LINC01133 (Fig. 8K). The results demonstrated that downregulation of LINC01133 significantly decreased the survival rate of AsPC-1 and BxPC-3 cells after ferroptosis induction using erastin (Fig. 8L), which suggested that LINC01133 may have been involved in the ferroptosis of pancreatic cancer cells.

Figure 8. Evaluation of ferroptosis-related lncRNA expression in cells and tissues. RT-qPCR analysis of (A) LINC01133, (B) CASC8, (C) SLC7A11 and (D) GPX4 expression levels in 30 paired pancreatic cancer tissues and ANT tissues. Comparison of the differential expression of (E) LINC01133 and (F) CASC8 in the Gene Expression Profiling Interactive Analysis database. Kaplan-Meier survival curves for the ferroptosis-related lncRNAs (G) LINC01133 and (H) CASC8. *P<0.05, **P<0.01 and ***P<0.001. (I and J) RT-qPCR analysis of LINC01133 and CASC8 expression levels in pancreatic cancer and HPDE normal pancreatic ductal epithelial cell lines. *P<0.05, **P<0.01 and ***P<0.001 vs HPDE6-C7. (K) AsPC-1 and BxPC-3 cells were transiently transfected with si-control and si-LINC01133 using Lipofectamine. LINC01133 expression was confirmed to be knocked down 24 h after transfection using RT-qPCR. ***P<0.001 vs si-control. (L) AsPC-1 and BxPC-3 cells transfected with si-control and si-LINC01133 were treated with 10 µM erastin for 24 h. Then, cell viability was assessed using a Cell Counting Kit-8 assay. Unpaired or paired two-tailed Student's t-test was used to compare two groups; one-way ANOVA followed by Tukey's post hoc test was used for multiple group comparison. lncRNA, long non-coding RNA; RT-qPCR, reverse-transcription quantitative PCR; T, pancreatic cancer tissues; ANT, adjacent non-tumor tissues; PAAD, pancreatic adenocarcinoma; si, small interfering RNAs; HR, hazard ratio.

Discussion

LncRNAs serve an important role in the ferroptotic process in cancers (24). In the present study, the relationship between lncRNAs and ferroptosis-related genes in PAAD, and the expression of ferroptosis-related lncRNAs and their relationship with the prognosis of patients with PAAD were assessed. A novel prognostic FRLS containing nine ferroptosis-related lncRNAs was then constructed to predict the prognosis of patients with PAAD. In both the training cohort from the TCGA database and the validation cohort from the ICGC database, the FRLS demonstrated robust prognostic predictive ability for patients with PAAD. Univariate and multivariate Cox regression analyses demonstrated that the FRLS was an independent prognostic factor for PAAD. Furthermore, combining the prognostic FRLS with clinicopathological factors, a nomogram with improved predictive ability for OS was constructed. Calibration curves demonstrated that the nomogram predictions for 1-year, 3-year and 5-year survival were close to the actual OS rates. Functional enrichment analyses demonstrated the difference in immune-related biological processes and classical tumor pathways between the high and low-risk groups. The differential immune infiltration between the two risk groups was further analyzed by comparison of the abundance of immune cells and immune functions and the mRNA expression levels of immune checkpoint genes. These results demonstrated that the nine-gene FRLS was a novel prognostic biomarker and potential therapeutic target for PAAD.

The prognostic signature constructed in the present study was derived from nine ferroptosis-related lncRNAs. Among them, LINC01133 had been previously reported to serve important roles in gastric cancer (45), ovarian cancer (46), breast cancer (47) and cervical cancer (48,49). In PAAD, LINC01133 was reported to be upregulated, and loss of LINC01133 was reported to suppress the development of PAAD cells via the Wnt signaling pathway (50). In another study, Huang et al (51) reported that overexpression of LINC01133 was associated with poorer prognosis in patients with PAAD and that the C/EBPβ-LINC01133 axis served an oncogenic function through upregulation of CCNG1. CASC8 was reported to be correlated with the progression of non-small cell lung cancer (52), retinoblastoma (53) and several other cancers (54,55). Upregulation of CASC8 was reported to be associated with the poor prognosis of patients with PAAD (56). In the prognostic signature identified in the present study, LINC01133 and CASC8 were upregulated and associated with poor prognosis in PAAD, which was consistent with the findings of the aforementioned studies. Furthermore, LINC02535 was demonstrated to regulate DNA damage repair in cervical cancer (57) and to promote cell growth in poorly differentiated gastric cancer (58). However, the role and prognostic value of the remaining ferroptosis-related lncRNAs in cancer were unknown. The specific mechanisms of ferroptosis-related lncRNAs and their contributions to ferroptosis in PAAD are complex and still obscure, and require further study.

The results of KEGG enrichment analyses demonstrated that several pathways were enriched when compared between the high and low-risk groups, including the p53, proteasome, glycolysis, calcium and chemokine signaling pathways. Previous studies reported that p53 signaling pathway was closely related to ferroptosis (59,60). Ferroptosis can be regulated by p53 signaling pathway through transcriptional or post-translational mechanisms. p53 can facilitate ferroptosis through the inhibition of SLC7A11 expression or by the induction of spermidine/spermine N (1)-acetyltransferase 1 and glutaminase 2 activity. However, p53 can also inhibit ferroptosis by the inhibition of dipeptidyl peptidase 4 or through the enhancement of cyclin dependent kinase inhibitor 1 A (61). It has been reported that proteasome function is often inhibited during ferroptosis (62). Numerous metabolic and degradation pathways, including the ubiquitin-proteasome system, orchestrate the complex ferroptotic response through the regulation of iron accumulation or lipid peroxidation (63). Recent studies reported that cellular energy metabolism activities such as glycolysis, pentose phosphate pathway and the tricarboxylic acid cycle are associated with ferroptosis by the regulation of essential ferroptosis molecules, including nicotinamide adenine dinucleotide phosphate, reactive oxygen species and glutathione (64). It was also reported that calcium-iron crosstalk regulates iron-induced ferroptosis (65). Martin-Sanchez et al (66) reported that ferrostatin-1, an inhibitor of ferroptosis, could prevent the upregulation of IL-33 and other chemokines and cytokines. In general, KEGG analyses suggested the probable mechanisms of the high-risk and low-risk groups, and these signaling pathways were generally related to ferroptosis. However, the specific mechanisms of ferroptosis-related lncRNAs in PAAD are complex and still unclear and require further evaluation in future work.

Immune cells that infiltrate tumors exert diverse effects on PAAD progression (67). The number and proportion of infiltrating immune cells serve important roles in regulation of the immunotherapy response and affect cancer progression (68). However, the relationship between the ferroptotic process and immune infiltration levels in PAAD remains unclear. In the present study, the numbers of B cells, CD8+ T cells, mast cells, neutrophils, pDCs, Th cells, Tfhs and TILs were significantly associated with risk score, which suggested that these immune cells may serve a substantial role in the development of PAAD. Compared with patients with PAAD in the low-risk group, those in the high-risk group tended to have fewer infiltrated tumor-killing immune cells, such as B cells, CD8+ T cells, Th cells and TILs, which suggested that the antitumor ability of the high-risk group was weaker. Therefore, we hypothesized that ferroptosis is markedly associated with the number and proportion of tumor-infiltrating immune cells in PAAD. Recently, certain studies have begun to report the effect of combining cancer immunotherapy with ferroptosis inducers. Wang et al (69) first reported that immunotherapy-activated CD8+ T cells could promote tumor cell lipid peroxidation and ferroptosis through the release of IFNγ to decrease the expression of SLC3A2 and SLC7A11. Given the potential role of the ferroptosis-related lncRNAs in the signature identified in the present study in the process of immune infiltration, the underlying mechanism of these lncRNAs in the regulation of the immune response of PAAD should be evaluated in further studies. Checkpoint blockade immunotherapy has been reported to have improved the prognosis of patients multiple advanced malignancies; however, not all patients respond to this treatment (70). The results of the present study demonstrated that patients with PAAD in the high-risk group, based on the FRLS, tended to have lower mRNA expression levels of immune checkpoints. Significant differences in these molecules between the two groups indicated differences in sensitivity to immunotherapies. However, further studies are needed to confirm the feasibility and effectiveness of this combined regimen of ferroptosis inducers and immune checkpoint inhibitors for the treatment of PAAD.

There have been numerous systematic bioinformatic analyses of PAAD; however, the present study proposed a novel nine-gene FRLS to predict prognosis for patients with PAAD. Furthermore, a prognostic nomogram including the FRLS was constructed that demonstrated improved accuracy for survival risk stratification. Importantly, the FRLS indicated in the present study may be involved in immune infiltration and have an effect on cancer immunotherapy. In addition, the AUC values for the training and test groups were 0.785 and 0.733, respectively. In previous studies, several prognostic models have also been constructed based on lncRNAs to assess the prognosis of patients with PAAD. Chen et al (71) constructed an evaluation model based on ferroptosis-related lncRNAs, and the AUC values were 0.70 at 1 year, 0.71 at 3 years and 0.75 at 5 years, which were lower than those for the model developed in the present study. In another study, Shi et al (72) reported a prognostic three lncRNA signature, and the AUC of this model was 0.716. Wei et al (73) built a nine lncRNA signature model for PAAD, and the AUC for this model for the prediction of the 2-year survival rate was 0.703. These results suggested that the prognostic model based on nine ferroptosis-related lncRNAs constructed in our study may improve the prognostic power for PAAD. The AUC values at 2 and 3 years slightly decreased compared with 1 year. Several other prognostic models also reported a decrease in the AUC values at 2 and 3 years (71,74,75). Numerous factors affect the AUC value and one of the important factors is the uncertainty of the sample. The so-called uncertainty of the sample refers to whether the corresponding label is uncertain for the same sample; that is, the samples with the same feature values. In the present study, the ICGC cohort had only six patients with a survival time of >3 years; however, five of them had a survival status of death, which may have affected the AUC value for 3-year overall survival. Furthermore, in a finite sample, the probability is commonly obtained by estimating the frequency of the sample. This estimate gradually approaches the true value as the sample size increases. As the ICGC cohort only had 76 PAAD samples (survival time >30 days), this may have also affected the AUC value for overall survival.

There were several limitations to the present study. First, the FRLS was only tested and validated in the TCGA and ICGC cohorts. If possible, it should be verified in other, larger cohorts to improve the accuracy and robustness of the prognostic prediction model. Furthermore, this was a retrospective study and it would be more convincing to evaluate the clinical application of the signature in prospective cohorts. Finally, LINC01133 and CASC8 were highly expressed in PAAD cells and tissues; however, they were not confirmed targets in PAAD treatment. Additional in vitro and in vivo experiments are required to fully characterize the role and mechanism of these lncRNAs in PAAD ferroptosis. In short, the prognostic model constructed in the present study requires further validation.

In conclusion, the present study constructed an FRLS for PAAD to predict prognosis. The prognostic signature and nine ferroptosis-related lncRNAs may be molecular biomarkers and therapeutic targets for patients with PAAD.

Supplementary Material

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by a grant from The Postdoctoral Program of Chongqing (grant no. 2021XM2017).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

JL, BN and YL designed the study. JL and WL acquired and analyzed the data. HW collected clinical samples and revised the manuscript. BN and YL wrote and revised the manuscript. All authors reviewed and approved the final manuscript. BN and YL confirm the authenticity of all the raw data.

Ethics approval and consent to participate

The experimental procedures for human tissues were approved by the ethics committee of Southwest Hospital (Chongqing, China; approval no. 20210312).

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References