Identification of a prognostic ferroptosis-related lncRNA signature in the tumor microenvironment of lung adenocarcinoma

Ferroptosis is closely linked to various cancers, including lung adenocarcinoma (LUAD); however, the factors involved in the regulation of ferroptosis-related genes are not well established. In this study, we identified and characterized ferroptosis-related long noncoding RNAs (lncRNAs) in LUAD. In particular, a coexpression network of ferroptosis-related mRNAs and lncRNAs from The Cancer Genome Atlas (TCGA) was constructed. Univariate and multivariate Cox proportional hazards analyses were performed to establish a prognostic ferroptosis-related lncRNA signature (FerRLSig). We obtained a prognostic risk model consisting of 10 ferroptosis-related lncRNAs: AL606489.1, AC106047.1, LINC02081, AC090559.1, AC026355.1, FAM83A-AS1, AL034397.3, AC092171.5, AC010980.2, and AC123595.1. High risk scores according to the FerRLSig were significantly associated with poor overall survival (hazard ratio (HR) = 1.412, 95% CI = 1.271–1.568; P < 0.001). Receiver operating characteristic (ROC) curves and a principal component analysis further supported the accuracy of the model. Next, a prognostic nomogram combining FerRLSig with clinical features was established and showed favorable predictive efficacy for survival risk stratification. In addition, gene set enrichment analysis (GSEA) revealed that FerRLSig is involved in many malignancy-associated immunoregulatory pathways. Based on the risk model, we found that the immune status and response to immunotherapy, chemotherapy, and targeted therapy differed significantly between the high-risk and low-risk groups. These results offer novel insights into the pathogenesis of LUAD, including the contribution of ferroptosis-related lncRNAs, and reveal a prognostic indicator with the potential to inform immunological research and treatment.


INTRODUCTION
Lung cancer is one of the leading causes of cancer-related deaths worldwide [1,2]. Although integrative therapies such as molecular targeted therapy, chemotherapy, and radiotherapy have been developed for lung adenocarcinoma (LUAD), the 5-year overall survival (OS) rate is only 15% [3,4]. Thus, it is essential to identify effective diagnostic markers, therapeutic targets, and prognostic factors.
Ferroptosis is a newly identified type of programmed cell death that is distinguished from apoptosis, necrosis, pyroptosis, and autophagic cell death [5,6]. It induces cell injury or death via the iron-dependent lipid peroxidation process [7][8][9]. Ferroptosisrelated genes have been identified as diagnostic markers and could even be potential drug targets, especially for the development and pharmacological study of anticancer agents and cancer therapies [10]. Ferroptosis plays an important role in the development of lung cancer [11][12][13][14][15]. However, its precise roles in the origin and progression of LUAD remain unclear.
Furthermore, the role of ferroptosis in the tumor immune microenvironment is not well established. Recent studies have demonstrated that immunotherapy-activated CD8+ T cells enhance ferroptosis-specific lipid peroxidation in tumor cells and increase the efficacy of cancer immunotherapy [16,17]. Intriguingly, IFNγ released from CD8+ T cells downregulates the expression of SLC7A11 and SLC3A2, thereby limiting cystine uptake by tumor cells and promoting tumor cell lipid peroxidation and ferroptosis [17]. This indicates that IFNγ is a vital mediator of tumor immune evasion and a potential target for improving clinical responses to immunotherapy. Furthermore, Lai et al. [11] identified a relationship between ferroptosis and non-small cell lung carcinoma prognosis. However, a comprehensive understanding of the mechanisms underlying the effects of ferroptosisrelated genes in LUAD and their functions is lacking.
Long noncoding RNAs (lncRNAs), which are >200 nucleotides in length and have limited protein coding ability, play essential roles in the development and progression of diverse cancers, including LUAD. The dysregulation of lncRNAs has been implicated in tumor cell proliferation, invasion, metastasis, and chemoresistance in lung cancer [18][19][20][21]. In addition, lncRNAs have been a major focus of research into ferroptosis [20,22,23]. For example, the p53related lncRNA P53RRA is involved in the regulation of SLC7A11, which is critical for the sensitivity of ferroptotic responses [23].
In the present study, we constructed a ferroptosis-related lncRNA signature and systematically assessed the correlations of ferroptosis-related lncRNAs with the prognosis and clinicopathological features of LUAD patients. We then established a nomogram that incorporates the ferroptosis-related lncRNA signature (FerRLSig) and clinical factors to predict the survival of these patients. The high-risk and low-risk groups were compared with respect to various factors, including immune status and response to immunotherapy. Our results provide valuable insights into the regulatory mechanisms by which ferroptosis contributes to LUAD and may help to improve the effectiveness of individualized treatment and assessments of prognosis.
Evaluation of FerRLSig as an independent prognostic factor for LUAD Next, univariate and multivariate Cox regression analyses were performed to determine whether FerRLSig is an independent prognostic model for OS in patients with LUAD. The HRs (95% CI) for the risk score were 1.456 (1.323-1.602) in the univariate Cox regression analysis (P < 0.001, Fig. 4A) and 1.412 (1.271-1.568) in the multivariate Cox regression analysis (P < 0.001, Fig. 4B), showing that FerRLSig was an independent prognostic indicator. In addition, the predictive accuracy of the model was assessed by a time-dependent receiver operating characteristic ROC analysis at 1, 3, and 5 years, with area under the curve (AUC) values of 0.781, 0.712, and 0.761, respectively (Fig. 4C). Then, we performed a PCA to compare the low-and high-risk groups based on the whole genome, ferroptosisrelated lncRNAs, and the risk model. As shown in Fig. 4D and E, the high-and low-risk groups could not be effectively discriminated using the whole genome or ferroptosis-related lncRNAs; however, using FerRLSig, the high-and low-risk patients could be clearly distinguished, further supporting the accuracy of the model. These results indicated that FerRLSig is a significant independent prognostic risk factor for patients with LUAD. Correlations between the risk score and clinicopathological factors To further assess the roles of FerRLSig in the development of LUAD, we assessed the correlations between the risk score and clinicopathological factors. As shown in Fig. 5 and Table S3, there were significant correlations between the risk score and pathologic stage (P < 0.05). In particular, the risk score was significantly higher for stages II-IV than for stage I ( Fig. 5A, P < 0.001), and the signature was associated with tumor stage (Fig. 5B, P < 0.001) and for patients with lymph node metastasis than for those without lymph node metastasis ( Fig. 5C, P < 0.001). Figure 5D illustrates that patients with a high-risk score (i.e., >1.17) had a significantly poorer prognosis in terms of survival status than patients with low-risk scores (P < 0.0001). These results indicate that FerRLSig is closely associated with LUAD progression and prognosis.

Construction of a predictive nomogram
We constructed a clinically adaptable nomogram to estimate the 1-, 3-, and 5-year survival probabilities of patients with LUAD using FerRLSig combined with other clinicopathological factors (Fig. 6A). The calibration plots of the nomogram for 1-year, 3-years, and 5years ( Fig. 6B-D) indicated that the mortality estimated by the nomogram was close to the actual mortality. Time-dependent ROC curves of 5-year OS were generated, and the AUC value for the clinical prognostic nomogram was 0.766, which was significantly higher than those for age, sex, tumor stage, tumor T stage, tumor M stage, and tumor N stage, further supporting the discriminative ability of FerRLSig in conjunction with clinicopathological factors for predicting survival in LUAD (Fig. 6E).  (Table S5, Fig. 7B). These results indicate that the lncRNAs in the newly developed signature may be associated with the tumor immune microenvironment. Differences in the response to immunotherapy, chemotherapy, and targeted therapy between high-risk and low-risk groups To assess the correlations between the risk score and tumor immune cell infiltration, the CIBERSORT algorithm was applied to compare the proportions of 22 immune cell types between the high-risk and low-risk groups. Differences in the infiltration of 22 immune cell types among patients with LUAD from TCGA are shown in Fig. S1A, showing that this is an intrinsic feature reflecting individual differences. The proportions of immune cells   < 0.001). B Ferroptosis-related lncRNAs in the cohorts stratified by tumor size (P < 0.001). C Ferroptosis-related lncRNAs in the cohorts stratified by regional lymph nodes (P < 0.001). D Ferroptosis-related lncRNAs in the cohorts stratified by survival outcome (P < 0.001). T, tumor size; N, regional lymph node metastasis.

Identification of FerRLSig-associated biological pathways
Y. Guo et al. differed between the two risk groups (Fig. S1B, where colors ranging from green to red represent the infiltration density from low to high). Patients with LUAD in the high-risk group had higher ratios of activated memory CD4 T cells (P < 0.001), M0 macrophages (P < 0.001), M1 macrophages (P = 0.002), and activated mast cells (P = 0.005). In contrast, memory B cells (P = 0.047), plasma cells (P = 0.023), monocytes (P = 0.001), resting dendritic cells (P < 0.001), resting mast cells (P < 0.001), and monocytes (P = 0.002) were negatively correlated with the risk score (Fig. 8A). The substantial immune infiltration observed in the low-risk group partly reflected the reduction in malignancy and effects of various treatments, indicating that our signature is not only a prognostic marker but also reflects levels of immune cell infiltration.
In addition, we evaluated the associations between the risk score and the efficacy of targeted therapeutics and chemotherapeutics for LUAD. Patients with low-risk scores were highly sensitive to the targeted therapeutics gefitinib (P = 0.008) and erlotinib (P < 0.001). In addition, patients with low-risk scores were highly sensitive to the chemotherapeutics cisplatin (P = 0.0068), etoposide (P < 0.001), paclitaxel (P < 0.001), docetaxel (P < 0.001), and gemcitabine (P = 0.001), suggesting that FerRLSig is a potential predictor of chemosensitivity (Fig. S3). These results could also explain the poor prognosis of patients with high risk scores, as these patients are likely to be more resistant to chemotherapy and targeted therapy.

DISCUSSION
Lung cancer is the leading cause of cancer-related mortality worldwide. The pathogenesis of LUAD is unclear, contributing to the lack of effective treatments [24,25]. Compared with single clinical biomarkers, the integration of multiple biomarkers into a single model can improve predictive accuracy and assist in the development of individualized treatment plans.
Ferroptosis is a novel iron-dependent and nonapoptotic form of cell death that is mediated by oxidative modifications of membrane polyunsaturated phospholipids. As an adaptive mechanism to eliminate malignant cells, ferroptosis is a promising Y. Guo et al.
new target for cancer treatment [26,27]. Battaglia et al. [28] reported that ferroptosis plays an important role in cancer. Various lncRNAs and microRNAs (miRNAs) involved in the regulation of ferroptosis have been identified [18,29]. For example, the lncRNA P53RRA binds to Ras GTPase-activating protein-binding protein 1 (G3BP1) and prevents its interaction with p53 to induce cell cycle arrest, apoptosis, and ferroptosis in lung cancer [20].
Although several lncRNA signatures associated with lung cancer have been reported, previous studies have focused on relapse-free survival or OS for patients with early-stage lung cancer [30,31], and general prognostic models for lung cancer remain warranted. In the present study, we constructed a prognostic model using 10 ferroptosis-related lncRNAs (AL606489.1, AC106047.1, LINC02081, AC090559.1, AC026355.1, FAM83A-AS1, AL034397.3, AC092171.5, AC010980.2, and AC123595.1). Based on the ROC curve, the lncRNA signature had a moderate predictive performance for OS. In addition, our newly developed nomogram is expected to improve clinical decision-making and guide the development of treatment strategies.
The most significant contribution of this research is the demonstration of the relationship between FerRLSig and the tumor immune microenvironment. Notably, the complex interplay between tumor cells and the tumor microenvironment not only plays a pivotal role in tumor development but also has significant effects on immunotherapeutic efficacy and overall survival [36,37]. In this study, a functional enrichment analysis indicated that ferroptosis-related lncRNAs are primarily involved in immune pathways. Patients with high-risk scores had higher proportions of activated memory CD4 T cells, M0 macrophages, M1 macrophages, and mast cells, confirming the roles of ferroptosis-related lncRNAs in the regulation of tumor immune infiltration. Since our results link FerRLSig to immune infiltration in LUAD, these ferroptosis-related lncRNAs may be targets for combined treatments with immune checkpoint inhibitors.
The combination of immune checkpoint blockade with immunotherapies such as CTLA-4, PD-1, and PD-L1 inhibitors is a promising approach to treat a variety of malignancies, and an activated tumor immune microenvironment is correlated with good outcomes of immune checkpoint inhibitor treatment [38,39]. Interestingly, PD-L1 was highly expressed in the highrisk group, indicating that patients with high-risk scores may benefit more from anti-PD-L1 immunotherapy, while CTLA-4 was highly expressed in the low-risk group, indicating that patients in the low-risk group may benefit more from anti-CTLA-4 immunotherapy. This provides new insights into tumor immunotherapy.
The current study had several limitations. First, we used a single data source. Second, it was a retrospective study. Third, wellknown prognostic factors such as chemotherapy data and tumor markers were not included in the nomogram since data for these parameters were incomplete. Thus, additional prospective studies are needed to verify the prognostic value of FerRLSig. In addition, functional experiments should be conducted to further elucidate the molecular mechanisms underlying the effects of ferroptosisrelated lncRNAs.
In summary, we constructed a ferroptosis-related lncRNA signature for LUAD to predict prognosis. We established an effective nomogram including FerRLSig. Importantly, the FerRLSig generated in our study might be associated with immune infiltration levels and even the efficacy of tumor immunotherapy.

MATERIALS AND METHODS Datasets and sample extraction
RNA-sequencing (RNA-seq) data for LUAD were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). A total of 535 patients with LUAD were included in the study. In addition, complete clinical information was downloaded from TCGA. A total of 477 patients with complete follow-up information and survival time >30 days and 325 patients with complete clinicopathological data were included in subsequent analyses.

Bioinformatics analysis
Coexpression networks were visualized using Cytoscape 3.7.2, while principal component analysis (PCA) was conducted to evaluate the distribution of patients with different risk scores. The R package scatterplot3D was used to generate the PCA plots. The nomogram was developed and validated using the rms package in R version 4.0.3 (http:// www.r-project.org/); the R package ggalluvial was used to obtain the Sankey map.

Identification of a prognostic ferroptosis-related lncRNA model for LUAD
A total of 1621 ferroptosis-related lncRNAs were identified by constructing a ferroptosis-related mRNA-lncRNA coexpression network using |Pearson correlation coefficient | > 0.3 and P < 0.001 as thresholds. The coexpression networks were visualized using Cytoscape 3.7.2. Then, based on univariate analyses with P < 0.01 as the threshold, significant prognostic ferroptosisrelated lncRNAs were identified and subsequently incorporated into a Fig. 7 Enrichment plots from gene set enrichment analysis (GSEA). A GSEA results showing differential enrichment of genes in GO with ferroptosis-related lncRNAs (5 GO items, namely, positive regulation of activation of innate immune response, innate immune response activating signal transduction, interleukin 1 mediated signaling pathway, positive regulation of innate immune response, and regulation of apoptotic signaling pathway), showed significantly differential enrichment in the high expression phenotype. Five GO items, namely, CD8positive alpha beta T cell activation, mast cell-mediated immunity, regulation of leukocyte-mediated immunity, regulation of lymphocytemediated immunity, and T cell-mediated immunity, showed significantly differential enrichment in the low expression phenotype. B GSEA results showing differential enrichment of genes in KEGG with ferroptosis-related lncRNA expression. (5 KEGG items namely, positive regulation of pancreatic cancer, cell cycle, p53 signaling pathway, pathogenic Escherichia coli infection, and small cell lung cancer were significantly differential enrichment in high expression phenotype; 5 KEGG items namely, ferroptosis -related lncRNA, FC epsilon RI signaling pathway, autoimmune thyroid disease, and allograft rejection and graft versus host disease showed significantly differential enrichment in the ferroptosis-related lncRNA low-expression phenotype based on the normalized enrichment score (NES), nominal p value (NOM P value), and FDR value).
multivariate Cox regression analysis to establish risk scores. The risk score for each patient was calculated using the following formula: risk score = ∑i coefficient(lncRNA1) × expression(lncRNA1) + coefficient(lncRNA2) × expression(lncRNA2) + …… + coefficient(lncRNAn) × expression(lncRNAn). A linear regression analysis was used to evaluate the relationship between survival and high-risk and low-risk groups. The consistency index (c-index), calibration curve, and receiver operating characteristic (ROC) curves were used to explore the accuracy of the model. Demographic data were included in multivariate Cox regression analysis to confirm whether the risk score is an independent prognostic indicator. The distribution of survival statuses according to risk score levels was evaluated.

Statistical analysis
All computational and statistical analyses were conducted using R (version 4.0.3). Survival analysis was performed using the Kaplan-Meier method.
The Kruskal-Wallis test was used to compare differences between groups. The chi-squared test or Fisher's exact test was used for analyses of clinical information. Spearman or Pearson correlation coefficients were used to evaluate the relationships among lncRNA expression, immune infiltration, and immune checkpoint gene expression.

Significance of the model for the prediction of clinical treatment response
Levels of immune cell infiltration were quantified in each LUAD sample using CIBERSORT (http://cibersort.stanford.edu/). The relationships between the FerRLSig score and the expression levels of ICI genes were determined using the ggstatsplot package in R, and violin plots were generated for visualization. IC 50 values for various antitumor drugs recommended for lung cancer treatment, such as cisplatin, etoposide, Fig. 8 Comparison of the immune microenvironment and immune checkpoints with LUAD between the high-risk and low-risk groups. A Violin plot of immune-infiltrating lymphocytes between the low-risk and high-risk groups, in which orange indicates high-risk samples and blue indicates low-risk samples. B Volcano plots for the enrichment of expression of immune checkpoints between the low-risk and high-risk groups. The differential expression of six immune checkpoints, (C) CD4, (D) CD44, (E) CD27, (F) CD48, (G) CTLA4, and (H) PD-L1, between the high-risk group and the low-risk group. Orange indicates high-risk patients, and blue indicates low-risk patients.
docetaxel, gefitinib, erlotinib, gemcitabine, and paclitaxel, were compared between groups using the Wilcoxon signed-rank test using pRRophetic and ggplot2 in R.

Gene set enrichment analysis
Gene set enrichment analysis (GSEA) was performed for the high-risk and low-risk groups according to the prognostic model. Significant predefined biological processes and pathways were enriched when NOM P < 0.05 and false discovery rate (FDR) < 0.25. c5.go. v7.1 symbols.gmt and c2.cp.kegg. v7.1.symbols.gmt were selected from the Molecular Signatures Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp) as the reference files.

DATA AVAILABILITY
The datasets generated and analyzed during the current study are available in the Figshare (https://figshare.com/s/f34a12c9a77534612ea7) repository.