In Silico Analysis of Ion Channels and Their Correlation with Epithelial to Mesenchymal Transition in Breast Cancer

Simple Summary Breast cancer involves changes in the healthy cells of the breast resulting in rapid and abnormal division of cells that later spread to other parts of the body through the process of metastasis, which involves epithelial mesenchymal transition (EMT). Ion channels play a significant role in the switch from epithelial to mesenchymal transition through their contributions to cellular motility, cell volume regulation and cell cycle progression. Comprehensive computational analyses were performed to understand the role of ion channels in tumor/metastatic samples of breast cancer and their correlation with EMT. Abstract Uncontrolled growth of breast cells due to altered gene expression is a key feature of breast cancer. Alterations in the expression of ion channels lead to variations in cellular activities, thus contributing to attributes of cancer hallmarks. Changes in the expression levels of ion channels were observed as a consequence of EMT. Additionally, ion channels were reported in the activation of EMT and maintenance of a mesenchymal phenotype. Here, to identify altered ion channels in breast cancer patients, differential gene expression and weighted gene co-expression network analyses were performed using transcriptomic data. Protein–protein interactions network analysis was carried out to determine the ion channels interacting with hub EMT-related genes in breast cancer. Thirty-two ion channels were found interacting with twenty-six hub EMT-related genes. The identified ion channels were further correlated with EMT scores, indicating mesenchymal phenotype. Further, the pathway map was generated to represent a snapshot of deregulated cellular processes by altered ion channels and EMT-related genes. Kaplan–Meier five-year survival analysis and Cox regressions indicated the expression of CACNA1B, ANO6, TRPV3, VDAC1 and VDAC2 to be potentially associated with poor survival. Deregulated ion channels correlate with EMT-related genes and have a crucial role in breast cancer-associated tumorigenesis. Most likely, they are potential candidates for the determination of prognosis in patients with breast cancer.


Introduction
Breast cancer is a life-threatening disease and is one of the most common types of cancer prevalent in individuals of all age groups, with the majority of cases in females. It is also one of the most aggressive tumors with multifaceted gene expression levels at different stages of tumor progression. Though breast cancer starts as a local disease, mutations occur at distinct junctures of tumor differentiation, facilitating tumor cells to metastasize [1]. The process of metastasis starts locally by invading the host tissues surrounding the primary tumor and into the blood or lymphatic vessels through a mechanism known as epithelialmesenchymal transition (EMT) [2].
Ion channels are multimeric proteins located in the plasma membrane of the cell. They form a passageway extending from one side of the membrane to the other, thereby allowing the flow of ions into and out of the cell based on the electrochemical gradient. This generates a membrane potential that mediates a large number of biological functions within cells and across cell membranes [3]. Recent evidence has shown that ion channels are involved in the progression and pathology of various cancers [4][5][6]. Eag1 (KCNH1), a potassium channel, was reported to be involved in the proliferation and cell cycle of liposarcoma cells [7]. TRPM1, a cation-permeable channel, was found to be a prognostic marker for metastasis in melanoma [8]. TRPV6 was implicated in prostate adenocarcinoma and colorectal cancer cell lines [8]. Additionally, studies carried out by Ko et al. reported the association of several ion channels with various pathological features in breast cancer [3]. The activation of EMT programs is an expository mechanism for the gain of malignant phenotypes by epithelial cancer cells [9]. Often, this activation relies on signaling events between cancer cells and neighboring stromal cells [10]. Specific ion channels in various aspects of EMT induction have been reported, including the downregulation of CFTR, an ion channel that promoted EMT, migration and invasion in breast cancer [11]. Knockdown of hERG1 (KCNH2) expression led to reversion of the EMT profile in colorectal cancer cell lines, leading to reacquisition of the epithelial-like profile [12]. Expression of EMT transcription factors together with overexpression of TRPV1 cation channels led to hepatocarcinogenesis. Further inhibition of TRPV1 inhibited the growth of hepatocellular carcinoma cells [13]. The upregulation of ASIC1 and ASIC3 led to acidity-induced EMT through elevation of intracellular Ca 2+ concentration in pancreatic cancer cells. It was also shown that ASIC1 and ASIC3 positively correlated with mesenchymal marker vimentin and inversely correlated with epithelial marker E-cadherin in those cells [14].
The current study aims to identify altered ion channels and ion channels co-expressed with EMT-related genes in patients with breast cancer. Several bioinformatic analysis of ion channels and EMT-related genes led to the identification of altered ion channels in tumor/metastatic samples along with ion channels co-expressed with EMT-related genes in breast cancer patients. EMT scoring is a promising and versatile tool for the systematic investigation of EMT dynamics in cancer progression [15]. EMT scores for the expression profiles were calculated using 76-gene EMT signature based (GS76), multinomial logistic regression-based (MLR) and Kolmogorov-Smirnov test-based (KS) metrics [16]. Protein-protein interactions (PPIs) network analysis revealed 32 ion channels that interacted with 26 hub EMT-related genes in breast cancer. The correlation of potential ion channels with the EMT scores was used to statistically evaluate their potential in the EMT program. Furthermore, data mining of identified ion channels revealed the possible repercussions of altered expression of ion channels in cellular processes. A pathway map depicting various reactions, including the role of ion channels in serotonin signaling, insulin signaling, calcium signaling, adipocyte metabolism, nitric oxide signaling, glutamatergic signaling and osmotic stress, was generated. CACNA1B, ANO6, TRPV3, VDAC1 and VDAC2 were found to be prognostically significant. These analyses could thus identify ion channels that can be further studied to validate their potential role as molecular markers of breast cancer. Targeting these ion channels might lead to suppression of tumor growth. Understanding the role of these ion channels and their detailed mechanism with EMT programs may pave the way for developing new therapeutic strategies to improve the clinical outcomes of patients with breast cancer.

Material and Methods
The workflow of this study is depicted in Figure 1. EMT programs may pave the way for developing new therapeutic strategies to improve the clinical outcomes of patients with breast cancer.

Material and Methods
The workflow of this study is depicted in Figure 1.

Data Collection
Primarily, a list of 328 ion channels was taken from the HGNC database [17] and 1184 EMT-related genes from dbEMT (version 2.0) [18]. Publicly available transcriptomic datasets were analyzed to investigate gene expression in breast cancer. RNA sequencing data (RNA-Seq) of breast cancer patients and normal individuals corresponding to the ion channels and EMT-related genes were downloaded from the UCSC Xena portal, including the Cancer Genome Atlas (TCGA) and Genotype Tissue Expression (GTEx) data [19]. The number of metastatic samples obtained in the RNA-Seq dataset was limited; thus, the inclusion of microarray samples was taken into consideration. Microarray datasets GSE42568 and GSE52604 were retrieved from the Gene Expression Omnibus (GEO) database. The detailed description regarding the datasets have been mentioned in Tables S1-S3.

Data Collection
Primarily, a list of 328 ion channels was taken from the HGNC database [17] and 1184 EMT-related genes from dbEMT (version 2.0) [18]. Publicly available transcriptomic datasets were analyzed to investigate gene expression in breast cancer. RNA sequencing data (RNA-Seq) of breast cancer patients and normal individuals corresponding to the ion channels and EMT-related genes were downloaded from the UCSC Xena portal, including the Cancer Genome Atlas (TCGA) and Genotype Tissue Expression (GTEx) data [19]. The number of metastatic samples obtained in the RNA-Seq dataset was limited; thus, the inclusion of microarray samples was taken into consideration. Microarray datasets GSE42568 and GSE52604 were retrieved from the Gene Expression Omnibus (GEO) database. The detailed description regarding the datasets have been mentioned in Tables S1-S3.

Data Pre-Processing and Identification of Differentially Expressed Genes
The log (read count + 1) normalized RNA-Seq data were transformed to raw reads for input in DESeq2 [20]. The probe IDs of microarray datasets were converted to gene symbols using the annotation files GPL570 ((HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array) and GPL6480 (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F) obtained from the GEO database. Gene symbols were subsequently superimposed with the list of ion channels and EMT-related genes. Microarray expression profiles corresponding to those genes were set apart for the detection of differentially expressed genes (DEGs). DEGs for ion channel and EMT gene sets were identified individually. R Bioconductor package DESeq2 was used for RNA-Seq data, and the limma (v3.28.14) package [21] was used for the identification of DEGs in microarray datasets. In RNA-Seq and GSE42568 datasets, DEGs were identified for 3 subgroups-Normal vs. Tumor (HT), Tumor vs. Metastatic (TM) and Normal vs. Metastatic (HM). GSE52604 contained only metastatic and normal samples. The genes with an adjusted p-value (padj) ≤ 0.05 and log 2 fold change (log 2 FC) > 0.6 were selected as upregulated and the genes with a padj ≤ 0.05 and log 2 FC < −0.6 were selected as downregulated genes. A non-redundant list consisting of all DEGs for the three subgroups was prepared (Tables S4 and S5).

Construction of Co-Expression Networks of Altered Ion Channels and EMT-Related Genes
The expression profiles of the DEGs were further employed for the construction of scale-free co-expression networks using the weighted gene co-expression network analysis (WGCNA) R package (v4.0.0) [22,23]. Since the datasets belonged to different platforms, they were analyzed independent of each other. First, co-expression networks for DEGs in ion channels were constructed on the three datasets. Later, a combined gene set comprising both ion channels and EMT-related genes was formed for co-expression network construction. The genes in each of the datasets were filtered using goodSamplesGenes function in R. The adjacency correlation matrix was calculated based on the scale-free network model [23]. A suitable soft threshold power was selected as the soft-thresholding parameter to ensure a scale-free network using the pickSoftThreshold function. The obtained adjacency matrix was further used to derive the signed Topological Overlap Matrix (TOM) and the corresponding dissimilarity matrix (1-TOM). Based on the dissimilarity matrix, hierarchical clustering was performed to group genes with similar expression profiles into the same gene modules using the DynamicTreeCut algorithm [24].
For the selection of candidate modules, expression profiles of each module were summarized by module eigengenes (MEs) and correlated with the binary trait (normal, tumor, metastatic) values. Thus, the module-trait relationship was obtained, and the p-value was calculated as a confidence measure. The association of individual genes in the module with the binary trait was quantified by the Gene Significance (GS) value. Module membership (MM) was calculated as the correlation of gene expression profiles with the MEs. Further, intramodular gene connectivity with the binary trait was established using GS vs. MM plots [23]. Finally, the best correlated modules with a significant p-value were selected from each subgroup in each dataset. Parameters adjusted at each step of WGCNA are provided in Tables S7 and S8. The obtained networks of combined gene sets were exported to Cytoscape with the weight threshold value set to 0.02 for visualization.

Identification of Altered Ion Channels in Tumor and Metastatic States
The modules chosen subsequent to WGCNA of ion channels were examined for the identification of ion channels in tumor and metastatic states. The overlapping ion channels in the selected modules of three datasets belonging to HT were shortlisted as ion channels involved in the process of transition from normal to tumor. Similarly, overlapping ion channels in TM were selected as ion channels in the transition between the tumor to metastatic phenotype of breast cancer.

EMT Score Calculation
EMT scores of the transcriptomic profiles of samples from patients with breast cancer were calculated using previously developed methods, GS76, MLR and KS [16,25,26]. A well-defined set of gene signatures along with a classification algorithm forms the basis for EMT score calculation in these methods. A higher GS76 score indicates a more epithelial sample, but a higher MLR and KS score indicates a more mesenchymal sample. MLR quantifies the extent of EMT on a scale of [0, 2], and KS scores are defined on a scale of [−1, +1] [16].
EMT scoring methods aid in quantifying the extent of EMT phenotype of a sample. Different methods based on different gene signatures and algorithms have been developed for this purpose, due to which the scores obtained from these methods may vary for the same dataset. However, computing the scores with more than one method may capture the extent of EMT potential in all the genes involved in EMT process, thus expanding the search space for further analysis.

Protein-Protein Interaction Networks of Altered Ion Channels Co-Expressed with EMT-Related Genes
Protein-protein interaction networks (PPINs) of the selected modules obtained through WGCNA of the combined gene set were constructed using the STRING (v11.0) database [27] to deduce the association through known and predicted interactions between the genes in the modules. The interactions with a medium confidence and FDR of 5% were visualized in Cytoscape (v3.8) [28] and Gephi (v0.9.2) software. To further identify ion channels interacting with hub EMT-related genes, the networks were analyzed using NetworkAnalyser [29], a Cytoscape plugin that calculates the various properties of a network, including degree centrality, betweenness centrality and closeness centrality. The top 15 genes with the highest degree centrality measure were selected as hub genes from each module. Thereafter, the ion channels interacting with these genes were shortlisted.

Computation of Correlation of EMT Scores with Identified Ion Channels
The EMT scores of each sample estimated by the GS76, MLR and KS methods (mentioned in Section 2.5) were correlated with the expression profiles of the potential ion channels interacting with EMT-related genes and in tumor/metastatic states of breast cancer. Correlating EMT scores with expression profiles of genes of interest provides statistical evidence that supports the selection of those genes for further assessment.

Generation of Pathway Map of Altered Ion Channels and EMT-Related Genes
A data mining approach was used to annotate reactions involving identified ion channels and EMT-related genes in various cellular processes in patients with breast cancer. The articles were searched and screened from the PubMed database pertaining to possible effects of alterations in ion channels and EMT-related genes in breast cancer patients. Furthermore, various reactions, such as activation, inhibition and transportation, were annotated [30][31][32][33]. These reactions describe the translocation of ions and small molecules between subcellular compartments through ion channels and the role of ion channels with EMT-related genes in several cellular processes. Thereafter, a pathway map depicting the annotated reactions involving identified ion channels and EMT-related genes was generated in Graphical Pathway Markup Language (GPML) format using PathVisio (version 3.3.0) [34]-an opensource pathway drawing software. Nodes describe the entities (proteins, genes) and edges represent the relationship between the nodes in the map.

Survival Analysis of Identified Altered Ion Channels
The R Bioconductor "RTCGA" package was utilized to obtain the clinical data for survival analysis. Cox proportional hazard regression analysis was performed to determine the correlation between gene expression and 5-year survival rate of patients with breast cancer. The "survival" R Bioconductor package was used to calculate the log rank p-values and the hazard ratios (HR) with a confidence interval of 95%. The survival differences between high and low expressions of the putative ion channels were visualized by generating Kaplan-Meier (KM) survival plots using the R Bioconductor "survminer" package. The group cut-off criteria were set to median value. Further, the relationship between overall survival (OS) and expression profiles of putative ion channels was determined by KM OS plots using Gene Expression Profiling Interactive Analysis (GEPIA) [35]. GEPIA is a web server for interactive analysis of cancer and normal gene expression profiles. It analyses the dependency of OS of patients on high and low expression of genes. The calculation of HR was set to the Cox PH model with a confidence interval of 95%, and group cut-off criteria were again set to the median value.

Altered Ion Channels in Tumor and Metastatic States of Breast Cancer
Identification of altered ion channel modules through WGCNA analysis of the differentially expressed ion channels and further selection of ion channels based on the obtained overlaps among the datasets led to the identification of altered ion channels in tumor and metastatic states of breast cancer.

Identification of Altered Ion Channel Modules by Weighted Gene Co-Expression Network Analysis
The soft-thresholding powers chosen to ensure a scale-free network model and the total number of modules obtained using the DynamicTreeCut algorithm of the hierarchical clustering method are provided in Table S7. The heatmap revealed the association of each module with normal, tumor and metastatic phenotypes along with p-values as confidence measures. Assessment of the modules using GS vs. MM plots revealed the intramodular gene connectivity with respect to a particular binary trait. Based on these criteria, selecting one module, each having a better co-relation and p-value resulted in 8 modules for further analysis (Figures 3 and S1).

Altered Ion Channels in Tumor and Metastatic States of Breast Cancer
Identification of altered ion channel modules through WGCNA analysis of the differentially expressed ion channels and further selection of ion channels based on the obtained overlaps among the datasets led to the identification of altered ion channels in tumor and metastatic states of breast cancer.

Identification of Altered Ion Channel Modules by Weighted Gene Co-Expression Network Analysis
The soft-thresholding powers chosen to ensure a scale-free network model and the total number of modules obtained using the DynamicTreeCut algorithm of the hierarchical clustering method are provided in Table S7. The heatmap revealed the association of each module with normal, tumor and metastatic phenotypes along with p-values as confidence measures. Assessment of the modules using GS vs. MM plots revealed the intramodular gene connectivity with respect to a particular binary trait. Based on these criteria, selecting one module, each having a better co-relation and p-value resulted in 8 modules for further analysis (Figures 3 and S1).

Overlapping Ion Channels in GSE42568, GSE52604 and TCGA Datasets
Twenty-two ion channels overlapped in the selected modules corresponding to the three datasets in the HT subgroup (Table 2). From the 22 HT ion channels, 14 were reported previously as significant in tumor initiation and growth in various tumor types. Among those, four ion channels were studied in breast cancer tumor growth. Another set of 22 ion channels were similarly identified in the TM subgroup (Table 3). Of these, 15 were reported to be involved in metastasis in various cancers. This also included 6 ion channels that were reported to play a role in breast cancer progression and metastasis.

Identification of Altered Ion Channels Interacting with EMT-Related Genes
WGCNA of DEGs obtained from the combined gene set dataset and further PPIN analysis of significantly co-expressed modules led to the identification of altered ion channels interacting with EMT-related genes.

Weighted Gene Co-Expression Network Analysis for the Identification of Altered Combined Gene Set Modules
The soft-thresholding powers chosen, the number of modules obtained and the selected significant combined gene set modules are provided in Table S8 and Figure S2. On the basis of the module-trait relationship and GS vs. MM plots, 8 modules with a better co-relation and p-value were chosen for further analysis.

Network Analysis of Altered Ion Channels Co-Expressed with EMT-Related Genes
The modules obtained upon performing WGCNA for the combined gene set represent groups of highly co-expressed genes, although a large number of ion channels were found to be differentially expressed and co-expressed in individual analyses of the expression profiles. Only a subset of ion channels clustered together with EMT-related genes when co-expression analysis of the combined gene set was performed ( Figure S3). The altered ion channels clustered with EMT-related genes were the ion channels that were co-expressed with EMT-related genes and could be involved in the EMT process. The role of TRPM7 in EMT through associations with EGF reported by Davis et al. or the induction of EMT by TGF-β involving CFTR reported by Zhang et al. expounds such notable contributions of ion channels in the EMT program [36,37].
PPI networks through STRING analysis of the combined gene set selected modules revealed the predicted interactions between the genes in the module (Figures 4 and S4).
Further, the PPINs of combined gene set analyzed using NetworkAnalyser resulted in the identification of hub genes in the network. Hub genes in a co-expression network are the highly interconnected nodes in a module. Out of the top 15 selected hub genes in HT of the microarray dataset (GSE42568), ten EMT-related genes were found to interact with ion channels. Of these, nine EMT-related genes were connected to GJA1 (Table 4). Eight genes from HT of RNA-Seq (TCGA) were found to interact with ion channels (Table 5). Similarly, in HM of GSE42568, eight genes were linked to ion channels and CFTR was connected to four of them ( Table 6). HM of the microarray (GSE52604) showed only five genes that interacted with ion channels and RNA-Seq (TCGA) showed seven genes as connected to ion channels (Tables 7 and 8). Six ion channels interacting with EMT-related genes were common to ion channel identified in normal to tumor state of breast cancer and two ion channels were common to the ion channels identified in tumor to metastatic state of breast cancer. Thus, overall, 32 ion channels were found interacting with 26 hub EMT-related genes.
Several ion channels altered in breast cancer patients were found to interact with multiple EMT-related genes. GJA1 and GJB2 belong to the gap junction family. GJA1 was found to interact with the hub EMT-related genes JUN, MYC, FGF2, PTEN, KDR, RHOA, CAV1, ITGB1 and CXCL12. GJB2 was found in interactions with AKT1 and CDH1. TRPC6 and TRPC1 belong to the family of transient receptor cation channels and control the flow of calcium ions. TRPC6 was found to be interacting with JUN and RHOA and TRPC1 was found to be connected to RHOA and CAV1. VDAC1 is a voltage-dependent anion channel and was linked to GAPDH, HSP90AA1 and HSPA4. KCNH2, a voltage-gated potassium channel, was found to be linked with SRC, HSPA4 and HSP90AA1. CFTR is one of the most widely studied ion channels, dysregulation of which has been reported in various pathophysiological conditions [37,38]. It interacted with TP53, GAPDH, BRCA1 and DNMT1. ANO1, a calcium-activated chloride channel, was connected to CDH1 and ERBB2. AQP5 is an aquaporin involved in the movement of water across cell membranes [39]. It interacted with GAPDH and CDH1. CLIC1, a chloride channel, was found interacting with AKT1 and GAPDH.

Correlation of Identified Ion Channels with EMT Scores
The EMT scores obtained using GS76, MLR and KS methods for each sample in RNA-Seq and microarray datasets are provided in Tables S9-S11. The samples with a negative GS76 score, a positive KS score and a higher MLR score can be interpreted as mesenchymal samples. Subsequently, the samples with a positive GS76 score, a negative KS score and a lower MLR score could be an epithelial sample. The correlation between the scores relative to the specific method in each dataset was calculated and is depicted in Figure 5 and provided in Table 9. KS and MLR scores correlated positively with each other and GS76 scores correlated negatively with both KS and MLR scores.   Further, the correlations between the obtained EMT scores from each method and the expression levels of the identified ion channels in tumor, metastatic and ion channels interacting with EMT-related genes were computed (Table S12).
The ion channels CLIC2, GJA4, HTR3C, CLCN6, SCN3A, ANO3, LRRC8C and GJA5 were found to be negatively correlated to GS76 and positively correlated to both KS and MLR scores in all three datasets, indicating that these ion channels might have a higher probability of contributing to acquiring and/or stabilizing a mesenchymal phenotype.
Nevertheless, all ion channels other than SCNN1G, CLCN3, KCNJ10 and KCNH2 had a similar correlation trend in the expression values corresponding to one/more datasets. Figure 6 depicts the correlation of each of the identified ion channels with the EMT scores. The bubble plot represents the variation in the correlation of the genes with the three scoring methods. It landscapes the fact that the ion channels that may have a higher probability of contributing to the mesenchymal phenotype, correlating positively with the MLR and KS methods and negatively with the GS76 method. Further, the correlations between the obtained EMT scores from each method and the expression levels of the identified ion channels in tumor, metastatic and ion channels interacting with EMT-related genes were computed (Table S12).

Pathway Map of Identified Ion Channels and Their Associations with EMT
Generation of pathway maps depicting potential ion channels in several cellular processes may aid in understanding their significance in various aspects of tumorigenesis in breast cancer patients. Figure 7 summarizes the overall identified ion channels and the events likely to occur in breast cancer patients when dysregulated ion channels function together with EMT-related genes. The pathway map consists of 66 molecules and 58 reactions. Several ion channels were differentially expressed and were found to have significant roles in key cellular processes.

Pathway Map of Identified Ion Channels and Their Associations with EMT
Generation of pathway maps depicting potential ion channels in several cellular processes may aid in understanding their significance in various aspects of tumorigenesis in breast cancer patients. Figure 7 summarizes the overall identified ion channels and the events likely to occur in breast cancer patients when dysregulated ion channels function together with EMT-related genes. The pathway map consists of 66 molecules and 58 reactions. Several ion channels were differentially expressed and were found to have significant roles in key cellular processes. Figure 7. Depiction of events that may occur in breast cancer patients upon dysregulation of ion channels with EMT-related genes. The pathway map was generated using PathVisio (v3.3.0). Cancer phenotypes may appear due to alterations in important cellular processes, such as calcium signaling, insulin secretion, adipocyte metabolism, nitric oxide signaling and glutamatergic signaling.

Association of Identified Putative Ion Channels in Survival of Breast Cancer Patients
The ion channels with HR > 1 and a p(HR) < 0.05 were selected as significant ion channels in survival of patients with breast cancer (Figure 8). Among the ion channels identified in association with EMT-related genes and tumor/ metastatic states of breast cancer, five ion channels-CACNA1B (HR: 1.14,  . Depiction of events that may occur in breast cancer patients upon dysregulation of ion channels with EMT-related genes. The pathway map was generated using PathVisio (v3.3.0). Cancer phenotypes may appear due to alterations in important cellular processes, such as calcium signaling, insulin secretion, adipocyte metabolism, nitric oxide signaling and glutamatergic signaling.

Association of Identified Putative Ion Channels in Survival of Breast Cancer Patients
The ion channels with HR > 1 and a p(HR) < 0.05 were selected as significant ion channels in survival of patients with breast cancer (Figure 8). Among the ion channels identified in association with EMT-related genes and tumor/ metastatic states of breast cancer, five ion channels-CACNA1B (HR: 1.14, p

Discussion
A deeper understanding of the role of ion channels in cells expressing cancerous phenotypes needs to be elucidated. Blockade of ion channels have been demonstrated to influence various pathophysiological conditions [40], making ion channels potential biomarkers in cancer diagnosis and therapeutics. This study focused on the identification of putative ion channels in tumor and metastatic states and their association with EMT programs of breast cancer through various computational analyses.
In the present study, analysis of transcriptomic data belonging to different platforms revealed the presence of several altered ion channels and EMT-related genes in patients with breast cancer. WGCNA of the altered ion channels led to the identification of significantly relevant co-expressed ion channels that were clustered into modules. Each of the selected modules consisted of a distinct set of ion channels. Although a few overlapping ion channels were noticed in the modules representing normal to tumor state and tumor to metastatic state of breast cancer. The normal to metastatic state modules appeared to be populated with a manifold of co-expressed ion channels as compared to the normal to tumor state modules. This could be an indication that as cancer progresses from normal

Discussion
A deeper understanding of the role of ion channels in cells expressing cancerous phenotypes needs to be elucidated. Blockade of ion channels have been demonstrated to influence various pathophysiological conditions [40], making ion channels potential biomarkers in cancer diagnosis and therapeutics. This study focused on the identification of putative ion channels in tumor and metastatic states and their association with EMT programs of breast cancer through various computational analyses.
In the present study, analysis of transcriptomic data belonging to different platforms revealed the presence of several altered ion channels and EMT-related genes in patients with breast cancer. WGCNA of the altered ion channels led to the identification of significantly relevant co-expressed ion channels that were clustered into modules. Each of the selected modules consisted of a distinct set of ion channels. Although a few overlapping ion channels were noticed in the modules representing normal to tumor state and tumor to metastatic state of breast cancer. The normal to metastatic state modules appeared to be populated with a manifold of co-expressed ion channels as compared to the normal to tumor state modules. This could be an indication that as cancer progresses from normal to tumor and tumor to metastatic, the number of ion channels being deregulated aggravates.

Ion Channels Identified as Putative Ion Channels in the Tumor State of Breast Cancer
• HTR3C, a ligand-gated ion channel, is one of the receptors for serotonin, a hormone that functions as a neurotransmitter and mitogen. Activation of this receptor causes fast, depolarizing responses. Additionally, SCN3A, a voltage gated sodium channel, is known to maintain depolarization in the enterochromaffin cells, which results in the regulation of serotonin release [57]. Serotonin functions as a tumor-suppressant in non-transformed breast cells and early-stage breast cancers. During tumor progression, cells acquire genetic or epigenetic alterations in serotonin signaling. This makes them resistant to suppressive actions of serotonin and favors tumor-promoting actions [58]. • CLCN6 is a voltage-dependent chloride channel and has a role in regulating blood pressure levels and hypertension [59]. cMyc a protooncogene responsible for cell proliferation in various cancers, transcriptionally regulates GRK4 protein that was reported to be overexpressed in breast cancer tissues [60]. GRK4 has been demonstrated to be associated with an increased risk of hypertension, indicating hypertension as an important factor in breast cancer [60]. • GLRB, a glycine receptor, is a ligand-gated ion channel that mediates the inhibitory effects of glycine. Vascular endothelial growth factor (VEGF) has a crucial role in cancer progression as it promotes the formation of new blood vessels. Activation of VEGF receptor results in activation of phospholipase C-gamma and increases intracellular Ca 2+ concentration. VEGF-induced cell proliferation is dependent on intracellular Ca 2+ concentration. Hyperpolarization of the cell membrane due to glycine-gated chloride channels blocks the influx of Ca 2+ , thereby minimizing VEGF-mediated signaling. Thus, changes in the functioning of GLRB may promote tumor growth [61]. • Anoctamins (ANOs) are Ca 2+ activated chloride channels. The activation of receptors of growth hormone signaling takes place through RAS-RAF-ERK, PI3K-AKT and DAG-IP3 pathways. It is known that anoctamin-controlled calcium channels are relevant for the activation of ERK-1,2. Additionally, the rise in intracellular Ca 2+ together with activation of RAS-RAF/ERK pathway is a major aspect of cell proliferation [62]. ANO1 has been widely studied in various tumors in this respect. However, very few studies have reported the involvement of other members of anoctamins. ANO3, ANO5 and ANO6 were identified in this study. • LRRC8C is a volume-regulated anion channel (VRAC). VRACs are activated on cell swelling and play a critical role in cell volume regulation [63]. The Fad158 gene, which has a crucial role in adipocyte differentiation, belongs to the LRRC8 family and is also known as LRRC8C [64]. As adipocytes are the major components of breast tissue, abnormal adipocyte metabolism leads to accumulation of tumor-supporting cells and other effects, such as insulin resistance, dyslipidemia and oxidative stress [65]. These effects may further lead to the aggressive nature of the tumor. Additionally, VRACs are known to mediate cellular uptake of drugs, such as cisplatin and carboplatin, which are widely used in the treatment of cancer [66]. However, LRRC8C was not reported to have a direct effect on the uptake of these drugs but was stated to have a possible role in mediating the action of LRRC8A [66].

Ion Channels Identified as Putative Ion Channels in the Metastatic State of Breast Cancer
• TRPV3 is not a well-understood ion channel and can be activated by temperature.
Although it has been reported as overexpressed in non-small cell lung cancer [84]. • GABRG3 belongs to the GABA-A receptor gene family of hetero-pentameric ligandgated ion channels. GABA-A receptor GABRP is required for maintaining basal-like cytokeratin expression, ERK1/2 phosphorylation and pro-migratory phenotype of breast cancer cells [85]. GABA-B receptors promote metastasis by enhancing ERK phosphorylation and thus activating metalloproteins that enable tumor cells to penetrate the basement membrane [86]. These studies indicate the importance of GABA receptors in the metastatic progression of breast cancer. • KCNT1 is a Na + activated K + channel with diverse functions, including insulin secretion and cell volume regulation. With the help of insulin receptors, insulin regulates endothelial cell migration, proliferation and production of VEGF. It also activates PI3K/Akt signaling that promotes nitric oxide (NO) release. NO increases endothelial survival, migration, proliferation and vascular permeability [87]. • KCNN1 is a Ca 2+ dependent potassium channel involved in regulating cell volume. Increased permeability of K + due to activation of these channels results in membrane hyperpolarization. This enhances Ca 2+ entry and thus contributes to a decrease in the regulatory volume of a cell [88]. • CLCNKB is a voltage-gated chloride channel. Chloride channels participate in cell volume regulation, membrane potential stabilization, signal transduction and transepithelial transport. CLCNKB has been widely studied in Bartter syndrome. Its role in breast cancer is not well established. These channels aid in maintaining intracellular ion concentration lower than extracellular ion concentration, preventing osmotic cell swelling. Altered expression in these channels may result in osmotic stress contributing to cell migration [89]. • GJA5 is a member of the connexin family. Connexins are involved in the formation of heterologous gap junctions between tumor and endothelial cells, which facilitate intravasation and extravasation [90]. Various gap junctions have been reported in this regard. Very few studies have reported GJA5 [91]. • GRIK5 belongs to the glutamate-gated ionic family. Cells with breast cancer phenotype secrete high levels of glutamate and metastasize to bones. These cells with excess glutamate result in cancer-induced bone pain, which is a significant co-morbidity in advanced stage breast cancer patients [92].
A subset of ion channels, including GJA1, TRPC6, VDAC1 and AQP5, clustered together with hub EMT-related genes when PPI network analysis of the combined gene set was performed.

Ion Channels Interacting with Multiple Hub EMT-Related Genes in Breast Cancer
• GJA1 has roles in the cAMP pathway, Wnt signaling pathway and activation of p38. It was also reported as a tumor suppressor and a potent molecule in the acceleration of cancer progression [93]. • TRPC6 is present ubiquitously in human tissues and is involved in Ca 2+ dependent pathways. It is also known to be upregulated in various other disease conditions [94]. • TRPC1, another ion channel involved in modulation of Ca 2+ , is an established biomarker in certain cancers [95]. • VDAC1 is a multifunctional channel involved in controlling the communications between mitochondria and the rest of the cell [96].
• ANO1 channel is involved in cell proliferation, survival, migration, contraction, secretion and neuronal excitation [97]. • AQP5, an aquaporin, is involved in the movement of water across cellular membranes. This process is one of the most important phenomena resulting in cell movement, cellular viscosity and signal transductions [39]. • CLIC1 was identified in maintenance of cell volume, ion homeostasis, trans-epithelial transport, pH regulation and cell cycle regulation [98]. • GJB2 regulates cell migration and colonization and thus aggressive phenotype in breast cancer [99]. • KCNH2 was reported to be involved in cell proliferation and cell migration processes implicating MAP kinase and c-fos pathways through a cell line study [100]. • CFTR, an anion channel, is known for its role in ion and acid-base homeostasis. It was also reported as a tumor suppressor [101].
An increasing number of studies are currently aiming to gain insights into the EMT program in cancer, which includes identifying and understanding the role of specific ion channels in the induction and maintenance of various aspects of EMT programs. The outcomes of these studies provide a preliminary perception that ion channels as therapeutic targets may be useful to control EMT in cancer. The involvement of putative ion channels in various cellular processes further substantiates its potential as a target. Thus, ion channels correlating with EMT-related genes identified through this study may aid in modulating several processes involved in tumor growth and progression and could act as promising targets. Kaplan-Meier plots and Cox regression analysis determined the prognostically significant ion channels in breast cancer among the putative ion channels. Overall, this study could identify various ion channels in breast cancer. Nevertheless, due to the data mining approach, the outcomes of the study may appear in the form of overfitting or underfitting. Thus, further theoretical and experimental studies need to be carried out to validate the obtained findings. Together, our analysis is a relevant data source for ion channels in breast cancer and may help in its management.

Conclusions
A systems biology-based approach was used to identify putative ion channels in tumor growth and metastasis and their correlation with EMT-related genes through analysis of RNA-Seq and microarray-based expression profiles of patients with breast cancer. Amid the identified ion channels, also present were ion channels already established as prognostic markers in breast cancer. Functional annotations of the altered ion channels led to the identification of processes contributing to cell proliferation, cell migration and cell volume regulation in breast cancer. However, detailed mechanisms underlying the possible effects of the identified ion channels and their associations with EMT need to be further characterized in vitro and in vivo.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/cancers14061444/s1, Table S1: Clinical Information of breast cancer patients from TCGA and GTEx datasets, Table S2: Clinical information of breast cancer patients in GSE52604 dataset, Table S3: Clinical information of breast cancer patients in GSE42568 dataset, Table S4: non-redundant list of differentially expressed ion channels in HM, NT and TM corresponding to patients with breast cancer, Table S5: non-redundant list of differentially expressed EMT-related genes in HM, NT and TM corresponding to patients with breast cancer, Table S6: Differential expression of ion channels overlapped in EMT gene-set, Table S7: Parameters used for WGCNA analysis of ion channel gene-set in HT, TM and HM, Table S8: Parameters used for WGCNA analysis of ion channels integrated with EMT gene-set in HT, TM and HM, Table S9: A list of samples present in RNA-Seq dataset and their corresponding EMT score estimated using GS76, MLR and KS EMT scoring methods, Table S10: A list of samples present in GSE42568 dataset and their corresponding EMT score estimated using GS76, MLR and KS EMT scoring methods, Table S11: List of samples present in GSE52604 dataset and their corresponding EMT score estimated using GS76, MLR and KS EMT scoring methods, Table S12: Correlation of expression values of ion channels identified as interacting with EMT and ICs identified in tumor and metastatic states belonging to RNA-Seq, GSE42568 and GSE52604 with GS76, MLR and KS EMT scoring methods, Figure