Yuanmei Chen1†Xinyi Huang2†Kunshou Zhu1†Changkun Li3Haiyan Peng3Lin Chen2Zhengrong Huang4Yangfan Zhang3Guibin Weng1Tianya Xiao1
Junqiang Chen2*
Yuanji Xu2*- 1Department of Thoracic Surgery, Fujian Medical University Cancer Hospital, Fujian Cancer Hospital, Fuzhou, China
- 2Department of Radiation Oncology, Fujian Medical University Cancer Hospital, Fujian Cancer Hospital, Fuzhou, China
- 3Fujian Key Laboratory of Innate Immune Biology, Biomedical Research Center of South China, Fujian Normal University Qishan Campus, Fuzhou, China
- 4Department of Integrative Traditional Chinese and Western Medicine, Fujian Medical University Cancer Hospital, Fujian Cancer Hospital, Fuzhou, China
Globally, esophageal cancer (ECA) is the seventh most common cancer and sixth most common cause of cancer-associated mortality. However, there are no reliable prognostic and predictive molecular markers for ECA; in addition, the pathogenesis of ECA is not fully elucidated. The expressions of circular RNAs (circRNAs), microRNAs (miRNAs), and messenger RNAs (mRNAs) of ECA and control groups were obtained from the RNA-sequencing (RNA-seq) data of our hospital, the Gene Expression Omnibus (GEO), and The Cancer Genome Atlas (TCGA) datasets. Analyses of differentially expressed genes, the circRNA–miRNA–mRNA–competing endogenous RNA (ceRNA) network, and functional/pathway enrichment were conducted. The key targets in the ceRNA network that showed significant results in survival Cox regression analyses were selected. Furthermore, analyses of immune infiltration and autophagy genes related to the key targets were performed. Seven circRNAs, 22 miRNAs, and 34 mRNAs were identified as vital genes in ECA; the nuclear factor-κ-gene binding (NF-κB) and phosphatidylinositol-3 kinase/protein kinase B (PI3K-Akt) signaling were identified as the most enriched pathways. In addition, the LIM domain containing 2 (LIMD2) was an independent predictor of prognosis in ECA patients and closely associated with immunity and autophagy. Moreover, quantitative reverse-transcription polymerase chain reaction (qRT-PCR) revealed significant upregulation of LIMD2 expression in ECA tissues. ECA may be closely correlated with NF-κB and PI3K/Akt signaling. In addition, LIMD2 could be a potential prognostic and predictive marker of ECA.
Introduction
Esophageal cancer (ECA), including esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC), is the seventh most common malignant tumor in the world and the sixth most common primary cause of tumor-associated death, especially in Asia (Bray et al., 2018). In 2018, approximately 572,034 new ECA cases were diagnosed (accounting for 3.2% of all cancer-related deaths) and 508,585 deaths were caused by ECA (accounting for 3.3% of all malignancies) (Bray et al., 2018). Surgical removal of tumors is the conventional treatment for most types of ECA. Although advances in surgical procedures, chemoradiotherapy, and immunotherapy have helped improve the overall survival (OS) of patients, the current 5-year survival rates of ECA patients still range from 15 to 25% (Pennathur et al., 2013). Therefore, in-depth exploration of potential prognostic and predictive markers, therapeutic targets, and mechanisms of ECA is a key imperative to improve the treatment outcomes and prognosis of ECA patients.
Due to the rapid advances in microarray and RNA-sequencing (RNA-seq) technology, multiple circular RNAs (circRNAs), microRNAs (miRNAs), and messenger RNAs (mRNAs) have been identified as important genes in ECA (Ritchie et al., 2015). CircRNAs are endogenous non-coding RNAs with closed-loop structures without 5′caps and 3′ tails (Kristensen et al., 2019). Recent studies have found that the mRNAs targeted by miRNAs in the circRNA–miRNA–mRNA–competing endogenous RNA (ceRNA) regulatory network could serve as a key therapeutic target in cancer. In addition, circRNAs were shown to play a vital role in the occurrence and development of cancers by positively or negatively regulating the circRNA–miRNA–mRNA–ceRNA axis or by acting as protein “sponges” in cancer cells. CircRNAs in the ceRNA network may indirectly regulate those target mRNAs by serving as miRNA “sponges” (Zhang et al., 2017; Wu et al., 2019). In a study, low expression of large tumor suppressor kinase 1 (LATS1) was found to be related to the tumor stage and poor prognosis of gastric cancer patients; in addition, circLARP4 was found to inhibit the malignant biological behavior of gastric cancer as a tumor suppressor via the modulation of the circLARP4/miR-424-5p/LATS1 axis (Yang et al., 2018).
Genome-wide analyses of ECA have indicated complex mutation situations and discovered significant gene mutations (including TP53, MLL2, NOTCH1, and PTEN), repeated copy number amplifications in SOX2, TERT, and FGFR1, as well as frequent deletion of RB1 (Chang et al., 2017; (Cancer Genome Atlas Research Network, 2017). However, previous ECA genome studies had some limitations, including relatively small sample size, low-coverage whole-genome sequencing, and one-sided analysis based only on chip data or sequencing data. Therefore, more comprehensive studies are required to overcome the shortcomings of previous studies.
The purpose of our study was to explore the potential prognostic and predictive markers and to unravel the mechanisms of ECA from our RNA-seq data, the Gene Expression Omnibus (GEO), and The Cancer Genome Atlas (TCGA) database. Differentially expressed circRNAs (DEcircRNAs), differentially expressed miRNAs (DEmiRNAs), and differentially expressed mRNAs (DEmRNAs) were identified using the Limma R package, and the circRNA–miRNA–mRNA–ceRNA regulatory network was constructed based on these differentially expressed genes. In addition, the enriched pathways were investigated by performing multiple functional and pathway enrichment analyses. Furthermore, the key target mRNAs in the ceRNA regulatory network were determined using survival Cox regression analyses. Finally, analyses of immune infiltration autophagy genes related to the key targets were performed. Our study may provide novel insights into the prognosis and treatment of ECA based on the pathogenetic mechanism.
Methods
Sample Collection
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.
RNA Sequencing
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.
Datasets
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.
Differential Expression Analyses
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.
Construction of the ceRNA and circRNA-RBP Network
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.
Functional and Pathway Enrichment Analyses
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.
Gene Set Variation Analysis
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.
Survival Analysis
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.
Immune Infiltration Related Analysis
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.
Autophagy Correlation Analysis
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.
Gene Set Enrichment Analysis
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.
Single-Factor and Multi-Factor Cox Regression Analyses
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.
qRT-PCR Validation for the Expression of LIMD2
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 Analysis
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.
Results
Data Collection and Preprocessing
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.
FIGURE 1. Identification of differentially expressed circRNAs (DEcircRNAs) in esophageal cancer (ECA). (A, D) Principal component analysis (PCA) of circRNA expression between ECA and normal groups in GSE131969 and the present RNA-sequencing (RNA-seq) data. (B, E) Volcano plot of the distributions of all DEcircRNAs in GSE131969 and the current RNA-seq data. (C, F) Hierarchical clustering heatmap of dysregulated circRNAs between ECA and the adjacent normal tissues in GSE131969 and the present RNA-seq data. (G, H) Comparison of differentially expressed circRNAs identified in GSE131969 and our RNA-seq data.
FIGURE 2. Identification of differentially expressed mRNAs (DEmRNAs) in ECA. (A, D) PCA of mRNA expression between ECA and normal groups in TCGA cohort and the present RNA-seq data. (B, E) Volcano plot of the distributions of all differentially expressed mRNAs in TCGA cohort and the current RNA-seq data. (C, F) Hierarchical clustering heatmap of dysregulated mRNAs between ECA and the adjacent normal tissues in TCGA cohort and the present RNA-seq data. (G, H) Comparison of differentially expressed mRNAs identified in TCGA cohort and our RNA-seq data.
FIGURE 3. Identification of differentially expressed miRNAs (DEmiRNAs) in ECA, the ceRNA network, and the circRNA–RNA binding protein (RBP) network. (A) PCA of miRNA expression between ECA and normal groups in TCGA cohort. (B) Volcano plot of the distributions of all DEmiRNAs in TCGA cohort. (C) Hierarchical clustering heatmap of dysregulated miRNAs between ECA and the adjacent normal tissues in TCGA cohort. (D) CeRNA network of circRNA–miRNA–mRNA in ECA; (E) circRNA–RBP binding network.
Identification of DEcircRNAs, DEmiRNAs, and DEmRNAs
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.
Construction of the CeRNA Network and CircRNA–RBP Network
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).
FIGURE 4. Biological enrichment analyses of target mRNAs. (A) Enrichment results based on the DisGeNET database. (B) Enrichment results based on the gene ontology (GO) and biological process (BP) databases. (C) Enrichment results based on the TRRUST database. (D) Comparison of target mRNAs predicted by miRNAs and differentially expressed mRNAs. (E) GSVA analysis of key targets.
Functional and Pathway Enrichment Analyses
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).
Survival Analysis
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.
FIGURE 5. Effect of target mRNA expression on the OS of ECA patients in TCGA database. Kaplan–Meier (KM) curves demonstrating the distinct outcomes of ECA patients. (A–L) KM curve of LIMD2 (A), ARL11 (B), CARD11 (C), ENKUR (D), HS3ST3A1 (E), KIT (F), LCK (G), LPCAT1 (H), NABP1 (I), NRXN1 (J), PPP1R18 (K), and SERPINB7 (L).
Immune Infiltration Analysis
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 R2 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.
FIGURE 6. Immune infiltration correlation analysis in ECA based on TCGA cohort. (A) Fraction of the immune tumor microenvironment (TME) cells in the ECA and normal groups, and TME cells between the ECA and normal groups are analyzed in the figure given below. Within each group, the scattered dots represent the expression values of TME cells. The bottom and top of the boxes are the 25th and 75th percentiles (interquartile range). The whiskers encompass 1.5 times the interquartile range. (B) Immune infiltration correlation analysis; (C) Plots showing Spearman’s correlation of LIMD2 and memory activated CD4 T cells.
Analysis of LIMD2 and Autophagy-Related Genes
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.
FIGURE 7. Plots showing Spearman’s correlation of LIMD2 and autophagy-related genes in TCGA cohort. In every plot, the X-axis represents LIMD2 and the Y-axis represents autophagy-related genes. (A–I) Spearman’s correlation of LIMD2 and BCL-2 (A), CCL2 (B), CXCR4 (C), DRAM1 (D), GRID1 (E), NLRC4 (F), PELP1 (G), PRKCQ (H), and RGS19 (I). Dashed line in each plot is the regression line.
GSEA
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).
Cox Regression Analysis
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.
TABLE 1. Results of univariate analysis and multivariate analysis showing prognostic factors for OS.
Validation of LIMD2 by QRT-PCR
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).
TABLE 2. Primers for LIMD2 and GAPDH.
FIGURE 8. Validation of the expression of target LIMD2. The LIMD2 expression in ECA (n = 20) and adjacent normal tissues (n = 20) was evaluated by qRT-PCR; the results were analyzed using the paired sample t-test (A) and Wilcoxon rank sum test (B). Results expressed as mean ± standard deviation (SD). *p < 0.01 and **p < 0.001.
Discussion
In this study, seven circRNAs, 22 miRNAs, and 34 target mRNAs in the circRNA–miRNA–mRNA–ceRNA regulatory network were identified as crucial genes in ECA; in addition, the nuclear factor-k-gene binding (NF-κB) and phosphatidylinositol-3 kinase/protein kinase B (PI3K-Akt) signaling were identified as the most enriched pathways. Survival Cox regression analyses indicated that the target LIMD2 in the ceRNA network may act as an independent predictor of OS in ECA patients. On further analyses of immune infiltration and autophagy genes related to target LIMD2, LIMD2 showed a close linkage with immunity and autophagy in ECA. Furthermore, we verified the high expression level of LIMD2 by qRT-PCR based on 20 pairs of ECA and normal samples. Hence, LIMD2 is a potential molecular marker of prognostic and predictive significance in ECA.
The ceRNA networks may play a vital role in the development of cancer (Liu et al., 2021). Seven circRNAs, 22 miRNAs, and 34 mRNAs were identified in the ceRNA network. The increase of has_circ_0000154 (circDCAF6) identified in our study was related to tumor invasion, positive lymph node metastasis, and a higher TNM stage in GC patients. It could serve as an independent prognostic indicator (Wu et al., 2019). A total of 12 identified mRNAs (HMGA2, CCL22, CDKN2D, CEACAM6, DKK2, KIT, KRT14, LEF1, LPCAT1, NCCRP1, NRXN1, and PPP1R18) were found to be related to ECA. In a previous study, high-mobility group AT-hook 2 (HMGA2) was shown to regulate transcription by inducing structural alterations in the chromatin (Mansoori et al., 2021). Several studies have shown that HMGA2 is re-expressed in most tumors and plays a vital role in tumorigenesis (Lu et al., 2021). The stability of HMGA2 may be regulated by hepatitis B X-interacting protein (HBXIP) via the Akt-PCAF pathway, thereby promoting the growth of ECA cells (Wu et al., 2020). All 22 DEmiRNAs in the ceRNA network have been reported to be related to cancer. MiR-141-3p plays an important role in various carcinomas (Huang et al., 2017), and it was found to be highly upregulated in ECA cells (Phatak et al., 2021). A recent study reported that miR-141-3p may inhibit the expression of pleckstrin homology domain leucine-rich repeat protein phosphatase-2 (PHLPP2), a negative regulator of the PI3K/AKT signaling, and could serve as an biomarker in ECA (Ishibashi et al., 2018).
CircRNAs with RBP binding sites can act as sponges for RBPs and may indirectly modulate their functions (Zhang et al., 2017). In this study, the six identified RBPs (SRSF1, SRSF2, PCBP2, TIA1, FUS, and FMR1) were closely related to ECA. PCBP2 performs multiple functions, such as stabilization of mRNAs and silencing or promotion of translation (Ye et al., 2016). Several studies have shown that PCBP2 may promote tumor growth. Ye et al. found that PCBP2 regulates the proliferation and apoptosis of ESCC cells and may serve as a novel therapeutic target in ESCC (Ye et al., 2016). Our findings are consistent with those of the previous study.
In our study, multiple diseases were enriched in ECA, including adult T cell lymphoma (ATL), childhood Langerhans cell histiocytosis, and classical Hodgkin’s lymphoma that are strongly associated with immune cells. ATL is a T cell lymphoproliferative tumor of mature CD4+ CD25 + T cells (Mehta-Shah et al., 2017), and Langerhans cell histiocytosis is characterized by the accumulation of Langerhans cells and antigen-presenting cells (Rayamajhi et al., 2020). Our results also revealed a close relation of many enriched GO and KEGG terms with immune response, including positive regulation of T cell activation and immune system process. GSEA showed that immunity was strongly associated with ECA, including systemic lupus erythematosus and the T cell receptor signaling pathway. The NF-κB signaling and PI3K/AKT signaling pathways were the most enriched pathways. Zheng et al. depicted the entire immune landscape, including the innate and acquired immune cell map, in ESCC and adjacent tissues; their work revealed that ESCC is enriched in immune-suppressive cell mass (Zheng et al., 2020). Tong et al. revealed that 14-3-3ζ may enhance the invasion and growth of ESCC cells by inhibiting the S1PR2 protein expression via the NF-κB pathway (Tong et al., 2018). GSVA indicated enrichment of ether lipid metabolism in ECA. Disorders of ether lipid metabolism are vital signs of tumors, which serve as the basis of tumor pathogenicity (Piano et al., 2015). Cao et al. found that lncRNAs may interact with their adjacent coding RNAs to modulate ether lipid metabolism (Cao et al., 2013).
Furthermore, we found that 12 genes were strongly related to OS. We included these 12 mRNAs in Cox regression analyses and identified LIMD2 as a vital prognostic factor. In our study, high expression of LIMD2 was associated with better OS in ECA. However, Zhang et al. found that high expression of LIMD2 may enhance the progression of non–small cell lung carcinoma (NSCLC); in addition, the overexpression of LIMD2 was closely related to lymph node metastasis, distant metastasis, and advanced stage. The survival time of NSCLC patients with the overexpression of LIMD2 was shorter than that of patients with lower expression, indicating that LIMD2 may serve as a therapeutic target in NSCLC (Zhang et al., 2019). Our study indicated that LIMD2 may serve as an independent predictor of prognosis in ECA patients. In the TME cell composition analysis, Yu et al. found that the high-immunity group of cutaneous melanoma specimens had the highest level of memory activated CD4 T cells, and the OS rate was poor (Yu et al., 2020). A recent research unraveled the connection between clinical information and immune signatures in GC. They noticed that the high-risk group showed greater proportion of memory activated CD4 T cells and M1 macrophages (Liu et al., 2020). Our study has also demonstrated positive relationship of LIMD2 with nine autophagy-related genes. In our study, BCL-2, an autophagy-related gene showed a positive relationship with LIMD2, and BCL-2–associated athanogene 3 (BAG3) was shown to be involved in multiple biological processes, including cell proliferation, cell vitality, and apoptosis (Rosati et al., 2011). Compared with Barrett’s metaplasia (BE) and normal samples, the amount of T lymphocytes with downregulation of BCL-2 was notably increased in EAC (Berndt et al., 2010). In short, LIMD2 is a potential marker of prognostic and predictive significance in ECA.
Nevertheless, some limitations of our study should be considered while interpreting the results. First, owing to the use of several datasets in our analyses, the effect of inter-batch differences on our results cannot be ruled out. In addition, further studies with a larger sample size are required to verify the high expression of LIMD2 in ECA. Moreover, the sample size of squamous carcinoma and adenocarcinoma was too small to be analyzed separately. To avoid unreliable results, we analyzed both histological subtypes. We will further expand the sample size to explore the differences between different subtypes. Finally, further in vitro and in vivo biological experiments on immunity and autophagy are required for in-depth characterization of the functions of the identified targets and potential mechanisms in ECA.
In conclusion, we identified seven circRNAs, 22 miRNAs, and 34 mRNAs in the ceRNA network in ECA. The NF-κB and PI3K/Akt signaling pathways were the most enriched pathways in ECA. Furthermore, LIMD2 was notably upregulated in ECA tissue samples and may serve as a potential prognostic and predictive marker that is closely associated with immunity and autophagy. Our study provides insights into the pathogenesis, prognosis, and therapeutic strategy for ECA.
Data Availability Statement
The datasets generated and analyzed for this study can be found in the National Center for Biotechnology Information Gene Expression Omnibus (Barrett et al., 2013) repository under the accession number GSE131969.
Ethics Statement
The studies involving human participants were reviewed and approved by the Ethics Committee of the Fujian Cancer Hospital approved the use of human tissues (Project Ethics Number: K2021-027-01). The patients/participants provided their written informed consent to participate in this study.
Author Contributions
YC and YX conceived and designed the whole project. YC and XH drafted the manuscript. XH and KZ analyzed the data and prepared figures. XH and CL performed the analyses. HP and LC computed the tables. ZH and YZ collected the data. GW, TX, and JC provided specialized expertise. YX designed the whole project and revised the manuscript. All authors read and approved the final manuscript.
Funding
This study was supported in part by grants from Fujian Provincial Health Technology Project, China (#2020GGA017), the Natural Science Foundation of Fujian Province, China (#2021J01428, #2020J011123), and China Association of Gerontology and Geriatrics (#2021-CSGOR-1), Fujian Provincial Clinical Research Center for Cancer Radiotherapy and Immunotherapy (2020Y2012), and the National Clinical Key Specialty Construction Program.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We would like to thank Medjaden Bioscience Limited (MJD) for the paper revision support.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.774432/full#supplementary-material
References
Alboukadel, K., Marcin, K., Przemyslaw, B., and Scheipl, F. (2020). Survminer: Drawing Survival Curves Using 'ggplot2'. R package version 0.4.8 Available at:https://CRAN.R-project.org/package=survminer/ (Accessed March 20, 2021).
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene Ontology: Tool for the Unification of Biology. Nat. Genet. 25, 25–29. doi:10.1038/75556
Barrett, T., Wilhite, S. E., Ledoux, P., Evangelista, C., Kim, I. F., Tomashevsky, M., et al. (2013). NCBI GEO: Archive for Functional Genomics Data Sets-Update. Nucleic Acids Res. 41, D991–D995. doi:10.1093/nar/gks1193
Berndt, U., Philipsen, L., Bartsch, S., Hu, Y., Röcken, C., Bertram, W., et al. (2010). Comparative Multi-Epitope-Ligand-Cartography Reveals Essential Immunological Alterations in Barrett's Metaplasia and Esophageal Adenocarcinoma. Mol. Cancer 9, 177. doi:10.1186/1476-4598-9-177
Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: A Cancer J. Clinicians 68, 394–424. doi:10.3322/caac.21492
Cao, W., Wu, W., Shi, F., Chen, X., Wu, L., Yang, K., et al. (2013). Integrated Analysis of Long Noncoding RNA and Coding RNA Expression in Esophageal Squamous Cell Carcinoma. Int. J. Genomics 2013, 1–10. doi:10.1155/2013/480534
Cancer Genome Atlas Research Network (2017). Integrated Genomic Characterization of Oesophageal Carcinoma. Nature 541, 169–175. doi:10.1038/nature20805
Chang, J., Tan, W., Ling, Z., Xi, R., Shao, M., Chen, M., et al. (2017). Genomic Analysis of Oesophageal Squamous-Cell Carcinoma Identifies Alcohol Drinking-Related Mutation Signature and Genomic Alterations. Nat. Commun. 8, 15290. doi:10.1038/ncomms15290
Chen, C., Chen, S., Cao, H., Wang, J., Wen, T., Hu, X., et al. (2020). Prognostic Significance of Autophagy-Related Genes within Esophageal Carcinoma. BMC Cancer 20, 797. doi:10.1186/s12885-020-07303-4
Dennis, G., Sherman, B. T., Hosack, D. A., Yang, J., Gao, W., Lane, H. C., et al. (2003). DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome Biol. 4, P3. doi:10.1186/gb-2003-4-5-p3
Dweep, H., and Gretz, N. (2015). MiRWalk2.0: a Comprehensive Atlas of microRNA-Target Interactions. Nat. Methods 12, 697. doi:10.1038/nmeth.3485
Han, H., Cho, J.-W., Lee, S., Yun, A., Kim, H., Bae, D., et al. (2018). TRRUST V2: an Expanded Reference Database of Human and Mouse Transcriptional Regulatory Interactions. Nucleic Acids Res. 46, D380–D386. doi:10.1093/nar/gkx1013
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7
Huang, S., Wa, Q., Pan, J., Peng, X., Ren, D., Huang, Y., et al. (2017). Downregulation of miR-141-3p Promotes Bone Metastasis via Activating NF-Κb Signaling in Prostate Cancer. J. Exp. Clin. Cancer Res. 36, 173. doi:10.1186/s13046-017-0645-7
Ishibashi, O., Akagi, I., Ogawa, Y., and Inui, T. (2018). MiR-141-3p Is Upregulated in Esophageal Squamous Cell Carcinoma and Targets Pleckstrin Homology Domain Leucine-Rich Repeat Protein Phosphatase-2, a Negative Regulator of the PI3K/AKT Pathway. Biochem. Biophysical Res. Commun. 501, 507–513. doi:10.1016/j.bbrc.2018.05.025
Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., and Morishima, K. (2017). KEGG: New Perspectives on Genomes, Pathways, Diseases and Drugs. Nucleic Acids Res. 45, D353–D361. doi:10.1093/nar/gkw1092
Kristensen, L. S., Andersen, M. S., Stagsted, L. V. W., Ebbesen, K. K., Hansen, T. B., and Kjems, J. (2019). The Biogenesis, Biology and Characterization of Circular RNAs. Nat. Rev. Genet. 20, 675–691. doi:10.1038/s41576-019-0158-7
Li, Z., Dong, K., Guo, P., Tan, Z., Zhang, F., Tian, Y., et al. (2020). Identification of Autophagy-Related Genes and Small-Molecule Drugs in Esophageal Carcinoma. Med. Sci. Monit. 26, e921855. doi:10.12659/MSM.921855
Liu, C., Chen, B., Huang, Z., Hu, C., Jiang, L., and Zhao, C. (2020). Comprehensive Analysis of a 14 Immune-Related Gene Pair Signature to Predict the Prognosis and Immune Features of Gastric Cancer. Int. Immunopharmacology 89, 107074. doi:10.1016/j.intimp.2020.107074
Liu, G., Sun, J., Yang, Z.-F., Zhou, C., Zhou, P.-Y., Guan, R.-Y., et al. (2021). Cancer-associated Fibroblast-Derived CXCL11 Modulates Hepatocellular Carcinoma Cell Migration and Tumor Metastasis through the circUBAP2/miR-4756/IFIT1/3 axis. Cell Death Dis 12, 260. doi:10.1038/s41419-021-03545-7
Lu, W., Huang, Z., Wang, J., and Liu, H. (2021). Long Non-coding RNA DANCR Accelerates Colorectal Cancer Progression via Regulating the miR-185-5p/HMGA2 axis. J. Biochem. mvab011, undefined, undefined. doi:10.1093/jb/mvab011
Mansoori, B., Mohammadi, A., Ditzel, H. J., Duijf, P. H. G., Khaze, V., Gjerstorff, M. F., et al. (2021). HMGA2 as a Critical Regulator in Cancer Development. Genes 12, 269. undefined. doi:10.3390/genes12020269
Mehta-Shah, N., Ratner, L., and Horwitz, S. M. (2017). Adult T-Cell Leukemia/Lymphoma. Jop 13, 487–492. doi:10.1200/JOP.2017.021907
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12, 453–457. doi:10.1038/nmeth.3337
Pennathur, A., Gibson, M. K., Jobe, B. A., and Luketich, J. D. (2013). Oesophageal Carcinoma. The Lancet 381, 400–412. doi:10.1016/S0140-6736(12)60643-6
Phatak, P., Noe, M., Asrani, K., Chesnick, I. E., Greenwald, B. D., and Donahue, J. M. (2021). MicroRNA‐141‐3p Regulates Cellular Proliferation, Migration, and Invasion in Esophageal Cancer by Targeting Tuberous Sclerosis Complex 1. Mol. Carcinogenesis 60, 125–137. doi:10.1002/mc.23274
Piano, V., Benjamin, D. I., Valente, S., Nenci, S., Marrocco, B., Mai, A., et al. (2015). Discovery of Inhibitors for the Ether Lipid-Generating Enzyme AGPS as Anti-cancer Agents. ACS Chem. Biol. 10, 2589–2597. doi:10.1021/acschembio.5b00466
Piñero, J., Ramírez-Anguita, J. M., Saüch-Pitarch, J., Ronzano, F., Centeno, E., Sanz, F., et al. (2020). The DisGeNET Knowledge Platform for Disease Genomics: 2019 Update. Nucleic Acids Res. 48, D845–D855. doi:10.1093/nar/gkz1021
Rayamajhi, S., Shahi, R., Maharjan, S., Sharma, S., and Sudhir Suman, K. (2020). Muscle Infiltrative Adult Multisystem Langerhans Cell Histiocytosis Detected on Fluorodeoxyglucose Positron Emission Tomography/computed Tomography - A Rare Case. Indian J. Nucl. Med. 35, 342–344. doi:10.4103/ijnm.IJNM_88_20
Reimand, J., Isserlin, R., Voisin, V., Kucera, M., Tannus-Lopes, C., Rostamianfar, A., et al. (2019). Pathway Enrichment Analysis and Visualization of Omics Data Using g:Profiler, GSEA, Cytoscape and EnrichmentMap. Nat. Protoc. 14, 482–517. doi:10.1038/s41596-018-0103-9
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). EdgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26, 139–140. doi:10.1093/bioinformatics/btp616
Rosati, A., Graziano, V., De Laurenzi, V., Pascale, M., and Turco, M. C. (2011). BAG3: a Multifaceted Protein that Regulates Major Cell Pathways. Cel Death Dis 2, e141. doi:10.1038/cddis.2011.24
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13, 2498–2504. doi:10.1101/gr.1239303
Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING V11: Protein-Protein Association Networks with Increased Coverage, Supporting Functional Discovery in Genome-wide Experimental Datasets. Nucleic Acids Res. 47, D607–D613. doi:10.1093/nar/gky1131
Therneau, T. (2020). A Package for Survival Analysis in R. R Package Version 3.2-7. Available at: https://CRAN.R-project.org/package=survival/ (Accessed March 20, 2021).
Tong, S., Chen, S.-C., Xu, K.-Y., Fang, B., Wang, S.-H., and Wang, J.-J. (2018). 14-3-3ζ Promotes Esophageal Squamous Cell Carcinoma Invasion by Repressing S1PR2 Protein Expression through NF-Κb Signaling. Arch. Biochem. Biophys. 643, 7–13. doi:10.1016/j.abb.2018.02.009
Wu, L., Liu, D., and Yang, Y. (2019a). Enhanced Expression of Circular RNA Circ-DCAF6 Predicts Adverse Prognosis and Promotes Cell Progression via Sponging miR-1231 and miR-1256 in Gastric Cancer. Exp. Mol. Pathol. 110, 104273. doi:10.1016/j.yexmp.2019.104273
Wu, Y., Wang, X., Xu, F., Zhang, L., Wang, T., Fu, X., et al. (2020). The Regulation of Acetylation and Stability of HMGA2 via the HBXIP-Activated Akt-PCAF Pathway in Promotion of Esophageal Squamous Cell Carcinoma Growth. Nucleic Acids Res. 48, 4858–4876. doi:10.1093/nar/gkaa232
Wu, Y., Xie, Z., Chen, J., Chen, J., Ni, W., Ma, Y., et al. (2019b). Circular RNA circTADA2A Promotes Osteosarcoma Progression and Metastasis by Sponging miR-203a-3p and Regulating CREB3 Expression. Mol. Cancer 18, 73. doi:10.1186/s12943-019-1007-1
Yang, C., Yuan, W., Yang, X., Li, P., Wang, J., Han, J., et al. (2018). Circular RNA Circ-ITCH Inhibits Bladder Cancer Progression by Sponging miR-17/miR-224 and Regulating P21, PTEN Expression. Mol. Cancer 17, 19. doi:10.1186/s12943-018-0771-7
Ye, J., Zhou, G., Zhang, Z., Sun, L., He, X., and Zhou, J. (2016). Poly (C)-binding Protein 2 (PCBP2) Promotes the Progression of Esophageal Squamous Cell Carcinoma (ESCC) through Regulating Cellular Proliferation and Apoptosis. Pathol. - Res. Pract. 212, 717–725. doi:10.1016/j.prp.2016.05.008
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). ClusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118
Yu, J., Xie, M., Ge, S., Chai, P., Zhou, Y., and Ruan, J. (2020). Hierarchical Clustering of Cutaneous Melanoma Based on Immunogenomic Profiling. Front. Oncol. 10, 580029. doi:10.3389/fonc.2020.580029
Zhang, F., Qin, S., Xiao, X., Tan, Y., Hao, P., and Xu, Y. (2019). Overexpression of LIMD2 Promotes the Progression of Non-small C-ell L-ung C-ancer. Oncol. Lett. 18, 2073–2081. doi:10.3892/ol.2019.10473
Zhang, J., Liu, H., Hou, L., Wang, G., Zhang, R., Huang, Y., et al. (2017). Circular RNA_LARP4 Inhibits Cell Proliferation and Invasion of Gastric Cancer by Sponging miR-424-5p and Regulating LATS1 Expression. Mol. Cancer 16, 151. doi:10.1186/s12943-017-0719-3
Zheng, Y., Chen, Z., Han, Y., Han, L., Zou, X., Zhou, B., et al. (2020). Immune Suppressive Landscape in the Human Esophageal Squamous Cell Carcinoma Microenvironment. Nat. Commun. 11, 6268. doi:10.1038/s41467-020-20019-0
Keywords: esophageal cancer, competing endogenous RNA, LIM domain containing, prognostic, predictive
Citation: Chen Y, Huang X, Zhu K, Li C, Peng H, Chen L, Huang Z, Zhang Y, Weng G, Xiao T, Chen J and Xu Y (2021) LIMD2 is a Prognostic and Predictive Marker in Patients With Esophageal Cancer Based on a ceRNA Network Analysis. Front. Genet. 12:774432. doi: 10.3389/fgene.2021.774432
Received: 12 September 2021; Accepted: 22 October 2021;
Published: 18 November 2021.
Edited by:
Graziano Pesole, University of Bari Aldo Moro, ItalyReviewed by:
YM Tsui, The University of Hong Kong, Hong Kong SAR, ChinaMaria La Mantia, University of Palermo, Italy
Copyright © 2021 Chen, Huang, Zhu, Li, Peng, Chen, Huang, Zhang, Weng, Xiao, Chen and Xu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yuanji Xu, xuyuanji@fjmu.edu.cn; Junqiang Chen, junqiangc@163.com
†These authors have contributed equally to this work and share first authorship