Stemness-related gene signature for predicting therapeutic response in patients with esophageal cancer.

Background: Extensive research has indicated that tumor stemness promotes tumor progression. However, the underlying role of stemness-related genes (SRGs) in esophageal cancer (ESCA) remains unclear.
Methods: This study identified differentially expressed stemness-related (DESR) messenger RNAs (mRNAs), microRNAs (miRNAs), and long non-coding RNAs (lncRNAs) in ESCA, and correlated them with the clinical features of patients with ESCA to develop a prognostic risk assessment model. Functional analysis, protein-protein interaction (PPI) analysis, competing endogenous RNA (ceRNA) networks, and tumor-infiltrating immune cell analyses were performed to corroborate the results obtained from the model.
Results: Correlation analysis of the stemness enrichment scores revealed 1,106 DESR genes (DESRGs), 84 DESRmiRNAs, and 320 DESRlncRNAs were identified from The Cancer Genome Atlas Esophageal Carcinoma (TCGA-ESCA) dataset. Network clustering was performed and the top 20 connection points were identified, including CDC20 that connects to 136 adjacent nodes. A ceRNA network was constructed, including 17 DESRmiRNAs, 44 DESRlncRNAs, and 55 DESRGs. Conclusions: NCAPG [log2fold change (FC) =1.81; q value =2.68×10-11] was significantly upregulated in ESCA and positively correlated with resting natural killer (NK) cells, suggesting that human NK cells rest via the overexpression of NCAPG in ESCA. hsa-miR-1269a is significantly upregulated in ESCA patients with poor prognostic features. CD4+ resting memory T cells (P<0.01) were significantly negatively correlated with hsa-miR-1269a. The insights presented in this study will contribute to the development of innovative therapeutics for the treatment of patients with ESCA. 2022 Translational Cancer Research. All rights reserved.

Introduction

Esophageal cancer (ESCA) is the second deadliest gastrointestinal cancer. Patients generally have poor prognosis, with a 5-year survival rate of 40% for early-stage patients (1). There are two pathological types of ESCA that show regional specificities because of differences in living habits and genetic backgrounds. Esophageal squamous cell carcinoma is more prevalent in Southeast Asia, while esophageal adenocarcinoma is more frequent in Europe and America (2). Patients with ESCA can be divided into resectable patients and those who are inoperable (3). Traditional treatments such as surgery, radiotherapy, chemotherapy, immunotherapy, and targeted therapy, have limited efficacy in both patients with resectable and those with inoperable ESCA, leading to poor outcomes (4). In patients with advanced metastatic ESCA, nanoparticles AbraxaneTM and DoxilTM can be used as treatment options (3). Carboplatin and paclitaxel are the most common preoperative treatment agents for ESCA (5), while capecitabine and fluorouracil are administered postoperatively (6). Minimally invasive esophagectomy is a common surgical method for patients with operable ESCA to reduce postoperative morbidity. For patients with inoperable ESCA, cebox and cisplatin are used in combination to prolong survival (3). PembrolizumabTM, which targets specific genetic features, has been approved by the Food and Drug Administration (FDA) for the treatment of advanced ESCA (7). Therefore, current research on ESCA should be continuously adjusted to develop more effective treatment strategies.

Cancer stem cells (CSCs) are highly relevant for tumorigenesis, progression, recurrence, and metastasis of various cancers (8). The messenger RNA (mRNA) expression-based stemness index (mRNAsi) is used to evaluate the unique characteristics of CSCs and to assess tumor development (9). There is mounting evidence to support the existence of CSCs in lung adenocarcinoma (10), breast cancer (11), colorectal cancer (12), and prostate cancer (13), as well as their roles in metastasis, drug resistance, and cancer adaptation to changing microenvironments. For instance, single-cell transcriptome analysis revealed the landscape of intra-tumoral heterogeneity and stemness-related subpopulations of liver cancer cells (14).

Stemness-related genes (SRGs) have been identified in numerous cancers, such as cervical squamous cell carcinoma (15), colorectal cancer (16), renal clear cell carcinoma (17), lung cancer (18), gastric cancer (19), and ESCA (20). Moreover, SRGs have been shown to promote the formation, progression, and metastasis of different cancers (21). The AGR3 gene has been shown to promote the stemness of CSCs by upregulating SRGs via Frizzled 4 (FZD4), which is involved in the Wnt/β-catenin signaling pathway (22). Study examining ESCA stemness and clinical characteristics have indicated that higher-stage and metastatic tumors feature great phenotypic dedifferentiation (20). This study examined how SRGs promote ESCA formation, progression, and metastasis. The differentially expressed stemness-related (DESR) genes (DESRGs), DESR microRNAs (DESRmiRNAs), and DESR long non-coding RNAs (DESRlncRNAs) were firstly identified and functional analyses were performed. CSCs could change tumor microenvironment into immunosuppressive (23). Tumors with significant positive DESRGs, DESRlncRNAs, and DESRmiRNAs expression may exist higher immune infiltration in the tumor microenvironment (21). The DESRGs, DESRmiRNAs, and DESRlncRNAs were used to develop a prognostic risk assessment model and to construct a protein-protein interaction (PPI) network and a competing endogenous RNA (ceRNA) network, which together with the tumor-infiltrating immune cell analyses, corroborated the model results (). Potential stemness-related prognostic markers for ESCA were identified. We present the following article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-1723/rc).

Figure 1

Flowchart showing the process involved in this study. mRNA, messenger RNA; lncRNA long non-coding RNA; miRNA, microRNA; TCGA, The Cancer Genome Atlas; GSEA, gene set enrichment analysis; ssGSEA, single-sample GSEA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ceRNA, competing endogenous RNA.

Methods

Data collection and analysis of differential gene expression

The ESCA gene expression profile was obtained from The Cancer Genome Atlas-ESCA (TCGA-ESCA) database containing 11 normal tissues and 152 tumor tissues, and analyzed using the R package TCGA biolinks (24). Among these, 145 tumor tissues had complete clinical information, including age, sex, survival status, and tumor, node, metastasis (TNM) stage. The read counts were converted to transcripts per kilobase million (TPM) to compare the relative gene expression level between samples, while the miRNA data were converted to reads per kilobase million (RPKM) (25).

Differential gene expression between tumor and normal tissue samples was analysis using the R package DESeq2 (26). The significant threshold for the mRNA and lncRNA differential expression was defined as the log2 of the absolute fold change (FC) value greater than 1 (log2FC >1) and adjusted P value less than 0.05 (q value <0.05). The threshold to determine differentially expressed miRNA was set as log2FC >0.5 and q value <0.05. The differential analysis results are displayed as heat maps using the R package pheatmap (https://CRAN.R-project.org/package=pheatmap) and as volcano maps using the ggplot2 R package (27). The DESRGs were obtained using the gene set enrichment analysis (GSEA) (https://www.gsea-msigdb.org/gsea/index.jsp) tool of the GSEA software (28). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Construction of the PPI network and the ceRNA network

The single-sample GSEA (ssGSEA) method (29) was used to obtain the enrichment score for cell stemness of each sample based on the ESCA expression data and the DESRG set. To identify the DESRmiRNAs and DESRlncRNAs, correlation analysis was performed between the differentially expressed mRNA and the lncRNA or miRNA levels. Significant correlation coefficient greater than 0.2 (ρ>0.2; P<0.05) indicated strong positive correlation.

To construct the PPI network, the STRING database (https://string-db.org/), containing 2,031 species, 9.6 million proteins, and 1.38 million PPIs was used. Our PPI analysis included experimental results, text-mined results from PubMed abstracts, and synthesized and predicted results from bioinformatics analyses (30). A correlation coefficient above 0.7 (ρ>0.7; P<0.05) indicated a strong association between the differentially expressed genes and cell stemness. The PPI results were exported from the STRING database and visualized using the Cytoscape StringApp (31). In addition, the molecular complex detection (MCODE) plugin from Cytoscape was used to analyze hub genes in the PPI network.

MiRcode (https://bio.tools/miRcode) provides transcriptome-wide human miRNA target predictions based on comprehensive GENCODE gene annotations (https://www.gencodegenes.org/), including 10,000 lncRNA genes (32). Therefore, MiRcode was used to identify the miRNAs and lncRNAs in our dataset. The TargetScan database was used to predict the biological targets of the identified miRNAs by searching for the presence of conserved 6mer-8mer sites on the mRNA sequences matching the seed region of each miRNA (33). The miRDB database (http://mirdb.org/) was used to confirm the predicted miRNA targets and functional annotations by implementing its MirTarget tool, which uses thousands of miRNA-target interactions (MTIs) from high-throughput sequencing experiments (34). miRTarBase (https://mirtarbase.cuhk.edu.cn) is an experimentally validated MTI database that has accumulated more than 360,000 MTIs obtained from natural language processing after a manual survey of relevant literature (35). For the analysis using miRTarBase, research articles were systematically screened to collect those related to miRNA functional studies. The MTIs from TargetScan, miRDB, and miRTarBase analysis were combined with the identified DESRGs, DESRmiRNAs, and DESRlncRNAs to construct a ceRNA network, which was visualized using Cytoscape.

Functional analysis

Gene Ontology (GO) analysis is a standard method for large-scale functional enrichment studies that includes terms related to biological processes (BPs), molecular functions (MFs), and cellular components (CCs) (36). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database for storing information on genomes, biological pathways, diseases, and drugs (37). GO and KEGG pathway enrichment analyses of DESRGs were performed using the clusterProfiler (38) R package, with a cutoff value of false discovery rate (FDR) <0.05 for statistical significance.

To investigate the differences in enriched BPs between normal and tumor tissues, GSEA was performed based on the TCGA gene expression profiling dataset of patients with ESCA (39). For GSEA analysis, the control “c2.cp.kegg.v7.1.entrez” dataset used was downloaded from the MsigDB database and compared with the gene expression profiling dataset of patients with ESCA (FDR <0.25).

Development and validation of a prognostic risk assessment model

Univariate Cox (uni-Cox) regression analysis (P<0.05) was performed for the DESRGs, DESRmiRNAs, and DESRlncRNAs to estimate ESCA prognosis. The uni-Cox results were visualized using forest plot. Survival analysis was performed on the DESRGs, DESRmiRNAs, and DESRlncRNAs and the results were displayed using the R package survminer (40).

Least absolute shrinkage and selection operator (LASSO) regression analysis was performed to filter the DESRGs, DESRmiRNAs, and DESRlncRNAs with a 10-fold cross-validation. Furthermore, the DESRGs, DESRmiRNAs, and DESRlncRNAs screened by the LASSO method were used for multivariate Cox (multi-Cox) proportional hazards regression and risk model construction. The LASSO glmnet package (41) was used to divide the TCGA samples with complete clinical information into a training set and a validation set at a ratio of 1:1. The training dataset was composed of the DESRGs, DESRmiRNAs, and DESRlncRNAs (single factor P<0.05). The accuracy of the final model was verified through the validation set. The R package ggrisk (https://CRAN.R-project.org/package=ggrisk) was used to display the risk factor and perform the survival analysis on the risk scores of the training and validation sets. The R package pROC was used to obtain the receiver operating characteristic (ROC) curve, nomogram, and standard curve, showing the accuracy of the model. Finally, multi-Cox analysis with clinical traits was used to determine whether the risk score is independent of other traits.

Tumor-infiltrating immune cells analyses

CIBERSORT (https://cibersort.stanford.edu/index.php) was used to analysis 22 tumor-infiltrating immune cell types in each tissue, using the LM22 signature matrix with 1,000 permutations (42). The gene expression matrix data (TPM) was uploaded to CIBERSORT combined with the LM22 eigengene matrix, and the immune cell infiltration matrix was derived. A histogram was drawn showing the distribution of the 22 infiltrating immune cells in each sample. The R package corrplot (https://github.com/taiyun/corrplot) was used to calculate and display the correlation between the prognosis-related mRNA, lncRNA, miRNA, and the 22 tumor-infiltrating immune cell types.

Statistical analyses

All statistical analyses were performed using R software version 4.1.3 (https://www.r-project.org/). The ROC curves assessed the accuracy of the risk score in estimating the prognosis by calculating the area under the ROC curve (AUC). Spearman’s correlation (ρ) was used to determine all correlations between two samples. All statistical analyses were two-sided, with a significance of P<0.05.

Results

Differential expression analyses

DESeq2 was used to analyze the differential expression between TCGA-ESCA tumors [152] and normal tissues [11]. Baseline data of TCGA patients with ESCA was detected in . A total of 3,504 differentially expressed mRNAs (including 1,900 upregulated and 1,604 downregulated mRNAs; ), 1,564 differentially expressed lncRNAs (including 993 upregulated and 571 downregulated lncRNAs; ), and 209 differentially expressed miRNAs (including 130 upregulated and 79 downregulated miRNAs; ) were identified. Heatmaps were used to identify the most significant differentially expressed mRNAs, lncRNAs, and miRNAs ().

Table 1

Baseline data of TCGA patients with ESCA

Type	Alive (n=98)	Deceased (n=46)	P value
Tumor grade, n (%)			0.238
G1	16 (16.3)	2 (4.3)
G2	38 (38.8)	19 (41.3)
G3	21 (21.4)	12 (26.1)
GX	23 (23.5)	13 (28.3)
Gender, n (%)			0.448
Female	17 (17.3)	5 (10.9)
Male	81 (82.7)	41 (89.1)
Age (years), n (%)			0.912
<60	43 (43.9)	19 (41.3)
≥60	55 (56.1)	27 (58.7)
Stage, n (%)			0.0011
I	17 (17.3)	3 (6.5)
II	48 (49.0)	14 (30.4)
III	31 (31.6)	22 (47.8)
IV	2 (2.0)	7 (15.2)
T staging, n (%)			0.191
T1	18 (18.4)	9 (19.6)
T2	25 (25.5)	9 (19.6)
T3	55 (56.1)	26 (56.5)
T4	0 (0.0)	2 (4.3)
M staging, n (%)			0.00238
M0	86 (87.8)	34 (73.9)
M1	1 (1.0)	7 (15.2)
MX	11 (11.2)	5 (10.9)
N staging, n (%)			0.00909
N0	48 (49.0)	11 (23.9)
N1	33 (33.7)	28 (60.9)
N2	6 (6.1)	4 (8.7)
N3	6 (6.1)	0 (0.0)
NX	5 (5.1)	3 (6.5)

TCGA, The Cancer Genome Atlas; ESCA, esophageal cancer.

Figure 2

The DEGs. (A-C) The volcano plots show the DEmRNAs, DElncRNAs, and DEmiRNAs, respectively. The red nodes represent upregulated DEGs, the blue nodes represent downregulated DEGs, and the gray nodes represent DEGs that were not significantly different. (D-F) The heatmaps show the gene expression in TCGA-ESCA tumors samples and normal tissues. DEGs, differentially expressed genes; DEmRNAs, differentially expressed messenger RNAs; DElncRNAs, differentially expressed long non-coding RNAs; DEmiRNAs, differentially expressed microRNAs.

Stemness-related ceRNA network

Correlation analysis of the stemness enrichment scores revealed 1,106 DESRGs, 84 DESRmiRNAs, and 320 DESRlncRNAs (). Network clustering was performed () and the top 20 connection points were identified, including CDC20 that connects to 136 adjacent nodes (). A ceRNA network was constructed, including 17 DESRmiRNAs, 44 DESRlncRNAs, and 55 DESRGs ().

Figure 3

The stemness-related ceRNAs. (A) A PPI network diagram of the DESRGs. Blue indicates high expression in tumor tissue, while red indicates low expression in tumor tissue. (B-D) Gene network maps of the hubs screened with the cytoscape plugin MCODE. (E) Histogram depicting the PPI network nodes for the DESRGs. (F) A diagram of the stemness-related ceRNA network. Green presents miRNA, purple presents lncRNA, and orange presents mRNA. ceRNA, competing endogenous RNA; PPI, protein-protein interaction; DESRGs, differentially expressed stemness-related genes; MCODE, molecular complex detection; miRNA, microRNA; lncRNA, long non-coding RNA; mRNA, messenger RNA.

Functional analysis of the DESRGs

GO enrichment analysis for the DESRGs () identified nuclear division, organelle fission, mitotic nuclear division, DNA-dependent DNA replication, and DNA replication as enriched initiation terms. KEGG pathway enrichment analysis revealed that the DESRGs were mainly enriched in biological pathways such as cell cycle, DNA replication, human papillomavirus infection, Cushing syndrome, and protein digestion and absorption (). GSEA functional analysis revealed that the DESRGs were enriched in the following KEGG pathways: cytokine-cytokine receptor interaction, cell cycle, calcium signaling pathway, neuroactive ligand-receptor interaction, fatty acid metabolism, DNA replication, spliceosome, NOD-like receptor signaling pathway, vascular smooth muscle contraction pathway, and graft-versus-host disease (). The most enhanced pathways were cytokine-cytokine receptor interaction, cell cycle, and DNA replication (), whereas the most inhibited pathways were related to calcium signaling pathway, neuroactive ligand-receptor interaction, and fatty acid metabolism ().

Figure 4

Enrichment functional analysis of the DESRGs. (A,B) GO analysis showing enrichment of the DESRGs in the BP, CC, and MF. The network graph shows the overall enrichment. (C) KEGG analysis of the DESRGs. BP, biological process; CC, cellular component; MF, molecular function; DESRGs, differentially expressed stemness-related genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 5

GSEA analysis. (A) The general enrichment pathway of the Mountain Map Exhibition. (B) The most significant top three pathways activated in ESCA. (C) The top three pathways most significantly inhibited in ESCA. NES, normalize enrichment score; GSEA, gene set enrichment analysis; ESCA, esophageal cancer.

Survival analysis of the DESRGs, DESRmiRNAs, and DESRlncRNAs

Uni-Cox and survival analysis of the DESRGs, DESRmiRNAs, and DESRlncRNAs was performed to identify the correlation on the survival of patients with ESCA. The uni-Cox model (P<0.05) used 73 DESRGs, 7 DESRmiRNAs, and 23 DESRlncRNAs (Table S1). The pseudogene MTRNR2L13 and the chromosome 15 open reading frame 54 (C15orf54) were identified as high-risk factors associated with patient mortality. Notably, the centromere protein I (CENPI), PRICKLE2 antisense RNA 1 (PRICKLE2-AS1), and human miRNAs hsa-miR-101-1 and hsa-miR-101-2 were found to negatively correlate with patient mortality. Survival analysis revealed six highly significant DESRGs, DESRmiRNAs, and DESRlncRNAs (Figure S1).

Construction and validation of the prognostic model

The TCGA-ESCA datasets were randomly divided into training and validation sets, and univariate analysis was performed to detect prognostic-related DESRGs. LASSO regression was used to identify meaningful genes in the uni-Cox analysis, including six DESRGs, namely, NCAPG, HIST1H4J, KIAA1456, SCUBE2, PCCA, and ARID3A, that were selected to build the prediction model (). The expression profiles of these selected genes were used to calculate the risk scores in the training and validation sets separately. The total risk score was calculated as follows: expression level of (NCAPG × 0.0245 + HIST1H4J ×1.2293 + KIAA1456 × 0.0358 + SCUBE2 × 0.1375 + PCCA × 0.0738 + ARID3A × 0.0134). The distribution of risk scores in patients with ESCA and the correlation between survival time and risk score in the training set revealed that patients with high scores were associated with a higher mortality (), indicating that high-risk patients generally had poor overall survival (OS). A similar trend was observed in the validation set (). Based on the median value of the risk scores, the patients were divided into a high-risk group and a low-risk group. The survival rate of patients in the high-risk group was significantly lower than that in the low-risk group in both the training (P<0.0001; ) and validation sets (P=0.01; ). The AUC for 1-, 2-, and 3-year survival rates in the training group were 0.995, 0.937, and 0.818, respectively (), while those for the validation group were 0.732, 0.838, and 0.955, respectively (). These results indicated that the DESRGs selected as prognostic indicators have robust potential for survival prediction in patients with ESCA.

Figure 6

Construction and validation of the DESRG prognostic model. (A) LASSO regression analysis plot of meaningful genes in uni-Cox analysis. (B) A risk factor plot for the training set. (C) A risk factor plot for the validation set. The Kaplan-Meier survival curve analysis revealed an OS in the high-risk group from the training set (D) and the validation set (E). The red line represents the high-risk group. The green line represents the low-risk group. Survival-dependent ROC curves validated the prognostic significance of DESRGs as prognostic indicators in the training set (F) and the validation set (G). The green line represents the 1-year AUC; the blue line represents the 2-year AUC; and the red line represents the 3-year AUC. A forest plot showing the uni-Cox (H) and multi-Cox (I) regression analyses. (J) A nomogram for predicting the survival probability of ESCA patients for predicting 1-, 2-, and 3-year survival outcomes of ESCA. (K) Calibration curves of the OS-nomogram. *, P<0.05; **, P<0.01; ***, P<0.001. AUC, area under the ROC curve; ROC, receiver operating characteristic; CI, confidence interval; OS, overall survival; DESRG, differentially expressed stemness-related gene; LASSO, least absolute shrinkage and selection operator; uni-Cox, univariate Cox; multi-Cox, multivariate Cox; ESCA, esophageal cancer.

Uni-Cox () and multi-Cox () regression analyses with clinical variables (risk score, sex, age, grade, and stage) revealed that the risk score and stage were independent prognostic indicators for patients with ESCA. A nomogram combining risk score and clinicopathological parameters was selected to predict ESCA patient prognosis (). Calibration curves were used to compare the actual probabilities of survival and the predicted survival rates. The results showed a significant correlation between the probability of OS and the actual survival rates for 1, 2, and 3 years ().

A prognostic model for the DESRlncRNAs and DESRmiRNAs was constructed and validated. The regression analysis identified SCOC-AS1, RP11-626E13.1, RP1-136B1.1, and RP11-16E23.3 as meaningful lncRNAs for the prediction model. The risk score for DESRlncRNAs was calculated as follows: expression level of SCOC-AS1 × 0.7209 + expression level of RP11-626E13.1 × 2.2361 + expression level of RP1-136B1.1 × 1.4229 + expression level of RP11-16E23.3 × 6.7802. Furthermore, the miRNAs selected for the risk analysis were hsa-miR-269a, hsa-miR-214, hsa-miR-3652, and hsa-miR-503. The DESRmiRNAs risk score was calculated as follows: expression level of hsa-miR-1269a × 0.003 + expression level of hsa-miR-214 × 0.03181 + expression level of hsa-miR-3652 × 0.1195 + expression level of hsa-miR-503 × 0.04297. The risk score results for the DESRlncRNAs (Figure S2) and DESRmiRNAs (Figure S3) showed that the model provided high accuracy in predicting the prognosis of patients with ESCA.

Tumor-infiltrating immune cells

To analyze the effect of tumor-infiltrating immune cells on the model, the CIBERSORT algorithm was used to calculate the degree of infiltration of 22 types of infiltrating cells and display their overall distribution scores () and correlations (). The correlation heatmap of the 22 types of infiltrating cells revealed that activated and resting natural killer (NK) cells, activated and resting mast cells, resting M0 macrophages, and CD4+ memory T cells had significant negative correlations. Monocytes were positively correlated with resting dendritic cells.

Figure 7

The landscape of immune infiltration in ESCA. (A) The difference in immune infiltration. (B) A correlation heat map of 22 types of immune cells. Red represents a positive correlation, blue represents a negative correlation. Candidate DESRGs (C), DESRlncRNAs (D), and DESRmiRNAs (E) that show significant positive correlations with immune infiltration cells are shown in red, whereas that show negative correlations with immune infiltration cells are shown in blue. NK, natural killer; ESCA, esophageal cancer; DESRG, differentially expressed stemness-related genes; DESRlncRNA, differentially expressed stemness-related long non-coding RNA; DESRmiRNA, differentially expressed stemness-related microRNA.

The association between infiltrating levels of the immune cells and the expression of DESRGs (), DESRlncRNAs (), and DESRmiRNAs () from patients with ESCA was analyzed. ARID3A was found to be positively correlated with infiltrating levels of immune cells, including follicular helper T cells (P<0.001) and activated CD4+ memory T cells (P<0.01). KIAA1456 expression was positively correlated with gamma delta (γδ) T cells (P<0.01) and negatively correlated with activated mast cells (P<0.01). The expression of NCAPG was positively correlated with resting NK cells (P<0.01). The expression of PCCA was positively correlated with resting CD4+ memory T cells (P<0.01) and plasma cells (P<0.01), and negatively correlated with activated CD4+ memory T cells (P<0.01). SCUBE2 expression was negatively correlated with activated dendritic cells (P<0.01). RP1-136B1.1 levels were positively correlated with activated NK cells (P<0.001) and CD8+ T cells (P < 0.01), and negatively correlated with follicular helper T cells (P<0.01), resting NK cells (P<0.001), and M0 macrophage (P<0.001). Finally, hsa-miR-1269a expression was negatively correlated with resting CD4+ memory T cells (P<0.01).

Discussion

To the best of our knowledge, this is the first study to perform a comprehensive characterization of ESCA stemness based on combined analyses of mRNA, miRNA, and lncRNA expression in ESCA tumors and normal tissues. A total of 1,106 DESRGs, 84 DESRmiRNAs, and 320 DESRlncRNAs were identified using the TCGA-ESCA dataset.

Non-coding RNAs, including miRNAs and lncRNAs, play important roles in transcriptional and translational regulation in several cancers (43). lncRNAs act as miRNA sponges to reduce their regulatory abilities (44). Interactions between lncRNAs, miRNAs, and mRNAs form ceRNA networks (45). Understanding the role of ceRNA networks in pan-cancer analysis is crucial (46). However, to date, no such study has been conducted on stemness-related ceRNA networks in ESCA. In this present investigation, we constructed a stemness-related ceRNA network based on 55 DESRGs, 17 DESRmiRNAs, and 44 DESRlncRNAs, and analyzed its correlation with the survival data of patients with ESCA. PRICKLE2-AS1, FAM13A-AS1, C15orf54, and LINC00304 were identified as high-risk DESRlncRNAs that correlated with high ESCA severity. Previous study has demonstrated that LINC00304 directly promotes cell proliferation and cell cycle progression in prostate cancer (47). PRICKLE2-AS1 has been shown to be differentially methylated in vulvar squamous cell carcinoma (48) and FAM13A-AS1 is involved in regulating RNA splicing and mRNA processing in thyroid cancer (49). Herein, PRICKLE2-AS1, FAM13A-AS1, and LINC00304 were shown to be downregulated in ESCA tumors, whereas C15orf54 was upregulated (log2FC =1.51; q value =0.02). C15orf54 is known to be positively associated with high-grade gastric cancer (50) and indeed, C15orf54 may be a potential biomarker for ESCA prognosis.

NCAPG is one of the DESRGs selected to build the prediction model. Increased expression of NCAPG promotes oncogenesis in non-small cell lung cancer (51), glioma (52), lung adenocarcinoma (53), prostate cancer (54), and esophageal squamous cell carcinoma (55). NCAPG (log2FC =1.81; q value =2.68×10−11) was significantly upregulated in ESCA and positively correlated with resting NK cells, suggesting that human NK cells rest via the overexpression of NCAPG in ESCA. Candidate DESRGs, DESRlncRNAs, and DESRmiRNAs that show significant positive correlations with immune infiltration cells, increased expression of SRGs and increased immune infiltration may be related to poor prognosis. Notably, the oncogene NCAPG in prostate cancer is directly targeted by the antitumor hsa-miR-145-3p, which is associated with cancer pathogenesis (54). In the present study, we found that hsa-miR-1269a was significantly upregulated in ESCA tumors (log2FC =3.98; q value =2.92×10−4) with poor prognostic features. Additionally, the number of resting CD4+ memory T cell as negatively corrected with hsa-miR-1269a (P<0.01). hsa-miR-145-3p is known as a diagnostic and oncogenic marker in hepatocellular carcinoma (56), colorectal cancer (57), and non-small cell lung cancer (58), and it may be involved in the molecular pathogenesis of various cancers, including ESCA (55).

ARID3A acts as an oncogene, and its overexpression affects the development of colorectal cancer (59), ovarian cancer (60), gastric cancer (61), hepatocellular carcinoma (62), and rectal cancer (63). This report revealed that ARID3A was upregulated and positively associated with infiltrating follicular helper T cells and activated CD4+ memory T cells, and thus, may play an immunological role in ESCA. Follicular helper T cells are a subset of T-helper cells, which provide critical support to B cells and are essential for maintaining humoral immune responses (64), such as germinal center formation, affinity maturation, and the development of high-affinity antibodies (65). Therefore, ARID3A is a promising independent prognostic biomarker corresponding to follicular helper T cells and activated CD4+ memory T cells.

A previous report showed that a de novo variant in the human HIST1H4J gene causes a syndrome analogous to HIST1H4C-associated neurodevelopmental disorder (66). Furthermore, HIST1H4J is the strongest negative predictor of ovarian cancer (67) and an important prognostic and chemoresistance predictor of ESCA (68). Herein, upregulation expression of HIST1H4J (log2FC =1.39; q value =1.87×10−3) was observed in patients with ESCA.

Hsa-miR-1269a serves as a diagnostic and oncogenic marker in hepatocellular carcinoma (56), colorectal cancer (57) and non‐small cell lung cancer (58). Our investigations showed that hsa-miR-1269a is significantly upregulated in ESCA patients with poor prognostic features. CD4+ resting memory T cells (P<0.01) were significantly negatively correlated with hsa-miR-1269a.

This is the first study to perform a comprehensive characterization of ESCA stemness based on a combined analysis of mRNA, miRNA, and lncRNA expression, and correlation with prognosis prediction. SRG network will facilitate improved treatment regimens for patients with ESCA and may promote the development of innovative therapeutics in the future.

The article’s supplementary files as