Whole-transcriptome analysis reveals a potential hsa_circ_0001955/hsa_circ_0000977-mediated miRNA-mRNA regulatory sub-network in colorectal cancer

Background: Circular RNAs (circRNAs), a novel class of non-coding RNAs, have been found to act as microRNA (miRNA) sponges and thus play key roles in biological processes and pathogenesis. However, studies regarding circRNAs in colorectal cancer (CRC) remain inadequate. Results: By differential expression analysis, 10 candidate circRNAs (6 upregulated and 4 downregulated circRNAs) were chosen. 9 of 10 circRNAs were available on CSCD and their structure showed the binding potential of miRNA. Intersection analysis revealed that miR-145-5p, miR-3127-5p, miR-761, miR-4766-3p, miR-135a-5p, miR-135b-5p, miR-374a-3p and miR-330-3p were 8 miRNAs with the most potential in binding circRNAs. Further expression validation and correlation analysis demonstrated hsa_circ_0001955/miR-145-5p and hsa_circ_0000977/miR-135b-5p axes as key pathways in CRC. Subsequently, target gene prediction, differential expression analysis, intersection analysis and correlation analysis showed that CDK6, MMP12 and RAB3IP were the three potential downstream targets of hsa_circ_0001955/miR-145-5p axis and FOXO1, MBNL1, MEF2C, RECK, PPM1E, TTLL7 and PCP4L1 were the seven potential downstream targets of hsa_circ_0000977/miR-135b-5p axis in CRC. Finally, we also confirmed that expression of hsa_circ_0001955 or hsa_circ_0000977 was significantly positively correlated with their individual targets in CRC. Conclusions: In the present work, we constructed a potential hsa_circ_0001955/hsa_circ_0000977-mediated circRNA-miRNA-mRNA regulatory network in CRC by a series of in silico analysis and experimental validation. Methods: Whole-transcriptome microarrays from CRC and matched normal samples were obtained from GEO. The structure of circRNA was identified by CSCD. starBase and miRNet were successively used to predict miRNA of circRNA and target gene of miRNA. Expression correlation between RNA-RNA interactions was assessed using GEO and TCGA data. Finally, a potential circRNA-miRNA-mRNA network was established based on competing endogenous RNA (ceRNA) hypothesis.


INTRODUCTION
Colorectal cancer (CRC) ranks as the third most common malignancy and is the third leading cause of cancer-associated deaths in the United States [1]. A lot of factors, including diet, polyps and chronic inflammation, have been documented to link to CRC pathogenesis, consisting of multiple processes. However, the alterations and roles of genomics and epigenetics in CRC initiation and progression are still rarely known. Moreover, despite received standard therapies, like surgery combined with chemotherapy and radiotherapy, many CRC patients will eventually recur or metastasize [2]. Therefore, it is urgent to further investigate molecular mechanism of CRC carcinogenesis, and seek and develop effective therapeutic targets for patients with CRC. Circular RNAs (circRNAs) are a class of novel endogenous non-coding RNAs (ncRNAs) with a covalently closed continuous loop, and they are generated by the selective cleavage of premature RNAs [3]. The loop structure produced by linking 5' and 3' end renders circRNAs resistant to exonucleases and more stable than their corresponding linear transcripts [4]. In 2011, Salmena L et al. proposed a competing endogenous RNA (ceRNA) hypothesis, in which ncRNAs could positively regulate mRNA expression through competitively binding to shared microRNAs (miRNAs) [5]. Recent studies have well-demonstrated that circRNAs increase oncogene or tumor suppressor expression by ceRNA mechanism, and thereby participate in cancer development. For example, Xue D et al. found that circ-AKT3 inhibited clear cell renal cell carcinoma metastasis by modulating miR-296-3p/Ecadherin axis [6]; Han D et al. suggested that circMTO1 acted as the sponge of miR-9 to suppress progression of hepatocellular carcinoma [7]; Wang C et al. showed that circFOXO3 induced cell apoptosis in urothelial carcinoma by interaction with miR-191-5p [8]. circRNAs have also been confirmed to play critical roles in various biological behaviors of CRC, including proliferation, migration and invasion [9,10]. Moreover, circRNAs can serve as diagnostic or prognostic biomarkers for CRC patients [11,12].
Establishment of circRNA-miRNA-mRNA regulatory networks may provide key clues for us to explore the detailed molecular mechanism of human disorders including malignancies. For example, Liu K et al. uncovered the possible mechanism of Hsp90 inhibitorinduced cell death of CRC by identification of a circRNA-miRNA-mRNA regulatory network [13]. Yuan W et al. identified circRNAs as ceRNAs for miRNA-mRNA in CRC [14]; Song W et al. also constructed a circRNA-associated ceRNA network in CRC [15]. The two studies constructed circRNA-miRNA-mRNA network using RNA data from different CRC patients.
In this study, we obtained and analyzed the circRNA-miRNA-mRNA ceRNA regulatory network using 10 CRC patients' circRNA, miRNA and mRNA expression profile data. Furthermore, TCGA data were employed to validate the analytic results. A hsa_circ_ 0001955/hsa_circ_0000977-mediated miRNA-mRNA regulatory sub-network in colorectal cancer was successfully established. Finally, experimental validation was given. Those findings from the present work may provide novel insight into the molecular mechanism of CRC carcinogenesis.

Selection of 9 potential circRNAs in CRC
To find the potential circRNAs in CRC, circRNA dataset of 10 CRC patients (GSE126094) from GEO database was selected. GEO2R tool provided by GEO database was used to perform differential expression analysis for GSE126094. As listed in Figure 1A and Supplementary Table 1, a total of 1850 significant  DECs, containing 1831 upregulated DECs and 19 downregulated DECs, were screened out. Next, we further selected some more potential DECs from the 1850 significant DECs based on the selection criterion of |log 2 FC| > 4, and 10 circRNAs, including 6 upregulated and 4 downregulated circRNAs, were finally identified ( Figure 1B and Table 1). The circBase names and parental genes of the 10 circRNAs were presented in Table 2. Subsequently, the structure of 9 circRNAs, hsa_circ_0072088 ( Figure 2H) and hsa_circ_0043278 ( Figure 2I) (hsa_circ_0000981 not available) was described in Figure 2 using the data from CSCD, indicating that all the 9 circRNAs had microRNA response elements (MREs). Taken together, the 9 circRNAs may be the key circRNAs in CRC by acting as miRNA sponges.

Identification of downstream target genes of miRNAs in CRC
It has been widely acknowledged that miRNAs exert their biological functions by negatively regulated expression of downstream target genes. Therefore, the next step was to ascertain the downstream target genes of hsa_circ_0001955/miR-145-5p and hsa_circ_ 0000977/miR-135b-5p axes. Target genes of miR-145-5p and miR-135b-5p were predicted through miRNet database, and 238 and 83 target genes were forecasted to potentially bind to miR-145-5p and miR-135b-5p, respectively (Supplementary Table 4). DEGs between  CRC and normal controls were obtained by GEO2R analysis using GSE126092 from the same 10 CRC patients with GSE126093 and GSE126094. As presented in Figure 6A and Supplementary Table 5, a total of 2787 DEGs, including 1194 upregulated DEGs and 1593 downregulated DEGs, were finally found. According to action mechanism of miRNA, miRNA expression should be negatively associated with expression of target gene. Thus, we performed intersection analysis for target genes set and DEGs set. As shown in Figure 6B, 19 genes (CBFB, CD44, CDK4, CDK6, CFTR, MMP1, MMP12, MYC, PIGF, POU5F1, ROBO2, SNTB1, TGFBI, VEGFA, PSAT1, LMNB2, RAB3IP, C11orf65 and CRNDE) were commonly appeared in miR-145-5p's target genes set and upregulated DEGs set. And 10 genes (APC, FOXO1, MBNL1, MEF2C, RECK, KLF4, PPM1E, ARC, TTLL7 and PCP4L1) were commonly appeared in miR-135b-5p's target genes set and downregulated DEGs set ( Figure 6C). Next, expression levels of the 19 upregulated genes and 10 downregulated genes were further validated using TCGA data by starBase database. The results demonstrated that 17 of 19 upregulated genes (except CFTR and POU5F1) and 10 of 10 downregulated genes fitted the previous analytic findings as shown in Figure 6D and 6E, respectively. The 27 genes were considered as the downstream target genes of miR-145-5p and miR-135b-5p in CRC and selected for subsequent analysis.

DISCUSSION
Increasing evidences have demonstrated that circRNAs serve as crucial players in initiation and progression of human cancers, including CRC. In spite of some studies have reported their constructed circRNA-miRNA-mRNA regulatory network in CRC [14,15]. However, the current knowledge of circRNA-associated ceRNA network in CRC remain inadequate and need to be further investigated. In the present study, by a series of in silico analysis, we identified a potential circRNA-miRNA-mRNA ceRNA regulatory network in CRC pathogenesis.
After performing differential expression analysis, 10 circRNAs were selected for further research. Most of the 10 circRNAs have been found to be closely linked to cancer development and progression. For example, hsa_circRNA_103809 regulates CRC proliferation and migration via miR-532-3p/FOXO4 axis [16]; hsa_circRNA_101555 is involved in carcinogenesis of breast cancer [17] hsa_circRNA_102619 played a suppressive role in progression of pancreatic ductal adenocarcinoma [18]. By using circBase and CSCD databases, 9 of 10 circRNAs were finally identified.
Multiple studies have confirmed that circRNAs positively regulate expression of downstream genes by acting as miRNA sponges [19,20]. Therefore, potential miRNAs of the 9 circRNAs were predicted. Followed by a series of analyses including miRNA expression analysis and correlation analysis between circRNA and miRNA in CRC, two pathways (hsa_circ_0001955/miR-145-5p and hsa_circ_0000977/miR-135b-5p) were regarded as the key axes in carcinogenesis of CRC. The functions of miR-145-5p and miR-135b-5p have been investigated, which were in accordance with ceRNA hypothesis. For example, miR-145-5p inhibited growth, migration and invasion of CRC [21,22]. Silencing miR-135b-5p sensitized CRC cells to oxaliplatin-induced apoptosis by upregulating FOXO1 expression [23].
As is known to all, miRNAs exert their functions by negatively modulating gene expression [24][25][26][27]. Firstly, target genes of miR-145-5p and miR-135b-5p were predicted by a comprehensive database, namely miRNet. Then, the DEGs between CRC tissues and normal tissues were screened. 29 targets were selected by intersection of target gene set and DEG set, after which expression of the 29 targets were validated using TCGA data by starBase. Subsequently, correlation analysis of target gene with miRNA was performed through both GEO and TCGA data. Finally, we found that three targets of miR-145-5p (CDK4, MMP12 and RAB3IP) and seven targets of miR-135b-5p (FOXO1, MBNL1, MEF2C, RECK, PPM1E, TTLL7 and PCP4L1) possessed the most potential. These targets have been well-documented to closely link to cancer progression, including CRC. For instance, miR-143-3psuppressed CDK4 enhanced CRC growth [28].
According to ceRNA mechanism, there should be positive expression relationship between circRNA and target gene. Intriguingly, hsa_circ_0001955 expression was positively correlated with expression of CDK4, MMP12 and RAB3IP, and hsa_circ_0000977 was also positively associated with expression of FOXO1, MBNL1, MEF2C, RECK, PPM1E, TTLL7 and PCP4L1 in CRC.
Taken all these data together, we successfully established a circRNA-miRNA-mRNA triple ceRNA sub-network mediated by two circRNAs (hsa_circ _0001955 and hsa_circ_0000977) in CRC. Moreover, the expression levels of RNAs in the network and expression correlation were followingly experimentally validated using clinical CRC samples, partially supporting the analytic accuracy of in silico analysis.
The current work provides novel insight into the detailed molecular mechanism of CRC carcinogenesis,  AGING and the research route may also be applied to explore the ceRNA molecular interaction mechanism of circRNA-miRNA-mRNA network in initiation and progression of other types of human cancer. However, more efforts should be given to elucidate the function of the established network in CRC with in vitro and in vivo experiments. Meanwhile, the diagnostic and prognostic values of each RNA in the established network should be evaluated using a large number of clinical CRC data in the future, which will help to develop promising diagnostic and prognostic biomarkers for patients with CRC.

CONCLUSIONS
In summary, this study reveals a potential hsa_circ_ 0001955/hsa_circ_0000977-mediated circRNA-miRNA-mRNA ceRNA sub-network in CRC by wholetranscriptome analysis, including differential expression analysis, intersection analysis and correlation analysis. Despite need more experimental and clinical validations, targeting the molecules in this sub-network may represent a promising treatment for patients with CRC.

Inclusion of datasets
In this study, we intended to construct a key circRNA-miRNA-mRNA regulatory network in CRC. Firstly, we selected the datasets for analysis using NCBI GEO database (http://www.ncbi.nlm.nih.gov/geo/). Selection criteria: only datasets that focused on human CRC tissue samples were included; datasets regarding CRC cell lines or animal CRC tissue samples were excluded; only datasets with sample count more than ten were included; only datasets containing circRNA expression, miRNA expression and mRNA expression were included. Finally, only dataset GSE126095 met all these criteria mentioned above. GSE126095 had three sub-datasets (GSE126094, GSE126093 and GSE126092) from 10 CRC tissues and 10 matched normal tissues. GSE126094, based on the platform of GPL19978 Agilent-069978 Arraystar Human CircRNA microarray V1, investigated circRNA expression profile; GSE126093, based on the platform of GPL18058 Exiqon miRCURY LNA microRNA array G7, investigated miRNA expression profile; and GSE126092, based on the platform of GPL21047 Agilent-074348 Human LncRNA v6 4X180K.
circBase and CSCD (Cancer-Specific CircRNA Database) analysis circBase, a database for circRNA-associated studies, provides scripts to acquire known and novel circRNAs in sequencing data [30]. circBase was employed to identify parental genes of potential circRNAs. CSCD (http://gb.whu.edu.cn/CSCD/), which is an online tool to investigate cancer-specific circRNAs, was utilized to obtain structure of potential circRNAs.

miRNA prediction
miRNAs that might bind to circRNAs were predicted by starBase database (http://starbase.sysu.edu.cn/). starBase database is an open-source platform that provides the ncRNA interactions from CLIP-seq, degradome-seq and RNA-RNA interactome data [31,32]. The circBase ID of circRNA was entered into starBase database, after which potential miRNAs of circRNA were directly presented on the webpage.

Validation of miRNA expression
The miRNAs that commonly appeared in target miRNAs set and DEmiRNAs set were regarded as potential miRNAs of circRNAs in CRC. Expression of these miRNAs was further validated using TCGA data by starBase database [31,32]. P < 0.05 was considered as statistically significant.
Target gene prediction miRNet (https://www.mirnet.ca/), an integrated platform linking miRNAs, targets and functions, was introduced to perform miRNA's target gene prediction [33]. The miRNA-gene interactions were directly downloaded from the webpage.

Validation of gene expression
The genes that commonly appeared in target genes set and DEGs set were considered as potential target genes of miRNAs in CRC. Expression of these target genes were further confirmed using TCGA data using starBase database [31,32]. P < 0.05 was considered as statistically significant.

Correlation analysis for RNA-RNA interactions
The expression correlation of hsa_circ_0001955 or hsa_circ_0000977 with miRNA in CRC was evaluated using the data from GSE126094 and GSE126093; the expression correlation of miR-145-5p or miR-135b-5p with target gene in CRC was assessed using the data from GSE126093 and GSE126092; the expression correlation of hsa_circ_0001955 or hsa_circ_0000977 with target gene in CRC was determined using the data from GSE126094 and GSE126092. Then, starbase database was also employed to validate the expression correlation of miR-145-5p or miR-135b-5p with target gene. |R| > 0.1 and P-value < 0.05 were set as the thresholds for identifying significant RNA-RNA interactions.

Clinical tissues and qRT-PCR
Fresh frozen CRC cancer tissues and adjacent normal tissues were obtained from 20 CRC patients who received surgical treatment at the First Affiliated Hospital of Zhejiang University (Hangzhou, China), and were frozen and stored in liquid nitrogen. This study was approved by the Ethical Committee of First Affiliated Hospital of Zhejiang University. Each informed consent was acquired from every patient. The expression of potential circRNAs, miRNAs and genes in these clinical samples was detected by the method of qRT-PCR as previously described [26]. The Primers used in this study were designed and purchased from RiboBio Co. Ltd (Guangzhou, China).

Statistical analysis
Most of statistical analyses have been calculated by the online databases or tools. Expression correlation for RNA-RNA interactions was perform by Person correlation coefficient using Graphpad Prism software (version 7). Paired Student's t-test was used to evaluate expression differences of circRNA, miRNA or gene between normal group and cancer group. Results with P-value < 0.05 were considered as statistically significant.

Ethics approval
This study was approved by the Ethical Committee of First Affiliated Hospital of Zhejiang University. Each informed consent was acquired from every patient.

AUTHOR CONTRIBUTIONS
All authors contributed to the conception and design of this manuscript, collection and assembly of data, manuscript writing, and final approval.