Cell-type specific features of circular RNA expression

Abstract

Thousands of loci in the human and mouse genomes give rise to circular RNA transcripts; at many of these loci, the predominant RNA isoform is a circle. Using an improved computational approach for circular RNA identification, we found widespread circular RNA expression in Drosophila melanogaster and estimate that in humans, circular RNA may account for 1% as many molecules as poly(A) RNA. Analysis of data from the ENCODE consortium revealed that the repertoire of genes expressing circular RNA, the ratio of circular to linear transcripts for each gene, and even the pattern of splice isoforms of circular RNAs from each gene were cell-type specific. These results suggest that biogenesis of circular RNA is an integral, conserved, and regulated feature of the gene expression program.

Figure 1. Bioinformatic and statistical method for detecting circular isoforms.

A) We created a custom database of all UCSC known-gene annotated scrambled exon-exon junctions. By mapping paired end 76 nt sequencing reads from poly(A) depleted RNA, we detected thousands of distinct circular RNA isoforms, including many cases where multiple circular isoforms are transcribed from the same locus. Our informatic pipeline required that one read (read 1) map to a diagnostic exon x - exon y junction (y< = x) and the other read map within the inferred circular isoform. B) Statistical scores improve filtering: We modeled the distribution of alignment statistics for reads from under an empirical null. Estimating the empirical null distribution of alignment quality for read 1 (required to map to a diagnostic circular exon-exon junction) and read 2 (which need not be junctional) allows us to compute a per-circular isoform FDR and statistically identify artifacts. The bulk of detected circles (illustrated at right) have alignment profiles that distinguish them from those detected under the null model. We used an estimated FDR threshold of .025, shown on plot.

Figure 2. Predicted circular isoforms are resistant to RNase R.

HeLa RNA was treated with RNase R or a mock treatment, and then subjected to qPCR with isoform-specific primers. The fraction of linear and circular isoforms was normalized to the value measured in the mock treatment. All tested circular isoforms resisted RNase R, including CYP24A1 (1106 nt), FAT1 (3283 nt), HIPK3 (1099 nt), RNF220 (742 nt), PVT1 (410 nt) and ABTB1 (130 nt). The depletion of the FAT1 circle (the largest circle tested) by RNase R may be due to occasional nicking by contaminating endonuclease activity. We hypothesize that the apparent increase in abundance of some upon RNase R treatment is due to more efficient priming in the RT after linear and ribosomal RNA depletion.

Figure 3. Northern blot shows the dominant isoform of CAMSAP1 is circular.

Northern blot on total HeLa RNA probed for exons 2 and 3 reveals three distinct bands: the largest, the canonical linear isoform of CAMSAP1 (∼7800 nt), a 1446 nt representing a circular isoform of CAMSAP1 containing exon 2, exon 3 and the intervening intron; a 425 nt band representing the fully spliced circular isoform of CAMSAP1 consisting of exons 2 and 3.

Figure 4. Intron length is enriched around exons defining circular RNA, but alone not explanatory of circular RNA expression.

Intron lengths flanking circular isoforms are calculated as described in the main text. A) and B) show the genome-wide distributions of flanking intron length, normalized by their quantile rank within a gene (shortest = 0; longest = 100); A) weights each isoform by total reads summed over all replicates and samples; B) counts each isoform once, provided it has at least 20 distinct read counts.

Figure 5. qPCR validation of relative circular RNA expression across cell type.

Total RNA from A549, AG04450 and HeLa cells was probed by qPCR using primers specific for circular isoforms of the indicated genes, and abundances were normalized using primer efficiencies estimated with a dilution series. Sequencing-based estimates are shown by comparison; sequencing values are depicted as a log fraction of total circle counts per experiment. Each qPCR and sequence value is calculated from the average of two biological replicates. Expression of LINC00340 and LPAR1 in HeLa is not detectable with the sequencing depth in this data, and these values were pinned at −60 on the log scale.

Figure 6. Quantitative regulation of circular to linear isoform ratios.

A) Examples of circular RNAs with cell-type dependent expression as predicted by a genome-wide statistical model. Circular isoform abundance was estimated as a fraction of total circular RNA expression per replicate, and error bars represent statistical variation (3.5 sd of the mean); replicates are depicted separately. Two circular isoforms of AMBRA1 are shown. For each gene, cell types are ordered left to right by monotonic increasing expression of the linear isoform as measured by RPKM, with RPKM value overlaid as a solid dot. Bar plot colors are consistent across cell types: red representing cancer cell lines, blue H1-HESC and greens are non-cancers; shading from dark to light representing endoderm, mesoderm and ectoderm. B) Across cell lines, no genome-wide trend between circle expression and linear transcript expression as measured by log RPKM.

Figure 7. Circular isoform expression patterns involves a variety of splicing patterns including proximal pairing and combinatorial expression.

Gene structures are represented along the axes with tick marks at splice site boundaries. Each circle is centered at the genomic coordinates corresponding to the donor and acceptor splice sites of the detected circular isoform. The length of the line segment is proportional to the log of the expression level of the circular isoform; the ring represents the maximum expression of the circular isoform across cell types. CYP24A1 and MCU exhibit striking expression preference for a single circular isoform. MBOAT2 exhibits a strong preference for a single splice acceptor site paired with multiple donors. ABCC1 exhibits preference for proximal pairs of splice donors and acceptors; CAMSAP1 exhibits a strong preference for either a particular single acceptor and or a particular single donor, whereas PICALM is an example of a gene with high combinatorial use of splice sites.

Figure 8. Circular isoform expression is regulated within individual genes.

Circular isoform expression across cell types is shown for LPAR1, RNF19B, and ZFAND6 as in Fig. 6. The barplots depict expression of the corresponding linear isoforms in RPKM units. Cell types are colored as red or green to highlight distinct patterns of circular isoform expression. Differential circular RNA isoform expression in LPAR1, RNF19B and ZFAND6 cannot be explained by differences in expression level of polyadenylated transcripts of these genes or sampling depth.