Identification of novel survival-related lncRNA-miRNA- mRNA competing endogenous RNA network associated with immune infiltration in colorectal cancer

Increasing studies have reported that long noncoding RNAs (lncRNAs) play critical roles in the initiation and progression of carcinogenesis. However, the underlying regulatory mechanisms of lncRNA-related competing endogenous RNA (ceRNA) network in colorectal cancer (CRC) are not fully understood. In the present study, we systematically analyzed the expression levels and prognostic values of dysregulated microRNAs (miRNAs) in human CRC to identify novel survival-related lncRNA-miRNA-mRNA ceRNA regulatory network. As a result, 28 dysregulated miRNAs were obtained, and hsa-miR-195-5p was identified as a key oncogene in human CRC based on analyses of expression levels and prognostic values. By means of stepwise prediction and validation, two upstream lncRNAs (NEAT1, XIST) and eight downstream mRNAs (ACOX1, CYP26B1, IRF4, ITPR1, LITAF, PHLPP2, RECK, and TPM2) were identified as key genes that interact with hsa-miR-195-5p. A ceRNA regulatory network consisted of these key genes was constructed, and Gene Set Enrichment Analysis (GSEA) indicated the possible association of key mRNAs with CRC onset and progression. Importantly, immune infiltration analysis revealed that the ceRNA network was remarkably associated with infiltration abundance of multiple immune cells and expression levels of immune checkpoints. These findings indicate that NEAT1 and XIST are potential prognostic factors that affect CRC onset and progression by targeting miR-195-5p.


Introduction
Malignant tumor, with more than 18.0 million newly diagnosed cases and 9.0 million deaths in 2018, has become the main threat to public health around the world [1]. Among them, colorectal cancer (CRC) is one of the most frequently diagnosed digestive malignancies, which remains the third leading cause of cancer-related death worldwide [1,2]. Early diagnosis and timely radical resection are critical for CRC. The 5-year survival rate is approximately 90% when the localized disease was diagnosed at an early stage. However, the 5-year survival rate of patients with distant metastasis dropped sharply to only 12% [3]. Metastasis and recurrence are the major causes of death in CRC patients. It is estimated that more than 50% of CRC patients will die from metastasis-related complications, which might be due to the lack of specific and sensitive biomarkers [4]. To make matters worse, available epidemiological data indicate that the incidence rates of human CRC have increased faster than any other types of cancer around the world [5]. Despite advances in surgical technology and targeted therapy, the prognosis of patients with CRC remains dismal because of recurrence and distant metastasis. Therefore, exploring the molecular mechanisms underlying colorectal carcinogenesis is imperative and meaningful to provide novel prognostic and therapeutic biomarkers for clinical CRC treatment.
The occurrence of CRC is a complex multi-step process involving accumulation of various genetic and epigenetic alternations [6]. In recent years, increasing evidence has indicat-ed that long noncoding RNA (lncRNA) plays a significant role in the development of multiple types of malignant tumor, including CRC [7][8][9][10]. LncRNA is a class of non-coding transcripts with more than 200 nucleotides in length [11]. Salmena et al [12] first proposed the competing endogenous RNA (ceRNA) hypothesis that genes could achieve crosstalk among each other by forming a regulatory network. The ceRNA hypothesis postulates that lncRNAs could act as endogenous miRNA sponges and inhibit miRNA function, thereby modulating miRNA-mediated target inhibition. In other words, the expression level of miRNA could be negatively correlated with its upstream lncRNA and downstream mRNA simultaneously [13]. Therefore, ceRNA crosstalk is a vital mechanism underlying the complex pathogenesis and multistage progression of human CRC, but our knowledge about the lncRNA related ceRNA network in human CRC remains limited, exploring the ceRNA regulatory network might be a mately, the expression levels and prognostic values of the upstream lncRNAs and downstream mRNAs in the ceRNA regulatory network were assessed, and a novel ceRNA network in which each RNA was significantly correlated with the prognosis of patients with CRC was constructed successfully. Furthermore, we used the Gene Set Enrichment Analysis (GSEA) software package to explore the pathway enrichment of this survival-related ceRNA regulatory network. Finally, Tumor Immune Estimation Resource (TIMER) database was used to estimate the correlations between ceRNA network and immune infiltration level in CRC microenvironment ( Figure 1). Through our efforts, it is hoped that the results of our research will contribute to enhancing our understanding of the regulatory mechanisms underlying the pathogenesis of human CRC and finding new therapeutics for clinical CRC treatment. potential pointcut for identifying novel early diagnostic and therapeutic targets of CRC and deserves further research.
In the present study, dysregulated miRNAs in human CRC were identified by mining the gene expression profiles obtained from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) database. Subsequently, two miRNAs were identified as key miRNA for further study through comprehensive consideration of the expression levels and prognostic values of these miRNAs in human CRC. Then, the upstream lncRNAs and downstream mRNAs of the key miRNAs were predicted by using multiple bioinformatic tools, and a ceRNA regulatory network was constructed and visualized by using Cytoscape software. In addition, the Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the ceRNA regulatory network were performed by using the DAVID database. Ulti-

Data download and preprocessing
To screen out the dysregulated miRNAs in human CRC, human miRNA expression profile GSE108153 (including 21 cancer samples and 21 adjacent normal samples) [14] was downloaded from the Gene Expression Omnibus database (GEO, http://www.ncbi.nlm.nih.gov/ geo) [15], which was based on the GPL197-30 platform (Agilent-046064 Unrestricted_Hu-man_miRNA_V19.0_Microarray). Meanwhile, to enhance the reliability of our results, miRNA expression profiles (including 9 normal samples and 539 cancer samples) and clinical information (including 547 cases) of patients with CRC were obtained from The Cancer Genome Atlas (TCGA, https://gdc-portal.nci. nih.gov/) database.

Identification of dysregulated miRNAs
First, we converted the miRNA probe ID in the matrix files to mature miRNA name based on the annotation document of its platform files through the Perl language (https://www.perl. org/). Subsequently, the "Limma" package in R software (V 3.6.2; http://www.r-project.org/) was used to identify differentially expressed miRNAs between CRC samples and normal samples both in miRNA expression profiles obtained from GEO and TCGA databases. Adjusted P-value of < 0.05 and absolute value of log2 fold change (log2 FC) ≥1 were considered as the cut-off criteria [16]. Finally, the differentially expressed miRNAs dysregulated in both GEO and TCGA datasets were chosen for subsequent analysis by using the intersect function of the Venn Diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Identification of key miRNAs
To identify the key miRNAs associated with the prognosis of CRC patients, the expression levels and prognostic values of the dysregulated miRNAs were assessed by using clinical information obtained from the TCGA database. The dysregulated miRNAs not only differentially expressed between the CRC cancer samples and normal samples, but also associated with the prognosis of CRC patients were identified as key miRNAs in human CRC and were chosen for subsequent research, according to the criteria of a P-value < 0.05. In addition, we further validated the expression pattern of key miRNAs in multiple types of cancer by using the Gene Expression Display Server (GEDS, http://bioinfo.life.hust.edu.cn/web/GEDS/) database, a web-based tool for quantification, comparison, and visualization of gene expression data based on The Cancer Genome Atlas, Genotype-Tissue Expression (GETx), Cancer Cell Line Encyclopedia (CCLE), and MD Anderson Cell Lines Project (MCLP) databases [17].

Construction of ceRNA regulatory network and functional enrichment analysis
The starBase database (http://starbase.sysu. edu.cn/) is an online tool that allows researchers to perform Pan-Cancer analysis on RNA-RNA and RBP-RNA interactions [18]. The Lnc-Base V2 database (http://carolina.imis.athenainnovation.gr/) is a web-based tool which caters information regarding cell type specific miRNA-lncRNA interaction and enables users to easily identify interactions in 66 different cell types, spanning 36 tissues for human and mouse [19]. The upstream lncRNAs which interacted with the key miRNAs were predicted by using these two databases. Only those lncRNAs predicted by both two databases and interacted with all key miRNAs simultaneously were identified as upstream lncRNAs. Subsequently, the downstream mRNAs which interacted with the key miRNAs were predicted by using miRTar-Base (http://mirtarbase.mbc.nctu.edu.tw/) [20], targetScan (http://www.targetscan.org) [21], and miRDB (http://www.mirdb.org/) [22]. To enhance the reliability of the predicted interaction relationships, only the predicted mRNAs presented in all three databases were identified as the downstream mRNAs. Ultimately, the lncRNA-miRNA-mRNA regulatory network was constructed and visualized by Cytoscape software (v3.7.2, http://www.cytoscape.org/) [23]. Furthermore, the functional enrichment analysis of this ceRNA regulatory was performed by using the Database for Annotation, Visualization, and Integrated Discovery (DAVID v6.8 https://david.ncifcrf.gov/home.jsp) [24] bioinformatics database. Only the terms with a P-value < 0.05 were considered as statistically significant.

Gene expression validation and survival analysis
To further enhance the reliability of the predicted interaction relationships, the expression lev-els of upstream lncRNAs were assessed by using the Gene Expression Profiling Interactive Analysis (GEPIA; http://gepia.cancer-pku.cn/ detail.php), a web-based tool that provides fast and customizable functionalities including differential expression analysis, profiling plotting, correlation analysis, and patient survival analysis based on TCGA and Genotype-Tissue Expression data [25]. The expression levels of downstream mRNAs were assessed by using the gene expression profiles in human CRC obtained from the TCGA database. Subsequently, the prognostic values of the upstream lncRNAs and downstream mRNAs were validated in human CRC by using the PrognoScan (http://dna00.bio.kyutech.ac.jp/PrognoScan/) database, a web-based tool for assessing the biological relationship between gene expression and prognosis of multiple types of cancer [26]. Only genes that complied with the ceRNA hypothesis and associated with the prognosis of patients with CRC were chosen for subsequent analysis.

Identification of survival related ceRNA regulatory network and functional annotation
Based on the prediction and evaluation above, a novel lncRNA-miRNA-mRNA regulatory network associated with the prognosis of CRC patients was constructed and visualized by Cytoscape software. In addition, to further explore the underlying role of this survival-related ceRNA regulatory network in the process of human CRC, functional annotation of genes in the ceRNA regulatory network was analyzed by the Gene Set Enrichment Analysis (GSEA, https://www.gsea-msigdb.org/) database, with an initial database of 1,325 biologically defined gene sets [27]. The enriched terms with P-value of < 0.05 were considered statistically significant.

Correlation of ceRNA network with immune infiltration level in human CRC
The Tumor Immune Estimation Resource (TIMER, https://cistrome.shinyapps.io/timer/) is an online database established for systematical analysis of immune infiltrates across multiple cancer types [28]. We applied the TIMER database to assess the correlations of gene expression in the ceRNA regulatory network with the level of immune infiltrates in human CRC. In addition, the correlations between the ceRNA regulatory network and abundance of T-cell infiltration in human CRC were confirmed by using the GEPIA database.

Statistical analysis
Most of the statistical analyses were carried out using the bioinformatic tools mentioned above, while the rest were carried out using R software (v 3.6.2). Identification of differentially expressed genes was assessed by Student's t-test, and Benjamini-Hochberg False Discovery Rate (FDR) method was applied to adjust the P-value. Kaplan-Meier curve analysis with log-rank test was used to evaluate the prognostic values of miRNAs in human CRC. Spearman correlation coefficients were used to evaluate the correlations. A P-value < 0.05 was considered statistically different.

Identification of dysregulated miRNAs in human CRC
As shown in the volcano plot and heatmap, differentially expressed miRNAs in GEO and TCGA datasets were identified separately according to the predefined threshold ( Figure 2A-D). In the GSE108153 dataset, a total of 49 dysregulated miRNAs were identified. Meanwhile, a total of 504 dysregulated miRNAs were also screened out in the TCGA dataset. A total of 28 dysregulated miRNAs were selected for subsequent analyses by using the intersect function of the Venn Diagram ( Figure 2E).

Identification of key miRNAs
To evaluate the prognostic values of dysregulated miRNAs for overall survival in human CRC, clinical information, including survival status and survival time of 547 CRC patients, were obtained from the TCGA database. The results showed that five of the dysregulated miRNAs (miR-195-5p, miR-10b-5p, miR-21-3p, miR-125b-5p, and miR-145-5p) were significantly associated with the prognosis of CRC patients (P < 0.05; Figure 3A-E). Meanwhile, we compared the expression levels of these survivalrelated miRNAs between cancer samples and normal samples in the TCGA dataset. After combining the prognostic value and expression level, we found that only two of them (miR-195-5p and miR-10b-5p) showed both high expression and poor prognosis in human CRC (P < 0.05; Figure 3F-J). To further enhance the reliability of our screening results, we verified the two dysregulated miRNAs in the GEDS database. As shown in Figure 3K and 3L, miR-195-5p and miR-10b-5p were dysregulated in multiple cancer types, including CRC. Thus, those two miRNAs were chosen for subsequent analyses.

Construction of ceRNA regulatory network
The upstream lncRNAs that could potentially bind to the key miRNAs were predicted through starBase and LncBase V2 databases, and a total of three upstream lncRNAs were predicted to regulate the two key miRNAs simultaneously ( Figure 4A). Subsequently, to identify downstream mRNAs that could potentially be targeted by those two key miRNAs, we applied miR-TarBase, targetScan, and miRDB to improve the credibility of the predicted results. As a result, a total of 287 downstream mRNAs were predicted to be targeted by miR-195-5p, while there were 22 downstream mRNAs could potentially be targeted by miR-10b-5p ( Figure  4B, 4C). Based on the prediction above, a ceRNA regulatory network comprised of six lncRNA-miRNA pairs and 310 miRNA-mRNA pairs was established and visualized by Cytoscape software ( Figure 4D).

Functional analysis for the ceRNA regulatory network
We used the DAVID database to further explore the biological functions of the ceRNA regulatory network. A total of 222 GO terms and 45 KEGG pathways were enriched, and the top 10 enriched GO terms of each group and the top   A. The intersection of upstream lncRNAs targeted by hsa-miR-9-5p and hsa-miR-10b-5p in LncBase v2 and starBase databases. B, C. The intersection of downstream mRNAs targeted by hsa-miR-9-5p and hsa-miR-10b-5p in miRTarBase, Targetscan, and miRDB databases. D. The lncRNA-miRNA-mRNA regulatory network constructed by Cytoscape software. The green round rectangle in the network represents miRNA. The pink ellipse in the network represents mRNA. The red diamond in the network represents lncRNA. 20 enriched KEGG pathways according to P-value were displayed (Tables 1, 2). In the biological processes (BP) group, the ceRNA regulatory network is mainly enriched in the process of regulation of transcription from RNA polymerase II promoter, regulation of cell migration, regulation of cell adhesion, and regulation of the apoptotic process. In the cellular component (CC) group, the ceRNA regulatory network is mainly enriched in nucleoplasm, cytoplasm, and nucleus. In the molecular function (MF) group, the ceRNA regulatory network is mainly enriched in protein binding, beta-catenin binding, and ubiquitin-protein ligase activity. These results indicated that the ceRNA network was significantly enriched in terms correlated with cell proliferation, cell migration, and cell apop-tosis. The KEGG enrichment analysis revealed that the ceRNA regulatory network was significantly enriched in multiple cancer-related pathways such as PI3K-Akt signaling pathway, FoxO signaling pathway, signaling pathways regulating pluripotency of stem cells, and AMPK signaling pathway.

Identification and validation of key lncRNA and mRNA
Based on the ceRNA hypothesis, the miRNA could be negatively correlated with its upstream lncRNA and downstream mRNA simultaneously. Therefore, we assessed the expression levels and prognostic values of the predicted lncRNAs for overall survival in patients with CRC by using GEPIA and PrognoScan databases, respectively. Only two lncRNAs (NEAT1 and XIST) showed both low expression and improved prognosis in patients with CRC ( Figure 5A, 5B). As for downstream mRNA, we collected gene expression profiles in human CRC from the TCGA database to assess the expression levels of the predicted mRNA. Meanwhile, the prognostic values of those mRNAs for overall survival in patients with CRC were also evaluated using the PrognoScan database. A total of 1964 downregulated genes were identified in human CRC from the TCGA database, and we then integrated the 1964 downregulated genes with the downstream mRNAs in the ceRNA regulatory network. As a result, no specific gene was identified to be targeted by miR-10b-5p, while there were 16 specific genes could potentially be targeted by miR-195-5p ( Figure 5C,  5D).
Subsequently, survival analysis found that eight of the 16 specific genes (ACOX1, CYP26B1, IRF4, ITPR1, LITAF, PHLPP2, RECK, and TPM2) were significantly associated with the prognosis of patients with CRC. The expression boxplot and survival curve of those mRNAs are displayed in Figure 5E-L. Those two lncRNAs and eight mRNAs were chosen for further analysis.

Integration of survival-related ceRNA regulatory network and GSEA enrichment analysis
Based on the validation above, a survival-related ceRNA regulatory network including two lncRNA-miRNA pairs (NEATI-miR-195-5p and XIST-miR-195-5p) and eight miRNA-mRNA pairs (miR-195-5p-ACOX1/CYP26B1/IRF4/ITPR1/ LITAF/PHLPP2/RECK/TPM2) was established ( Figure 6A). Each component in the ceRNA regulatory network had significant prognostic values in human CRC, and fully complied with the rule of the ceRNA hypothesis. After that, we used the GSEA software package to further explore the biological function of the survival- GOTERM_MF_DIRECT GO: 0043565~sequence-specific DNA binding 20 0.001444 related ceRNA regulatory network. As shown in Figure 6B, GSEA enrichment analysis indicated that the survival-related ceRNA regulatory network was significantly enriched in multiple cancer-related pathways, such as the JAK-STAT signaling pathway, mTOR signaling pathway, TGF-β signaling pathway, Wnt signaling pathway, ERBB signaling pathway, and VEGF signaling pathway.

Correlation of ceRNA network with immune infiltration level in human CRC
The correlation between the survival-related ceRNA regulatory network and tumor-infiltrating immune cells in the CRC microenvironment was assessed by using TIMER. As shown in Figure  7A-H, we found that gene expression in the ceRNA regulatory network was significantly associated with infiltrating immune cells, including B cells, CD4 + /CD8 + T cells, macrophages, meutrophils, and dendritic cells, in the human CRC microenvironment. Subsequently, the correlations between the ceRNA regulatory network and the abundance of T cells infiltration in human CRC were further evaluated by using the GEPIA database, and the results indicated that the survival-related ceRNA regulato-  Figure 8I-L). Therefore, these findings further support that this survival-related ceRNA regulatory network was closely correlated with immune infiltration level in human CRC, suggesting that the ceRNA regulatory network might be an important participant in immune escape in the CRC microenvironment.

Discussion
Accumulating studies have reported that lncRNA acts as a post-transcriptional regulator   that regulates downstream mRNA expression by sponging miRNA, and participates in multiple processes of tumorigenesis, including CRC. For instance, lncRNA CRNDE has been reported to function as a sponge for miR-181a-5p, thereby promoting the progression and chemoresistance of CRC cells through modulating the expression levels of TCF4 and the activity of Wnt/β-catenin signaling [29]. Previous research reported that lncRNA SNHG1 promotes CRC cell growth by sponging miR-154-5p to regulate the expression level of CCND2, which plays a critical role in cell cycle progression [30]. LncRNA MIR17HG was reported to increase the expression level of NF-κB/ RELA through competitively sponging miR-375, thereby promoting tumorigenesis and metastasis in CRC cells. In addition, MIR17HG was reported to upregulate the expression of PD-L1, suggesting that it may participant in immune escape in the CRC microenvironment [31]. Nevertheless, the clinical utility of lncRNA-related ceRNA regulatory network is limited at present. Exploring the lncRNA related ceRNA regulatory network will help to provide a comprehensive perspective of the molecular mechanisms of gene interaction and regulation in human CRC.
In the present study, dysregulated miRNAs with prognostic significance in human CRC were identified in cancer samples compared with adjacent non-cancer samples from TCGA and GEO databases. Then, a novel lncRNA-miRNA-mRNA regulatory network involved in CRC was constructed through stepwise prediction and validation. Functional enrichment analysis revealed that this novel ceRNA regulatory network was significantly enriched in some cancerrelated pathways, such as PI3K-Akt signaling pathway [32], FoxO signaling pathway [33], Wnt signaling pathway [34], and AMPK signaling pathway [35]. To enhance the reliability of our study, we further assessed the expression levels and prognostic values of genes in the ceRNA regulatory network and constructed a survivalrelated ceRNA regulatory network; each constituent in the network has prognostic significance in human CRC, which might provide novel diagnostic and therapeutic targets for clinical treatment of CRC.
LncRNAs are defined as transcripts that are > 200 nucleotides in length and have no ability to code protein [36]. Previous studies have reported that lncRNA participates in multiple biological processes and tumorigenesis by sponging miRNA [37]. In the present study, we found that lncRNA NEAT1 and XIST were not only dysregulated in CRC, but also significantly associated with poor prognosis, and these two lncRNAs were identified as key lncRNA in CRC. Studies have shown that expression of lncRNA XIST was dysregulated in CRC cell lines and tissues. XIST knockdown inhibited the cell migration and invasion of CRC cells in vitro and in vivo by regulating the processes of epithelial-mesenchymal transition (EMT), an important step in tumor progression and metastasis [38,39]. In addition, dysregulated XIST has been found to be involved in chemoresistance of CRC to Doxorubicin, with the silencing of XIST significantly enhancing the anti-tumor effect of Doxorubicin in CRC [40]. There is also increasing evidence that lncRNA NEAT1 was dysregulated in human CRC tissues and was significantly correlated with worse overall survival in CRC patients. NEAT1 regulates the Wnt/βcatenin signaling pathway by targeting DDX5, thereby facilitating CRC cell migration and invasion in vitro and in vivo [41]. Furthermore, NEAT1 also has been found to participate in chemoresistance of CRC, and NEAT1 knockdown increased the sensitivity of 5-fluorouracil and promoted apoptosis in CRC cells [42]. These studies above improved the reliability of our research finding that dysregulation of lncRNA XIST and NEAT1 was closely correlated with the pathogenesis of CRC, which is valuable for further study.
MiRNAs are defined as endogenous non-coding RNAs that are 21-25 nucleotides in length and suppress protein-coding gene expression by binding to the 3'-untranslated regions (UTRs) of mRNAs [43]. Dysregulation of miRNA has been reported to be associated with the onset and progression of human CRC. For instance, expression of miR-330 was downregulated in CRC samples compared with normal control, and CRC patients with low miR-330 expression showed poorer overall survival than those with high miR-330 expression. As of mechanism, miR-330 suppresses cell proliferation, migration, invasion, and angiogenesis by reducing the phosphorylation of AKT and STAT3 via targeting HMGA2 in vitro [44]. Moreover, expression of miR-128-3p in tumor samples was closely associated with oxaliplatin sensitivity in human CRC, and forced-expression of miR-128-3p significantly improved the sensitivity of CRC cells to oxaliplatin [45]. In the present study, we found that miR-195-5p was sponged by XIST and NEAT1, it was significantly associated with the prognosis of patients with CRC, and miR-195-5p was identified as a potential biomarker in CRC. Studies have shown that expression of miR-195-5p was dysregulated in CRC, and miR-195-5p acts as an independent risk factor for overall survival in CRC [46]. Importantly, miR-195-5p has been proved to regulate the stemness and chemoresistance of CRC cells through inhibiting Notch signaling pathway, and forced-expression of miR-195-5p significantly improved the sensitivity of CRC cells to 5-fluorouracil [47]. Intriguingly, we also found that the role of NEAT1/XIST-miR-195-5p regulatory axes in human cancer has been proved by previous studies, which partially enhanced the reliability of our results [48,49]. However, studies focused on the significance of NEAT1/XIST-miR-195-5p regulatory axes in CRC remains largely limited, and they are worthy of further exploration.
We further analyzed and proved the downstream mRNAs that were targeted by the key miRNA, and eight survival-related mRNAs were identified as key mRNAs in CRC that were regulated by miR-195-5p. In agreement with our results, previous studies have reported that these key mRNAs participate in the initiation and progression of multiple types of cancers, including CRC. IRF4, as a member of the interferon regulating factor family, is downregulated by PIP5K1A in human CRC [50]. PHLPP2 belongs to the phosphokinase family that has been found to be down-regulated in CRC, and CRC patients with high PHLPP2 expression showed a better prognosis than those with low PHLPP2 expression [51]. TPM2 has been found to be down-regulated in CRC cell lines and tis-sues, and reduction of TPM2 was associated with increased migration and proliferation of CRC cells in vitro [52]. RECK is a unique matrix metalloproteinase (MMP) regulator that has been found to be lowly expressed in CRC tissues and cell lines. RECK knockdown enhanced the migratory and invasive rates of CRC cells [53]. ITPR1 is an isoform of inositol 1,4,5-trisphosphate receptor (ITPR), and there is increasing evidence that ITPR expression was associated with aggressiveness of CRC cells [54]. CYP26B1 belongs to the CYP26 enzyme hydroxylate retinoic acid family that has been reported to be dysregulated in CRC, and CYP26B1 expression was significantly correlated with the clinical outcome of patients with CRC [55]. ACOX1, which catalyzes the initial step for peroxisomal β-oxidation, has been reported to associate with peroxisomal disorders and carcinogenesis in liver cancer [56]. LITAF has been reported to function as a tumor suppressor and is frequently down-regulated in multiple cancer types because its expression is regulated by the tumor suppressor protein, p53 [57]. However, there is currently no research that focused on the significance of ACOX1 and LITAF in human CRC, so they are worthy of further study. Furthermore, GSEA analysis revealed that all of these key mRNAs were significantly enriched in well-known tumorigenesis-related pathways, such as the JAK-STAT signaling pathway [58], TGF-β signaling pathway [59], ERBB signaling pathway [60], and VEGF signaling pathway [61]. These reports partially improved the credibility of our research, and further research on the significance of the novel ceRNA network in CRC is valuable.
Tumor microenvironment (TME), a complex system composed of tumor cells, stromal cells, immune cells, fibroblasts, blood vessels, and cytokines, has been reported to regulate some well-known oncogenic process including apoptosis, angiogenesis, hypoxia, and immune escape [62]. Immune cells are important components of TME, which induced the host immune response by secreting cytokines, thereby inhibiting or promoting the progression of tumor cells [63]. For example, the regulatory T cells (Tregs) secrete immunosuppressive cytokines TGF-β and interleukin-10 (IL-10) to protect cancer cells from cytotoxic T cell response and induce immune tolerance, while CD8 + T cells secrete interleukin-2 (IL-2), tumor necrosis factor α (TNF-α) and interferon-γ (IFN-γ) to kill cancer cells in the TME [64]. Immunotherapy has shown encouraging results as an emerging therapy for some cancers. However, existing evidence suggested that most CRC patients were insensitive to immunotherapy compared with esophageal cancer and lung cancer. Therefore, exploring the infiltration abundance and function of immune cells is important to enhance the efficacy of immunotherapy in CRC. In the present study, we found that the ceRNA regulatory network is significantly associated with the infiltration of immune cells, including CD4 + /CD8 + T cells and macrophages. Importantly, we found that the ceRNA regulatory network is significantly correlated with the expression of some well-known immune checkpoints. Accumulating studies have reported that tumor-infiltrating lymphocytes, including CD4 + and CD8 + T cells, are important components of immune cells with anti-tumor functions. Furthermore, immune cells are highly heterogeneous in TME compared with noncancerous areas, and infiltration of tumor-infiltrating lymphocytes was associated with improved survival in some cases [63]. Studies also have indicated that ligands of immune checkpoints are highly expressed in cancer cells and TME, leading to T cells' functional exhaustion and T cell failure, thereby inducing immune escape of CRC cells in the TME [65]. Based on previous studies and our results above, we speculate that the survivalrelated ceRNA regulatory network identified by our study may play an important role in immune escape of cancer cells in CRC TME. Further study on this survival-related ceRNA network will help to improve the efficacy of immunotherapy in CRC.
Some limitations inevitably existed in the present study. First and foremost, all conclusions of our study were based on online databases such as TCGA and GEO. Further studies will be needed to validate our findings, and to elucidate whether and how the survival-related ceRNA regulatory network regulates CRC. Second, we did not evaluate the subtype of CRC, such as mismatch repair (MMR) and microsatellite instability (MSI) status, which may affect the expression profiles and prognosis of CRC patients. Third, we cannot obtain all the clinicopathologic data of each patient, and the survival analysis was performed by online database automatically, which may lead to bias.
In conclusion, we successfully constructed a survival-related ceRNA regulatory network in CRC by means of stepwise prediction and validation, and each component of the ceRNA regulatory network was significantly correlated with the prognosis of patients with CRC. More importantly, immune infiltration analysis revealed that the survival-related ceRNA regulatory network was remarkably associated with infiltration abundance of multiple immune cells and expression levels of immune checkpoints. Our findings provide valuable clues into improving the efficacy of targeted therapy and immunotherapy for human CRC.

Disclosure of conflict of interest
None.