Circular RNA in Chemonaive Lymph Node Negative Colon Cancer Patients
Journal: 2021/April - Cancers
Abstract:
Circular RNAs (circRNAs) appear important in tumor progression of colon cancer (CC). We identified an extensive catalog of circRNAs in 181 chemonaive stage I/II colon tumors, who underwent curative surgery between 2007 and 2014. We identified circRNAs from RNAseq data, investigated common biology related to circRNA expression, and studied the association between circRNAs and relapse status, tumor stage, consensus molecular subtypes (CMS), tumor localization and microsatellite instability (MSI). We identified 2606 unique circRNAs. 277 circRNAs (derived from 260 genes) were repeatedly occurring in at least 20 patients of which 153 showed a poor or even negative (R < 0.3) correlation with the expression level of their linear gene. The circular junctions for circSATB2, circFGD6, circKMT2C and circPLEKHM3 were validated by Sanger sequencing. Multiple correspondence analysis showed that circRNAs were often co-expressed and that high diversity in circRNAs was associated with favorable disease-free survival (DFS), which was confirmed by Cox regression analysis (Hazard Ratio (HR) 0.60, 95% CI 0.38-0.97, p = 0.036). Considering individual circRNAs, absence of circMGA was significantly associated with relapse, whereas circSATB2, circNAB1, and circCEP192 were associated with both MSI and CMS. This study represents a showcase of the potential clinical utility of circRNAs for prognostic stratification in patients with stage I-II colon cancer and demonstrated that high diversity in circRNAs is associated with favorable DFS.
Keywords: CMS; circular RNA; colon cancer; prognostic stratification.
Relations:
Content
References
(57)
Diseases
(2)
Conditions
(2)
Anatomy
(1)
Similar articles
Articles by the same authors
Discussion board

Circular RNA in Chemonaive Lymph Node Negative Colon Cancer Patients

1. Introduction

Colon cancer is one of the most common types of cancer with over 1 million new cases worldwide and around 9800 new cases in the Netherlands in 2018 [1]. Up to 21% of patients with stage I–II colon cancer and 40% of patients with stage III colon cancer will develop metastatic disease after curative surgery [2]. As much of the disparity in prognosis for clinically comparable patients remains unexplained, efforts are being directed at finding so far unknown factors that may play a role in the development and progression of colon cancer.

Transcriptome sequencing studies have identified many short and long RNAs with non-protein-coding ability [3,4,5,6]. These non-coding RNAs (ncRNAs) have received increasing attention in recent years due to their aberrant expression features associated with colorectal cancer (CRC) carcinogenesis [7,8,9]. Recent studies have shown that non-coding microRNAs can function as promising biomarkers for stage II [10,11] and stage I–II colon cancer patients [12]. Circular RNAs (circRNAs) represent a re-discovered, abundant class of non-coding RNA molecules [13]. Altered expression of circRNAs is observed in cancer tissue compared to normal tissue [14,15,16,17,18,19], and particularly in CRC [20]. circRNA biogenesis derives from back-splicing, but the regulation and the frequency of this event are under investigation [20]. circRNAs form covalently closed, continuous loop structures produced through an end-to-end formation during transcription [21,22,23]. Increasing evidence shows that circRNAs can function as miRNA sponges, transcription regulators, and interfere with splicing, as well [20]. They are conserved, abundant and often exhibit tissue-, developmental-, and stage-specific expression [24,25].

Due to their special circular structure, circRNAs are usually more stable than linear RNAs and are not easily degraded by exonucleases. They have been proven to remain stable in saliva, blood, and exosomes, which makes them promising biomarkers for the diagnosis, prognosis, and therapeutic assessment of cancer patients [26]. Recently, two novel circRNAs, both derived from the gene BCL2L12, were identified as biomarkers for stage II CRC patients [27]. Compared to conventional available cancer biomarkers (e.g., PSA and CEA), circRNAs are expected to have higher sensitivity and specificity in diagnosis and prognosis [28].

Taken together, these characteristics indicate that circRNAs could represent new clinical diagnostic and prognostic markers, and possibly provide new leads for the treatment of diseases. In this study we describe the identification of an extensive catalog of circRNAs in a large cohort of 181 chemonaive, stage I/II primary colon tumors and related these to tumor stage, localization, Consensus Molecular subtypes (CMS), microsatellite instability (MSI) status and clinical outcome.

2. Materials and Methods

2.1. Study Population and Patient Selection

Fresh-frozen tumor tissue was collected from 181 patients with stage I–II colon cancer undergoing curative surgery. These patients had been enrolled in the MATCH-study—a prospective multicenter cohort study in seven hospitals in the region of Rotterdam, the Netherlands-between 2007 and 2014. Patients have given informed consent on the storage and use of tissue samples, and the collection of clinical data for research purposes. The MATCH study was approved by the Erasmus MC IRB (MEC-2007-088). Inclusion criteria and additional clinical characteristics have been described [29].

Disease-free survival (DFS) was defined as the time elapsed between the date of surgery and either the date of any recurrence of disease or the date of the last follow-up visit at which a patient was considered to have no recurrence.

2.2. Sample Collection and Processing

Sample collection and processing, as well as RNA isolation and RNA sequencing have been described in detail previously [30,31,32]. All samples were reviewed by a pathologist (CHMvD) to ensure the presence of sufficient tumor cells (≥40%). Only samples with an RNA integrity number of at least 7.0 were selected for RNAseq analysis. RNA integrity numbers were assessed using the MultiNA Microchip Electrophoresis system (Shimadzu, Kyoto, Japan) [33].

2.2.1. Microsatellite Instability

MSI analyses have previously been performed and described [32]. In short, the MSI analyses made use of the MSI Analysis System from Promega©, which is a fluorescent PCR-based assay for detection of microsatellite instability in seven markers, including five mononucleotide repeat markers (BAT-25, BAT-26, NR-21, NR-24 and MONO-27) and two pentanucleotide repeat markers (Penta C and Penta D) [34].

2.2.2. Consensus Molecular Subtypes

The CMS classification was performed using the “CMSclassifier” package (https://github.com/Sage-Bionetworks/CMSclassifier, accessed on 23 May 2019), using the single-sample prediction parameter [35].

2.2.3. Identification of circRNAs

The methodology used to identify circRNA reads has previously been described in detail [36]. In short, the method developed by Smid et al. uses sequence reads that have a “secondary alignment” (SA) tag. When using paired-end sequence data, and assuming a circRNA molecule is present, the sequence read that aligns over the crossing of the junction would “point toward” its read-mate somewhere in the circle. Aligning these reads to the linear reference, the junction read will get an SA tag which will be assigned to two locations if and only if this is the one and unique alignment configuration the STAR software can find [37,38]. The read-mate aligns somewhere in between these two locations. Finding additional read pairs showing this configuration with a breakpoint at the exact same location strengthens the evidence for circular transcripts. We included only regions with at least five reads crossing the circular junction. After filtering, GENCODE annotation was used to obtain the exon locations of genes that exactly matched to the circular region. For each sample, STAR also gives the raw read counts for all genes. These were normalized (Trimmed Mean of M-values implemented in edgeR [39]), and the normalized read counts were used to correlate with the number of junction reads of the circular transcripts. The script is also available at https://bitbucket.org/snippets/MSmid/Le949d/identify-circularrna-reads (accessed on 30 October 2018).

2.2.4. Multiple Correspondence Analysis (MCA)

For a substantial number of genes, only a linear transcript is detected in the majority of samples, which results in many missing values per circRNA. This in turn, complicates the use of standard cluster analysis for the identification of sample groups with similar circRNA-related biology. Therefore, circRNA data were considered categorical, i.e., a circRNA was scored as either “present” or “absent” in a sample. These categorical data are suitable for a multiple correspondence analysis (MCA), which is a generalized principle component analysis. An MCA generates a combined plot that shows both patients and circRNAs in such a way that patients and circRNAs that have similar patterns are closer together. Thus, the colon cancer tumor samples and circRNAs are projected onto the same plane, in which the relative distance to either the samples or the circRNAs is meaningful. The 0,0 point corresponds to a sample or circRNA with an average profile. The R-package “ade4” was used to perform the MCA in R version 3.4.1. Custom functions to plot the MCA results are available upon request of the authors.

2.2.5. Reverse Transcription, Quantitative PCR, and Sanger Sequencing

Candidate circular RNAs were selected, and divergent primers were designed that are only able to amplify and detect the circular and not the corresponding linear mRNA (Table S1). Total RNA, isolated with RNA-Bee according to the manufacturer’s instructions (CS105B, TEL TEST), was reverse transcribed into cDNA with the H-minus RevertAid First Strand cDNA Synthesis Kit (K1632, ThermoFisher Scientific, Waltham, MA, USA), followed by an RNAse-H step (AM2293; Ambion). For Sanger sequencing, cDNA from five individual patient samples were used for PCR. cDNA was amplified for 35 cycles using Phusion high-fidelity DNA polymerase (ThermoFisher Scientific) in a total reaction volume of 25 μL, containing 400 nM of each primer and 160 µM dNTPs. For every circRNA two resulting PCR amplicons were purified from gel using the Qiaquick Gel Extraction Kit from Qiagen (Hilden, Germany) according to the manufacturer’s protocol and subjected to Quick Shot Sanger Sequencing by BaseClear BV (Leiden, The Netherlands).

2.3. Statistical Analyses

As indicated above for a substantial number of genes, only a linear transcript was detected in the majority of samples, which results in many missing values per circRNA. This hampers statistical analyses of circRNA expression levels and therefore we categorized the circRNA data into “present” or “absent” for statistical evaluation as well. We used circRNAs present in at least 20 samples to ensure a sufficient number of events for subsequent statistical analyses. STATA version 14 and SPSS Version 24.0 (SPSS, Inc., Chicago, IL, USA) were used to perform the statistical tests that are also indicated in the text. Cox’s proportional-hazards regression was used to evaluate the (log-transformed) number of uniquely present circRNAs per sample, hereafter called circRNA diversity, with DFS, or as “present”/“absent” when evaluating individual circRNAs. Survival curves were evaluated using the logrank test (for individual circRNAs) or with the logrank test for trend (for circRNA diversity, after dividing into three equal quantiles). Pearson’s correlation was used to correlate the circRNA expression with the expression of the linear gene it was derived from. Analyses between categorical variables (like present/absent of a circRNA versus MSI yes/no) were analyzed using Fisher’s exact test. Reported p-values are two-sided and considered significant at p ≤ 0.05. p-values were corrected for multiple testing using Benjamini–Hochberg’s FDR correction when evaluating multiple circRNAs, which were considered significant at p < 0.10.

2.1. Study Population and Patient Selection

Fresh-frozen tumor tissue was collected from 181 patients with stage I–II colon cancer undergoing curative surgery. These patients had been enrolled in the MATCH-study—a prospective multicenter cohort study in seven hospitals in the region of Rotterdam, the Netherlands-between 2007 and 2014. Patients have given informed consent on the storage and use of tissue samples, and the collection of clinical data for research purposes. The MATCH study was approved by the Erasmus MC IRB (MEC-2007-088). Inclusion criteria and additional clinical characteristics have been described [29].

Disease-free survival (DFS) was defined as the time elapsed between the date of surgery and either the date of any recurrence of disease or the date of the last follow-up visit at which a patient was considered to have no recurrence.

2.2. Sample Collection and Processing

Sample collection and processing, as well as RNA isolation and RNA sequencing have been described in detail previously [30,31,32]. All samples were reviewed by a pathologist (CHMvD) to ensure the presence of sufficient tumor cells (≥40%). Only samples with an RNA integrity number of at least 7.0 were selected for RNAseq analysis. RNA integrity numbers were assessed using the MultiNA Microchip Electrophoresis system (Shimadzu, Kyoto, Japan) [33].

2.2.1. Microsatellite Instability

MSI analyses have previously been performed and described [32]. In short, the MSI analyses made use of the MSI Analysis System from Promega©, which is a fluorescent PCR-based assay for detection of microsatellite instability in seven markers, including five mononucleotide repeat markers (BAT-25, BAT-26, NR-21, NR-24 and MONO-27) and two pentanucleotide repeat markers (Penta C and Penta D) [34].

2.2.2. Consensus Molecular Subtypes

The CMS classification was performed using the “CMSclassifier” package (https://github.com/Sage-Bionetworks/CMSclassifier, accessed on 23 May 2019), using the single-sample prediction parameter [35].

2.2.3. Identification of circRNAs

The methodology used to identify circRNA reads has previously been described in detail [36]. In short, the method developed by Smid et al. uses sequence reads that have a “secondary alignment” (SA) tag. When using paired-end sequence data, and assuming a circRNA molecule is present, the sequence read that aligns over the crossing of the junction would “point toward” its read-mate somewhere in the circle. Aligning these reads to the linear reference, the junction read will get an SA tag which will be assigned to two locations if and only if this is the one and unique alignment configuration the STAR software can find [37,38]. The read-mate aligns somewhere in between these two locations. Finding additional read pairs showing this configuration with a breakpoint at the exact same location strengthens the evidence for circular transcripts. We included only regions with at least five reads crossing the circular junction. After filtering, GENCODE annotation was used to obtain the exon locations of genes that exactly matched to the circular region. For each sample, STAR also gives the raw read counts for all genes. These were normalized (Trimmed Mean of M-values implemented in edgeR [39]), and the normalized read counts were used to correlate with the number of junction reads of the circular transcripts. The script is also available at https://bitbucket.org/snippets/MSmid/Le949d/identify-circularrna-reads (accessed on 30 October 2018).

2.2.4. Multiple Correspondence Analysis (MCA)

For a substantial number of genes, only a linear transcript is detected in the majority of samples, which results in many missing values per circRNA. This in turn, complicates the use of standard cluster analysis for the identification of sample groups with similar circRNA-related biology. Therefore, circRNA data were considered categorical, i.e., a circRNA was scored as either “present” or “absent” in a sample. These categorical data are suitable for a multiple correspondence analysis (MCA), which is a generalized principle component analysis. An MCA generates a combined plot that shows both patients and circRNAs in such a way that patients and circRNAs that have similar patterns are closer together. Thus, the colon cancer tumor samples and circRNAs are projected onto the same plane, in which the relative distance to either the samples or the circRNAs is meaningful. The 0,0 point corresponds to a sample or circRNA with an average profile. The R-package “ade4” was used to perform the MCA in R version 3.4.1. Custom functions to plot the MCA results are available upon request of the authors.

2.2.5. Reverse Transcription, Quantitative PCR, and Sanger Sequencing

Candidate circular RNAs were selected, and divergent primers were designed that are only able to amplify and detect the circular and not the corresponding linear mRNA (Table S1). Total RNA, isolated with RNA-Bee according to the manufacturer’s instructions (CS105B, TEL TEST), was reverse transcribed into cDNA with the H-minus RevertAid First Strand cDNA Synthesis Kit (K1632, ThermoFisher Scientific, Waltham, MA, USA), followed by an RNAse-H step (AM2293; Ambion). For Sanger sequencing, cDNA from five individual patient samples were used for PCR. cDNA was amplified for 35 cycles using Phusion high-fidelity DNA polymerase (ThermoFisher Scientific) in a total reaction volume of 25 μL, containing 400 nM of each primer and 160 µM dNTPs. For every circRNA two resulting PCR amplicons were purified from gel using the Qiaquick Gel Extraction Kit from Qiagen (Hilden, Germany) according to the manufacturer’s protocol and subjected to Quick Shot Sanger Sequencing by BaseClear BV (Leiden, The Netherlands).

2.2.1. Microsatellite Instability

MSI analyses have previously been performed and described [32]. In short, the MSI analyses made use of the MSI Analysis System from Promega©, which is a fluorescent PCR-based assay for detection of microsatellite instability in seven markers, including five mononucleotide repeat markers (BAT-25, BAT-26, NR-21, NR-24 and MONO-27) and two pentanucleotide repeat markers (Penta C and Penta D) [34].

2.2.2. Consensus Molecular Subtypes

The CMS classification was performed using the “CMSclassifier” package (https://github.com/Sage-Bionetworks/CMSclassifier, accessed on 23 May 2019), using the single-sample prediction parameter [35].

2.2.3. Identification of circRNAs

The methodology used to identify circRNA reads has previously been described in detail [36]. In short, the method developed by Smid et al. uses sequence reads that have a “secondary alignment” (SA) tag. When using paired-end sequence data, and assuming a circRNA molecule is present, the sequence read that aligns over the crossing of the junction would “point toward” its read-mate somewhere in the circle. Aligning these reads to the linear reference, the junction read will get an SA tag which will be assigned to two locations if and only if this is the one and unique alignment configuration the STAR software can find [37,38]. The read-mate aligns somewhere in between these two locations. Finding additional read pairs showing this configuration with a breakpoint at the exact same location strengthens the evidence for circular transcripts. We included only regions with at least five reads crossing the circular junction. After filtering, GENCODE annotation was used to obtain the exon locations of genes that exactly matched to the circular region. For each sample, STAR also gives the raw read counts for all genes. These were normalized (Trimmed Mean of M-values implemented in edgeR [39]), and the normalized read counts were used to correlate with the number of junction reads of the circular transcripts. The script is also available at https://bitbucket.org/snippets/MSmid/Le949d/identify-circularrna-reads (accessed on 30 October 2018).

2.2.4. Multiple Correspondence Analysis (MCA)

For a substantial number of genes, only a linear transcript is detected in the majority of samples, which results in many missing values per circRNA. This in turn, complicates the use of standard cluster analysis for the identification of sample groups with similar circRNA-related biology. Therefore, circRNA data were considered categorical, i.e., a circRNA was scored as either “present” or “absent” in a sample. These categorical data are suitable for a multiple correspondence analysis (MCA), which is a generalized principle component analysis. An MCA generates a combined plot that shows both patients and circRNAs in such a way that patients and circRNAs that have similar patterns are closer together. Thus, the colon cancer tumor samples and circRNAs are projected onto the same plane, in which the relative distance to either the samples or the circRNAs is meaningful. The 0,0 point corresponds to a sample or circRNA with an average profile. The R-package “ade4” was used to perform the MCA in R version 3.4.1. Custom functions to plot the MCA results are available upon request of the authors.

2.2.5. Reverse Transcription, Quantitative PCR, and Sanger Sequencing

Candidate circular RNAs were selected, and divergent primers were designed that are only able to amplify and detect the circular and not the corresponding linear mRNA (Table S1). Total RNA, isolated with RNA-Bee according to the manufacturer’s instructions (CS105B, TEL TEST), was reverse transcribed into cDNA with the H-minus RevertAid First Strand cDNA Synthesis Kit (K1632, ThermoFisher Scientific, Waltham, MA, USA), followed by an RNAse-H step (AM2293; Ambion). For Sanger sequencing, cDNA from five individual patient samples were used for PCR. cDNA was amplified for 35 cycles using Phusion high-fidelity DNA polymerase (ThermoFisher Scientific) in a total reaction volume of 25 μL, containing 400 nM of each primer and 160 µM dNTPs. For every circRNA two resulting PCR amplicons were purified from gel using the Qiaquick Gel Extraction Kit from Qiagen (Hilden, Germany) according to the manufacturer’s protocol and subjected to Quick Shot Sanger Sequencing by BaseClear BV (Leiden, The Netherlands).

2.3. Statistical Analyses

As indicated above for a substantial number of genes, only a linear transcript was detected in the majority of samples, which results in many missing values per circRNA. This hampers statistical analyses of circRNA expression levels and therefore we categorized the circRNA data into “present” or “absent” for statistical evaluation as well. We used circRNAs present in at least 20 samples to ensure a sufficient number of events for subsequent statistical analyses. STATA version 14 and SPSS Version 24.0 (SPSS, Inc., Chicago, IL, USA) were used to perform the statistical tests that are also indicated in the text. Cox’s proportional-hazards regression was used to evaluate the (log-transformed) number of uniquely present circRNAs per sample, hereafter called circRNA diversity, with DFS, or as “present”/“absent” when evaluating individual circRNAs. Survival curves were evaluated using the logrank test (for individual circRNAs) or with the logrank test for trend (for circRNA diversity, after dividing into three equal quantiles). Pearson’s correlation was used to correlate the circRNA expression with the expression of the linear gene it was derived from. Analyses between categorical variables (like present/absent of a circRNA versus MSI yes/no) were analyzed using Fisher’s exact test. Reported p-values are two-sided and considered significant at p ≤ 0.05. p-values were corrected for multiple testing using Benjamini–Hochberg’s FDR correction when evaluating multiple circRNAs, which were considered significant at p < 0.10.

3. Results

3.1. CircRNA Expression in Colon Cancer

We analyzed RNAseq data of 181 patients with chemonaive, stage I/II primary colon cancer. The median follow-up time was 53 months (IQR 37–59). Clinical and histopathological characteristics are listed in Table 1.

Table 1

Clinical and histopathological characteristics.

Clinical VariablesCategoriesn = 181%
GenderFemale92(50.9)
Male89(49.2)
Age (median, IQR)70(63–76)
Tumor stageStage I66(36.5)
Stage II115(63.5)
T statusT266(36.5)
T3110(60.8)
T45(2.8)
Nodal statusN0 ≥ 10 nodes assessed149(82.2)
N0 < 10 nodes assessed32(17.3)
Tumor gradeGood16(8.8)
Poor10(5.5)
Moderate152(84)
Unknown3(1.7)
LocationRight92(50.8)
Left89(49.2)
MSI statusMSI44(24.3)
MSS137(75.7)
RelapseNo152(84)
Yes29(16)

Abbreviations: MSI = Microsatellite instability, MSS = Microsatellite stable.

Circular RNAs were defined as present when at least five reads crossed the circular junction [36]. This resulted in the identification of 2606 distinct circRNAs in the entire cohort, of which 1860 were derived from known genes. Sixty-three percent of these were repeatedly occurring in at least two colon cancer samples (n = 1172) (Figure 1), whereas 277 (15%) were observed in 20 samples or more (Table S2). The most repeatedly occurring circRNAs were derived from SMARCA5, HIPK3, ZKSCAN1 and FBXW7 and were observed in 177 samples each (n.b. not the same 177 samples for all four circRNAs) (Table 2). For 29 genes we observed that more than one unique circRNA was derived from the linear sequence (Table S2). Relative to the total number of samples expressing at least one circRNA from the respective gene, varying levels of co-expression between circRNAs derived from the same gene were observed with a median of 26.19% (range: 4.40–82.30%).

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g001.jpg

Histogram showing the distribution of circRNA occurrence in the 181 stage I/II colon cancer samples. A total of 1860 distinct circRNAs were identified which were derived from known genes.

Table 2

Most frequently recurring circRNAs.

Circular RegionEnsembl Gene IDGeneExonsNr. of SamplesR *
chr4:143543509-143543973ENSG00000153147SMARCA515-161770.296
chr11:33286413-33287512ENSG00000110422HIPK321770.541
chr7:100023419-100024308ENSG00000106261ZKSCAN12-31770.539
chr4:152411303-152412530ENSG00000109670FBXW721770.437

* R indicates the Pearson correlation between the number of circRNA reads and mRNA reads for that gene.

We correlated the number of circRNA reads per circRNA with the expression of the linear gene from which the circRNA was derived. To avoid possible spurious correlations, only the 277 circRNAs found in at least 20 samples were analyzed. The vast majority of circRNAs showed a positive correlation with the linear gene from which the circRNA is derived (Figure 2). However, this correlation was poor (R < 0.3) for 126 circRNAs. Twenty-seven circRNAs showed a negative correlation with their corresponding linear gene, suggesting the circRNA may function independently from the linear transcript.

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g002.jpg

Histogram showing the distribution of the observed correlations between the number of circRNA and mRNA reads per gene.

We randomly selected four circRNAs representative of the entire list of identified circRNAs to get an unbiased validation of our circRNA identification pipeline. Our four candidates include two from the top-20 circRNAs showing the most positive and negative correlations to their linear counterparts respectively (SATB2_chr2:199368605-199433515, R = 0.61 and PLEKHM3_chr2:207976651-207977587, R = −0.17), and one novel circRNA (FGD6 (chr12:95208843-95211268), R = 0.39) not present in circBase (circbase.org (accessed on 23 May 2019) [40]). The identified junctions in the RNAseq data were verified for all four selected circRNAs by Sanger sequencing, thereby demonstrating the validity of our circRNA identification algorithm [36] (Figure 3).

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g003.jpg

Sanger sequencing was used to validate the circular exon junctions of circSATB2, circKMT2C, circFGD6, and circPLEKHM3 identified from the RNAseq data.

3.2. CircRNA Expression Patterns Are Associated with Relevant Clinical Factors

For a substantial number of genes, only a linear transcript is detected in the majority of samples, which results in many missing values per circRNA. This in turn complicates the use of standard cluster analysis for the identification of circRNA/sample groups with similar circRNA-related biology. To be able to investigate circRNA profiles, we categorized circRNAs as “present” or “absent” in a sample and used this in a multiple correspondence analysis (MCA) to find naturally occurring subgroups. An MCA plot projects the colon tumor samples and circRNAs onto the same plane, in which the relative distance to either the samples or the circRNAs is meaningful. As such, samples that group close together have more similar circRNA profiles. In addition, since circRNAs have two states (present/absent), both these states are used in the analysis. Thus, two circRNAs that are “present” frequently in the same samples (co-occurrence) will be placed at a short distance, but this is also true for circRNAs that are mutually exclusive (presence of a circRNA and absence of the other circRNA) across the samples. Coloring the circRNA states will reveal the co-occurrence/mutual-exclusivity.

For the MCA analysis, we used the 277 repeatedly occurring circRNAs (i.e., those which are present in at least 20 samples) and labelled these in each sample as “present” or “absent” as defined above. After MCA analysis, we first colored genes based on the circRNA state (Figure 4a). As shown by the clear separation of the “present” and “absent” state, circRNAs do not show mutual-exclusivity (which would show up as a red triangle among blue triangles or vice versa), but rather are often co-expressed in the same samples. Furthermore, the variability (spread among the x-axis) of the “present” profiles indicates different circRNA expression profiles are present among the samples, or, in other words, samples show a large diversity of expressed circRNAs.

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g004.jpg

MCA analysis plots in which samples (closed circles) and circRNAs (triangles) are projected onto the same plane. (a). Blue and red indicate circRNAs without (absent) or with (present) circRNA expression, where similarity of two circRNA expression profiles (either both present (both red), both absent (both blue), or mutually exclusive (one red and one blue)) over the samples results in a small relative distance between these circRNAs. (b). Samples are colored based on relapse status, red indicates patients who relapse (1) showing circRNA profiles that are close to the “circRNA absent” group. (c). Samples are colored based on tumor stage I (blue) or stage II (red) showing no clear distinction in circRNA profiles. (d). The consensus molecular subtype (CMS) of the samples are indicated. CMS3 samples are most closely located to the “circRNA present” group. (e). Samples colored based on tumor side (left = 0; in blue, right = 1; in red) and (f), based on microsatellite instability (MSI) (MSI = 1; in red, MSS = 0; in blue).

Next, we colored samples according to relapse status (Figure 4b), tumor stage (Figure 4c), CMS (Figure 4d), tumor localization (Figure 4e) and MSI (Figure 4f). Patients showing relapse (1) (Figure 4b) have profiles that are close to the “absent” group (group without circRNAs), which indicates that few genes give rise to circRNA expression in these samples. Indeed, when analyzing circRNA diversity (the number of distinct circRNA molecules in a sample) we found that a high diversity in circRNAs is associated with a favorable DFS: Cox regression using circRNA diversity as (log-transformed) continuous variable: Hazard Ratio (HR) 0.60, 95% CI 0.38–0.97, p = 0.036. Figure 5 shows Kaplan–Meier curves in which the levels of diversity of circRNAs were split into three equal quantiles to visualize the association between circRNA-diversity and DFS. The difference in DFS between the three quantiles was evaluated using the logrank test for trend, to account for the ordered structure of the sample groups (high, intermediate and low circRNA diversity; p = 0.050). High diversity was not associated with other clinical factors such as tumor stage, tumor side, MSI status, or CMS (diversity as continuous variable).

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g005.jpg

Kaplan–Meier survival curves of disease free survival in which patients were grouped in 3 equal quantiles based on their diversity in circRNA expression. The red line represents the quantile with the lowest diversity in circRNAs, the blue line represents the quantile with intermediate diversity in circRNAs and the green line represents the quantile with the highest diversity in circRNAs.

The MCA plot of tumor side (Figure 4e) shows that right-sided tumors are closer to the absent group (group of samples without circRNAs) than the left-sided tumors –therefore also closer to the relapse group, but this association was not significant. Combining CMS grouping (Figure 4d) and MSI (Figure 4f) shows that, as expected, samples that are CMS1 and those that are MSI tumors have a similar position. Combining CMS (Figure 4d) and circRNA diversity (Figure 4a) leads to the conclusion that CMS3 patients have the highest diversity in circRNAs, and CMS2 patients the lowest. As to tumor stage, there was no clear distinction between stage I and II tumors with regard to circRNAs (Figure 4e).

Next to this global analysis of overall circRNA profiles in the samples, we also investigated whether the presence/absence of specific circRNAs was associated with relapse status, tumor stage, CMS, tumor localization, and MSI. Whereas no specific circRNAs were significantly associated with tumor stage and localization, the presence/absence of nine circRNAs was associated with CMS and five with MSI (Table 3; Fisher exact test p < 0.0003; Benjamini-Hochberg corrected p-value < 0.1). circSATB2 (Special AT-rich sequence-binding protein 2), circNAP1 (Nucleosome assembly protein) and circCEP192 (Centrosomal Protein 192) each correlated with both MSI and CMS. Only absence of circMGA (MAX dimerization protein) was significantly associated with relapse (Table 3; Fisher exact test p = 0.0002; Benjamini–Hochberg corrected p-value = 0.06). Kaplan-Meier analysis showed that patients in whom circMGA was detected (n = 94) had a favorable DFS compared to patients in which circMGA was not detected (n = 87; log-rank test p < 0.001, Cox HR 0.22 95%CI 0.09–0.53, p < 0.001) (Figure 6).

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g006.jpg

Kaplan–Meier analysis of patients in whom circMGA was detected (n = 94) versus patients in whom circMGA was not detected (n = 87).

Table 3

Association of presence/absence of circRNAs with relapse status, tumor stage, CMS and MSI.

NameEnsemble Gene IDCircular RegionFisher p-ValueFDR *Comparison
circSATB2ENSG00000119042chr2:199368605-1994335152.47 × 10−86.84 × 10−6MSI
2.41 × 10−86.68 × 10−6CMS
circNAB1ENSG00000138386chr2:190659158-1906731536.48 × 10−70.000179MSI
0.000075840.02078CMS
circZBTB44ENSG00000196323chr11:130260856-1302619303.4 × 10−50.009347MSI
circCEP192ENSG00000101639chr18:12999421-130192070.0001520.041758MSI
0.000078080.021316CMS
chr18:12999421-130306090.0002330.063043CMS
circUBAP2ENSG00000137073chr9:33960826-339732380.0002380.064846MSI
circMGAENSG00000174197chr15:41668828-416699590.0002350.064979DFSI
circASPHENSG00000198363chr8:61618978-616536611.19 × 10−73.28 × 10−5CMS
circCCSER2ENSG00000107771chr10:84438512-844776653.9 × 10−50.010728CMS
circZNF609ENSG00000180357chr15:64499293-645001670.000170.046128CMS
circFUT8ENSG00000033170chr14:65561337-655617670.000260.070135CMS
circMRPS35ENSG00000061794chr12:27714780-277241870.0002990.08035CMS

* False discovery rate (Benjamini–Hochberg procedure). Abbreviations: MSI = microsatellite instability, CMS = consensus molecular subtype.

3.1. CircRNA Expression in Colon Cancer

We analyzed RNAseq data of 181 patients with chemonaive, stage I/II primary colon cancer. The median follow-up time was 53 months (IQR 37–59). Clinical and histopathological characteristics are listed in Table 1.

Table 1

Clinical and histopathological characteristics.

Clinical VariablesCategoriesn = 181%
GenderFemale92(50.9)
Male89(49.2)
Age (median, IQR)70(63–76)
Tumor stageStage I66(36.5)
Stage II115(63.5)
T statusT266(36.5)
T3110(60.8)
T45(2.8)
Nodal statusN0 ≥ 10 nodes assessed149(82.2)
N0 < 10 nodes assessed32(17.3)
Tumor gradeGood16(8.8)
Poor10(5.5)
Moderate152(84)
Unknown3(1.7)
LocationRight92(50.8)
Left89(49.2)
MSI statusMSI44(24.3)
MSS137(75.7)
RelapseNo152(84)
Yes29(16)

Abbreviations: MSI = Microsatellite instability, MSS = Microsatellite stable.

Circular RNAs were defined as present when at least five reads crossed the circular junction [36]. This resulted in the identification of 2606 distinct circRNAs in the entire cohort, of which 1860 were derived from known genes. Sixty-three percent of these were repeatedly occurring in at least two colon cancer samples (n = 1172) (Figure 1), whereas 277 (15%) were observed in 20 samples or more (Table S2). The most repeatedly occurring circRNAs were derived from SMARCA5, HIPK3, ZKSCAN1 and FBXW7 and were observed in 177 samples each (n.b. not the same 177 samples for all four circRNAs) (Table 2). For 29 genes we observed that more than one unique circRNA was derived from the linear sequence (Table S2). Relative to the total number of samples expressing at least one circRNA from the respective gene, varying levels of co-expression between circRNAs derived from the same gene were observed with a median of 26.19% (range: 4.40–82.30%).

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g001.jpg

Histogram showing the distribution of circRNA occurrence in the 181 stage I/II colon cancer samples. A total of 1860 distinct circRNAs were identified which were derived from known genes.

Table 2

Most frequently recurring circRNAs.

Circular RegionEnsembl Gene IDGeneExonsNr. of SamplesR *
chr4:143543509-143543973ENSG00000153147SMARCA515-161770.296
chr11:33286413-33287512ENSG00000110422HIPK321770.541
chr7:100023419-100024308ENSG00000106261ZKSCAN12-31770.539
chr4:152411303-152412530ENSG00000109670FBXW721770.437

* R indicates the Pearson correlation between the number of circRNA reads and mRNA reads for that gene.

We correlated the number of circRNA reads per circRNA with the expression of the linear gene from which the circRNA was derived. To avoid possible spurious correlations, only the 277 circRNAs found in at least 20 samples were analyzed. The vast majority of circRNAs showed a positive correlation with the linear gene from which the circRNA is derived (Figure 2). However, this correlation was poor (R < 0.3) for 126 circRNAs. Twenty-seven circRNAs showed a negative correlation with their corresponding linear gene, suggesting the circRNA may function independently from the linear transcript.

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g002.jpg

Histogram showing the distribution of the observed correlations between the number of circRNA and mRNA reads per gene.

We randomly selected four circRNAs representative of the entire list of identified circRNAs to get an unbiased validation of our circRNA identification pipeline. Our four candidates include two from the top-20 circRNAs showing the most positive and negative correlations to their linear counterparts respectively (SATB2_chr2:199368605-199433515, R = 0.61 and PLEKHM3_chr2:207976651-207977587, R = −0.17), and one novel circRNA (FGD6 (chr12:95208843-95211268), R = 0.39) not present in circBase (circbase.org (accessed on 23 May 2019) [40]). The identified junctions in the RNAseq data were verified for all four selected circRNAs by Sanger sequencing, thereby demonstrating the validity of our circRNA identification algorithm [36] (Figure 3).

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g003.jpg

Sanger sequencing was used to validate the circular exon junctions of circSATB2, circKMT2C, circFGD6, and circPLEKHM3 identified from the RNAseq data.

3.2. CircRNA Expression Patterns Are Associated with Relevant Clinical Factors

For a substantial number of genes, only a linear transcript is detected in the majority of samples, which results in many missing values per circRNA. This in turn complicates the use of standard cluster analysis for the identification of circRNA/sample groups with similar circRNA-related biology. To be able to investigate circRNA profiles, we categorized circRNAs as “present” or “absent” in a sample and used this in a multiple correspondence analysis (MCA) to find naturally occurring subgroups. An MCA plot projects the colon tumor samples and circRNAs onto the same plane, in which the relative distance to either the samples or the circRNAs is meaningful. As such, samples that group close together have more similar circRNA profiles. In addition, since circRNAs have two states (present/absent), both these states are used in the analysis. Thus, two circRNAs that are “present” frequently in the same samples (co-occurrence) will be placed at a short distance, but this is also true for circRNAs that are mutually exclusive (presence of a circRNA and absence of the other circRNA) across the samples. Coloring the circRNA states will reveal the co-occurrence/mutual-exclusivity.

For the MCA analysis, we used the 277 repeatedly occurring circRNAs (i.e., those which are present in at least 20 samples) and labelled these in each sample as “present” or “absent” as defined above. After MCA analysis, we first colored genes based on the circRNA state (Figure 4a). As shown by the clear separation of the “present” and “absent” state, circRNAs do not show mutual-exclusivity (which would show up as a red triangle among blue triangles or vice versa), but rather are often co-expressed in the same samples. Furthermore, the variability (spread among the x-axis) of the “present” profiles indicates different circRNA expression profiles are present among the samples, or, in other words, samples show a large diversity of expressed circRNAs.

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g004.jpg

MCA analysis plots in which samples (closed circles) and circRNAs (triangles) are projected onto the same plane. (a). Blue and red indicate circRNAs without (absent) or with (present) circRNA expression, where similarity of two circRNA expression profiles (either both present (both red), both absent (both blue), or mutually exclusive (one red and one blue)) over the samples results in a small relative distance between these circRNAs. (b). Samples are colored based on relapse status, red indicates patients who relapse (1) showing circRNA profiles that are close to the “circRNA absent” group. (c). Samples are colored based on tumor stage I (blue) or stage II (red) showing no clear distinction in circRNA profiles. (d). The consensus molecular subtype (CMS) of the samples are indicated. CMS3 samples are most closely located to the “circRNA present” group. (e). Samples colored based on tumor side (left = 0; in blue, right = 1; in red) and (f), based on microsatellite instability (MSI) (MSI = 1; in red, MSS = 0; in blue).

Next, we colored samples according to relapse status (Figure 4b), tumor stage (Figure 4c), CMS (Figure 4d), tumor localization (Figure 4e) and MSI (Figure 4f). Patients showing relapse (1) (Figure 4b) have profiles that are close to the “absent” group (group without circRNAs), which indicates that few genes give rise to circRNA expression in these samples. Indeed, when analyzing circRNA diversity (the number of distinct circRNA molecules in a sample) we found that a high diversity in circRNAs is associated with a favorable DFS: Cox regression using circRNA diversity as (log-transformed) continuous variable: Hazard Ratio (HR) 0.60, 95% CI 0.38–0.97, p = 0.036. Figure 5 shows Kaplan–Meier curves in which the levels of diversity of circRNAs were split into three equal quantiles to visualize the association between circRNA-diversity and DFS. The difference in DFS between the three quantiles was evaluated using the logrank test for trend, to account for the ordered structure of the sample groups (high, intermediate and low circRNA diversity; p = 0.050). High diversity was not associated with other clinical factors such as tumor stage, tumor side, MSI status, or CMS (diversity as continuous variable).

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g005.jpg

Kaplan–Meier survival curves of disease free survival in which patients were grouped in 3 equal quantiles based on their diversity in circRNA expression. The red line represents the quantile with the lowest diversity in circRNAs, the blue line represents the quantile with intermediate diversity in circRNAs and the green line represents the quantile with the highest diversity in circRNAs.

The MCA plot of tumor side (Figure 4e) shows that right-sided tumors are closer to the absent group (group of samples without circRNAs) than the left-sided tumors –therefore also closer to the relapse group, but this association was not significant. Combining CMS grouping (Figure 4d) and MSI (Figure 4f) shows that, as expected, samples that are CMS1 and those that are MSI tumors have a similar position. Combining CMS (Figure 4d) and circRNA diversity (Figure 4a) leads to the conclusion that CMS3 patients have the highest diversity in circRNAs, and CMS2 patients the lowest. As to tumor stage, there was no clear distinction between stage I and II tumors with regard to circRNAs (Figure 4e).

Next to this global analysis of overall circRNA profiles in the samples, we also investigated whether the presence/absence of specific circRNAs was associated with relapse status, tumor stage, CMS, tumor localization, and MSI. Whereas no specific circRNAs were significantly associated with tumor stage and localization, the presence/absence of nine circRNAs was associated with CMS and five with MSI (Table 3; Fisher exact test p < 0.0003; Benjamini-Hochberg corrected p-value < 0.1). circSATB2 (Special AT-rich sequence-binding protein 2), circNAP1 (Nucleosome assembly protein) and circCEP192 (Centrosomal Protein 192) each correlated with both MSI and CMS. Only absence of circMGA (MAX dimerization protein) was significantly associated with relapse (Table 3; Fisher exact test p = 0.0002; Benjamini–Hochberg corrected p-value = 0.06). Kaplan-Meier analysis showed that patients in whom circMGA was detected (n = 94) had a favorable DFS compared to patients in which circMGA was not detected (n = 87; log-rank test p < 0.001, Cox HR 0.22 95%CI 0.09–0.53, p < 0.001) (Figure 6).

An external file that holds a picture, illustration, etc.
Object name is cancers-13-01903-g006.jpg

Kaplan–Meier analysis of patients in whom circMGA was detected (n = 94) versus patients in whom circMGA was not detected (n = 87).

Table 3

Association of presence/absence of circRNAs with relapse status, tumor stage, CMS and MSI.

NameEnsemble Gene IDCircular RegionFisher p-ValueFDR *Comparison
circSATB2ENSG00000119042chr2:199368605-1994335152.47 × 10−86.84 × 10−6MSI
2.41 × 10−86.68 × 10−6CMS
circNAB1ENSG00000138386chr2:190659158-1906731536.48 × 10−70.000179MSI
0.000075840.02078CMS
circZBTB44ENSG00000196323chr11:130260856-1302619303.4 × 10−50.009347MSI
circCEP192ENSG00000101639chr18:12999421-130192070.0001520.041758MSI
0.000078080.021316CMS
chr18:12999421-130306090.0002330.063043CMS
circUBAP2ENSG00000137073chr9:33960826-339732380.0002380.064846MSI
circMGAENSG00000174197chr15:41668828-416699590.0002350.064979DFSI
circASPHENSG00000198363chr8:61618978-616536611.19 × 10−73.28 × 10−5CMS
circCCSER2ENSG00000107771chr10:84438512-844776653.9 × 10−50.010728CMS
circZNF609ENSG00000180357chr15:64499293-645001670.000170.046128CMS
circFUT8ENSG00000033170chr14:65561337-655617670.000260.070135CMS
circMRPS35ENSG00000061794chr12:27714780-277241870.0002990.08035CMS

* False discovery rate (Benjamini–Hochberg procedure). Abbreviations: MSI = microsatellite instability, CMS = consensus molecular subtype.

4. Discussion

With the use of RNAseq data, we could establish the presence of a wide variety of circRNAs in chemonaive lymph node negative, stage I/II primary colon tumors. Previous studies have been limited by the small number of circRNAs screened, the small sample size and retrospective data. Our study, however, concerned 181 patients included in a prospective, multicenter cohort study, and is therefore, to our knowledge, the largest circRNA-based biomarker discovery study done in stage I/II colon cancer.

The four most repeatedly occurring circRNAs that we found (177/181 samples), circSMARCA5, circHIPK3, circFBXW7 and circZKSCAN1, have also been described as such in previous studies. circSMARCA5 was reported to be induced during epithelial-to-mesenchymal transition, which is an important mechanism during the metastatic process that has been associated with the pathogenesis of several cancers [26,41,42,43,44]. circHIPK3 has been described to promote CRC growth and metastasis by sponging miR-7 [45]. Furthermore, previous research in CRC cell lines showed that circFBXW7 is conducive in controlling the progression of CRC through NEK2, mTOR, and PTEN signaling pathways [37]. The correspondence of our finding with previous results clearly underlines the validity of our approach in identifying circRNAs. In addition, we performed Sanger sequencing to verify four randomly selected circRNAs (circSATB2, circKMT2C circFGD6, and circPLEKHM3) and successfully validated the identified circular junctions for all four circRNAs.

In the studied cohort of chemonaive lymph node negative colon cancer patients, a first highlight was the finding that high diversity of circRNAs present in colon cancer tissue was associated with favorable DFS. Vo et al showed that across different cancer types, total circRNA abundance was lower in cancer compared to normal tissue, suggesting that the reduction of circRNA generation could be associated with loss of cellular differentiation [46]. More specifically, presence of circMGA was significantly associated with a favorable DFS. Together, these findings support the idea that circRNAs might play a functional role in cancer metastasis [26]. Recent studies provide evidence for a tumor suppressive role for the gene MGA (MAX dimerization protein) in colorectal cancer [47]. In lung adenocarcinoma, the molecular function of MGA appears to be antagonistic to that of MYC. To our knowledge, this is the first study associating the circRNA emanating from this gene with colon cancer or any other malignancies.

A second highlight of this study is the association between circRNAs and distinct colorectal cancer subtypes. Presence/absence of nine and five circRNAs was significantly associated with CMS and MSI, respectively, of which circSATB2, circNAB1, circCEP192 were overlapping. Although we were unable to find a suitable publicly available RNAseq dataset to validate the associations we found between circRNAs and clinical parameters in our cohort of stage I–II colon cancer patients, a number of the circRNAs we found to be associated with distinct subtypes of colon cancer were described before in cancer. circSATB2 has been described to play a notable role in the progression of lung cancer by binding to miR-326 [48], which in turn is associated with CRC [49]. The association between CEP192, NAB1 and CRC or other cancers, has not been described in previous studies. A role in CRC was proven for circZNF609 (Zinc Finger Protein 609), which is down-regulated in CRC tissue and promotes apoptosis in CRC by upregulating p53 [50]. circUBAP2 (ubiquitin associated protein 2) facilitates CRC progression by sponging miR-199a to upregulate VEGFA which implies that circUBAP2 may be a potential therapeutic biomarker for CRC [51]. Furthermore, circZBTB44 (Zinc Finger and BTB Domain Containing 44) and circZNF609 are both upregulated in acute lymphoblastic leukemia [52] of which especially circZNF609 has a known oncogenic potential in multiple other cancers as well [53,54,55,56,57]. CircASPH (Aspartate Beta-Hydroxylase) expression is upregulated in lung adenocarcinoma [58] and, finally, circFUT8 (Fucosyltransferase 8) functions as a tumor suppressor in bladder cancer cells where low circFUT8 was associated with poor prognosis, high histological grade, and lymph node metastasis [59]. The largest strength of this study is its prospective, multicenter study design and that it is, to our knowledge, the largest circRNA-based biomarker discovery study performed in stage I/II colon cancer. However, as mentioned before, a limitation of this study is that we were unable to find a suitable publicly available dataset to validate the associations we found between circRNAs and clinical parameters. Furthermore, some of the subgroup analyses, such as MSI, resulted in rather small sample sizes in outcome, increasing the chance of type II errors.

5. Conclusions

In conclusion, this study generated a comprehensive catalog of circRNAs in colon cancer and demonstrated the potential biological and clinical relevance of circRNAs in patients with stage I–II colon cancer. We demonstrated high diversity in circRNAs is associated with favorable DFS. As such, circRNAs represent a promising addition to the biomarker repertoire for colon cancer.

Department of Surgery, Erasmus MC-University Medical Center Rotterdam, 3015 GD Rotterdam, The Netherlands; ln.cmsumsare@grebnednav.i (I.v.d.B.); ln.cmsumsare@kaarbdvhgrebeoc.r (R.R.J.C.v.d.B.); ln.cmsumsare@snamrezji.j (J.N.M.I.)
Department of Medical Oncology, Erasmus MC Cancer Institute, University Medical Center Rotterdam, 3015 GD Rotterdam, The Netherlands; ln.cmsumsare@dims.m (M.S.); ln.cmsumsare@dreewed.v (V.d.W.); ln.cmsumsare@snekeof.j (J.A.F.); ln.cmsumsare@snetram.j (J.W.M.M.)
Department of Pathology, Erasmus MC-University Medical Center Rotterdam, 3015 GD Rotterdam, The Netherlands; ln.cmsumsare@nezruednav.m.h.c
Correspondence: ln.cmsumsare@gnitliw.s
These authors contributed equally to this paper.
Received 2021 Mar 16; Accepted 2021 Apr 13.
Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Abstract

Simple Summary

Colon cancer (CC) is one of the most common types of cancer. Circular RNAs (circRNAs) appear to play an important role in tumor progression of CC. They are stably expressed in saliva, blood, and exosomes, potentially rendering them promising biomarkers for the diagnosis, prognosis, and treatment of CC. In this study we describe the identification of an extensive catalog of circRNAs in a large cohort of 181 chemonaive, stage I/II primary colon tumors and related circRNA expression to consensus molecular subtypes (CMS), microsatellite instability (MSI) status and clinical outcome. We observed that a high diversity in circRNAs was associated with favorable disease-free survival, and that several circRNAs were associated with MSI and CMS, demonstrating the potential clinical value of circRNAs in CC.

Abstract

Circular RNAs (circRNAs) appear important in tumor progression of colon cancer (CC). We identified an extensive catalog of circRNAs in 181 chemonaive stage I/II colon tumors, who underwent curative surgery between 2007 and 2014. We identified circRNAs from RNAseq data, investigated common biology related to circRNA expression, and studied the association between circRNAs and relapse status, tumor stage, consensus molecular subtypes (CMS), tumor localization and microsatellite instability (MSI). We identified 2606 unique circRNAs. 277 circRNAs (derived from 260 genes) were repeatedly occurring in at least 20 patients of which 153 showed a poor or even negative (R < 0.3) correlation with the expression level of their linear gene. The circular junctions for circSATB2, circFGD6, circKMT2C and circPLEKHM3 were validated by Sanger sequencing. Multiple correspondence analysis showed that circRNAs were often co-expressed and that high diversity in circRNAs was associated with favorable disease-free survival (DFS), which was confirmed by Cox regression analysis (Hazard Ratio (HR) 0.60, 95% CI 0.38–0.97, p = 0.036). Considering individual circRNAs, absence of circMGA was significantly associated with relapse, whereas circSATB2, circNAB1, and circCEP192 were associated with both MSI and CMS. This study represents a showcase of the potential clinical utility of circRNAs for prognostic stratification in patients with stage I–II colon cancer and demonstrated that high diversity in circRNAs is associated with favorable DFS.

Keywords: colon cancer, circular RNA, prognostic stratification, CMS
Abstract

Abbreviations: MSI = Microsatellite instability, MSS = Microsatellite stable.

* R indicates the Pearson correlation between the number of circRNA reads and mRNA reads for that gene.

* False discovery rate (Benjamini–Hochberg procedure). Abbreviations: MSI = microsatellite instability, CMS = consensus molecular subtype.

Click here for additional data file.(141K, zip)

Author Contributions

I.v.d.B.: conceptualization, formal analysis, investigation, writing—original draft, visualization. M.S.: conceptualization, methodology, formal analysis, investigation, writing—review and editing. R.R.J.C.v.d.B.: conceptualization, resources, writing—review and editing. C.H.M.v.D.: validation, investigation, writing—review and editing. V.d.W.: validation, investigation, writing—review and editing. J.A.F: conceptualization, resources. J.N.M.I.: conceptualization, resources, writing—review and editing, supervision. J.W.M.M.: conceptualization, resources, writing—review and editing, supervision. S.M.W.: conceptualization, methodology, validation, investigation, writing—original draft, visualization, supervision. All authors have read and agreed to the published version of the manuscript.

Author Contributions

Funding

This research received no external funding.

Funding

Institutional Review Board Statement

This research was funded by Erasmus MC IRB, grant number MEC-2007-088.

Institutional Review Board Statement

Informed Consent Statement

Patients had given informed consent on the storage and use of tissue samples and the collection of clinical data for research purposes.

Informed Consent Statement

Data Availability Statement

The data that support the findings of this study are available on request from the corresponding author.

Data Availability Statement

Conflicts of Interest

The authors declare no conflict of interest.

Conflicts of Interest

Footnotes

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Footnotes

References

  • 1. Cijfers over kanker [(accessed on 23 May 2019)]; Available online: .[PubMed]
  • 2. Sargent D.J., Patiyil S., Yothers G., Haller D.G., Gray R., Benedetti J., Buyse M., Labianca R., Seitz J.F., O’Callaghan C.J., et al End Points for Colon Cancer Adjuvant Trials: Observations and Recommendations Based on Individual Patient Data From 20,898 Patients Enrolled Onto 18 Randomized Trials From the ACCENT Group. J. Clin. Oncol. 2007;25:4569–4574. doi: 10.1200/JCO.2006.10.4323.] [[PubMed][Google Scholar]
  • 3. Wang Z., Gerstein M., Snyder MRNA-Seq: A revolutionary tool for transcriptomics. Nat. Rev. Genet. 2009;10:57–63. doi: 10.1038/nrg2484.] [[Google Scholar]
  • 4. Iyer M.K., Niknafs Y.S., Malik R., Singhal U., Sahu A., Hosono Y., Barrette T.R., Prensner J.R., Evans J.R., Zhao S., et al The landscape of long noncoding RNAs in the human transcriptome. Nat. Genet. 2015;47:199–208. doi: 10.1038/ng.3192.] [[Google Scholar]
  • 5. Muers M.RGenome-wide views of long non-coding rnas. Nat. Rev. Genet. 2011;12:742–743. doi: 10.1038/nrg3088.] [[PubMed][Google Scholar]
  • 6. Sanger H.L., Klotz G., Riesner D., Gross H.J., Kleinschmidt A.KViroids are single-stranded covalently closed circular RNA molecules existing as highly base-paired rod-like structures. Proc. Natl. Acad. Sci. USA. 1976;73:3852–3856. doi: 10.1073/pnas.73.11.3852.] [[Google Scholar]
  • 7. Yan B., Gu W., Yang Z., Gu Z., Yue X., Gu Q., Liu LDownregulation of a long noncoding rna-ncrupar contributes to tumor inhibition in colorectal cancer. Tumour Biol. 2014;35:11329–11335. doi: 10.1007/s13277-014-2465-0.] [[PubMed][Google Scholar]
  • 8. Cao Y., Lin M., Bu Y., Ling H., He Y., Huang C., Shen Y., Song B., Cao Dp53-inducible long non-coding RNA PICART1 mediates cancer cell proliferation and migration. Int. J. Oncol. 2017;50:1671–1682. doi: 10.3892/ijo.2017.3918.] [[PubMed][Google Scholar]
  • 9. Xu J., Zhang R., Zhao JThe Novel Long Noncoding RNA TUSC7 Inhibits Proliferation by Sponging MiR-211 in Colorectal Cancer. Cell. Physiol. Biochem. 2017;41:635–644. doi: 10.1159/000457938.] [[PubMed][Google Scholar]
  • 10. Jacob H., Stanisavljevic L., Storli K.E., Hestetun K.E., Dahl O., Myklebust M.PA four-microrna classifier as a novel prognostic marker for tumor recurrence in stage ii colon cancer. Sci. Rep. 2018;8:6157. doi: 10.1038/s41598-018-24519-4.] [[Google Scholar]
  • 11. Bullock M.D., Pickard K., Mitter R., Sayan A.E., Primrose J.N., Ivan C., Calin G.A., Thomas G.J., Packham G.K., Mirnezami A.HStratifying risk of recurrence in stage II colorectal cancer using deregulated stromal and epithelial microRNAs. Oncotarget. 2015;6:7262–7279. doi: 10.18632/oncotarget.3225.] [[Google Scholar]
  • 12. Grassi A., Perilli L., Albertoni L., Tessarollo S., Mescoli C., Urso E.D.L., Fassan M., Rugge M., Zanovello PA coordinate deregulation of microRNAs expressed in mucosa adjacent to tumor predicts relapse after resection in localized colon cancer. Mol. Cancer. 2018;17 doi: 10.1186/s12943-018-0770-8.] [[Google Scholar]
  • 13. Su M., Xiao Y., Ma J., Tang Y., Tian B., Zhang Y., Li X., Wu Z., Yang D., Zhou Y., et al Circular RNAs in Cancer: Emerging functions in hallmarks, stemness, resistance and roles as potential biomarkers. Mol. Cancer. 2019;18:1–17. doi: 10.1186/s12943-019-1002-6.] [[Google Scholar]
  • 14. Memczak S., Jens M., Elefsinioti A., Torti F., Krueger J., Rybak A., Maier L., Mackowiak S.D., Gregersen L.H., Munschauer M., et al Circular RNAs are a large class of animal RNAs with regulatory potency. Nat. Cell Biol. 2013;495:333–338. doi: 10.1038/nature11928.] [[PubMed][Google Scholar]
  • 15. Zhang Y., Zhang X.-O., Chen T., Xiang J.-F., Yin Q.-F., Xing Y.-H., Zhu S., Yang L., Chen L.-LCircular Intronic Long Noncoding RNAs. Mol. Cell. 2013;51:792–806. doi: 10.1016/j.molcel.2013.08.017.] [[PubMed][Google Scholar]
  • 16. Ashwal-Fluss R., Meyer M., Pamudurti N.R., Ivanov A., Bartok O., Hanan M., Evantal N., Memczak S., Rajewsky N., Kadener SCircrna biogenesis competes with pre-mrna splicing. Mol. Cell. 2014;56:55–66. doi: 10.1016/j.molcel.2014.08.019.] [[PubMed][Google Scholar]
  • 17. Hansen T.B., Jensen T.I., Clausen B.H., Bramsen J.B., Finsen B., Damgaard C.K., Kjems JNatural RNA circles function as efficient microRNA sponges. Nat. Cell Biol. 2013;495:384–388. doi: 10.1038/nature11993.] [[PubMed][Google Scholar]
  • 18. Qu S., Yang X., Li X., Wang J., Gao Y., Shang R., Sun W., Dou K., Li HCircular RNA: A new star of noncoding RNAs. Cancer Lett. 2015;365:141–148. doi: 10.1016/j.canlet.2015.06.003.] [[PubMed][Google Scholar]
  • 19. Zhang LCircular RNA: The main regulator of energy metabolic reprogramming in cancer cells. Thorac. Cancer. 2019;11:6–7. doi: 10.1111/1759-7714.13251.] [[Google Scholar]
  • 20. Artemaki P.I., Scorilas A., Kontos C.KCircular RNAs: A New Piece in the Colorectal Cancer Puzzle. Cancers. 2020;12:2464. doi: 10.3390/cancers12092464.] [[Google Scholar]
  • 21. Hsu M.-T., Coca-Prados MElectron microscopic evidence for the circular form of RNA in the cytoplasm of eukaryotic cells. Nat. Cell Biol. 1979;280:339–340. doi: 10.1038/280339a0.] [[PubMed][Google Scholar]
  • 22. Wilusz J.E., Sharp P.AA Circuitous Route to Noncoding RNA. Science. 2013;340:440–441. doi: 10.1126/science.1238522.] [[Google Scholar]
  • 23. Cocquerelle C., Mascrez B., Hétuin D., Bailleul BMis-splicing yields circular RNA molecules. FASEB J. 1993;7:155–160. doi: 10.1096/fasebj.7.1.7678559.] [[PubMed][Google Scholar]
  • 24. Salzman J., Chen R.E., Olsen M.N., Wang P.L., Brown P.OCell-type specific features of circular rna expression. PLoS Genet. 2013;9:e1003777. doi: 10.1371/annotation/f782282b-eefa-4c8d-985c-b1484e845855.] [[Google Scholar]
  • 25. Jeck W.R., Sorrentino J.A., Wang K., Slevin M.K., Burd C.E., Liu J., Marzluff W.F., Sharpless N.ECircular RNAs are abundant, conserved, and associated with ALU repeats. RNA. 2012;19:141–157. doi: 10.1261/rna.035667.112.] [[Google Scholar]
  • 26. Zhang M., Xin YCircular RNAs: A new frontier for cancer diagnosis and therapy. J. Hematol. Oncol. 2018;11 doi: 10.1186/s13045-018-0569-5.] [[Google Scholar]
  • 27. Karousi P., Artemaki P.I., Sotiropoulou C.D., Christodoulou S., Scorilas A., Kontos C.KIdentification of Two Novel Circular RNAs Deriving from BCL2L12 and Investigation of Their Potential Value as a Molecular Signature in Colorectal Cancer. Int. J. Mol. Sci. 2020;21:8867. doi: 10.3390/ijms21228867.] [[Google Scholar]
  • 28. Wang S., Zhang K., Tan S., Xin J., Yuan Q., Xu H., Xu X., Liang Q., Christiani D.C., Wang M., et al Circular RNAs in body fluids as cancer biomarkers: The new frontier of liquid biopsies. Mol. Cancer. 2021;20:1–10. doi: 10.1186/s12943-020-01298-z.] [[Google Scholar]
  • 29. van Vugt J.L., Braak R.R.C.V.D., Lalmahomed Z.S., Vrijland W.W., Dekker J.W., Zimmerman D.D., Vles W.J., Coene P.-P.L., Ijzermans J.NImpact of low skeletal muscle mass and density on short and long-term outcome after resection of stage I-III colorectal cancer. Eur. J. Surg. Oncol. 2018;44:1354–1360. doi: 10.1016/j.ejso.2018.05.029.] [[PubMed][Google Scholar]
  • 30. Sieuwerts A.M., Lyng M.B., Gelder M.E.M.-V., De Weerd V., Sweep F.C., Foekens J.A., Span P.N., Martens J.W., Ditzel H.JEvaluation of the ability of adjuvant tamoxifen-benefit gene signatures to predict outcome of hormone-naive estrogen receptor-positive breast cancer patients treated with tamoxifen in the advanced setting. Mol. Oncol. 2014;8:1679–1689. doi: 10.1016/j.molonc.2014.07.003.] [[Google Scholar]
  • 31. Sieuwerts A.M., Gelder M.E.M.-V., Timmermans M., Trapman A.M., Garcia R.R., Arnold M., Goedheer A.J., Portengen H., Klijn J.G., Foekens J.AHow ADAM-9 and ADAM-11 Differentially from Estrogen Receptor Predict Response to Tamoxifen Treatment in Patients with Recurrent Breast Cancer: A Retrospective Study. Clin. Cancer Res. 2005;11:7311–7321. doi: 10.1158/1078-0432.CCR-05-0560.] [[PubMed][Google Scholar]
  • 32. Coebergh van den Braak R.R.J., Sieuwerts A.M., Kandimalla R., Lalmahomed Z.S., Bril S.I., Van Galen A., Smid M., Biermann K., Van Krieken J.H.J.M., Kloosterman W.P., et al High mRNA expression of splice variant SYK short correlates with hepatic disease progression in chemonaive lymph node negative colon cancer patients. PLoS ONE. 2017;12:e0185607. doi: 10.1371/journal.pone.0185607.] [[Google Scholar]
  • 33. Smid M., Coebergh van den Braak R.R.J., van de Werken H.J.G., van Riet J., van Galen A., de Weerd V., van der Vlugt-Daane M., Bril S.I., Lalmahomed Z.S., Kloosterman W.P., et al Gene length corrected trimmed mean of m-values (getmm) processing of rna-seq data performs similarly in intersample analyses while improving intrasample comparisons. BMC Bioinform. 2018;19 doi: 10.1186/s12859-018-2246-7.] [[Google Scholar]
  • 34. Bacher J.W., Flanagan L.A., Smalley R.L., Nassif N.A., Burgart L.J., Halberg R.B., Megid W.M.A., Thibodeau S.NDevelopment of a Fluorescent Multiplex Assay for Detection of MSI-High Tumors. Dis. Markers. 2004;20:237–250. doi: 10.1155/2004/136734.] [[Google Scholar]
  • 35. Guinney J., Dienstmann R., Wang X., De Reyniès A., Schlicker A., Soneson C., Marisa L., Roepman P., Nyamundanda G., Angelino P., et al The consensus molecular subtypes of colorectal cancer. Nat. Med. 2015;21:1350–1356. doi: 10.1038/nm.3967.] [[Google Scholar]
  • 36. Smid M., Wilting S.M., Uhr K., Rodríguez-González F.G., De Weerd V., Der Smissen W.J.P.-V., Van Der Vlugt-Daane M., Van Galen A., Nik-Zainal S., Butler A., et al The circular RNome of primary breast cancer. Genome Res. 2019;29:356–366. doi: 10.1101/gr.238121.118.] [[Google Scholar]
  • 37. Lu H., Yao B., Wen X., Jia BFBXW7 circular RNA regulates proliferation, migration and invasion of colorectal carcinoma through NEK2, mTOR, and PTEN signaling pathways in vitro and in vivo. BMC Cancer. 2019;19:1–8. doi: 10.1186/s12885-019-6028-z.] [[Google Scholar]
  • 38. Dobin A., Davis C.A., Schlesinger F., Drenkow J., Zaleski C., Jha S., Batut P., Chaisson M., Gingeras T.RSTAR: Ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21. doi: 10.1093/bioinformatics/bts635.] [[Google Scholar]
  • 39. Robinson M.D., Oshlack AA scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25. doi: 10.1186/gb-2010-11-3-r25.] [[Google Scholar]
  • 40. Glažar P., Papavasileiou P., Rajewsky NcircBase: A database for circular RNAs. RNA. 2014;20:1666–1670. doi: 10.1261/rna.043687.113.] [[Google Scholar]
  • 41. Gurzu S., Turdean S., Kovecsi A., Contac A.O., Jung IEpithelial-mesenchymal, mesenchymal-epithelial, and endothelial-mesenchymal transitions in malignant tumors: An update. World J. Clin. Cases. 2015;3:393–404. doi: 10.12998/wjcc.v3.i5.393.] [[Google Scholar]
  • 42. Cai J., Chen Z., Zuo XCircsmarca5 functions as a diagnostic and prognostic biomarker for gastric cancer. Dis. Markers. 2019;2019 doi: 10.1155/2019/2473652.] [[Google Scholar]
  • 43. Barbagallo D., Caponnetto A., Brex D., Mirabella F., Barbagallo C., Lauretta G., Morrone A., Certo F., Broggi G., Caltabiano R., et al CircSMARCA5 Regulates VEGFA mRNA Splicing and Angiogenesis in Glioblastoma Multiforme Through the Binding of SRSF. Cancers. 2019;11:194. doi: 10.3390/cancers11020194.] [[Google Scholar]
  • 44. Kong Z., Wan X., Zhang Y., Zhang P., Zhang Y., Zhang X., Qi X., Wu H., Huang J., Li YAndrogen-responsive circular RNA circSMARCA5 is up-regulated and promotes cell proliferation in prostate cancer. Biochem. Biophys. Res. Commun. 2017;493:1217–1223. doi: 10.1016/j.bbrc.2017.07.162.] [[PubMed][Google Scholar]
  • 45. Zeng K., Chen X., Xu M., Liu X., Hu X., Xu T., Sun H., Pan Y., He B., Wang SCircHIPK3 promotes colorectal cancer growth and metastasis by sponging miR. Cell Death Dis. 2018;9:1–15. doi: 10.1038/s41419-018-0454-8.] [[Google Scholar]
  • 46. Vo J.N., Cieslik M., Zhang Y., Shukla S., Xiao L., Zhang Y., Wu Y.-M., Dhanasekaran S.M., Engelke C.G., Cao X., et al The Landscape of Circular RNA in Cancer. Cell. 2019;176:869–881.e13. doi: 10.1016/j.cell.2018.12.021.] [[Google Scholar]
  • 47. Mathsyaraja H., Catchpole J., Eastwood E., Babaeva E., Geuenich M., Cheng P.F., Freie B., Ayers J., Yu M., Wu N., et al Loss of mga mediated polycomb repression promotes tumor progression and invasiveness. bioRxiv. 2020 doi: 10.1101/2020.10.16.334714.[PubMed][Google Scholar]
  • 48. Zhang N., Nan A., Chen L., Li X., Jia Y., Qiu M., Dai X., Zhou H., Zhu J., Zhang H., et al Circular RNA circSATB2 promotes progression of non-small cell lung cancer cells. Mol. Cancer. 2020;19:1–16. doi: 10.1186/s12943-020-01221-6.] [[Google Scholar]
  • 49. Wu L., Hui H., Wang L.-J., Wang H., Liu Q.-F., Han S.-XMicroRNA-326 functions as a tumor suppressor in colorectal cancer by targeting the nin one binding protein. Oncol. Rep. 2015;33:2309–2318. doi: 10.3892/or.2015.3840.] [[PubMed][Google Scholar]
  • 50. Zhang X., Zhao Y., Kong P., Han M., Li BExpression of circZNF609 is Down-Regulated in Colorectal Cancer Tissue and Promotes Apoptosis in Colorectal Cancer Cells by Upregulating p. Med. Sci. Monit. 2019;25:5977–5985. doi: 10.12659/MSM.915926.] [[Google Scholar]
  • 51. Dai J., Zhuang Y., Tang M., Qian Q., Chen J.PCircrna ubap2 facilitates the progression of colorectal cancer by regulating mir-199a/vegfa pathway. Eur. Rev. Med. Pharmacol. Sci. 2020;24:7963–7971.[PubMed][Google Scholar]
  • 52. Buratin A., Paganin M., Gaffo E., Molin A.D., Roels J., Germano G., Siddi M.T., Serafin V., De Decker M., Gachet S., et al Large-scale circular RNA deregulation in T-ALL: Unlocking unique ectopic expression of molecular subtypes. Blood Adv. 2020;4:5902–5914. doi: 10.1182/bloodadvances.2020002337.] [[Google Scholar]
  • 53. He Y., Huang H., Jin L., Zhang F., Zeng M., Wei L., Tang S., Chen D., Wang WCircZNF609 enhances hepatocellular carcinoma cell proliferation, metastasis, and stemness by activating the Hedgehog pathway through the regulation of miR-15a-5p/15b-5p and GLI2 expressions. Cell Death Dis. 2020;11:1–12. doi: 10.1038/s41419-020-2441-0.] [[Google Scholar]
  • 54. Li M., Li Y., Yu MCircRNA ZNF609 Knockdown Suppresses Cell Growth via Modulating miR-188/ELF2 Axis in Nasopharyngeal Carcinoma. Oncotargets. 2020;13:2399–2409. doi: 10.2147/OTT.S234230.] [[Google Scholar]
  • 55. Liu Z., Pan H.-M., Xin L., Zhang Y., Zhang W.-M., Cao P., Xu H.-WCirc-ZNF609 promotes carcinogenesis of gastric cancer cells by inhibiting miRNA-145-5p expression. Eur. Rev. Med. Pharm. Sci. 2019;23:9411–9417.[PubMed][Google Scholar]
  • 56. Zhu L., Liu Y., Yang Y., Mao X.-M., Yin Z.-DCircRNA ZNF609 promotes growth and metastasis of nasopharyngeal carcinoma by competing with microRNA-150-5p. Eur. Rev. Med. Pharm. Sci. 2019;23:2817–2826.[PubMed][Google Scholar]
  • 57. Rossi F., Legnini I., Megiorni F., Colantoni A., Santini T., Morlando M., Di Timoteo G., Dattilo D., Dominici C., Bozzoni ICirc-ZNF609 regulates G1-S progression in rhabdomyosarcoma. Oncogene. 2019;38:3843–3854. doi: 10.1038/s41388-019-0699-4.] [[Google Scholar]
  • 58. Xu L., Ma Y., Zhang H., Lu Q.-J., Yang L., Jiang G.-N., Liao W.-LHMGA2 regulates circular RNA ASPH to promote tumor growth in lung adenocarcinoma. Cell Death Dis. 2020;11:1–10. doi: 10.1038/s41419-020-2726-3.] [[Google Scholar]
  • 59. He Q., Yan D., Dong W., Bi J., Huang L., Yang M., Huang J., Qin H., Lin TcircRNA circFUT8 Upregulates Krüpple-like Factor 10 to Inhibit the Metastasis of Bladder Cancer via Sponging miR-570-3p. Mol. Oncol. 2020;16:172–187. doi: 10.1016/j.omto.2019.12.014.] [[Google Scholar]
Collaboration tool especially designed for Life Science professionals.Drag-and-drop any entity to your messages.