In order to confirm the clonal status and differential phenotype of the SW620 and SW480 cell lines used in this study, we performed a preliminary experiment using nude mice by intrasplenic (IS) injections to induce tumor development (see section Materials and Methods). Tumors were visible in the liver, spleen, and pancreas (Figure 1A, Table S1) following injection of the CRC cells. Numerous liver metastatic nodules were observed in mice injected with SW620 cells after 4 weeks. In contrast, no or very few nodules were detected in mice injected with NCM460 or SW480 cells (Figure 1A). In the SW620 group, tumors were visible in 7/7 livers, and all of the tumors weighed more than 2.0 g. The mean weight of liver tumors was significantly higher than that of other tumor groups at the same developmental stage (Figure 1B). Taken together, our findings indicate that the tumorigenic and metastatic potential of SW620 cells is higher than that of other cell lines, consistent with previous descriptions of these cells (Bergmann-Leitner and Abrams, 2000; Hewitt et al., 2000; Provenzani et al., 2006).

Circular RNA and mRNA were extracted from the metastatic SW620 cells, their primary SW480 colon carcinoma cells, and normal colorectal cells (NCM460), and sequenced (details in section Materials and Methods). Three biological replicates were performed. For each biological replicate, about 25 and 9 Gb circRNA and mRNA data were generated, respectively (Table S2).

FastQC and fastx-toolkit (see section Materials and Methods) were used to check and control the quality of the sequencing reads, respectively. Principle component analysis (PCA, details in section Materials and Methods) of gene expression showed that the three individual replicates of circRNA-Seq and RNA-seq data for each cell line clustered together, indicating that the samples were highly reproducible (Figure S1). In general, based upon our analysis, the data generated by this study are of high quality and therefore we undertook for following analysis.

Using the previously described CIRI2 circRNA prediction tool (Gao et al., 2017), we were able to identify 33,962 unique circRNAs from our three cell lines (Tables S3, S4). Of the 33,962 unique circRNAs, 25,042 (73.7%) were identified previously in other studies (Chen et al., 2008, 2016; Ghosal et al., 2013; Glazar et al., 2014; Table S3) and 8,920 were newly identified by this study. Of the 33,962 unique circRNAs, 28,340 (83.4%) appear to be derived from protein-coding exons, which is similar to the results of previous studies (Zheng et al., 2016).

Approximately the same number (~14,000) of circRNAs were identified from both of the two CRC lines, while a significantly larger number (25,329) of circRNAs were identified from the normal cell line (Figure 2A, Table S3) despite the fact that it had the lowest raw data amount (66.7 Gb) among three cell lines (Table S2). Of the circRNAs identified from two CRC cell lines, only 20,552 unique circRNAs were found.

Of the 25,329 circRNAs from normal cell lines (NCM460), 13,410 circRNAs were only found in the normal cell lines but not the two CRC cell lines. In contrast, much fewer, about 3,800 specific circRNAs, were identified in each CRC cell line suggesting that circRNAs may be expressed significantly lower in CRC cells relative to their normal cells.

The circRNAs from two CRC cell lines have some unique characteristics when compared to those identified in the normal cell lines. For example, the sizes (genomic lengths) of the circRNAs are significantly smaller than those from the normal line (Figure 2B, P-value < 2.2e-16 by Wilcoxon rank-sum test). We also found that in the two CRC lines, the genomic lengths of circRNAs from metastasis line are much shorter than those from their primary cell line (Figure 2B). Circularization score H (Ivanov et al., 2015) was used to measure the reverse complementary matches of the flanking introns of circRNAs identified from three cell lines. A significantly lower value of score H in the flanking introns of circRNAs from CRC lines than in normal lines was observed (Figure 2C, P-value < 2.2e-16 by Wilcoxon rank-sum test). It has been reported that reverse complementary sequences of flanking introns, contribute to circRNA biogenesis (Ivanov et al., 2015). Usually, the circularization score H of reverse complementary sequences are much higher than normal sequences. Therefore, compared to the normal line, circRNAs may be more likely to be generated by other biogenesis mechanisms rather than reverse complementary sequences in CRC lines, in which predicted circRNAs have significant lower score H of flanking introns.

A total of 2,919 unique DECs were identified between the CRC cells (SW480 or SW620) and their normal cells NCM460 (termed as ncDECs) in this study (Figure 3A, Table S4). Among these unique DECs, there are 2,056 DECs between the SW480 and NCM460, and 1,758 DECs between the SW620 and NCM460 (Figure 3A).

Some interesting results could be found from the expression patterns of 2,919 ncDECs (Figure 3B, Table 1). For example, 840 have similar expression levels in two CRC lines but are significantly decreased (709, Table 1, pattern NO. 1) or increased (131, Table 1, pattern NO. 7) compared with normal cell lines. Fifty-five (Table 1, pattern NO. 11–16) of 2,919 ncDECs not only differentially expressed between cancer and normal cell lines, but also differentially expressed between primary and metastasis cancer cell lines. The 895 mentioned above are common between the two ncDEC sets (i.e., the overlap of blue and green circles in Figure 3A, the expression of which were shown in Figure 3C). In addition, pattern NO. 2 and NO. 4 showed 32.6% circRNAs differentially expressed between primary cancer and normal cell lines; pattern NO. 3 and NO. 5 showed 22.0% differentially expressed between metastasis cancer and normal cell lines.

A total of 623 DECs were identified (FDR < 0.05 and fold-change >= 2) between the metastatic SW620 and its primary SW480 cells (termed as pmDECs; Figure 3A, Table S4). Among these 623 pmDECs, there are 487 (78.2%) that are common with the 2,919 ncDECs, of which 277 are also differentially expressed between NCM460 and SW620, and 265 are differentially expressed between NCM460 and SW480. Only 55 ncDECs (Table 1, pattern NO. 11–16) are differentially expressed in all three cell lines, which indicates that those circRNAs may involve in the progression of cancer.

Of the 623 pmDECs, 348 are down-regulated and 275 up-regulated in the SW620 line (Figures 3D,E). Many circRNAs showed no detectable level of expression in SW620 cells despite being highly expressed in SW480, or no expression in SW480 but highly expressed in SW620 (Figure 3E). For example, circRNAs from two genes, GLI3 (10 circRNAs) and RAPGEF5 (8 circRNAs), are down-regulated in SW620 (or with no detectable level of expression) relative to SW480. Meanwhile, some circRNAs are expressed in SW620 such as circRNAs from SLC22A3 locus (6 circRNAs) but not in SW480.

Twenty-eight of the 623 pmDECs were generated from 24 genes that are included in the cancer Gene Census (http://cancer.sanger.ac.uk/census; Table S5). Both significantly up- and down-regulated expression of these 28 circRNAs in SW620 cell lines relative to SW480 cell lines were observed. Similarly, the expression patterns of the genes from which they derive also have the same expression patterns in the SW620 relative to the SW480 cell lines, with the exception of two genes (GPHN and SPECC1). Each of the 24 genes had multiple loci generating circRNAs, including four genes (NFIB, ETNK1, SPECC1, ARHGAP26) which gave rise to two pmDECs, respectively (Table S5).

CircRNAs that are significantly differentially expressed between primary tumor cell and metastatic tumor cell, and also differentially expressed between normal cell and metastatic tumor cell, might be related to metastasis processes. A total of 277 DECs fell into this category and were used to investigate the relationship between levels of circRNA expression and expression of their host genes. We found a positive correlation between the expression of the 277 pmDECs and their host genes (Figure 4A; i.e., more highly expressed circRNAs have more highly expressed host genes in SW620 relative to SW480). The 277 pmDECs between the SW620 and SW480 were further used to plot expression fold-changes with their host genes (Figure 4B). The results indicate that the DECs between the metastatic SW620 and its primary SW480 cells have a highly positive expression correlation with their host genes (R² = 0.62, P < 2.2e-16, Figure 4B). Based on the expression fold-changes of the ncDECs identified from the two CRC cell lines and their host genes, a positive correlation is also observed (Figure S2).

Many circRNAs were identified from metastasis-related genes (Ramaswamy et al., 2003; Jin H. et al., 2017; Weyden et al., 2017) or from genes that are involve in regulating miRNAs that target metastasis-related genes (Table S6). For example, it is reported by Li J. et al. (2015) that miR-96 can significantly lower the expression of Foxo3 which is associated with lymph node metastasis. We found at least 49 pmDECs that could serve as sponges for miR-96 or other miRNAs that target Foxo3 (Presentation 1, Table S7). The down-regulation of such circRNAs in SW620 would be expected to result in an up-regulation of miR-96 and down-regulation of Foxo3 (1.64-folds, P = 8.0e-27). As in the case of Foxo3 regulation, many complex circRNA-miRNA-mRNA interactions were predicted (Presentation 1, Table S7).

Another example is the SMURF2 network. SMURF2 has a key role in Ubiquitin E3 ligase regulation of the TGF-beta pathway (Jin K. et al., 2017) and was found significantly down-regulated (four-fold, P = 1.5e-246 by Wald test in DESeq2 package) in SW620 metastatic cells relative to SW480 primary tumor cells (Table S4). Two pmDECs from the SMURF2 were detected in our studies, both of which were also significantly down-regulated in SW620 (13.5-fold, P = 8.3e-7 and 11.9-fold, P = 4.4e-7 by Wald test in DESeq2 package) based on their junction spanning reads. These two pmDECs are up-regulated in SW480 compared with NCM460 but not differentially expressed in NCM460 and SW620.

We have analyzed the circRNAs in two CRC cell lines, primary tumor cell lines and their metastasis derivative cell lines. We showed that a large number of circRNAs are present in tumor and metastatic cells that have significantly different expression levels relative to normal cells (2,919 ncDECs), including 895 ncDECs that are significantly differentially expressed in both NCM460/SW480 and NCM460/SW620 pairings. We also found a group of 623 pmDECs that are differentially expressed in primary cells compared to its metastasis cell. Among these are 277 pmDECs that are also differentially expressed between NCM460 and SW620. The expression patterns of these 277 pmDECs are unique among the three cell lines and could potentially be involved in the metastasis processes. In the sets, only about 180 circRNAs are in common with CRC-related circRNAs identified in previous studies (Table 2) and therefore, these circRNAs constitute a new large pool of candidate circRNAs that could serve as biomarkers for CRC development and metastasis. To our knowledge this is the first large-scale identification of metastasis-related circRNAs for CRC.

Interestingly, circRNA circ-001569 previously reported to be involved in CRC (Xie et al., 2016) was not found among our circRNAs. However, we did identify at least 40 different circRNAs (3 DECs) that arise from the same ABCC1 locus that gives rise to circ-001569. Alternative splicing is a major contributing factor to the formation of circRNAs, and we believe that the different isoforms of circ-001569 identified in our experiments are of importance to CRC development and metastasis.

In this study, we compared the circRNAs from normal and CRC cells and found several features which appear to be in common for CRC circRNAs. First, our findings suggest that, in general, circRNAs are expressed significantly lower in CRC cells compared to their normal cells counterparts (P = 8.7e-5; Figure S3). This observation is also in agreement with the findings of Zheng et al. (2016) and Bachmayr-Heyda et al. (2015) and suggest that this expression pattern indicates the role of circRNAs in CRC development. A positive expression relationship (although rare cases with negative relationship) was observed between circRNAs and their host genes based on this study (P < 2.2e-16). It is reasonable to suggest is that the lower expression of circRNAs in CRC cells may be responsive to lower expression of genes involving tumor depression.

Second, CRC circRNA sizes become smaller than that from their normal cell (Figure 2B). Similar trends were also observed in two prior studies on CRC circRNAs in which different CRC lines were used (Bachmayr-Heyda et al., 2015; Zheng et al., 2016). However, we also discovered the circularization score H of flanking introns of CRC circRNAs are significantly lower than that from their normal cell lines (Figure 2C), which was not reported previously. Therefore, we measured the circRNAs from both CRC and normal cells identified by Bachmayr-Heyda et al. (2015) and Zheng et al. (2016), respectively (Wilcoxon rank-sum test), and found same results in circRNA sizes (P < 2.2e-16 and P = 3.94e-7, respectively) but not in score H (Table S8). The lower score H may imply the main biosynthesis mechanisms of circRNAs in CRC cells are not based on reverse complementary sequences. And the changes of score H indicated that it is also a new aspect to study the mechanism of circRNA formation.

Particularly, in terms of sizes and score H of flanking introns of circRNAs, CRC metastasis circRNA sizes further become smaller than that from their primary tumor cells while score H of flanking introns of CRC metastasis circRNAs are also significantly lower than that from their primary tumor cells. This is the first observation for metastasis circRNAs and clearly additional studies in other CRC lines or other tumor tissues are merited to confirm the generality of this observation.

The characteristics of circRNAs in CRC cells and normal cells such as expression level, size and score H of flanking introns, indicates environmental inhibition of circRNAs biosynthesis in CRC cells in multiple perspectives.

Previous studies have shown that circRNAs can be used as good biomarkers for tumor tissues (Wang X. et al., 2015; Ahmed et al., 2016; Kulcheski et al., 2016; Taborda et al., 2017; Zhang et al., 2017). In this study, we identified some circRNAs highly expressed only in CRC lines but not in normal lines or conversely. These circRNAs could serve as good measurable indicators for CRC development and metastasis. Based on expression patterns of circRNAs identified by this study, we can find some circRNAs highly expressed only in CRC SW480 and SW620 lines but not in normal NCM460 line (putative CRC biomarker) or only SW620 but not in SW480 and NCM460 (putative CRC metastasis biomarkers). By this way, top 12 circRNA biomarkers for CRC cells and top 15 biomarkers for CRC metastasis cells were selected and could be used for biomarker development in diagnosis of CRC in future (Table S9).

Normal human colon mucosal epithelium cells NCM460 (INCELL, San Antonio, Texas, USA), primary colorectal cancer cells SW480 and the lymph node metastatic variant cells SW620 were obtained from State Key Lab of Diagnostic and Treatment of Infectious Diseases, the First Affiliated Hospital, Zhejiang University School of Medicine. All media and supplements were purchased from Hyclone (Logan, UT, United States) and Sigma. Anti-CD63 (25682-1-AP) antisera was purchased from Proteintech. NCM460, SW480, and SW620 cells were cultured in DMEM medium supplemented with 10% FBS, 100 IU/mL penicillin-streptomycin at 37°C with 5% CO₂.

This study was carried out in accordance with the recommendations of Institutional Animal Care and Use Committee (IACUC) guidelines. The protocol was approved by the Laboratory Animal Welfare Ethics Committee, Zhejiang University (ZJU20170983). Six-week-old female Balb/C nu-nu mice were obtained from SHANGHAI SLAC LABORATORY ANIMAL CO. LTD (ShangHai, China), and maintained under specific pathogen-free conditions. Mice were injected as previously described (Hewitt et al., 2000). For intrasplenic injection (IS), 5 × 10⁶ cells were obtained and resuspended in 100 uL of PBS. Mice were anesthetized and a midline abdominal incision was made. The spleen was exteriorized via an abdominal midline incision and tumor cell suspension was injected slowly into the spleen. For assays of tumorigenesis and liver metastasis, mice were euthanized 4 weeks after injection and organs were weighed.

Total RNAs from these three cell lines were extracted using RNeasy Mini Kit (QIAGEN). For preparation of RNA-Seq libraries enriched for circRNAs, total RNA was first treated with the Ribo-Zero rRNA Removal Kit (Epicenter) to remove rRNA according to the manufacturer's instructions, and then treated by RNase R (Epicenter) to deplete linear RNA. RNase R treatment was done by incubating 1 μg of rRNA-depleted total RNA with 5 U RNase R in 1 × RNase R buffer at 37°C for 30 min. RNA sequencing libraries (NCM460, SW480, and SW620) for Illumina Hiseq 3000 platform were constructed according to the manufacturer's instructions (Illumina). Three biological replicates were performed.

FastQC and fastx toolkit (Andrews, 2010; HannonLab, 2014) were used to check the quality of the sequencing reads and remove the adaptors and low-quality reads, respectively. RNA-seq and circRNA-seq data by this study have been deposited in GenBank under the accession number or Project PRJNA393626.

In order to evaluate the reproducibility of replicates, data of gene expression were used to do PCA. For each circRNA-seq and mRNA-seq sample, FASTQ reads were first mapped to human genome (GRCh38/hg38) obtained from Ensembl by STAR (Dobin et al., 2013). Then we used HTSeq-count (Anders et al., 2015) to count supporting reads of annotated genes and finally DESeq2 (Love et al., 2014) to draw PCA plots.

For each circRNA-seq sample, FASTQ reads were first mapped to human genome (GRCh38/hg38) by BWA (Li and Durbin, 2009) and then CIRI2 (Gao et al., 2015) was used to identify circRNAs. The number of reads spanning back-spliced junctions was used as an absolute measure of circRNA abundance. While running CIRI2, default parameters were used for samples of cell lines.

According to the annotation of human genome (GRCh38/hg38) obtained from Ensembl (http://www.ensembl.org), all identified circRNAs were annotated. For each circRNA, we searched for transcript fragments that have overlap with the genomic position of circRNA and then defined the corresponding gene of this transcript fragment as the host gene of this circRNA.

In order to evaluate the characteristics of circRNAs from three cell lines (NCM460, SW480, and SW620). We used the algorithm presented by Ivanov et al. (2015) to calculate the score H of flanking introns of circRNAs. Only exonic circRNAs that come from middle exons of transcripts were used, based on the calculation principles of this algorithm. Therefore, in order to reduce the errors caused by total circRNA numbers, all score H have been weighted average adjusted by the total number of circRNAs identified in each cell lines.

For each circRNA-seq sample, DESeq2 (Love et al., 2014) was used to analyze the significant differences of gene expression among three cell lines (NCM460, SW480, and SW620). To estimate the relative expression of a circRNA, we normalized the number of reads spanning the back-spliced junction to the total number of mapped reads (units in billion) and read length. Therefore, the formula of a circRNA's relative expression is “number of circular reads/number of mapped reads (units in billion)/read length.” Here, we defined circRNAs with |log2 FC| > 1 and padj < 0.05 (P-value adjusted for multiple testing using Benjamini-Hochberg) as differential expressed circRNAs (DECs).

For each RNA-seq sample, FASTQ reads were first mapped to human genome (GRCh38/hg38) obtained from Ensembl by STAR (Dobin et al., 2013). Then we used HTSeq-count (Anders et al., 2015) to count supporting reads of genes and DESeq2 (Love et al., 2014) to analyze the significant differences of gene expression among three groups (NCM460, SW480, and SW620). We defined genes with |log2FC| > 1 and padj < 0.05 (P-value adjusted for multiple testing using Benjamini-Hochberg) as differential expressed genes (DEGs).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.