Simultaneous Study of Circular RNAs and Messenger RNAs in Colorectal Cancer: The Unbalanced Fate of a Couple?

1. Introduction

Known for decades [1,2], the messenger RNA (mRNA) is the centerpiece of the transcriptome. This intermediary between DNA and protein is regulated by various mechanisms, including transcription [3,4,5], splicing [6], stability [7], capping [8], transport [9] and translation [10]. Among the elements involved in this regulation, many types of RNA have been documented. MicroRNAs, but also pseudogene transcripts and long non-coding RNAs, have been shown to play a role in regulating mRNA [11]. These different elements are grouped together in a network known as competing endogenous RNA (ceRNA), and accumulating evidence points to its preponderant role in various pathologies such as cancer [12], neurodegenerative [13], cardiovascular [14] or autoimmune diseases [15]. Recently, with growing interest, circular RNAs (circRNAs) have been included in this ceRNA network [16].

circRNAs are single-stranded, closed-loop covalent structures, devoid of 5' caps nor 3' poly-A tails. They result from a process known as backsplicing, whereby a 5' splice site binds to an upstream 3' splice site [17]. Interest in circRNAs has increased greatly, mainly because the levels of certain circRNAs have proved important for better diagnosis, prognosis, disease monitoring in Alzheimer disease [18], cardiovascular disease [19] or cancer [20]. These observations are often attributed to their functions as microRNA and RNA binding protein sponges [21] but other mechanisms are emerging e.g. cap-independent translation that enables the generation of peptides involved in different tumorigenesis pathways [22,23] or transcription modulation. This last mechanism is particularly interesting, as it has been shown that circRNAs can regulate the levels of various mRNAs, including those of the parental gene [24,25] and a physiological balance between these two forms of transcript has been demonstrated [26,27]. Consequently, it may be hypothesized that a disruption of this physiological balance represents a novel, hidden mechanism in genetic diseases [28,29,30,31]. The question is particularly relevant in diseases where known mechanisms explain only a fraction of cases, the rest are referred to as “missing heritability”. Colorectal cancer (CRC) exemplifies this issue as the majority of non-polyposis familial and early onset microsatellite stable (MSS) cases have no identified genetic cause [32]. To investigate this hypothesis, we performed a simultaneous analysis of circRNA and mRNA expression in blood samples for 23 CRC predisposition genes in a cohort of 712 CRC patients (i.e. the “missing heritability” cohort) and 249 matched controls using a novel, dedicated technique named SEALigHTS (Splice and Expression Analyses by exon Ligation and High-Throughput Sequencing) [29,33]. We described the backsplicing landscape for these genes, studied the circRNA-mRNA balance between patients and controls and searched for a competitive mechanism between circRNA and mRNA expression.

2. Materials and Methods

Patients and controls

The collection of unexplained CRC cases was previously described [32] and includes 1029 patients with a personal and/or familial history suggestive of hereditary predisposition to CRC, Lynch syndrome and polyposis excluded (Table 1). All patients were recruited by the French network of Cancer Genetics Departments. In addition, 500 healthy volunteers without any personal or familial CRC history among their first-degree relatives were recruited by the Clinical Investigation Center of Rouen University Hospital. Among them, 712 patients and 500 controls had their blood sampled on PAXgene tubes for RNA extraction. To ensure quality, PAXgene tubes were transported within two days after sampling to the investigating center and extracted within 2 weeks. All of the 712 patients were included in the study. Among 500 healthy controls with available RNA samples, matched by age and sex to patients, 249 were randomly included. Written informed consent was obtained for all (ethics approval: DC-2013-1759).

Inclusion criteria

• CRC in two first-degree relatives, one being diagnosed before 61 years of age
• CRC diagnosed before 51 years of age or advanced colorectal adenoma (diameter over 1 cm, and/or tubulovillous or villous and/or with high-grade dysplasia) before 41 years of age
• Multiple primary CRCs in the same individual, the first one being diagnosed before 61 years of age

Exclusion criteria

• Lynch syndrome, as defined by the presence in the patient of a germline MMR gene pathogenic variant and/or a MSI tumour, in the context of a suggestive presentation (familial history, early age of onset)
• Adenomatous polyposis, as defined by more than 10 histologically proven adenomas
• Hamartomatous polyposis, as defined by the presence of histologically proven hamartomas

Table 1. Description of the criteria used to select patients for inclusion in the collection.

RNA extraction

RNA was extracted from peripheral blood with the PAXgene Blood RNA kit (PreAnalytiX, Switzerland) and concentration was determined with the NanoDrop 1000 Spectrophotometer (ThermoFisher Scientific, MA, USA), with a mean of 88 ng/µL (min: 4; max: 446).

SEALigHTS

Principle

Initially used as a multiplex technique for fusion transcript detection [34] and for measuring gene expression [35], SEALigHTS has been adapted and validated to study and quantify splicing and backsplicing [29]. Briefly, SEALigHTS allows for the simultaneous exploration of all exon-exon junctions on a panel of genes of interest, thanks to probes designed at exon extremities. Following reverse transcription and probe hybridization on cDNA, nearby probes are ligated if splicing and/or backsplicing occurs and the number of ligations is quantified using unique molecular identifiers and high throughput sequencing. All possible combinations of exons, i.e. splicing and backsplicing, are detectable.

Protocol

Oligonucleotides probes contained (i) specific sequences for each exon (19 – 30 bases length to obtain an optimal melting temperature of 70°C), (ii) Unique Molecular Identifiers (UMI), consisting of 7 random bases, to count the number of ligations (iii) complementary sequences of universal PCR primers. In accordance with the transcripts present in Ensembl and/or GTex Portal and literature, 788 probes (Supplementary Table S1) were designed at exon boundaries for the 23 CRC predisposition genes listed by the French experts from the Groupe Génétique et Cancer (GGC) [36] (Table 2). For exon ends where a single nucleotide polymorphism (SNP) is present with a frequency higher than 1% in the Caucasian population, the corresponding probe was designed with a combination of the 2 nucleotides involved to avoid any hybridization problems. If the alternative and canonical splice sites to be explored are close, partial probe overlap occurs. In this case, and when the sequence overlap was greater than 50%, the alternative probe was not selected. To control for DNA contamination, intronic probes were designed for each gene and added in the probe mix (Supplementary Table S1). Following quantification, 14−500 ng of total RNA were converted into cDNA using a SuperScript™ VILO™ cDNA Synthesis Kit (Invitrogen, Carlsbad, CA, USA). cDNAs were incubated for 1 h at 60 °C with our mix of 788 oligonucleotide probes (Supplementary Table S1) in 1× SALSA MLPA buffer (MRC Holland, Amsterdam, The Netherlands). Following hybridization, neighboring probes were ligated using the thermostable SALSA DNA ligase (MRC Holland, Amsterdam, The Netherlands).

Among the 23 genes studied (Table 2), the gene expression given by GTex Portal is up to 10,000 times higher for some genes than for others. After an initial trial, it appeared that PTEN and RPS20 high expression mandated a separate amplification mix in combination with the Q5® Hot Start High-Fidelity 2X Master Mix (NEB, Ipswich, MA, USA), which is why 2 independent PCRs were performed using different universal PCR primers. PCR products were purified using AMPure XP magnetic beads (Beckman Coulter, Brea, CA, USA). The library containing PTEN and RPS20 was diluted 1:20 in the library with the other genes. Sequencing of amplicons was carried out using a NextSeq® system with 75 cycles (Illumina, San Diego, CA, USA).

Sequencing reads, comprising a maximum of 75 bases, including UMIs, left and right probes, and barcodes, were demultiplexed using the barcodes and aligned with probe sequences. Exon junctions were counted for quantification purposes, thanks to the 7 random bases of the UMI, allowing 16 384 different combinations of unique molecules. A minimum of 65% mapping of left and right probes was deemed necessary to avoid primer dimers and ensure sufficient ligation and reliable results. A custom Python script, available on request, generated schematic backsplicing and splicing profiles (Supplementary Figure S1).

Bioinformatics and statistics analyses

Ratio circRNAs/mRNAs

Linear mRNAs (splicing) were distinguished from circRNAs (backsplicing) thanks to the order of the probes. If a probe located at the 3′ boundary of an exon is ligated with a probe located at the 5′ boundary of the following exon, the transcript is linear. Conversely, if this 3′ probe is ligated to a 5′ probe of a preceding exon, the transcript is circular (Supplementary Figure S1). Each circular RNA is a unique molecule, so all circular junctions must be considered, unlike an mRNA molecule made up of successive linear junctions. Hence, for each gene, the ratio between circRNAs and mRNAs was obtained by dividing the total number of UMIs counts of all circular junctions by the median number of UMIs counts of all canonical linear junctions of that gene i.e. ${\displaystyle \frac{UMIallcircularjunctions}{medianUMIallcanonicallinearjunctions}}\times 100$. Importantly, some of the UMIs attributed to linear junctions may actually stem from circRNAs. This occurs when a circRNA contains multiple exons, allowing probes to hybridize to internal exons. Therefore, the signal may be misinterpreted as linear junctions, impacting mRNA expression calculation. To mitigate this, only linear junction UMIs with no potential influence from circRNAs expressed at more than 10% of the corresponding mRNA were counted for mRNA calculation. A Student t-test was used to assess differences between patient and control groups, and p-values were corrected with a Bonferroni method to guarantee an overall control at 5% of type-I error rate.

mRNA and circRNAs expression

For each sample and junction, the UMI counts for each junction were normalized according to the library size i.e. ${\displaystyle \frac{UMIjunction}{UMIalljunctions}}=UMInormalizedjunction$. Hence, for each gene and sample, mRNA expression was calculated as the mean of normalized “circular-free” canonical linear junctions. In contrast, as circRNAs represented distinct transcripts, expression per gene was obtained by summing the UMIs of normalized circular junctions.

Circular RNAs characterization

To characterize the full circRNA sequences, divergent primers were designed (Supplementary Figure S2A and B). From 500ng of RNA, cDNA was generated with the SuperScript™ VILO™ cDNA Synthesis kit (Invitrogen, Carlsbad, CA, USA), followed by a PCR step with the Q5® Hot Start High-Fidelity 2X Master Mix kit (NEB, Ipswich, MA, USA). Gel-purified PCR products were sequenced using PCR primers and the BigDye v3 kit (Applied Biosystems, South San Francisco, CA, USA) and analyzed by capillary electrophoresis on a "3500 Genetic Analyzer" (Applied Biosystems, South San Francisco, CA, USA).

3. Results

3.1. General considerations

Of the 961 samples analyzed, 657 (68.4%) i.e. 433 patients and 224 controls, met the quality criteria (65% mapping) with a mean UMI count of 1,183,735 ensuring homogeneous samples and results, as shown by principal component analysis (PCA) (Supplementary Figure S3).

GREM1 and EPCAM genes are poorly expressed in whole blood. As the number of UMI is gene expression-dependent we were not able to detect all the canonical junctions and therefore GREM1 and EPCAM were excluded from further analyses.

For the 21 genes analyzed, 559 linear junctions i.e. canonical and alternatives were identified (Supplementary Table S1), with diversity in their numbers between genes (Figure 1 and Table 2). Although the number of linear junctions depends on the number of exons, the number of expressed alternative junctions differed between genes, e.g. APC and MSH2 both have 16 exons but 43 and 25 total linear junctions were expressed in blood, respectively. As the probes hybridize on both PMS2 and its pseudogene PMS2CL, their respective mRNAs and circRNAs cannot be distinguished. Therefore, PMS2 was not considered in the analyses.

3.2. Backsplicing landscape

Of the 20 genes studied (PMS2 excluded), SEALigHTS detected 220 circular junctions (Supplementary Table S1) of which 47 (21.36%) were novel (Table 2), not reported in RJunBase [37], a database compiling data on linear junctions, circular junctions, and fusion transcripts. A wide diversity in the number and level of expression of circRNAs (Table 2 and Figure 1) was found. Three genes (AXIN2, GALNT12 and NTHL1) did not show circular junctions. For the remaining 17, the median number of circular junctions per gene was 7 and ranged from 1 to 42 for RPS20 and MSH3, respectively. The number of circular junctions was not linked to the number of linear junctions e.g. MUTYH had 35 linear junctions and only 2 circular junctions, unlike SMAD4 which has 25 and 28 respectively (Table 2 and Figure 1). While circRNAs are expressed at low levels for most of the genes (17/20), POLD1, BUB1 and SMAD4 had higher or similar levels of circRNA compared to mRNA molecules. The relative expression of circRNA to mRNA varied between genes and was not related to the number of circRNAs produced by the gene. For example, POLD1 had the highest circRNA/mRNA expression level 547.28%, i.e. an relative abundance of 547 circRNAs for 100 mRNA molecules, with only 4 circRNAs produced, while APC had a ratio of 9.65% with 31 circRNAs produced (Table 2).

Following circRNA detection, we wanted to describe the entire sequence of the most frequent circRNAs i.e. POLD1_circRNA_3-2, BUB1_circRNA_19-10, BUB1_circRNA_9-2, BUB1_circRNA_9-6, SMAD4_circRNA_8-6, SMAD4_circRNA_10-5, SMAD4_circRNA_8-5 and found alternative backsplicing events in 2 cases. Sequencing the circRNA joining exons 3 to 2 of POLD1 (POLD1_circRNA_3-2) characterized a POLD1_circRNA_3-2p (Supplementary Figure S2C), i.e., with the use of an alternative splice site in exon 2, reported in RJunBase for both linear and circular transcripts. For the circRNA joining exons 10 to 5 of SMAD4 (SMAD4_circRNA_10-5), sequence analysis showed a full-length backsplice transcript and an alternative one with exons 6 and 7 skipping (Supplementary Figure S2D). As these alternative splicing events are also present in mRNAs, this suggests the use of the same (back)splicing machinery to produce messenger and circular RNAs to increase circRNAs diversity.

3.3. Ratio circRNA/mRNA between CRC patients and controls

The ratio was calculated for the 17 genes harboring circular junctions. We observed a significant 1.93-fold increase in the mean circRNA/mRNA ratio in patients compared to controls (433 vs 224; p<0.001), (Figure 2A). This result was confirmed independently on POLD1, BUB1 and SMAD4, which exhibited the highest circRNA/mRNA ratios. This increase was predominantly driven by POLD1 (Figure 2B), with a fold change of 3.11, rather than BUB1 and SMAD4 (fold changes 1.36 and 1.33, respectively) (Figure 2C and D). We identified 31 outlier patients i.e. with a cirRNA/mRNA ratio above mean +/-2 DS. Among these patients, 4 exhibited outlier ratios for all three genes (POLD1, BUB1, and SMAD4), 14 and 13 patients displayed outlier ratios for two and one of these genes, respectively. To explore potential genotype-phenotype correlations, we investigated whether this subgroup of 31 patients was associated with specific clinical phenotypes. No significant association was found between this subgroup and the presence of adenocarcinoma (p = 0.43; Fisher’s Exact test), adenoma (p = 0.75; Pearson's Chi-squared test), familial cancer (p = 0.51; Pearson's Chi-squared test), early sporadic cancer (p = 0.41; Fisher’s exact test), or multiple primary tumors (p = 0.64; Fisher’s exact test).

To shed light on a possible mechanism, we then investigated whether this circRNA/mRNA increase was linked to a drop or an increased level of mRNA or circRNA, respectively. Regression analysis between circRNA and mRNA expression showed that this imbalance probably reflected an increase in circRNA production or accumulation rather than a lower POLD1, BUB1 and SMAD4 expression, thereby suggesting an independent regulation (Figure 3).

4. Discussion

CircRNAs are now established as key players in human pathology, but original data remain scarce. We investigated whether disruption to the coregulation of mRNA–circRNA could represent a novel disease mechanism in a cohort of 712 unexplained CRC cases and 249 matched controls. While our data do not support a competitive regulatory mechanism, they reveal a significant imbalance in the mRNA–circRNA couple in blood between CRC patients and controls.

These results were obtained using SEALigHTS, a novel targeted high-throughput approach that, to our knowledge, is the only technology enabling the simultaneous quantitative analysis of mRNAs and circRNAs, including from FFPE-compatible material. By relying on probe hybridization at exon boundaries and UMI-based quantification prior to amplification, SEALigHTS circumvents the major limitations of conventional RNA-seq strategies for circRNA detection [38], without requiring RNA enrichment or specific sample treatment [39].

Our study provides new insights into i) the mechanism of backsplicing ii) the landscape of circRNAs for 20 CRC genes iii) the involvement of circRNAs from these 20 genes in CRC. These different aspects will be discussed successively. Linear splicing and backsplicing of exons are mediated by the spliceosome and occur simultaneously during the transcription of most human genes. Briefly, two main models for exonic circRNAs biogenesis coexist. Following the first model, known as "lariat-driven circularization" or "exon skipping", the skipping of one or more exons of the pre-mRNA results in a spliced mRNA and a lariat containing the exon(s) that will be circularized. In the second "Intron-Pairing-Driven Circularization" [17] model the presence of intronic repeated inverted complementary sequences, such as Alu sequences, enables the pairing of two introns that border exons that will be circularized. These sequences promote the spatial connection of the 5' splice site junction to that of the upstream 3' site [40] and leads to the formation of a circular transcript. Although both models may explain our results, we counted 220 circular junctions and 138 alternative linear splicing (i.e. exon skipping) junctions, not related to circRNAs for the most part. As an example, BUB1 has 28 circular junctions but 7 alternative linear splicing junctions. Consequently, and based on these observations, the “Intron-Pairing-Driven Circularization” model is more likely used.

Albeit counterintuitive, the number of circRNAs doesn’t increase with the number of exons/introns (Table 2). This might be explained by independent regulation and/or the fact that circRNA biogenesis also depends on several parameters such as the size of flanking introns and the presence of Alu sequences (see the “Intron-Pairing-Driven Circularization” model), the transcription speed [41] in particular through interaction with DNA forming circRNA:DNA hybrids (circR loops) [42], and the strength of splice sites [26]. POLD1_circRNA_3-2, BUB1_circRNA_9-2 and BUB1_circRNA_4-3 were more highly expressed than their parental mRNA, but the vast majority of circRNAs (217/220, 94%) represented less than 10% of the corresponding mRNA (Supplementary Table S1).

These results confirm that backsplicing is less efficient than splicing [43] and suggest a transcriptional backsplicing background as already known for linear splicing. More surprisingly, at the gene level, the overall expression of circRNAs compared with the corresponding mRNAs is often (9/20 genes) greater than 10% (Table 2). The values observed for POLD1, BUB1 and SMAD4, 5 times higher for the patients’ outlier group for POLD1 (Figure 2), can hardly be considered as background noise and could have a biological impact (Table 2). This suggests that each gene has its own circRNA repertoire with different, finely but uncoupled circRNA/mRNA expression levels [44,45]. Although a competitive mechanism has been ruled out, we demonstrate that the couple may be unbalanced in certain pathological situations [20,31].

Patients suspected of hereditary predisposition to colorectal cancer had a 1.93-fold higher circRNA/mRNA ratio as compared to controls, an imbalance especially driven by POLD1 with a 3.11-fold change. No correlation with the phenotype could be made for the subgroup of outlier patients with the highest circRNA/mRNA ratio. Our results suggest an increased production or a lack of clearance of POLD1 circRNA rather than a low POLD1 expression, a finding in line with known POLD1 contribution to hereditary CRC predisposition that is not related to decreased expression but to missense pathogenic variants in the exonuclease domain [46]. Interestingly, circPOLD1 itself has recently been shown to play a functional role in tumorigenesis. In a recent study, circPOLD1 expression was found to increase with lesion severity in cervical cancer and to promote oncogenic signaling through interaction with RNA-binding proteins such as YBX1, leading to activation of pathways involved in cell proliferation and tumor progression [47]. These data support the biological relevance of circPOLD1 and reinforce the interest of the imbalance observed in our study.

In any case, the circRNA/mRNA imbalance found in CRC patients suggests that the circRNAs studied play a role in some CRC, either as a cause or a consequence of other hidden defect(s). As circRNAs are expected to be finely regulated with cell-type dependent levels [48], this imbalance could also reflect cellular composition changes following tumorigenesis or chemotherapy, and may serve as a biomarker linked to the onset of cancer and progression [20]. Given the phenotype of the patients under study i.e. a cohort of patients suspected of hereditary predisposition to colon cancer, the small subgroup of 31 outlier patients is of special interest and requires further investigations e.g. cosegregation analyses of this imbalance in families, circRNA-mRNA tumor analyses and the search for newly formed peptide(s) to shed light on novel mechanism(s).

5. Conclusions

Overall, splicing and backsplicing analyses with SEALigHTS highlight the RNA architecture of each gene, indicating no evidence of a competitive mechanism or coordinated regulation between circular and linear transcripts. The observed imbalance is largely driven by a single circRNA with oncogenic potential (circPOLD1), suggesting that specific circRNAs may play a pivotal role in disease susceptibility. These findings open new avenues for the study of circRNAs in blood as potential biomarkers or functional contributors to cancer predisposition.