The study was approved by the Institutional Review Board Medical Ethics Committee of the Kidwai Memorial Institute of Oncology, Bangalore, India (KMIO/MEC/021/24 dated 30 November 2016), and conducted according to the guidelines of the Declaration of Helsinki.

Informed consent was obtained from all subjects involved in the study.

The raw reads from IOB-KMIO and GSE32424 datasets were analyzed using the workflow depicted in Supplementary Figure S1. The raw reads were quality checked using FastQC, quality trimmed using Trimmomatic (v0.39), and clean reads were obtained. Furthermore, clean reads were mapped and aligned to the human genome GRCh38.105 using the STAR RNA-seq aligner (v2.7.10a) (Supplementary Tables S6, S7). Furthermore, the aligned reads were quantified using HTSeq (v1.99.2). R (v4.1) scripts were used to merge the HTSeq output files, and the R Bioconductor package, biomaRt (Durinck et al., 2009), was used to convert the Ensembl IDs to HGNC gene symbols. The ESCC RNA-seq log (read count +1) normalized data from TCGA were transformed to raw read counts using R scripts. Furthermore, the read counts corresponding to ion channels, lipid metabolism- and EMT-related genes were parsed using R scripts. Supplementary Figure S2 depicts the overall workflow used for the identification of ion channels and their interactions with lipid metabolism and EMT.

Differentially expressed ion channels, lipid metabolism- and EMT-related genes in the IOB-KMIO and GSE32424 datasets were identified using the R Bioconductor package DESeq2 (v4.1.2) (Love et al., 2014). An adjusted p-value (padj) ≤ 0.05 and log2-fold change (log2FC) > 0.6 were set as the criteria for the selection of upregulated genes, and the genes with padj ≤ 0.05 and log2FC < −0.6 were selected as downregulated genes. A non-redundant list with all DEGs obtained from the two datasets was generated.

The expression profiles of ESCC patients from IOB-KMIO and GSE32424 datasets corresponding to DEGs were parsed. Furthermore, these profiles were provided as input to the R Bioconductor package weighted gene co-expression network analysis (WGCNA) (Zhang and Horvath, 2005; Langfelder and Horvath, 2008) to construct the scale-free co-expression networks. The read counts provided as input were DESeq2-normalized. On the basis of a scale-free network model, an adjacency correlation matrix was calculated (Langfelder and Horvath, 2008). A suitable soft threshold power that ensured a scale-free network was selected based on the pickSoftThreshold function. Thereafter, a signed topological overlap matrix (TOM) and the corresponding dissimilarity matrix (1-TOM) were derived from the obtained adjacency matrix. The dissimilarity matrix was used for performing hierarchical clustering with the DynamicTreeCut algorithm (Langfelder et al., 2008) in order to group the genes with similar expression profiles into the same gene modules.

The expression profiles of each obtained modules were epitomized by module eigengenes (MEs) and correlated with binary trait (normal and tumor) values, and a relationship between the module and traits was obtained. The p-value was calculated as a confidence measure for the module–trait relationship. Furthermore, gene significance (GS) and module membership (MM) values were calculated to estimate the association of individual genes in the module with the binary trait and the correlation of gene expression profiles with MEs, respectively. The intramodular gene connectivity, with the binary trait, was estimated on the basis of GS vs. MM plots (Langfelder and Horvath, 2008). Thereafter, the most correlated modules with a significant p-value were selected from each dataset. The obtained weighted gene co-expression networks were exported to Cytoscape with the weight threshold value set to 0.02 for visualization.

In this study, 76-gene EMT signature-based (GS76), multinomial logistic regression-based (MLR), and Kolmogorov–Smirnov test-based (KS) methods developed in the past (Chakraborty et al., 2020; George et al., 2017; Mandal et al., 2021) were used for the calculation of EMT scores of transcriptomic profiles of samples from patients with ESCC. A higher KS and MLR score indicated a more mesenchymal sample, and a higher GS76 score indicated a more epithelial sample. MLR quantifies the extent of EMT on a scale of [0, 2], and KS scores are defined on a scale of [−1, +1] (Chakraborty et al., 2020). The scatter plots between different EMT scores were plotted using the R package “ggplot2” (Wickham, 2016) with correlation coefficients and p-values calculated using the R package “corrplot” (Wei et al., 2017) by Pearson’s method.

The estimated EMT scores of IOB-KMIO samples by the GS76, MLR, and KS methods were correlated with the expression profiles of the module ion channels and lipid metabolism-related genes interacting with EMT-related genes in ESCC.

To infer the interactions between the genes in the modules, protein–protein interaction networks (PPINs) were constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (v11.0) database (von Mering et al., 2003; Szklarczyk et al., 2021). The interactions with a medium confidence and false discovery rate of 1% were further visualized in Cytoscape (v3.9) (Shannon et al., 2003). The unconnected nodes in the network were discarded. Thereafter, NetworkAnalyser, a Cytoscape plugin, was used to calculate the various properties of the network, including degree centrality and betweenness centrality. The top 15 nodes in the network with the highest degree centrality were selected as the hub nodes in the networks.

The Molecular Complex Detection (MCODE) (Bader and Hogue, 2003) tool, a Cytoscape plugin, was used to identify the dense clusters representing the highly interacting nodes within the PPI networks. The parameter maximum depth was set to 100 degree, k-core was maintained as 2, and the node score was set to 0.2. Thereafter, the clusters consisting of ion channels, and lipid metabolism- and EMT-related nodes were selected to identify ion channels interacting with both lipid metabolism- and EMT-related nodes.

An extensive literature search was carried out using PubMed, and reactions like activations, inhibitions, and translocations that describe the interactions between the clustered ion channels, lipid metabolism- and EMT-related genes were annotated. Thereafter, a signaling pathway map consisting of the deregulated ion channels, lipid metabolism- and EMT-related genes was generated using PathVisio (version 3.3.0) in the Graphical Pathway Markup Language (GPML) format.

TCGA-ESCA dataset was used to compare the alterations in the ion channels obtained from IOB-KMIO and GSE32424 datasets. Differentially expressed ion channels, lipid metabolism- and EMT-related genes in the RNA-seq TCGA-GTEx ESCC dataset were identified using the R Bioconductor package DESeq2 (v4.1.2). An adjusted p-value (padj) ≤ 0.05 and log2FC > 0.6 were set as the criteria for the selection of upregulated genes, and the genes with padj ≤ 0.05 and log2FC < −0.6 were selected as downregulated genes.

Furthermore, the clinical data on patients with ESCC were extracted using the R Bioconductor package “RTCGA” for survival analysis. To determine the correlation between ion channel expression and a 5-year survival rate of patients with ESCC, Cox proportional hazard regression analysis was performed. The R Bioconductor package “survival” was utilized to calculate the hazard ratio (HR) with a confidence interval of 95% and the log rank p-values. Kaplan–Meier (KM) survival plots were generated to visualize the differences between the high and low expressions of the putative ion channels using the R Bioconductor package “survminer”. The median value was used as the group cut-off criteria.

Ion channels are potential clinical targets for restorative intervention due to their localization to the plasma membrane and different regulation mechanisms (Anderson et al., 2019). The list of known drugs interacting with all the human genes was collected from The Drug–Gene Interaction Database (DGIdb) (Freshour et al., 2021). This list was parsed with the list of potential ion channels interacting with lipid metabolism- and EMT-related genes.

A total of 282 genes, including 39 ion channels, 108 lipid metabolism genes, and 135 EMT-related genes, were found as differentially expressed in the IOB-KMIO dataset. In total, 1,013 genes, including 89 ion channels, 433 lipid metabolism genes, and 491 EMT-related genes, were found differentially expressed in the GSE32424 dataset (Table 1 and Supplementary Table S8).

Figure 1(i) represents the overlaps between the upregulated and downregulated genes in IOB-KMIO and GSE32424 datasets. A total of 1,081 DEGs were identified. Figure 1(ii) represents the overlaps between ion channels, and lipid metabolism- and EMT-related genes across IOB-KMIO and GSE32424 datasets.

WGCNA of the identified DEGs led to the detection of gene modules consisting of genes co-expressed in each of the dataset. For the IOB-KMIO dataset, a soft thresholding power of 14 was chosen and the minimum module size of 50 with a deepSplit parameter set to 0 was used to construct co-expression networks. Similarly, for the GSE32424 dataset, a soft thresholding power of 18 was chosen and the minimum module size of 50 with a deepSplit parameter set to 0 was used to construct co-expression networks. The IOB-KMIO dataset was split into five modules. Based on the p-value and correlation values, the blue and brown modules were selected (Figure 2). GSE32424 was split into two modules—turquoise and blue. Both the modules showed a good correlation among the genes within the modules (Supplementary Figure S3).

The EMT scores obtained using GS76, MLR, and KS metrics for the samples in IOB-KMIO are provided in Supplementary Table S9. A negative GS76 score with a positive KS score and a higher MLR score indicates the sample to be a mesenchymal sample. Similarly, the sample with a positive GS76 score and a negative KS score with a lower MLR score could be an epithelial sample. The correlation between the scores corresponding to the respective method in the dataset was calculated and is depicted in Figure 3 and Supplementary Table 10. MLR and KS scores correlated positively with each other. GS76 scores correlated negatively with both MLR and KS scores.

PPI networks obtained through the STRING database for the co-expressed gene modules revealed the predicted PPIs between co-expressed DEGs (Supplementary Figure S4). The lipid metabolism- and EMT-related proteins with multiple interactions were identified as hub nodes in the network (Figure 4; Supplementary Figure S5). Several ion channels were found interacting with the hub nodes in the network. In the network generated from the blue module of the KMIO dataset, the ion channels GJA1, P2RX4, and GJC1 interacted with the hub nodes AKT1, CCL2, CDH2, and FN1. TRPV3 and ITPR3 were found interacting with EGFR in the network generated from the KMIO brown module dataset. In the network from the GSE32424 turquoise module dataset, 13 ion channels, TRPV3, GJB1, GJC1, P2RX4, ITPR3, TRPC6, RYR2, CLIC4, LRRC8A, TRPC1, GJA5, TRPV2, and KCNMA1, interacted with seven hub nodes that included EGFR, CDH2, IL6, AKT1, CAV1, EP300, and IGF1. VDAC2, ITPR2, VDAC1, GJB2, and CFTR were found interacting with HSPA4, MAPK3, and CDH1 in the network from the GSE32424 blue module (Table 2).

Furthermore, PPI networks were analyzed using the MCODE algorithm to identify dense networks within the PPI networks. In total, 11 clusters in the IOB-KMIO blue module PPI network, 10 clusters in the IOB-KMIO brown module PPI network, 26 clusters in the GSE32424 turquoise module PPI network, and 20 clusters in the GSE32424 blue module PPI network were identified. These network clusters represented the highly interacting protein clusters within the PPINs. Out of these, five clusters—one from the IOB-KMIO blue module, one from the IOB-KMIO brown module, two from the GSE32424 turquoise module, and one from the GSE32424 blue module—comprised ion channels together with lipid metabolism- and EMT-related proteins (Figures 5A–E) (Supplementary Table S11). Several ion channels present in the cluster including GABRR2, ANO1, and TRPV3 were found to be negatively correlated with GS76 and positively correlated with both MLR and KS scores, indicating that these ion channels have a higher probability toward stabilizing/acquiring a mesenchymal phenotype. Similarly, GABRQ, GABRA3, and TRPM7 negatively correlated with GS76 and positively correlated with MLR and KS scores, indicating that they may have a higher affinity toward maintaining the epithelial phenotype (Supplementary Table S9).

The significance of the interactions between ion channels, lipid metabolism- and EMT-related genes in various aspects of ESCC tumorigenesis can be laid down by understanding their role in cellular processes. Pathway maps enable the depiction of these cellular processes that may aid in better understanding the disease etiology (Sharma et al., 2020). Figure 6 summarizes the reaction events that are likely to occur in ESCC patients through interactions between deregulated ion channels, lipid metabolism- and EMT-related genes. The pathway map comprises of 77 molecules and 46 reactions. Key cellular processes, including insulin signaling and estrogen receptor-calcium signaling pathways deregulated in ESCC, are depicted in the pathway map.

The deregulated ion channels in IOB-KMIO and GSE32424 were compared with the deregulated ion channels in TCGA-ESCC dataset. In total, 240 ion channels were found deregulated in TCGA-ESCC dataset, of which 18 were found deregulated in IOB-KMIO and 49 were present in the deregulated list of ion channels from the GSE32424 dataset. In total, 11 ion channels, GJB2, GRIN2D, TPCN1, KCTD6, RYR2, HTR3B, SCNN1B, CLIC3, KCNE3, CNGA1, and CLCN1, were identified as deregulated in all three datasets (Supplementary Figure S6). Out of 16 ion channels that interacted with EMT- and lipid metabolism-related genes, 14 were found deregulated in at least two datasets. TRPM2 and CLCN1 were found deregulated in all three datasets (Supplementary Table S12).

The interconnection of ion channels with the survival of patients with ESCC was estimated through KM plots and Cox proportional hazard regression analysis. Ion channels with HR > 1 were selected as major ion channels for the survival of ESCC patients (Figure 7). Six ion channels—GJA1 [HR: 3.24, p (HR): 0.07], GABRQ [HR: 1.68, p (HR): 0.19], GABRA3 [HR: 2.42, p (HR): 0.11], KCNN4 [HR: 3.88, p (HR): 0.04], TRPC1 [HR: 2.72, p (HR): 0.09], and ITPR3 [HR: 4.09, p(HR): 0.04]—exhibited HR > 1. Patients with high expression of ITPR3 (p-value: 0.03) had a better survival rate. Additionally, 782 gene–drug interactions were identified between the 16 ion channels—ANO1, CLCN1, CLCN2, GABRA3, GABRA4, GABRE, GABRP, GABRQ, GABRR2, GJA1, ITPR3, KCNN4, TRPC1, TRPM2, TRPM7, and TRPV3 drugs (Supplementary Table 13).

GJA1 was found interacting with a suppressor of cytokine signaling 3 (SOCS3) and other proteins involved in the mechanism leading to tumor progression through insulin resistance-related pathways. Obesity is one of the major risk factors contributing toward the onset of ESCC (Etemadi et al., 2012). Elevated levels of IL-6 observed in obese individuals could lead to an increase in the levels of suppressor of cytokine signaling 1 (SOCS1) and SOCS3 in white adipose tissue and the liver. At the molecular level, this increase could impair the insulin functionality by binding of the SOCS proteins to insulin receptors, IRS1 and IRS2. Another mechanism leading to the development of insulin resistance by overexpression of the SOCS protein was through reduction in the tyrosine phosphorylation of IRS1 and IRS2 required for insulin transduction (Wunderlich et al., 2013). This signal transduction is carried out through phosphoinositide 3-kinases (PI3K) and the MAPK pathway accounting for mitogenesis, transformation, and apoptosis protection in tumor cells (Arcidiacono et al., 2021). Wunderlich and others have shown the involvement of SOCS1 and SOCS3 in insulin resistance through the elevation of interleukin-6 (IL-6), thus leading to the activation of the JAK-STAT pathway.

GJA1 was found interacting with the signal transducer and activator of transcription (STAT), and other proteins involved in this pathway. STAT3 aberrant activation could lead to esophageal cellular transformation, angiogenesis, EMT, and metastasis. Sugase and colleagues have identified the association of the JAK/STAT pathway in the prognosis of patients with ESCC through immunohistochemistry studies (Sugase et al., 2017). The correlation between the phosphorylation levels of STAT3 and JAK1 was significant in predicting the poor prognosis of ESCC patients, indicating their importance in better management of ESCC (Rah et al., 2022).

ALDH1A1 was found to interact with the ion channels GABRE, GABRR2, GABRQ, GABRA3, ANO1, and GJA1 (Figure 5B). The PI3K-AKT pathway known to be deregulated in ESCC is activated by the overexpression of ALDH1A1. ALDH1A1 increases the expression of AKT serine/threonine kinase 1 (AKT1) and p-AKT. The phosphorylation of AKT activates several substrates that regulate cellular proliferation, survival, and motility (Shi et al., 2019). ALDH1A1 plays a significant role in maintaining stem cell properties in patients with ESCC through the AKT pathway. ALDH1A1 could be a potential target in patients with ESCC (Wang et al., 2020).

The estrogen receptor-Ca²⁺ signaling pathway was found to be associated with male predominance of ESCC. The estrogen receptors present on the plasma membrane activate multiple intracellular secondary messengers, including Ca²⁺ signaling. An increase in the cytoplasmic Ca²⁺ levels occur due to the release of Ca²⁺ from the endoplasmic reticulum through ion channels, such as ITPR3 (Zhang et al., 2017).

Cox proportional regression analysis with high HR indicated the alterations in the ion channels GJA1, GABRQ, GABRA3, KCNN4, TRPC1, and ITPR3 to possibly be associated with poor survival in patients with ESCC. Patients with low expression of GJA1 (p-value: 0.14), KCNN4 (p-value: 0.56), and TRPC1 (p-value: 0.55) and higher expression of GABRQ (p-value: 0.14), GABRA3 (p-value: 0.57), and ITPR3 (p-value: 0.032) may show better survival. However, the number of patients with ESCC was limited to attain statistical significance to support the expression difference using a log rank p-value <0.05 in survival analysis except ITPR3. Of these, GJA1 and ITPR3 were found interacting with hub proteins in the interactome and with other proteins in the highly interacting network clusters. Several ion channels found deregulated in the present study, including GJA1, and ITPR3, interacted with various known drugs. The presence of ion channels on the plasma membrane makes them accessible targets to therapeutic drugs. Moreover, their putative significance as diagnostic and prognostic markers and the availability of multiple drugs interacting with the potential ion channels make them ideal targets for clinical studies. In vivo mouse model studies showed that SKF96365, a phenylethylimidazole, and an inhibitor of store operated calcium channels (SOCE) resulted in decreased tumor volume by 50% in 4 weeks (Leanza et al., 2016). The inhibition of TRPC6 in ESCC cells leads to the suppression of proliferation and induction of the G2/M phase arrest (Zhang et al., 2018). The inhibition of TRPV2 by tranilast; inhibition of ANO1 by T16Ainh-A01, digallic acid, and tannic acids; and the inhibition of voltage-gated sodium channels by valproic acid were reported as potential targeted therapeutic agents in esophageal squamous cell carcinoma (Shiozaki et al., 2018; Capatina et al., 2022; Bill et al., 2014; Wheler et al., 2014). However, these channels are ubiquitously present, and toxicity toward normal cells on targeting these channels has not been addressed in depth. Further in vivo and in vitro studies could be carried out to elucidate the role of identified potentially altered ion channels interacting with genes of lipid metabolism- and EMT-related genes.