Circular RNA repertoires are associated with evolutionarily young transposable elements

Essential revisions:

1. The authors collected three organs from five species for RNA sequencing to annotate circRNAs for their subsequent analysis, and obtained about 4,500 circRNAs in human and much fewer in other four species. However, significantly more circRNAs have been previously reported. for example, in some publicly-available circRNA databases, such as circAtlas (Wu et al., Genome Biol 2020, PMID: 32345360), circBase (Glažar et al., RNA 2014, PMID: 25234927), CIRCpedia (Zhang et al., Genome Res 2016, PMID: 27365365), more than 400,000 circRNAs are annotated from human samples. Very recently, two publications by using third generation sequencing technologies further confirmed a large number of circRNAs in human (including Xin et al., Nat Commun 2021, PMID: 33436621 and Zhang et al., Nat Biotechnol 2021, PMID: 33707777). It is not clear why the authors only used their own circRNA annotation for the current study, as this might lead to a significant bias in the downstream analysis (e.g. affecting the parental genes of circRNAs to be compared to genes that do not generate circRNAs).

We thank the Reviewers for raising the issue – this is indeed an important point to clarify. It is true that the publications and circRNA databases cited by the Reviewers have reported on very high numbers of molecularly distinct circRNA molecules, whereas we base our study on circRNA counts that are in the low thousands.

1) The most important answer to why we “only used [our] own circRNA annotation for the current study” lies in the unavailability of comparable data from previous studies. It is correct that we could have retrieved hundreds of thousands of circRNA sequences that have been reported from human samples and cell lines in the literature, yet we would not have been able to obtain matched data from, for example, opossum or macaque organs. In our experience, it is key for clean, unbiased evolutionary analyses to generate datasets in parallel and with the same methodology. We therefore do not think that it would have been a good option to use the large published human datasets and compare them with data from the other species that we would anyway have to generate ourselves. Thus, the most straightforward solution consisted in acquiring all biological samples ourselves, subsequently allowing us to ensure the best possible consistency of organs and of the data production and analysis pipelines. In particular, (i) all organs originated from young, sexually mature male individuals; (ii) our highquality biological samples, RNA and library preparations are fully traceable; (iii) the data was produced for all samples in parallel using the most appropriate experimental protocols.

We also note, however, that we actually did analyse complete published datasets, too, in order to independently validate some of our findings – please see the analysis of the human-mouse-rat datasets from Werfel et al., 2016; now in Figure 3 – Figure supplement 4.

2) Second, why do we “only” report on a few thousand circRNAs, as compared to the published hundreds of thousands? To understand this difference, we need to elaborate slightly more. The first important point to make is that the publications cited by the reviewers – for example the two most recent ones by Xin et al., and by Zhang et al., both from 2021 – actually used relatively permissive conditions for circRNA detection and annotation. We will explain this just below using the publication by Xin et al., in addition, at the end of this reply to Essential Revision 1 (see “Supplementary Rebuttal Information”), we provide a table that analyzes data from the three databases mentioned by the reviewers (circAtlas, circBase, CIRCpedia) and that summarizes why we think that an inclusion of the reported circRNAs would not have been helpful for our study.

With regard to the example of Xin et al., one should first note that this particular study is based on RNase R-treated samples, without the use of additional mock-treated libraries. We, however, used both library types and the enrichment factor between treated vs. mock treated as one of the criteria to define high-confidence circRNAs (see Essential revision 3). Moreover, the study by Xin et al., is based exclusively on a cancer-derived cell line, HEK293. HEK293 cells are hypotriploid (64 chromosomes) and are cytogenetically not very stable (as compared to organ samples that are the basis of our study). Chromosomal instability can easily be envisioned to lead to exon rearrangements that would lead to de novo creation of BSJs and abnormal splicing events that are not biologically relevant in vivo. Next, Xin et al., prepared 6 libraries - i.e. 3 technical library replicates sequenced from 2 biological HEK293 cell samples - and reported on around 40’000-60’000 BSJs/circRNA isoforms per library. Of note, before applying our enrichment/filtering steps to remove noise/sporadic BSJ detection in our samples, the numbers of BSJs that we detect are within the same order of magnitude (see Supplementary File 2). Please also note that in Xin et al., even for the “high confidence BSJs” (defined as BSJs with absolute read count = 2 across all libraries), the similarity even between technical replicates appears to be relatively modest (see Figure 2a-b in Xin et al.,), which indicates that there may be an issue with rather lenient thresholds for circRNA calling (“noise”).

In summary, we think that estimates such as the “more than 400’000 circRNAs from human samples” are the result of (i) the inclusion of many cancer cell line samples in the databases that contribute with abnormal gene expression and splicing events and of (ii) relaxed thresholds – in some studies a single reported BSJ read in a single sample is sufficient for entry into the database. We are not convinced that these (in part stochastically occurring?) BSJs are suitable for downstream analyses like those that we carry out in our study.

3) As mentioned in the above paragraphs and shown in Supplementary File 2, we actually do obtain relatively high numbers of BSJs/circRNAs per library (i.e. in the tens of thousands); why, then, do we base our downstream analyses nevertheless only on a few thousand circRNAs?

The reason lies in various filtering steps, which are described in Materials and Methods and also shown in Supplementary File 2 for the example of the human libraries. For human liver, cerebellum and testis we thus begin with approximately 25,000, 50,000 and 40,000 detected BSJs, respectively, per tissue sample (=per biological replicate). After removing poorly mapped events and merging across biological replicates and RNase-treated and -untreated samples, these correspond to a total number of 66,405 (liver), 137,615 (cerebellum) and 94,831 (testis) unique BSJs (RNase-treated and -untreated combined). Of the various filtering steps indicated in Supplementary File 2, the ones that decrease the numbers most strongly are the following: Filtering step 2 leads to a reduction by about 75%; briefly, in this step we remove those circRNAs that can be detected in RNase R-treated samples, but not at all in any of the untreated ones. The rationale behind this filtering step is to remove the very low expression BSJs that would be considered “absent” if only based on RNase R-untreated RNA samples. In Filtering step 3 we analyse enrichment of BSJs in RNase R-treated vs. untreated samples; we remove the BSJs for which the log2-enrichment of treated vs. untreated libraries is less than 1.5. This step removes around 15% of the original read count. Finally, in Filtering step 4, we calculate the mean RPM value for each BSJ across untreated replicate libraries, and keep the BSJs with an RPM > 0.05. We consider these loci strong circRNA candidates and use them for our subsequent analyses.

Taken together, we think we understand the Reviewers’ concern that “it is not clear why the authors only used their own circRNA annotation for the current study, as this might lead to a significant bias in the downstream analysis (e.g. affecting the parental genes of circRNAs to be compared to genes that do not generate circRNAs)”. However, we would like to argue that our approach is actually suited to precisely reduce bias that could occur when combining other studies’ circRNA annotations where available (i.e. using circRNAs that were discovered in other human or mouse samples/studies), and compare them to our data from the more unconventional species such as opossum or macaque, for which only little data is publicly available. Second, our filtering pipeline is designed to provide a pool of robust parental genes and circRNAs rather than a collection of inconsistently or sporadically expressed circRNAs (e.g. circRNAs that may be specific only to some cancer cell lines).

Supplementary File 2 has been newly added to the manuscript; in addition, a new section in Materials and methods, entitled “Generation of high confidence circRNA candidates from the comparison of RNaseR-treated vs. -untreated samples” now explains in detail the various filtering steps.

2. A complementary question is whether the model put forward by the authors applies to highly abundant circRNAs, which are likely to be enriched in functional circRNAs, as the sense in the field is that only highly abundant circRNAs -particularly in the nervous system- might be functional. Can the authors partition their assessment of conservation for the different quantiles of expression? It is indeed possible that only circRNAs in the top quantile of expression will provide insights into the potential functional importance of circRNAs. The authors could also take the opposite approach: for those few circRNAs for which some functionality has been proven, are the pattern of evolution and elements flanking the circularizable exons similar to the ones they describe generally?

We thank the Reviewers for this excellent suggestion. We have taken the proposed approach and compared properties of the top 10% most highly expressed circRNAs with the lower 90%. We have included these interesting analyses in the revised manuscript – mainly in Supplementary File 9, with further information in Supplementary File 10 and Figure 3—figure supplement 5.

First, we compared the general properties of the highly expressed vs. other circRNAs. In short, this analysis revealed that highly expressed circRNA are more likely to be expressed in all three tissues of a species, to be generated from a hotspot, and to be shared between species. Second, we went to the generalised linear model (Table 1) that we used in our study to distinguish parental vs. non-parental genes. Within the model, we compared the prediction values for the highly vs. lowly expressed circRNAs. For opossum, rhesus macaque and human, we found that the highly expressed parental genes were associated with higher prediction values. This effect was mainly driven by GC content (opossum, rhesus macaque, human) and genomic length (rhesus macaque, human).

In summary, our analyses therefore suggest that the overall model that we put forward in our study also applies to highly expressed circRNAs. The overall higher shared-ness of this group of circRNAs makes it probable that potentially functional circRNAs are more likely to be found in this group as well. These findings have been integrated in the manuscript.

Finally, we have analysed known functional circRNAs in our model – please see Figure 3-Figure supplement 3, “Properties of ‘functional circRNAs’ from literature.”

3. The authors identified circRNAs using RNAseR-treated samples. While RNAseR enriches generally for circRNAs, the authors don't have data reporting the abundance of the circRNAs in mock treated samples. Without a mock control is impossible for the authors to determine the actual resistance to the degradation by the exonuclease and hence determine which potential backsplicing junctions are really coming from circular molecules. The authors could do at least for some of the tissues/species show that the circRNAs they identified have been previously shown to be real (e.g., by a canonical comparison mock vs RNAseR treated samples).

We apologize for not having described more explicitly the experimental protocol for library generation in the main text, but only in the Materials and methods section and in the Supplemental Material. We fully agree with the Reviewers’ assessment that to conclude actual resistance to RNase R treatment (and thus, likely, circularity), a mock control is required. All RNase R-treated libraries of our study were accompanied by a parallel library from the same RNA preparation, but without prior RNase R treatment. In our previously submitted version, this information could be found, for example, in Figure 1—figure supplement 1. Moreover, we had described the use of the mock-treated samples in the materials as methods in two separate paragraphs. We therefore feel that this “Essential Revision” item does not require additional analyses (or experiments), but rather an improvement in presentation in the main text, in order to avoid misunderstandings. We have thus adapted the manuscript now giving all essential information at the beginning of the Results section. We have also moved the figure panel that shows the RNase R treatment vs. control strategy from Figure 1—figure supplement 1 to Figure 1 to show the protocol details more prominently. Finally, in the Materials and methods section, we have expanded the explanation of the various steps for filtering and RNase R enrichment analyses – see new section entitled “Generation of high confidence circRNA candidates from the comparison of RNase R-treated vs. -untreated samples”.

4. For the analysis of enriched TEs of flanking introns of circRNAs or circRNAs shared across species, except the ages of TEs, other factors should be also included in these analyses. For example, whether pairing abilities of RVCs alone can explain the phenomenon? Whether the number of RVCs in the flanking introns is related to the conserved or species specific circRNAs? In addition, In this case, the general conclusion about circRNA biogenesis and TE insertion events can be biased.

We thank the Reviewers for the opportunity to clarify these issues. To address them, we have conducted a more refined analysis of the dimer composition in the flanking introns of shared and species-specific circRNA loci. The obtained results have led to a change in the presentation and layout of the final part of the Results section (“Flanking introns of shared circRNA loci are enriched in evolutionarily young TEs”) and of Figure 5.

We would like to point out that our comparison of TEs in flanking introns in shared and species-specific circRNAs was – already in the first version of the manuscript – based on multiple variables, notably: the phylogenetic distance (=age) of the TE, the number of observed dimers (=frequency), their degradation rates (=milliDiv), as well as the binding affinity of a given dimer (=minimal free energy, MFE). However, after having received the Reviewers’ comments and upon reexamining our original manuscript version, we realised that we had indeed not made this sufficiently clear in text and figures.

We have therefore made the following changes. First, we already clarified some of the conceptual points in Figure 4 (where we analyse parental vs. non-parental genes). We now make it clearer where and how age and binding affinity are assessed in Figure 4C-F (mouse) and Figure 4—figure supplement 2 (other species). In particular, it should now be clearer that Figure 4D refers to age and Figure 4E to binding energy. The former use of the phrase “binding score” was changed in the new version to “TE pairing score (age + MFE)”, to indicate that it integrates both the components age and binding energy (Figure 4F).

The major changes relating to Essential revision 4 we implemented in Figure 5, where shared vs. species-specific circRNAs are compared. Current Figure 5A and Figure 5-Figure supplement 1 now show the dimer count per gene on the x- and y-axis (the corresponding panel in the previous manuscript version was based on total count); the “per gene” analysis facilitates the interpretation of enrichment in shared vs. species-specific dimer counts. Furthermore, we have included a new supplementary figure – Figure 5—figure supplement 2 – to portray the increase of possible dimers that can be formed (“dimer interaction landscape”) due to increased number of dimers in the flanking introns of shared vs. species specific circRNA loci.

Next, we reconsidered how we could calculate more precisely the actual degradation rates of the relevant dimers; the outcome is found in the new Figure 5C and Figure 5-Figure supplement 3 (left panel). Briefly, we now show a comparison of milliDiv values between the individual predicted dimer-forming repeats – rather than the mean milliDiv values from repeats of the same dimer, as in the old version – should conceptually be a more precise approach. With this analysis, we now conclude that degradation levels between top-5 dimers do not differ between the flanking introns of shared and species-specific circRNA loci. Because degradation rates alone may not fully capture the actual decline in pairing of the dimer, we next added the analysis of the actual binding affinity for all top-5 dimers (Figure 5, right panel) calculated from the individual genomic sequences (rather than from the reference sequence). In the next analysis, we selected for each gene the least degraded dimer, i.e. the one that is strongest and that is likely to form (see Material and Methods). Again these analyses showed no significant difference between the flanking introns of shared and species-specific circRNA loci, fully in line with our model.

We believe that these changes have significantly improved our analysis. However, our deep analyses of dimer landscapes are rather novel; in this context, we would like to note that we are quite aware that dimer formation will be the result of a complex interplay between TE frequency, degradation levels, and sequence affinity – how these individual parameters should be weighted and how they can be integrated into a quantitative measure is not simple to evaluate. Nevertheless, all our analyses suggest that degradation rates and sequence affinities are fairly similar between shared and species-specific circRNA loci; we therefore also assume that TE frequency is one of the driving factors for circRNA formation.

Suggested text revisions:

a. The discussion should further elaborate the concept that at least for some circRNAs, their main functionality is to regulate/limit the expression of the host gene, hence the exact exons that form circRNAs might not be important. Therefore it would not so important which exon in particular makes the cirRNA or even which elements in the flanking introns are driving this process but rather the existence of this during evolution. This is relevant especially for those circRNAs for which the circular molecules are abundant in comparison with the linear mRNAs generated from the locus and there are at least some cases in the literature in which has been shown that production of the circRNA comes at the expense of the linear transcript.

We now more clearly include this point towards the end of the discussion.

b. Discuss the evolutionary implications of the observation that circRNAs in fruitfly is not driven by intronic pairing of repetitive/transposable elements (Westholm, et al., Cell Rep 2014, PMID: 25544350).

We have looked at the cited publication in some detail in order to see how we could best accommodate this suggested text revision. First, we noticed that the publication is relatively careful in its wording and does not explicitly state that “circRNAs in fruitfly is not driven by intronic pairing of repetitive/transposable elements”. A telling statement from the paper is, for example: “However, overall, it appears that Drosophila RNA circularization does not appear to be driven by flanking sequence or structural complementarity, as in mammals.” - the authors are thus very careful in their wording in order to avoid overstating this particular point. Closer examination of their data shows, for example, in Figure S3A that their approach using MEME recovers one short substring of the ALU element described by Jeck et al., (2012) in human circRNA flanking introns (although, of note, this Figure unfortunately does not show the same analysis in control introns). In Figure S3B, a comparison of Drosophila introns (500 nt window) recovers some repetitive sequences that are of lower complexity and found in both circRNA and non-circRNA introns, and that are therefore dismissed as possible transposable element substrings. Still, their Figure S4B actually does show that circRNA flanking introns contain significantly better stem-loop-potential than non-circRNA introns, especially for short sequence window sizes (interestingly, the effect is clear until window sizes of 100 nt, but absent when using a 500 nt window, possibly indicating that the MEME analysis could have resulted in other outcomes when using smaller windows as well). Westholm et al., do not extensively develop the lack of evidence for intron pairing or repetitive/transposable elements in their study, but concentrate on other main findings. As it stands and only based on this study, therefore, we feel it would be misleading to develop and discuss evolutionary scenarios in our mammalian-centric study.

The situation is also further complicated by the fact that Kramer et al., (2015) (PMID: 26450910) actually have suggested a role for intronic repeats in circRNA biogenesis in Drosophila. Also Ashwal-Fluss et al., (2014) (PMID 25242144) identify numerous reverse complementary repeats between the two flanking introns of the mbl circRNA (see Figure S2B in their publication). Despite the specific role of MBL in this circularisation event, this finding suggests that sequence-dependent pre-mRNA folding may have an additional, and possibly more ancestral, role in circRNA biogenesis in this particular example (reviewed in Patop et al., (2019) PMID: 31343080). Finally, there is good evidence for reverse-complement sequence mediated circRNA biogenesis in C. elegans (Ivanov et al., 2015; PMID 25558066). With all these caveats to be navigated – and with the idea of simplifying/shortening our discussion in view of other reviewers’ comments, we hope the editor and reviewers will accept that we prefer not to dedicate a part of the discussion to conservation of circRNA biogenesis beyond the subject of our mammal-based study.

c. On page 6, the authors stated "Coding genes in rhesus macaque and human are characterised by a bimodal GC content distribution". The H3 category contains all genes with GC content > 52% while L2 category contains genes with 37-42% GC contents (Figure 3A). It is not correct to define it as bimodal GC content distribution when the group standard has a large difference.

This has been changed.

Former sentences:

“Coding genes in rhesus macaque and human are characterised by a bimodal GC content distribution (see peaks in L2 and H3 for non-parental genes). By contrast, the two rodents displayed a unimodal distribution (peak in H1), whereas opossum coding genes were generally GC-poor (in agreement with Galtier and Mouchiroud, 1998; Mikkelsen et al., 2007).”

have thus been changed in the revised version to:

“Non-parental genes displayed a unimodal distribution in the two rodents (peak in H1), were generally GC-poor in opossum (peak in L1), and showed a more complex isochor structure in rhesus macaque and human (peaks in L2 and H3), in agreement with previous findings (Galtier and Mouchiroud, 1998).”

d. More data are needed to sustain the statement "circRNA parental genes are characterised by low GC content and high sequence repetitiveness". Low GC content and high sequence repetitiveness may be the consequences of the long intron length of parental genes and/or the enrichment of repeat elements in long introns. More stringent conditions, like comparable intron lengths between circRNA and non-circRNA parent genes would need to be set up for this type of comparison.

To address this “text revision”, we have now rephrased the cited statement, e.g. in the subsection title (now: “CircRNA parental genes are associated with low GC content and high sequence repetitiveness.”).

In terms of the need for additional data to support the statement, we would like to point out that in particular the modeling approach that we applied actually is designed to isolate the individual predictors, i.e. it calculates the association of, for example, GC content, over the other effects, such as intron length. We are therefore not fully sure what is meant with the comment that “more stringent conditions are needed”.

Reviewer #1 (Recommendations for the authors):

I would recommend the authors to consider the following revisions:

1. There seems to be some contradiction between the argument that circRNAs tend to be generated from constitutive, widely expressed exons and the statement on page 3 that circRNAs "showed considerable tissue-specificity".

We have now rephrased this sentence to be clearer:

“Identified circRNAs were generally small in size, overlapped with protein-coding exons, frequently detectable only in one of the tissues, and were flanked by long introns (Figure 1—figure supplement 3).”

This wording is better compatible with the main underlying reason for “tissue-specificity”, namely that the stochastic splicing noise around the constitutive exons will not lead to reproducible circle formation across all tissues. Second, similar to linear alternative splicing, one needs to differentiate between the transcript and exon level: alternative RNA isoforms can show tissue-specific expression patterns, while being partially composed of the same (constitutive) exons – the exon itself can be considered expressed in multiple tissues (because it is used in multiple, but different isoforms), yet the transcript is not.

2. The authors may want to elaborate further on the possible links between the preferred genomic architectures of circRNAs hotspots and the concept of exon definition (e.g. enrichment in long flanking introns and differences in GC content between exons and introns, as reported by Ast and colleagues).

We agree that this would be a very interesting idea to address in the future. However, we felt that in its current state the manuscript was already quite loaded with information analyses, including comments that we should shorten the discussion. We therefore did not extend our analyses in this additional direction at this point. We hope the Reviewer will find this decision acceptable.

3. The predictive model generated by the authors produces very impressive results. In addition to sequence complementarity, a number of protein factors have been shown to support the formation of stem-loop structures associated with back splicing, including hnRNP L, MBNL1 or PTB. The authors may be in an excellent position to assess whether the presence of potential binding sites for these proteins contributes evolutionarily to the generation of circRNAs (or subtypes of them). More specifically, does the inclusion of binding sites for these factors improve the performance of their predictive model?

The reviewer raises an excellent point and opportunity for further in-depth analyses.

It is indeed possible that circRNAs become first expressed by the formation of hairpin-like structures driven by TEs in their flanking introns. In a second step (and in case that the circRNA has conferred a benefit for the organism), circRNA expression may be stabilised by the evolution of new binding sites for proteins such as hnRNP proteins, MBNL1 and PTB. At the same time, this reduces the selective pressure on the TE, leading to its degradation and circRNA candidates without detectable dimers. While we focus on the “first steps” towards the expression of a circRNA, the integration of binding site information from external datasets to understand “what happens next” can be considered an excellent starting point for a follow up study. Yes, as for the previous point, we found this interesting request too complex to accommodate at this point, and thus beyond the scope of this first study. We do, however, expand a bit more on RBPs in the discussion.

4. Figure 4A: it would be good to provide some measurement of the statistical relevance of the enrichment scores observed.

Measurements of statistical relevance have been integrated in the figure and the statistical tests used are mentioned in the figure legend.

5. Figure 5 and related discussion about evolutionary age of TEs associated with circRNA biogenesis: is it possible that dimers co-evolving over long periods of time accumulate mutations but retain -through compensatory changes- their capacity to form base pairing interactions leading to looping out and back-splicing?

The Reviewer raises an interesting point. Our modified Figure 5 (see above, Essential revision 4) now looks at aspects of this idea by calculating the actual binding energies of the local TE copies predicted to form dimers.

In general, however, we would think that most sequence changes – even the compensatory ones – are rather random in nature than selected for. The most relevant TEs that we identify are SINE elements, which are rather GC-rich, but integrate into GC-poor environments. Our analyses (not shown in the manuscript) indicate that many of the mutations we observe are G/C -> A/T changes, for which we would tend to hypothesize that they occur as a consequence of the GC-poor environment to which the integrated TE is adapting. This GC adaptation to the local GC content will happen in both TEs of the dimer in parallel, which can lead on average to more seemingly “compensatory” mutations than expected. However, these compensatory changes would not be driven by selective pressure, but by the local GC content. Therefore, although it is possible that a few of them are real compensatory mutations, it would represent a major project to properly identify, classify and analyse these signals.

Reviewer #2 (Recommendations for the authors):

[…]

The regulation of circRNA biogenesis can be complex. It should be careful to evaluate different conditions to draw a convincing conclusion. The major pitfall of this study is to use relatively small population of circRNAs for the evolutionary study, which may lead to biased findings. In addition, more stringent controls are required for analyses throughout the study. The Discussion part seems to be unusually long, containing lots of assumptions that are needed to be proven. Please shorten.

The first part of Reviewer #2’s recommendations (“The major pitfall of this study is to use relatively small population of circRNAs […]”) is largely overlapping with Essential revision 1 at the beginning of this rebuttal document. Please refer to this section for a detailed response.

The context of the second statement “In addition, more stringent controls are required for analyses throughout the study.“ is not fully clear from the above paragraph. However, in the more detailed “Public evaluation summary”, we find as point 4 the statement: “More data are needed for the statement "circRNA parental genes are characterised by low GC content and high sequence repetitiveness". Low GC content and high sequence repetitiveness may be the consequences of the long intron length of parental genes and/or the enrichment of repeat elements in long introns. More stringent conditions, like comparable intron lengths between circRNA and non-circRNA parent genes, should be set up for this type of comparison.”

We agree with the Reviewer that it is important to account for a possible dependence of variables, as shown by the example that the Reviewer cites: GC content and sequence repetitiveness will clearly not be independent from intron length. This is precisely why we used generalised linear models, which is the state-of-the-art approach when multiple correlated dependent variables are predicted. To avoid possible biases by correlated variables, we have also assessed all variables for multicollinearity before including them into the analysis (see Material and Methods, section Generalised linear models). We prefer this approach of a global model that includes all potential predictive parameters, over the Reviewer’s suggestion of individual comparisons of single features, precisely because such comparisons are prone to biases like the one described by the Reviewer.

The third point – shortening of discussion and removal of assumptions – is well taken. We have substantially restructured and rewritten the discussion accordingly.

Reviewer #3 (Recommendations for the authors):

[…]

The manuscript is interesting and timely. However, the conclusions are too broad and the specific issues need to be addressed experimentally and/or by doing additional analysis.

More specifically:

– The approach for circRNA detection is at least problematic. Without a mock control is impossible for the authors to determine the actual resistance to the degradation by the exonuclease and hence determine which potential backsplicing junctions are really coming from circular molecules. One thing the authors could do is at least for some of the tissues/species show that the circRNAs they identified have been previously showed to be real (e.g., by a canonical comparison mock vs RNAseR treated samples).

This point has been answered under Essential Revision 3.

– Authors seems to ignore that in vivo studies focused on highly expressed circRNAs. Can the authors partition their assessment of conservation for the different quantiles of expression? Indeed, I would believe that only circRNAs in the top quantile of expression are the ones that will really give insights into the potential functional importance of circRNAs.

This point has been answered under Essential Revision 2.

– Notion of noise for lowly expressed circRNAs is prevalent in the field and the key question is whether the circRNAs expressed in the neural system at high levels have function. So only on those circRNAs is relevant to determine evolutionary conservation. Indeed, the authors could take the opposite approach; for those few circRNAs for which some functionality been proven, are the pattern of evolution and elements flanking the circularizable exons similar to the ones they describe generally?

The largest part of this point, too, is picked up by several of the Essential revisions - please refer to the beginning of this document. Moreover, please see Figure 3-Figure supplement 3 for an analysis of known, functional circRNAs.

– The manuscript ignores potential cis regulatory mechanisms of circRNA biogenesis. One could argue that an important feature of a given locus is to produce circRNAs as a way to limit the amount of linear RNA. In this sense, it would not so important which exon in particular makes the cirRNA or even which elements in the flanking introns are driving this process but rather the existence of this during elevation. This is relevant especially for those circRNAs for which the circular molecules are abundant in comparison with the linear mRNAs generated from the locus and there are at least some cases in the literature in which has been shown that production of the circRNA comes at the expense of the linear transcript. The authors could include this also as a criterion for evolution and not rule it out based on some theoretical assumptions.

We now dedicate more space to this important concept in the discussion of the revised manuscript (see also above our response to suggested text revisions).