Non-coding effects of circular RNA CCDC66 promote colon cancer growth and metastasis

Circular RNA (circRNA) is a class of non-coding RNA whose functions remain mostly unknown. Recent studies indicate circRNA may be involved in disease pathogenesis, but direct evidence is scarce. Here we characterize the functional role of a novel circRNA, circCCDC66, in colorectal cancer (CRC). RNA-Seq data from matched normal and tumor colon tissue samples identified numerous circRNAs specifically elevated in cancer cells, several of which were verified by quantitative RT-PCR. CircCCDC66 expression was elevated in polyps and colon cancer and was associated with poor prognosis. Gain-of-function and loss-of-function studies in CRC cell-lines demonstrated that circCCDC66 controlled multiple pathological processes, including cell proliferation, migration, invasion, and anchorage-independent growth. In-depth characterization revealed that circCCDC66 exerts its function via regulation of a subset of oncogenes, and knockdown of circCCDC66 inhibited tumor growth and cancer invasion in xenograft and orthotopic mouse models, respectively. Taken together, these findings highlight a novel oncogenic function of circRNA in cancer progression and metastasis.


Introduction
Colorectal cancer was the third most common cancer with approximately 1.4 million new cases in 2012 (1) and was ranked as the second leading cause of cancer-related death in United States of America (2). Despite the intensive investigations and the therapeutic improvements, over 50% of patients with colon cancer ultimately die from this disease, which highlights the necessity to unravel other unknown mechanisms contributing to colon cancer malignancy. Circular RNA (circRNA), a class of RNA molecules with circular configuration formed by either typical spliceosome-mediated or lariat-typed splicing between an upstream splice acceptor and a downstream splice donor, was recently found to be widely spread in eukaryotes (3,4). Circular RNA arises from exonic, intronic and intergenic regions (3). The exonic circRNAs (exclusively having exon sequences) typically reside in the cytoplasm while exon intronic circRNAs (intron-retaining circRNAs) remain in nuclei (5). The exoncontaining circRNAs are the end products of splicing and are the most studied groups to date.
The biogenesis of circRNAs involves several distinct mechanisms. An early study proposes that the accumulation of circRNAs is caused by the retardation of cell proliferation because circRNAs are passively diluted in proliferative cells (6). Others propose that the biogenesis of circRNAs is facilitated by an RNA binding protein like Quaking (7). The formation of exonic circRNAs is facilitated by complementary sequences in the flanking introns (8)(9)(10) and by proteins modulating the proper interaction between upstream and downstream introns (7,(11)(12)(13). However, the roles of these proteins in circRNA dysregulation in cancer have not been fully characterized.
In addition to the biogenesis of circRNA, genome-wide analyses revealed that circRNAs are differentially expressed in various cancerous tissues/cell lines. It has been reported that cancerous cell lines tend to express a more diverse pattern of circRNAs compared to noncancer cell lines (14). Bachmayr-Heyda, et al. reported that the overall levels of circRNAs were globally downregulated in colorectal and ovarian cancers, and were negatively correlated with the disease status and proliferation (6). In adherent to Bachmayr-Heyda's study, two recent studies also demonstrated that circITCH and circFOXO3 are downregulated in esophageal squamous cell carcinoma and in a breast cancer cell line, respectively (15,16). In contrast, the levels of several circRNAs were elevated during the process of epithelial-mesenchymal transition (7) and in cultured primary endothelial cells (17). These studies demonstrate that the regulation of circRNA expression is tightly regulated under distinct circumstances, and that the investigations into circRNA are still in its infancy, and yet the roles of circRNA in cancer are not fully characterized.
Despite the negative correlation between the overall circRNA expression levels and normal/ cancer status, the biological or pathological functions of individual circRNA in a particular cancer remain largely unknown. CircRNAs are suggested to serve as microRNA (miRNA) sponges, regulators of splicing, and as transcription regulators (5,13,18,19); nonetheless, the experimental data supporting these hypotheses are scarce. Herein, we aimed to identify circRNAs that are clinically relevant to colorectal cancer and to investigate the roles of these circRNAs in colorectal cancer pathogenesis.
Surgery at the National Cheng Kung University Hospital, and the non-tumorous tissues were collected at the site 10 cm away from the affected region. The stages and classification for collected tumors were histologically verified by pathologists.

RNA sequencing
Total RNA isolated from 12 paired tumor/adjacent non-tumor specimens were subjected to ribosomal RNA (rRNA) removal using the RiboMinus™ technology (Invitrogen). The rRNA-depleted RNA sample was quantified by the Agilent 2100 Bioanalyzer, and used to construct cDNA template by the SOLiD™ Total RNA-Seq Kit according to the manufacturer's instructions. In brief, the RNA was fragmented using RNase III hydrolysis followed by ligation with strand-specific adapters and reverse transcribed to generate cDNA. Fragments of cDNA between 150 to 250 nucleotides were subsequently isolated by electrophoresis in 6% urea-TBE acrylamide gel. The gel-isolated cDNA was amplified with 15 cycles of PCR to produce enough templates for the SOLiD™ EZ Bead™ system to generate template bead library for ligation-based sequencing on the SOLiD™ Sequencer System. The sequencing raw data ran through the SOLiD Accuracy Enhancement Tool (Life Technologies) for color space correction (20). The color space reads were aligned to the human reference genome (hg19, http://genome.ucsc.edu) using the mapping algorithm implemented in the RNA-Seq pipeline of LifeScope Genomic Analysis Software (Life Technologies). The mapped reads were reported in BAM files and used in the following bioinformatic analysis developed in house for circular RNA identification.

Bioinformatic identification of circRNAs
The steps for identification of circRNAs were illustrated in the Fig. 1A. Briefly, the first 25 bp and last 25 bp were extracted from the BAM files of hard-clipped reads (with H in CIGAR string) and unmapped reads and converted to fastq format, followed by second round of mapping to human genome (hg19) using CLC Genomics Workbench (QIAGEN). The pairs of half-reads mapped to the same genes were used for further analyses. To identify the backsplice junctions in the full-length reads, the exon positions for the pairs of half-reads having the same gene ID with reversed order in the genomic coordinate were used to preassemble the reference sequences as backspliced exons. Full-length read sequences matched to the pre-assembled reference sequences, either one alignment in the pre-assembled exons or two alignments in a single exon (Fig. S1), using BLAST2 (parameters: bl2seq -F F -U Tm T -T F -q -2 -W 28 -G 0 -p blastn -D 1) were considered as the circRNA candidates. The rest of full-reads either with zero or multiple (>2) matches were re-aligned to full length of mRNA, and considered as the circRNA candidates if the matched positions have reversed order.

Analysis of miRNA targeting sites
The miRNA binding sites in circRNAs were analyzed by miRanda (Version 3.3a, parameter: -strict) using the miRNA list from miRBase 21 (downloaded in Feb 2015) (21,22), while miRTarBase (Release 6.0) (23) was used to identify the miRNA targeting sites in 3′-UTR for genes.

Gene Set Enrichment Analysis
The Gene Set Enrichment Analysis was performed according to a previous study (24). The CRC RNA-Seq data were used as expression profile, and the differentially expressed genes with miRNA sites identified in circCCDC66 (designated circCCDC66-regulated genes) were inputted as the gene set for enrichment analysis. Similarly, the gene set of differential expressed genes without miRNA sites identified in circCCDC66 was used as non-circCCDC66-regulated ones.

RNA interference knockdown
Specific knockdown of circCCDC66 was achieved using two siRNA oligonucleotides designed and synthesized by Dharmacon to target the backsplice junction (siJCT1 and siJCT2, see sequences in Fig. 4A and Table S1). HCT116 or HT-29 cells were transfected with desired siRNA oligonucleotides at a final concentration of 50 nM, using Lipofectamine 2000 (Life Technologies) reagent, according to the manufacturer's instructions. Expression of both linear and circular CCDC66 transcripts was quantified using RT-qPCR.

Molecular cloning of circRNA for overexpression and validation
To validate the backsplice junction, a set of primers spanning the splice site was design for each circRNA candidate (Table S2). Amplicon with expected size was cloned to pCRII (Invitrogen) and sequenced. pCIRC2, the circular RNA expression plasmid, was developed by the Yen lab (detail to be presented elsewhere). For overexpression of circCCDC66_10_8, we first used pCIRC2 vector to include part of the intron sequences from CiRS-7 for efficient circularization of transcript. The sequence of exon 8 to 10 was amplified using PCR and cloned to pCIRC2 through MfeI and XmaI sites, followed by reconstitution of the intron/exon junction using 3′ SS F/R and 5′ SS F/R primers (see Table S3 for primer sequences). The full length of CCDC66 exon 8-10 and the flanking intronic sequence were verified using Sanger sequencing.

Xenograft and orthotopic models of colorectal cancer in mice
NOD/SCID mice (NOD.CB17-Prkdc scid /NcrCrl) were purchased from the Animal Center at the College of Medicine, National Cheng Kung University, and approval for the experiments was obtained from the institutional animal care and use committee. HCT116 cells (1 × 10 6 cells) transfected with control oligonucleotides or oligonucleotides against circCCDC66 were resuspended in matrigel (BD, 354230; cells in medium:matrigel = 1:1) and injected either subcutaneously into the dorsal flank or into the surface of cecum of 8 weeks old mice. After 4 weeks, the mice were sacrificed, xenografted tumors, cecum, colon and liver were collected for inspection using haematoxylin and eosin staining, and immunohistochemical staining.

Statistical analyses
The data were expressed as means ± standard deviation of the mean in the bar charts and were analyzed by either student's t test (two samples) or one-way analysis of variance (ANOVA) using GraphPad Prism 5.0 (GraphPad Software, Inc. La Jolla, CA, USA). The categorical data were analyzed by Chi-square test. Results with p values less than 0.05 were considered statistical significance. The receiver operating characteristic (ROC) analysis was applied to evaluate the power of expression of circCCDC66 as prognostic marker of CRC, and the area under curve (AUC) was calculated in Prism 5.0.

Differential expression of circRNAs in colorectal cancer tissues
To identify the expression of circular RNAs in clinical samples from patients with colorectal cancer, we applied an in-house developed pipeline to consider that backspliced exons typically caused the problems to the reads aligner -typically unmapped or partially mapped (clipped reads), thus we focused on those unmapped reads and reads with hard clipping (which indicates a fragment of sequence was missing in the reference sequence). To facilitate the mapping of those reads and identify the potential backsplice events, the first and last 25 bp of the reads were extracted and remapped to the human genome (Fig. 1A). Pairs of half-reads were grouped 'exonic' if they fall in the exon regions, while the rest were grouped 'intergenic/intronic'. Among the 1645 exonic half-reads, the half-reads with reversed orders within the same gene were further manually curated and collapsed to 74 circular RNAs with canonical splice signals (Table S4). The majority (71.62%) of the 74 circRNA candidates sat within the coding region, while 21.62% covered the 5′-UTR and partial CDS, 5.40% covered the 3′-UTR and partial CDS, and 1.35% spanned from the beginning to the end of transcript (Fig. 1B).
As a proof-of-concept, several circRNAs were selected for further investigation. They were chosen because their parental genes are involved in cell cycle regulation (CCNB1 and CDK13), signal transduction (PLCE1 and USP3), epigenetic modulation (KDM1A and KMT2E), or transcriptional regulation (FOXK2 and BMP2K). In addition, we selected the circRNA derived from CCDC66 because there is no known function associated with its parental gene. The presence of these candidate circRNAs were first validated with RT-PCR using outward primer pairs in two colon cancer cell lines, HT-29 and HCT116, and in a set of clinical CRC specimens. Circular RNAs derived from the genes BMP2K, KMT2E and USP3 showed multiple backsplice events, while the rest of circular RNAs yielded a single amplicon ( Fig. 1C and S2A). Each selected circRNA candidate was verified by omitting reverse transcriptase to exclude the potential DNA contamination and treated with RNase R prior to reverse transcription to confirm the lack of free 5′-and 3′-ends ( Fig. 1C, top) and the levels of linear transcript were used to demonstrate the efficacy of RNase R treatment (Fig. 1C, bottom). The backsplice junctions were further confirmed by Sanger sequencing of the excised bands (Fig. S3). Next, differential expression of the circRNA candidates was quantified in 48 pairs of colorectal tumors and corresponding adjacent non-tumor tissues. Among the investigated circRNA candidates, circCCDC66_10_8, circCCNB1_7_6, and circCDK13_2_2 were significantly elevated in the colorectal tumor tissues ( Fig. 1D and 1E) while circFOXK2_3_2, circKDM1A_4_2, circKDM1A_9_2 and circPLCE1_2_2 showed no significant difference (Fig. S2B). Note that we used a set of inward primers to amplify exon 6 to 11 of CCDC66 and did not detect the exon-skipped product (Fig. 1E, bottom), which indicates circCCDC66 is not a product of exon skipping.
To explore the novel functions of circular RNA, we chose circCCDC66 (circCCDC66_10_8) for further investigation as its parental transcript possesses no known biological function. RT-PCR results from multiple cell lines demonstrated that circCCDC66 were expressed in all tested colorectal (Caco-2, HCT116, HT-29, LS123), breast (MCF-7, MDA-MB-231, MDA-MB-468), pancreatic (BxPC-3, MIA PaCa-2) and cervical (HeLa) cancer cell lines except the cell line derived from normal colon epithelial cells (CCD 841 CoN), suggesting that the elevated level of circCCDC66 may be general to tumorigenesis of these cancers (Fig. 1F). Cellular localization analysis revealed that circRNAs mainly resided in the cytoplasm (RNU6B as positive control for nuclear fraction; 18S rRNA for cytoplasmic fraction), a phenomenon that is consistent with the property of exonic circRNAs (Fig. 1G). Next, we quantified the level of the linear transcript of CCDC66 by RT-qPCR to determine if the elevated level of circCCDC66 is due to the higher levels of parental transcripts. Surprisingly, the level of linear CCDC66 was lower in cancer cells compared to adjacent normal tissues (Fig. 1H). To answer whether the low level of linear transcript is due to being processed to generate circular RNA, we analyzed the correlation of linear and circular CCDC66 levels of each patient. The results showed that there is no negative correlation between the levels of circular and linear transcripts (Fig. 1I), which suggests that the low levels of linear CCDC66 in cancer is not caused by a high converting rate to generate circCCDC66. Intriguingly, we found that the level of linear CCDC66 transcript was dramatically decreased in cancer cells treated with hypoxia while the level of circular CCDC66 transcript remained no change (Fig. S4), together caused the fraction of circular CCDC66 enriched (Fig. 1J). The nature of high stability (Fig. 1K) of circular RNA might partly contribute to this differential response to hypoxia between linear and circular transcripts. These data imply that the elevated/accumulated levels of circCCDC66 could be the consequence of recurring intratumoral hypoxic stress, which is further enhanced by the stable nature of circular RNA.
We next addressed the potential clinicopathologic implications of elevated circCCDC66 during the development of colon cancer. We found that the level of circCCDC66 was elevated in the pre-cancerous polyp tissues and became higher in the tumorous tissues ( Fig.  2A). The level of circCCDC66 was increased in all stages of colorectal cancer tissues compared to their paired noncancerous counterparts (Fig. 2B). Of noted, patients with a higher level of circCCDC66 have a poor prognosis as evidenced by lower overall survival rates (Fig. 2C). This phenomenon is specific to circCCDC66 as the level of linear CCDC66 mRNA was not associated with any prognostic markers (Fig. 2D, linear CCDC66 transcript). To better demonstrate the circCCDC66 may be used as a diagnostic marker for CRC, the area under receiver operating characteristic curve (AUC) using the expression levels of circCCDC66 was calculated. The AUC of circCCDC66 indicated that 88% of randomly chosen CRC patient will have higher levels of circCCDC66 (ranked higher) than the randomly chosen normal subjects (Fig. 2E) (25), suggesting the expression level of circCCDC66 is a good predictive biomarker for CRC detection and prognosis.

CircCCDC66 promotes cancer cell proliferation, migration, and metastasis
To investigate the role of circCCDC66 in CRC, we designed two siRNA oligonucleotides to target the unique backsplice junction. The backsplice junction-specific siRNA (siJCT1 and siJCT2) successfully knocked down circCCDC66 but had no effect on the level of linear CCDC66 mRNA (Fig. 3A). The specific knockdowns led to significantly reduced cell number in both HCT116 and HT-29 ( Fig. 3B and 3C), which was mediated through suppression of cell proliferation (Fig. 3D), but not through the induction of apoptosis (Fig.  3E). In addition to cell proliferation, our data demonstrated that knockdown of circCCDC66 significantly reduced cell migration as assayed by the number of cells migrating through the transwell membrane, and also reduced invasion as assayed by the number of cells traveling through the matrigel in both HCT116 and HT-29 cells (Fig. 3F and 3G).
To further characterize the role of circCCDC66, we used pCIRC2, a circRNA expression vector, to express circRNA with exons 8 to 10 of CCDC66 (Fig. 4A). The plasmid was first transfected to 3T3L1 cells, a mouse cell line without detectable circCCDC66, for demonstrate the circularization of designed construct (Fig. 4A and 4B). The specificity of pCIRC2 was further confirmed by mutating splice signals (mutated splice acceptor/mSA and mutated splice donor/mSD), which abolished the expression of circCCDC66 (Fig. 4A and  4B). Transfection of pCIRC2-CCDC66 e8-10 led to the expression of circCCDC66 as confirmed by RT-PCR and Sanger sequencing of the backspliced RNA junction (Fig. 4B). Next, pCIRC2-CCDC66 e8-10 was transfected to HCT116 cells and the expression level of circCCDC66 was quantified ( Fig. 4C and 4D). Overexpression of circCCDC66 significantly increased the number of cell colonies in the soft agar as compared to the circular mCherry RNA that was used as a control (Fig. 4E, upper). In addition, the BrdU incorporation assay showed an increased number of BrdU-positive cells overexpressing circCCDC66 (Fig. 4E,  bottom), suggesting that circCCDC66 stimulates DNA synthesis. Importantly, the overexpression of circCCDC66 had no effect on the levels of the endogenous parental linear CCDC66 transcript, which was measured using two independent sets of inward primers designed to amplify the linear transcripts ( Fig. 4C and 4D, center and right).

CircCCDC66 may serve as a miRNA sponge
To explore the molecular function of circCCDC66, we performed a miRNA target sequence analysis using CCDC66 exon 8 to 10. The analysis revealed that this region was potentially targeted by 99 miRNAs (Fig. 5A and supplementary Table S5). Unlike miR-7 which has multiple binding sites in CDR1as (19), most of miRNAs targeting circCCDC66 have single binding site except miR-670 and miR-4326 (Supplementary Table S5). The distribution of potential miRNA binding sites were illustrated and summarized as cumulative counts at positions where miRNAs bind (Fig. 5B). This prompted us to hypothesize that circCCDC66 may function as a miRNA sponge to protect a specific group of oncogenes from attacks by miRNAs in colorectal cancer. Among these 99 miRNAs, sixty-two are expressed in colorectal cancer tissues according to three independent genome-wide studies -GSE63119, GSE46622 and GSE43793 (26-28) (Fig. 5A, left panel). Of special note, there were many miRNAs (high-numbered) registered after these three genome-wide studies were conducted; thus, the CRC-expressed miRNAs may be underestimated in these three datasets. To test the miRNA sponge hypothesis, we filtered oncogenes and tumor suppressor genes with significant changes in colon cancer (fold change > 1.5, FDR < 0.4 and p < 0.05) and inspected whether these genes could be predicted as targets of miRNAs that have matching sequences in circCCDC66 (Fig. 5A, right panel). The results revealed that 37.5% of oncogenes have the targeting sites of those CRC-expressed miRNAs (circCCDC66.CRC.miRNA+), while only 23.5% of the tumor suppressor genes contain such sites ( Fig. 5C and Fig. S5A, upper panels). This phenomenon of preferential regulation on oncogenes is not restricted to our RNA-Seq expression profile because a similar result was observed when using an independent dataset of CRC from Gene Expression Omnibus (GSE32323) (Fig. 5C and Fig. S5A, lower panels). In contrast, when applying the same analyses using 5000 randomly selected miRNAs and the aforementioned CRC expression profiles, no difference was found between oncogenes and tumor suppressor genes, suggesting that randomly selected miRNAs have no regulatory preference on oncogenes (Fig. 5D). To test whether the predicted circCCDC66-regulated genes are enriched in the tumor transcriptome, gene set enrichment analysis (GSEA) was performed and the results indicated that predicted circCCDC66-regulated genes (by using CRC expressing miRNAs) were enriched in the tumor transcriptome (Fig. 5E, left), while non-circCCDC66-regulated genes showed no enrichment (Fig. 5E, right). Consistent results were observed when using circCCDC66-regulated genes (predicted by using full set of miRNAs) (Fig. S5B).
The circCCDC66-miRNAs regulatory network is a critical component for oncogene activation. Using this regulatory axis of circCCDC66-miRNAs-oncogenes, we revealed that many circCCDC66-regulated oncogenes are upregulated in the CRC tumor specimens (Fig.  6A). As a proof-of-concept, we manipulated the levels of circCCDC66 using overexpression and knockdown approaches and investigated the expression levels of several candidate oncogenes that are known to play pivotal roles in CRC. The levels of DNMT3B, EZH2, MYC, and YAP1, which are the targets of the predicted miRNAs, were found to be significantly upregulated by circCCDC66 overexpression (Fig. 6B) or downregulated by knockdown (Fig. S6A). In contrast, EFNA4 and FGF9, which are not targeted by the predicted miRNAs, were not affected ( Fig. 6B and S6A). Coherently, we found that MYC protein level was increased by the overexpression of circCCDC66 in HCT116 while decreased by circCCDC66 depletion (Fig. 6C and S6B). To further consolidate the role of circCCDC66 as miRNA sponge, we developed a novel in vitro transcription combined with in vivo circularization protocol to test circCCDC66's miRNA sponge effect (Fig. 6D, upper panel). The in vitro transcribed biotinylated linear RNA was efficiently circularized in cells ( Fig. 6D, lower panel), and was used to pull down potential miRNAs. As a proof-of-concept, miR-33b, miR-93, and miR-185, three miRNAs predicted to be sponged by circCCDC66 were tested. Results demonstrated that significant levels of predicted miRNAs were pulled down by biotinylated circCCDC66 (Fig. 6E, upper panel). In contrast, the levels of unrelated miRNA, miR-22, miR-148a and miR-152 were not enriched (Fig. 6E, lower panel). In addition, we cloned the 3′-UTR of human MYC gene to the downstream of renilla luciferase (Fig. 6F). The MYC 3′-UTR activity was significantly increased by the overexpression of circCCDC66 (Fig. 6G) and suppressed by circCCDC66 knockdown (Fig. 6H). In contrast, the non-targeted FGF9 3′-UTR reporter activity was not affected by the overexpression (Fig  6G) or knockdown of circCCDC66 (Fig. S6C), further consolidating the hypothesis that circCCDC66 regulates gene expression as miRNA sponge. To investigate whether this protective effect of circCCDC66 on MYC 3′-UTR activity is mediated through sponging out miR-33b and/or miR-93, inhibitors of two miRNAs were administered to the cells with circCCDC66 knockdown. Knockdown of circCCDC66 released miR-33b and miR-93, which then targeted the 3′-UTR of MYC but not FGF9 transcript and caused the reduction of luciferase activity. Administration of inhibitors of miR-33b or miR-93 relieved the suppressive effect of circCCDC66 knockdown on MYC 3′-UTR activity but had no effect on FGF9 3′-UTR activity ( Fig. 6H and Fig. S6C). Lastly, levels of MYC transcripts were quantified by RT-qPCR to test whether similar phenomenon can be observed in endogenous mRNA. As shown in Fig. 6I, knockdown of circCCDC66 (release miR-33b and miR-93) inhibited MYC mRNA expression while administration of miR-33b and miR-93 inhibitors abolished circCCDC66 knockdown effect. Taken together, these data support the idea that circCCDC66 functions as miRNA sponge to protect the MYC mRNA from the attack of miRNA-33b and miR-93.

CircCCDC66 promotes cancer cell proliferation, migration, and metastasis
To investigate the potential clinical relevance of circCCDC66 in vivo, HCT116 cells with or without transient circCCDC66 knockdown were first injected subcutaneously into dorsal flanks of NOD/SCID mice and allowed to proliferate for four weeks. The growth of tumors inoculated with circCCDC66 knocked-down HCT116 was significantly inhibited (Fig. 7A  and 7B). Interestingly, the expression of circCCDC66 in xenografted HCT116 tumor cells has returned to the pre-knockdown level as examined at the end of week 4 (Fig. S7) whereas the level of MYC remained low (Fig 7C), suggesting that suppression of circCCDC66 at early stage may have great impact on tumor growth. Furthermore, HCT116 cells with or without transient circCCDC66 knockdown were inoculated orthotopically in the cecum of NOD/SCID mice and allowed to proliferate for four weeks (Fig. 7D, top). The tumors formed by the HCT116 cells with the knockdown control aggressively disrupted the lining of smooth muscles and intercalated with the muscle cells (Fig. 7D, bottom-left). In contrast, tumors derived from the cells with the circCCDC66 knockdown typically retained a distinct boundary along the edge of the smooth muscle layer (Fig. 7D, bottom-right). The ability of HCT116 cells to form nodules in the liver, presumably due to the metastasis of the tumor from the cecum to the liver, was significantly reduced with circCCDC66 knockdown (Fig.  7E). Immunohistochemical staining using an antibody specific for pan-keratin showed that the metastasized cells were indeed epithelial cells but not hepatocytes (Fig. 7F).

Discussion
CircRNAs are a group of non-coding RNAs with largely unknown biological functions. In this study, we found that the expression of several circRNAs is enriched in colorectal tumor specimens. Among them, circCCDC66, a circRNA derived from a gene with unknown molecular functions, was upregulated in all stages of colon cancer and negatively correlated with prognosis. In vitro and in vivo studies showed that circCCDC66 accumulated under hypoxia and possesses oncogenic capability in terms of promoting cell proliferation, migration, and metastasis. Bioinformatic and genome-wide analyses revealed that circCCDC66 may preferentially protect a group of oncogenes from being attacked by panels of miRNAs and these circCCDC66-protected genes are positively correlated with colon cancer status. To our knowledge, this is the first report that thoroughly investigates the expression, regulation, function, and clinical implication of circRNA, circCCDC66, in human cancer.
The expression of circRNAs in human diseases is an emerging research topic. One pioneer study showed that the molecular signatures of circRNA among different types of cancer cell lines are highly diverse compared to normal cell lines (14). A few earlier studies reported that circRNAs are downregulated in cancer cells (29)(30)(31) while more and more recent studies demonstrated an enrichment of certain circRNAs in several types of human cancer (7,(32)(33)(34)(35). To investigate the expression and function of circRNAs in colon cancer pathogenesis, we performed RNA-Seq analysis of 12 pairs of CRC and adjacent non-tumor samples. Seventy-four novel exon-containing circRNAs were identified. By employing RT-qPCR quantification, we showed that the expression of a panel of novel circRNAs including circCCDC66, circCCNB1, and circCDK13 was elevated in colorectal cancer cell lines as well as in clinical specimens. RNase R treatment, sequencing verification, and cellular localization analysis confirmed that they are bona fide exonic circRNAs. Our finding is consistent with the results observed in TGF-β-treated human mammary epithelial cells (7) and in hypoxia-treated human umbilical venous endothelial cells (17), which suggest that circRNAs are regulated by complicated mechanisms in disease-and microenvironmentspecific manners.
Although the biogenesis and expression of circRNAs have been intensively investigated in the past couple years, the functions of circRNAs in human cancers are not well understood. In esophageal squamous cell carcinoma, circITCH protects its parental transcript, ITCH, which functions as a tumor suppressor by inhibiting the Wnt/β-catenin pathway-mediated cell proliferation (15). In addition, a recent study discovered that circFOXO3 works as a miRNA sponge to protect Foxo3 mRNA from miRNA attack (16). These two studies exemplify the tumor suppressive function of circRNAs and how the loss of these circRNAs contribute to tumor progression. In contrast, our current study reported that circCCDC66 possesses oncogenic capacity through protecting multiple oncogenes from being attacked by a group of miRNAs. Thus, overexpression of circCCDC66 potentiates multiple tumor characteristics including proliferation, migration, and invasion. These functions greatly contribute to cancer metastasis and malignancy, which supports the observation that high levels of circCCDC66 are associated with poor prognosis of CRC patients. It is worth noting that the relationship of high expression levels and poor prognosis was not seen in linear CCDC66 mRNA, suggesting that it is likely mediated by the non-coding function of circCCDC66. Indeed, gene set enrichment analysis revealed that most circCCDC66protected genes are upregulated in cancer but not in normal colon cells. In vivo study using orthotopic injection of circCCDC66-knockdown colon cancer cells showed that knockdown of circCCDC66 impairs tumor nodule establishment in the liver. Taken together, our study provides multiple lines of evidence to illustrate a novel mechanism of circRNA in cancer progression.
Our finding that circCCDC66 functions as an oncogenic circRNA is different from the reported tumor suppressive circRNAs in several ways (Fig. 7G). First, the expression of linear CCDC66 in tumor specimens is not correlated with circCCDC66, implying that the regulation of CCDC66 transcription and the processing of the circular form are two independent events. In contrast, positive correlation was observed between circular and linear transcripts in the cases of circITCH and circFOXO3 (15,16). Second, the expression of circCCDC66 does not affect the level of the linear transcript in either overexpression or circCCDC66-specific knockdown experiments. This is different from the previously identified function of intron-exonic circRNA, which promotes transcription of parental genes (5). The lack of exon-skipped transcript (exon 7 to exon 11) suggests that it is unlikely a modulator of splicing event or a product of exon-skipping (13,36). One may expect that the expression of circCCDC66 would protect the parental linear transcript from miRNAmediated mRNA cleavage. However, there is growing evidence that miRNAs targeting the coding region and the 3′-UTR have distinct effects on mRNA levels. In contrast to the miRNA binding sites in the 3′-UTR, miRNA-mediated cleavage in CDS is typically hindered by translation machinery (37,38). Nevertheless, the lack of association with the active translation machinery may also make circRNA a better miRNA sponge as compared to its linear transcript. Third, unlike ciRS-7 which has an extremely high density of targeting sites for a single miRNA, circCCDC66 suppresses multiple miRNAs (Supplementary Table  S5) (18). Fourth, circCCDC66 impacts the expression of multiple genes simultaneously while circITCH and circFOXO3 exert their influence on the single mRNA target. Intriguingly, many of the genes that harbor the same miRNA binding sites found in circCCDC66 are oncogenes, while only a few are tumor suppressor genes, implicating that circCCDC66-associated regulatory network benefits the tumor development, and thus is preserved in the tumor cells (Fig. 1F). Similarly, during the preparation of this manuscript, a report by Zheng et al., showed that circHIPK3 sponges multiple miRNAs and promotes cell proliferation in vitro (39), suggesting that the regulation of multiple miRNAs by circRNA could be a general mechanism in cells. However, whether circHIPK3 is also elevated in the tumor tissues, preferentially regulates oncogenes or promotes tumorigenesis in vivo like circCCDC66 remains unsolved and warrants further investigation. Lastly, it is worth noting that the transient knockdown of circCCDC66 had a great suppressive effect on the tumor growth in vivo as evidenced by the return of circCCDC66 but not MYC expression in tumor cells. We reason that this might be due to the complicated gene-gene and genemicroenvironment interactions during the early stage of tumor development. However, detailed mechanisms remain unknown. It would be interesting to investigate whether circCCDC66 exerts such effects by targeting multiple oncogenes or epigenetic regulators and how to translate such findings to clinical application.
In summary, we identified a group of novel circRNAs expressed in colorectal cancer and among these, circCCDC66 negatively correlated with treatment outcomes. CircCCDC66 serves as a sponge to protect multiple oncogenes from being attacked by miRNAs. Our finding that one circRNA contains various binding sites for different miRNAs reveals the complex role of circRNA in cancer malignancy and highlights the importance of investigating the complicated circRNA-miRNA regulatory gene network in cancer progression and treatment efficacy.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.