Identification and validation of calcium extrusion-related genes prognostic signature in colon adenocarcinoma

Colon adenocarcinoma data from the TCGA and GEO databases

The transcriptomic data was obtained from the TCGA (https://gdc.cancer.gov/). A total of 524 samples were downloaded from TCGA-COAD, and the samples with incomplete clinical information, including age, gender, T stage, N stage, M stage, and TNM stage, and the overall survival time less than 30 days were excluded from this study. We included 471 samples, comprising 430 tumor samples and 41 healthy samples in the present study.The transcriptional expressions of nine calcium extrusion-related genes were examined in both tumor tissues and normal tissues (Tomczak, Czerwińska & Wiznerowicz, 2015). Table S1 displayed comprehensive details of nine genes. GSE17536, GSE29623, and GSE39582 datasets were used as validation cohorts in the present study (https://www.ncbi.nlm.nih.gov/geo/)  (Zheng et al., 2021; Liu et al., 2021). A total of 585 samples were downloaded from GSE39582, and 19 non-tumor samples were excluded in this study. Finally, 566 tumor samples were analyzed in the present study. 177 and 65 tumor samples were analyzed in GSE17536 and GSE29623, respectively. The details were shown in Table S2.

Clinical-pathological analysis

The clinical-pathological features of COAD patients were obtained from the TCGA database (https://gdc.cancer.gov/). The data is available to the public online and its usage does not violate the rights of individuals or institutions. After dividing the numeric values according to the median, a comparison was made between the low-risk and high-risk groups. A Pearson’s chi-square (χ2) test was used to assess the association between these categorical variables. Additionally, the correlation between the mRNA expression of calcium extrusion-related genes and the pathological stage of COAD patients was examined using Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/) (Cai et al., 2021).

Survival analysis in colon adenocarcinoma patients

The relationship between nine calcium extrusion-related genes expressions and the overall survival (OS) and disease-free survival (DFS) of COAD patients was investigated through the Kaplan–Meier Plotter tool (http://kmplot.com/analysis/). Survival analyses were conducted based on the median distribution, with a cutoff of 50% for high expression and 50% for low expression. The hazard ratio value and p value for each gene were uploaded to bioinformatic online software to visualize the prognosis (https://www.bioinformatics.com.cn/).

The prognosis value of the risk signature

A risk model was built using R software based on the nine calcium extrusion-related genes and the TCGA database was utilized as the training set. The GEO datasets (including GSE39582, GSE17536, and GSE29623) were employed as validation datasets. Patients were categorized into high-risk and low-risk groups based on the median value of their risk score. The disparities in patient survival between these subgroups were assessed using the R packages “survival” and “survminer”, and the receiver operating characteristic (ROC) curve was generated using the R package “survival-ROC” to calculate the area under curve (AUC) value. The transcriptional expressions of calcium extrusion-related genes in the high-risk and low-risk groups were visualized using GraphPad Prism 8.0 software. The heatmap was plotted by online tool.

Analysis of gene mutations, interaction networks, and functional enrichment

The mutation data of patients with COAD were acquired from the cBioPortal (https://www.cbioportal.org/). The tumor mutation burden (TMB) for each sample was determined using the R packages (“maftools”). TCGA-COAD samples were divided into microsatellite stability (MSS), low-frequency microsatellite instability (MSI-L), and high-frequency MSI(MSI-H) using the R packages (“TCGA biolinks”). Protein-protein interaction network analysis was performed using the STRING database (https://cn.string-db.org/), while gene-gene interaction network analysis was carried out using Gene-MANIA (http://genemania.org/). We conducted the functional enrichment of SLC8B1 and its co-expression genes using R packages (“clusterProfiler”, “ggplot2”, and “enrichplot”).

Immune cell infiltration and responsiveness to immunotherapy

The immune checkpoint genes used in this study were downloaded from a previous study (Sun et al., 2020). The relationship between the expressions of immune checkpoint genes and the mRNA expressions of seven calcium extrusion-related genes in COAD was examined. The estimation of immune cell composition was performed using CIBERSORT R package. The list of immune cell signatures was obtained from TISIDB (http://cis.hku.hk/TISIDB/download.php). Sanger box software was used to explore the correlation between the risk score and the expressions of immune cell signature genes (http://www.sangerbox.com/home.html). Additionally, the ESTIMATE algorithm, available in the R package “estimate”, was utilized to calculate the estimated score, immune score, stromal score and tumor purity. The anti-PD1 and anti-CTLA4 immunotherapy response in patients with COAD was determined using the TIDE algorithm (http://tide.dfci.harvard.edu/login/) (Jiang et al., 2018). PRJEB25780 was used as an immunotherapy validation cohort and was obtained from the European Nucleotide Archive (ENA, https://www.ebi.ac.uk/). The samples were divided into low-risk and high-risk groups based on the median risk score as previously described. For each of the 198 FDA-approved drugs, the IC50 was calculated by through R package (“pRRophetic”). The Genomics of Drug Sensitivity in Cancer database (GDSC, https://www.cancerrxgene.org/) provided information on drugs.

Cell culture

We purchased human colorectal cancer cell (HCT-116, SW620, SW480, HT-29 and RKO) from the American Type Culture Collection (ATCC, Manassas, VA, USA). SW620 and SW480 were cultured in L-15 supplemented with 10% fetal bovine serum (FBS). HCT-116 and RKO was cultured in DMEM supplemented with 10% FBS. HT-29 was cultured in RPMI-1640 supplemented with 10% FBS.

Cell transfection, virus packaging and stable cells lines

Cells were transfected with siRNAs using Lipofectamine® 2000 reagent provided by Thermo Fisher Scientific following the manufacturer’s instructions. The cells were seeded 50–60% confluent and transfected with 50 nM siRNAs. NC-siRNA was used as the control. The control was NC-siRNA. Shanghai Generay Biotech (China) synthesized all the siRNAs and following is a list of siRNA sequences. SLC8A3 siRNA: sense, 5′-GAAGAGUCCUAUGAGUUCAAGACUA-3′, SLC8A3 siRNA: antisense, 5′-UAGUCUUGAACUCAUAGGACUCUUC-3′; SLC24A2 siRNA: sense, 5′-CAACGUUGGCAUAGGCACAAUUGUA-3′, SLC24A2 siRNA: antisense, 5′-UACAAUUGUGCCUAUGCCAACGUUG-3′; SLC24A3 siRNA: sense, 5′-GGUGCCUGCUGAGGGAUUCUAUUUA-3′, SLC24A3 siRNA: antisense, 5′-UAAAUAGAAUCCCUCAGCAGGCACC-3′; SLC24A4 siRNA: sense, 5′-GAGGCUGGUAAUGAUUUCUAUGACG-3′, SLC24A4 siRNA: antisense, 5′-CGUCAUAGAAAUCAUUACCAGCCUC-3′.

Cells were co-transfected with lentiviral expression vectors and packaging plasmids (psPAX2 and pMD2.G) using Lipofectamine 2000 (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s protocols. Lentiviral supernatants were collected 48 h after transfection and filtered through a 0.45-micron filter. After transducing the lentiviruses into CRC cells, lentiviruses were added to medium containing polybrene (8 µg/ml). Stable cells were selected using puromycin for 1 week and detected with immunoblotting.

Total RNA extraction and quantitative real-time PCR

Total RNA extraction was extracted using Trizol reagent (Invitrogen, Waltham, MA, USA) and reverse transcribed with a reverse transcription kit (Takara, Shiga, Japan). Subsequently, qRT-PCR was performed with SYBR Green (Bio-Rad, Hercules, CA, USA) on an ABI PRISM 7900HT detection system (Applied Biosystems, Waltham, MA, USA). The primer sequences used in this study were produced by Generay Biotech (Shanghai, China). Relative mRNA expression levels were quantified by the 2^−ΔΔCt method. Experiments were performed in triplicate and GAPDH was used as an internal control. The primer sequences were as follows: SLC8A3 forward sequence 5′–3′: GAGGCCAAGAGGATAGCAGAG and reverse sequence 5′–3′: TGCTGCACTGACGGTGATG; SLC24A2 forward sequence 5′–3′: CCAAGGAGACTACCCGAAAGA and reverse sequence 5′–3′: CAGACAATGGCTAAGGCTATGAA; SLC24A3 forward sequence 5′–3′: CCCGTGCTGTCCCTTGTATT and reverse sequence 5′–3′: CCCCATCCGTTAATGCTGGT; SLC24A4 forward sequence 5′–3′: TGTGCCGAGACTCCGTGTA and reverse sequence 5′–3′: TGTTGACCGGGTTACCGTTTG; GAPDH forward sequence 5′–3′: GTCTCCTCTGACTTCAACAGCG and reverse sequence 5′–3′:ACCACCCTGTTGCTGTAGCCA.

Plasmids, reagents, antibodies

SLC8A3, SLC24A2, SLC24A3 and SLC24A4 were subcloned into PLVX. A DNA sequencing test was performed on all vectors to confirm their validity. Polybrene and puromycin were purchased from Sigma. This research used the following antibodies: anti-GAPDH (60004-1-lg, Proteintech, San Diego, CA, USA), anti-SLC8A3 (ANX-013, Alomone, Jerusalem, Israel), anti-SLC24A2 (ab192419, Abcam), anti-SLC24A3 (ab129776, Abcam, Cambridge, UK), anti- SLC24A4 (118992-1-AP, Proteintech, San Diego, CA, USA). All antibodies were obtained from commercial suppliers.

Colony formation assay

In each well, 1  × 10³ cells were placed in 6-well plates and left to incubate for 9 days in a culture medium containing 10% FBS. Following staining with 1% crystal violet, colonies per well were counted.

Cell counting kit-8 assay

We performed a cell counting kit-8 (CCK8) assay using a CCK8 kit according to the manufacturer’s instructions (Beyotime). A total of 5,000 RKO cells were evenly distributed in a 96-well plate with a total volume of 100 µL. After incubating for 24 h, 48 h, 72 h or 96 h, each well was supplemented with 10 µL of CCK8 reagent and incubated at 37 °C for 1 h. The plates were analyzed using a 450 nm wavelength.

Scratch wound healing assay

Cells were seeded in six-well plates and allowed to reach near-complete confluence within a 24-hour period. Subsequently, a wound was created in each well using a 10 µL pipette tip. Following this, the cells were subjected to multiple washes with PBS and subsequently incubated in serum-free culturing medium. The process of wound closure was observed at both 0 and 48 h using a Leica light microscope. The extent of cell migration post-wounding was assessed at 0 and 48 h.

Statistical analyses

Statistical analyses were conducted using GraphPad Prism 8.0 software. Multiple comparisons were performed utilizing the ANOVA test, while the assessment of the correlation between two variables was carried out using Spearman’s rank correlation coefficient. * P < 0.05, ** P < 0.01, *** P < 0.001.

Bioinformatic analyses revealed mRNA differential expression patterns of calcium extrusion-related genes in COAD

This study process is illustrated in Fig. 1. Firstly, using the RNA-Seq data from TCGA, we found the mRNA expression level of SLC8B1, SLC8A1, SLC8A2, SLC8A3, SLC24A3, and SLC24A4 was significantly decreased, while SLC24A2 was significantly increased in COAD compared with that in normal colon tissues (Fig. 2). Moreover, the similar results of SLC8B1, SLC8A1, SLC8A3, SLC24A4 and SLC24A2 were observed in GSE39582 (Fig. S1). These results suggest that the abnormal expressions of calcium extrusion-related genes might be associated with COAD progression.

Next, we analyzed the relationship between the expressions of calcium extrusion-related genes and tumor stage. In COAD, the expression of SLC24A1 gradually decreased from stage I to IV, while SLC8A2 expression increased from stage I to IV in COAD, and the expressions of SLC8B1, SLC8A1, SLC8A3, SLC24A2, SLC24A3, SLC24A4 and SLC24A5 did not show significant differences among different tumor stage (Fig. 3). Taken together, these results indicated that calcium extrusion-related genes, especially SLC24A1 and SLC8A2, may be involved in the development of COAD.

Meanwhile, the relationship between expression of calcium extrusion-related genes and clinicopathological features in COAD was systematically analyzed. Results indicated that mRNA expressions of the calcium extrusion-related genes were not associated with age, gender, T stage, N stage and M stage in COAD (Figs. S2A–S2G). Then, the expression of calcium extrusion-related genes was evaluated using the Kaplan–Meier Plotter database to predict the prognosis of COAD patients. However, the mRNA expression of calcium extrusion-related gene had no prognostic values for COAD patients (Figs. S3A–S3D). Overall, these results suggest the individual calcium extrusion-related gene had not a good predictive value in COAD.

Gene mutation signature, protein-protein interaction network and correlation analysis of calcium extrusion-related genes in COAD

To further investigate the molecular characteristics of calcium extrusion-related genes, gene mutation, interaction network between proteins and correlation analysis were performed. In cBioPortal database, the calcium extrusion-related genes were rarely mutated; SLC8A1 showed the highest mutation rate of 6%, following by SLC24A2 and SLC24A3 with the same mutation rate of 4% (Fig. S4A). In addition, the expression of calcium extrusion-related genes had no significant correlation with TMB in COAD (Fig. S4B). We analyzed the correlation between the calcium extrusion-related genes with the protein-protein interaction network formed by SLC8B1, SLC8A1, SLC8A2 and SLC8A3. The expression of them was significantly correlated in COAD (Figs. S4C–S4D). The co-expression gene interactions opined the function of the calcium extrusion-related genes mainly focused on cellular calcium ion homeostasis, active transmembrane transporter activity and calcium ion transport (Fig. S4E).

The association between mRNA expression of calcium extrusion- related genes and immune cell infiltration in COAD

We also studied the association between the expression of calcium extrusion-related genes and the infiltration of immune cells in COAD. CIBERSORT was used to determine the cellular composition of 22 immune cells, as shown in Fig. 4A and Fig. S5. Data indicated that the proportion of M0 macrophage cells were significantly different between the low-expression and high-expression groups divided by the median expression of SLC8A1, SLC8A2, SLC8A3, SLC24A2 and SLC24A3, respectively. There was a notable disparity in the percentage of M2 macrophage cells between the low and high categories of SLC8B1, SLC8A1, SLC24A2, SLC24A3, and SLC24A4. Additionally, there was a notable disparity in the proportion of resting CD4 memory T cells between the low-expression and high-expression categories of SLC8A2, SLC24A2, and SLC24A4. Furthermore, the distribution of activated CD4 memory T cells exhibited a significant contrast between the two groups of SLC8A1 and SLC24A3, respectively. Taken together, these results suggested that the immune cell infiltration was closely related with the expressions of calcium extrusion-related genes, which needs to be clarified in the future.

The association between the mRNA expression of calcium extrusion-related genes and immune checkpoint and immunotherapy response

Numerous studies reported that immune checkpoint molecules were closely associated with the response to immunotherapy in cancers. Therefore, we examined the expression levels of immune checkpoint molecules in different groups. Our results indicated that PD-1, PD-L1, PD-L2 and CTLA-4 expression level have significant differences between the low-and high-expression of almost all calcium extrusion-related genes, except for SLC8A2 (Fig. 4B). The potential of anti-PD1 and anti-CTLA4 immunotherapy response in the low-and high-expression groups of calcium extrusion-related genes was examined using TIDE. Our results showed that the SLC8B1-high expression group might respond better to immunotherapy, while the low expression group of SLC8A2, SLC8A3, SLC24A2 or SLC24A3 might respond better to immunotherapy (Fig. 4C). Additionally, the analysis of functional enrichment revealed that SLC8B1 and its co-expressed genes were involved in T cell activation, leukocyte intercellular adhesion, and the positive regulation of cell adhesion (Fig. S6). Altogether, these results demonstrated that the expressions of calcium extrusion-related genes might be associated with immunotherapy responses in COAD.

Construction of calcium extrusion-related genes risk signature

Due to the failure prediction for prognosis of COAD using the individual calcium extrusion-related gene, we established a novel risk model based on expression level of calcium extrusion-related genes to predict prognostic outcomes of COAD patients. The TCGA-COAD dataset served as the training dataset, while three GEO datasets (GSE29623, GSE39582, and GSE17536) were utilized as the validation datasets. The results indicated that high-risk group had a worse prognosis for patients with COAD in training dataset (Figs. 5A and 5C) and validation datasets (Figs. 5B, 5D and Figs. S7A and S7B). Then, prognostic risk models were used in predicting 1-, 3-and 5-year overall survival, ROC curves were plotted with AUC values of 0.584, 0.623 and 0.636, respectively (Fig. 5E). Corresponding AUC values of 0.706, 0.675 and 0.704 were observed in the GSE17536 cohort; AUC values of 0.827, 0.637 and 0.792 were observed in the GSE29623 cohort, and AUC values of 0.564, 0.587 and 0.594 were observed in GSE39582 cohort (Fig. 5F and Fig. S7C). The clinical characteristics of the low-risk and high-risk groups were also analyzed and we found that age and T stage among the two risk groups were significantly different (Table 1). Altogether, these results indicated that the risk score model is effective in predicting the prognosis of COAD patients and patients with high risk have worse survival.

Mutation profile of COAD patients in high- and low-risk groups

Next, we mapped the mutation landscapes of COAD in both the high-risk and low-risk groups. In the high-risk group, APC (70%), TP53 (56%), TIN (53%), KRAS (49%) and PIK3CA (35%) were the top five mutated genes, while APC (70%), TP53 (54%), TIN (52%), KRAS (37%) and MUC16 (34%) were the most frequently mutated genes in low-risk group (Fig. 6A). These gene mutations are shown in Fig. S8. Both groups shared APC, TP53, TTN, KRAS, PIK3CA, MUC16, SYNE1, FAT4, OBSCN, RYR2, ZFHX4, DNAH5, CSMD3, FLG and DNAH11 mutations at high frequencies (Fig. 6B). Interestingly, the frequency of KRAS mutation were higher in patients of high-risk group, whereas MUC16 had a lower frequency of mutation in the high-risk group (Fig. 6C). The tumor mutation burden (TMB) was an important prognostic indicator for COAD patients. Therefore, we evaluated the relationship between risk score and TMB. The results demonstrated a significant negative correlation between risk score and TMB in COAD (Fig. 6D). The overall survival rate of high TMB patients was lower than that of low TMB patients (Fig. 6E). Moreover, in high TMB subgroup, high-risk group had worse prognosis than low-risk group (Fig. 6F). In conclusion, the above results indicated combination of risk score and TMB might be a more effective biomarker in predicting COAD prognosis.

Calcium extrusion-related risk score was associated with immune signatures and immunotherapy responses in COAD

Using the CIBERSORT algorithm, we analyzed the proportion of 22 immune cells in the high-and low-risk groups and assessed relationship between risk score and tumor immune microenvironment (TIME). A significant increase in T cells regulatory (Tregs) and M0 macrophages was observed in the high-risk group, while the proportion of B cells naive, T cell CD4 memory resting, and M2 macrophages were remarkably decreased in the high-risk group (Fig. 7A, Fig. S9A). Similarly, risk score was positively correlated with majority signatures of Tregs (Fig. 7B). The relationship between the risk score and other immune cells signatures have been shown in Fig. S9B. Consisted with previous studies, we found that high-risk groups had lower immune score than low-risk group (Fig. 7C). Collectively, these findings suggest that the increased infiltration of Tregs-involved tumor immunosuppressive microenvironment might contribute to the poor survival for patients with COAD.

In addition, our data showed that six inhibitory immune checkpoints and eight stimulatory immune checkpoints were altered in high-risk group. The expressions of inhibitor immune checkpoints including CD160, CD244, IL10, KIR2DL1 and KIR2DL3, were declined in high-risk groups. The expressions of stimulatory immune checkpoints including ICOSLG, PVR, TNFSF13 and TNFSF15 were increased in high-risk groups (Fig. 7D). In patients with COAD, the effectiveness of risk scoring in predicting the outcomes of immunotherapy was evaluated by the TIDE algorithm. Our findings indicated that the high-risk group exhibited a lower response rate to immunotherapy in comparison to the low-risk group, but this disparity did not achieve statistical significance (Fig. 7E). Subsequently, the PRJEB25780 cohort, comprising of 45 individuals diagnosed with advanced stomach cancer receiving PD-L1 inhibitor therapy, was employed to determine the effectiveness of the risk signature associated with calcium extrusion in accurately predicting the outcomes of immunotherapy. The low-risk group showed a considerably higher immunotherapy response rate compared to the high-risk group (Fig. 7F). As shown in Fig. 7G, the immunotherapy response rate was significantly higher in the subgroup with a high-immune score than those with a low-immune score. Notably, within the high-immune score subgroup, the low-risk group demonstrated a higher immunotherapy response rate (63.00%) compared to the high-risk group (58.00%) (Fig. 7H). These findings indicated that risk score alone or combined with risk score and immune score could predict the immunotherapy responses in patients with COAD.

TMB has been shown to serve as a valuable biomarker for the purpose of selecting immune checkpoint blockade in various cancer types. The subgroup with a high TMB exhibited a slightly higher immunotherapy response rate (51.00%) compared to the subgroup with a low TMB (44.00%), although this disparity did not reach statistical significance (Fig. S9D). Moreover, our findings indicated that the rate of immunotherapy response in the low-risk group (56.00%) was comparatively higher than that observed in the high-risk group (46.00%) within the high tumor mutational burden (TMB) subgroup, despite the lack of statistical significance (Fig. S9E). The phenotype linked to microsatellite instability-high (MSI-H) is a unique category of tumors that display an increased vulnerability to immunotherapy. Notably, A significantly higher immunotherapy response rate was observed in the MSI-H subgroup compared to the MSS subgroup (Fig. S9F). Moreover, our findings indicated that the immunotherapy response rate in the low-risk group (68.00%) was comparatively higher than that observed in the high-risk group (61.00%) within the MSI-H subgroup, although this difference did not reach statistical significance (Fig. S9G). Taken together, these findings suggested that combined with risk score, immune score, TMB, or microsatellite status might be a potential strategy to predict immunotherapy response in patients with COAD.

Risk score predicts therapeutic benefits in COAD

In order to assess the efficacy of risk score as a biomarker for prognosticating drug response in COAD, the IC50 value of 198 drugs was deduced for TCGA-COAD patients. The results showed that the IC50 values of BMS-754807, Dabrafenib, Daporinad, PLX-4720, SB216763, UMI-77, MIM1, LJI308, PD0325901, and Telomerase inhibitor IX in the low-risk group were significantly lower than the high-risk group (Fig. 8A and Fig. S10). These results demonstrated that patients in the low-risk group were relatively sensitive to these agents. However, patients in high-risk group exhibited a higher sensitivity to AZD6482, KU-55933, PF-4708671, Acetalax, GSK2699962A, RO-3306, NU7441, Ribociclib, Sapitinib, and Taselisib (Fig. 8B and Fig. S10). In conclusion, the results indicated that patients in different group might be sensitive to various drugs, which was of great importance for the development of individualized treatment and clinical medication guidance. Tables 2 and 3 provided comprehensive details about the top 10 vulnerable medications for both the low-risk and high-risk subgroup.

Function analysis of the previously unreported model genes SLC8A3, SLC24A2, SLC24A3 and SLC24A4 in CRC

Previous studies have shown that SLC8A1, SLC8A2, and SLC8B1 contribute to the initiation and progression of cancer (Muñoz et al., 2015; Qu et al., 2017; Pathak et al., 2020). Therefore, we examined previously unreported model genes for their functional roles (SLC8A3, SLC24A2, SLC24A3 and SLC24A4) in CRC cells. We first investigated the mRNA expression level of SLC8A3, SLC24A2, SLC24A3 and SLC24A4 in different types of CRC cells. Then, we assessed cell viability and colony formation in those overexpression or knockdown RKO cell, which exhibited a relative median expression level of SLC8A3, SLC24A2, SLC24A3 and SLC24A4 (Fig. S11A). To validate the efficiency of cell transfection, western blots were conducted, and the results were presented in Fig. 9A and Fig. S11B. The results indicated that the downregulation of SLC8A3 or SLC24A4 expression promoted RKO cells growth in vitro (Figs. 9B–9C), while the overexpression of SLC8A3 or SLC24A4 inhibited RKO cells growth in vitro (Figs. S11C–S11D). The knockdown or overexpression of SLC24A2 and SLC24A3 unexpectedly exhibited a weak effect on CRC growth (Figs. 9B–9C and Figs. S11C–S11D). Moreover, we assessed CRC cell migration in those overexpression or knockdown RKO cell by scratch wound healing assay. The results revealed that the migration ability of RKO cells was markedly increased when SLC8A3 or SLC24A4 was downregulated. By contrast, there was a significant decrease in the migration capability of CRC cells in SLC8A3 or SLC24A4 overexpression cells (Fig. 9D and Fig. S11E). However, cell migration was not affected by the knockdown or overexpression of SLC24A2 and SLC24A3 (Fig. 9D and Fig. S11E). The results showed that SLC24A2 and SLC24A3 did not affect CRC cell proliferation and migration. Altogether, these results suggest that SLC8A3 and SLC24A4 may be involved in CRC growth and metastasis as tumor suppressor genes.