Transcriptomic analysis of calcium remodeling in colorectal cancer

International Journal of Molecular Sciences Article Transcriptomic Analysis of Calcium Remodeling in Colorectal Cancer Enrique Pérez-Riesgo 1,2, Lucía G. Gutiérrez 1,2, Daniel Ubierna 1, Alberto Acedo 3, Mary P. Moyer 4, Lucía Núñez 1,2 and Carlos Villalobos 1,* 1Institute of Molecular Biology and Genetics (IBGM), National Research Council (CSIC), 47003 Valladolid, Spain; [email protected] (E.P.-R.); gonzalez.gutierr[email protected] (L.G.G.); [email protected] (D.U.); [email protected] (L.N.) 2Department of Biochemistry and Molecular Biology and Physiology, University of Valladolid, 47005 Valladolid, Spain 3AC-Gen Reading Life, 47011 Valladolid, Spain; [email protected] 4INCELL Corporation, San Antonio, TX 78249, USA; [email protected] *Correspondence: [email protected]; Tel.: +34-983-184-821 Academic Editors: Peter J. K. Kuppen and Dario Marchetti Received: 8 March 2017; Accepted: 13 April 2017; Published: 27 April 2017 Abstract: Colorectal cancer (CRC) cells undergo the remodeling of intracellular Ca 2+ homeostasis, which contributes to cancer hallmarks such as enhanced proliferation, invasion and survival. Ca 2+ remodeling includes critical changes in store-operated Ca 2+ entry (SOCE) and Ca 2+ store content. Some changes have been investigated at the molecular level. However, since nearly 100 genes are involved in intracellular Ca 2+ transport, a comprehensive view of Ca 2+ remodeling in CRC is lacking. We have used Next Generation Sequencing (NGS) to investigate differences in expression of 77 selected gene transcripts involved in intracellular Ca 2+ transport in CRC. To this end, mRNA from normal human colonic NCM460 cells and human colon cancer HT29 cells was isolated and used as a template for transcriptomic sequencing and expression analysis using Ion Torrent technology. After data transformation and filtering, exploratory analysis revealed that both cell types were well segregated. In addition, differential gene expression using R and bioconductor packages show significant differences in expression of selected voltage-operated Ca 2+ channels and store-operated Ca 2+ entry players, transient receptor potential (TRP) channels, Ca 2+ release channels, Ca 2+ pumps, Na + /Ca 2+ exchanger isoforms and genes involved in mitochondrial Ca 2+ transport. These data provide the first comprehensive transcriptomic analysis of Ca2+ remodeling in CRC. Keywords: colorectal cancer; Ca2+ remodeling; transcriptomics; RNA-sequencing 1. Introduction Ca 2+ signaling is involved in the control of a large series of cell and physiological functions in health and disease from exocytosis and muscle contraction to gene expression, cell proliferation, migration and cell death [ 1 ]. Over the last few years, evidence has been accumulating that dis-homeostasis of intracellular Ca 2+ may be involved in many different disorders including cancer [2,3]. In support of this view, several authors have reported that remodeling of Ca2+ signaling in different forms of cancer contributes to cancer hallmarks, including exaggerated cell proliferation, acquisition of cell migration and invasion capabilities, and enhanced resistance to cell death [ 2 , 3 ]. We have recently compared intracellular Ca 2+ homeostasis in normal human colonic cells and human colon carcinoma cells [ 4 ]. We found critical differences that may contribute to several of the above cancer hallmarks. The most salient characteristics of Ca2+ remodeling in colon cancer cells relative to their normal counterparts are the dramatic rise in store-operated Ca 2+ entry (SOCE) and the partial Int. J. Mol. Sci. 2017,18, 922; doi:10.3390/ijms18050922 www.mdpi.com/journal/ijms Int. J. Mol. Sci. 2017,18, 922 2 of 23 depletion of intracellular Ca 2+ stores [ 4 ]. Enhanced entry of Ca 2+ in colon cancer cells correlates with enhanced cell proliferation and all SOCE antagonists tested and/or targeting molecular players involved in this pathway result in inhibition of cell proliferation and cell migration [ 4 – 6 ]. In addition, the partial depletion of Ca 2+ stores provides resistance to cell death, probably because of impaired mitochondrial Ca2+ overload and apoptosis [4–6]. SOCE, first reported by James W Putney in 1986 [ 7 , 8 ] is a ubiquitous Ca 2+ entry pathway activated after agonist-induced release of Ca 2+ from intracellular stores at the endoplasmic reticulum (ER). This is due to the opening of intracellular Ca 2+ release channels gated by phospholipase C-dependent synthesis of inositol trisphosphate IP 3 receptors (IP 3 Rs) and/or Ca 2+ -induced Ca 2+ release mediated by ryanodine receptors (RYRs). At the molecular level, it is well established that Stromal Interaction Molecules (STIM) sense ER Ca 2+ levels. Emptying of Ca 2+ stores promotes STIM oligomerization at the ER endomembranes [ 9 ] and their interaction with cation channels in the plasma membrane of the Orai1 [ 10 ] and transient receptor potential (TRP) families of Ca 2+ channels [ 11 ]. Activation of SOCE not only permits the refilling of Ca 2+ stores but also promotes a sustained increase in intracellular Ca 2+ concentration that is critical for the sustained activation of effector proteins including, for instance, calcineurin and the ensuing activation of the nuclear factor of activated T cells [10]. After stimulation, the increased cytosolic free Ca 2+ concentrations ([Ca 2+ ] cyt ) must return to resting levels in the low nM range. This is achieved by systems that transport Ca 2+ back from cytosol to the extracellular space or into the ER against a large electrochemical gradient for Ca 2+ . This is carried out by different adenosine triphosphate (ATP)ase Ca 2+ pumps including plasma membrane Ca 2+ ATPases (PMCAs) [ 12 ], sarcoplasmic and/or endoplasmic reticulum Ca 2+ ATPases (SERCAs) [ 13 ] and the secretory pathway Ca 2+ ATPases (SPCAs) [ 14 ]. An additional cotransporter, the Na + /Ca 2+ exchanger, uses the passive transport of Na + into the cell to extrude Ca 2+ out of the cell [ 15 ]. Finally, as has become increasingly evident in recent years, mitochondria contribute to clear Ca 2+ loads by means of the mitochondrial Ca 2+ uniporter (MCU) and regulatory proteins [ 16 ]. In this case, Ca 2+ is transported down a huge electromotive force, the mitochondrial potential ( ∆Ψ ) of about − 180 mV, negative inside the mitochondrial matrix, that enables Ca 2+ influx into mitochondria [ 17 , 18 ]. This potential is particularly high in cancer cells due to the Warburg effect, the metabolic signature of cancer cells [ 19 ]. Interestingly, mitochondria are also involved in control of SOCE in different cell types [ 20 , 21 ], including colorectal cancer cells [ 18 , 22 ] where, store-operated Ca 2+ channels may undergo slow, Ca 2+ -dependent inactivation unless this process is prevented by Ca 2+ removal by surrounding mitochondria [20,21]. Our recent analysis of Ca 2+ remodeling in colon cancer cells revealed increased expression of several molecular players involved in SOCE [ 4 – 6 ]. However, Ca 2+ remodeling is complex and may involve changes in expression and/or activity of hundreds of genes related to Ca 2+ signaling [ 23 ], with nearly 80 of them directly involved in Ca 2+ transport across biological membranes. By taking advantage of Next Generation Sequencing technologies, particularly RNA-sequencing (RNA-seq) by Ion Torrent methodologies, we were able to study in detail intracellular Ca 2+ homeostasis in NCM460 cells reflecting normal human colonic cells [ 24 , 25 ] and HT29 cells [ 26 ] representing human colorectal cancer cells, to investigate, for the first time, the transcriptomics of Ca 2+ remodeling in CRC. To this end, mRNA was isolated from four independent cultures of NCM460 and four independent cultures of HT29 cells to be used as templates for transcriptomic analysis of the selected 77 genes shown in Table 1. The gene transcripts studied include all ten voltage-operated Ca 2+ channels (VOCCs); all seven well-established players involved in SOCE; all 27 TRP channels; all six Ca 2+ release channels—including IP 3 receptors and ryanodine receptor isoforms—all nine Ca 2+ pumps including PMCAs, SERCAs and SPCAs; the three Na + /Ca 2+ exchanger isoforms; all known genes involved in mitochondrial Ca 2+ transport; and other related proteins. Clustering analysis segregated data in two different cell populations and genes were segregated into six independent gene families according to Int. J. Mol. Sci. 2017,18, 922 3 of 23 expression behavior. Overall, we found that 30 genes are expressed differentially in CRC cells relative to normal cells, thus providing a first comprehensive view of Ca2+ remodeling in CRC. Table 1. Genes involved in calcium transport analyzed in this study pooled by gene families. Gene Group Protein Name (Gene Name) Voltage Operated Calcium Channels Cav1.1 (CACNA1S); Cav1.2 (CACNA1C); Cav1.3 (CACNA1D); Cav1.4 (CACNA1F); Cav2.1 (CACNA1A); Cav2.2 (CACNA1B); Cav2.3 (CACNA1E); Cav3.1 (CACNA1G); Cav3.2 (CACNA1H); Cav3.3 (CACNA1I) Store-Operated Calcium Entry Player Orai1; Orai2; Orai3; STIM1; STIM2; CRACR2A (EFCAB4B); MS4A12 TRP Channels TRPC1; TRPC3; TRPC4; TRPC5; TRPC6; TRPC7; TRPV1; TRPV2; TRPV3; TRPV4; TRPV5; TRPV6; TRPM1; TRPM2; TRPM3; TRPM4; TRPM5; TRPM6; TRPM7; TRPM8; TRPA1; TRPML1 (MCOLN1); TRPML2 (MCOLN2); TRPML3 (MCOLN3); TRPP1 (PKD2); TRPP2 (PKD2L1); TRPP3 (PKD2L2) Calcium Release Channels IP3R1 (ITPR1); IP3R3 (ITPR2); IP3R3 (ITPR3); RYR1; RYR2; RYR3 Calcium Pumps PMCA1 (ATP2B1); PMCA2 (ATP2B2); PMCA3 (ATP2B3); PMCA4 (ATP2B4); SERCA1 (ATP2A1); SERCA2 (ATP2A2); SERCA3 (ATP2A3); SPCA1 (ATP2C1); SPCA2 (ATP2C2) Sodium Calcium Exchangers NCX1 (SLC8A1); NCX2 (SLC8A2); NCX3 (SLC8A3) Mitochondrial Calcium Transport Proteins MCU; MICU1; MICU2 (EFHA1); MICU3 (EFHA2); MCUR1 (CCDC90A); EMRE (C22orf32); MCUb (CCDC109B); VDAC1; VDAC2; VDAC3 Other Proteins Bcl-2 (BCL2); Calsequestrin 1 (CASQ1); Calsequestrin 2 (CASQ2); PICALM; Phospholamban (PLN) Cav: Ca 2+ channel; STIM: Stromal Interaction Molecules; TRPC: canonical TRP channels; TRPV: vaniloid family of transient receptor potential channels; TRPM: melastatin family of transient receptor potential channels; TRPP: polycystine family of transient receptor potential channels; IP 3 R: inositol trisphosphate receptors; RYR: ryanodine receptors; PMCA: plasma membrane Ca 2+ ATPases; SERCA: sarcoplasmic and/or endoplasmic reticulum Ca 2+ ATPases; SPCA: secretory pathway Ca 2+ ATPases; NCX: sodium-calcium exchanger; MCU: mitochondrial Ca 2+ uniporter; EMRE: MCU regulator; VDAC: voltage-dependent anion channels; PICALM: phosphatidylinositol binding clathrin assembly protein. 2. Results and Discussion 2.1. Data Set Transformation and Filtering mRNA from four independent cultures of normal (NCM460) human colonic cells and four independent cultures of human colorectal (HT29) cells were isolated and used as templates for transcriptomic analysis of 77 genes involved in intracellular Ca 2+ transport (Table 1) as detailed in the Methods section. RNA-seq data sets consist of aligned reads from gene fragmentations against a reference genome. Accordingly, the longer a gene is, the higher the amount of fragments that come from it. In addition, the longer the gene is, the more likely it is that one of the reads aligns against this gene. As we are concerned with analyzing not only differential expression between normal and tumor cells, but also with differential expression among genes of the same family, raw gene data was normalized for the corresponding gene length (l) and the size of the library samples (n), thus making expression values among different genes readily comparable. As expression values after normalization are very small, they were multiplied by 10 9 to obtain expression values as Reads per Kilobase per Millions of mapped reads (RPKM) that can be obtained using the following expression (1). RPKM =reads l·n·109(1) Another important issue relative to RNA-seq data sets is that they contain usually a large amount of low and disparate values. Thus, it is very important to be careful with the scale employed. If we Int. J. Mol. Sci. 2017,18, 922 4 of 23 use a raw data set, almost all data will plot on low values and only a few on higher values, thus occluding large parts of data (Figure 1A,B). This problem could affect techniques of Exploratory Data Analysis such as Principal Components Analysis or Cluster Hierarchical Analysis, since highly expressed genes could dominate them. To solve this issue, raw data has to be transformed numerically using, for example, Box-Cox Transformations or Logarithmic Transformations. In the present study, we used a data transformation similar to logarithmic transformation in base 2. Transformation shrinks data variation across samples. However, low-expression genes tend to over-dominate when this transformation is used. This is due, on one hand, to Poisson noise, related to low values of reads and, on the other hand, to the fact that logarithm transformation amplifies differences among low values. Thus, low-expression genes will show relatively larger differences among samples (Figure 1C,D). To solve this problem, we have used Regularized Logarithm Transformation (rlog). In this case, the outcomes of highly expressed genes are similar to outcomes obtained by logarithmic transformation in base two. In addition, differences among low values essentially disappear (Figure 1E,F). Int. J. Mol. Sci. 2017, 18, 922 4 of 22 occluding large parts of data (Figure 1A,B). This problem could affect techniques of Exploratory Data Analysis such as Principal Components Analysis or Cluster Hierarchical Analysis, since highly expressed genes could dominate them. To solve this issue, raw data has to be transformed numerically using, for example, Box-Cox Transformations or Logarithmic Transformations. In the present study, we used a data transformation similar to logarithmic transformation in base 2. Transformation shrinks data variation across samples. However, low-expression genes tend to overdominate when this transformation is used. This is due, on one hand, to Poisson noise, related to low values of reads and, on the other hand, to the fact that logarithm transformation amplifies differences among low values. Thus, low-expression genes will show relatively larger differences among samples (Figure 1C,D). To solve this problem, we have used Regularized Logarithm Transformation (rlog). In this case, the outcomes of highly expressed genes are similar to outcomes obtained by logarithmic transformation in base two. In addition, differences among low values essentially disappear (Figure 1E,F). Figure 1. Gene expression level data transformation. Raw data set values (A,B). Log2 transformation (C,D); rlog transformation (E,F). RNA-seq: RNA-sequencing. Another key point of RNA-seq data analysis is filtering in order to minimize hypothesis testing. We have to take into account that RNA-seq data sets contain a large amount of values equal to zero, since not all genes are expressed by every kind of cell, and that is very important relative to filtering. Thus, those genes that may not play any role or show no relationship to phenotype may be removed from the data expression set. Nevertheless, we need to be careful in this step to avoid removing genes which marginally do not show any kind of activity in a cell but may act jointly with other genes. Thus, we have considered taking into account the criteria set up by Hackett [27] and Seyednasrollah [28], who suggested filtering out those genes that are not readily expressed in any of the samples and whose median expression profile is lower than 0.125 RPKM [29]. Data transformed in this manner were filtered in order to remove those genes considered as non-expressed, whose median expression profile is lower than 0.125 RPKM (Figure 2). Interestingly, three types of gene populations according Figure 1. Gene expression level data transformation. Raw data set values ( A , B ). Log2 transformation (C,D); rlog transformation (E,F). RNA-seq: RNA-sequencing. Another key point of RNA-seq data analysis is filtering in order to minimize hypothesis testing. We have to take into account that RNA-seq data sets contain a large amount of values equal to zero, since not all genes are expressed by every kind of cell, and that is very important relative to filtering. Thus, those genes that may not play any role or show no relationship to phenotype may be removed from the data expression set. Nevertheless, we need to be careful in this step to avoid removing genes which marginally do not show any kind of activity in a cell but may act jointly with other genes. Thus, we have considered taking into account the criteria set up by Hackett [ 27 ] and Seyednasrollah [ 28 ], who suggested filtering out those genes that are not readily expressed in any of the samples and whose median expression profile is lower than 0.125 RPKM [ 29 ]. Data transformed in this manner Int. J. Mol. Sci. 2017,18, 922 5 of 23 were filtered in order to remove those genes considered as non-expressed, whose median expression profile is lower than 0.125 RPKM (Figure 2). Interestingly, three types of gene populations according to expression values emerge after filtering: genes expressed at low, intermediate and high expression levels (Figure 2). Int. J. Mol. Sci. 2017, 18, 922 5 of 22 to expression values emerge after filtering: genes expressed at low, intermediate and high expression levels (Figure 2). Figure 2. Gene expression level filtering. Histograms of data before (A) and after (B) the data set filtration process. Filtering was carried out by removing genes with a profile expression lower than 0.125 Reads per Kilobase per Millions of mapped reads (RPKM). Data are shown as rlog transformed. 2.2. Exploratory Data Analysis: Hierarchical Clustering After data transformation and filtration, two different exploratory data analyses were carried out; although they are not hypothesis testing, they are useful for generating hypotheses. The first one, known as Hierarchical Cluster Analysis, is carried out both in samples and genes. This analysis, using Euclidean distance as dissimilarity measurement, allows the breaking up of samples or genes into groups of members sharing common characteristics. It also establishes which genes are co-regulated. The second exploratory data analysis is known as Principal Components Analysis, which is very useful in order to reduce the dimension of the data or to find patterns among genes or samples, thus getting new variables which are a linear combination of the original ones. Hierarchical Clustering Analysis forms observation clusters whose members share common characteristics. For example, it is possible to detect different populations, such as healthy and tumor phenotypes, thereby clustering each sample into a healthy or tumor group. Furthermore, if the existence of different groups is known, it could be possible to identify a rule to sort the observations. For example, we can sort samples with Discriminant Analysis [30]. In the present study, we clustered both samples and genes. Thus, if we classify samples, we can discover the different populations they came from and, if we classify genes, we can use the information to carry out, for example, an Enrichment Gene Set Analysis that identifies which gene clusters are related to the tumoral phenotype. However, this analysis was not undertaken here because of the low number of genes studied and for the sake of brevity. We first carried out the Hierarchical Cluster Analysis considering the eight samples used as observations and genes as variables. The dissimilarity measurement chosen is Euclidean distance, and the analysis is shown as a heatmap with dendograms in both axes of Figure 3. According to the map, it is clear that two groups have been formed, where the first one contains all four samples belonging to the normal, healthy phenotype (A, B, C, D), and the second one contains all four samples corresponding to the tumor phenotype (E, F ,G ,H) (Figure 3). The distance among samples belonging to the same phenotype is much smaller than the distance among samples belonging to different phenotypes. Each value represented in the heatmap corresponds with the distance between each sample couple. A second Hierarchical Clustering Analysis was carried out in a similar way, except that now the 77 genes are considered as observations and the eight samples as variables. The result of the analysis is shown in Figure 4. The Euclidian distances map shows that the 77 genes in the pool break into six different groups. All the genes in each group behave similarly in the sense that, regardless of whether they are high or low expressed in normal and tumor cells, expression levels are similar for all genes within each group and different from expression levels of the other gene families according to the Euclidian distances shown in Figure 4. Figure 2. Gene expression level filtering. Histograms of data before ( A ) and after ( B ) the data set filtration process. Filtering was carried out by removing genes with a profile expression lower than 0.125 Reads per Kilobase per Millions of mapped reads (RPKM). Data are shown as rlog transformed. 2.2. Exploratory Data Analysis: Hierarchical Clustering After data transformation and filtration, two different exploratory data analyses were carried out; although they are not hypothesis testing, they are useful for generating hypotheses. The first one, known as Hierarchical Cluster Analysis, is carried out both in samples and genes. This analysis, using Euclidean distance as dissimilarity measurement, allows the breaking up of samples or genes into groups of members sharing common characteristics. It also establishes which genes are co-regulated. The second exploratory data analysis is known as Principal Components Analysis, which is very useful in order to reduce the dimension of the data or to find patterns among genes or samples, thus getting new variables which are a linear combination of the original ones. Hierarchical Clustering Analysis forms observation clusters whose members share common characteristics. For example, it is possible to detect different populations, such as healthy and tumor phenotypes, thereby clustering each sample into a healthy or tumor group. Furthermore, if the existence of different groups is known, it could be possible to identify a rule to sort the observations. For example, we can sort samples with Discriminant Analysis [ 30 ]. In the present study, we clustered both samples and genes. Thus, if we classify samples, we can discover the different populations they came from and, if we classify genes, we can use the information to carry out, for example, an Enrichment Gene Set Analysis that identifies which gene clusters are related to the tumoral phenotype. However, this analysis was not undertaken here because of the low number of genes studied and for the sake of brevity. We first carried out the Hierarchical Cluster Analysis considering the eight samples used as observations and genes as variables. The dissimilarity measurement chosen is Euclidean distance, and the analysis is shown as a heatmap with dendograms in both axes of Figure 3. According to the map, it is clear that two groups have been formed, where the first one contains all four samples belonging to the normal, healthy phenotype (A, B, C, D), and the second one contains all four samples corresponding to the tumor phenotype (E, F ,G ,H) (Figure 3). The distance among samples belonging to the same phenotype is much smaller than the distance among samples belonging to different phenotypes. Each value represented in the heatmap corresponds with the distance between each sample couple. A second Hierarchical Clustering Analysis was carried out in a similar way, except that now the 77 genes are considered as observations and the eight samples as variables. The result of the analysis Int. J. Mol. Sci. 2017,18, 922 6 of 23 is shown in Figure 4. The Euclidian distances map shows that the 77 genes in the pool break into six different groups. All the genes in each group behave similarly in the sense that, regardless of whether they are high or low expressed in normal and tumor cells, expression levels are similar for all genes within each group and different from expression levels of the other gene families according to the Euclidian distances shown in Figure 4. Int. J. Mol. Sci. 2017, 18, 922 6 of 22 Figure 3. Heatmap Hierarchical Cluster of samples with healthy phenotype (A–D) and tumor phenotype (E–H). Pseudocolor scale shows hierarchical distance from minimum (0, blue) to maximum (6, red). Figure 4. Hierarchical clustering of genes. Genes have been divided into six different groups according to their Euclidean distances relative to their expression values (in RPKM). 2.3. Gene Correlations Next, we questioned possible correlations between gene couples. In this way, correlation values for each gene couple were plotted in the heatmap shown in Figure 5. Thus, in the case that two genes are positively and strongly correlated, if the expression of one of them is low in tumor samples, the expression of the other one will also be low in the normal sample. Conversely, if the pair of genes is negatively and strongly correlated and the expression of one of them is low in tumor samples, the expression of the other one would be very high in normal cells (Figure 5). Thus, for correlated genes, expression of one of them will allow prediction of the behavior of the correlated genes. For instance, MCU and MICU1 are very positively co-regulated with Orai2 and RYR1 and are thus extremely and negatively co-regulated with RYR2. Thus, when the expression of MCU and MICU1 increases, the expression of Orai2 and RYR1 will be also high, whereas the expression of RYR2 would be much lower. Nevertheless, not all genes are co-regulated, either positively or negatively, and there are many genes that are not co-regulated. For example, MCU and MICU1 do not co-regulate with ATP2A3, and they are almost uncorrelated with TRPM3 and TRPM4. Therefore, while the expression Figure 3. Heatmap Hierarchical Cluster of samples with healthy phenotype ( A – D ) and tumor phenotype ( E – H ). Pseudocolor scale shows hierarchical distance from minimum (0, blue) to maximum (6, red). Int. J. Mol. Sci. 2017, 18, 922 6 of 22 Figure 3. Heatmap Hierarchical Cluster of samples with healthy phenotype (A–D) and tumor phenotype (E–H). Pseudocolor scale shows hierarchical distance from minimum (0, blue) to maximum (6, red). Figure 4. Hierarchical clustering of genes. Genes have been divided into six different groups according to their Euclidean distances relative to their expression values (in RPKM). 2.3. Gene Correlations Next, we questioned possible correlations between gene couples. In this way, correlation values for each gene couple were plotted in the heatmap shown in Figure 5. Thus, in the case that two genes are positively and strongly correlated, if the expression of one of them is low in tumor samples, the expression of the other one will also be low in the normal sample. Conversely, if the pair of genes is negatively and strongly correlated and the expression of one of them is low in tumor samples, the expression of the other one would be very high in normal cells (Figure 5). Thus, for correlated genes, expression of one of them will allow prediction of the behavior of the correlated genes. For instance, MCU and MICU1 are very positively co-regulated with Orai2 and RYR1 and are thus extremely and negatively co-regulated with RYR2. Thus, when the expression of MCU and MICU1 increases, the expression of Orai2 and RYR1 will be also high, whereas the expression of RYR2 would be much lower. Nevertheless, not all genes are co-regulated, either positively or negatively, and there are many genes that are not co-regulated. For example, MCU and MICU1 do not co-regulate with ATP2A3, and they are almost uncorrelated with TRPM3 and TRPM4. Therefore, while the expression Figure 4. Hierarchical clustering of genes. Genes have been divided into six different groups according to their Euclidean distances relative to their expression values (in RPKM). 2.3. Gene Correlations Next, we questioned possible correlations between gene couples. In this way, correlation values for each gene couple were plotted in the heatmap shown in Figure 5. Thus, in the case that two genes are positively and strongly correlated, if the expression of one of them is low in tumor samples, the expression of the other one will also be low in the normal sample. Conversely, if the pair of genes is negatively and strongly correlated and the expression of one of them is low in tumor samples, the expression of the other one would be very high in normal cells (Figure 5). Thus, for correlated genes, expression of one of them will allow prediction of the behavior of the correlated genes. For instance, Int. J. Mol. Sci. 2017,18, 922 7 of 23 MCU and MICU1 are very positively co-regulated with Orai2 and RYR1 and are thus extremely and negatively co-regulated with RYR2. Thus, when the expression of MCU and MICU1 increases, the expression of Orai2 and RYR1 will be also high, whereas the expression of RYR2 would be much lower. Nevertheless, not all genes are co-regulated, either positively or negatively, and there are many genes that are not co-regulated. For example, MCU and MICU1 do not co-regulate with ATP2A3, and they are almost uncorrelated with TRPM3 and TRPM4. Therefore, while the expression of MICU1 and MCU is enhanced in the tumor phenotype, the expression of ATP2A3,TRPM3 and TRPM4 seems not to vary much. Therefore, this analysis enabled us to uncover which genes behave in the same way and which do not when comparing healthy and tumor phenotypes. In addition, it can predict whether or not the behaviors are similar. Int. J. Mol. Sci. 2017, 18, 922 7 of 22 of MICU1 and MCU is enhanced in the tumor phenotype, the expression of ATP2A3, TRPM3 and TRPM4 seems not to vary much. Therefore, this analysis enabled us to uncover which genes behave in the same way and which do not when comparing healthy and tumor phenotypes. In addition, it can predict whether or not the behaviors are similar. Figure 5. Correlation or co-regulation between couples of genes. The larger the circle and the darker the color, the higher the correlation (either positive in blue or negative in red) between each pair of genes. 2.4. Principal Component Analysis The Principal Component Analysis describes the variation produced by a multivariate observation, such that new variables are made from linear combinations of the original variables. These new variables are known as Principal Components (PCs). Thus, if the observation has p original variables, up to p PCs can be made, which are sorted by the amount of explained variance by each of them, where PC1 explains the largest amount of variance, followed by PC2, and so on. Therefore, this analysis is intended to reduce the dimension of the observations. With the new dimensions selected, PCs explain the largest possible amount of variance. Indeed, a criterion for deciding how many PCs to keep is that the proportion of variance explained for all PCs selected is larger than 70%. Other criteria are taken from the Decay graph, which represents the explained variance by each PC against the corresponding PC. Thus, the number of PCs located in the Decay graph before the slope of the graph changes drastically (Figure 6) shows that the variance explained does not increase much despite considering more PCs. As it is really difficult for multivariate data to verify the assumption that they fit a normal distribution, the Principal Component Analysis is considered as a kind of exploratory data analysis. Figure 5. Correlation or co-regulation between couples of genes. The larger the circle and the darker the color, the higher the correlation (either positive in blue or negative in red) between each pair of genes. 2.4. Principal Component Analysis The Principal Component Analysis describes the variation produced by a multivariate observation, such that new variables are made from linear combinations of the original variables. These new variables are known as Principal Components (PCs). Thus, if the observation has p original variables, up to p PCs can be made, which are sorted by the amount of explained variance by each of them, where PC1 explains the largest amount of variance, followed by PC2, and so on. Therefore, this analysis is intended to reduce the dimension of the observations. With the new dimensions selected, PCs explain the largest possible amount of variance. Indeed, a criterion for deciding how many PCs to keep is that the proportion of variance explained for all PCs selected is larger than 70%. Other criteria are taken from the Decay graph, which represents the explained variance by each PC against the corresponding PC. Thus, the number of PCs located in the Decay graph before the slope of the graph changes drastically (Figure 6) shows that the variance explained does not increase much despite considering more PCs. Int. J. Mol. Sci. 2017,18, 922 8 of 23 Int. J. Mol. Sci. 2017, 18, 922 8 of 22 This is why data have been filtered and transformed previously—to fit normal distribution as much as possible. Furthermore, data have been centered with their mean, and standardized with their variance. Figure 6. Principal Component Analysis (PCA), where genes are variables and samples are observations. The proportion of variance explained by PC1 is equal to 59.82%, and 16.19% for PC2. (A) PC2 vs. PC1; (B) Decay Variance Explained graph; (C) Correlation coefficients. In the present study, a Principal Component Analysis between samples as a function of the expression profile of p genes has been carried out and the results are shown in Figure 6. We found that PC1 clearly explains the difference between phenotypes, since the projections of the values for each sample over PC1 show how the samples belong to a healthy phenotype and are well separated from the samples belonging to the tumor phenotype. Since the variance explained by PC1 is 59.82%, and the one explained by the two first PCs is 76.79%, together with the fact that the slope of the Decay graph changes drastically from PC2, it is a good decision to keep only the first two PCs. Given that PC1 clearly reproduces the two different groups related with the phenotype, it is interesting to evaluate the influence of each original variable gene, which is proportional to a coefficient associated with the linear combination for each gene. The way in which this is evaluated in the present study is by estimating the correlation between each gene and PC1. Genes with correlations larger than |0.7| are shown in Table 2. Accordingly, it is possible to obtain the following expression (2): PC1= 󰇛.󰇜· .   (2) where PC1 is the value of the Principal Component 1, rlog(expression of gene.i) is the value of expression for each gene after transformation (see supplementary data for raw expression data of individual genes) and βgene.i is the coefficient value for the same gene (Table 2). This enables the sorting of a given sample to the normal or tumor phenotype. It is important to take into account that the values of the original variables have been transformed, centered and standardized as reported. Figure 6. Principal Component Analysis (PCA), where genes are variables and samples are observations. The proportion of variance explained by PC1 is equal to 59.82%, and 16.19% for PC2. ( A ) PC2 vs. PC1; (B) Decay Variance Explained graph; (C) Correlation coefficients. As it is really difficult for multivariate data to verify the assumption that they fit a normal distribution, the Principal Component Analysis is considered as a kind of exploratory data analysis. This is why data have been filtered and transformed previously—to fit normal distribution as much as possible. Furthermore, data have been centered with their mean, and standardized with their variance. In the present study, a Principal Component Analysis between samples as a function of the expression profile of p genes has been carried out and the results are shown in Figure 6. We found that PC1 clearly explains the difference between phenotypes, since the projections of the values for each sample over PC1 show how the samples belong to a healthy phenotype and are well separated from the samples belonging to the tumor phenotype. Since the variance explained by PC1 is 59.82%, and the one explained by the two first PCs is 76.79%, together with the fact that the slope of the Decay graph changes drastically from PC2, it is a good decision to keep only the first two PCs. Given that PC1 clearly reproduces the two different groups related with the phenotype, it is interesting to evaluate the influence of each original variable gene, which is proportional to a coefficient associated with the linear combination for each gene. The way in which this is evaluated in the present study is by estimating the correlation between each gene and PC1. Genes with correlations larger than |0.7| are shown in Table 2. Accordingly, it is possible to obtain the following expression (2): PC1 = p ∑ i=1 rlog(expression o f gene.i)·βgene.i(2) where PC1 is the value of the Principal Component 1, rlog(expression of gene.i) is the value of expression for each gene after transformation (see supplementary data for raw expression data of individual genes) and β gene.i is the coefficient value for the same gene (Table 2). This enables the sorting of a given sample to the normal or tumor phenotype. It is important to take into account that the values of the original variables have been transformed, centered and standardized as reported. Int. J. Mol. Sci. 2017,18, 922 9 of 23 Table 2. Coefficient values from Principal Component (PC1) corresponding to each gene. Raw data were transformed, filtered, centered and standardized. To obtain PC1 value for a given sample, use coefficients in this table in expression 2. Gene Expression Gene Expression Gene Expression Gene Expression Gene Expression ATP2A1 −0.166 CACNA1C −0.174 ITPR1 0.143 ORAI3 0.090 TRPA1 −0.147 ATP2A2 0.148 CACNA1D 0.162 ITPR2 −0.173 PICALM −0.089 TRPM3 −0.070 ATP2A3 −0.088 CACNA1G 0.072 ITPR3 0.152 PKD2 0.020 TRPM4 −0.113 ATP2B1 0.172 CACNA1H −0.173 MCOLN1 −0.175 PKD2L1 −0.169 TRPM5 −0.168 ATP2B4 −0.173 CACNA1I −0.084 MCOLN2 −0.174 RYR1 0.103 TRPM7 0.105 ATP2C1 −0.121 CASQ1 0.093 MCOLN3 −0.148 RYR2 −0.174 TRPV1 −0.141 ATP2C2 0.142 CCDC109B −0.167 MCU 0.149 SLC8A1 −0.010 TRPV2 0.096 BCL2 −0.141 CCDC90A 0.140 MICU1 0.131 SLC8A2 0.165 TRPV3 −0.060 C22orf32 0.131 EFCAB4B −0.169 ORAI1 0.061 STIM1 0.174 TRPV4 0.026 CACNA1B −0.141 EFHA1 −0.149 ORAI2 0.168 STIM2 0.057 TRPV6 0.133 VDAC1 −0.175 VDAC2 0.147 VDAC3 −0.170 TRPC1 0.041 Int. J. Mol. Sci. 2017,18, 922 16 of 23 in SERCA activity could lead to possible differences in Ca 2+ store content between normal and CRC cells. In fact, we have reported recently that CRC cells display the partial depletion of Ca 2+ stores [ 4 ]. This effect has been partially attributed to loss of ER Ca 2+ sensor STIM2 [ 4 ]. Whether changes in SERCA isoform contribute to this remodeling as well remains to be established. Expression of Na + /Ca 2+ exchanger isoforms is shown in Figure 12. These co-transporters normally extrude Ca 2+ back to the external medium in exchange for Na + using the electrochemically favorable gradient for Na + . There are three isoforms named NCX1, NCX2 and NCX3, but only NCX1 and NCX2 are expressed in normal and colon cancer cells. Interestingly, while expression of NCX1 is similar in normal and CRC cells, expression of NCX2 is dramatically enhanced in CRC cells, actually showing the largest fold increase in expression in tumor cells of all upregulated genes tested (Figure 12). Int. J. Mol. Sci. 2017, 18, 922 17 of 22 stores [4]. This effect has been partially attributed to loss of ER Ca 2+ sensor STIM2 [4]. Whether changes in SERCA isoform contribute to this remodeling as well remains to be established. Expression of Na+/Ca2+ exchanger isoforms is shown in Figure 12. These co-transporters normally extrude Ca2+ back to the external medium in exchange for Na+ using the electrochemically favorable gradient for Na+. There are three isoforms named NCX1, NCX2 and NCX3, but only NCX1 and NCX2 are expressed in normal and colon cancer cells. Interestingly, while expression of NCX1 is similar in normal and CRC cells, expression of NCX2 is dramatically enhanced in CRC cells, actually showing the largest fold increase in expression in tumor cells of all upregulated genes tested (Figure 12). Figure 12. Expression of Na+/Ca2+ exchanger isoforms in normal and colon cancer cells. Expression levels of NCXs in normal colonic cells (grey bars) and colon cancer cells (red bars). * p statistically significant by three independent methods. Gene expression levels for molecular players involved in mitochondrial Ca2+ transport are shown in Figure 13. Ca2+ enters mitochondria down its electrochemical gradient through the mitochondrial Ca2+ uniporter (MCU) that is modulated by several regulatory proteins including mitochondrial Ca2+ uptake (MICU) isoforms 1 to 3, mitochondrial calcium uniporter regulator 1 (MICUR1), mitochondrial calcium uniporter dominant negative beta (MCUb) and MCU regulator (EMRE). Voltage-dependent anion channels (VDAC) 1 to 3 in the outer mitochondrial membrane have been involved in apoptosis. Interestingly, messengers involved in mitochondrial Ca2+ transport are very highly expressed in both normal and CRC cells (Figure 13). In addition, all genes but MICU3 are expressed in both normal and tumor cells (Figure 13). Expression of the channel MCU and its positive modulator MICU1 is enhanced in tumor cells (Figure 13A,B), while expression of negative modulator MICU2 is decreased. Expression of other MCU modulators including MICUR1 and EMRE is similar (Figure 13B,C,E), while MCUb is significantly downregulated in colon cancer cells (Figure 13D). Finally, all three VDACs are expressed at very high levels in both normal and tumor cells (Figure 13F). Interestingly, expression of all of them is significantly different in colon cancer cells. VDAC1 and 3 are downregulated in colon cancer cells while VDAC2 is significantly overexpressed in cancer cells. Finally, expression of other related proteins is shown in Figure 14. Bcl-2 modulates Ca2+ release and apoptosis. This gene is significantly downregulated in cancer cells (Figure 14) while expression of two other proteins including calsequestrin 1 and phosphatidylinositol binding clathrin assembly protein (PICALM), is similar (Figure 14B,C). No expression of calsequestrin2 and phospholamban is observed in normal and colon cancer cells. Figure 12. Expression of Na + /Ca 2+ exchanger isoforms in normal and colon cancer cells. Expression levels of NCXs in normal colonic cells (grey bars) and colon cancer cells (red bars). * pstatistically significant by three independent methods. Gene expression levels for molecular players involved in mitochondrial Ca 2+ transport are shown in Figure 13. Ca 2+ enters mitochondria down its electrochemical gradient through the mitochondrial Ca 2+ uniporter (MCU) that is modulated by several regulatory proteins including mitochondrial Ca 2+ uptake (MICU) isoforms 1 to 3, mitochondrial calcium uniporter regulator 1 (MICUR1), mitochondrial calcium uniporter dominant negative beta (MCUb) and MCU regulator (EMRE). Voltage-dependent anion channels (VDAC) 1 to 3 in the outer mitochondrial membrane have been involved in apoptosis. Interestingly, messengers involved in mitochondrial Ca 2+ transport are very highly expressed in both normal and CRC cells (Figure 13). In addition, all genes but MICU3 are expressed in both normal and tumor cells (Figure 13). Expression of the channel MCU and its positive modulator MICU1 is enhanced in tumor cells (Figure 13A,B), while expression of negative modulator MICU2 is decreased. Expression of other MCU modulators including MICUR1 and EMRE is similar (Figure 13B,C,E), while MCUb is significantly downregulated in colon cancer cells (Figure 13D). Finally, all three VDACs are expressed at very high levels in both normal and tumor cells (Figure 13F). Interestingly, expression of all of them is significantly different in colon cancer cells. VDAC1 and 3 are downregulated in colon cancer cells while VDAC2 is significantly overexpressed in cancer cells. Finally, expression of other related proteins is shown in Figure 14. Bcl-2 modulates Ca 2+ release and apoptosis. This gene is significantly downregulated in cancer cells (Figure 14) while expression of two other proteins including calsequestrin 1 and phosphatidylinositol binding clathrin assembly protein (PICALM), is similar (Figure 14B,C). No expression of calsequestrin2 and phospholamban is observed in normal and colon cancer cells. Int. J. Mol. Sci. 2017,18, 922 17 of 23 Int. J. Mol. Sci. 2017, 18, 922 18 of 22 Figure 13. Expression of genes involved in mitochondrial Ca2+ transport in normal and colon cancer cells. Expression levels of Mitochondrial Calcium Uniporter (MCU) (A), Mitochondrial Calcium Uptake (MICU) isoforms (B), Mitochondrial Calcium Uptake Regulator 1 (MICUR1) (C), Mitochondrial Calcium Uniporter Dominant Negative Beta Subunit (MCUb) (D), Essential MCU regulator (EMRE) (E) and voltage-dependent anion channels (VDACs) (F) in normal colonic cells (grey bars) and colon cancer cells (red bars). * p statistically significant by three independent methods. Figure 14. Expression of selected genes coding for other selected proteins. Expression levels of BCL2, (A) calsequestrins CASQ1 and CASQ2 (B), PICALM (C) and PLN (Phospholamban) (D) in normal colonic cells (grey bars) and colon cancer cells (red bars). * p statistically significant by three independent methods. Figure 13. Expression of genes involved in mitochondrial Ca 2+ transport in normal and colon cancer cells. Expression levels of Mitochondrial Calcium Uniporter (MCU) ( A ), Mitochondrial Calcium Uptake (MICU) isoforms ( B ), Mitochondrial Calcium Uptake Regulator 1 (MICUR1) ( C ), Mitochondrial Calcium Uniporter Dominant Negative Beta Subunit (MCUb) ( D ), Essential MCU regulator (EMRE) ( E ) and voltage-dependent anion channels (VDACs) ( F ) in normal colonic cells (grey bars) and colon cancer cells (red bars). * pstatistically significant by three independent methods. Int. J. Mol. Sci. 2017, 18, 922 18 of 22 Figure 13. Expression of genes involved in mitochondrial Ca2+ transport in normal and colon cancer cells. Expression levels of Mitochondrial Calcium Uniporter (MCU) (A), Mitochondrial Calcium Uptake (MICU) isoforms (B), Mitochondrial Calcium Uptake Regulator 1 (MICUR1) (C), Mitochondrial Calcium Uniporter Dominant Negative Beta Subunit (MCUb) (D), Essential MCU regulator (EMRE) (E) and voltage-dependent anion channels (VDACs) (F) in normal colonic cells (grey bars) and colon cancer cells (red bars). * p statistically significant by three independent methods. Figure 14. Expression of selected genes coding for other selected proteins. Expression levels of BCL2, (A) calsequestrins CASQ1 and CASQ2 (B), PICALM (C) and PLN (Phospholamban) (D) in normal colonic cells (grey bars) and colon cancer cells (red bars). * p statistically significant by three independent methods. Figure 14. Expression of selected genes coding for other selected proteins. Expression levels of BCL2, ( A ) calsequestrins CASQ1 and CASQ2 ( B ), PICALM ( C ) and PLN (Phospholamban) ( D ) in normal colonic cells (grey bars) and colon cancer cells (red bars). * pstatistically significant by three independent methods. Int. J. Mol. Sci. 2017,18, 922 18 of 23 Changes in expression of all significantly expressed transcripts are shown in Figure 15A. To best appreciate the differences, all genes that are significantly downregulated in tumor cells are removed for the genes shown in normal cells (Figure 15B) but show up in the tumor cell (Figure 15C). In summary, we report here a comprehensive transcriptomic analysis of Ca 2+ remodeling in colorectal cancer cells. A similar analysis has been recently carried out in glioblastoma using data gathered from repository data sets of a large series of glioblastomas [ 23 ]. Data shown here were generated using Ion Torrent RNA-sequencing from mRNA isolated from normal human colonic NCM460 cells and colon cancer HT29 cells. Therefore, while the study has the limitation that both cell lines may not reflect entirely normal human colonic and colon cancer cells, respectively, it does profit from the fact that changes in intracellular Ca 2+ handling between the two cell lines have been analyzed in detail [4]. The results provide a strong molecular basis for Ca2+ remodeling in CRC. Int. J. Mol. Sci. 2017, 18, 922 19 of 22 Changes in expression of all significantly expressed transcripts are shown in Figure 15A. To best appreciate the differences, all genes that are significantly downregulated in tumor cells are removed for the genes shown in normal cells (Figure 15B) but show up in the tumor cell (Figure 15C). In summary, we report here a comprehensive transcriptomic analysis of Ca2+ remodeling in colorectal cancer cells. A similar analysis has been recently carried out in glioblastoma using data gathered from repository data sets of a large series of glioblastomas [23]. Data shown here were generated using Ion Torrent RNA-sequencing from mRNA isolated from normal human colonic NCM460 cells and colon cancer HT29 cells. Therefore, while the study has the limitation that both cell lines may not reflect entirely normal human colonic and colon cancer cells, respectively, it does profit from the fact that changes in intracellular Ca2+ handling between the two cell lines have been analyzed in detail [4]. The results provide a strong molecular basis for Ca2+ remodeling in CRC. Figure 15. Molecular players differentially expressed in NCM460 normal colonic and HT29 colorectal cancer cells. Genes significantly upregulated (↑) and downregulated (↓) in colorectal cancer cells vs. normal colonic cells are shown in specific locations in the plasma membrane, the endoplasmic reticulum (ER) and mitochondria (A). To best appreciate the differences, molecular players downregulated in cancer cells have been removed in the prototypical, normal colonic cell (B), whereas only genes upregulated in CRC are shown in the prototypic CRC cell (C). Molecular players in red are those with the largest fold changes in expression from the normal to the tumor phenotype. PMCAs, plasma membrane calcium ATPases; VOCCS, voltage-operated calcium channels; SOCE, store-operated calcium entry; TRPs, Transient Receptor Potential channels. Figure 15. Molecular players differentially expressed in NCM460 normal colonic and HT29 colorectal cancer cells. Genes significantly upregulated ( ↑ ) and downregulated ( ↓ ) in colorectal cancer cells vs. normal colonic cells are shown in specific locations in the plasma membrane, the endoplasmic reticulum (ER) and mitochondria ( A ). To best appreciate the differences, molecular players downregulated in cancer cells have been removed in the prototypical, normal colonic cell ( B ), whereas only genes upregulated in CRC are shown in the prototypic CRC cell ( C ). Molecular players in red are those with the largest fold changes in expression from the normal to the tumor phenotype. PMCAs, plasma membrane calcium ATPases; VOCCS, voltage-operated calcium channels; SOCE, store-operated calcium entry; TRPs, Transient Receptor Potential channels. Int. J. Mol. Sci. 2017,18, 922 19 of 23 3. Materials and Methods 3.1. Materials HT29 cells were donated by JoséCarlos Fernández-Checa (CSIC, Barcelona, Spain). NCM460 cells were obtained after a material transfer agreement with INCELL Corporation (San Antonio, TX, USA). Dulbecco’s Modified Eagle’s Medium (DMEM), penicillin, streptomycin, L-glutamine and fetal bovine serum came from Lonza (Basel, Switzerland). M3:10TM medium is from INCELL Corporation, San Antonio, TX, USA). 3.2. Cell Culture Cells are cultured in DMEM 1 g/L glucose or in M3:10TM medium as reported previously [ 4 ] and supplemented with 1% penicillin-streptomycin, 1% L-glutamine and 10% fetal bovine serum. Cells are maintained in standard conditions (37 ◦ C, 10% CO 2 ) and cultured once a week. All cells were used at passages 3 to 10. 3.3. mRNA Isolation and Ion Torrent Reading Total cellular RNA was isolated from NCM460 and HT29 cells using Trizol reagent (Invitrogen, Carlsbad, CA, USA). Four totally independent samples from each cell line were used for molecular analysis. Extracted RNA integrity was tested by electrophoresis on agarose gels and the purity and concentration were determined by spectrophotometry. RNA was reverse transcribed using a High Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific, Waltham, MA, USA) and the cDNA diluted prior to PCR amplification. A targeted, multiplex PCR primer panel was designed using the custom Ion Ampliseq Designer v1.2 (Thermo Fisher Scientific), where target regions were selected for 77 genes directly or indirectly involved in Ca 2+ transport and listed in Table 1. The panel was designed to amplify PCR products appropriate for use with RNA from eight different cellular cultures, where four of them belonged to one cellular line and the other four to another. Sequencing libraries were prepared with Ion AmpliSeq RNA Library kit (Thermo Fisher Scientific) using the custom primers. After primer digestion, adapters and Ion Xpress Barcode Adapters (Thermo Fisher Scientific) were ligated to the amplicons followed by magnetic bead purification. Amplicon size and DNA concentration were measured using an Agilent High Sensitivity DNA Kit (Agilent Technologies, Santa Clara, CA, USA) according to the manufacturer’s recommendation. A total of eight picomols of each of the eight samples were pooled for emulsion PCR (ePCR) on Ion Sphere Particles (ISPs) using the Ion PGM Hi-Q OT2 Kit (Thermo Fisher Scientific) using the Ion OneTouch 2 Instrument. Following automated Ion OneTouch, enrichment of template-positive ISP samples were loaded on an Ion 316 V2 Chip and sequenced on a Personal Genome Machine (PGM-Thermo Fisher Scientific) System with Ion PGM Hi-Q Sequencing Kit (Thermo Fisher Scientific). Raw sequencing reads were initially filtered for high quality reads, and the adaptors were removed using the Ion Torrent Suite 4.4.3; reads were then aligned with the hg19_rna reference sequence by TMAP [ 42 ] using default parameters. Resulting binary format files (BAM) storing sequence data were processed through an in-house quality control (QC) filter in order to keep those reads whose quality is larger or equal than Q20 and remove polyclonal and primer dimers. After that, outcomes show that 253 Mbp were amplified: the length read median is 108 bp. Reads were aligned against hg19_rna reference genome using GATK LeftAlignIndels module, for 77 genes of interest, where average coverage was 93%, and mean raw accuracy was 99.3%. BAM files were aligned. Amplicon primers were trimmed from aligned reads by Torrent Suite 4.4.3. 3.4. Transformation and Filtration of Raw Data Set Raw data were transformed with an R/Bioconductor package called DESeq2 [ 33 ], applying rlog function on data, and normalized against their gene length, obtained from R package EDASeq [ 43 ]. The R code used and the raw data employed in this study is provided as a supplementary file data (esetCalcium file) so that anyone can reproduce the results. Subsequently, using R software, Int. J. Mol. Sci. 2017,18, 922 20 of 23 the transformed data set was filtered in order to remove those genes that are not expressed by samples; thus genes whose median of expression profile was lower than 0.125 RPKMs were removed [29]. 3.5. Exploratory Data Analysis A transformed and filtered data set was analyzed through Exploratory Data Analysis. At first, for both samples and genes, a Hierarchical Clustering Analysis using the Euclidean distance was carried out as well as a correlation between different genes, using R software and stats, heatmap [ 44 ] and Corrplot [ 45 ] R packages. A second method, the Principal Components Analysis, was applied in order to reveal patterns and reduce the dimensionality of the data, using R software and stats, EDASeq [ 43 ] and IRanges R packages. 3.6. Differential Expression Analysis In order to carry out a differential expression analysis, three different methods were used to assess significance. These methods are EdgeR [ 32 ], DESeq2 [ 33 ] and Limma [ 34 ], and they were implemented using R software and EnrichmentBrowser R package [ 35 ] because this package allows the performing of differential expression analyses either based on Limma, EdgeR or DESeq2, and adjusted pvalue using the Benjamini-Hochberg method [ 46 ]. A critical issue is normalization of data. We have used the default method of normalization in each package as suggested by Seyednasrollah [ 27 ] which has a low influence on the analysis outcomes. In addition, we have controlled False Discovery Rate (FDR), proposed by Benjamini-Hochberg [ 46 ] where the percent of false positives among whole test where null hypothesis has been rejected is controlled. Supplementary Materials: Supplementary materials can be found at www.mdpi.com/1422-0067/18/5/922/s1. Acknowledgments: We thank David del Bosque for technical assistance. This work was supported by grants BFU2015-70131-R from Ministry of Economy and Competitivity of Spain to Carlos Villalobos and VA145U13. from Regional Government of Castilla y León, Spain to Lucía Núñez. Enrique Pérez-Riesgo is supported by a fellowship from the Spanish Association Against Cancer (AECC). Daniel Ubierna is supported by predoctoral fellowships from Junta de Castilla y León, Spain and the EU Social Fund. Lucía G. Gutiérrez is supported by a predoctoral fellowship from the University of Valladolid, Spain. Author Contributions: Enrique Pérez-Riesgo carried out all the statistical analyses. Lucía G. Gutiérrez and Daniel Ubierna cultured the cells, isolated quality mRNA for analysis and carried out preliminary analysis. Alberto Acedo carried out the Ion Torrent transcriptomic sequencing. Mary P Moyer developed and provided NCM460 cells. Lucía Núñez and Carlos Villalobos the idea, directed the study and developed the cartoon model. Carlos Villalobos and Enrique Pérez-Riesgo wrote the manuscript. Conflicts of Interest: The authors declare no conflict of interest. The funding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results. Abbreviations CRC Colorectal cancer VOCC Voltage-operated Ca2+ channels SOCE Store-operated Ca2+ entry TRP Transient receptor potential channels IP3R Inositol trisphosphate receptor RyR Ryanodine receptor PMCA Plasma membrane Ca2+/ATPase SERCA Sarcoplasmic and endoplasmic reticulum Ca2+/ATPase SPCA Secretory pathway Ca2+/ATPase MCU Mitochondrial calcium uniporter Int. J. Mol. Sci. 2017,18, 922 21 of 23 References 1. Berridge, M.J. The Inositol trisphosphate/calcium signaling pathway in health and disease. Physiol. Rev. 2016,96, 1261–1296. [CrossRef] [PubMed] 2. Prevarskaya, N.; Ouadid-Ahidouch, H.; Skryma, R.; Shuba, Y. Remodelling of Ca 2+ transport in cancer: How it contributes to cancer hallmarks? Philos. Trans. R. Soc. Lond. B Biol. Sci. 2014 ,369, 20130097. [CrossRef] [PubMed] 3. Stewart, T.A.; Yapa, K.T.; Monteith, G.R. Altered calcium signaling in cancer cells. Biochim. Biophys. Acta 2015,1848, 2502–2511. [CrossRef] [PubMed] 4. Sobradillo, D.; Hernández-Morales, M.; Ubierna, D.; Moyer, M.P.; Núñez, L.; Villalobos, C. A reciprocal shift in transient receptor potential channel 1 (TRPC1) and stromal interaction molecule 2 (STIM2) contributes to Ca 2+ remodeling and cancer hallmarks in colorectal carcinoma cells. J. Biol. Chem. 2014 ,289, 28765–28782. [CrossRef] [PubMed] 5. Villalobos, C.; Sobradillo, D.; Hernández-Morales, M.; Núñez, L. Remodeling of calcium entry pathways in cancer. Adv. Exp. Med. Biol. 2016,898, 449–466. [PubMed] 6. Villalobos, C.; Sobradillo, D.; Hernández-Morales, M.; Núñez, L. Calcium remodeling in colorectal cancer. Biochim. Biophys. Acta 2017. [CrossRef] [PubMed] 7. Putney, J.W., Jr. A model for receptor-regulated calcium entry. Cell Calcium 1986,7, 1–12. [CrossRef] 8. Parekh, A.B.; Putney, J.W., Jr. Store-operated calcium channels. Physiol. Rev. 2005 ,85, 757–810. [CrossRef] [PubMed] 9. Liou, J.; Kim, M.L.; Heo, W.D.; Jones, J.T.; Myers, J.W.; Ferrell, J.E., Jr.; Meyer, T. STIM is a Ca 2+ sensor essential for Ca 2+ -store-depletion-triggered Ca 2+ influx. Curr. Biol. 2005 ,15, 1235–1241. [CrossRef] [PubMed] 10. Feske, S.; Gwack, Y.; Prakriya, M.; Srikanth, S.; Puppel, S.H.; Tanasa, B.; Hogan, P.G.; Lewis, R.S.; Daly, M.; Rao, A. A mutation in Orai1 causes immune deficiency by abrogating CRAC channel function. Nature 2006 , 441, 179–185. [CrossRef] [PubMed] 11. Cheng, K.T.; Ong, H.L.; Liu, X.; Ambudkar, I.S. Contribution and regulation of TRPC channels in store-operated Ca2+ entry. Curr. Top. Membr. 2013,71, 149–179. [PubMed] 12. Strehler, E.E. Plasma membrane calcium ATPases: From generic Ca 2+ sump pumps to versatile systems for fine-tuning cellular Ca2+.Biochem. Biophys. Res. Commun. 2015,460, 26–33. [CrossRef] [PubMed] 13. Stammers, A.N.; Susser, S.E.; Hamm, N.C.; Hlynsky, M.W.; Kimber, D.E.; Kehler, D.S.; Duhamel, T.A. The regulation of sarco(endo)plasmic reticulum calcium-ATPases (SERCA). Can. J. Physiol. Pharmacol. 2015 , 93, 843–854. [CrossRef] [PubMed] 14. Brini, M.; Calì, T.; Ottolini, D.; Carafoli, E. Calcium pumps: Why so many? Compr. Physiol. 2012 ,2, 1045–1060. [PubMed] 15. Giladi, M.; Shor, R.; Lisnyansky, M.; Khananshvili, D. Structure-functional basis of ion transport in sodium-calcium exchanger (NCX) proteins. Int. J. Mol. Sci. 2016,17, 1949. [CrossRef] [PubMed] 16. De Stefani, D.; Patron, M.; Rizzuto, R. Structure and function of the mitochondrial calcium uniporter complex. Biochim. Biophys. Acta 2015,1853, 2006–2011. [CrossRef] [PubMed] 17. Villalobos, C.; Núñez, L.; Montero, M.; García, A.G.; Alonso, M.T.; Chamero, P.; Alvarez, J.; García-Sancho, J. Redistribution of Ca 2+ among cytosol and organella during stimulation of bovine chromaffin cells. FASEB J. 2002,16, 343–353. [CrossRef] [PubMed] 18. Valero, R.A.; Senovilla, L.; Núñez, L.; Villalobos, C. The role of mitochondrial potential in control of calcium signals involved in cell proliferation. Cell Calcium 2008,44, 259–269. [CrossRef] [PubMed] 19. Bonnet, S.; Archer, S.L.; Allalunis-Turner, J.; Haromy, A.; Beaulieu, C.; Thompson, R.; Lee, C.T.; Lopaschuk, G.D.; Puttagunta, L.; Bonnet, S.; et al. A mitochondria-K + channel axis is suppressed in cancer and its normalization promotes apoptosis and inhibits cancer growth. Cancer Cell 2007 ,11, 37–51. [CrossRef] [PubMed] 20. Hoth, M.; Fanger, C.M.; Lewis, R.S. Mitochondrial regulation of store-operated calcium signaling in T lymphocytes. J. Cell Biol. 1997,137, 633–648. [CrossRef] [PubMed] 21. Gilabert, J.A.; Parekh, A.B. Respiring mitochondria determine the pattern of activation and inactivation of the store-operated Ca2+ current ICRAC. EMBO J. 2000,19, 6401–6407. [CrossRef] [PubMed] Int. J. Mol. Sci. 2017,18, 922 22 of 23 22. Núñez, L.; Valero, R.A.; Senovilla, L.; Sanz-Blasco, S.; García-Sancho, J.; Villalobos, C. Cell proliferation depends on mitochondrial Ca 2+ uptake: Inhibition by salicylate. J. Physiol. 2006 ,571, 57–73. [CrossRef] [PubMed] 23. Robil, N.; Petel, F.; Kilhoffer, M.C.; Haiech, J. Glioblastoma and calcium signaling: Analysis of calcium toolbox expression. Int. J. Dev. Biol. 2015,59, 407–415. [CrossRef] [PubMed] 24. Moyer, M.P.; Manzano, L.A.; Merriman, R.L.; Stauffer, J.S.; Tanzer, L.R. NCM460, a normal human colon mucosal epithelial cell line. In Vitro Cell. Dev. Biol. Anim. 1996,32, 315–317. [CrossRef] [PubMed] 25. Alcarraz-Vizan, G.; Sánchez-Tena, S.; Moyer, M.P.; Cascante, M. Validation of NCM460 cell model as control in antitumor strategies targeting colon adenocarcinoma metabolic reprogramming: Trichostatin A as a case study. Biochim. Biophys. Acta 2013,1840, 1634–1639. [CrossRef] [PubMed] 26. Marshall, C.J.; Franks, L.M.; Carbonell, A.W. Markers of neoplastic transformation in epithelial cell lines derived from human carcinomas. J. Natl. Cancer Inst. 1977,58, 1743–1751. [CrossRef] [PubMed] 27. Hackett, N.R. RNA-Seq quantification of the human small airway epithelium transcriptome. BMC Genom. 2012,13, 82. [CrossRef] [PubMed] 28. Seyednasrollah, F. Comparison of software packages for detecting differential expression in RNA-seq studies. Brief. Bioinform. 2013,16, 59–70. [CrossRef] [PubMed] 29. Mortazavi, A. Mapping and quantifying mammalian transcriptomes by RNA-seq. Nat. Methods 2008 ,5, 621–628. [CrossRef] [PubMed] 30. Pérez, A.G. Métodos Avanzados de Estadística Aplicada; Técnicas Avanzadas; Universidad Nacional de Educación a Distancia: Madrid, Spain, 2014. (In Spanish) 31. Zhang, Z.H. A Comparative Study of Techniques for Differential Expression Analysis on RNA-Seq Data. PLoS ONE 2014,9, e103207. [CrossRef] [PubMed] 32. Robinson, M.M. EdgeR: A Bioconductor package for. Bioinformatics 2010 ,26, 139–140. [CrossRef] [PubMed] 33. Anders, S. Differencial expression analysis for sequence count data. Genome Biol. 2010 ,11, R106. [CrossRef] [PubMed] 34. Smyth, G.K. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 2004,3. [CrossRef] [PubMed] 35. Geistlinger, L.; Csaba, G.; Zimmer, R. Bioconductor’s Enrichment Browser: Seamless navigation through combined results of set- & network-based enrichment analysis. BMC Bioinform. 2016,17, 45. 36. Dziegielewska, B.; Brautigan, D.L.; Larner, J.M.; Dziegielewski, J. T-type Ca 2+ channel inhibition induces p53-dependent cell growth arrest and apoptosis through activation of p38-MAPK in colon cancer cells. Mol. Cancer Res. 2014,12, 348–358. [CrossRef] [PubMed] 37. Koslowski, M.; Sahin, U.; Dhaene, K.; Huber, C.; Türeci, O. MS4A12 is a colon-selective store-operated calcium channel promoting malignant cell processes. Cancer Res. 2008 ,68, 3458–3466. [CrossRef] [PubMed] 38. Lopez, J.J.; Albarran, L.; Gómez, L.J.; Smani, T.; Salido, G.M.; Rosado, J.A. Molecular modulators of store-operated calcium entry. Biochim. Biophys. Acta 2016,1863, 2037–2043. [CrossRef] [PubMed] 39. Shibao, K.; Fiedler, M.J.; Nagata, J.; Minagawa, N.; Hirata, K.; Nakayama, Y.; Iwakiri, Y.; Nathanson, M.H.; Yamaguchi, K. The type III inositol 1,4,5-trisphosphate receptor is associated with aggressiveness of colorectal carcinoma. Cell Calcium 2010,48, 315–323. [CrossRef] [PubMed] 40. Pierro, C.; Cook, S.J.; Foets, T.C.; Bootman, M.D.; Roderick, H.L. Oncogenic K-Ras suppresses IP3-dependent Ca 2+ release through remodeling of the isoform composition of IP3Rs and ER luminal Ca 2+ levels in colorectal cancer cell lines. J Cell Sci. 2014,127, 1607–1619. [CrossRef] [PubMed] 41. Ribiczey, P.; Tordai, A.; Andrikovics, H.; Filoteo, A.G.; Penniston, J.T.; Enouf, J.; Enyedi, A.; Papp, B.; Kovács, T. Isoform-specific up-regulation of plasma membrane Ca 2+ ATPase expression during colon and gastric cancer cell differentiation. Cell Calcium 2007,42, 590–605. [CrossRef] [PubMed] 42. Ion Torrent Analysis. Available online: https://github.com/iontorrent/TS/tree/master/Analysis/TMAP (accessed on 3 June 2015). 43. Dudoit, S.; van der Laan, M.J. Multiple Testing Procedures with Applications to Genomics; Springer: New York, NY, USA, 2008. 44. Kolde, R. Pheatmap: Pretty Heatmaps; Package Manual; The Comprehensive R Archive Network: Vienna, Austria, 2015. Int. J. Mol. Sci. 2017,18, 922 23 of 23 45. Simko, T.W. Corrplot: Visualization of a Correlation; Package Manual; The Comprehensive R Archive Network: Vienna, Austria, 2016. 46. Benjamini, Y. Controling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. B Methodol. 1995,57, 289–300. © 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).