The Midkine-related ceRNA Network is a Prognostic Marker of Colon Cancer With Perineural Invasion


 BackgroundPerineural invasion (PNI) is a prominent characteristic of multiple solid tumors and indicates poor prognosis[1, 2]. A growing body of evidence emphasizes the critical role of competitive endogenous RNA(ceRNA) regulatory networks in various human cancers. However, the complexity and behavior characteristics of the ceRNA network in colon cancer with perineural invasion were still unclear. In this study, we aimed to clarify a midkine(MDK)-related ceRNA regulatory network and identify potential prognostic markers associated with colon cancer with perineural invasion.Material and MethodsWe extract the expression profiles of three RNAs (long non-coding RNAs [lncRNAs], microRNAs [miRNAs], and mRNAs) from The Cancer Genome Atlas (TCGA) database. In order to construct the lncRNA-miRNA-mRNA triple regulatory network in colon cancer with perineural invasion, we performed joint analysis in the MDKhigh and MDKlow expression groups, as well as the colon cancer and paracancerous tissue groups. Kaplan-Meier method was used to plot the Overall survival (OS) curves, and the log-rank test was used for testing. Univariate and multivariate Cox regression analysis were used to determine OS-related characteristics.ResultsKCNQ1OT1-miR-454-3p/miR-204-5p-MAP1B ceRNA network related to the prognosis of colon cancer was obtained through bioinformatics analysis. Importantly, we determined the KCNQ1OT1/MAP1B axis in the ceRNA through correlation analysis, and through Cox regression analysis it seems to be a clinical prognostic model. ConclusionIn conclusion, the current study constructing a ceRNA based KCNQ1OT1/MAP1B axis may be a novel important prognostic factor for the diagnosis and prognosis of colon cancer with perineural invasion.

obtained through bioinformatics analysis. Importantly, we determined the KCNQ1OT1/MAP1B axis in the ceRNA through correlation analysis, and through Cox regression analysis it seems to be a clinical prognostic model.

Conclusion
In conclusion, the current study constructing a ceRNA based KCNQ1OT1/MAP1B axis may be a novel important prognostic factor for the diagnosis and prognosis of colon cancer with perineural invasion.

Background
Perineural invasion (PNI), malignant invasion of nervous structures and nerve sheaths that could be a source of metastatic spread beyond the extend of any local invasion, indicating an increased risk for poor prognosis and decreased survival in colon cancer. [3][4][5] Then, Midkine(MDK) is a secretory alkaline heparin-binding protein and a neural growth factor. The MDK family is comprised of two members, MDK and pleiotrophin(PTN) [6]. MDK and PTN share a number of common receptors and biological characteristics, including interaction with heparin [7] and development of the nervous system [6]. Studies have showed that MDK is overexpressed in colon cancer and levels of MDK in tumor samples of patients with colon cancer are signi cantly higher compared to control samples. [8] Additionally, the expression of MDK is being closely related to PNI, while levels of MDK and syndecan-3 are signi cantly higher in patients with PNI than in those without PNI [9]. Thus, these observations provide evidence that MDK may promote tumorigenesis and peripheral nerve in ltration of colon cancer. Additional studies are needed to further understand the role of MDK in colon cancer development and whether it is a viable therapeutic target.
In 2011, miscellaneous researchers token the ceRNA fray, which can de ne the transcriptome and form a large-scale regulatory RNA network. They avowed a detailed post-transcriptional regultory network, in which Long non-coding RNAs (lncRNAs) and other RNAs base competing with microRNAs (miRNAs) and acting as natural miRNA sponges by virtue of sharing no less than one miRNA response element (MRE). Moreover,these lncRNAs and protein-coding mRNAs can compete with miRNAs and participate in the regulation of this complex network by combining with miRNAs through MRE [10,11]. Many researchers have found that ceRNA regulation can participate in the occurrence and development of colon cancer with perineural invasion [12][13][14]. Thus, these networks should attach importance into complex gene interactions and identify potential biomarkers to diagnose and treat colon cancer with perineural invasion. We can acquire the RNA sequencing data from TCGA or GEO database which provide lncRNA, miRNA and mRNA data of various cancers and act as an excellent resource for data mining and biological discovery [15]. Based on these public databases, comprehensive ceRNA regulatory networks were built to explore more accurate prognostic biomarkers in more studies.
In our study, a ceRNA network related to MDK in colon cancer with perineural invasion was established in Figure 1. In order to construct the colon cancer-related lncRNA-miRNA-mRNA triple regulatory networks ,we set two groups of MDK high and MDK low expression in 449 colon cancer samples by differential expression analysis. Next, Functional enrichment analysis is used to evaluate the functional role and potential mechanisms of the network in colon cancer with perineural invasion. Then, a key ceRNA network was identi ed by expression analysis, survival analysis, and nuclear-cytoplasmic localization analysis, and RNAs of the network came from hub triple regulatory networks. Moreover, the correlation of four predictive genes and MDK can be performed and the vital role of KCNQ1OT1/MAP1B axis can be found in colon cancer with perineural invasion. Finally, the diagnostic and prognostic function of MAP1B in colon cancer with perineural invasion was identi ed by Cox regression analysis, and the possible values of MAP1B were obtained by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses.

Data Collection and processing
RNAseq data, miRNAseq data, and corresponding COAD clinical data were downloaded from the TCGA data portal. Meanwhile, to further validate our results ,we also downloaded the gene expression pro les from GEO: GSE44076 (including 98 colon cancer samples and 98 adjacent nontumor samples; platform, GEO: GPL13667). The expression of MDK and MAP1B in colon cancer cell lines and at the protein level was veri ed by searching CCLE (https://portals.broadinstitute.org/ccle) and HPA (http://www.proteinatlas.org/) .The mutation status of MDK and MAP1B were obtained through the cBioPortal for Cancer Genomics.
DElncRNAs, DEmiRNAs, and DEmRNAs in COAD MDK high -MDK low in colon cancer samples and colon cancer-adjacent nontumorous samples were set to construct four groups. Differentially expressed lncRNAs (DElncRNAs), differentially expressed miRNAs (DEmiRNAs), and differentially expressed mRNAs (DEmRNAs) between COAD samples and corresponding paracancerous tissues were analyzed and normalized by using the edgeR package in R.
GraphPad Prism 9 software (version 9.0.2) were used to visualize Volcano plots of the DERNAs and the heatmap clustering.

Construction of the ceRNA Regulatory network
The miRcode database was used to match DElncRNAs and DEmiRNAs. MiRNA target genes were predicted on the basis of two databases: miRDB and TargetScan (version 7.2). Subsequently, we integrated the interaction between DEmiRNAs and DElncRNAs or DEmRNAs to construct a ceRNA regulatory network. Cytoscape (version 3.8.2) was used to visualize the ceRNA network.
We obtain DElncRNA sequences from LNCipedia, and identify the DElncRNA cellular localization based on its sequence by the lncLocator database. The hub triple regulatory network was identi ed by The Cytoscape plug-in cytoHubba.
Functional enrichment analysis "ClusterPro ler" package in R software was used for functional enrichment analysis, and GO biological processes and KEGG pathways at the signi cant level (q-value < 0.05) were employed.

Survival analysis and establishment of a speci c prognosis model
To determine the prognostic characteristics of DERNAs, combining the clinical data the survival curves of these samples with differentially expressed mRNA, lncRNA and miRNA were plotted by using GraphPad Prism 9 software (version 9.0.2) based on Kaplan-Meier curve analysis. P values < 0.05 were regarded as signi cant.

MDK overexpression act as the perineural invasion suppressor role and prognostic biomarker in Colon Cancer
We found that MDK was downexpressed in normal colon tissue, but overexpressed in colon cancer tissues based on TCGA database( Figure 2A).To further investigate the possible role of MDK in colon cancer, overexpression of MDK was demonstrated by immunohistochemistry (IHC) staining obtained from the Human Protein Atlas (HPA) .Compared with MDK in tumor tissues is mainly expressed in epithelial cells and neurons, normal tissues are mainly expressed in glandular cells .The patient data from the IHC are listed in( Figure 2B, Table S1).Then, MDK expression level in perineural invasion positive patients is higher than patients without perineural invasion based on TCGA clinical database.( Figure 2C) Owing to overexpression of MDK in tumor specimens with perineural invasion, we then studied the clinical signi cance of MDK expression in colon cancer patients. We found that poor overall survival (OS) in patients of the colon cancer cohort was associated with higher expression of MDK according to the Kaplan-Meier survival curves ( Figure 2D).
In a word, our data demonstrate that MDK expression is overregulated in colon cancer .The high expression of MDK in tumor tissues may be related to colon cancer perineural invasion and lead to a decrease in overall survival ability.

Identi cation of Signi cant Differential Genes
Based on the above analysis, a potential prognostic model for colon cancer with perineural invasion patients can be constructed by the ceRNA network related to MDK. Also, we set MDK high and MDK low expression groups according to the expression levels in colon cancer samples ,which are different from those in cancer and paracancerous groups. Construction of the lncRNA-miRNA-mRNA triple regulatory network We conducted a joint analysis in the MDK high and MDK low expression groups as well as in colon cancer and normal colon tissue groups .First, we put the remaining DElncRNAs into the mircode and TarBase database to identify potential miRNAs targeting lncRNAs. Six predicted miRNAs were selected from 105 DEmiRNAs. We then identi ed target genes of the common 6 miRNAs by selecting mRNAs shared by TargetScan and miRDB database. 708 mRNAs were predicted by all two databases. Finally, these candidate target mRNAs and above-mentioned 1255 DEmRNAs were further intersected to obtain 36 hub mRNAs. Figure 4A presents the relationship including 36 mRNAs ,17 lncRNAs and 6 miRNAs.
To further seek the potential roles associated with the triple regulatory network, functional enrichment analysis (including GO and KEGG) was used by Metascape. The results suggested that the DEmRNAs involved in the triple regulatory network were particularly enriched in the "GPCR ligand binding," "axoneme assembly," and "Neuroactive ligand-receptor interaction"( Figure 4C).

Validation of a model with perineural invasion-speci c prognostic value
To construct a accurate ceRNA of great prognostic value in colon cancer with perineural invasion, we rst analyzed the expression levels of RNAs from the hub triple regulatory network.
Next, the target sites in the KCNQ1OT1 and MAP1B 3'UTRs were forecasted to pair with miR-454-3p and miR-204-5p by TarBase and TargetScan ( Figure 7C). Moreover, correlation analysis was utilized to indicate a positive relationship between KCNQ1OT1 expression and MAP1B expression ( Figure 7D). Therefore, the KCNQ1OT1/MAP1B axis in the ceRNA network was selected as a potential prognostic model for the following step analysis.

Clinical relevance of the KCNQ1OT1/MAP1B axis in Colon Cancer with perineural invasion Patients
To explore whether the expression levels of KCNQ1OT1 and MAP1B were affected by clinical factors, we analyzed the association of two RNAs with clinical characteristics. The results demonstrated that both KCNQ1OT1 and MAP1B expression levels were positively correlated with the age(KCNQ1OT1, p = 0.035; MAP1B, p = 0.048) and TNM stage(KCNQ1OT1, p = 0.030; MAP1B, p = 0.017) as shown in TableS2. However, there is no remarkeble correlation in sex, lymph node metastasis, distant metastasis and BMI (p > 0.05; Table S2).
Moreover ,we use Univariate and multivariate Cox regression analyses to further analyze the prognostic signi cance of clinical features. In univariate Cox regression analysis models of KCNQ1OT1 and MAP1B, some prognostic factors (TNM stage, lymph node metastasis, and distant metastasis) were signi cantly associated with OS (p < 0.05; Table S3-S4) in colon cancer patients in TCGA cohorts. Importantly, both KCNQ1OT1 (hazard ratio [HR] = 1.665, p = 0.046) and MAP1B (HR = 1.249, p = 0.031) highexpression levels were closely relevant to a worse prognosis(TableS4-S4). In multivariate Cox regression analysis of KCNQ1OT1, age (HR = 0.669, p = 0.000), distant metastasis (HR = 1.529, p = 0.000), and Lymph-node metastasis (HR = 1.741, p = 0.012) were still relevant to OS in colon cancer patients. However, increased KCNQ1OT1 expression was not associated with a poor prognosis (HR = 1.089, p = 0.411; Table S3) by multivariate Cox regression analysis of KCNQ1OT1. Thus, MAP1B may act as an independent prognostic factor for colon cancer with perineural invasion patients.

Con rmation of MAP1B abnormally high expression
For better comprehension, 98 paired colon samples from the GEO: GSE44076 colon cancer cohorts in the GEO database was performed the difference analysis of MAP1B in colon cancer patients. Our data revealed that the expression of MAP1B in the colon cancer tissues was signi cantly higher than that of the paired nontumorous colon tissues, which was consistent with our results from TCGA data( Figure  S2D) . Furthermore, to analyze whether the genetic mutation act as a potential mechanism to promote the abnormal high-expression of MAP1B, we use an OncoPrint plot to show the ampli cation of MAP1B gene in TCGA colon cancer dataset. However, MAP1B expression is not signi cantly associated with copy number values among colon cancer samples( Figure S2B-2C). Hence, the results revealed that there is no remarkable correlation between MAP1B abnormal expression and copy number alterations in colon cancer.

MAP1B-related functional enrichment analysis in colon cancer with perineural invasion
To further explore the possible function of MAP1B in colon cancer with perineural invasion, we performed GO and KEGG pathway enrichment analysis of MAP1B in colon cancer (the top 200 correlated genes).

Discussion
Colon Cancer is a refractory disease with high morbidity and mortality worldwide. [16][17][18] Although surgical resection can improve the prognosis of some colon cancer patients, the prognosis of colon cancer with perineural invasion is still poor. [19][20][21] Therefore, exploring the biology and novel prognostic biomarkers of Colon Cancer with perineural invasion may provide clinicians new tools to treat the disease.
As far as the hypothesis of a ceRNA network postulates, lncRNAs, mRNAs, and other types of RNAs regulate each other's expression and affect tumorigenesis and cancer progression through competing for binding to miRNAs by sharing MREs [22]. The ceRNA regulatory network can better explain the interaction among a variety of RNA types at the genetic level. Therefore, we attempted to construct a MDK-related ceRNA triple network in colon cancer with perineural invasion, as well as correlate it with prognosis. Since increasing evidence has demonstrated that MDK acts as a neurogenesis suppressor functionally involved in many types of cancers. [9,[23][24][25] Interestingly, we veri ed similar results in this study through correlation analysis, survival analysis and IHC analysis.
In our study, rst, 17 lncRNAs, 6 miRNAs, and 36 mRNAs was obtained to establish the lncRNA-miRNA-mRNA triple regulatory network through differential analysis. Enrichment analysis identi ed that the network was mainly associated with "extracellular structure organization," "cell-cell junction," and "glycosaminoglycan binding." Next, a key triple regulatory network, including six lncRNAs, six miRNAs, and three mRNAs ,was obtained by hub analysis based on a score >2. Finally, we performed expression analysis and survival analysis on this hub regulatory network. Moreover, To con rm the localization of the ceRNA network, we also performed subcellular localization analysis of the six lncRNAs. In a word, the KCNQ1OT1-miR-454-3p/miR-204-5p-MAP1B highexpressed ceRNA network related to the prognosis of colon cancer with perineural invasion was obtained.
To further explore these genes in PubMed, we observed that miR-454-3p, miR-204-5p, and KCNQ1OT1 have been researched for their functions in cancer or in the relationship with cancer. Li et al. [26] provided results suggesting that miR-454-3p is overexpressed in the tumor tissues in patients with colorectal cancer and participates in the progression of colorectal cancer by promoting cancer cell proliferation, invasion, and liver metastasis. However, we found that low expression of miR-454-3p executed the tumor suppressive activity via downregulating NFATc2 in Glioblastoma [27] .Wa et al. [28] study showed that miR-204-5p expression is reduced in PCa tissues and serum sample with bone metastasis compared with that in PCa tissues and serum sample without bone metastasis, which is associated with advanced clinicopathological characteristics and poor bone metastasis-free survival in PCa patients Additionally, another study suggested that KCNQ1OT1 promotes cell proliferation and autophagy and inhibits cell apoptosis via regulating miR-204-5p/ATG3 axis in Non-small cell lung cancer. [29] However, few studies have investigated the role of MAP1B in cancer with perineural invaion according to searches in PubMed.
Microtubule-associated protein 1B (MAP1B) contains both actin and MT binding domains. MAP1B is synthesized as a polyprotein which is proteolytically cleaved to generate a light chain (LC1) and a heavy chain (HC), and is predominantly expressed in the nervous system, where it has been described as the rst microtubule-associated protein (MAP) to be expressed during early stages of development. Most MAP1B functions are related to events occurring during development of the nervous system. Thus, MAP1B regulates axonal elongation and guidance, and neuronal migration. Nonetheless, MAP1B expression persists during adult stages, speci cally in brain areas showing high synaptic plasticity like the olfactory bulb, hippocampus and cerebellum. [30] Moreover, The acidic extracellular microenvironment downregulate MAP1B expression level of the colorectal cancer cell SW620 that was cultured in acidic medium (pH 6.5) for more than 3 months to be acidosis-adapted. The study reveal that the tumorigenicity ,increased invasion and liver metastasis capacity of SW620-AA cells was dramatically enhanced compared with SW620-NA cells. [31] In this study, KCNQ1OT1 and MAP1B has been proven overexpressed in colon cancer, and high expression level of KCNQ1OT1 and MAP1B indicated worse outcome by survival analysis. However, miR-454-3p and miR-204-5p expression were downregulated in colon cancer samples, and the low expression level of miR-454-3p and miR-204-5p was correlated with poor OS prognosis in colon cancer.
By further analyzing the biological role of MAP1B, GO and KEGG enrichment analysis for MAP1B was performed. Axoneme assembly act as the enrichment term related to MAP1B. Extracellular structure organization, cell-cell junction and glycosaminoglycan binding are mainly related to MAP1B by GO enrichment analysis (including BP, CC, and MF). Therefore, MAP1B may be a key gene that promotes nerve growth in the process of colon cancer with perineural invasion.
Therefore, we have demonstrated a potential prognostic biomarker in clinical applications by constructing the ceRNA-based KCNQ1OT1/ MAP1B axis. However, several limitations must also be concerned. First, further experimental investigation should be performed to verify the binding a nity of lncRNAs, miRNAs, and mRNAs obtained from the database. Second, further study for the function and mechanism of the KCNQ1OT1/ MAP1B axis need to be carried out in colon cancer with perineural invasion by experiments.

Conclusion
In a word, we constructed a MDK-related ceRNA (KCNQ1OT1-miR-454-3p/miR-204-5p-MAP1B) highexpressed network and veri ed potential prognostic markers associated with colon cancer with perineural invasion. One hub gene in the prognostic signature, MAP1B, was identi ed as a novel potential biomarker for colon cancer with perineural invasion through data mining and experiments.

Declarations
Availability of data and materials All the data referred to this study are available from TCGA and GEO databases.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.