Nav1.6 drives colorectal cancer proliferation and invasion through MAPK signaling pathway

Abstract

BACKGROUND

Voltage-gated sodium channels (VGSCs, or Navs) are highly expressed in various tumors and play a critical role in tumor metastasis and invasion.

AIM

To identify Nav1.6-associated cancer genes through bioinformatics analysis and experimental validation, with the goal of determining the role of Nav1.6 in colorectal cancer (CRC) metastasis.

METHODS

The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) data were analyzed using weighted correlation network analysis (WGCNA) and Venn analysis to identify Nav1.6-associated genes in CRC. siRNA, real-time PCR, and western blotting were employed to validate the Nav1.6-associated cancer genes and signaling pathways identified in CRC. Cell counting kit-8 and Transwell migration assays were used to assess the proliferation and migration of CRC cells.

RESULTS

The analysis of TCGA and GEO datasets, along with WGCNA, identified 575 differentially expressed genes associated with SCN8A (Nav1.6) in CRC, which were particularly enriched in MAPK signaling pathways. Tissue microarray analysis of surgical samples revealed elevated Nav1.6 levels in CRC tissues, which were predominantly in the cytoplasm and nucleus rather than in the membrane. Cytoplasmic Nav1.6 expression increased with T stage increases, consistent with the TCGA findings. SCN8A knockdown in colon tumor cells significantly reduced cell proliferation and invasion and downregulated key proteins in the RAF-MAPK pathway.

CONCLUSION

These findings suggest that Nav1.6 promotes CRC cell proliferation and invasion which is related to the MAPK signaling pathway.

Key Words: Nav1.6; Colorectal cancer; Weighted correlation network analysis; MAPK signaling; Cancer proliferation

Core Tip: Voltage-gated sodium channels are prominently expressed in various tumors and play a crucial role in tumor metastasis and invasion. Our research identified 575 differentially expressed genes related to Nav1.6 in colorectal cancer (CRC), particularly enriched in MAPK signaling pathways. Tissue microarray revealed heightened Nav1.6 levels in CRC tissues. The cytoplasmic expression of Nav1.6 increases with T staging. Knockdown of SCN8A in colon tumor cells significantly decreased cell proliferation and invasion, along with downregulation of key proteins in the RAF-MAPK pathway. These findings suggest that Nav1.6 enhances CRC cell proliferation and invasion through modulation of the MAPK signaling pathway.

INTRODUCTION

Colorectal cancer (CRC) is the third most common malignancy globally, with the second highest mortality rate[1]. Its incidence has steadily increased in recent years[1]. In 2023, it is estimated that approximately 153020 individuals will be diagnosed with CRC, resulting in 52550 deaths. This includes 19550 new cases and 3750 deaths among individuals under 50 years of age[1]. CRC metastasis refers to the spread of cancerous cells from the colon or rectum to distant organs. This complex process plays a pivotal role in CRC progression and significantly impacts patient prognosis and treatment outcomes[2]. Metastasis typically involves the detachment of cancer cells from the primary tumor, the invasion of adjacent tissues, and entry into the blood or lymphatic vessels, where they may establish secondary tumors in organs such as the liver, lungs, and peritoneum[2].

CRC is staged according to the tumor-node-metastasis classification, which takes into account the extent of the primary tumor (T), regional lymph node involvement (N), and the presence of distant metastasis (M). The "T" stage refers to the depth of invasion of the primary tumor into the surrounding tissues and organs. It is classified as follows: T1: The tumor invades the submucosa (the second layer of the colon wall). T2: The tumor invades the muscularis propria (the muscle layer). T3: The tumor invades through the muscularis propria into the subserosa or into the non-peritonealized pericolic or perirectal tissues. T4: The tumor invades the serosa (outermost layer of the colon) or adjacent structures, such as nearby organs or tissues[3].

The metastatic potential of CRC is influenced by various factors, including the molecular characteristics of the tumor, the tumor microenvironment, and the host’s immune response. A comprehensive understanding of the mechanisms driving CRC metastasis is essential for developing effective therapeutic strategies, particularly for advanced-stage disease, with the ultimate goal of improving patient survival outcomes[4,5].

Voltage-gated sodium channels (VGSCs or Navs) are large molecular complexes embedded in the cell membrane, composed of an α subunit that forms the ion-conducting pore and one or more associated β subunits[6]. The VGSC α subunit family includes nine isoforms-Nav1.1 through Nav1.9-encoded by the SCN1A-SCN11A genes[6]. VGSCs are predominantly expressed in excitable cells, such as neurons and cardiac myocytes, where they are essential for the initiation of action potentials and the transmission of electrical signals[7]. Specifically, Nav1.1, Nav1.2, and Nav1.3 are primarily expressed in the central nervous system; Nav1.4 is found in skeletal muscle; Nav1.5 in cardiac muscle; and Nav1.6, Nav1.7, Nav1.8, and Nav1.9 in the peripheral nervous system, including the enteric and sensory systems[6]. Recent studies have demonstrated that VGSCs are expressed in various cancer cells, including breast, cervical, colorectal, and prostate cancers, as well as in melanoma, neuroblastoma, and non-small cell lung cancer[8]. The aberrant overexpression of VGSC α subunits in cancer cells has been implicated in promoting tumor cell migration and invasion[9]. Notably, Nav1.1, Nav1.2, Nav1.3, Nav1.4, and Nav1.9 show heightened expression in ovarian cancer, non-small cell lung cancer, and prostate cancer[8]. The Nav1.7 α subunit has also been shown to facilitate the progression of gastric cancer through the MACC1-mediated upregulation of NHE1[10]. Furthermore, the expression of Nav1.5 is closely associated with poor prognosis in several cancers, including breast cancer[11], non-small cell lung cancer[12], ovarian cancer[13], and prostate cancer[14].

Nav1.6, encoded by the SCN8A gene, is expressed in a variety of solid tumors, including breast cancer, cervical cancer, non-small cell lung cancer, ovarian cancer, and prostate cancer[8]. The elevated expression of Nav1.6 has been shown to enhance the metastasis and invasiveness of cervical cancer cells, and the specific inhibition of Nav1.6 currents significantly reduces their invasive potential[15]. In our previous study, we demonstrated for the first time that Nav1.6 is aberrantly overexpressed in more than 50% of CRC tumor tissues[16]. Furthermore, the high expression of Nav1.6 correlated positively with vascular invasion and lymph node metastasis in patients with CRC[16]. These findings suggest that Nav1.6 may contribute to CRC metastasis; however, the underlying mechanisms remain to be fully elucidated. A deeper understanding of the role of Nav1.6 in CRC lymph node metastasis and the molecular pathways through which elevated Nav1.6 expression promotes this process is essential for advancing our comprehension of CRC metastasis.

With advancements in sequencing technologies, large-scale gene co-expression analysis has become a cornerstone in understanding disease mechanisms and has become increasingly pivotal in life sciences research. This study aimed to investigate how Nav1.6 potentially promotes lymph node metastasis by leveraging data from The Cancer Genome Atlas (TCGA) database. The analysis revealed a positive correlation between elevated Nav1.6 expression and the CRC T stage. This finding was further validated through tumor tissue microarray (TMA) analysis of CRC patients. Using weighted correlation network analysis (WGCNA) of TCGA data, we identified several differentially expressed genes (DEGs) associated with Nav1.6, including members of the MAPK pathway. Further examination of CRC single-cell sequencing data from the Gene Expression Omnibus (GEO) database identified epithelial cells as the predominant cell type in tumor tissues, where Nav1.6, MAP3K10, and other tumor-related genes were prominently co-expressed. The siRNA-mediated knockdown of Nav1.6 in colorectal cells confirmed these findings. Collectively, these results provide novel insight into the molecular mechanisms by which Nav1.6 may contribute to colorectal carcinogenesis. They not only underscore the role of VGSCs in tumor metastasis but also open new avenues for exploring the molecular basis of CRC metastasis.

MATERIALS AND METHODS

Initial screening of microarray datasets

The microarray datasets utilized in this study were obtained from the NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/) and the TCGA database. Human CRC single-cell sequencing data (GSE201348) were selected for analysis. For each dataset, the differential expression analysis of tumor tissue versus normal tissue was performed using the DESeq2 program to generate a list of DEGs. A threshold of P < 0.05 and a fold change of > 2 were applied. Several web-based tools were employed to identify SCN8A-related gene signatures, including Cluster Heatmap, Principal Component Analysis, GO/KEGG analysis, Venn Diagram Generator, Co-expression Correlation Analysis, and Protein-Protein Interaction Network (https://cloud.oebiotech.com), to identify Nav1.6-related DEGs.

Single-cell sequencing data analysis

Single-cell transcriptomic data from GSE201348 were analyzed using the Seurat package in R (version 4.2.2). Batch effects across samples were mitigated using the batch-balanced K-nearest neighbors (BBKNN) method from the BBKNN package. Non-linear dimensionality reduction was performed using uniform manifold approximation and projection. Cell types were classified through unsupervised clustering based on the expression of canonical marker genes. The expression levels of specific genes across cell clusters were visualized with bubble plots, where the size of each bubble represented the proportion of cells expressing the gene, and the color indicated the scaled average expression level.

Weighted gene co-expression network analysis

WGCNA was conducted using the WGCNA package in R[17,18]. An unsigned co-expression network was constructed based on an adjacency matrix, and genes with highly similar co-expression patterns were grouped using average linkage hierarchical clustering of the topological overlap. The Dynamic Hybrid Tree Cut algorithm was applied to segment the hierarchical clustering tree, with modules defined as the branches resulting from tree cutting. Moderately large and distinct modules were identified by setting the minimum module size to 30 genes and applying a merging height threshold of 0.25. The resulting modules were assigned random color labels. Each module was characterized by its first principal component, termed the module eigengene, which represents the scaled (standardized) expression profile of the module. Correlations between module eigengenes and clinical phenotypes, such as estimated glomerular filtration rate and interstitial fibrosis, were then analyzed.

Clinical sample collection

The study enrolled patients diagnosed with stage II-III colon cancer (American Joint Committee on Cancer 8^th edition), with tumor lesions greater than 2 cm in size, who had not received radiotherapy or chemotherapy prior to laparoscopic surgery. A total of 90 patients aged 34-81 years participated, including 61 males and 29 females. Among these, 19 were classified as having stage II, 69 as having stage III, and two as having stage I. Lymph node metastasis was present in 55 cases and absent in 25 cases, while perineural invasion was observed in 22 cases and absent in 68 cases. Following surgery, specimens were promptly collected for in vitro analysis. A 0.5 cm × 0.5 cm segment of CRC tissue, isolated from the luminal surface of the intestine, was used for further examination. A 0.5 cm × 0.5 cm section of intestinal wall tissue (primarily the mucosa and submucosa), located 5 cm from the tumor margin, served as the normal tissue control. All clinically obtained samples were stored at -86 °C within 2 minutes of isolation. The remaining tumor tissue was preserved in 4% formalin for subsequent hematoxylin and eosin (H&E) staining.

Immunohistochemistry

A TMA was constructed using 1 mm² tissue cores from 90 patients, utilizing a TMArrayer (Pathology Devices). Areas with high tumor content were selected for inclusion in the TMA block based on H&E-stained sections from the entire tissue block, ensuring the presence of CRC cells in each sample. TMA sections, 4 µm thick, were cut using a microtome (Leica RM2125 RTS), mounted on SuperFrost™ Plus slides (Thermo Fisher Scientific), and air-dried overnight at room temperature. Nav1.6 staining was performed following previously established protocols[16]. After staining, the slides were scanned using a Vectra 3.0 scanner, and images were visualized with Phenochart v.1.0.8 (AKOYA Bioscience). Optimal scanning protocols were determined by adjusting the exposure times for each L cube at × 10 and × 20 magnifications. Whole-slide scans were conducted using a × 10 objective lens, while high-resolution images of individual core tissue regions were captured with a × 20 objective lens. TMA core images were analyzed using HALO v2.2 software (Indica Labs). The random forest tissue classifier was trained on multiple tissue areas to identify tissue and slide regions. Positivity for each marker was assessed based on staining intensity, ensuring a minimum 10-fold signal-to-noise ratio, staining patterns that are consistent with the manufacturer’s specifications and previous literature, and comparable staining with control tissues. The analysis settings specific to each staining category were optimized by adjusting the thresholds, and all samples were analyzed.

Immunofluorescence

Cut the paraffin embedded tissue sample into 3-5 μm thickness using a slicer and adhere it to a glass slide. According to the manufacturer's instructions, use TSA fluorescence kit (Melady)^® Biosciences is used for target protein labeling. These slices underwent a one-hour baking process at 60 °C, followed by dewaxing in xylene and graded ethanol series (100%, 95%, 80%, and 70%), and rinsing in ddH₂O. Use a microwave oven to perform thermally induced antigen repair in EDTA buffer at pH 9.0. Quench endogenous peroxidase activity with 3% hydrogen peroxide, then seal the slices with a blocking solution. After 20 minutes, shake off the blocking solution and incubate the slices with primary antibody overnight at 4 °C. The slices were incubated with HRP secondary antibodies suitable for the species at room temperature for 20 minutes, washed off with PBS, and then incubated with appropriate fluorescent dyes for TSA staining. The first antibodies we use are CD31 (Abcam, ab231436, 1:50), KI67 (Abcam, ab92742, 1:500), DCX (Santa Cruz Biotechnology Company, SC-271390, 1:100). After completing TSA staining of the first antibody, repeat the process for the second antibody. After staining all targets, wash the slices with PBS and counterstain with DAPI. Finally, install the slice using a suitable mounting medium, perform full slide scanning imaging using a multi-channel fluorescence scanner (3DHISTECH, Pannoramic MIDI), and view the image using SlideViewer software.

Cell culture and transfection

All cell lines were obtained from the Shanghai Institute of Cell Biology. CaCo2 cells were cultured in MEM supplemented with 10% FBS and 1% penicillin-streptomycin; DLD-1 cells were cultured in 1640 supplemented with 10% FBS and 1% penicillin-streptomycin; and HT-29 cells were cultured in McCoy’s 5A medium supplemented with 10% FBS and 1% penicillin-streptomycin. All cell culture reagents were purchased from Gibco (Thermo Fisher Scientific). Cells were transfected with the si-hNav1.6 plasmid (GenePharma, China) using Lipofectamine 3000 (Thermo Fisher Scientific) diluted in Opti-MEM (Gibco) to knock down the expression of Nav1.6. The plasmid and transfection reagent were diluted separately, incubated for 5 minutes, mixed, and allowed to stand for 20 minutes at room temperature. The transfection mixture was then added to DLD-1 cells seeded in 6-well plates and incubated for 5 hours. Afterward, the medium was replaced with a complete growth medium. Cells were collected for further analysis after 24 hours.

Quantitative RT-PCR assay

Total RNA was extracted using Trizol (Vazyme) following the manufacturer’s instructions. Next, 1 μg of RNA was reverse-transcribed into complementary DNA (cDNA) using HiScript III RT SuperMix for quantitative PCR (qPCR) (Vazyme). qPCR was performed on an ABI 7500 System (Invitrogen Life Technologies) using human-specific primers (Table 1). The PCR cycling conditions were as follows: Initial denaturation at 95 °C for 30 seconds, followed by 40 cycles of denaturation at 95 °C for 5 seconds, and annealing/extension at 60 °C for 30 seconds. A melting curve analysis was conducted using SYBR Green (Vazyme) to confirm the specificity of the PCR products. The samples were analyzed in duplicate to ensure the reliability of the qRT-PCR data, and the results were based on the average of the duplicate values for each sample. Relative gene expression was calculated using the 2^-ΔΔCt method, using GAPDH as the internal control and untreated samples from control mice as the calibrator. PCR data were normalized to the mean values of the control group to minimize variability. The primers used are listed in Table 1.

Table 1 Primers for quantitative RT-PCR assay.

Primer name	Sequences (5’-3’)
hGAPDH-F	GTCTCCTCTGACTTCAACAGCG
hGAPDH-R	ACCACCCTGTTGCTGTAGCCAA
hNav1.6-F	GGATTGAGACCATGTGGGACTG
hNav1.6-R	ATCTGTGGCAGCCAGGTTGTCT

Western blots

After transfection, cells were harvested by centrifugation, washed with ice-cold 1 × PBS, and lysed using radioimmunoprecipitation assay (RIPA) buffer (Beyotime). The lysed proteins in the supernatant were collected following centrifugation and mixed with 5 × SDS loading buffer (Beyotime). The mixture was then boiled and loaded onto SDS-polyacrylamide gels. Protein concentrations were quantified using a BCA assay (Beyotime) to ensure equal loading. Proteins were separated by electrophoresis, transferred to PVDF membranes, and blocked with 10% BSA. The membranes were incubated overnight at 4 °C with primary antibodies against Nav1.6 (Abcam ab65166), MAP3K10 (Bioworld BS61034), phospho-c-Raf (Cell Signaling Technology 9427), ERK1/2 (Abcam ab184699), MEK1/2 (Abcam ab178876), Raf1 (Abcam ab137435), MEK1/2 (phospho-S221) (Abcam ab278564), GAPDH (Abcam ab8245), and β-actin (AiFang Biological AF10678). After washing, the membranes were incubated with HRP-conjugated secondary antibodies and subjected to chemiluminescent detection. Image analysis for quantification was performed using Image Lab software, with background subtraction applied.

Transwell assay

Cell invasion ability was assessed using Transwell assays with 24-well chambers containing 8 μm polycarbonate membranes (Millipore, MA, United States). The filters were pre-coated with 500 ng/mL Matrigel solution (BD Biosciences, United States) for invasion assays and incubated for 4 hours at 37 °C. The lower chamber was placed in a 24-well plate containing 500 μL of medium with 10% FBS, and 1 × 10⁵ cells in 200 μL of serum-free medium were seeded into the upper chamber. After 18 hours of incubation at 37 °C, non-invading cells on the upper membrane surface were gently removed by scraping. The invading cells on the membrane were fixed with methanol, stained with 0.5% crystal violet (Beyotime), and photographed. Invasive cells were counted in fewer than five random 200 × microscopic fields per well using a Nikon Eclipse Ti inverted research microscope. A cell migration assay was conducted in parallel under the same conditions, except that the chambers were not pre-coated with Matrigel.

Wound healing assays

Cells (1 × 10⁵) were seeded in 6-well plates. After 16 hours, the complete medium was replaced with a low-serum medium (2%). A uniform wound was created in each well by scratching the surface with a 10 μL pipette tip once the cells reached 90% confluence. The cells were gently washed twice with PBS to remove any detached cells, and a serum-free medium was added. Multiple reference marks were made at the center of the denuded surface to ensure consistency in the wound area. Scratch wounds were photographed using an inverted microscope at 0, 12, and 24 hours. Image analysis was performed using AxioVision Rel. 4.8 software to assess the migratory capacity of the cancer cells. The experiment was repeated three times to ensure reproducibility.

Cell counting kit-8 assay

DLD-1 cells (5 × 10³ cells/well) were seeded in a 96-well plate and incubated with complete medium containing the specified concentration of CNB at 37 °C for 12, 24, and 48 hours. Cell viability was assessed using the cell counting kit-8 (CCK-8) assay (TargetMOL, C0005), following the manufacturer's instructions. The optical density at 450 nm was measured using a microplate reader (Thermo Fisher).

Statistical analysis

The Student's t-test, one-way ANOVA, and two-way ANOVA without correction were used to assess statistical significance between groups (^aP < 0.05, ^bP < 0.01, ^cP < 0.001, and ^dP < 0.0001). Data are presented as mean ± SEM. All analyses were conducted using GraphPad Prism 10. The statistical review of the study was performed by a biomedical statistician.

RESULTS

Nav1.6 is highly expressed in CRC and is positively correlated with T stage

In our previous study, we analyzed the mRNA expression of VGSCs in CRC tissue using collected CRC samples and found that SCN8A was significantly overexpressed in both tumor tissues and metastatic lymph nodes, with a strong correlation with the lymph node stage[18]. We extracted transcriptomic sequencing data for rectal adenocarcinoma (READ) and colon adenocarcinoma (COAD) from the TCGA database and analyzed SCN8A expression to further explore the relationship between Nav1.6 expression and tumor metastasis in CRC. Our results showed a significant upregulation of SCN8A in both READ (Figure 1A) and COAD tumor tissues (Figure 1B). Additionally, SCN8A expression increased with advancing T stages (Figure 1C), further supporting our earlier findings that the aberrant overexpression of Nav1.6 may be linked to CRC metastasis.

Figure 1 Expression of SCN8A (Nav1.6) in colorectal cancer cells. A and B: Expression of Nav1.6 in rectal adenocarcinoma and colon adenocarcinoma. Each dot represents a patient; C: Expression of Nav1.6 in different stages of colorectal cancer; D: Representative images of immunohistochemical staining of Nav1.6 in cancer tissues and adjacent normal tissues from 90 colorectal cancer (CRC) patients; E: Heatmap showing the expression level of Nav1.6 in tumor and normal colon tissues after staining with Nav1.6 antibody. The orange-red color indicates a high quantitative value, while white represents a low quantitative value; F: H-Score of Nav1.6 signal intensity in 97 cancer and adjacent tissues; G-I: Proportion of DAP-positive cells (G), positive nuclei (H), and positive cytoplasm (I) in cancer and normal colon tissues; J: Cytoplasmic expression of Nav1.6 in different stages of CRC. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001; ^dP < 0.0001; READ: Rectal adenocarcinoma; COAD: Colon adenocarcinoma; TCGA: The Cancer Genome Atlas.

We performed TMA (immunohistochemical staining) analysis of surgical resection samples from 90 CRC patients to validate these findings, which included adjacent normal tissues and tumor tissues (2 T1, 19 T2, and 69 T3). CRC typically progresses through several stages, reflecting the extent of tumor growth and spread within the colon or rectum as well as to other parts of the body. This study focused on stages T1, T2, and T3 of stage II. The results showed significantly elevated Nav1.6 expression in tumor tissues (Figure 1D). Subsequent analysis using HALO TMA analysis software included heatmap generation, H-SCORE analysis, and quantification of the proportion of positive cells, as well as positive nuclear and cytoplasmic staining. These analyses revealed a notably higher proportion of Nav1.6-positive cells in tumor tissues compared to normal colorectal tissues (Figure 1D). Importantly, Nav1.6 expression was significantly elevated in both the cytoplasm and nucleus of CRC tissues (Figure 1E-I), suggesting not only overexpression but also the aberrant subcellular localization of Nav1.6. Further analysis showed a positive association between cytoplasmic Nav1.6 expression and local T staging (Figure 1J), confirming the results from the TCGA database analysis and highlighting the link between cytoplasmic Nav1.6 expression and local staging in CRC.

Identification and function enrichment of Nav1.6-related DEGs in CRC tissue

We performed WGCNA on DEGs (fold change ≥ 1.5) derived from the RNA-seq data of CRC tumor and normal samples in TCGA to identify Nav1.6-related tumor genes in CRC. A total of 575 DEGs, comprising 230 upregulated genes and 345 downregulated genes, were selected from the intersection of SCN8A-correlated genes in COAD, SCN8A-correlated genes in READ, and DEGs from tumors versus normal comparisons. These genes were used to construct a co-expression network (Figure 2A and B). Genes with low expression variability (standard deviation ≤ 0.5) were excluded, resulting in 293 genes for further analysis. A soft-thresholding power was applied to the network topology to assess scale independence and mean connectivity. After screening, a soft-threshold power of 20 was selected, achieving an approximate scale-free topology with a scale-free topology fit index ≥ 0.9, as shown in Figure 2C. This selected power was used to construct a weighted co-expression network model, dividing the 293 genes into five modules, with the gray module representing genes that could not be assigned to any specific module (Figure 2D). A heatmap of the selected 293 genes demonstrated relatively high modular independence (Figure 2E). Among these, 167 genes (turquoise module in Figure 2F) were positively correlated with the T stage, while 98 genes (blue module in Figure 2F) were negatively correlated with the T stage. Therefore, these two modules were identified as key modules in SCN8A-related CRC progression. Further analysis revealed significant correlations between the molecular module and gene significance (GS) for specific traits. Notably, significant correlations were observed between the green and magenta modules and GS for the T stage, as shown in Figure 2G.

Figure 2 Weighted correlation network analysis of the gene set associated with Nav1.6 expression. A: Screening of SCN8A-correlated genes from tumors vs normal differentially expressed genes (DEGs) in the colon adenocarcinoma and rectal adenocarcinoma databases of The Cancer Genome Atlas; B: Volcano plot showing the upregulated genes (red) and downregulated genes in the 575 intersecting DEGs; C: Selection of soft-threshold power. Panel C shows the scale-free fit index (left) and average connectivity (right) vs soft-threshold power. A power of 12 was selected as the fit index curve flattens at high values (> 0.9); D: Cluster dendrogram and module assignment from weighted correlation network analysis. Genes were clustered based on dissimilarity measures (1-TOM). The branches represent modules of highly interconnected genes. The dynamic tree cut colors indicate modules identified by the Dynamic TreeCut method. Some modules were merged based on similarity, and the final merged modules were used for subsequent analysis; E: Heatmap representing the topological overlap matrix (TOM) of all genes in the analysis. The degree of overlap is indicated by color shading, with darker colors representing higher overlap and lighter colors indicating lower overlap. The gene dendrogram is shown on the left, and module assignments are shown on the top; F: Correlation of module eigengenes with clinical and pathological traits. Each row corresponds to a module eigengene, and the columns represent clinical traits. The values in the cells are presented as "Pearson r (P value)" and are color-coded based on the direction and magnitude of the correlation (red = positive correlation; blue = negative correlation). The 16 co-expressed transcript modules are shown along with their respective correlation parameters with clinical and pathological traits; G: Distribution of genes positively correlated with T stage in colorectal cancer (blue). READ: Rectal adenocarcinoma; COAD: Colon adenocarcinoma.

We performed KEGG enrichment analysis on the 167 T stage-related genes identified through WGCNA to further identify Nav1.6-related functional genes. The bubble chart from the KEGG analysis revealed that the upregulated genes were predominantly involved in MAPK signaling, nicotine addiction, and morphine addiction (Figure 3A and B). In contrast, the downregulated genes were mainly associated with pathways in cancer, renin secretion, and protein digestion and absorption. Given the well-established role of MAPK signaling pathways in CRC metastasis, we further examined the expression of KEGG-enriched MAPK-related genes in both normal and CRC tissues using data from TCGA and GEO. The results showed that the expression of MAP2K1, MAP3K8, and MAP3K10 was significantly elevated in colon cancer tissues in the READ and COAD datasets, with MAP3K10 showing high expression in both colon and rectal cancer tissues (Figure 3C-E).

Figure 3 KEGG enrichment analysis of Nav1.6-related differentially expressed genes from weighted correlation network analysis-enriched molecules. A and B: KEGG pathway analysis of differentially expressed genes (DEGs). Upregulated DEGs are shown in red (A), while downregulated DEGs are shown in blue (B). The advanced bubble chart shows the top 20 enriched signaling pathways for DEGs. The X-axis represents the rich factor (rich factor = number of DEGs enriched in the pathway/total number of genes in the background gene set), and the Y-axis represents the enriched pathway. The color indicates the significance of the enrichment, and the size of the bubble corresponds to the number of DEGs enriched in the pathway; C-E: Expression of KEGG-enriched MAPK genes in colon and rectal cancer tissues. ^bP < 0.01; ^cP < 0.001; READ: Rectal adenocarcinoma; COAD: Colon adenocarcinoma.

Expression of Nav1.6-related MAPKs in different types of cells in CRC tissue

We analyzed CRC single-cell sequencing data to assess the expression and distribution of SCN8A and SCN8A-related functional genes across different cell types during CRC progression and predict the potential molecular mechanisms underlying SCN8A involvement in CRC. The analysis was performed using the GEO database GSE201348, which includes data from 48 polyps, 27 normal tissues, and six CRC samples. A total of 12 distinct cell lineages were identified in the non-paired controls (Figure 4A-D), with each corresponding to specific cell markers (Figure 4E). A bar graph was used to depict the proportion of cell types in healthy and CRC tissues. Compared to the healthy group, the proportion of functional cell types, such as pericytes, goblet cells, and glial cells, was significantly reduced in CRC, whereas the proportion of epithelial cells and fibroblasts was significantly increased. This suggests that as CRC progresses, the proliferation of epithelial cells in the tumor becomes more prominent (Figure 4F).

Figure 4 Expression distribution of top differentially expressed genes enriched in different cell types by KEGG analysis. A-D: UMAP visualization of single cells in normal colon tissue and tumor tissue. Each dot represents a cell, color-coded by cell type; E: Dot plot showing the top three marker genes for each identified cell type. The size of each dot indicates the proportion of the cell population expressing each gene, while the color reflects the expression level; F: Bar graph representing the average proportion of different cell types in normal and tumor tissues; G: UMAP distribution of Nav1.6 expression in different cell types in normal colon tissue and tumor tissue; H: Bubble plot showing the relative expression of differentially expressed genes enriched according to weighted correlation network analysis. The color and size of the dot represent the average expression and the percentage of positively expressed hepatocytes in each cluster, respectively; I-M: Expression of Nav1.6-related MAPK genes in epithelial cells in different stages of colorectal cancer; N: Tumor tissue microarray analysis of expression of MAP3K10 in tumor tissues and paratumor tissues; O: Immunofluorescence analysis of expression of MAP3K10 and Nav1.6 in tumor tissues and paratumor tissues. CRC: Colorectal cancer.

We further examined the expression of SCN8A and functional, T-stage-related genes across different cell types. Figure 4G shows that SCN8A was predominantly expressed in goblet and epithelial cells in healthy colon tissue. However, as CRC progressed, goblet cells gradually disappeared, and the proportion of epithelial cells increased, leading to SCN8A being predominantly expressed in epithelial cells in CRC. In contrast, SCN8A-related MAPK genes were widely expressed across various cell types, with genes such as MAP2K1 and MAP3K8 showing higher average expression levels and a greater proportion of positive monocytes (Figure 4H).

Next, we analyzed the difference in MAPK gene expression most strongly correlated with SCN8A expression in epithelial cells. As shown in Figure 4I-M, SCN8A and MAP3K10 were significantly upregulated as CRC progressed, whereas MAP2K1 and ANGPT2 were downregulated. The decreased expression of these genes, which are positively co-expressed with SCN8A at the tissue level, suggests that changes in gene expression within CRC tissues may arise from alterations in the expression and distribution of these genes across various cell types (Figure 4I-M). Among the selected genes, only MAP3K10 showed a trend toward upregulated expression during CRC progression, similar to that of Nav1.6. Therefore, we further investigated the expression of MAP3K10 in CRC using TMA analysis. The results indicated that MAP3K10 expression was significantly higher in CRC tissues compared to that in adjacent non-tumor tissues (Figure 4N). In addition, immunofluorescence showed that MAP3K10 had significant co-localization with Nav1.6 in tumor which suggested that MAP3K10 could interact with Nav1.6 in CRC tissues (Figure 4O).

Knockdown of Nav1.6 downregulates the MAPK signaling pathway in CRC cells

We established an SCN8A knockdown model using a colon cancer cell line to validate the correlation between SCN8A expression and MAPK signaling-related genes identified in TCGA and GEO databases and elucidate the molecular role of Nav1.6 in CRC metastasis. Initially, Nav1.6 expression was assessed in CaCo2, DLD-1, and HT-29 CRC cell lines. Western blot analysis revealed that Nav1.6 expression was highest in DLD-1 cells (Figure 5A). However, the scratch assay showed that the expression of Nav1.6 in different CRC cell lines did not significantly correlate with the wound repair ability of tumor cells (Figure 5B and C). Based on these results, DLD-1 cells were selected for the Nav1.6 knockdown experiments, where siRNA-mediated interference reduced Nav1.6 expression by approximately 50% (Figure 5D and E). Fluorescence photos of cells after transfection showed that the transfection efficiency of SiNav1.6 was around 50% (Figure 5F). Subsequently, the Transwell migration assays demonstrated a significant decrease in cell migration 24 h after SCN8A knockdown in DLD-1 cells (Figure 5G and H). CCK-8 assays were performed to assess the effect of Nav1.6 on cell viability, which showed that the viability of DLD-1 cells significantly decreased 24 and 48 hours post-SCN8A knockdown (Figure 5I). These findings suggest that Nav1.6 may be involved in the proliferation and metastasis of CRC cells.

Figure 5 Knockdown of Nav1.6 inhibits the invasion and proliferation of colorectal cancer cells. A: Expression of Nav1.6 in three different human colorectal cancer cell lines; B and C: Wound healing assays of cell migration in CaCo2, HT-29, and DLD-1 cells. Images of wound closure are shown at the indicated time points (0, 12, and 24 hours) after scratching; D and E: Expression of Nav1.6 in DLD-1 cells before and after siRNA-mediated interference; F: Images of DLD-1 cells after transfection with Control siRNA or siNav1.6; G: Transwell invasion assay before and after Nav1.6 knockdown in DLD-1 cells; H: Statistical analysis of the Transwell invasion assay results; I: Cell viability of DLD1 cells at different time points after transfection with SiNav1.6. ^bP < 0.01; ^cP < 0.001.

We performed RNA-sequencing (RNA-seq) analysis on DLD-1 cells before and after SCN8A knockdown to further investigate whether Nav1.6 contributes to the proliferation and invasion of CRC cells by regulating the MAPK signaling pathway. As shown in the volcano plot (Figure 6A), SCN8A knockdown resulted in the downregulation of 77 genes and the upregulation of 71 genes (|log2 fold change| > 1, P < 0.05). Among the top 20 DEGs, MAP3K10 and RAF1 were significantly downregulated after SCN8A knockdown, which corroborates the findings from TCGA data analysis (Figure 6B). This suggests that the MAPK signaling pathway may play a role in regulating CRC cell proliferation and invasion through Nav1.6. Given that MAP3K10 is associated with both the ERK signaling and TGF-β pathways, we analyzed alterations in the MAPK pathway before and after Nav1.6 knockdown. Western blot analysis revealed a marked decrease in the expression of MAP3K10, as well as in the ratio of p-ERK1/2 to total ERK1/2 and p-c-Raf to total c-Raf, in DLD-1 cells following Nav1.6 knockdown compared to the control group (Figure 6C-F). Additionally, the levels of MEK1/2 were reduced in DLD-1 cells after siNav1.6 treatment (Figure 6G). These findings suggest that Nav1.6 plays a role in CRC metastasis by regulating the MAPK signaling pathway.

Figure 6 Knockdown of Nav1.6 inhibits the RAF-MAPK signaling pathway. A: Volcano plot of differentially expressed genes (DEGs) before and after Nav1.6 knockdown; B: Heatmap of the top DEGs before and after Nav1.6 knockdown. Red indicates upregulated DEGs, and blue indicates downregulated DEGs; C: Expression of MAPK-related proteins MAP3K10, c-RAF, p-cRAF, ERK, p-ERK, and MEK1/2 before and after Nav1.6 knockdown in DLD-1 cells; D-G: Ratios of MAP3K10 to β-actin, p-ERK1/2 to total ERK1/2, p-cRAF to total c-RAF, and MEK1/2 to GAPDH in DLD-1 cells before and after Nav1.6 knockdown. ^aP < 0.05; ^cP < 0.001; ^dP < 0.0001; NS: Not significant.

DISCUSSION

CRC is the leading cause of cancer-related morbidity and mortality worldwide. Understanding the molecular mechanisms that drive CRC progression is essential for developing more effective therapeutic strategies[19]. In this study, we focused on the VGSC Nav1.6 (SCN8A), which has been implicated in various cancers, including CRC. Our results provide strong evidence that Nav1.6 is overexpressed in CRC tissues, particularly in the cytoplasm, and is positively correlated with tumor progression, particularly the T stage. Additionally, our functional studies highlight a potential role for Nav1.6 in promoting CRC cell migration and survival, suggesting that it could be a valuable therapeutic target.

We first confirmed the elevated expression of SCN8A in CRC through transcriptomic data analysis of TCGA data. Both READ and COAD samples exhibited significantly higher levels of SCN8A expression in tumor tissues compared to that in normal tissues, consistent with our previous findings in CRC and metastatic lymph nodes. We observed that SCN8A expression increased with advancing CRC T stages, a result that was further validated by TMA analysis in a cohort of 90 CRC patients. The positive correlation between SCN8A expression and local T stage strongly suggests that Nav1.6 is involved in tumor progression, particularly in the early stages of metastasis, a key determinant of CRC prognosis.

The aberrant subcellular localization of Nav1.6, with significant elevation in both the cytoplasm and nucleus of CRC tissues, adds a novel aspect to its role in CRC. The cytoplasmic expression of Nav1.6 was correlated with local T staging, providing further evidence that it may contribute to CRC progression[16]. The shift in localization from normal colon tissue (where Nav1.6 is primarily expressed in goblet and epithelial cells) to CRC tissue (where it is predominantly expressed in epithelial cells) suggests that Nav1.6 may play an evolving role in tumorigenesis as CRC progresses.

We conducted WGCNA and identified several SCN8A-associated DEGs in CRC tissues to further explore the molecular pathways involved. These genes were clustered into distinct modules, with one module positively correlated with advancing T stages and another negatively correlated. KEGG pathway analysis of these DEGs revealed that the upregulated genes were primarily involved in MAPK signaling, a well-known pathway implicated in cancer progression. Specifically, MAPK-related genes, such as MAP3K10, MAP2K1, and MAP3K8, were found to be overexpressed in CRC tissues and significantly correlated with SCN8A. Given the central role of MAPK signaling in cell proliferation, migration, and survival, our findings suggest that SCN8A may regulate CRC progression through its interaction with the MAPK pathway[20].

Single-cell RNA-seq data further supported this hypothesis. In healthy colon tissue, SCN8A was mainly expressed in goblet and epithelial cells. However, as CRC progressed, the proportion of goblet cells decreased, and SCN8A expression became more localized in epithelial cells, reinforcing the idea that Nav1.6 is primarily involved in epithelial cell proliferation during CRC progression. Moreover, MAPK-related genes such as MAP3K10 were expressed at higher levels in epithelial cells, highlighting the potential link between SCN8A and MAPK signaling in CRC.

Functional experiments in CRC cell lines, specifically DLD-1 cells, confirmed the role of Nav1.6 in regulating CRC cell migration and survival. SCN8A knockdown resulted in a significant reduction in cell migration and viability, which was accompanied by the downregulation of MAPK signaling components, including MAP3K10, c-RAF, and ERK1/2. These findings suggest that Nav1.6 may facilitate CRC progression by modulating the MAPK pathway, with potential implications for therapeutic targeting.

While these findings strongly suggest a role for Nav1.6 in CRC metastasis, some limitations of our study should be acknowledged. First, the transfection efficiency in our knockdown experiments was not optimal, with around 50% knockdown efficiency. This could have impacted the strength of the observed phenotypic effects, and further optimization of transfection methods may yield stronger results. Additionally, while our data suggest that Nav1.6 plays a role in promoting CRC cell migration and survival, further in vivo studies are needed to confirm the therapeutic potential of targeting Nav1.6 in CRC.

CONCLUSION

Our study provides compelling evidence that Nav1.6 is overexpressed in CRC tissues and is positively correlated with tumor progression. We propose that Nav1.6 contributes to CRC metastasis, at least in part, by regulating the MAPK signaling pathway. Given its critical role in tumor progression, Nav1.6 could serve as a potential biomarker and therapeutic target for CRC. However, additional studies are necessary to fully elucidate its functional role and therapeutic potential.

ACKNOWLEDGEMENTS

We express our deepest gratitude to all of the doctors and nurses at the Department of Gastrointestinal Surgery at Quzhou People's Hospital for their kind support.