All samples were collected from patients between January, 2013 and December, 2015 in the Department of Gastrointestinal and Hernia Surgery, the First Affiliated Hospital of Kunming Medical University (Kunming, China). The project was reviewed and approved by the Ethics Committee of the First Affiliated Hospital of Kunming Medical University. All patients included in the present study provided written informed consent prior to surgery.

Primary CRC samples were obtained from 3 male patients (age, 66.33±7.37 years) with CRC and lung metastasis (Exp A) and 3 male patients (age, 64.67±11.93 years) with CRC without metastasis (Ctrl A), respectively. All patients were histologically confirmed to have CRC based on colonoscopy (26) and did not receive any other forms of therapy at the time of enrollment. At the time of surgery, all tissue samples were immediately frozen in liquid nitrogen and stored at −80°C for further use.

According to the manufacturer's protocol, total RNA was extracted from each sample using a homogenizer and TRIzol regent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA). Subsequently, the quantification and quality of purified RNA were assessed using a NanoDrop ND-1000 (NanoDrop Technologies; Thermo Fisher Scientific, Inc.). The integrity of RNA was assessed by electrophoresis on a denaturing 1.5% agarose gel.

Sample labeling and human circRNA array hybridization were performed according to the manufacturer's protocol (Arraystar, Inc., Rockville, MD, USA). Briefly, total RNA were digested with RNase R (Epicentre; Illumina, Inc., San Diego, CA, USA) to remove linear RNA and enrich circRNA. Subsequently, the enriched circRNA were amplified and transcribed into fluorescent cRNA utilizing a random priming method (Arraystar Super RNA Labeling kit; Arraystar, Inc.). The labeled cRNA were purified using an RNeasy mini kit (Qiagen GmbH, Hilden, Germany). The concentration and specific activity of the labeled cRNA (pmol Cy3/µg cRNA) were measured using a NanoDrop ND-1000. A total of 1 µg of each labeled cRNA was fragmented by adding 5 µl 10X blocking agent and 1 µl 25X fragmentation buffer (Arraystar Super RNA Labeling kit; Arraystar, Inc.). Subsequently, the mixture was heated at 60°C for 30 min, and finally 25 µl 2X hybridization buffer was added to dilute the labeled cRNA. Following this, 50 µl hybridization solution was dispensed into the gasket slide and assembled to the circRNA expression microarray slide. The slides were incubated for 17 h at 65°C in an Agilent Hybridization Oven (Agilent Technologies, Inc., Santa Clara, CA, USA). The hybridized arrays were washed using Agilent wash buffer 1 and wash buffer 2, fixed and scanned using the Agilent Scanner G2505C (both from Agilent Technologies, Inc.).

Scanned images were imported into Agilent Feature Extraction software (version 11.0.1.1; Agilent Technologies, Inc.) for raw data extraction. When comparing two groups of profile differences (such as CRC with lung metastasis vs. without lung metastasis), the fold change (i.e. the ratio of the group averages) between the groups for each circRNA is computed. The statistical significance of the difference may be conveniently estimated by a t-test. circRNA with fold changes ≥2 and P-values of <0.05 were selected as being significantly differentially expressed.

GO analysis (geneontology.org) may be used to construct meaningful annotation of gene products in a wide variety of organisms. GO contains three domains, including biological process (BP), cellular components (CC) and molecular function (MF). The -log₁₀ (P-value) indicates the enrichment score representing the significance of GO term enrichment among genes producing differentially expressed circRNA. Fold enrichment represents the proportion of the changes in genes as proportions in the GO database in this function (such as, the greater the proportion, the more reliability, which is seen as the more significant changes in BP, CC or MF).

KEGG (genome.jp/kegg/) was utilized to harvest pathway clusters covering our knowledge of the molecular interaction and reaction networks in genes producing differentially expressed circRNA.

In the present study, circRNA/miRNA interactions were predicted with Arraystar's home-made miRNA target prediction software (Arraystar, Inc.) based on TargetScan (targetscan.org/vert_71/) and miRanda (microrna.org/microrna/home.do). The differentially expressed circRNA within all the comparisons were annotated in detail with the circRNA/miRNA interaction information.

High-throughput microarray is an efficient approach for studying the biological function of RNA. As demonstrated in Fig. 1A, a total of 431 circRNA were detected to be differentially expressed with fold change ≥2.0, P<0.05 and false discovery rate <0.05 in the Exp A group compared with the Ctrl A group. Among them, 192 and 239 circRNA were upregulated and downregulated, respectively. Among the dysregulated circRNA, in the light of their relation with their coding genes, the circRNA were classified into five categories: 331 were exonic, 45 were intronic, 32 were sense overlapping, eight were intergenic and 15 were antisense (Fig. 1B). According to the log fold change, the top 10 upregulated and downregulated circRNA were selected (Fig. 1C).

According to the number of linear counterparts of differentially overexpressed circRNA involved in BPs, the top 10 BPs were classified (Fig. 2A). Among them, 21 linear counterparts of differentially overexpressed circRNA were involved in transcription, DNA-templated, 15 were involved in regulation of transcription, DNA-templated, 15 were involved in small molecule metabolic process, 14 were involved in viral process, 11 were involved in the regulation of transcription from RNA polymerase II promoter, 10 were involved in signal transduction, 9 were involved in cellular protein metabolic process, 9 were involved in gene expression, 9 were involved in positive regulation of transcription from RNA polymerase II promoter and 8 were involved in DNA repair. According to the enrichment score, the top 10 BPs were selected (Fig. 2B). These included regulation of transcription from RNA polymerase II promoter, viral process, protein refolding, protein desumoylation, chromatin assembly or disassembly, response to ethanol, regulation of actin filament polymerization, regulation of smooth muscle cell proliferation, positive regulation of mitochondrial depolarization and double-strand break repair. Additionally, according to the fold enrichment, the top 10 BPs were selected (Fig. 2C). These included protein desumoylation, chromatin assembly or disassembly, regulation of actin filament polymerization, regulation of smooth muscle cell proliferation, positive regulation of mitochondrial depolarization, cardiovascular system development, cellular response to indole-3-methanol, sarcoplasmic reticulum calcium ion transport, valine metabolic process and histone-serine phosphorylation.

According to the number of linear counterparts of down-regulated circRNA involved in BPs, the top 10 BPs were classified (Fig. 3A). Among them, 25 linear counterparts of differentially downregulated circRNA were involved in transcription, DNA-templated, 24 were involved in small molecule metabolic process, 15 were involved in regulation of transcription, DNA-templated, 14 were involved in gene expression, 13 were involved in negative regulation of transcription from RNA polymerase II promoter, 13 were involved in positive regulation of transcription from RNA polymerase II promoter, 12 were involved in viral process, 12 were involved in signal transduction, 11 were involved in protein transport and 11 were involved in small GTPase-mediated signal transduction. According to the enrichment score, the top 10 BPs were selected (Fig. 3B). These included cellular response to cholesterol, regulation of cell motility, protein transport, response to purine-containing compound, cardiac muscle cell differentiation, neural tube closure, positive regulation of cytoplasmic mRNA processing body assembly, apoptotic process involved in patterning of blood vessels, mammary gland duct morphogenesis and positive regulation of B cell receptor signaling pathway. Furthermore, according to the fold enrichment, the top 13 BPs were selected (Fig. 3C). These included response to purine-containing compound, positive regulation of cytoplasmic mRNA processing body assembly, apoptotic process involved in patterning of blood vessels, mammary gland duct morphogenesis, positive regulation of B cell receptor signaling pathway, B cell apoptotic process, cellular response to cholesterol, regulation of cell motility, activated T cell proliferation, mitotic nuclear envelope reassembly, positive regulation of receptor recycling, cellular response to nutrient levels and positive regulation of osteoblast proliferation.

According to the number of linear counterparts of differentially overexpressed circRNA identified in CC, the top 10 CC categories were classified (Fig. 4A). Among them, 51 linear counterparts of differentially overexpressed circRNA were involved in the cytoplasm, 50 were involved in the nucleus, 37 were involved in the cytosol, 35 were involved in the nucleoplasm, 31 were involved in plasma membrane, 30 were involved in the membrane, 24 were involved in the extracellular exosome, 23 were involved in the integral component of the membrane, 15 were involved in the mitochondrion and 14 were involved in the nucleolus. According to enrichment score, the top 10 CCs were selected (Fig. 4B). These included cytoplasmic membrane-bounded vesicle, membrane, cyclin-dependent protein kinase activating kinase holoenzyme complex, cell, pericentric heterochromatin, midbody, endosome lumen, chromosome (telomeric region), nucleoplasm and lysosomal membrane. Furthermore, according to the fold enrichment, the top 13 CCs were selected (Fig. 4C). These included cyclin-dependent protein kinase activating kinase holoenzyme complex, pericentric heterochromatin, GATOR2 complex, U2-type prespliceosome, transferase complex, Sec61 translocon complex, endosome lumen, pyruvate dehydrogenase complex, telomerase holoenzyme complex, synaptic cleft, nonhomologous end joining complex, elongator holoenzyme complex and lipopolysaccharide receptor complex.

According to the number of linear counterparts of differentially downregulated circRNA identified in CCs, the top 10 CCs were classified (Fig. 5A). Among them, 78 linear counterparts of differentially downregulated circRNA were involved in the cytoplasm, 69 were involved in the nucleus, 56 were involved in the cytosol, 43 were involved in the nucleoplasm, 41 were involved in the plasma membrane, 37 were involved in the extracellular exosome, 33 were involved in the membrane, 31 were involved in the integral component of the membrane, 26 were involved in mitochondrion and 19 were involved in the nucleolus. According to the enrichment score, the top 10 CC categories were selected (Fig. 5B). These included focal adhesion, cytosol, cytoplasm, mitochondrion, cell-cell adherens junction, nucleolus, extracellular vesicle, nuclear euchromatin, stress fiber and dendritic spine head. Additionally, according to the fold enrichment, the top 10 CC categories were selected (Fig. 5C). These included dendritic spine head, dynactin complex, spindle pole centrosome, Fuse Toes complex, U2-type prespliceosome, cerebellar mossy fiber, nonhomologous end joining complex, telomerase holoenzyme complex, intracellular organelle and euchromatin.

According to the number of linear counterparts of differentially overexpressed circRNA identified in MF, the top 10 MFs were classified (Fig. 6A). Among them, 90 linear counterparts of differentially overexpressed circRNA were involved in protein binding, 23 were involved in ATP binding, 21 were involved in zinc ion binding, 19 were involved in poly(A) RNA binding, 19 were involved in DNA binding, 17 were involved in metal ion binding, 13 were involved in transcription factor activity, sequence-specific DNA binding, 8 were involved in protein complex binding, 8 were involved in ubiquitin protein ligase binding and 6 were involved in transcription regulatory region DNA binding. According to the enrichment score, the following top 10 MFs were selected (Fig. 6B): Protein complex binding, SUMO-specific protease activity, C-C chemokine receptor type 5 (CCR5) chemokine receptor binding, ubiquitin protein ligase binding, zinc ion binding, DNA replication origin binding, calcium-transporting ATPase activity, oxidoreductase activity [acting on paired donors, with incorporation or reduction of molecular oxygen, NAD(P)H as one donor, and incorporation of one atom of oxygen], telomeric DNA binding and androgen receptor binding. Additionally, according to the fold enrichment, the following top 10 MFs were identified (Fig. 6C): Androgen receptor binding, CCR5 chemokine receptor binding, DNA replication origin binding, calcium-transporting ATPase activity, oxidoreductase activity [acting on paired donors, with incorporation or reduction of molecular oxygen, NAD(P)H as one donor, and incorporation of one atom of oxygen], phosphatase inhibitor activity, nuclear factor (NF) of activated T cells protein binding, calcium, potassium: sodium antiporter activity, protein methyltransferase activity and glycosylphosphatidylinositol-linked ephrin receptor activity.

According to the number of linear counterparts of differentially downregulated circRNA identified in MF, the top 10 MFs were classified (Fig. 7A). Among them, 120 linear counterparts of differentially downregulated circRNA were involved in protein binding, 29 were involved in poly(A) RNA binding, 27 were involved in ATP binding, 23 were involved in metal ion binding, 18 were involved in DNA binding, 16 were involved in zinc ion binding, 11 were involved in identical protein binding, 10 were involved in transcription factor activity, sequence-specific DNA binding, 9 were involved in ubiquitin-protein transferase activity and 9 were involved in protein homodimerization activity. According to the enrichment score, the following top 10 MF categories were selected (Fig. 7B): Poly(A) RNA binding, protein binding, ligand-dependent nuclear receptor transcription coactivator activity, transcription regulatory region DNA binding, translation elongation factor activity, protein complex binding, ubiquitin-protein transferase activity, transferase activity, carboxy-lyase activity and type 1 angiotensin receptor binding. According to the fold enrichment, the following top 10 MFs were selected (Fig. 7C): Carboxy-lyase activity, type 1 angiotensin receptor binding, 5′-deoxyribose-5-phosphate lyase activity, biotin carboxylase activity, protein methyltransferase activity, mitogen-activated protein kinase p38 binding, death effector domain binding, protein-arginine N-methyltransferase activity, sodium channel inhibitor activity and neuroligin family protein binding.

In order to understand the pathways that differentially expressed circRNA involved in lung metastasis in CRC, KEGG pathway analysis was performed to analyze the genes producing differentially expressed circRNA. According to the enrichment score, the top 10 pathways in which the genes producing upregulated circRNA were involved in are demonstrated in Fig. 8A. While, based on the enrichment score, the top 10 pathways in which the genes producing downregulated circRNA were involved are demonstrated in Fig. 8B.

To investigate the circRNA/microRNA interactions, the binding sites of circRNA were predicted. As demonstrated in Table I, it was identified that hsa_circRNA_105055 (upregulated), hsa_circRNA_086376 (downregulated) and hsa_circRNA_102761 (downregulated) could bind with miR-7 regulating target genes PRKCB, EPHA3, BRCA1 and ABCC1. Additionally, hsa_circRNA_102348 and hsa_circRNA_102562 could bind with miR-410 regulating target genes PDHX, ATM and ABCC1.