Establishment and Analysis of a Prognostic Model of Autophagy-Related lncRNAs in ESCA.

Esophageal cancer (ESCA) is a malignant tumor of the upper gastrointestinal tract, with a high mortality rate and poor prognosis. Long noncoding RNAs (lncRNAs) play a role in the malignant progression of tumors by regulating autophagy. This study is aimed at establishing a prognostic model of autophagy-related lncRNAs in ESCA and provide a theoretical basis to determine potential therapeutic targets for ESCA. The transcriptome expression profiles were downloaded from The Cancer Genome Atlas (TCGA). We identified autophagy-related mRNAs and lncRNAs in ESCA using differential expression analysis and the Human Autophagy Database (HADb). Four differentially expressed autophagy-related lncRNAs with a prognostic value were identified using Cox regression and survival analyses. Furthermore, the combination of the selected lncRNAs was able to predict the prognosis of patients with ESCA more accurately than any of the four lncRNAs individually. Finally, we constructed a coexpression network of autophagy-related mRNAs and lncRNAs. This study showed that autophagy-related lncRNAs play an important role in the occurrence and development of ESCA and could become a new target for the diagnosis and treatment of this disease.

1. Introduction

Esophageal cancer (ESCA) is an aggressive malignant gastrointestinal tumor that affects the epithelial tissue of the esophagus [1], with a high mortality rate and poor prognosis, causing more than 500,000 deaths annually [2]. At present, the diagnosis of ESCA relies mainly on gastroscopy and pathological evidence. However, ESCA symptoms are often discrete in the early stages of the disease and are mostly diagnosed during the late and advanced stages. Despite treatments such as surgery and neoadjuvant chemotherapy, the survival rate of ESCA remains low [2]. Therefore, it is important to understand the mechanism underlying the occurrence and development of ESCA and identify biomarkers for its diagnosis and prognosis.

Autophagy is a survival-promoting pathway that functions in the capture, degradation, and circulation of intracellular proteins and organelles in the lysosomes [3]. Autophagy retains the functions of organelles, prevents the toxic accumulation of cellular waste products, and provides a substrate for maintaining metabolism during starvation [4]. The role of autophagy in cancer depends on the availability of nutrients, microenvironmental pressure, and the immune system. Although autophagy inhibits tumorigenesis in some cancers, it promotes this process in most cancers [5, 6]. Wu et al. demonstrated that the tight junction protein CLDN1 activates AMPK/STAT1/ULK1 signaling in the esophageal squamous cell carcinoma cell lines TE10 and TE11, which induces autophagy in esophageal squamous cell carcinoma cells, and enhances cancer cell proliferation and metastasis [7]. Several studies reported that autophagy is upregulated in hypoxic tumor areas, inhibits tumor-induced inflammation, promotes tumor cell survival, and increases growth and invasiveness [8, 9]. Another study reported altered expression of autophagy markers in patients with ESCA and a significant association between microtubule-associated light chain 3 (LC3), the most characteristic autophagy marker in ESCA, and poor survival in these patients [10]. Therefore, autophagy-related genes are potential targets for the treatment of ESCA and have broad prospects for clinical applications.

At present, the Human Genome Project has deciphered approximately 25,000 genes; however, only approximately 2% of these genes encode proteins, and many are noncoding genes that are transcribed to noncoding RNAs. Among these, those with a sequence length of >200 bp are called long noncoding RNAs (lncRNAs) [11]. Although these genes are not translated into proteins, they play crucial biological roles, such as regulating gene transcription, translation, and shearing processes and regulating microRNA and protein folding [12]. lncRNAs promote the occurrence and development of tumors by regulating tumor cell proliferation, migration, and invasion, modulating the cell cycle, and inhibiting apoptosis [13]. H19 is a classic cancer-promoting lncRNA that mediates the metastasis of ESCA in vitro and in vivo through the STAT3/EZH2/β-catenin axis. It is negatively regulated by let-7 at the posttranscriptional level, and its downregulation inhibits the proliferation, migration, and invasion of ESCA cells and promotes apoptosis [14]. UCA1 is another lncRNA that is overexpressed in gastrointestinal cancers and has an important regulatory role in cancer progression by acting as a competing endogenous RNA to regulate the expression of the target protein SOX4, thereby promoting ESCA cell proliferation [15]. UCA1 can also promote the glycolytic process in ESCA cells by sequestering miR-203 and alleviating its inhibitory effect on hexokinase 2 (HK2), thereby increasing HK2 levels and promoting the Warburg effect, cell proliferation, and metastasis [16]. Although most lncRNAs have tumor-promoting effects, some are tumor suppressors that inhibit tumor proliferation and migration and promote tumor cell apoptosis through multiple molecular mechanisms [13]. For example, the NEF overexpression reduces the expression of Wnt/β-catenin pathway-related proteins in ESCA cells, thereby inhibiting tumor cell proliferation, migration, and invasion [17]. Uc061hsf.1, a lncRNA, is a tumor suppressor gene and a direct transcriptional target of p53 [18]. Additionally, it regulates the expression of the downstream transcription factor FOXA1 and inhibits the proliferation and migration of ESCA cells [18].

Chen et al. reported that autophagy-related lncRNA prognostic markers, such as AL355574.1, are associated with the immune microenvironment and survival outcomes in patients with gastric cancers [19]. Increasing studies have found that lncRNAs affect the malignant progression of tumors by regulating autophagy [20-24]; however, the specific mechanism by which this regulation occurs remains unclear. Therefore, this study used coexpression analysis of autophagy-related mRNAs and lncRNAs in ESCA to identify autophagy-related lncRNAs, establish a prognostic model, and analyze its correlation with the clinical characteristics and survival of patients with ESCA.

2. Materials and Methods

2.1. Datasets and Sample Extraction

Transcriptome sequencing data and clinical-related data of patients with ESCA were downloaded from the Cancer Gene Atlas database [TCGA, Repository (http://cancer.gov/)]. The information obtained included the age, sex, survival time, survival status, tumor stage, and grade of patients with ESCA and other clinical data, as well as the transcriptome data of the tumors of 160 patients with ESCA and control samples comprising 11 adjacent tissues. The prcomp function of R was used to perform a principal component analysis (PCA).

2.2. Screening of Autophagy-Related Differentially Expressed Genes

The “limma” R package was used to sort and screen out differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNA) in ESCA using the following criteria: |log2 FC > 1| and P < 0.05. The list of autophagy genes was obtained from the Human Autophagy Database (HADb, http://autophagy.lu/clustering/index.htm). Autophagy-related DEmRNAs were obtained using the intersection of DEmRNAs and autophagy-associated genes listed on the website (http://bioinformatics.psb.ugent.be/webtools/Venn/). The correlation between lncRNAs and autophagy-related DEmRNAs was calculated using Pearson's correlation. lncRNAs with |r| > 0.3 and P < 0.001 were considered autophagy-related lncRNAs. The coexpression network was visualized using Cytoscape 3.7.2.

2.3. Identification of Prognostic Autophagy-Related lncRNAs

The prognostic value of autophagy-related lncRNAs was calculated using multiple Cox regression with P < 0.05. We constructed a risk score based on the linear combination of autophagy-related lncRNA expression levels multiplied by a regression coefficient. Based on this, participants were divided into high-risk and low-risk groups. Differences in survival between the two groups were compared using the log-rank test.

2.4. Development of the Prognostic Model

Univariate and multivariate Cox regression analyses were performed to explore the correlation between the autophagy-related lncRNA risk score and the age, sex, tumor grade, and stage of patients with ESCA. A nomogram was used to predict patient survival. The concordance index (C-index), calibration curves, and receiver operating characteristic (ROC) curves were used to determine the accuracy of the model.

2.5. Functional Analysis

Gene set enrichment analysis (GSEA) was conducted using the “clusterProfiler” package in R to determine the functional enrichment of autophagy-related lncRNAs. The enrichment plot package was used to visualize the results of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. StarBase (https://starbase.sysu.edu.cn/) and miRWalk (http://mirwalk.umm.uni-heidelberg.de/) website were used to predict the possible downstream targets of candidate lncRNAs.

2.6. Statistical Analysis

Statistical analyses were conducted using R language (version 3.6). Statistical tests were bilateral, and statistical significance was set at P < 0.05.

3. Results

3.1. Construction of a Coexpression Network

We identified 14142 recognizable lncRNAs, eight of which were downregulated and 53 were upregulated in ESCA (Figure 1(a)). Similarly, 638 DEmRNAs, including 472 with high expression and 166 with the low expression, were identified using TCGA dataset (Figure 1(b)). In total, 257 genes were obtained from the HADb, of which five (PINK1, BIRC5, ITPR1, PRKAB1, and GABARAPL1) were differentially expressed in ESCA simultaneously; therefore, we defined them as autophagy-related DEmRNAs. We used Pearson's correlation analysis to identify DElncRNAs that have a coexpression relationship with autophagy-related DEmRNAs (|r| > 0.3 and P < 0.001; Table 1). Finally, we identified five mRNAs (PINK1, BIRC5, ITPR1, PRKAB1, and GABARAPL1) and 11 lncRNAs (MAFG-DT, AC007637, AC091563, AC004982, FENDRR, AC037198, SNHG1, SEMA3B, ZNF710, and LINC02381) that might regulate autophagy in ESCA (Figures 1(c) and 1(d)).

Figure 1

(a) Differentially expressed lncRNAs in esophageal cancer. (b) Differentially expressed mRNAs in esophageal cancer. (c) Venn diagram showing the intersection of DEGs and HADb. (d) Venn diagram showing the intersection of DElncRNAs and ARlncRNAs (autophagy-related lncRNAs).

Table 1

Autophagy-related mRNAs (ARGmRNAs) and lncRNAs using correlation analysis.

ARGgene	lncRNA	cor	P value
BIRC5	MAFG-DT	0.46477	5.97E-10
ITPR1	AC007637.1	0.332063	1.79E-05
PRKAB1	AC007637.1	0.61085	9.73E-18
PRKAB1	AC091563.1	0.377207	8.80E-07
PRKAB1	AC004982.1	0.0318618	4.02E-05
ITPR1	FENDRR	0.790606	1.76E-35
GABARAPL1	FENDRR	0.308751	7.12E-05
ITPR1	AC037198.2	0.675193	1.23E-22
PINK1	AC092718.4	-0.33651	1.36E-05
PINK1	SNHG1	-0.35644	3.37E-06
BIRC5	SNHG1	0.341363	9.98E-06
PRKAB1	SEMA3B-AS1	0.312406	5.78E-05
ITPR1	ZNF710-AS1	0.413656	5.42E-08
PRKAB1	LINC02381	0.311697	6.02E-05

3.2. Identification of a Prognostic Autophagy-Related lncRNA Signature

According to the results of multivariate Cox regression analysis, four autophagy-related lncRNAs (AC092718, SEMA3B-AS1, FENDRR, and LINC02381) had prognostic value for patients with ESCA (P < 0.05). Of these, FENDRR, LINC02381, and SEMA3B-AS1 were prognostic risk factors, and AC092718 was a favorable prognostic factor (Figure 2). The four lncRNAs were used to establish an autophagy-related lncRNA signature. To further evaluate the prognostic value of the four-gene combination in patients with ESCA, the patients were divided into high- and low-risk groups according to the median score of the Cox regression model. The patients' risk scores continued to increase from left to right (Figure 3(b)). According to the survival diagram, patients in the high-risk group had shorter survival times and higher mortality rates than those in the low-risk group (Figure 3(c)). Kaplan–Meier survival analysis also confirmed that the survival time of patients in the high-risk group was significantly shorter than that of patients in the low-risk group (P < 0.05; Figure 3(a)). The expression of the four genes in the high- and low-risk groups were presented as a heat map in Figure 3(d). These results suggested that the combination of these four genes can be used as a specific prognostic index for patients with ESCA.

Figure 2

Survival analysis based on the four lncRNAs.

Figure 3

Multivariate Cox regression model analysis. (a) Overall survival curves for the four-gene combination in ESCA. (b) Risk scores: red indicates high risk, and green indicates low risk. (c) Survival diagram: red nodes indicate death, and green nodes indicate survival. (d) Heat map for the four-gene combination.

3.3. Construction of a Coexpression Network

Cytoscape was used to visualize the coexpression network, which included four lncRNAs and four mRNAs associated with autophagy in ESCA (Figure 4(a)). A hazard ratio (HR) > 1 was considered a risk factor, whereas HR < 1 was considered a protective factor. The Sankey diagram confirmed that FENDRR, LINC02381, and SEMA3B-AS1 were risk factors for patients with ESCA, whereas AC092718 was a protective factor (Figure 4(b)).

Figure 4

Coexpression network and Sankey diagram of the four prognostic autophagy-related lncRNAs. (a) Coexpression network of autophagy in Cytoscape; red indicates lncRNAs, and green indicates mRNAs. (b) Sankey diagram.

3.4. Clinical Value of the Autophagy-Related lncRNA Signature

The forest plot of the univariate analysis showed that stage, M, N, and risk scores could predict the prognosis of patients with ESCA (Figure 5(a)). Multivariate Cox analysis showed that only the risk score could be used as an independent prognostic factor (Figure 5(b)). However, ROC curve analysis showed that the area under the curve (AUC) values of the risk model, M, and N were 0.834, 0.544, and 0.646, respectively (Figure 5(d)). Therefore, compared to that with TNM staging, the risk model constructed herein was more accurate in predicting patient survival. As indicated in the nomogram, the risk score was the largest contributor to the 3-and 5-year overall survival rates of patients with ESCA (Figure 5(c)). The C-index of the prognostic model was 0.796 (95% CI: 0.739–0.853), and the 5-year survival rate AUC of the risk score was 0.834, indicating its reliable predictive ability (Figure 5(d)). To investigate the differences and similarities among the grouped samples, PCA was performed based on the expression profiles (Figure 6(a)) and autophagy-related lncRNA prognostic signatures (Figure 6(b)). Collectively, our results indicated that the autophagy-related lncRNA signature might be closely related to the progression and prognosis of ESCA.

Figure 5

Evaluation of prognostic models based on four autophagy-related lncRNAs. (a) Forest plots for (a) univariate and (b) multivariate Cox regression analyses in ESCA. (c) Nomogram of 3-year or 5-year overall survival based on the risk score, age, and TNM stage. (d) Receiver operating characteristic (ROC) curve analysis based on risk score and the clinicopathologic parameters.

Figure 6

PCA analysis based on the (a) expression profiles and (b) autophagy-related lncRNA prognostic signature.

3.5. Functional Analysis

To investigate the biological characteristics of the proposed lncRNAs, we analyzed GO enrichment analysis using GSEA. GO analysis showed that autophagy-related biological processes such as autophagosome assembly and autophagosome organization were significantly enriched in the ESCA group (Figure 7(a)). Subsequently, we identified 585 possible downstream targets of candidate lncRNAs. The KEGG enrichment analysis of target genes showed that autophagy and mitophagy were significantly enriched and identified pathways that also included pathways in cancer, gastric cancer, endocytosis, and p53 signaling pathways (Figure 7(b)).

Figure 7

Functional analysis. GO analysis using GSEA (a) and KEGG pathway (b) enrichment analysis.

4. Discussion

ESCA is one of the most malignant cancers worldwide, and China has the highest incidence of ESCA in the world [25]. The incidence of ESCA is associated with age, sex, obesity, eating habits, genetic susceptibility, and other risk factors. Its occurrence and development are multifactorial and complex processes [26]. Therefore, it is necessary to identify new biomarkers to improve its prognosis.

Autophagy-related pathological processes are increasingly being recognized as important disease mechanisms for nonmalignant (neurodegeneration and diffuse lung parenchymal disease) and malignant diseases [27]. Although autophagy has a dichotomous role in the regulation of cancer, increasing studies have shown the prosurvival role of autophagy in cancer progression and metastasis [28]. lncRNAs participate in the regulation of tumors and can be used as a molecular marker to predict the prognosis of patients with cancers [29]. Studies have shown that lncRNAs play a vital role in the development of ESCA [30-32]; however, there are no reports on autophagy-related lncRNA models that can predict the survival of patients with ESCA.

In this study, we constructed a coexpression network of lncRNAs and autophagy-related genes and used Cox regression analysis to identify four autophagy-related lncRNAs that might affect the prognosis of ESCA, specifically AC092718, FENDRR, LINC02381, and SEMA3B-AS1. We constructed an ESCA autophagy-related lncRNA risk assessment model and analyzed the correlation between the lncRNA model and the clinical characteristics and survival rate in patients with ESCA (Figure 8).

Figure 8

Schematic diagram of the experimental design.

AC092718.4 is a positive predictor of ovarian cancer and can play a role in cancer-related biological processes through a lncRNA–miRNA–mRNA regulatory network, such as cell cycle regulation, chromatin binding, modification, and remodeling the mTOR signaling pathway and the ovarian-specific BRCA1-related genome surveillance complex [33]. FENDRR, a lethal noncoding developmental regulatory RNA in the fetus, has been explored for its role in various cancers. In gastric cancer, FENDRR inhibits cell invasion and migration by downregulating fibronectin 1, and its low expression is associated with poor prognosis [34]. In breast cancer, FENDRR inhibits cell proliferation and is associated with a good prognosis [35]. In colorectal adenocarcinoma, the high FENDRR expression indicates a poor survival rate and promotes autophagy, apoptosis, and aging [36]. LINC02381 can be used as an oncogenic lncRNA and a tumor suppressor in cancer. It can aggravate the malignant phenotype and behavior of glioma cells by regulating the CBX5 expression. LINC02381 knockdown inhibits the malignant behavior of glioma cells and increases their proliferation [37]. In addition, LINC02381 is persistently expressed at low levels in breast cancer tissues, which might be related to immunity [38]. In gastric cancer tissues, the LINC02381 expression is downregulated, and its increased expression inhibits the activity of the Wnt pathway and cell cycle, increases cell apoptosis and caspase activity, and reduces the cell survival rate and proliferation rate of the human gastric cancer cell lines AGS and MKN45 [39]. SEMA3B-AS1 is a tumor suppressor lncRNA in cardia adenocarcinoma [40]. In hepatocellular carcinoma, the overexpression of SEMA3B-AS1 reduces hepatocellular carcinoma cell proliferation [41]. Additionally, SEMA3B-AS1 is significantly downregulated in ESCA, and its expression is related to TNM staging and lymph node metastasis. The SEMA3B-AS1 overexpression inhibits ESCA cell viability and invasiveness [42].

Our autophagy-related lncRNA signature could accurately predict the prognosis of patients with ESCA. This signature could be an independent indicator of ESCA, as revealed by univariate and multivariate Cox analyses. C-index and ROC curve analyses indicated that the model showed good discrimination and accuracy, suggesting that it might serve as a potential predictive method for patients with ESCA. GO analysis suggested that the selected lncRNAs were involved in the biological process of autophagosome assembly and autophagosome organization. Pathway enrichment analysis of downstream genes was also reported to be associated with autophagy and pathways in cancer. The results of functional enrichment analysis further confirmed that these prognostic autophagy-related lncRNAs play a role by regulating autophagy (Figure 8).

This study had some limitations. First, the data used in this study were based on TCGA and HADb databases which are limited; therefore, the analysis results might be biased. Second, we did not verify the expression of these four lncRNAs in ESCA. We plan to prospectively collect patients and conduct follow-up to construct a test cohort to further validate the accuracy of the model in the next stage. Third, functional experiments have not been performed, and the potential molecular mechanism associated with the predictive effect of autophagy-related lncRNAs remains unclear.

In conclusion, the prognostic model established in this study provided an effective basis to reveal the function of lncRNAs related to autophagy in ESCA. We identified four autophagy-related lncRNAs that were significantly related to the prognosis of patients with ESCA and can be used to distinguish among patients with different risk statuses. Therefore, these four autophagy-related lncRNAs and their markers could serve as molecular biomarkers and therapeutic targets for patients with ESCA.

5. Conclusions

This study established an autophagy-related coexpression network of mRNAs and lncRNAs in ESCA and identified four genes with potential applications for the diagnostic and prognostic analysis of patients with ESCA. These data could provide new insights into the diagnosis and treatment of ESCA. We aim to investigate the specific mechanisms by which these genes regulate autophagy in future studies.