Between July and December 2020, a total of 20 pairs of tumor and adjacent tissues were collected from patients with pathologically confirmed ECA at the Department of Thoracic Surgery, Fujian Cancer Hospital. The Ethics Committee of the Fujian Cancer Hospital approved the use of human tissues (Project Ethics Number: SQ 2020-063-01), and informed consent was acquired. The study protocol conformed to the principles enshrined in the Declaration of Helsinki.

Six pairs of esophageal cancer and adjacent normal tissues from Fujian Cancer Hospital between July 1st and July 17th, 2020 were used for RNA-sequencing. First, we characterized circRNA transcripts by sequencing analysis of ribosomal RNA and linear RNA. The total RNA was extracted by the Trizol method. Then, every sample was sequenced on Illumina HiSeq yielding an average of 42.38 million reads, which were mapped to the human reference genome (GRCh38/hg38) by TopHat2 (v2.1.1). The CIRC explorer program (v2.2.3) was used with the fusion junctions obtained from TopHat2 to identify both the circularizing junction and the spliced sequence of circRNAs. The whole step of library construction and sequencing was performed at Shanghai Life genes Technology Co. Ltd.

The chip data of circRNA expression profiles, that is, GSE131969 (GPL 19978) were downloaded from the GEO (Barrett et al., 2013). The miRNA and mRNA expression and clinical information in TCGA were obtained from the University of California Santa Cruz (USCS) Xena. The circRNA microarray data included three ECA tissues and three adjacent normal esophageal tissues. The miRNA sequencing data included 185 ECA samples and 13 adjacent samples, and the mRNA sequencing data included 162 ECA samples and 11 adjacent samples.

The Limma R package (version 4.0.2) (Ritchie et al., 2015) and the edge R package (version 3.14.0) (Robinson et al., 2010) were used for differential analyses. The expressions of all differentially expressed genes (DEGs) were visualized by volcano maps and two-way clustering heat maps. Principal component analysis (PCA) was conducted with the “pca3d” R package to explore gene expression patterns of ECA and normal groups.

MiRanda was used to predict the target miRNAs of the DEcircRNAs. Target mRNAs of DEmiRNAs were predicted using miRWalk (version 2.0) (Dweep, Gretz, 2015). The databases used for targeted mRNA prediction included miRWalk, Microt4, miRanda, mirbridge, miRDB, miRMap, and Targetscan. Overlapping mRNAs in six or more databases were considered the target mRNAs. Moreover, the ceRNA network was constructed and visualized by Cytoscape (version 3.6.1) (Shannon et al., 2003). The potential RNA binding proteins (RBPs) were analyzed by the CatRAPID database. A star rating score > 2 was considered indicative of strong binding effect. Cytoscape software was used to construct the circRNA–RBP binding network.

The mRNAs in the ceRNA network were analyzed using Search Tool for the Retrieval of Interacting Genes (STRING) (Szklarczyk et al., 2019), Database for Annotation, Visualization, and Integrated Discovery (DAVID) (Dennis et al., 2003), and Metascape (Zhou et al., 2019). Metascape, DAVID, and STRING databases were used to conduct analyses based on DisGeNET (version 7.0) (Piñero et al., 2020), Gene Ontology (GO) (Ashburner et al., 2000), biological processes (BPs), Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) database (version 2) (Han et al., 2018), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2017) analyses.

The enrichment scores of pathways were calculated using the GSVA (Hänzelmann et al., 2013), and a scoring matrix was obtained. Based on the matrix, the difference analysis was performed using Limma and visualized with a two-way clustering heat map.

Based on the optimal cutoff calculated by the survminer R package (Alboukadel et al., 2020), the patients were divided into low and high gene expression groups. Between-group differences were evaluated using the log-rank test with the survival (version 3.2-7) R package (Therneau, 2020), and target genes that showed significant association with survival were chosen (p < 0.05). Furthermore, single-factor and multi-factor Cox regression analyses were performed on mRNAs with significant results. Finally, LIM domain containing 2 (LIMD2) with significant results was selected as a key target.

Based on TCGA database, CIBERSORT (Newman et al., 2015) was used to perform immune infiltration–related analysis, and immune cells were visualized between groups. Based on the ratio of various immune cells, Spearman’s correlation analysis was performed to assess the correlation between immune cells and between targets and immune cells.

The autophagy-related genes were obtained from the Human Autophagy Database (HADb). The expressions of key target LIMD2 and autophagy-related genes were analyzed by Spearman’s correlation analysis, and the p-values were adjusted using the BH method.

GSEA (Reimand et al., 2019), performed by R package cluster Profiler (version 3.8.0) (Yu et al., 2012), was used to analyze the significant function and pathway difference between high- and low- LIMD2 groups. P-value <0.05 and FDR q-value < 0.25 were set as the cutoff criteria.

Among the clinical variables, age, M stage, N stage, T stage, sex, tumor stage, and LIMD2 expression were selected for univariate Cox regression analysis. Multivariate Cox regression analyses were conducted using variables that showed significant results in univariate Cox regression analysis.

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) was applied to verify the expression of the target LIMD2 in 20 pairs of ECA and adjacent normal esophageal tissues. The primers of LIMD2 were purchased from BioSune (Shanghai, China) (Table 2). An RT2162; All-in-One Mix with dsDNase (Monad Biotech Co. Ltd., Shanghai, China) was used to synthesize cDNA from 1 µg of total RNA. The qRT-PCR analyses were conducted on the Quant Studio 6 Flex qRT-PCR system (Applied Biosystems, Thermo Fisher Scientific Co. Ltd., United States) using the Hieff ®qPCR SYBR® Green Master Mix, Low Rox (Yeasen, Biotechnology Co., Ltd, Shanghai, China). The reaction was: 95°C for 10 min, then 40 cycles of 95°C for 15 s and 6°C for 1 min. The reference gene was GAPDH, and the relative gene expression levels were calculated using the 2−ΔΔCt method.

Statistical analyses were performed using R software (version 3.6.1), GraphPad Prism (version 8.0.1.244), SPSS (version 24.0), and the bioinformatics tools mentioned above. Differential expressions of genes were obtained by two-tailed Student’s t-test. The Benjamini and Hochberg FDR method was conducted to adjust the p-values. Enrichment analyses were analyzed using the hypergeometric test and Bonferroni correction. The outcomes were expressed as the mean ± SD and between-group differences were assessed using the paired Student’s t-test and Wilcoxon rank sum test. P values < 0.05 were considered indicative of statistical significance.

The differences between the clustering of the circRNAs (Figures 1A,D), mRNAs (Figures 2A,D), and miRNAs (Figure 3A) of the ECA and adjacent samples were analyzed. The results showed that most principal components could separate or tended to separate ECA samples from normal esophageal tissues.

CircRNAs in GSE131969 and RNA-seq data (p- value <0.05 and |log2 FC| > 1.5), mRNAs in RNA-seq data (p- value <0.05 and |log2 FC| > 1) and TCGA dataset (p value <0.05 and |log2 FC| > 1.5), as well as miRNAs inTCGA (p value <0.05 and |log2 FC| > 1.5) were analyzed. A total of 368 DEcircRNAs were identified in GSE131969 (Figure 1B), 493 DEcircRNAs were identified in RNA-seq data (Figure 1E), 1576 DEmRNAs were identified in RNA-seq data (Figure 2B), 2450 DEmRNAs were identified in TCGA (Figure 2E), and 86 DEmRNAs were identified in TCGA (Figure 3B). The DEcircRNAs in the GEO (Figure 1C) and RNA-seq data (Figure 1F) and DEmRNAs in RNA-seq data (Figure 2C) and TCGA (Figure 2F), as well as DEmiRNAs in TCGA (Figure 3C), were visualized using the hierarchical cluster maps. Seven upregulated circRNAs (Figure 1G), zero downregulated circRNAs (Figure 1H) and 61 upregulated mRNAs (Figure 2G) and six downregulated mRNAs (Figure 2H) were significantly differentially expressed in both datasets. Finally, seven DEcircRNAs in the GSE131969 (Supplementary Figure S1A) and RNA-seq data (Supplementary Figure S1B) were visualized.

The miRNAs targeted by the seven DEcircRNAs were predicted and 22 circRNA–miRNA pairs were identified. The mRNAs targeting the miRNAs were identified, and 34 mRNAs were obtained through the intersection with DEmRNAs (Figure 4D). Subsequently, the differential expression of the 34 mRNAs TCGA database (Supplementary Figure S3A) and RNA-seq data were visualized (Supplementary Figure S3B). A ceRNA network was established based on seven circRNAs, 22 miRNAs, and 34 mRNAs (Figure 3D). The circRNA–RBP network was constructed, including seven circRNAs and 37 RBPs (Figure 3E).

The enriched terms analyzed based on the DisGeNET database are shown in Figure 4A, including adult T-Cell lymphoma and childhood Langerhans cell histiocytosis. The enriched GO terms are shown in Figure 4B, including positive regulation of cell adhesion and positive regulation of protein kinase B signaling. The enriched terms analyzed based on the TRRUST database were SP1, RELA, and NFKB1 (Figure 4C). Moreover, GSVA showed enrichment of ether lipid metabolism or acute myeloid leukemia in the cancer group (Figure 4E). Based on the DAVID database, the outcomes showed that mRNAs were enriched in positive regulation of T cell activation and cell migration (BP), extracellular space, plasma membrane (cellular components, CCs), DNA binding (Molecular Functions, MFs), and NF-kappa B signaling pathway (KEGG, Supplementary Figure S2A). Based on the STRING database, the results showed that targets were enriched in positive regulation of lymphocyte activation and immune system process (BPs), extracellular space (CC), protein binding and DNA binding (MFs), NF-kappa B pathway and ether lipid metabolism (KEGG), and regulation of PI3K/AKT signaling (record control the memory, RCTM, Supplementary Figure S2B).

Targets that showed a significant association with survival (p < 0.05) included LIMD2, ARL11, and SERPINB7 (Figures 5A–L). LIMD2 with significant regression results was selected as a key target.

The composition of immune cells in the tumor microenvironment (TME) in the ECA and normal groups was analyzed. We observed significant differences in TME cell infiltration and composition in the two groups, including memory activated CD4 T cells and M1 macrophages (Figure 6A). The R ² of TME cells was also calculated, and the results showed that memory activated CD4 T cells had statistical significance (Figure 6B). A positive relationship was observed between LIMD2 and memory activated CD4 T cells (Figure 6C). These analyses indicated a positive correlation of the ECA group with immune-relevant signatures.

Several studies have indicated the importance of autophagy in the development of ECA (Chen et al., 2020; Li et al., 2020). We investigated the relationship between LIMD2 and autophagy-related genes (Figures 7A–I). The top nine autophagy-related genes that showed a positive relationship with LIMD2 were identified, including BCL-2, CCL2, RGS19, DRAM1, GRID1, NLRC4, PELP1, PRKCQ, and RGS19.

The results demonstrated that the ECA patients with high expression of LIMD2 showed enrichment of the following pathways: systemic lupus erythematosus, DNA replication, and T cell receptor signaling pathway (Supplementary Figures S4A–I).

Univariate and multivariate Cox regression analyses were performed to explore independent predictors of survival in ECA patients. A total of 162 ECA patients were included in the regression analysis. In univariate Cox regression analysis, LIMD2, M stage, N stage, and tumor stage showed a strong association with OS. Based on the results of univariate Cox analysis, multivariate Cox regression analysis was performed to analyze the effect of LIMD2 and other related clinical phenotypes on the prognosis of ECA. In multivariate Cox regression analysis, LIMD2 showed a significant correlation with OS (Table 1, Supplementary Figure S5). The results showed that LIMD2 is an independent predictor of OS in ECA patients.

To verify the RNA-seq and bioinformatics results, the target LIMD2 was chosen for validation by quantitative reverse-transcription polymerase chain reaction (qRT-PCR) with primers (Table 2). The results showed significantly higher expression of LIMD2 in ECA tissues compared with that of normal tissues (Figure 8).