Global accumulation of circRNAs during aging in Caenorhabditis elegans

Background Circular RNAs (CircRNAs) are a newly appreciated class of RNAs that lack free 5′ and 3′ ends, are expressed by the thousands in diverse forms of life, and are mostly of enigmatic function. Ostensibly due to their resistance to exonucleases, circRNAs are known to be exceptionally stable. Previous work in Drosophila and mice have shown that circRNAs increase during aging in neural tissues. Results Here, we examined the global profile of circRNAs in C. elegans during aging by performing ribo-depleted total RNA-seq from the fourth larval stage (L4) through 10-day old adults. Using stringent bioinformatic criteria and experimental validation, we annotated a high-confidence set of 1166 circRNAs, including 575 newly discovered circRNAs. These circRNAs were derived from 797 genes with diverse functions, including genes involved in the determination of lifespan. A massive accumulation of circRNAs during aging was uncovered. Many hundreds of circRNAs were significantly increased among the aging time-points and increases of select circRNAs by over 40-fold during aging were quantified by RT-qPCR. The expression of 459 circRNAs was determined to be distinct from the expression of linear RNAs from the same host genes, demonstrating host gene independence of circRNA age-accumulation. Conclusions We attribute the global scale of circRNA age-accumulation to the high composition of post-mitotic cells in adult C. elegans, coupled with the high resistance of circRNAs to decay. These findings suggest that the exceptional stability of circRNAs might explain age-accumulation trends observed from neural tissues of other organisms, which also have a high composition of post-mitotic cells. Given the suitability of C. elegans for aging research, it is now poised as an excellent model system to determine whether there are functional consequences of circRNA accumulation during aging. Electronic supplementary material The online version of this article (10.1186/s12864-017-4386-y) contains supplementary material, which is available to authorized users.


Background
Circular RNAs (circRNAs) have recently been identified as a natural occurring family of widespread and diverse endogenous RNAs [1,2]. They are highly stable molecules mostly generated by backsplicing events from proteincoding genes. The expression trends of circRNAs are only recently emerging thanks to RNA-seq library preparation methods that deplete ribosomal RNA (ribo-depletion) rather than enrich for polyadenylated RNA. Most cir-cRNAs are derived from protein-coding genes, and thus one challenge in mapping and quantifying circRNAs is to distinguish reads that can be uniquely ascribed to circular molecules versus linear RNAs emanating from the same gene. Elements located within introns flanking circularizing exons play a role in promoting circRNA biogenesis [3][4][5][6], and several RNA binding proteins and splicing factors have been shown to influence circRNA expression [4,[7][8][9][10].
Despite the current interest in circRNAs, their functions are only beginning to emerge [2]. Recent reports have identified roles for circRNAs in regulating transcription, protein binding, and sequestration of microRNAs [11][12][13][14]. Some circRNAs can be translated via cap-independent mechanisms to generate proteins [15][16][17]. Moreover, cir-cRNAs have been implicated in antiviral immunity [18,19], and expression patterns of circRNAs in the brain suggest that they might serve important functions in the nervous system [20].
Several RNA-seq studies have found that circRNAs are differentially expressed during aging. Over 250 circRNAs increased in expression within Drosophila head tissue between 1 and 20 days of age [21]. Trends for increased circRNA expression have also been identified during embryonic/postnatal mouse development [10,22,23], suggesting that circRNA accumulation might begin early in development. We recently reported that circRNAs were biased for age-accumulation in the mouse brain [24]. In hippocampus and cortex,~5% of expressed circRNAs were found to increase from 1 month to 22 months of age, whereas~1% decreased [24]. This accumulation trend was independent of linear RNA changes from cognate genes and thus was not attributed to transcriptional regulation. CircRNA accumulation during aging might be a result of the enhanced stability of circRNAs compared to linear RNAs [13,25]. Age-related deregulation of alternative splicing [26,27] leading to increased circRNA biogenesis might also play a role.
C. elegans is a powerful model organism for studying aging. Previously, thousands of circRNAs were annotated from RNA-seq data obtained from C. elegans sperm, oocytes, embryos, and unsynchronized young adults [4,13]. Here, we annotated circRNAs from very deep total RNA-seq data obtained from C. elegans at different aging time points and uncovered 575 novel circRNAs. A massive trend for increased circRNA levels with age was identified. This age-accumulation was independent of linear RNA changes from shared host genes. Our findings suggest that circRNA resistance to degradation in post-mitotic cells is largely responsible for the age-upregulation trends identified both here in C. elegans, and possibly in neural tissues of other animals.

Genomic features of circRNAs in C. elegans
We set out to map C. elegans circRNAs genome-wide and quantify their expression at different ages using RNA-seq. We performed an aging paradigm of wild-type Bristol N2 worms that involved treatment with 5-fluoro-2′-deoxyuridine (FUdR) to inhibit DNA synthesis during embryogenesis and prevent egg-hatching. RNA from whole worms from three independent biological replicates corresponding to four aging time-points were collected: L4-larval stage (L4), Day-1 (D-1), Day-7 (D-7), and Day-10 (D-10) adults (Additional file 1: Figure S1). Ribo-depleted total RNA-seq library preparation was performed, followed by sequencing using paired-end 125 nt reads. De novo mapping of circRNAs from these total RNA-seq datasets was performed using the find_circ algorithm [13] (see Methods), with the added restriction of only annotating circRNAs that shared known exonic splice sites (ce11 UCSC genome). This strategy for mapping circRNAs requires the use of only "back-spliced reads" (Fig. 1a), which represent a very low percentage of a typical RNA-seq run (Additional file 2: Table S1). This approach was required to distinguish reads corresponding to circRNAs versus their linear counterparts that share the same exons. Of thẽ 1.9 billion paired-end reads generated, only 111,895 reads (0.006% of total) mapped to circRNA junctions after removal of PCR duplicate reads (Additional file 2: Table  S1). To annotate circRNAs with high confidence, we used a cut-off of 12 unique aligned reads per circRNA across the 12 libraries. This minimum read cut-off was more stringent compared to previous circRNA annotations in C. elegans [4,13]. Using our annotation pipeline (Additional file 3: Figure S2), we confidently identified a total of 1166 circRNAs. In this high confidence list, 591 circRNAs were previously annotated [4,13], and 575 were novel (Fig. 1b).
Most of the 1166 circRNAs mapped to codingsequence (CDS) regions of exons (82.1%), followed by circRNAs mapping to exons encompassing 5′ UTR regions and CDS (14.8%) (Fig. 1b). We found that 797 genes express at least one circRNA. As shown for haf-4 ( Fig. 1c), most genes that produced circRNAs expressed a single circRNA (576/797 genes). On the other hand, some genes were found to generate a large number of circRNAs. For instance, the afd-1 gene was found to generate 8 different circRNAs (Fig. 1d). Overall, 221 out of the 797 genes generated two or more circRNAs (Fig. 1e). The number of exons within circRNAs ranged from 1 to 13, but it was most common for them to harbor 2 exons (Fig. 1f ). Only 6.6% of the 1166 circRNAs contained 5 or more exons. The reliance of this analysis on back-spliced reads precludes the determination of whether these multi-exon circRNAs have introns or particular exons removed. As previously found for Drosophila and mice [21,24] there was a bias for circRNAs to emanate from the 5′ end of genes (Fig. 1g).
Base pairing between introns that flank circularizing exons is thought to bring 5′ and 3′ splice sites in close proximity to promote circRNA biogenesis over linearsplicing (Fig. 1h). We used Basic Local Alignment Search Tool (BLAST) alignment of introns that flank circRNAforming exons to identify reverse complementary matches (RCMs). We found that RCMs flanking circRNA loci were strongly enriched (P < 0.0001, Kruskal-Wallis test with Dunn's post-hoc test for multiple comparisons) compared to analogous introns flanking non-circularizing exons (Fig. 1h). Thus, consistent with previous reports [4], our analysis shows that C. elegans circRNAs tend to be flanked by introns that pair with one another.

Experimental validation of circRNAs
We next performed experimental validation of individual circRNAs annotated from our pipeline. One validation method was to prepare complementary DNA (cDNA) using random hexamers from total RNA, and then perform PCR using outward facing primers that should only amplify a back-spliced circRNA (Fig. 2a). The presence of a back-spliced junction was confirmed for 10/10 cir-cRNAs tested by Sanger sequencing of RT-PCR products (Additional file 4: Table S2). In addition, we confirmed a subset of circRNAs by treating total RNA with the exoribonuclease RNase R, which is known to preferentially degrade linear RNAs over circRNAs [5,13]. RT-qPCR experiments show that linear RNA cdc-42 was susceptible to degradation by RNase R, whereas 4/4 circRNAs tested were enriched upon RNase R treatment (Fig. 2b).
Although most circRNAs were of lower abundance compared to linear RNAs from the same gene, some cir-cRNAs annotated were of relatively high abundance. We  set out to confirm two of these high abundance cir-cRNAs using Northern blot analysis. In the case of the crh-1 gene, an increased abundance of reads aligning to an exon harboring a circRNA was clearly evident from linear RNA alignment tracks on Integrated Genomics Viewer [28] (Fig. 2c). We performed Northern analysis using a probe targeting this exon. This probe should detect both linear and circular transcripts of the crh-1 gene. As expected, an abundant circRNA migrating at the predicted size was detected (Fig. 2d). The expression of higher molecular weight linear RNAs was found to be diminished by RNase R treatment, whereas the circRNA bands were unaffected by RNase R treatment (Fig. 2d). We prepared polyA+ RNA from a column-based preparation, and collected and precipitated the unselected RNA (polyA-depleted). We found that polyA+ RNA had depleted circRNA levels relative to linear RNA. In contrast, polyA-samples showed enhanced levels of cir-cRNA relative to linear RNA (Fig. 2d). Analogous results were obtained using a probe for an abundant circRNA from the afd-1 gene (Fig. 2e). Together, these validations provide experimental support that our annotation pipeline detected bona fide circRNAs in C. elegans.

Global circRNA levels dramatically increase during aging
We next quantified the abundance of circRNAs from the different aging time-points. CircRNA read counts were normalized to their corresponding library size to obtain Transcripts Per Million reads (TPM) (Additional file 5: Table S3). Principal Component Analysis on the circRNA TPM values was performed on the 12 RNA-seq libraries (Fig. 3a). A close clustering of L4 to D-1, and of D-7 to D-10 was observed, suggesting that global circRNA expression levels reflect the age of the animal.
To further investigate trends in circRNA levels during aging, we plotted circRNA log 2 fold changes in TPM for all pairwise comparisons of the aging time-points (Fig. 3b). We found that 1052 circRNAs (90.2%) were at least 1.5-fold greater in D-10 versus L4 time-points, whereas only 37 circRNAs were 1.5-fold greater in L4. Similar trends were found in other pairwise comparisons between older (D-7, D-10) versus younger (L4, D-1) time-points (Fig. 3b). For instance, when comparing D-7 versus D-1, 80.8% of circRNAs were 1.5-fold or higher in the older time-point.
To gain statistical support for these dramatic aging trends we performed several additional analyses. The global expression of circRNA TPM values was compared across ages by non-parametrical Kruskal-Wallis test with Nemenyi post-hoc test for multiple comparisons (Fig. 3c).  Table S2). b RT-qPCR experiments on RNase R treated mixed age adult worms. Equal amounts of mock-treated and RNase R treated RNA were used for cDNA preparation prior to qPCR. Note the enrichment of circRNAs with RNase R treatment, whereas linear cdc-42 mRNA is not enriched. c RNA-seq track visualized using Integrated Genomics Viewer from D-7 worms showing read pileup at the crh-1 gene. Note the increased read number overlapping the circularized exon. cel_circ_0000438 and cel_circ_0000439 differ by 6 nucleotides in length at the 5′ end of the exon. d Northern blot using a probe overlapping the circularized exon of crh-1 (see panel c) detects bands corresponding to circRNA and mRNA from mixed age adult RNA. Relative circRNA to mRNA abundance is enriched in RNase R treated and polyA-samples compared to polyA+ samples. e Northern blot performed using a probe that detects afd-1 circRNA and mRNA. Red arrows denote circRNA bands. Black arrows denote likely linear RNA bands 2.5e-9 and 5.6e-4, respectively), which might reflect the ages being closer together. In order to identify the individual circRNAs with statistically significant changes in expression during C. elegans aging, we performed t-tests on TPM values between each time-point (P < 0.05, > 1.5 Fold Change (FC), False Discovery Rate (FDR) < 0.2). An overwhelming bias for upregulation of circRNAs during aging was uncovered. For instance, in the comparison of D-7 versus D-1, a total of 196 circRNAs were upregulated whereas only 2 were downregulated (Fig. 3d). Comparing D-10 versus L4 age time-points, 342 circRNAs were upregulated, whereas only 1 was downregulated (Fig. 3d). We next sought to identify circRNAs exclusively detected in one aging time-point. Using a 6 unique back-spliced read minimum cutoff, we plotted age-specific circRNAs (Fig. 3e). We detected 28 and 51 agespecific circRNAs in D-7 and D-10 libraries, respectively. In contrast, no age-specific circRNAs were found for the L4 and D-1 time-points. Note that the D-7 and D-10 specific circRNAs might be expressed in the younger worms, but at levels below the detection limits of the RNA-seq analysis. Together, these data demonstrate that circRNAs show an overwhelming bias for ageaccumulation in C. elegans.

Experimental validation of circRNA age-accumulation trends
We next set out to confirm circRNA expression trends for particular circRNAs that are generated from genes with interesting functions. We performed RT-qPCR validation for 11 individual circRNAs, including circRNAs generated from genes that are involved in lifespan determination (akt-1, crh-1, daf-16, daf-2). We selected cir-cRNAs with a variable range of overall expression levels for validation (Fig. 4a). Of these 11 circRNAs, 7 met the statistical threshold for increased expression between at least one old versus young time-point. On the other hand, two of the circRNAs were predicted to be significantly decreased during aging. Three of these circRNAs did not meet statistical significance for differential expression (D-10 versus L4) from the RNA-seq data, including two circRNAs of low abundance (Fig. 4a). Quantification of these same circRNAs by RT-qPCR revealed that 7/7 circRNAs predicted to increase during aging were significantly upregulated in D-10 versus L4 (Fig. 4b). Interestingly, the three circRNAs predicted to not increase during aging from the RNA-seq data (including two low abundance circRNAs from the daf-2 and daf-16 genes) were found to robustly increase during aging by RT-qPCR. Finally, the two circRNAs predicted to decrease during aging (from ddx-19 and nhr-65 genes) were found to be unchanged when tested by RT-qPCR analysis (Fig. 4b).
RT-qPCR analysis also showed that the tested circRNAs did not continue to increase in D-10 animals, consistent with the RNA-seq results (Fig. 4b). Notably, RT-qPCR validation for most circRNAs showed fold-changes greater than those detected by RNA-seq differential expression analysis. Most of these RT-qPCR quantified changes were >10-fold between L4 and D-7. Remarkably, changes in circRNAs from the gld-2 and daf-16 genes were > 40-fold increased between L4 and D-7 (Fig. 4b). Overall, these RT-qPCR validations strongly support the trend of greater circRNA abundance in old (D-7, D-10) versus young (L4, D-1) animals. These confirmations also suggest that the actual number of circRNAs that increase during aging might be greater than what was found to significantly change from the RNA-seq analysis (Fig. 3d).
Host gene-independent circRNA accumulation during aging Next, we performed differential expression analysis on linear RNAs among the different age time-points. Linear RNAs previously found to be differentially regulated during aging displayed similar expression trends in our datasets. For example, between D-7 and D-1, hsp-70 and cht-1 were upregulated during aging, whereas fat-7, ifp-1, and ifd-1 were downregulated (Additional file 6: Table S4). In contrast to circRNA trends, a global bias for linear RNA differential expression was not evident. Similar numbers of upregulated and downregulated linear RNAs between aging time points were identified using Cuffdiff [29] (Additional file 6: Table S4). Scatterplots comparing old versus young time-points for linear RNA levels ( Fig. 5a-d), and circRNA levels ( Fig. 5e-h) exemplify the stark contrast in the age-related trends.
Although linear RNAs lacked a global bias for increased levels during aging, it was still possible that increased transcription of circRNA-hosting genes could contribute to the circRNA expression trends. Thus, we analyzed whether cir-cRNA accumulation was independent of host gene expression. Density plots were generated to contrast circRNA total read count fold-changes versus their counterpart linear RNA changes from the same host gene (Fig. 6a-f). For this analysis, we used an expanded list of circRNAs, requiring at least 6 reads per age time-point. Using the CircTest algorithm [30], we tested for statistically significant changes in circRNA expression independent of host gene expression. From 1239 circRNAs tested, 459 showed a host gene independent significant expression (ANOVA P value corrected <0.05, Additional file 7: Table S6). In the old versus young time-point comparisons, a clear upward vertical shift was evident in the density plots (reflecting increased circRNA expression), and only a minor horizontal shift to the right (reflecting increased linear RNA expression) (Fig. 6b-e). This suggests that circRNA accumulation trends are largely independent of linear RNA changes. For comparisons between closer time points (D-1 versus L4 and D-10 versus D-7) the density plots lacked clear vertical or horizontal shifts (Fig. 6a, f). We next examined the host gene independent circRNAs for significant changes between time-points. We found strong trends for these circRNAs to be increased in older versus younger aging time-points (Fig. 6g). For example, 194 circRNAs were significantly increased in a host-independent manner between D-10 and L4, whereas no circRNAs were found to Relative Quantity Fig. 4 Validation of circRNA age-accumulation. a RNA-seq quantification of select circRNAs during aging (TPM fold-change with L4 set at 1). Total number of reads across all libraries for each circRNA is noted above graph. Labels display circRNA names with host gene in brackets. b RT-qPCR data for the same selected circRNAs as in A). Data is normalized to cdc-42 mRNA. Note the greater magnitude of age-related expression changes reported by RT-qPCR versus RNA-seq for all circRNAs. Note that for cel_circ_0001331, there was a significant reduction (P < 0.05) between D-1 and D-10 (Additional file 5: Table S3). Error bars represent SEM. *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001; ****, P ≤ 0.0001, two-tailed t-test compared to L4 decrease. Together, these results show that circRNAs globally accumulate during aging in C. elegans independently of host gene expression.

Discussion
This study is the first to report age-accumulation of cir-cRNAs in C. elegans. Previous studies have documented the bias for circRNAs to be increased during aging in neural tissues of Drosophila and mice [21,24]. Interestingly, the trends uncovered here during C. elegans aging are much more dramatic. Many confirmed expression trends were >10-fold increased between L4 and D-7. Of the hundreds of differentially expressed circRNAs, the vast majority increased with age (Figs. 3, 4, 5). Experimental confirmations of the age-accumulation trends by RT-qPCR were of greater magnitude than those reported from the RNA-seq analysis. Although we have higher confidence in the accuracy of circRNA quantification by RT-qPCR over RNA-seq due to the low depth of junction-spanning sequencing reads, it has been discussed that RT-qPCR might overestimate cir-cRNA levels due to rolling circle reverse transcription [1]. It is notable that RT-qPCR experiments revealed two circRNAs to accumulate with age that were not significantly increased in the RNA-seq data, and two circRNAs predicted by RNA-seq to decrease with age were found to not change during aging (Fig. 4). Together, the evidence presented here suggests that many circRNAs not passing statistical thresholds for increased expression during the C. elegans aging paradigm are of too low abundance to be accurately quantified by the limited number of back-spliced RNA-seq reads. Thus, our reported numbers of significant age-accumulated circRNAs are most likely gross underestimations.
Why is the age-accumulation trend of circRNAs in C. elegans much stronger than in other organisms tested so far? After completing development and a brief reproductive period, C. elegans spends the remainder of its adult life comprised almost exclusively of post-mitotic cells [31]. The FUdR treatment employed in this study inhibits DNA synthesis and is commonly used to prevent egg-hatching in aging and life-span studies of C. elegans [32], and also reduces the presence of proliferating cells. We have previously proposed that either the high stability of circRNAs or alterations in alternative splicing with age could both contribute to the increase of circRNAs during aging [2]. The trends observed here in C. elegans argue for a low rate of circRNA decay being the responsible mechanism. We propose that the dramatic genome-wide increase of circRNA levels during aging are a consequence of the dominance of postmitotic cells in adult C. elegans combined with the high stability of circRNAs. As neurons are post-mitotic, perhaps this can also explain why age-accumulation is most notable in Drosophila heads (which are rich in neurons) [21], and in brain regions of mice [24]. To investigate the potential functional significance of age-accumulated circRNAs, we performed Gene Ontology (GO) analysis on the host genes of expressed C. elegans circRNAs, and found many significantly enriched categories, including an enrichment in the Biological Process category of "determination of adult lifespan" (Additional file 8: Figure S3, Additional file 9: Table S5). A clear bias for particular GO categories, however, not uncovered by these efforts. It is certainly possible that trans functions of circRNAs are completely distinct from the curated roles of their host genes.
Future work can now take advantage of the powerful genetics of C. elegans to delineate aging functions of individual circRNAs. Generating loss-of-circRNA mutants in C. elegans by disrupting base-pairing of flanking introns could be a fruitful approach. Various RNAs found to be differentially regulated during aging were subsequently found to impact lifespan in C. elegans mutant analysis, including linear RNAs [33], microRNAs [34] and long non-coding RNAs [35]. Given that these mutant studies on lifespan were based on comparatively modest fold-changes during aging, the massive upregulation trends for circRNAs provide solid rationale for disrupting or overexpressing circRNAs in C. elegans and testing for effects on lifespan and healthspan. However, one should also consider that the aging process might be impacted generally by the total compendium of hundreds of circRNAs accumulating in cells, as opposed to individual circRNAs. Thus, non-conventional approaches to alter the expression of many circRNAs simultaneously might be required to uncover age-related functions of circRNAs.

Conclusion
We have shown that circRNAs accumulate during aging on a genome-wide scale in C. elegans. Given that cir-cRNAs are also increased during aging in Drosophila [21] and mice [24], age-accumulation of circRNAs in postmitotic cells appears to be a universal phenomenon.

C. elegans maintenance and culturing
The C. elegans Bristol N2 wild type strain was grown and maintained as previously described [36]. To synchronize populations, gravid adults were bleached, eggs were collected and left overnight in 1X M9 buffer with rocking. Starvation arrested L1 larvae were placed on 150x100mm NGM plates with 10X concentrated E. coli OP50 as a primary food source, and kept at 15°C. At the L4 larval stage, animals were collected using a 25 μm nylon mesh (Sefar) and either harvested for RNA extraction, or placed on 150x100mm NGM plates containing 75 mM FUdR (Sigma Aldrich) with 10X concentrated E. coli OP50 and kept at 15°C until animals were harvested at D-1, D-7 and D-10 of adulthood (see Additional file 1: Figure S1). Three biological replicates were collected for each timepoint and used for RNA extraction and RNA-seq library preparation.

RNA extraction
A 250 μl mixture of animals in 1X M9 buffer was added to 750 μl of TRIzol LS (ThermoFisher Scientific) and immediately frozen with liquid N 2 . Lysates were freeze/ thawed at −80°C, disrupted with Mixer Mill 400 (Retsch) and Dounce homogenizer (Corning) to break apart the cuticle of animals. Any cellular debris was removed by low-speed centrifugation. RNA was extracted using the Purelink RNA mini-kit with DNAse I treatment (Ambion). RNA quality was assessed by Bioanalyzer (Agilent) and quantified using Quant-iT RiboGreen RNA Assay kit (ThermoFisher Scientific).

Library preparation and high-throughput sequencing
Libraries were prepared using the Illumina TruSeq Stranded Total RNA Library Prep Kit as recommended by the manufacturer (Illumina) with modified conditions to increase the size of the cloned fragments (fragmentation at 85°C × 5 min). Barcoded libraries were sequenced at New York Genome Center (New York, NY) using the Illumina HiSeq 2500 system to obtain pairedend 125 nt reads. Raw FASTQ files from the RNA-seq data were deposited at the NCBI Sequence Read Archive (BioProject: PRJNA357503, individual accession numbers are listed in Additional file 2: Table S1). PolyA+ RNA and polyA-RNA were obtained from total RNA using the NucleoTrap mRNA kit (Machery-Nagel). RNA bound to oligo(dT) beads was carried through the complete polyA+ enrichment according to the manufacturer's protocol while RNA that remained in the column flow-through (unbound to the oligo(dT) beads) was precipitated with isopropanol and washed with 70% ethanol. Northern analysis was performed as previously described [24], with probe hybridization taking place overnight at 42°C, and all blot washing steps at 50°C.

CircRNA prediction and mapping
For de novo identification of circRNAs, a computational pipeline was carried out as previously described with filtering for duplicate reads and removing circRNA annotations spanning multiple genes [24] (Additional file 3: Figure S2). We obtained circRNA junction spanning FASTA sequence templates of 200 nucleotides using the C. elegans genome (from WBcel235/ce11 UCSC genome) as a reference. Assignment of circRNAs to their corresponding parental genes was performed using custom R scripts based on the library GenomicFeatures [37].

CircRNA normalization and quantification
A minimum of 12 reads summed from all 12 libraries was required for each circRNA in order to be considered for downstream analyses. To account for variability due to differences in library size, counts attributed to individual circRNAs were normalized to total read count to obtain circRNA TPM values. We calculated the corresponding fold changes by pair-wise comparisons of the average TPMs between time-points. Individual unpaired t-tests were performed for these pairwise comparisons. P values were corrected for multiple hypothesis testing with FDR < 0.2. To identify stage-specific circRNAs (Fig. 3e), and for pair-wise comparisons of host gene expression independent changes using CircTest [30] (Fig. 6), we required a minimum of 6 reads per time-point.

CircRNA host gene independence test
The CircTest pipeline was implemented to test for host gene independent circRNAs [30]. DCC 0.0.4 (https:// github.com/dieterich-lab/DCC) was used to quantify host gene read counts. In accordance with the DCC pipeline, FASTQ files were mapped with STAR 2.5 [38] using the recommended parameters. Picard 2.2.4 (http:// broadinstitute.github.io/picard/) was used to remove duplicated reads from the BAM files. The circRNA raw read counts (using an annotation list of circRNAs meeting a 6 read minimum requirement) together with the DCC linear counts were used as input for the Circ.test function which was run with default parameters. To define a circRNA as host gene independently expressed, we set a cutoff of P < 0.05, FDR < 0.05.

Quantification of linear expression
Linear RNA-seq reads were mapped to the C. elegans ce11 annotation using TopHat [40] with default settings. To quantify the differential expression across time-points, we used Cuffdiff [40]. Genes with fold changes > 1.5 and a Benjamini-Hochberg corrected P value < 0.05 (default parameter used in the Cuffdiff algorithm) were considered to be differentially expressed.

GO analysis
We used the Cytoscape plugin ClueGO [41] (together with the GO Annotation (GOA EMBL-EBI) (http:// www.ebi.ac.uk/GOA) released on 11/17/2016. The background set for GO was all genes with annotated GO entries for the GOA annotation released on 11/17/2016. Network specificity was set to "medium". The enrichment statistic used was a two-side hypergeometric test, correcting for multiple testing with the Bonferroni method. The cutoff for considering a term as enriched was set at P < 0.05. To reduce the number of redundant terms we used the GO term grouping option, which uses a Kappa score to collapse terms that share elements. We set the minimal number of elements in a group at 3. Bar graphs of significant GO terms were created using the -log (P-value).

RCM and motif analysis
To identify RCMs, pairs of intron sequences flanking the circRNAs were extracted from the C. elegans genome (WBcel235/ce11) using custom scripts available at: https://github.com/alexandruioanvoda/IntronPicker. The corresponding sequences were used as input for the RCM analysis using custom scripts (https://github.com/ alexandruioanvoda/autoBLAST) that employed BLAST (parameters: blastn, word size 7, output format 5) to identify matches. Exons 2 and 8 from non-circRNA generating exons were used as controls to account for the possibility that intron pairing might be influenced by exon location within genes.
OD010440). None of these funding bodies were involved in the design, data collection, analysis, interpretation of data, or in writing the manuscript.

Availability of data and materials
The dataset supporting the conclusions of this article is available in the NCBI Sequence Read Archive repository. This information can be accessed from http://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA357503. See Additional file 2: Table S1 for accession numbers. Ethics approval and consent to participate Not applicable Consent for publication Not applicable