A Novel Autophagy-Related Long Non-Coding RNA Signature to Predict Prognosis and Therapeutic Response in Esophageal Squamous Cell Carcinoma

Introduction

Esophageal carcinoma is the most common malignant tumor of the digestive system, which includes two principal pathological subtypes: squamous cell carcinoma and adenocarcinoma.¹ In China, the main pathological type of esophageal cancer is esophageal squamous cell carcinoma (ESCC), accounting for more than 90% of the total cases.² In the past few decades, although the diagnosis and treatment strategies of ESCC have made great progress, the overall survival rate is still poor, with a 5-year survival rate of only about 15–25%.^3–5 Hence, the development of new effective prognostic biomarkers to guide ESCC treatment is imperative.

Autophagy is a physiological process during which cells use lysosomes to degrade damaged organelles and macromolecular substances.⁶ Currently, autophagy appears to play dual roles in tumorigenesis, maintenance, and tumor progression. In the process of tumor initiation and malignant transformation, autophagy exerted a tumor-suppressive role. However, in the process of tumor progression, autophagy can act as an oncogenic mechanism by promoting tumor cell metastasis and inhibiting tumor cell apoptosis.⁷ Current literature shows that autophagy is associated with the diagnosis, treatment and prognosis of ESCC patients.⁸

Long non-coding RNAs (lncRNAs) generally refer to the non-coding RNAs, which are at least 200 nucleotides in length.⁹ With the continuous deepening of research in recent years, lncRNAs have engaged much attention. There is a growing literature that lncRNAs play a key role in gene regulation, thus affecting a variety of biological processes, including tumor cell proliferation, differentiation, metastasis and apoptosis.^10,11 Recent studies also indicate that lncRNA can play a regulatory role in complex networks of autophagy by interacting with multiple autophagy-related DNA, RNA, or proteins.¹² For example, LINC00337 can induce autophagy and chemotherapeutic resistance to cisplatin in ESCC cells by upregulating the expression of TPX2.¹³ However, studies on the relationship between ESCC and autophagy-related lncRNAs (ARlncRNAs) are particularly scarce, which still need to be further elucidated.

Here, we comprehensively analyzed the expression profile of ARlncRNAs in 179 patients with ESCC to construct and validate an autophagy-related prognostic model. The present study aimed to enable patients with ESCC to benefit from individualized treatment and thus improve prognosis.

Materials and Methods

Data Preparation

We downloaded the RNA expression profiles and clinical data of ESCC from the GSE53625 dataset and The Cancer Genome Atlas (TCGA) database, respectively. The GSE53625 dataset contained normalized gene expression data of 179 ESCC patients, but the batch effect was still present, so the “limma” package was employed to remove the batch effect (Figure S1 and Table S1).^14,15 RNA-sequence data of TCGA-ESCC samples were analyzed as previously described.¹⁶ 222 autophagy-related genes (ARGs) were searched from the Human Autophagy Database (HADB, http://www.autophagy.lu/). Next, we adopted the Spearman correlation test to determine potential ARlncRNAs and the filter criteria for correlation coefficient and p-value were set at 0.3 and 0.001, respectively. Owing to the freely accessible resources TCGA and GEO databases, the study let off institutional review board approval.

Differentially Expressed ARlncRNAs

By using the “limma” package to compare the matched tumors and normal tissues of 179 ESCC patients, we screened differentially expressed autophagy-related lncRNAs (DEARlnRNAs) with both |log2FC| more than 1 and false discovery rate (FDR) less than 0.05.

Construction of Prognostic Signature

DEARlncRNAs were analyzed by univariate Cox regression analysis, and lncRNAs with p < 0.05 were defined as prognostic biomarkers. In order to avoid overfitting, we adopted the least absolute shrinkage and selection operator (LASSO) and multivariate Cox analysis method based on Akaike Information Criteria (AIC) to narrow the range of lncRNA, thus selecting the optimal prognostic risk model. Furthermore, a prognostic model was established by multiplying multivariate Cox regression coefficients by the expression level of each variable.

Validation and Assessment of the Prognostic Model

The risk score of each patient with ESCC was calculated according to the prognostic formula, and all ESCC patients were separated into two groups (high-risk group or low-risk group) based on the median risk score. To evaluate the prognostic value of the model, we plotted the Kaplan–Meier survival curves and receiver operating characteristic (ROC) curves. Meanwhile, we assessed model calibration using the Brier score and calibration curves. The model was internally validated using the Bootstrap method based on 1000 resampling, and AUC and Brier scores for the prognostic model were calculated by the “riskRegression” package.¹⁷ Clinical correlation analyses were performed to compare whether there were differences in risk score between clinical variables and the survival prediction ability of prognostic factors was further compared. To demonstrate whether the prognostic model could be assumed as a standalone prognostic variable, univariate and multivariate Cox regression analyses were implemented by combining the risk score with clinicopathological variables, and the results were drawn forest maps.

Gene Set Enrichment Analysis

To explore the possible mechanisms between high and low-risk groups of ESCC patients, we employed the “clusterProfiler” package to perform GSEA based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) (v7.2) gene set collections.

Immune Infiltration and Therapeutic Response

The relative proportion of 22 tumor-Infiltrating Immune Cells (TICs) in 179 ESCC samples was evaluated with CIBERSORT algorithm.^18,19 After a quality screening (p < 0.05), 107 cases of ESCC could be used for subsequent analysis. Wilcoxon test was implemented to compare the different abundances of 22 TICs types between high and low-risk groups. To assess the effect of immunotherapy, we detected the correlation of risk scores with immunotherapy-related targets, such as PD-1 (programmed cell death 1) and CD274 molecule (best known as PD-L1) using the Spearman correlation method. Meanwhile, the “pRRophetic” package was adopted to estimate the chemotherapeutic responses based on half maximal inhibitory concentration (IC50) of chemotherapy drugs for each ESCC patient available in the Genomics of Drug Sensitivity in Cancer (GDSC) database.

Relationship Between the Risk Score and Ferroptosis

Ferroptosis is a newly introduced form of programmed cell death discovered in recent years. Accumulating studies have revealed that ferroptosis can inhibit tumor growth or promote tumor proliferation in tumor development, and there is crosstalk with autophagy at the molecular level.^20,21 Then, Spearman correlation examined the correlation of risk scores with ferroptosis-associated genes that were identified from the FerrDb website (http://www.zhounan.org/ferrdb/index.html).

Construction of the ARlncRNA-mRNA Co-Expression Network

We analyzed the relationship between the nine ARlncRNAs and co-expression of ARGs through co-expression network and Sankey diagram. Furthermore, Gene ontology (GO) and KEGG pathway analysis of the co-expressed ARGs were performed using the “clusterProfiler” package.^22–24

Validation of Core lncRNAs

The expression levels of nine ARlncRNAs in tumor and normal tissues were verified in TCGA ESCC cohort. Furthermore, the ESCC samples from the TCGA cohort were categorized into two groups (the high and low-expression groups) based on the median expression level of a single lncRNA and K-M analysis was implemented to verify the survival prediction ability of core lncRNAs.

Statistical Analysis

All statistical analyses and visualization of results in this study were performed with R software version 4.0.2. Statistical significance was set at p-value <0.05. Wilcoxon test was used for two independent samples or the paired data. P < 0.05 and Benjamini-Hochberg adjusted p < 0.05 were considered as the cut-off criterion for enrichment analysis.

Results

Identification of 481 DEARlncRNAs

Figure 1 presents the overall procedures in this study. The expression profiles of 6615 lncRNAs and 189 ARGs were obtained from GSE53625. Then, Spearman correlation analysis between the lncRNAs and ARGs was performed and 5031 ARlncRNAs were identified. Next, a total of 481 DEARlncRNAs were identified between tumor tissues and normal tissues, of which 169 lncRNAs were upregulated and 312 lncRNAs were downregulated, respectively (Table S2). Volcano and heatmap plots are shown (Figure 2A and B).

Figure 1 Flow diagram of the overall procedures.

Figure 2 Selection of DEARlncRNAs and construction of prognostic signature. (A) Volcano plot of DEARlncRNAs between tumor tissues and normal tissues. Upregulation as red dot and downregulation as blue dot. (B) Heatmap of top 30 DEARlncRNAs between tumor tissues and normal tissues. (C) Distribution map of LASSO coefficients. (D) LASSO coefficient values and the selection of the best parameter (lambda). (E) Cox regression analysis of ARlncRNAs visualized as forest plot.

Construction and Validation of Prognostic Signature for ARlncRNAs

In total, 64 prognostic DEARlncRNAs in ESCC were extracted by univariate analysis of which the significance filtering condition was p < 0.05 (Table S3). Subsequently, the prognostic signature based on 9 ARlncRNAs was established using Lasso and stepwise multivariate Cox regression analysis (Figure 2C–E). The formula of the prognostic model was given as follows: Risk score = [0.21219 × LINC02024] + [0.46308 × LINC01711] + [- 0.10103 × LINC01419] + [- 0.11997 × LCAL1] + [0.19070 × FENDRR] + [0.15634 × ADAMTS9-AS1] + [- 0.12891 × AC025244.1] + [- 0.10636 × AC015908.6] + [0.49029 × AC011997.1]. Furthermore, we found that LINC02024, LINC01711, FENDRR, ADAMTS9-AS1 and AC011997.1 were risk factors for HR > 1, while LINC01419, LCAL1, AC025244.1 and AC015908.6 were favorable factors with HR < 1 (Table 1 and Figure 3).

Table 1 Multivariate Cox Results of Prognostic ARlncRNAs Based on GSE53625 Dataset

Figure 3 Kaplan-Meier survival curves of 9 ARlncRNAs identified in GSE53625 dataset. (A) LINC02024; (B) LINC01711; (C) LINC01419; (D) LCAL1; (E) FENDRR; (F) ADAMTS9-AS1; (G) AC025244.1; (H) AC015908.6; (I) AC011997.1.

We computed each patient’s risk score using the prognostic signature and separated all ESCC patients into a high-risk group (n = 89) and low-risk group (n = 90) based on the median risk score. According to the K-M survival curve, patients in the low-risk group had a longer survival time than those in the high-risk group (p < 0.001, Figure 4A). By plotting ROC curves, we found that AUC values of 1-, 3-, and 5-year survival rates in the original dataset were 0.755, 0.772, and 0.796, respectively, which illustrated the medium predictive accuracy of the prognostic model (Figure 4B). Figure 4C presents the distribution of risk scores, the survival status and the expression values of nine lncRNAs of patients. The calibration curves of the prognosis model for the prediction of 1-, 3- and 5-year OS rates demonstrated good agreement in the original dataset with Brier scores of 0.151, 0.192, and 0.18, respectively (Figure 4D). Through bootstrapping verification, the AUCs and Brier scores for predicting 1-, 3- and 5-year OS were 0.708, 0.715, 0.738, 0.166, 0.222 and 0.208, respectively, indicating that the model had good discrimination and calibration (Figure 4E).

Figure 4 Appraisal of prognostic signature in179 ESCC patients. (A) Kaplan-Meier survival analysis in the original dataset. (B) ROC curves for 1-, 3-, and 5-year survival in the original dataset. (C) Risk score analysis (risk curve; scatter plot of vital status by risk score; heatmap of the 9 lncRNAs expression). (D) Calibration curves of prognostic signature in the original dataset. (E) Calibration curves for internal validation by the Bootstrap method.

Independence of the Model

We conducted a correlation analysis between the risk score and clinicopathological features, and found that the risk score was related to survival outcomes (p < 0.001), T stage (p = 0.022, p = 0.013) and grade (p = 0.02) (Figures 5A–C and S2). Univariate Cox regression indicated that advanced age, higher tumor stage, N stage, grade, and higher risk score were risk factors for poor survival (Figure 5D). Furthermore, we performed a multivariate Cox regression analysis and found that the risk score could serve as an independent prognostic factor in ESCC patients, whereas the clinical parameters were not (Figure 5E). By means of ROC curves with multiple measures, we visualized that the AUC values of risk score were significantly higher than those of other clinical parameters, which furnished evidence to support that the foretelling effect of the prognostic signature had medium accuracy (Figure 5F).

Figure 5 Clinical significance and independent prognostic analysis of prognostic signature. (A) survival outcome; (B) T stages; (C) grade; (D) The univariate and (E) multivariate Cox regression analysis of risk score and clinical parameters by using overall survival (OS) of ESCC patients. Age was analyzed as a categorical variable: age≥60 =1, age<60=0; Tumor stage, T stage, N stage, grade, and risk score were analyzed as continuous variables. (F) ROC analysis of 5-year OS for the model and the clinical parameters.

Gene Set Enrichment Analysis

The GSEA results revealed that cancer- and autophagy-related signaling pathways were obviously enriched in the high-risk ESCC patients, including ECM receptor interaction, TGF- β signaling pathway and focal adhesion (Figure 6A). Meanwhile, metabolism-related signaling pathways were obviously enriched in the low-risk ESCC patients, including Arachidonic acid metabolism, Retinol metabolism, linoleic acid metabolism and drug metabolism cytochrome p450 signaling pathway (Figure 6B). Details of the enrichment analysis are given in Table S4.

Figure 6 Gene set enrichment analysis based on risk scores. (A) the result of high-risk group; (B) the result of low-risk group.

Immune Cell Infiltration and Therapeutic Response Analysis

Figure 7A shows the relative content distribution of 22 TICs in the GSE53625 cohort. After quality screening (p < 0.05), 107 samples (51 low-risk samples and 56 high-risk samples) of ESCC were then used to analyze the immune cells infiltration. We found that the infiltration level of plasma cells was upregulated, while the infiltration level of macrophages M2 cells was downregulated in the low-risk group compared with the high-risk group (Figure 7B). Then, correlation analysis revealed the mRNA expression levels of PD-1 (R = −0.25, p < 0.001) were negatively related to the risk score (Figure 7C), while no statistical difference was shown between the mRNA expression levels of PD-L1 and risk score (Figure 7D). Meanwhile, the results regarding the chemosensitivity of three commonly used drugs for ESCC treatment (cisplatin, paclitaxel and docetaxel) indicated that the estimated IC50 values of cisplatin were lower in the high-risk group (Figures 7E and S3). The above results revealed that the signature may be a valuable indicator to assess patients’ response to immunotherapy and chemotherapy.

Figure 7 Correlations between risk score and tumor-infiltrating immune Cell (TIC) and response to immunotherapy and chemotherapy. (A) The relative content distribution of 22 TICs of ESCC samples. (B) Differential analysis of immune cells proportions in the high and low-risk groups. (C and D) Correlation plots of the risk score and the expression levels of PD-1 and PD-L1. (E) Chemotherapeutic sensitivity of cisplatin in the high and low-risk groups.

Identification of the Ferroptosis Correlation with the Prognostic Signature

Among the 209 ferroptosis-related genes, 74 (35.41%) genes were significantly associated with risk scores, of which 32 and 42 genes separately were positively and negatively associated with risk scores. As shown in Figure 8, HIC1, ACSL3, NNMT, ANO6 and CDO1 are the leading five ferroptosis-related genes that have a positive association with the risk score, while the top five ferroptosis-related genes that are negatively associated with the risk score are HAMP, ALOX12B, MAFG, HILPDA and DUOX1.

Figure 8 Correlations between risk score and ferroptosis-related genes in ESCC. (A) The top five positive correlations. (B) The top five negative correlations.

ARlncRNA-mRNA Co-Expression Network

We developed a co-expression network containing nine ARlncRNAs and 42 ARGs and the Sankey plots indicated the correlation between them and risk types (protective or risk) (Figure 9A and B). GO biological process enrichment analysis of 42 ARGs showed that they were mainly involved in the autophagy, process utilizing autophagic mechanism, cellular response to starvation and macroautophagy (Figure 9C). KEGG analysis revealed that these ARGs mainly participated in the Autophagy-animal, Autophagy-other, FoxO and Cellular senescence signaling pathway (Figure 9D). Details of the enrichment analysis are given in Table S5.

Figure 9 ARlncRNA-mRNA co-expression network. (A) Co-expression network contained 9 ARlncRNAs and 42 autophagy mRNAs. LncRNA as red rectangle and mRNA as blue ellipse. (B) Sankey diagram of genes in co-expression network and risk types. (C and D) GO and KEGG enrichment analysis of 42 autophagy-related genes.

Validation of Core lncRNAs

We extracted the expression levels of 9 ARlncRNAs in tumor tissues and normal tissues from GSE53625 dataset and found that LINC01711, LINC01419, LCAL1, AC025244.1 and AC011997.1 were highly expressed in ESCC, while LINC02024, FENDRR, ADAMTS9-AS1 and AC015908.6 were lowly expressed in ESCC (Figure 10). Furthermore, we extracted the expression profiles of six of these lncRNAs from TCGA ESCC dataset. Among them, 4 lncRNAs (LINC01711, LINC01419, AC025244.1 and AC011997.1) were up-regulated in ESCC, and 2 lncRNAs (FENDRR and ADAMTS9-AS1) were down-regulated in ESCC, which was consistent with the differential expression pattern in the GSE53625 dataset (Figure 11A–F). In addition, K-M plots demonstrated that three lncRNAs (LINC01419, FENDRR, and AC011997.1) had similar prognostic values as those in the GSE53625 dataset (Figure 11G–I).

Figure 10 The expression profiles of 9 ARlncRNAs in tumor tissues and normal tissues from GSE53625 dataset. (A) LINC02024; (B) LINC01711; (C) LINC01419; (D) LCAL1; (E) FENDRR; (F) ADAMTS9-AS1; (G) AC025244.1; (H) AC015908.6; (I) AC011997.1.

Figure 11 The expression profiles and Kaplan-Meier survival curves of ARlncRNAs in TCGA ESCC dataset. The expression profiles of ARlncRNAs in tumor tissues and normal tissues ((A) LINC01419; (B) LINC01711; (C) FENDRR; (D) ADAMTS9-AS1; (E) AC025244.1; (F) AC011997.1). Kaplan-Meier survival curves of ARlncRNAs ((G) LINC01419; (H) FENDRR; (I) AC011997.1).

Discussion

In China, ESCC is the most common pathological subtype of esophageal cancer, and patients with ESCC often have a short survival due to unclear pathogenesis and insensitivity to treatment. Considering the significance of autophagy and lncRNA in the biology of ESCC and the fact that several studies have shown that models composed of ARlncRNAs could indicate prognosis in multiple cancer types,^25–28 the development of a robust prognostic biomarker composed of ARlncRNAs to guide ESCC treatment is imperative.

In our present study, we constructed a novel prognostic model composed of nine ARlncRNAs, which showed good prognostic performance and was a significant independent prognostic factor for predicting OS in ESCC patients. The model was subjected to bootstrap internal validation on the original GEO dataset, and the expression levels of core lncRNAs were verified on the TCGA ESCC dataset. Correlation analysis showed that the risk score of ARlncRNAs model was associated with survival outcomes, T stage, and tumor grade. The GSEA results revealed the prognostic signature may be involved in the cancer-related signaling pathways, autophagy-related signaling pathways and metabolism-related signaling pathways. Correlation analysis showed that the risk score was negatively correlated with the level of PD-1mRNA expression, and the risk score could predict the sensitivity of cisplatin therapy. In summary, our research demonstrated that autophagy dysfunction was important during the process of ESCC. Most importantly, our study was the first to demonstrate that prognostic signature composed of nine ARlncRNAs had predictive value for survival prognosis, tumor-infiltrating immune cells, immunotherapy response, and drug sensitivity of ESCC.

Among the nine autophagy-related lncRNAs that make up the prognostic model, only LINC01419, LCAL1, FENDRR and ADAMTS9-AS1 have been reported in ESCC or other cancers. Existing literatures have shown that LINC01419 is overexpressed and associated with malignant phenotypes in hepatocellular carcinoma, gastric cancer and ESCC.^29–31 Li et al revealed that LCAL1 was elevated in lung cancer tissues, and overexpressed LCAL1 may suppress autophagic cell death via AMPK/ULK1 pathway.³² Consistent with our results, the expression of FENDRR was significantly down-regulated in liver cancer, colon cancer and gastric cancer.^33–35 Xu et al indicated that low FENDRR expression in gastric cancer was related to poor prognosis and this study disclosed that low expression of FENDRR predicted good prognosis in ESCC, which may require further investigation.³⁶ ADAMTS9-AS1 has been reported as a prognostic risk factor in several tumors, including ESCC, prostate cancer and colorectal cancer.^37–39

CIBERSORT result indicated that the risk score was related to immune cell infiltration (plasma cells and macrophages M2 cells) in ESCC. Plasma cells can produce antibodies and thus improve the immune response, which may explain the upward trend of plasma cells in the infiltrating abundance of the low-risk group in this article.⁴⁰ M2 macrophages have been demonstrated to exhibit pro-tumoral activity, and have been related to poor prognosis of ESCC.^41,42 Since the efficacy of drugs is related to the degree of drug sensitivity of patients and individual differences, it is essential to detect the drug sensitivity for patients requiring immunotherapy and chemotherapy. Our study indicated that the mRNA expression level of PD-1 was negatively related to the risk score, and the estimated IC50 values of cisplatin were lower in the high-risk group. These data suggest that further studies will be interesting to explore specific mechanisms of risk signature and tumor-infiltrating Immune Cells, and the risk score may provide potential guidance for personalized treatments of ESCC.

Ferroptosis is a newly introduced form of programmed cell death discovered in recent years. Accumulating studies have revealed ferroptosis could inhibit tumor growth or promote tumor proliferation in tumor development, and there is crosstalk with autophagy at the molecular level.^20,21 In this study, we found that 35.41% (74/209) ferroptosis-related genes were significantly associated with risk scores, which may further provide more information for targeting ferroptosis therapy.

Nonetheless, there are still some limitations in this study. First, although we have used the bootstrap method to verify the robustness of the prognostic model established in 179 ESCC patients, and the expression levels of key lncRNAs have been verified on TCGA ESCC dataset, external datasets are still needed to further validate the prognostic signature. Second, the role and mechanism of lncRNAs that constitute the prognostic signature of esophageal squamous cell carcinoma needs to be explored in depth.

Conclusion

In brief, we constructed a novel prognostic signature composed of nine autophagy-related lncRNAs for ESCC and verified its reliability, which might provide a novel approach for effective prediction of clinical outcomes and individualized therapeutic strategies selection.