Identification of the circRNA-miRNA-mRNA regulatory network and its prognostic effect in colorectal cancer

Abstract

BACKGROUND

The high morbidity and mortality of colorectal cancer (CRC) have posed great threats to human health. Circular RNA (CircRNA) and microRNA (miRNA), acting as competing endogenous RNAs (ceRNAs), have been found to play vital roles in carcinogenesis. However, the biological function of ceRNAs in CRC pathogenesis and prognosis remains largely unexplored.

AIM

To identify the CRC-specific circRNA-miRNA-mRNA regulatory network and uncover the subnetwork associated with its prognosis.

METHODS

CircRNAs, miRNAs and mRNAs differentially expressed (DE) in CRC tissues were selected by expression file analysis in the Gene Expression Omnibus (GEO) database, and the downstream target molecules of circRNAs and miRNAs were predicted. Then, the intersection of differentially expressed RNA molecules with the predicted targets was determined to obtain a ceRNA network. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted to elucidate the possible mechanism of pathogenesis. A survival analysis using the gene profiles and clinical information in The Cancer Genome Atlas (TCGA) database was performed to identify the mRNAs associated with the clinical outcome of CRC patients and construct a prognostic subnetwork.

RESULTS

We downloaded three datasets (GSE126095, GSE41655 and GSE41657) of large-scale CRC samples from the GEO database. There were 55 DEcircRNAs, 114 DEmiRNAs and 267 DEmRNAs in CRC tissues compared with normal tissues. After intersecting these molecules with predicted targets, 19 circRNAs, 13 miRNAs and 28 mRNAs were chosen to develop a circRNA-miRNA-mRNA network. GO and KEGG functional enrichment analyses indicated that the retinol metabolic process, leukocyte chemotaxis, extracellular matrix remodeling, endoplasmic reticulum stress, alcohol dehydrogenase activity, gastric acid secretion, nitrogen metabolism and NOD-like receptor signaling pathway might participate in the tumorigenesis of CRC. After verifying the identified mRNA effect in the TCGA database, we finally recognized 3 mRNAs (CA2, ITLN1 and LRRC19) that were significantly associated with the overall survival of CRC patients and constructed a ceRNA subnetwork including 5 circRNAs (hsa_circ_0080210, hsa_circ_0007158, hsa_circ_0000375, hsa_circ_0018909 and hsa_circ_0011536) and 3 miRNAs (hsa-miR-601, hsa-miR-671-5p and hsa-miR-765), which could contain innovative and noninvasive indicators for the early screening and prognostic prediction of CRC.

CONCLUSION

We proposed a circRNA-miRNA-mRNA regulatory network closely associated with the progression and clinical outcome of CRC that might include promising biomarkers for carcinogenesis and therapeutic targets.

Key Words: CircRNA; miRNA; Network; Colorectal cancer; Prognosis; Biomarkers

Core Tip: The biological functions of circRNA and miRNA interactions and their potential as noninvasive biomarkers have not been well elucidated in colorectal cancer (CRC). In this study, we constructed a circRNA-miRNA-mRNA regulatory network with 19 circRNAs, 13 miRNAs and 28 mRNAs. GO and KEGG analyses indicated several signaling pathways probably involved in tumorigenesis. After being combined with survival analysis, a prognostic subnetwork was constructed including 5 circRNAs, 3 miRNAs and 3 mRNAs, which may represent novel diagnostic and prognostic candidate biomarkers, as well as therapeutic targets of CRC.

INTRODUCTION

Colorectal cancer (CRC) is one of the most common tumors worldwide. The global cancer burden report in 2018 showed that the incidence and mortality of CRC ranked third and second, respectively[1]. Conventional therapeutic options for CRC, such as chemotherapy and radiotherapy, cannot satisfy the ever-rising demand for overall and disease-free survival. At present, the high morbidity and mortality of CRC pose a great threat to the health of humans. Early screening and prevention should be actively and urgently carried out. The pathogenesis of CRC remains largely unexplored, and gene regulation disorders may play an important role in it.

In recent years, circular RNA (circRNA) and microRNA (miRNA), acting as noncoding RNAs, have attracted considerable research attention in a variety of diseases. These RNA molecules could regulate gene expression through complex mechanisms and interactions. CircRNA is a newly identified class of single-stranded circular, noncoding RNA molecules without 5′ poly-A and 3′ cap ends, which makes it resistant to degradation by RNA exonucleases and more stable than the linear RNA class[2]. Emerging evidence has proven that circRNAs are widely expressed in eukaryotic cells and can be implicated in physiological and pathological processes[3]. Salmena et al[4] first hypothesized that competing endogenous RNAs (ceRNAs) contain adequate miRNA response elements and may act as miRNA sponges to bind and compete with corresponding miRNAs, thereby sequestering miRNAs and regulating mRNA expression at the posttranslational level. MiRNAs are small, noncoding RNAs of approximately 20-22 nucleotides that can play a vital role in the regulation of gene expression, such as decreasing mRNA stability in various biological and pathological processes[5].

Accumulating evidence has revealed that the circRNA-miRNA-mRNA network could play a significantly important role in the development and progression of many diseases, especially cancer. For example, Song and Fu[6] discovered that the hsa_circ_00001666/hsa-mir-1229/CXCR5 axis could participate in the pathogenesis of CRC and act as a promising biomarker for targeted treatment. A related study indicated that hsa_circ_0005100 has pivotal value in the progression of CRC via the miR-1182/hTERT axis[7]. Hsa_circ_000984 could sequester miR-106b and consequently intensify the proliferation and migration of CRC cell lines[8]. These studies indicate that dysregulated circRNAs, miRNAs and mRNAs are closely related to the progression and prognosis of CRC and could be used as potential CRC-specific predictors, but the competitive regulatory pattern and biological function mechanism among circRNAs, miRNAs and mRNAs are still complicated and need further verification.

Along with the enormous advancement of RNA-sequencing technology, many public databases, such as Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA), have been established, and a surge of large-scale RNA sequence data are available. Recently, bioinformatics analysis has been widely used to help screen key genes, construct regulatory models, and select therapeutic targets that may participate in tumor development and prognosis as well as provide guidance for basic and clinical research, which is helpful to clarify the pathogenesis of tumors and guide clinical treatment options.

The objective of this study was to explore a competitive regulatory model among circRNAs, miRNAs and mRNAs and discover promising indicators for driving factors and mechanisms that induce the progression and affect the prognosis of CRC. To achieve this goal, we downloaded three datasets (GSE126095, GSE41655 and GSE41657) of large-scale CRC samples in the GEO database, selected differentially expressed (DE) circRNAs, DEmiRNAs and DEmRNAs in CRC tissues compared with normal controls, and predicted the downstream target molecules of circRNAs and miRNAs. After intersection with the differentially expressed RNA molecules and predicted targets, a circRNA-miRNA-mRNA network was identified. Functional enrichment analyses were conducted to identify the underlying mechanism involved in the pathogenesis of CRC. To verify the prognostic effect of the mRNAs found above, we performed survival analysis using the gene profiles and clinical information in the TCGA database. Finally, survival-related genes were determined, and a prognostic subnetwork was developed. Our results revealed the circRNA-miRNA-mRNA regulatory interaction and provided guidance for expanding the understanding of CRC progression, prognosis and therapeutic options. The flow chart of the procedure in our study is illustrated in Figure 1.

Figure 1 Flow chart of the procedure applied in this study. We identified the differentially expressed (DE) circRNAs, DEmiRNAs and DEmRNAs and then took the intersection of DEmiRNAs and DEmRNAs with the predicted targets of DEcircRNAs and intersected miRNA to obtain the competing endogenous RNA network in the Gene Expression Omnibus (GEO) database, performed GO and KEGG enrichment analyses, and constructed a prognostic subnetwork based on survival analysis in The Cancer Genome Atlas (TCGA) database. ceRNA: Competing endogenous RNA.

MATERIALS AND METHODS

Dataset retrieval

The expression profiles of three datasets (GSE126095, GSE41655 and GSE41657) were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo). Referring to the annotation information on the platform, probes were transformed to corresponding gene symbols. The GSE126095 dataset included 10 CRC and 10 normal tissues. The GSE41655 dataset contained 33 CRC tissues, 15 normal tissues and 59 colorectal adenomas, and GSE41657 included 25 CRC tissues, 12 normal tissues and 51 colorectal adenomas. CRC and normal tissue data in the two datasets above were chosen for comprehensive analysis. Moreover, clinical information and mRNA expression profiles were obtained from the TCGA database (https://portal.gdc.cancer.gov/). A total of 530 samples of patients (488 CRC tissues and 42 normal tissues), which had complete clinical characteristics, including age, sex, stage, survival time and survival state, and survival or follow-up time ≥ 30 d, were selected for further analysis. This study was approved by the Ethics Committee of China-Japan Friendship Hospital (No. 2018-116-K85-1). An informed consent statement was not necessary because all data were acquired from the GEO and TCGA databases and available to the public.

Identification of differentially expressed circRNAs, miRNAs and mRNAs

After downloading data from the GEO database, the “limma” package was utilized for background correction and normalization of the raw read counts of circRNA, miRNA and mRNA as well as identification of DEcircRNAs, DEmiRNAs and DEmRNAs between tumor samples and normal controls. The cut-off criteria of circRNA and mRNA were set at an adjusted P < 0.05 and |log₂ fold change (FC)| > 2, while the cut-off value of miRNA was set at an adjusted P < 0.05 and |log₂ FC|>1.

Target prediction and intersection for ceRNA network construction

In the circBase database (http://www.circbase.org/), specific information on the DEcircRNAs was available, and then the Cancer-Specific CircRNA database (CSCD) (http://gb.whu.edu.cn/CSCD/) was used to obtain structural patterns of circRNAs and predict the binding relationship between circRNAs and miRNAs. After intersection with the DEmiRNAs, promising miRNAs were finally identified, and their target mRNAs were predicted in the TargetScan and miRDB databases. When both databases supported mRNAs as candidate targets, these targets and DEmRNAs were intersected to obtain the final mRNA. After removing unconnected nodes, the circRNA-miRNA-mRNA network was developed based on the results above and visualized using Cytoscape 3.7.2.

GO and KEGG functional enrichment analysis

To reveal the pathophysiological processes and critical signaling pathways involved in the carcinogenesis of CRC, two widely used bioinformatics analysis methods, Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), were conducted. The “clusterProfiler” package in R/Bioconductor was used for analyzing GO term and KEGG pathway enrichment. An adjusted P < 0.05 was regarded as a statistical criterion.

Survival analysis and prognostic subnetwork construction

To assess the value of the identified ceRNA network and determine the mRNAs related to prognosis, the mRNA expression profile and clinical information for CRC patients were downloaded from the TCGA database. Kaplan–Meier curves were generated for survival analysis. The log-rank test was utilized for statistical analysis. The cut-off criterion was a P < 0.05. Finally, a ceRNA subnetwork was developed on the basis of the above verified mRNAs.

RESULTS

Identification of differentially expressed and intersecting circRNAs, miRNAs and mRNAs

The basic characteristics of the three GEO datasets (GSE126095, GSE41655 and GSE41657) and the TCGA is shown in Table 1. In the circRNA profile data of GSE126095, there were a total of 55 DEcircRNAs between the CRC samples and normal controls, among which 22 upregulated and 33 downregulated in CRC tissues (Figure 2). After removing 6 circRNAs that did not have details in the CSCD database, the structural models of 49 DEcircRNAs obtained from CSCD database are depicted in Figure 3, and the basic information of the these circRNAs is listed in Supplementary Table 1. Based on this, we successfully predicted that 1602 miRNAs might be the targets of these 49 circRNAs. For miRNAs with aberrant expression in GSE41655, 114 miRNAs were investigated to be differentially expressed in CRC, among which 58 were overexpressed and 56 were downregulated (Supplementary Table 2). A total of 25 miRNAs were identified from the intersection of 114 DEmiRNAs and 1602 circRNA targets predicted by the CSCD database (Supplementary Table 3). According to the TargetScan and miRDB databases, 7190 potential target genes for the 25 intersecting miRNAs were found. In the GSE41657 mRNA expression profile, 267 mRNAs were differentially expressed, among which 112 were highly expressed and 155 were expressed at lower levels in CRC compared with normal tissues (Supplementary Table 4). A total of 77 intersecting mRNAs were generated (Supplementary Table 5). We utilized Venn diagrams to illustrate the intersecting states of miRNA and mRNA (Figure 4A and B) and took the 20 miRNAs and 20 mRNAs expressing the most significant upregulation and downregulation, respectively, to draw heat maps (Figure 4C and D).

Figure 2 Heat map of circRNAs differentially expressed in GSE126095. By analyzing the circRNA profile, 55 differentially expressed circRNAs were obtained. Red indicates high expression, while blue indicates low expression. A deeper color indicates a more significant difference. The cut-off value was a |log₂ fold change| > 2 and an adjusted P < 0.05.

Figure 3 Structural pattern of 49 differentially expressed circRNAs obtained from Cancer-Specific CircRNA Database. Red indicates the miRNA response element, blue indicates the RNA binding protein, and green indicates the open reading frame. CSCD: Cancer-Specific CircRNA Database (http://gb.whu.edu.cn/CSCD/).

Figure 4 Venn diagrams and heat maps of miRNAs and mRNAs. A: Twenty-five miRNAs were identified from the intersection of 114 differentially expressed (DE) miRNAs with 1602 circRNA targets predicted by the Cancer-Specific CircRNA Database; B: 77 mRNAs were obtained from the intersection of 267 DEmRNAs and 7190 miRNA targets predicted using the TargetScan and miRDB databases; The expression level of each 20 DEmiRNAs (C) and 20 mRNAs (D) expressing the most significant upregulation and downregulation between control tissues (type C) and tumor tissues (type T), respectively. |log₂ fold change (FC)| > 1 and an adjusted P < 0.05 were considered the statistical criteria for DEmiRNA; adjusted P < 0.05 and |log₂ FC| > 2 were considered the statistical criteria for DEmRNA.

Table 1 Basic characteristics of 4 microarray datasets in Gene Expression Omnibus and The Cancer Genome Atlas databases.

Data source	Platform	Series	Sample size
			Tumor	Control
circRNA	GPL19978	GSE126095	10	10
miRNA	GPL11487	GSE41655	33	15
mRNA	GPL6480	GSE41657	25	12
mRNA	TCGA	None	488	42

TCGA: The Cancer Genome Atlas.

Construction of the CRC-specific ceRNA network

Based on the recognized 49 circRNAs, 25 miRNAs and 77 mRNAs, we removed unconnected nodes and chose 19 circRNAs (hsa_circ_0000520, hsa_circ_0000519, hsa_circ_0001955, hsa_circ_0028198, hsa_circ_0080210, hsa_circ_0007158, hsa_circ_0000375, hsa_circ_0000026, hsa_circ_0023685, hsa_circ_0000370, hsa_circ_0061817, hsa_circ_0005927, hsa_circ_0072088, hsa_circ_0018909, hsa_circ_0013912, hsa_circ_0071681, hsa_circ_0011536, hsa_circ_0043278, and hsa_circ_0006220), 13 miRNAs (hsa-miR-423-5p, hsa-miR-532-3p, hsa-miR-765, hsa-miR-1224-5p, hsa-miR-650, hsa-miR-769-5p, hsa-miR-671-5p, hsa-miR-1290, hsa-miR-125a-3p, hsa-miR-601, hsa-miR-198, hsa-miR-1202, and hsa-miR-1182) and 28 mRNAs (RNF43, DSG3, AZGP1, SST, DES, TCF21, MFAP4, EREG, BCAS1, C1QA, SPARCL1, CXCL3, EPHB3, TRAF3IP3, TRPM2, CA2, LRRC19, SCG2, C16orf89, ADH1A, MZB1, HAPLN1, S100A2, GPR34, MS4A12, ITLN1, DHRS9, and CHGB) in to construct a circRNA-miRNA-mRNA regulatory network utilizing Cytoscape 3.7.2 (Figure 5). The expression levels of these RNA molecules are shown in Figure 6.

Figure 5 CircRNA-miRNA-mRNA regulatory network of colorectal cancer. CircRNAs, miRNAs and mRNAs are represented by red diamonds, green triangles and blue ellipses, respectively. Dark colors represent overexpression, and light colors represent low expression.

Figure 6 The expression levels of circRNAs, miRNAs and mRNAs in the competing endogenous RNA (ceRNA) network. A and B: Expression levels of 19 circRNAs in the ceRNA network; C and D: Expression levels of 13 miRNAs in the ceRNA network; E and F: Expression levels of 25 mRNAs in the ceRNA network. ^aP < 0.001; ^bP < 0.01; ^cP < 0.05.

GO and KEGG functional enrichment analysis

GO analysis revealed that the enrichments of the identified mRNAs in the ceRNA network were mainly in the ‘retinoic acid metabolic process’, ‘retinol metabolic process’, ‘digestive tract morphogenesis’, and ‘leukocyte chemotaxis’ (biological processes; P < 0.005) (Figure 7A); ‘collagen-containing extracellular matrix’, ‘endoplasmic reticulum lumen’ and ‘endoplasmic reticulum chaperone complex’ (cellular components; P < 0.02) (Figure 7B); and ‘alcohol dehydrogenase (ADH) [NAD(P)⁺] activity’, ‘retinol dehydrogenase activity’ and ‘oxidoreductase activity’ (molecular functions; P < 0.02) (Figure 7C). The results of KEGG pathway analysis involving the ceRNA network indicated that the target genes were mainly enriched in ‘retinol metabolism’, ‘gastric acid secretion’, ‘nitrogen metabolism’, ‘NOD-like receptor signaling pathway’, ‘proximal tubule bicarbonate reclamation’, ‘collecting duct acid secretion’ and ‘tyrosine metabolism’ (P < 0.05) (Figure 7D). According to these results, the retinol metabolic process, leukocyte chemotaxis, extracellular matrix remodeling, endoplasmic reticulum stress, ADH activity, gastric acid secretion, nitrogen metabolism and NOD-like receptor signaling pathway might participate in the tumorigenesis of CRC.

Figure 7 GO and KEGG pathway analyses of the identified mRNAs in the competing endogenous RNA network. A: GO analysis of biological processes; B: GO analysis of cellular components; C: GO analysis of molecular functions; D: KEGG pathway analysis. P < 0.05 were considered significant.

Survival analysis and construction of the prognostic ceRNA subnetwork

We obtained mRNA profiles and clinical information of CRC patients in the TCGA database and performed survival analysis for each mRNA in the obtained ceRNA network. Finally, we recognized that 3 mRNAs (CA2, ITLN1, and LRRC19) were significantly correlated with the clinical outcome of CRC patients (Figure 8). The patients with the upregulation of CA2 (P = 0.002), ITLN1 (P = 0.001) and LRRC19 (P = 0.032) had a better prognosis than the corresponding group with low expression. Considering the 3 mRNAs identified above, we successfully constructed and visualized a ceRNA subnetwork including 5 circRNAs (hsa_circ_0080210, hsa_circ_0007158, hsa_circ_0000375, hsa_circ_0018909 and hsa_circ_0011536) and 3 miRNAs (hsa-miR-601, hsa-miR-671-5p and hsa-miR-765) (Figure 9).

Figure 8 Survival analysis of the identified mRNAs in the circRNA-miRNA-mRNA regulatory network using The Cancer Genome Atlas database . A: Comparison between patients with upregulation or downregulation of CA2 (P = 0.002); B: Comparison between patients with upregulation or downregulation of ITLN1 (P = 0.001); C: Comparison between patients with upregulation or downregulation of LRRC19 (P = 0.032).

Figure 9 The competing endogenous RNA subnetwork associated with colorectal cancer prognosis. Based on mRNAs with potential for the clinical outcome prediction of colorectal cancer, a prognostic competing endogenous RNA subnetwork was successfully developed. Diamonds, triangles and ellipses represent circRNAs, miRNAs and mRNAs, respectively.

DISCUSSION

We identified the intersection of CRC-specific DEcircRNAs, DEmiRNAs and DEmRNAs in the GEO database with the target molecules of circRNAs and miRNAs predicted by relevant databases. Then, 19 circRNAs, 13 miRNAs and 28 mRNAs were identified to develop a circRNA-miRNA-mRNA regulatory network that may play a pivotal role in the progression of CRC. Subsequently, we conducted GO and KEGG pathway analyses for the 28 mRNAs to expand our understanding of the vital pathophysiological process of CRC initiation and progression. Finally, we verified the differential expression of the identified mRNAs in the TCGA database, screened prognosis-related mRNAs by conducting a survival analysis, and constructed a prognostic subnetwork using 5 circRNAs, 3 miRNAs and 3 mRNAs.

CircRNAs, as stable, abundant and conserved ceRNAs that act as miRNA sponges, could be identified as valuable indicators for the diagnosis and pathogenesis of CRC. Similar studies also supported the evidence that DEcircRNAs identified in our study could be vital components in the ceRNA network, which modulate crucial gene expression in the initiation and progression of cancer, especially CRC. For example, the dysregulation of hsa_circ_0000520 could affect the tumorigenesis of cervical cancer and breast cancer through the miRNA-mRNA axis[9,10]. The hsa_circ_0072088 in the circRNA-miRNA-mRNA network was identified as being related to CRC and lung cancer progression[11,12]. The hsa_circ_0001955 was found to mediate a ceRNA network in CRC by bioinformatics analysis and experimental validation[13]. The hsa_circ_0005927 was verified in gastric cancer and could be a biomarker for gastric cancer screening[14]. The hsa_circ_0000026 was found to be expressed at low levels in gastric cancer and may correlate with the progression of CRC[15]. The hsa_circ_0000370 in plasma showed diagnostic value for CRC and might be involved in tumorigenesis[16]. All of these studies indicate that circRNAs may participate in CRC progression and could be vital biomarkers for diagnosis as well as therapeutic targets. The majority of 55 DEcircRNAs identified in our study were innovative biomarkers in CRC and still require further investigation in the future.

To further explore signaling pathways that might play an important role in the tumorigenesis and progression of CRC, we conducted GO and KEGG analyses of 28 identified genes in the ceRNA network. Growing studies have confirmed that the signaling pathways uncovered in our study participate in crucial pathological processes in many kinds of cancer. Bhattacharya et al[17] reported that the inhibition of retinoic acid signaling, a key regulator of intestinal immunity, could promote the tumorigenesis of CRC by cytotoxic T cells, and retinoic acid catabolizing enzyme was a promising negative predictor for the prognosis of CRC patients. Retinol dehydrogenase 16, one of the isoforms of the rate-limiting enzyme of the retinol cycle, was reported to increase the level of retinoic acid, and associate with the tumor size of hepatocellular carcinoma and poor overall survival of patients as well[18]. Mesenchymal cells in the intestine, called cancer-associated fibroblasts, exert critical functions to regulate a variety of activities, including intestinal inflammation, epithelial proliferation, extracellular matrix remodeling and metastasis, which could affect the microenvironment and promote CRC development and progression[19]. A recent study revealed that leukocyte chemotaxis and adhesion were distinctly reduced in the vasculature of CRC[20], and some miRNAs, including miR-15A and miR-16-1, were able to modulate the pathway of immune regulatory B cell chemotaxis in CRC, which could affect the tumor growth and survival time[21]. Endoplasmic reticulum stress is one of the pivotal processes in carcinogenesis and could represent an innovative therapeutic target in resistant tumors[22]. Cheng et al[23] found that endoplasmic reticulum stress participated in the apoptosis and autophagy in CRC induced by apatinib, a novel tyrosine kinase inhibitor. The consumption of alcohol was reported to increase the risk of colorectal adenomas, and ADH might modify the correlation between alcohol consumption and colorectal adenomas[24]. ADH expression has the potential to be a prognostic marker of pancreatic adenocarcinoma[25], and acetaldehyde is recognized to elevate the possibility of chemically induced rectal carcinogenesis[26]. A recent meta-analysis revealed that gastric acid suppressant use showed a significant correlation with poor survival for patients receiving oral chemotherapy for gastrointestinal tract cancer, supporting a possible negative impact of gastric acid suppressants on the survival outcome of CRC[27]. As one of the most fundamental requirements for biosynthesis, nitrogen metabolism is utilized and modulated to sustain the increased demand for nitrogen sources in cancer proliferation[28]. The regulation of nitrogen metabolism participates in the process of obesity-associated pancreatic cancer, small cell lung cancer and CRC metastasis[29-31]. The NOD-like receptor signaling pathway was one of the enriched pathways of genes with aberrant expression among pancreatic, thyroid, and renal cancer compared with healthy controls through bioinformatics analyses[32-34]. Furthermore, KEGG website (http://www.kegg.jp/) was utilized to analyze the specific steps of these biological processes affected by identified mRNAs in ceRNA network. It is well-known that CA2 participates in the processes of combining water and carbon dioxide to generate carbonic acid, acting as one of the key enzymes in the proximal tubule bicarbonate reclamation and collecting duct acid secretion. According to the pathway diagrams in KEGG website, CA2 also involves in gastric acid secretion and arginine biosynthesis in nitrogen metabolism. CXCXL3 is one of the downstream chemokines of NOD-nuclear factor-kappa B pathway in NOD-like receptor signaling pathway. ADH, as one of the key enzymes for the mutual transformation of all-trans-retinal and all-trans-retinol (vitamin A), participates in the final metabolic process of dopamine to 3-methoxy-4-hydroxy-phenylethylene-glycol in tyrosine metabolism. In summary, the retinol metabolic process, leukocyte chemotaxis, extracellular matrix remodeling, endoplasmic reticulum stress, ADH activity, gastric acid secretion, nitrogen metabolism and NOD-like receptor signaling pathway might represent essential signaling pathways involved in the pathogenesis of CRC. However, the molecular mechanism of tumorigenesis and progression is quite complicated and still requires further exploration.

We verified the differential expression of the mRNAs in the ceRNA network using the TCGA database, combined the results of the survival analysis, screened prognosis-related mRNAs, and finally used 5 circRNAs, 3 miRNAs and 3 mRNAs to construct a prognostic subnetwork. The roles of the 5 circRNAs (hsa_circ_0080210, hsa_circ_0007158, hsa_circ_0000375, hsa_circ_0018909 and hsa_circ_0011536) involved in tumorigenesis require investigation, which implies that these circRNAs might have the potential to become novel indicators for CRC diagnosis and targeted treatment. Some studies were consistent with our result that 3 miRNAs (hsa-miR-601, hsa-miR-671-5p and hsa-miR-765) might represent promising biomarkers of cancer progression and prognosis. MiR-601 was identified to hold diagnostic value for CRC with 69.2% sensitivity and 72.4% specificity for CRC diagnosis[35], and might suppress the proliferation and invasion of esophageal squamous cell carcinoma and breast cancer[36,37]. Some studies revealed that miR-671-5p could promote and maintain the oncogenesis and progression of various cancers, such as esophageal cancer, breast cancer, glioblastoma and melanoma[38-41], and act as a prognostic predictor of locally advanced rectal cancer due to the significant upregulation in pathological response to neoadjuvant chemoradiotherapy[42]. The expression level of hsa-miR-765 could be utilized to independently predict overall survival and disease-free survival and correlate with the tumor stage, clinical stage and lymph node metastasis in esophageal squamous cell carcinoma[43]. Hsa-miR-765 could also promote the aggressiveness of hepatocellular carcinoma and osteosarcoma[44,45]. To date, although some studies have focused attention on the roles of these 3 miRNAs in many types of cancer, hsa-miR-671-5p and hsa-miR-765 have not been recognized as promising biomarkers of CRC before our study.

We found that the upregulation of CA2, ITLN1 and LRRC19 might be related to better clinical outcomes in CRC patients. CA2, the gene that encodes carbonic anhydrase Ⅱ, was validated to be downregulated in CRC tissue and cell lines compared with healthy controls through experimental assays, and the overexpression of CA2 suppressed tumor cell growth in vitro and in vivo and elevated the sensitivity of CRC cells to chemotherapy drugs[46]. CA2 might also be useful for the survival prediction of CRC patients[47]. Other carbonic anhydrase isoforms, such as CA1, CA4, CA9 and CA12, could exert effects on favorable outcomes or rescue the tumor progression of CRC[48-50]. In summary, CA2 could be a beneficial predictor for CRC diagnosis and prognosis. The expression of ITLN1, which encodes intelectin-1 (also known as omentin-1), presented a sequentially descending trend with the mucosa-adenoma-carcinoma process[51]. The downregulation of ITLN1 was reported to be related to poor prognosis among CRC patients at advanced stages[52]. ITLN1 could also suppress tumorigenesis and correlate with a better prognosis probability of neuroblastoma and ovarian cancer[53,54]. The low expression of LRRC19 (leucine-rich repeat containing 19) was identified as an independent risk factor and prognostic biomarker of kidney renal clear cell carcinoma and was involved in selenium adjuvant therapy[55]. In addition, LRRC19 could have predictive power for sensitivity to AZD0530 in pancreatic tumor and promote the elimination of pathogenic bacteria from the kidneys[56,57]. However, whether LRRC19 could be regarded as an indicator of oncogenesis and the clinical outcome of CRC has not been illustrated before. In summary, CA2, ITLN1 and LRRC19 might have the potential to become novel biomarkers for CRC diagnosis, especially early screening, as well as patient prognosis. At present, the hsa_circ_0011536/hsa-miR-601/CA2 axis and the complex interaction among circRNAs (hsa_circ_0080210, hsa_circ_0007158, hsa_circ_0000375 and hsa_circ_0018909), miRNAs (hsa-miR-671-5p and hsa-miR-765) and mRNAs (ITLN1 and LRRC19) could provide novel insights into the pathogenesis of and therapeutic options for CRC and need further investigation.

Collectively, we have provided a deeper understanding of the circRNA-related ceRNA mechanism of CRC by developing a circRNA-miRNA-mRNA regulatory network. Through GO and KEGG functional enrichment analysis, the retinol metabolic process, leukocyte chemotaxis, extracellular matrix remodeling, endoplasmic reticulum stress, ADH activity, gastric acid secretion, nitrogen metabolism and the NOD-like receptor signaling pathway might play critical roles in the initiation and progression of CRC. After being combined with survival analysis, a prognostic subnetwork was constructed, including 5 circRNAs, 3 miRNAs and 3 mRNAs, which could be novel candidate biomarkers for the clinical outcome of CRC. Our research still has some shortcomings. First, molecular-level verification in clinical samples and CRC cell lines should be applied to validate biomarkers and clarify the actual significance of the regulatory ceRNA network and prognostic subnetwork. In addition, the lack of research on the downstream target molecules of the ceRNA network makes it difficult to completely elucidate the specific mechanism in the occurrence and development of CRC. Moreover, due to the sample size of the datasets, the sample size of this study is not too large.

CONCLUSION

In summary, we constructed a CRC-specific circRNA-miRNA-mRNA regulatory network and performed functional enrichment analysis, which may assist in revealing the mechanism of carcinogenesis. A prognostic ceRNA subnetwork was successfully developed, thereby identifying several RNA molecules that could serve as innovative and noninvasive indicators for the RNA-based early screening and prognostic prediction of CRC. Comprehensive experimental studies are still required to enrich the understanding of ceRNAs, which is essential to illustrate the pathogenesis and prognosis of CRC and provide new opportunities for targeted therapeutics.

ARTICLE HIGHLIGHTS

Research background

Colorectal cancer (CRC) is one of the most common malignant tumors worldwide. Circular RNA (circRNA) and microRNA (miRNA), acting as competing endogenous RNAs (ceRNAs), have been investigated to play vital roles in carcinogenesis. Accumulating evidence highlights that it is necessary to further explore the biological function of circRNA and miRNA in CRC pathogenesis and prognosis.

Research motivation

Dysregulated circRNAs, miRNAs and mRNAs could closely associate with the progression and prognosis of CRC and act as potential CRC-specific predictors, but the competitive regulatory pattern and biological function mechanism among circRNAs, miRNAs and mRNAs are still complicated and have not yet been elucidated.

Research objectives

This study aimed to uncover a CRC-specific competitive regulatory model among circRNAs, miRNAs and mRNAs and explore the subnetwork associated with CRC prognosis.

Research methods

Expression profiles of circRNAs, miRNAs and mRNAs were downloaded from the Gene Expression Omnibus (GEO) database. Differentially expressed (DE) circRNAs, miRNAs and mRNAs in CRC tissues and the predicted target molecules of circRNAs and miRNAs were intersected to obtain a CRC-specific ceRNA network. GO and KEGG pathway analyses were conducted to explore the mechanism of CRC pathogenesis. Based on the survival analysis using the gene profiles and clinical information in The Cancer Genome Atlas (TCGA) database, the mRNAs significantly associated with the clinical outcome of CRC patients were identified and a prognostic subnetwork was constructed.

Research results

There were 55 DEcircRNAs, 114 DEmiRNAs and 267 DEmRNAs of CRC in three datasets (GSE126095, GSE41655 and GSE41657) from GEO database. After intersected with predicted targets, 19 circRNAs, 13 miRNAs and 28 mRNAs were chosen to develop ceRNA network. Go and KEGG analyses indicated that several signaling pathways might participate in the tumorigenesis. After verifying effect in TCGA database by survival analysis, we finally recognized 3 mRNAs (CA2, ITLN1 and LRRC19) associated with prognosis, and constructed a ceRNA subnetwork including 5 circRNAs (hsa_circ_0080210, hsa_circ_0007158, hsa_circ_0000375, hsa_circ_0018909 and hsa_circ_0011536) and 3 miRNAs (hsa-miR-601, hsa-miR-671-5p and hsa-miR-765).

Research conclusions

A circRNA-miRNA-mRNA regulatory network closely associated with the progression and clinical outcome of CRC was identified, which might include promising biomarkers for carcinogenesis and therapeutic targets.

Research perspectives

We have provided a deeper understanding of the circRNA-related ceRNA mechanism of CRC by developing a circRNA-miRNA-mRNA regulatory network. Comprehensive experimental studies are still required to confirm our findings and provide new opportunities for targeted therapeutics.

ACKNOWLEDGEMENTS

We thank The Gene Expression Omnibus and The Cancer Genome Atlas project for providing invaluable datasets for statistical analyses.