Gene Expression Analysis of the Concomitant Existence of Lymphovascular and Perineural Invasion in Colorectal Cancer Nadiah Abu*, Nurul-Syakima Ab Mutalib, Rahman Jamal UKM Medical Molecular Biology Institute (UMBI), UKM Medical Center, Universiti Kebangsaan Malaysia, Jalan Yaacob Latiff, Cheras, 56000 Kuala Lumpur Received on 01/12/2017 / Accepted on 20/07/2018 ABSTRACT The invasion of cancer cells into the peritumoral, lymph node and perineural system could be detrimental on cancer patients. In colorectal cancer (CRC) patients, the presence of lymphovascular (LVI) and/or perineural (PNI) invasion could significantly influence on the survival rates, treatment options and recurrence tendencies. To date, no study has analyzed the molecular profile of the concomitant existence of LVI and PNI in CRC. Here, we reanalyzed The Cancer Genome Atlas (TCGA) CRC datasets and focused on cases where the information regarding LVI and PNI are available (n=176). We performed differential gene expression, methylation and microRNA analysis by comparing the groups having both or either LVI and PNI with the control group (LVI negative and PNI negative). Although there was no significant difference in the methylation and miRNA profiles, we identified a number of differentially expressed genes (DEGs). The comparison between the LVI+PNI+ and LVI-PNI- groups revealed key DEGs including SFTA2, PHACTR3, CRABP2, ODZ3, GRP, HAP1, CSDC2, TMEM59L and HDAC9. Meanwhile, in the LVI-PNI+ vs LVI-PNI- group, some of the DEGs found were PTPRR, EFNA2, FGF20, IGFL4, METRN and IGFBPL1. We believe that this study could be beneficial and add value to further understand the complex molecular profiles of CRC. Introduction Colorectal cancer (CRC) is the third most commonly diagnosed cancer in men and women worldwide (1). CRC has multiple clinical pathologies that allow for effective staging and classification when considering treatment. For instance, among the clinical signs assessed include the tumor volume, presence of lymph node metastasis, blood vessels, lymphovascular invasion and perineural invasion. Lymphovascular invasion (LVI) is one of the most critical steps in metastasis and the process involves the invasion of tumor cells in the lymphatic spaces, blood vessels or both in the peritumoral area (2). Moreover, LVI has been linked to other cellular pathways such as angiogenesis, cell proliferation and cell migration (3). LVI is associated with a wide variety of cancers including breast cancer (3), endometrial cancer (4), prostate cancer (5), gastric cancer (6) and CRC (7, 8). The usage of LVI as a prognostic marker for survival in a number of different cancers, especially CRC, has been widely studied (9-11). In fact, the presence of LVI has been a factor for clinicians to consider before making an informed decision regarding treatment (12, 13). For instance, in invasive breast carcinoma, the presence of LVI can be used to identify which adjuvant treatment is suitable (12, 14). Similarly in papillary thyroid carcinoma, patients with LVI are strongly recommended to receive post-operative radioactive iodine (RAI) therapy to have a better outcome (13). Besides LVI, perineural invasion (PNI) is also an important factor in determining treatment and survival rates (9). PNI is a process whereby cancer cells migrate and invade the surrounding nerves (15). PNI was first discovered in head and neck cancer, however, to date, multiple other types of cancer cells have been reported to be able to invade the perineural site as well (15, 16). For instance, in prostate cancer, it has been well established that PNI contributes to cancer aggressiveness and subsequently, lower survival rates (17-19). In CRC, the presence of PNI is associated with a wide range pathological features such as poor tumor differentiation, high grade tumor classification and aggressive tumor behavior (7, 20, 21). PNI assessment has not been a major aspect in the stratification of CRC patients. However, recently, it has been reported that PNI deserves extra attention especially as a prognostic factor (22). In terms of therapy in CRC patients, PNI should also be given considerable attention in influencing the treatment given (23). Adjuvant therapy for CRC patients are usually given based on the presence of clinical signs such as LVI, PNI or presence of secondary tumors (24, 25). There have been studies conducted that assessed the efficacy of a treatment based on different types of clinical pathology (25). Nevertheless, in this post-genomic era, besides assessing the clinical signs, we also need to correlate it with the molecular profile of each tumor. Since CRC, as with any other cancer, is known to be molecularly and clinically heterogeneous, combining both information would allow based on Illumina Methylation 450 beadchip in β value format was retrieved from cBioPortal (http://www.cbioportal.org/data_sets.jsp) while level 3 RNA sequencing dataset in RSEM (RNA-Seq by Expectation Maximization) format and microRNA sequencing in normalized reads per million (RPM) format were retrieved from Firebrowse (doi:10.7908/C1F76BX1). TCGA CRC patients’ information were obtained from GDC Data Portal (https://portal.gdc.cancer.gov/). From there, we filtered for the presence of lymphovascular invasion (LVI) and perineural invasion (PNI) status, where we excluded 280 datasets due to the lack of information regarding LVI and PNI invasion status. oncologists to execute a more comprehensive treatment approach. However, the concomitant existence of both LVI and PNI in CRC may make it more difficult to treat and compare. For instance, in a Swedish cohort of rectal cancer patients, it was discovered that the rate of recurrence was higher in patients having both LVI and PNI instead of just having either LVI or PNI (26). Similarly, it was reported that the overall survival rate was much lower in gastric cancer patients having both LVI and PNI (27). Therefore, in an effort to determine the concomitant existence of LVI and PNI in CRC, we re-analyzed The Cancer Genome Atlas (TCGA) data based on the molecular profile of LVI and/or PNI. We performed analysis related to gene expression, methylation expression by comparing CRC patients with LVI and/or PNI against patients without LVI and PNI. We identified differentially expressed genes that could be used to further understand the molecular profiles of tumors with invasive behavior. Materials and methods Data acquisition A total of 456 datasets corresponding to colon adenocarcinoma (COAD) cases (cite Cancer Genome Atlas Research Network, 2012) were downloaded from the TCGA portal in April 2017. Methylation data The final number of datasets used for the analysis was 176, where we divided it into 4 groups. The groups are LVI+PNI+ (presence of both LVI and PNI), LVI- PNI+ (absence of LVI and presence of PNI), LVI+PNI- (presence of LVI and absence of PNI) and LVI-PNI- (absence of both LVI and PNI). The patients’ characteristics and clinicopathological information are displayed in Table 1. Survival analyses Kaplan-Meier survival analyses were performed for the overall survival analysis (OS) and disease-free survival analysis (DFS) for patients with follow-up details. The OS analysis was conducted based on the duration from date of diagnosis to death, while for the DFS analysis, the duration was selected from the date of diagnosis to the event of relapse or recurrence (28). All the analyses were conducted on Graph Pad Prism software version 7.00 (Graphpad Prism, USA) using the Log-rank Mantel Cox Test. Bioinformatic analyses The miRNAseq and RNAseq level 3 datasets and methylation dataset from TCGA were used exclusively. For the RNAseq analysis, we used the RNAseq by Expectation-Maximization (RSEM) values to quantify the expression level of mRNAs. Additionally, for the miRNAseq analysis, we used the normalized reads per million (RPM) for quantification. Afterwards, we performed the Students unpaired t-test using the Benjamini Hochberg (HB) false discovery rate (FDR) multiple testing correction using the Bioconducter in R version 3.2.1. We then calculated the log2 fold change based on the mean expression of normalized log2 values for each category. Statistical significance was set at P≤0.05.Heatmaps were generated using Morpheus from the (https://software.broadinstitute.org/morpheus/). Venn diagrams were created using an online maker tool at (http://bioinformatics.psb.ugent.be/webtools/Venn/). The overall study design is displayed in Figure 1. Pathway enrichment analyses and gene ontology analyses Pathway enrichment and gene ontology analysis were performed using the online tool WEB-based GEne SeT AnaLysis Toolkit v. 2017 (WebGestalt) (www.webgestalt.org)(29). For the LVI+PNI+ comparison group, the input list contained 72 differentially expressed genes with P value ≤0.05, from which 69 genes were unambiguously mapped to the unique Entrez Gene IDs and 3 genes were mapped to multiple Entrez Gene IDs or could not be mapped to any Entrez Gene ID. For the LVI-PNI+ comparison group, the input list contained 43 differentially expressed genes with P value≤0.05, form which 39 genes were mapped to the unique Entrez Gene IDs and 4 genes were mapped to multiple Entrez Gene ID or could not be mapped to any Entrez Gene ID. All mapped Entrez Gene IDs were retrieved from the platform genome_protein-coding. The minimum number of Entrez Gene IDs in the category was set at 5. The FDR method used was BH. Moreover, we also performed the analyses using the Database for Annotation, Visualization and Integrated Discovery (DAVID) as a comparison. The pathway enrichment analysis was performed based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (30), and the gene ontology analysis was performed based on the biological process and molecular function. Analyses with Benjamini- adjusted P value ≤ 0.05 were considered statistically significant. Results Clinicopathological characteristics of CRC patients with LVI and/or PNI Overall, a total of 176 CRC patients were used for this analysis, where 24 were in the LVI+PNI+ group, 22 for the LVI-PNI+ group, 28 for the LVI+PNI- group and 102 for the LVI-PNI-group. Based on the demographic and clinicopathological features displayed in Table 1, there were considerable differences in the disease-free status, overall survival status, tumor stage and lymph node metastasis stage. The effects of LVI and PNI on the overall survival and disease-free survival analyses in CRC patients As illustrated in Table 1, among the four groups, LVI- PNI+ group had the lowest mean overall survival (OS) with 22.11 months, whereas the highest mean OS was LVI-PNI- group with 29.69 months. However, even though LVI-PNI+ patients had the lowest mean OS, the Kaplan-Meier survival curve was not significant against patients with LVI-PNI- (Log rank Mantel Cox Test, P=0.1792, Figure 2A). Meanwhile, for patients with LVI-PNI+, there was no statistically significant difference in terms of OS against patients with LVI- PNI- (Log rank Mantel-Cox Test, p=0.728, Figure 2B). Based on our Kaplan-Meier survival curve analysis as shown in Figure 2, patients who had both LVI and PNI (LVI+PNI+) had a lower OS rate as compared to patients who did not have LVI and PNI (LVI-PNI-) (Log rank Mantel-Cox Test, P=0.004, Figure 2C). Similarly, for the mean disease-free survival, LVI-PNI+ group had the lowest value with 17.89 months. However, interestingly, the LVI+PNI- group had the highest mean disease-free survival with 29.33 months compared to the LVI-PNI- group. As displayed in Figure 2F, the DFS for LVI+PNI+ against LVI-PNI- was not statistically significant (Log rank Mantel Cox Test, P=0.138). Likewise, when comparing LVI+PNI- against LVI-PNI-, no statistical significance was observed (Log rank Mantel Cox Test, P=0.104, Figure 2D). Interestingly, for the comparison of LVI-PNI+ against LVI-PNI-, the P value obtained was 0.0015, making the DFS curves significantly different (Figure 2E). Differentially expressed genes (DEGs) Based on the RNAseq reanalysis, only the LVI+PNI+/LVI-PNI- group and LVI-PNI+/LVI-PNI- group have differentially expressed genes (DEGs) with P values ≤ 0.05. According to Table 2, the LVI+PNI+/LVI-PNI- group had 76 differentially expressed genes, but 4 had to be eliminated from the analysis due to the low level of expression. The genes with log2 fold change more than 1.4 include SFTA2, PHACTR3, CRABP2, ODZ3, GRP, HAP1 CSDC2, HDAC9 and TMEM59L. The comparison of LVI- PNI+ against LVI-PNI- group produced 50 DEGs, but 8 had to be removed due to the low expression as shown in Table 3. Genes that had fold changes of more than log2 1.2 were EFNA2, PTPRR, METRN (up-regulated genes) and IGFL4, IGFBPL1 and FGF20 (down-regulated genes). The heatmap visualization for both sets of genes is displayed in Figure 3A and 3B. For the comparison of LVI+PNI- against LVI-PNI-, no significant DEGs were observed, so we reduced the stringency, and set the P value at p≤0.1. From there, only 32 DEGs were identified, where 5 genes were removed due to the low level of expression. Interestingly, as shown in Figure 3C, based on the two sets of genes, there were no overlapping genes observed in all two sets. Only 9 genes were found to be overlapping between the LVI- PNI+ set and LVI+PNI-set. The overlapped genes are ARL2, C11orf68, C11orf83, C17orf67, EIF2AK3, POMZP3, SLC25A39, SYK and THAP7. Gene Ontology and Pathway Enrichment Analysis The gene ontology analyses were performed using the DEGs discovered in the two sets of genes as shown in Table 4. We used both the WebGestalt and the DAVID program as a comparison in terms of the statistical significance. For the first set of comparison group (LVI+PNI+/LVI-PNI-), only 1 pathway was found to be enriched, which is the pathway related to microRNAs in cancer. However, only the WebGestalt application managed to identify this pathway. Next, we performed gene ontology analysis, for the biological process, where angiogenesis and regulation of calcium ion were significantly enriched. For the molecular function, most of the genes were involved in transporter-related activity. In the second set of genes (LVI-PNI+/LVI- PNI-), the enriched pathways were MAPK and PI3K- AKT pathways. For the gene ontology analysis, the enriched biological process were transmembrane receptor tyrosine kinase pathway and ammonium transport. Meanwhile for the molecular function, the phosphatase binding function was enriched. Methylation and miRNA analyses We also reanalyzed the global methylation and miRNA differential expression analysis for each of the comparison groups. However, no significant difference (P≤0.05) in the methylation status were found. Similarly, for the miRNA expression, there were no differentially expressed miRNAs (P≤0.05) discovered. Discussion In this study, we analyzed the gene expression data of the concomitant existence of lymphovascular invasion (LVI) and/or perineural invasion (PNI) in colorectal cancer (CRC) patients using the Cancer Genome Atlas (TCGA) dataset. For the overall survival analysis (OS), only the LVI+PNI+ group showed a statistically significant difference in terms of survival than the LVI-PNI- group. This indicates that CRC patients with both LVI and PNI have a shorter survival rate than patients with either LVI only or PNI only. Interestingly, even though the OS for the LVI+PNI+ CRC patients was more significantly affected, the patients with LVI- PNI+ showed a higher tendency of recurring/progressing according to the DFS analysis. In gastric cancer patients, it was reported that only LVI+PNI+ patients had a significant difference in DFS and OS rate than the other groups (27). This suggests that in different types of cancer, both LVI and PNI play ultimately different roles with regards to the OS and DFS. From the gene expression analysis, a total of 72 differentially expressed genes (DEGs) were discovered in the LVI+PNI+ vs LVI-PNI- comparison. Of all the deregulated genes, SFTA2 had the highest difference in terms of fold change. SFTA2 or also known as “surfactant associated 2” has been reported to be overexpressed in lung cancer (31). SFTA2 was predicted to interact with other mucin-like proteins that are highly abundant in the lung (31). It was only recently that SFTA2 was proposed as a biomarker to differentiate between lung adenocarcinoma and squamous cell carcinoma (32). However, not much research has been conducted to evaluate the molecular function of SFTA2 in regard to carcinogenesis especially in relation to LVI and PNI. Another gene that was upregulated in the LVI+PNI+ vs LVI-PNI- group is the phosphatase and actin regulator 3 (PHACTR3) gene. This gene is part of the PHACTR family which is highly involved in the reorganization of the cytoskeleton and actin (33). It was reported that PHACTR3 enhances cell motility and adhesion of HeLa cells by interacting with the actin cytoskeleton (34). This is in concordance with our hypothesis, since LVI+PNI+ patients have a higher rate of metastasis, proteins related to the movement of cancer cells should be elevated. A study conducted by Bosch et al found that the mRNA expression of PHACTR3 gene was low in CRC samples and significantly hypermethylated in advanced CRC than in normal colonic mucosa (35). In our study however, the expression of PHACTR3 was higher in advanced CRC (LVI+PNI+), and there was no significant difference in terms of the methylation status. This reflects the molecular heterogeneity in CRC and could explain why these findings contradict each other. Another interesting finding from the LVI+PNI+ vs LVI-PNI- comparison is the upregulation of CRABP2. CRABP2, known as Cellular Retinoic Acid Binding Protein 2, is a member of the retinoic acid and cytosolic fatty-acid binding protein family (36, 37). This gene functions by transporting vitamin A to retinoic acid receptors and shuffling the complexes between the cytoplasm and nucleus (36-39). In cancer however, there have been multiple contradicting reports as to whether CRABP2 is involved in the pro- tumorigenic or anti-tumorigenic process (36, 38, 39). In pancreatic cancer, a study conducted by Xiao et al found that CRABP2 could be used as a molecular marker for distinguishing the high grade pancreatic ductal adenocarcinoma cases from benign pancreatic conditions (38). This shows that CRABP2 is highly expressed in advanced pancreatic cancer, similar to what we have discovered in CRC. Additionally, in non-small lung cancer cells (NSCLC) patients, the expression of CRABP2 was also found to be upregulated (40). Favorskaya et al reported that the mRNA expression of CRABP2 was negatively correlated with the presence of lymph node metastasis in NSCLC patients (40). Interestingly, in our report, the LVI+PNI+ group had the highest presence of lymph node metastasis, and also a higher expression of CRABP2. This shows that the role of CRABP2 in different tumors is still not fully understood yet. Moreover, the expression of CRABP2 was shown to be significantly down regulated in esophageal cancer, confirming the unknown molecular dynamics of CRABP2 in different types of cancer (39). Besides CRABP2, the ODZ3 gene was also upregulated in LVI+PNI+ group. Little is known about ODZ3 except that it is part of the teneurins family, a group of highly conserved transmembrane proteins involved in intracellular signaling related to neural circuits and development (41, 42). In fact, ODZ3 was found to be heavily involved in microphthalmia in humans (43). In cancer, it was reported that ODZ3 was significantly mutated in colorectal cancer tissue (44-46). Apart from that, the molecular function of ODZ3 in relation to LVI and PNI is still elusive. The gastrin-releasing peptide (GRP) was also discovered to be overexpressed in the LVI+PNI+ group. This gene is known to be involved in multiple types of cancer including breast cancer (47), esophageal squamous cell carcinoma (48) and colon cancer (49). GRP acts as a growth factor by binding to GRP-receptors (GRPR) and is involved in other types of growth factor receptor signaling pathways such as EGFR and VEGFR (50). Molecules involved in growth factor signaling are highly associated with metastasis, angiogenesis and invasion. In fact, it was suggested that GRP and GRPR are involved in neoangiogenesis and microvascular perfusion in renal cell carcinoma (51). Moreover, a study conducted by Fleishcmann et al showed that GRPRs can be promising targets for vascular targeting applications since they are highly expressed in the vascular bed of various cancers (52). In CRC, a correlation was found between the high levels of mRNA of GRPR and lymphatic vessel invasion (45). This is similar to what we have discovered, where patients with more invasive CRC have a higher expression of GRP. The CSDC2 gene was also upregulated in regards to the LVI and PNI status in CRC. Nevertheless, there have been no reports on the function of CSDC2 in relation to CRC or any other types of cancer. Another gene that was upregulated was the TMEM59L gene. This gene is known to be involved in neuron-related pathways that modulates cellular oxidative stress (53). In terms of cancer, similar to the CSDC2 gene, there were limited information regarding its role in cancer. Since we did not find any significant difference in the methylation and miRNA profiles among all the groups, we hypothesized that there should be another epigenetic mechanism regulating the gene expression, such as long noncoding RNA or the recently discovered circular RNA. In the LVI+PNI+ vs LVI- PNI- gene set, we found that HDAC9, or histone deacetylase 9 gene was significantly upregulated. HDAC9 is part of the histone deacetylases that are commonly deregulated in tumorigenesis and are epigenetically involved in multiple gene regulation. In fact, HDAC9 has been reported to be overexpressed in different types of tumor including breast cancer (54) and oral squamous cell carcinoma (55). HDAC9 was found to be targeting several genes including SOX9 (54) and TRIM29 (56). In CRC, the inhibition of HDACs was shown to decrease the expression of EGFR. This could be particularly beneficial for patients who have the mutant KRAS where cetuximab is rendered inefficient (57). Thus knowing that HDAC9 is overexpressed in LVI+PNI+ patients could provide a more informed decision in treating these patients. We also found 42 DEGs when comparing the LVI- PNI+ group against the LVI-PNI- group. In this set of genes, it was discovered that the PTPRR gene is significantly upregulated. This gene codes for the protein tyrosine phosphatase receptor-type R enzyme and is part of the protein tyrosine phosphatase (PTP) family (58). The PTPs are often involved in multiple signaling pathways including MAPK, ERK, cell growth and cell cycle-related pathways. PTPRR has been primarily viewed as a neuronal phosphatase which regulates the ERK pathway and it has been associated with major depressive disorder (MDD) (59). In cancer, interestingly, the levels of PTPRR is inversely correlated with the tumor grade. For instance, in oral squamous cell carcinoma (OSCC), the expression of PTPRR is higher in lower grades of OSCC and vice versa (60). Similarly, in CRC, the expression of PTPRR is decreased in CRC tissue and cell lines as compared to normal samples (61). In fact, Menigetti et al also reported that PTPRR is epigenetically silenced in the early events of CRC carcinogenesis. However, in our findings, we found that patients who were PNI positive, had a high expression of PTPRR than patients who were negative for PNI. Since, PNI is a process involving the neural network, it could be postulated that PTPRR could be reactivated upon PNI. However, further studies need to be conducted to conclusively understand the dual roles of PTPRR in CRC. Besides PTPRR, the EFNA2 gene was also increased in the LVI-PNI+ group. EFNA2 is the gene encoding the Ephrin-A2 protein, a member of the ephrin family. The ephrin family can be divided into ephrin A and ephrin B depending on the type of ligands they bind to. In prostate cancer, it has been shown that ephrin A modulates the contact inhibition of locomotion (62). It has also been reported that in a Japanese patient diagnosed with adenoid cystic carcinoma, the expression of ephrin A2 was high and is associated with the clinical feature of perineural invasion (63). This was in agreement with our findings where PNI in CRC is correlated with the ephrin A2 expression. Another gene that was upregulated in this group comparison is the METRN gene. This gene codes for the meteorin protein, ans has been shown to be involved in the neurogenesis of glial cells (64). Apart from that, no other studies have reported on the function of this gene, much less in cancer. For the downregulated genes in LVI-PNI+ group, the three most significant ones are FGF20, IGFBPL1 and IGFL4. FGF20, also known as fibroblast growth factor 20, is part of the fibroblast growth factor (FGF) family (65). In normal human tissues, FGF20 is highly expressed in the brain, especially in the cerebellum (65). FGF20 was also found to be overexpressed in SW480 colon cancer cell line, LXI liver cancer cell line and NCI-N87 gastric cancer cell line (65). However, there are few studies which have been performed to understand the role of FGF20, especially in terms of its relationship with perineural invasion. Furthermore, IGFBPL1, or insulin-like growth factor binding protein-related protein 1, is a member of the IGFBP superfamily that affects the expression of IGF- 1 and IGF-2 (66). This protein has been reported to be significantly present in the developing forebrain of mice (67). IGFBPL-1 was found to be down-regulated in breast cancer (66) and that the methylation of IGFBPL-1 was associated with a lower OS and DFS in breast cancer patients. In our report, the expression of IGFBPL1 was indeed lower in the LVI-PNI+ group, but there were no significant methylation patterns observed. Furthermore, the IGFL4 gene, which is part of the small IGFL family, was also down-regulated in the LVI-PNI+ group (68). The specific function of IGFL4 has not been reported elsewhere, but it has been shown to be expressed in the cerebellum (68). Interestingly, there were no overlapping genes found in between the 2 sets of genes. This suggests that cancer cells undergoing PNI only, have a different regulation than cells undergoing both LVI and PNI. This suggests that the regulation of cancer cells in LVI and PNI is more complex than just the elevated expression of the predicted invasion-related genes. Moreover, in between the LVI-PNI+ and LVI+PNI- group, there are 9 overlapping genes, but since the P value set for the LVI+PNI- group was P≤0.1, there could only be a weak association between the genes. Furthermore, we performed gene ontology and pathway enrichment analyses for the DEGs found in the LVI+PNI+ vs LVI-PNI- group. We did not find any significantly enriched pathways which suggest that the LVI and PNI processes are dependent of each other. There is no particular pathway that could lead to both the concomitant existence of LVI and PNI in CRC. Targeting specific molecules that are deregulated in the presence of both LVI and PNI could be more feasible than targeting specific pathways. Meanwhile, for the LVI-PNI+ vs LVI-PNI- comparison, the pathways that were enriched include the MAPK and PI3K-AKT pathway. This is in concordance with the genes that were significantly deregulated, especially PTPRR and FGF20, where their expressions were reported to affect these pathways. Conclusion In summary, the comparison of LVI+PNI+ vs LVI- PNI- showed that most of the DEGs found were related to cell metastasis and angiogenesis. Interestingly, there were some contrasting evidence found between the expression of these genes in our analysis and previously reported studies. This indicates that CRC is a molecularly heterogeneous disease and different cases of cancer have different gene expression profiles. Nevertheless, the information retrieved from this part of the study could have further implications for LVI+PNI+ patients. Some of the most deregulated DEGs (PTPRR, FGF20, IGFBPL1 and IGFL4) found in the LVI-PNI+ vs LVI-PNI- comparison were related to neural network, which is highly probable since perineural invasion takes place mostly around the neural-related areas. These findings could be of further use to the treatment of PNI positive CRC patients. Collectively, our results show that though there were no significant difference in the methylation and miRNA profiles, there are some DEGs that could be useful for the LVI+PNI+ and LVI-PNI+ patients. However, further in-depth study is needed to actually address the effects of these genes to the overall outcome of therapy and prognosis in CRC.