Circular RNAs are long-lived and display only minimal early alterations in response to a growth factor

Abstract

Circular RNAs (circRNAs) are widespread circles of non-coding RNAs with largely unknown function. Because stimulation of mammary cells with the epidermal growth factor (EGF) leads to dynamic changes in the abundance of coding and non-coding RNA molecules, and culminates in the acquisition of a robust migratory phenotype, this cellular model might disclose functions of circRNAs. Here we show that circRNAs of EGF-stimulated mammary cells are stably expressed, while mRNAs and microRNAs change within minutes. In general, the circRNAs we detected are relatively long-lived and weakly expressed. Interestingly, they are almost ubiquitously co-expressed with the corresponding linear transcripts, and the respective, shared promoter regions are more active compared to genes producing linear isoforms with no detectable circRNAs. These findings imply that altered abundance of circRNAs, unlike changes in the levels of other RNAs, might not play critical roles in signaling cascades and downstream transcriptional networks that rapidly commit cells to specific outcomes.

Figure 1.

Unlike mRNAs and microRNAs, circular RNAs of mammary cells display minor changes in response to an extracellular cue. (A) Venn diagrams presenting genomic origins of circRNAs found in MCF10A mammary cells. The analysis comprises 1498 circRNA molecules we identified in MCF10A cells using RNA-sequencing. Of these, 1451 molecules (97%) overlap known transcripts. The remainder 47 circRNAs are either intergenic (14 transcripts) or antisense to known transcripts (N = 33). Inner circles represent Quantified fractions, meaning transcripts we followed also by using PCR; For circRNAs derived from known transcripts, the majority (>90%) of the Quantified fraction refers to circRNAs, the response of which to EGF was assayed using divergent and convergent sets of primers, while for the remainder of the Quantified fraction, including circRNAs derived from either antisense or intergenic regions, measurements were performed using only divergent sets of primers. (B) MCF10A human mammary epithelial cells were starved overnight for serum factors. Thereafter they were treated with EGF (10 ng/ml) for the indicated time intervals. High-throughput PCR and specific primers were applied on isolated RNA samples to amplify 241 of 1498 circRNA species previously identified using RNA sequencing. This group included >50% of the most abundant candidates, along with dozens of other circRNAs with varying expression levels. The presented heatmap (right panel) depicts time-dependent alterations in expression levels of specific circRNAs. These alterations are compared to mRNA and microRNA alterations we previously observed, using microarrays, while stimulating MCF10A under identical conditions (5,36). Note that all previously analyzed miRNAs are represented, but in order to match the size of the circRNA population, only randomly selected, 16.1% of all MCF10A's mRNA molecules, are depicted in the heatmap. Data were normalized to time zero and ordered according to the time point corresponding to the maximal change. Red squares represent an increase and blue squares represent a decrease, as shown in the scale bar on the right. CircRNA results represent biological duplicates performed in technical triplicates. (C) A histogram showing the range of abundance changes of mRNAs, miRNAs and circRNAs (N = 3608, N = 164 and N = 288, respectively) displayed by EGF-stimulated MCF10A cells. To construct the histogram, the maximal change value (induction or repression) along the stimulation interval (240 min) was found for each RNA molecule. Note that circRNAs exhibit narrower dynamic range (highlighted region) than mRNAs and miRNAs (P < 1e-100, F-test, Bonferroni corrected for multiple comparisons). Note that only 3 time points (30, 90 and 240 minutes) were available for 47 of the presented circRNAs.

Figure 2.

circRNAs are less dynamic and abundant relative to linear isoforms derived from the same host gene. (A) A linear diagram exemplifying primer pairs (arrows) for parallel measurements of levels of circular and linear RNA isoforms derived from the same host gene. The ERBB2 (chromosome 17) gene is shown as an example. The respective circRNA (red) spans five internal exons. The divergent primers (red arrows; circular transcript) flank the non-canonical splicing (circularization), whereas the convergent primers (blue arrows) were designed for measuring the corresponding linear isoform. (B) PCR and four pairs of primers (designed as in A) were used to follow the response of the indicated linear and circular isoforms of four EGF-regulated genes. Note that mRNAs corresponding to PUS7 and SLC25A32 are induced when MCF10A cells are stimulated with EGF, whereas PCMTD1 and PIK3C2G mRNAs undergo down-regulation under the same conditions. Shown are means ± S.E. of triplicates. P-values were calculated using two-way Anova with time and RNA type as categorical factors. (C) A histogram comparing the abundance of circRNAs and the corresponding linear transcripts (N = 203) of MCF10A mammary cells. Note that circRNAs are, on average, 36 times less abundant than the corresponding linear isoforms (P = 5.45e-78, t-test, two-tail, paired; A.U., arbitrary units). (D) A dot plot presenting the ratio between expression levels of each circRNA of MCF10A cells we characterized and the respective linear transcript, which is derived from the same host gene. The ratios are presented against the abundance of the latter isoform. Note that most circRNAs are expressed at lower levels than the corresponding linear transcripts; only five circRNAs (red dots; identified by names) exceeded the abundance of their respective linear isoforms, but the absolute expression levels of all five circRNAs are relatively low.

Figure 3.

Analysis of newly transcribed RNAs reveals that circRNAs are more stable and static than the linear isoforms derived from the same host genes. (A) MCF10A cells were treated with EGF as in Figure 1B and RNA was simultaneously metabolically labeled using 4-thiouridine (4sU), for the indicated time intervals. RNA was extracted (Total-RNA) with Trizol, biotinylated and purified on streptavidin magnetic beads (denoted 4sU-RNA). Flow-through RNA was also collected (denoted FT-RNA). Thereafter, RNA was reverse transcribed and quantified using high-throughput real time PCR (Fluidigm). The boxplot diagram shows enrichment of newly transcribed RNA (4sU labeled) in linear isoforms (blue dots) relative to the respective circRNA isoforms (red dots; P < 1e-20, t-test, two-tailed distribution, unequal variance, Bonferroni corrected for multiple comparisons, N = 61). Shown are the log₂ fold enrichments of 4sU-labeled RNA versus total RNA (Y axis). (B) MCF10A cells were metabolically labeled using 4sU, for 1 or 2 h. Thereafter, RNA was extracted, biotinylated and purified on streptavidin magnetic beads. Flow-through RNA was also collected. Next, RNA was reverse transcribed and quantified using high-throughput real time PCR (Fluidigm). Presented are the half-lives of 60 circRNAs and their corresponding linear counterparts. Half-life values were calculated from two samples, which were labeled with 4sU for 1 or 2 h and then averaged. All data were corrected for any bias introduced due to low uridine (short length) of RNA species (see ‘Materials and Methods’ section). The circRNAs (red dots) and their linear counterparts (blue dots) were sorted according to the difference between their half-lives from high to low. Error bars represent standard errors. The calculations of half life were performed using the HALO software (39) and the ratios between newly transcribed (RNA-4sU) and pre-existing RNAs (RNA-FT). (C) PCR and four pairs of primers were used to follow the response of newly transcribed RNA (solid lines) and total RNA (dashed lines) of both linear and circular isoforms of the indicated EGF-regulated genes. Shown are calculated fold changes of expression levels displayed by EGF-treated cells relative to basal (unstimulated) expression levels, measured in cells labeled with 4sU for the indicated intervals (Y axis). Shown are means ± S.E. of duplicates.

Figure 4.

Genes giving rise to circRNAs frequently express the respective mRNAs, and they are transcriptionally more active than genes transcribed into linear RNA only. (A) Basal (unstimulated) expression levels of mRNAs transcribed from genes that also give rise to circRNAs (red bars; N = 1051) are compared to mRNAs transcribed from genes that give rise only to linear RNAs (no detectable circRNAs; blue bars; N = 21 872). MCF10A RNA-seq data were obtained from a publicly available RNA-sequencing dataset (42). RPKM—reads per kilo basepairs per million reads. P = 6.13e-212, Mann–Whitney-U test. Note that only five circRNAs (leftmost quantiles) displayed no linear counterparts. (B) Basal H3K27 acetylation levels at promoter regions (transcription start site (TSS) ± 3000 bp) were determined in MCF10A cells. The acetylation levels displayed by genes transcribed into both circRNAs and linear RNAs (N = 781; red line) are compared to genes transcribed into linear RNAs with no indication of the respective circRNAs (N = 4992; blue line). The tag counts were normalized to 10 million sequencing tags for each sample. Shown are curves of tag densities of 25 nt bins. Tag density units are presented per basepair (bp) and per gene. Note that H3K27 acetylation levels are higher for genes producing circRNA compared to the control gene set (P = 1.26e-26, Mann–Whitney-U test). (C) Basal DNA methylation levels at the promoter region (TSS ± 1000 bp) of genes transcribed into both circRNAs and linear RNAs (N = 781) are compared to methylation levels displayed by genes that give rise to linear RNA only (N = 4992). DNA methylation data of MCF10A cells was derived from the ENCODE dataset. DNA methylation levels in the promoter region (−1000 to +1000 bp relative to the TSS) of each circRNA-producing host gene and control gene (no evidence for circRNA production) were tallied and a cumulative fraction was calculated and plotted for both groups. Note that DNA methylation levels are lower in genes producing also circRNAs compared to a group of genes giving rise to linear RNAs only (P = 3.15e-12, Mann–Whitney-U test). (D) Basal DNA methylation levels in the region between the TSS and TES (gene body) of genes transcribed into both circRNA and linear RNA (N = 399) are compared to levels of methylation displayed by genes that give rise to linear RNA with no detectable circRNAs (N = 2220). The methylation levels were tallied and normalized to gene length for each circRNA-producing host gene and for control genes (no evidence for circRNA production). A cumulative fraction was calculated and plotted for both groups. Note that DNA methylation levels in the body of genes producing circRNAs are lower compared to a control group of genes (P = 9.46e-57, Mann–Whitney-U test).

Figure 5.

Circular RNAs expressed by human mammary cells show no enrichment for binding sites of the microRNAs found in the same cells. (A) CircRNAs expressed in MCF10A cells were analyzed for the presence of binding sites (7 or 8 nt long) for miRNAs expressed by the same cells (33% uppermost percentile). The number of sites was tallied for each circRNA–miRNA pair, and the distribution of values is plotted. The solid and dotted curves indicate the averaged and 95% upper percentile, respectively, of results when repeating the analysis 1000 times using different permutations of site sequences (see ‘Materials and Methods’ section). (B) Shown are the numbers of hexamers within circRNAs of MCF10A cells that potentially pair to miRNAs expressed in the same cells (orange line), or to scrambled miRNA regions (black line). Error bars represent the distribution of values in 1000 control scrambled miRNA regions. Scrambled miRNA regions preserved mono-nucleotide and GC dinucleotide composition of the original 6 nt window. For each miRNA highly expressed in MCF10A cells (33% upper percentile), starting from the 5′ end, we checked 6-nt-long windows for complementarity to hexamers within the collection of circRNAs expressed in MCF10A cells (N = 1326). The total number of matches was normalized to the number of unique 6 nt windows analyzed. Highlighted is the central region of miRNAs found in a previous study to have statistically significant pairing to coding sequences (55).