Background {#Sec1}
==========

Honeybees (*Apis mellifera*) display a striking behavioral plasticity among their workers that is reflected in an age-dependent task allocation and thus represent a substantial model organism for phenotypic plasticity. Workers are able to execute varying specific behaviors in order to fulfill tasks that are essential for the viability of the colony, such as cleaning combs, feeding the larvae, guarding the nest entrance and foraging for food. To ensure economic efficiency and to prevent randomly performed tasks, the assignment of tasks has to be coordinated \[[@CR1]\]. Task allocation is predominantly dependent on the age of the worker bees, but is also flexible and can be adjusted to colony needs \[[@CR2]--[@CR4]\]. In experimental single cohort colonies (SCCs) that are solely composed of young bees, some colony members initiate foraging precociously irrespective of their age \[[@CR5]\]. Major differences in task-related behaviors exist between the typically younger nurse bees that feed the larvae inside the hive and the older foragers (≥18 days after emergence) that leave the hive to collect pollen, nectar and water \[[@CR6]--[@CR9]\]. This phenotypic plasticity is also reflected at the neuronal level. The overall brain volume is increased in forager bees compared to nurse bees \[[@CR10], [@CR11]\] especially in visually innervated brain structures \[[@CR12], [@CR13]\]. At synaptic levels, these changes involve for instance the density of synaptic complexes within mushroom body calyces caused by the growth of Kenyon cell dendrites and pruning of presynaptic boutons \[[@CR14]--[@CR16]\]. The regulation of these processes is poorly understood and seems to be highly complex. Various effectors are known which include the external environment, the colony state and internal stimuli such as (post-)transcriptional changes. Alterations in the expression ratio of hundreds of genes were detected, including some with synaptic functions \[[@CR17]--[@CR20]\]. Additionally, protein expression is affected as shown for instance for changes in neuropeptides \[[@CR21]\], the membrane proteome \[[@CR22]\] and the phosphoproteome \[[@CR23]\] in the worker bee brain, but also in peripheral tissues such as the hypopharyngeal gland that produces royal jelly proteins \[[@CR24]\]. Task- or age-related differences were also observed in the abundance of micro RNAs (miRNAs). Many of the identified miRNAs have a number of putative target genes that also exhibit functions in a neural context \[[@CR25]--[@CR27]\].

Circular RNAs (circRNAs) represent a class of RNA with considerable regulatory potential that was overlooked for decades and is currently under extensive research and discussion. An increasing number of studies show that circRNAs are abundant, differentially expressed and even have biological functions \[[@CR28], [@CR29]\]. In general, circRNAs arise from a back-splicing event. The 5'-end of a donor exon is joined to a 3'-end of an acceptor exon of the same molecule \[[@CR30]\]. This results in a so-called back-spliced junction (BSJ) which can be observed as junction-spanning reads (JSR) when mapping RNA-Seq data to a genome using a mapper that supports split reads. The abundance of circRNAs typically varies between tissues and is sometimes uncorrelated to the host mRNAs \[[@CR31]\]. This may indicate a specific function of circRNAs but might as well reflect distinct decay rates for linear compared to circular transcripts which lack accessible ends. Studies point out that circRNAs may act as regulators of alternative splicing \[[@CR32]\] or could feature miRNA sponges \[[@CR28], [@CR33]\].

Besides human and mice, the presence of circRNAs was verified and studied extensively in the fruit fly *Drosophila melanogaster* and recently in the silkworm *Bombyx mori* but no other insect so far \[[@CR28], [@CR30], [@CR31], [@CR34]\]. Important findings are the presence of numerous miRNA seed matches in line with a putative miRNA sponge function as well as the fact that circRNAs mainly derive from neural genes and accumulate in neural tissues in an age-depended manner \[[@CR32], [@CR35], [@CR36]\]. Following up on these findings, circular RNAs may contribute to regulating the age-related transition from nurse bees to foragers at the molecular level.

Results {#Sec2}
=======

Identification of circRNAs in the brain of honeybees {#Sec3}
----------------------------------------------------

As circRNAs do not feature 5'- or 3'-ends they are virtually resistant to RNase R treatment, which digests most linear RNAs. The enzyme can thus be used to enrich total RNA extracts for circRNAs \[[@CR37], [@CR38]\]. In order to identify these, we prepared RNA-Seq libraries from total RNA extracts of honeybee worker brains. The libraries were enriched for circular RNAs and compared to a non-enriched library. Each BSJ was considered as representative of a distinct circular RNA. We were able to detect a total of 3384 individual BSJs supported by at least three JSRs from the four libraries combining two different methods, see Fig. [1](#Fig1){ref-type="fig"}a. Based on these we provide two sets of circRNAs identified by applying different stringency thresholds (see Material and Methods for details). The low stringency set contains 1263 circRNAs found by both independent algorithmic methods (overlap). Only these BSJs were considered viable circRNA candidates because previous studies showed inconsistent results between different algorithms \[[@CR39], [@CR40]\]. Specifically, segemehl is known to produce very sensitive mapping results, potentially introducing false positives when solely relied upon. The high stringency set used for the following analyses is a subset containing 254 circRNAs with a higher amount of supporting reads along with a significant five-fold enrichment of the JSRs through RNase R treatment. The majority of the circular transcripts were even enriched by more than ten-fold (\>77*%*). Fig. 1Identified circRNAs by RNA-Seq. **a** Two independent algorithms were used to predict circRNAs: segemehl in conjunction with testrealign (green) reported 2725 BSJs while BWA with CIRI2 reported 1922 BSJs. The overlap was considered as low stringency set (yellow). The high stringency set additionally requires an enrichment through RNase R treatment and compelling read coverage from at least two independent sequencing libraries. See Materials and Methods for details. These 254 candidates are most likely genuine circRNAs and are therefore used for further analysis. **b** Strong evidence is found in three independent sequencing libraries for 139 candidates. The remaining 115 are well supported by two samples. Note that all candidates were also identified as enriched in *E*~+~ vs *E*~−~. **c** 122 host genes are orthologous to host genes of circRNAs identified in either *Drosophila* or *Bombyx* in previous studies

A vast majority of BSJs was flanked by a canonical GT/AG splice signal. Only five circRNAs did not show such a canonical splice site (see Additional file [1](#MOESM1){ref-type="media"}). In one case an annotation was not possible. The BSJ spans two exons that are (presumably) not spliced together. The coding exon of gene CG45167 (homolog of B52 in *D. melanogaster*) and its immediate downstream exon which starts with the 5'-UTR are not present in any currently annotated transcript variant. Details are illustrated in Additional file [1](#MOESM1){ref-type="media"}.

We picked some of the most significant circRNAs that were highly abundant or showed a particularly differential expression pattern between nurse or forager bee libraries. Presence and circularity of these selected circRNAs were verified further by additional PCR experiments, see Additional file [2](#MOESM2){ref-type="media"}. TaqMan based Real-Time PCR assays were used to examine the expression levels in nurse or forager bee for two salient circRNAs in an independent experimental approach, see "Quantification in nurses and foragers". A complete list of all 254 high confidence circRNAs including read levels and putative homologs in *D. melanogaster* and *B. mori* can be found in Additional file [3](#MOESM3){ref-type="media"}. An excerpt of the most prominent entities is shown in Table [1](#Tab1){ref-type="table"}. Table 1Excerpt of identified circRNAs in the brain of honeybee nurse and forager beescircRNA IDHost geneBeeBaseChr.JSRsEnrichedExpr.Homologyame_circ_0001970LOC413427GB43145LG114326.70.329⋆∘ame_circ_0000721LOC724885GB53835LG313910.20.344ame_circ_0002142LOC410393GB52063LG121245.00.628*†*ame_circ_0000163LOC408576GB42249LG11037.20.328⋆∘ame_circ_0000232Mup2GB49259LG19711.00.119ame_circ_0001780radGB49511LG10915.00.062⋆∘ame_circ_0001286LOC411534GB44365LG76393.00.731⋆∘ame_circ_0002579LOC409655GB47584LG164217.50.226ame_circ_0001822Rsmep2GB54272LG10345.00.022⋆∘ame_circ_0002577LOC409655GB47584LG161711.30.072ame_circ_0001852CoRestGB52614LG10105.00.013⋆∘ame_circ_0001099LOC411114GB44582LG530618.50.328*†*ame_circ_0000414LOC725294GB55364LG22167.60.312⋆∘ame_circ_0000397LOC408688GB49767LG21855.00.377⋆∘*†*ame_circ_0001712LOC408996GB42579LG91696.30.198⋆∘*†*ame_circ_0002576LOC409655GB47584LG1616814.90.644ame_circ_0001638LOC411347GB17597LG91599.40.400ame_circ_0001593LOC408991GB53310LG91487.90.105ame_circ_0001479LOC408957GB40504LG814718.40.339∘ame_circ_0000524LOC408718GB43446LG213036.20.313ame_circ_0001120sGC-alpha1GB52929LG612910.40.276*†*ame_circ_0000054LOC726544GB42188LG11247.50.480ame_circ_0001877LOC408309GB45167LG111219.40.085ame_circ_0000669LOC410044GB55791LG311813.00.370ame_circ_0000073LOC410717GB55293LG111111.00.290⋆∘ame_circ_0001340LOC411229GB42567LG71096.90.570∘All circRNAs were significantly enriched in *E*~+~ over the non-enriched set *E*~−~. *Set* refers to the RNA-Seq libraries in which the circRNA was enriched in addition (see Table [2](#Tab2){ref-type="table"}). The *host gene* is given according to the RefSeq GCF_000002195.4 annotation along with the corresponding *BeeBase* identifier \[[@CR74]\]. The respective chromosome is indicated in the *Chr.* column. The summarized number of *JSR*s is given along with the averaged normalized expression levels relative to the host gene expression *expr.* and fold enrichment *enriched*. The *homology* column indicates whether a *Drosophila* or *Bombyx* homolog was found in ⋆\[[@CR32]\] ∘\[[@CR35]\] or *†*\[[@CR34]\]. The first block corresponds to circRNAs that were particularly strongly expressed (many JSRs or high norm. expression) or showed signs of differential expression between nurse and forager bee libraries and were thus selected for verification of circularity and presence in further PCR experiments. The full list can be found in the Additional file [3](#MOESM3){ref-type="media"}

Homologs to fly and silkworm {#Sec4}
----------------------------

Honeybee circRNAs were compared to those found in fruit fly \[[@CR32], [@CR35]\] and silkworm \[[@CR34]\] based on the homology of their parental genes (Fig. [1](#Fig1){ref-type="fig"}c). Out of 254 honeybee circRNAs only 70 host gene homologs were found in silkworm (30%). In contrast, 203 homologous BSJs were identified for fruit fly (80%) which can be explained by the closer phylogenetic relationship to honeybee \[[@CR41]\]. Consistent with our results, circularized exons in fruit fly were found in 144 to 151 of these homologs (with respect to \[[@CR32], [@CR35]\], overlap 122 circRNAs). This finding is in line with a similar comparison of circRNAs in human and mouse. There, two-thirds of all host genes harboring back-splicing junctions could be correlated by homologies between the two species \[[@CR42]\]. A complete listing of the results can be found in Additional file [3](#MOESM3){ref-type="media"}.

GO term enrichment {#Sec5}
------------------

A GO term analysis (gene ontology term enrichment) was performed using all 203 circRNA host gene homologs correlated to fruit fly from which we extrapolated the functional annotation. High-level processes involved in synaptic development and regulation were significantly enriched. Given that the source samples were obtained from brain tissue, this is an expected result but it also resembles the finding that neurologically associated genes are a main source of circRNAs as found for *D. melanogaster* \[[@CR35]\]. The most enriched high level terms below a *p*-value of 10^−4^ were "anesthesia-resistant, medium- and long-term memory" (27 ×) "medium-term memory" (23 ×), "regulation of neuromuscular synaptic transmission" (21 ×) and "deactivation of rhodopsin mediated signaling" (21 ×). The former is especially remarkable. One representative of this group is the *radish* gene from which circAm*rad* (*ame_circ_0001780*) arises. We found that the abundance levels of circAm*rad* correlated with the acquired task of a bee (see "Abundance and task allocation" below). Consistent with this is also the enrichment of rhodopsin signaling and memory-related genes. Nurse bees take care of the brood inside the hive, where it is dark and the requirements to memory are different from those of foragers \[[@CR43]\]. After task transition to forager bees, they start to collect food from outside the hive, mostly at daylight, and need to find their way back to the hive afterward. A need for adaptation of rhodopsin signaling and a change in memory requirements is obvious. In fact, "positive phototaxis" showed the highest GO term enrichment (44 ×). The significance (*p*=1.87×10^−3^) however was above the applied threshold because the term only has four representatives in the reference set. A detailed overview of enriched GO terms can be found in Additional file [4](#MOESM4){ref-type="media"}.

Exon-intron structures {#Sec6}
----------------------

The majority of BSJs in honeybee corresponds exactly to exon boundaries of protein-coding regions (78%), see Fig. [2](#Fig2){ref-type="fig"}. Nearly all remaining cases are derived from 5'-UTR containing segments (17%). This is only slightly different from the set of (presumably) linearly spliced exons in the control but shows a trend towards 5'-UTRs. For both *D. melanogaster* datasets \[[@CR32], [@CR35]\] the overall proportion is similar but with a much stronger bias towards 5'-UTRs (∼30*%*) and non-canonical splice events, e.g. occurring in the middle of introns or exons in between genes (∼20*%*, other). The latter category was rarely found for honeybee circRNAs (\<2*%*). We note that this difference might be a result of different annotation qualities for honeybee (data from 2018) and fruit fly (data from before 2014) and should thus not be over-interpreted. Fig. 2Location of circRNAs in the original transcripts. BSJs are clustered into the following categories: part of the 5'-UTR, from the coding sequence exclusively, part of the 3'-UTR, spanning from 5' to 3'-UTR or other (in the middle of exons, introns, part of non-coding genes or in between genes). Most of the circRNAs in *A. mellifera* originate from exons in the coding portion of a transcript, which also represents the largest amount of exons transcriptome-wide. Compared to exons of a random control, circRNAs originate slightly more often from 5'-UTR exons as also described for fruit fly. However, we did not find as many BSJ including exons of the 3'-UTR as reported there. All high confidence circRNAs were mapped to annotated splice sites and we could not detect BSJ with sufficient read numbers outside of annotated regions

For fruit fly, it was reported that circRNAs mostly originate from the second exon of a transcript \[[@CR35]\]. This is also true for honeybee circRNAs. Figure [3](#Fig3){ref-type="fig"}, however, shows that this number is implied by the outstanding abundance of transcripts with only two exons. This is also visible in the randomized control distribution. Compared to this set, the observed starts at exon two are actually less than what would be expected. We identified two factors that correlate with back-splicing: The exon position and the number of exons. The further downstream an exon is located in a transcript and the more exons (and thereby splice junctions) it exhibits, the more likely circRNAs arise from the transcript. Fig. 3Exon position of circRNAs in the host gene. Exons of the parent transcripts were numbered. The number of the first an the last exon included in a circRNA is represented in the left (5'-end) and right (3'-end) panel, respectively. CircRNAs are stratified by (**a**) normalized circRNA expression into four groups (green to red with decreasing ratio) and (**b**) relative fold change between RNase R untreated and treated experiments (green to red with decreasing difference). As a reference, randomly selected exons in the control (gray) exhibit a higher number of second and third exons in a transcript just by chance

Another finding from fruit fly indicates that circRNAs with higher normalized expression tend to favor earlier exons than less expressed variants \[[@CR35]\]. We reproduced this by partitioning the BSJs according to their normalized expression levels, see Fig. [3](#Fig3){ref-type="fig"}a. A similar trend is visible in our data. However, the shift towards later exons for less expressed circRNAs (e.g. with relative expression levels \<0.05) is not as pronounced. Notably, our control exons exhibit a much stronger bias towards the second and third exon for the starts of circular junctions than any of the partitions with circularized exons (almost 60%), especially compared to the control used for fruit fly \[[@CR35]\]. An alternative stratification of BSJs by their relative fold change in RNA-Seq libraries enriched for circRNAs yields the same results, see Fig. [3](#Fig3){ref-type="fig"}b. The circRNAs presented here do not involve parts of the 5'-UTR more often than expected but transcripts with unusually long 5'-UTRs appear to be prone to circularization at an increased likelihood.

Intronic features {#Sec7}
-----------------

In honeybee, introns flanking circularized exons are significantly longer than those from linearly spliced exons, see Fig. [4](#Fig4){ref-type="fig"}a. They can span several thousand bases. This result is in line with findings from fruit fly and human \[[@CR35], [@CR44]\]. There, in addition, flanking introns showed increased levels of reverse complementarity compared to linearly spliced exons. Reverse complementary regions are thought to enhance the likelihood for base-pairing between the introns. This interaction likely guides back-splicing process \[[@CR36], [@CR45], [@CR46]\]. Following up on this assumption, introns were reciprocally scanned for reverse complementary matches at sequence-level using BLAST \[[@CR47]\], see Fig. [4](#Fig4){ref-type="fig"}b. While the result shows that introns flanking circularized exons are composed of regions with better complementary (represented by higher bitscores) in general, it is also obvious that complementarity is linked to the length of introns. Higher scores of complementarity matches are likely a result of the fact that introns flanking circularized exons are much longer than those from the control set. The most relevant regions for circularization are probably the end of the 5' flanking and the start of the 3' flanking intron, see Fig. [5](#Fig5){ref-type="fig"}a for a scheme. Even if the comparison is limited to these regions, the difference in complementary matches cannot explain why some exons are circularized and others are not. The median complementarity is about equal to the control introns that flank linear exons even though the latter show much higher variance especially towards introns with hardly any complementarity, see Fig. [5](#Fig5){ref-type="fig"}b. Fig. 4Introns flanking circRNAs. **a** Flanking introns upstream (left) and downstream (right) of circRNAs (yellow) are significantly longer (*t*-test with *p*\<0.001 for both) than those of a random control (gray). **b**BLAST bitscores of the best reverse complementary match of the shorter circRNA flanking intron to the other (yellow) compared to control intron pairs (gray). Introns flanking circRNAs are visibly longer than random introns but show similar bitscores for reverse complementary matches of equal intron lengths (correlation *ρ*=0.5). Notably, shorter introns in the control frequently exhibit relatively high bitscores, indicating that short complementarity regions are common among intron pairs Fig. 5Putative interactions and sequence properties of flanking introns. **a** Scheme indicating the relevant regions (red) of introns flanking circularized exons (yellow). **b** Screen of sequence complementarity using BLAST. The scores for circular and linear control intron pairs largely overlap. **c** Screen of intron interactions represented by the MFE of constrained RNAcofold calculations. The difference between both groups is significant. Still there is some overlap. **d** Relation of GC-content and MFEs. Both measures appear to be linked (correlation *ρ*=0.5, Spearman's rank correlation). **e** Screen of GC-content. A significant increase is found compared to the control. **f** Screen of CpG dinucleotides. A significant increase is found compared to the control which rarely exhibits CpG. **g** Screen of cytosine-content. A significant increase is found compared to the control. It explains the observed effects for GC-content and CpG frequency. **h** Screen of average base methylation. The number of methylated bases in introns flanking circRNAs is slightly but significantly increased compared to linear control introns. \**p*-value \<0.001

An RNA secondary structure prediction using RNAfold \[[@CR48]\] was used to investigate potential intron-intron interactions more specifically, see Fig. [5](#Fig5){ref-type="fig"}c. The difference is more obvious using this method. Cofolded complexes of the control introns exhibit much higher minimum free energy scores (MFE), indicating less base-pairing interaction. The difference is highly significant (*p*\<0.001). However, the MFE scores partly cover similar ranges, which does not allow for a clear distinction between circularized exons and linear splicing products. Figure [5](#Fig5){ref-type="fig"}d shows that the increase in folding-potential (represented by lower MFE scores) is linked to GC-content of the respective introns. Also the fact that the complementarity match, as well as the cofolding analysis, yielded similar results for all combinations of starts and ends of the flanking introns (e.g. pairing the end of the upstream intron with the end of the downstream intron) puts a direct effect of base-pairing in doubt. The GC-content in turn well discriminates circRNA introns from control introns, see Fig. [5](#Fig5){ref-type="fig"}e.

Methylation {#Sec8}
-----------

The intronic features raise the question, why the GC-content of circRNA flanking introns is elevated in such significant amounts (median shifted from 20 to 36%, *p*\<0.001). One reasonable explanation is an increase of potential DNA-methylation at these introns (CpG islands). Figure [5](#Fig5){ref-type="fig"}f illustrates that the CpG dinucleotide frequency is also significantly increased for circRNA flanking introns and nearly absent in the control group (∼1*%*). As CpG sites are preferentially methylated \[[@CR49], [@CR50]\], this indicates a significant increase of potential DNA-methylation sites. Moreover, cytosine methylation and hydroxymethylation at non-CG sites (CA, CT, CC) are reported to be enriched in introns of the honeybee \[[@CR51]\]. In line with this, Fig. [5](#Fig5){ref-type="fig"}f shows that also the cytosine mononucleotide frequency is significantly increased for circRNA flanking introns. While the genome comprises ∼16*%* cytosines, circRNA introns exhibit a median of ∼18*%* cytosines. Strikingly, the median cytosine-content of introns in the control group is as low as 10%. This can be translated into reduced methylation and hydroxymethylation potential and thereby fewer alternative splicing events for introns flanking canonically spliced RNAs compared to those that frequently result in circRNAs.

We evaluated publicly available whole-genome bisulfite sequencing data of worker bees from a previous study to comprehensively determine methylation levels \[[@CR52]\]. Figure [5](#Fig5){ref-type="fig"}h shows that the length-normalized accumulative DNA-methylation of introns flanking circular RNAs actually tends to be increased compared to those flanking random exons. Notably, the effect was not visible using only the closest 50 or 100 nucleotides of a flanking intron but became visible using a 200 nt window or full-length introns. This is probably due to the limited windows size which is likely too small for statistical assessment.

miRNA targets {#Sec9}
-------------

Potential miRNA target sites were annotated for all 254 circRNAs identified here. The results can be divided based on their degree of phylogenetic conservation. 3058 target sites were only conserved in *Apis* species. We argue that *Apis* species are too closely related to qualify as a reliable predictor for miRNA target sites. The sequence conservation in this set appears rather high in general. This is also reflected by a similar distribution of potential miRNA target sites compared to the control without any constraints on conservation, see Fig. [6](#Fig6){ref-type="fig"}a. Fig. 6Occurrence of putative miRNA target sites in circRNA exons. **a** Putative target sites normalized by exon length in differentially conserved sets. "All targets" represents no conservation at all. *Apis* represents conservation only in closely related species. Eusocial insects are sufficiently distant to consider conservation relevant. The number of miRNA binding sites conserved in this set is significantly higher (*p*\<0.001, Students *t*-test) than in the control. The absolute number of potential binding sites in the respective sets is given in parentheses below. **b** Illustration of miRNAs with potential target sites in circRNAs conserved even in eusocial insects. Only miRNAs with at least ten target sites and an at least five-fold enrichment over control are shown. Error bars indicate standard deviation in fold enrichment between different controls (*N*=42)

A set of 1076 sites is conserved in *Apis* and eusocial insects which are sufficiently distant to *A. mellifera* to reasonably infer conservation. With about 10.4 target sites per 1000 nt circRNAs have a 1.7 × increase in conserved, putative miRNA target sites compared to the median of the linear splice product control. Thus, in line with previous findings for *Drosophila* \[[@CR35]\], we report a general enrichment of conserved miRNA target sites in circRNAs over random linear counterparts. The most enriched miRNA target sites correspond to ame-miR-3748/ame-miR-3753 (∼10× enriched, same seed region) and ame-miR-3791 (∼9.2× enriched), see Fig. [6](#Fig6){ref-type="fig"}b. RNA expression studies show that the abundance levels of some miRNAs correlate with task or age of honeybees \[[@CR25]--[@CR27], [@CR53]\]. We did, however, not find a significant overlap of miRNAs corresponding to enriched target sites and miRNAs reported as differentially expressed in nurses and foragers. The complete list of potential target sites and their degree of conservation can be found in Additional file [5](#MOESM5){ref-type="media"}.

Quantification in nurses and foragers {#Sec10}
-------------------------------------

The circRNAs *ame_circ_0001780* and *ame_circ_0001822* showed a notable differential expression pattern in RNA-Seq results of nurse bees and foragers. For simplicity they will be termed according to their host genes in the further course of the study: circAm*rad* and circAm*rsmep2*, respectively. As the experimental setup is not suitable for any reliable quantitative assertions, we decided to perform a targeted quantitative Real-Time PCR for these candidates at different developmental stages. In addition, we compared the expression patterns in bees with age-related task allocation to those undergoing a task allocation due to colony needs (same-age, SCC), see Fig. [7](#Fig7){ref-type="fig"}. Fig. 7Expression analysis. TaqMen expression analysis of (**a**) circAm*rsmep2* and (**b**) circAm*rad* in brains of nurse bees and foragers from colonies with a typical age structure and SCCs consisting of bees of the same age. The expression is given relative to nurse bees. Bars show relative mean expression with standard error. The number of biological replicates is indicated in the bars. Significant differences are indicated (\* *p*\<0.05, \*\* *p*\<0.01, Two-tailed unpaired Students *t*-test). In contrast to circAm*rsmep2* the expression of circAm*rad* seems correlated with the allocated task rather than the bees' age

For circAm*rsmep2* we found that expression in the brain is higher in foragers than in nurse bees (Fig. [7](#Fig7){ref-type="fig"}a). This difference, however, does not seem to be directly task-related. In a SCC where the nurse and forager bees have exactly the same age, no expression differences are observed (Fig. [7](#Fig7){ref-type="fig"}a). Our interpretation is that this expression difference most likely depends on the bees' age but not on its task.

In contrast, circAm*rad* is higher expressed in brains of nurse bees than in brains of foragers (Fig. [7](#Fig7){ref-type="fig"}b, typical). Strikingly, this is inversely correlated with the expression of the linear product which is strongly increased in foragers (XM_393494.2/Am*rad*, *log*~2~ *ratio*∼6.1) \[[@CR18]\] and holds true independent of the age-related task transition. The expression levels in the SCC experiment (Fig. [7](#Fig7){ref-type="fig"}b, single cohort) are similar to that of typical colonies where tasks are allocated based on a bee's age. This data suggests a correlation of acquired task and circAm*rad* levels.

Discussion {#Sec11}
==========

We remark that all circRNAs reported here are expressed in the brain of nurse and forager bees. In contrast, 2513 circRNAs reported for *D. melanogaster* \[[@CR35]\] and 3916 for *B. mori* \[[@CR34]\] are based on samples of different developmental stages, tissues and even cultured cells and do not ensure RNase R enrichment. All analyses were performed using the high stringency set which requires a notable amount of supporting JSRs along with a significant enrichment through RNase R treatment. We argue that the enrichment control is necessary to discern genuine circRNAs from potential trans-splicing or exon-shuffling events. Otherwise independent experiments would be required to further support the sequencing-based evidence which is not feasible given the high number of involved loci. Focussing on this significantly expressed subset of circRNAs allows us to investigate genomic properties that are an inherent part of circRNA deriving loci. The inclusion of candidates with less confidence (less enrichment, fewer supporting JSRs) would introduce additional noise into statistical analyses. This observation was made e.g. regarding lower read numbers in a *D. melanogaster* study \[[@CR32]\]. For the majority of circRNAs identified here, the amount of canonically spliced transcripts (linear) is at least the same as the amount of back-spliced transcripts (circular). For this reason, it is unlikely that the circRNAs presented here arose from a mapping artifact, e.g. due to misalignment of reads or repeating gene copies.

Even though the exact mechanism is unknown so far, DNA-methylation is known to induce alternative splicing in honeybees \[[@CR54], [@CR55]\]. Our results indicate that this might affect the formation of circRNAs as well. A potential connection to the age and the allocated task is given by the fact that methylation patterns also vary depending on the age and allocated task of an individual bee \[[@CR49], [@CR50], [@CR56]\]. It was even shown that reverted nurse bees regain their original methylation patterns independent of their age \[[@CR52]\]. To follow up on this hypothesis, we correlated data from a methylation study \[[@CR52]\] with our annotation of circRNAs. While relevant social roles are comparable, we note that collection times and extraction methods differ from experiments done in this study. Ideally, the libraries used for circRNA detection and DNA-methylation analysis should be derived from the same biological sample. Without further experimental investigation, a strong conclusion cannot be drawn yet. We argue, however, that the data presented here provides first indications for a link of circularization and DNA-methylation in honeybees. On this basis, we speculate that the age-depended increase of circRNA abundance is not (only) due to potentially lower decay rates of circRNAs compared to linear products but also a result of increasing DNA-methylation that leads to alternative splicing accompanied by an increase of circRNA formation.

The quantification of circAm*rsmep2* in nurses and foragers indicates that its abundance most likely depends on the bees' age but not on its task. Supposedly, circAm*rsmep2* accumulates over time in the brain of worker bees, as shown for certain circRNAs in the nervous system from mammals to flies \[[@CR35], [@CR57]\]. On the other hand, a significant increase of the linear product in foragers was reported previously (XM_393489.3/Am*rsmep2*, *log*~2~ *ratio*∼2.8) \[[@CR18]\]. The observed increase of the circular product circAm*rsmep2* might thus be a consequence of generally increased expression of the host gene, which codes for a RIM-family (Rab3a-interacting molecule) protein. Studies in Tetrapoda species (human, mouse, chicken and so on) show that this family plays an important role in neuronal plasticity, especially in neurotransmitter release and in organizing active zones in plasma membranes \[[@CR58], [@CR59]\].

Contrary to this, the data for circAm*rad* suggests a correlation of acquired task and circAm*rad* levels. Either the task of the bee is influencing circAm*rad* expression or *vice versa*. Its host gene is orthologous to the *radish* gene in *D. melanogaster*, which is known to play a crucial role in the amnesia-resistant memory (ARM). Unlike the long-term memory ARM does not require protein *de novo* synthesis \[[@CR60]\] and thus represents a low costs memory form \[[@CR61], [@CR62]\]. *Rad* also exhibits circRNAs in the fly (Table [1](#Tab1){ref-type="table"}), but whether this circRNA is involved in ARM or whether ARM is also present in honeybees, has not yet been investigated.

Studying circRNAs in the context of synaptic plasticity and neuronal processes promises further insights into the mechanism of task allocation and behavioral regulation of honeybees and probably also of other insects. New evolving techniques such as genome editing using CRISPR/Cas9 which is also available in honeybees \[[@CR63]\] and the microinjection of short interfering RNAs into the *medial ocellus* \[[@CR64]\] will be promising approaches to study the physiological and behavioral effects of altered circRNA levels. The latter could be used to decrease circRNA levels in the brain by specifically targeting circular junctions and thereby promoting their decay. Genome editing, on the other hand, might provide means to induce changes to introns that alter the formation of circRNAs.

Conclusion {#Sec12}
==========

In total, 1263 RNAs were identified in the brain of honeybees (*A. mellifera*) using RNA-Seq. This social model organism complements previous studies in *D. melanogaster* and *B. mori*. Given the evolutionary relationship of the three species, the amount of homologous host genes with circular products, and the number of putative miRNA targets, it can be assumed that circRNAs are a common feature among insects.

Back-splicing in honeybee occurs preferentially towards the end of transcripts and in transcripts with a high number of exons. As reported for *D. melanogaster*, back-splicing is correlated with the length of the 5' and 3' flanking introns \[[@CR35]\]. Additionally, a correlation was found regarding the cofolding probability of these intronic regions as well as their CpG- and cytosine-contents which might be relevant for DNA-methylation. In fact, the methylation was found to be increased for circRNA-flanking introns.

A number of circRNAs identified here were confirmed in independent PCR experiments. Two circRNAs showed a differential expression in nurse and forager bees. While circAm*rsmep2* seems to accumulate with age, the expression of circAm*rad* correlates with the performed task. This finding is the first indication of a link between the circRNA and the social role of honeybees. In general, an increase in phototaxis- and memory-related host genes for circRNAs was found.

Methods {#Sec13}
=======

Collection of bees {#Sec14}
------------------

Bees were derived from colonies with normal age structure and with a naturally mated queen located on the grounds of the University of Würzburg. Bees were considered as nurse bees if they clearly poked their head into open brood cells containing young larvae. Foragers were captured when returning from a foraging flight and having huge pollen loads at their hind legs. Collected bees were frozen in liquid nitrogen immediately. A single cohort colony was established by transferring 2500 newly emerged bees (marked by the same color immediately after hatching) into a small hive together with one queen in one brood frame and one frame with pollen and honey. Single cohort colony bees were collected at the age of eleven days and controlled for their social task. For details regarding biological replicates in the RNA-sequencing analysis and quantitative Real-Time PCR, refer to the respective sections below.

RNA-Seq {#Sec15}
-------

We used a total of four RNA-Seq libraries to determine circular transcripts present in the brain of honeybees. First, an enrichment control was compiled from the brains of ten dissected nurse bees and ten dissected foragers. Total RNA was extracted with Isol-RNA lysis reagent (5PRIME, Hilden, Germany) and treated with DNase I. The sample was divided into two halves. One half (*E*~+~) was treated with 3 units RNase R (epicentre, Madison, USA) per *μ*g total RNA. Digestion was performed for 30 min at 37 °C. For the other half (*E*~−~) an equivalent volume of double distilled water was added. Afterward, both samples were purified using phenol--chloroform extraction. Efficacy of the RNase R treatment was verified in a control experiment shown in Additional file [6](#MOESM6){ref-type="media"}. Second, we took additional samples from ten nurses and ten foragers separately and treated both with RNase R as described above (samples *F*~+~ and *N*~+~, respectively) in order to distinguish task dependent expression levels. In sum, three enriched libraries were independently generated. Only results which were confirmed in at least two libraries (considered as biological replicates) are reported. The RNA quality was ensured using an Agilent 2100 Bioanalyzer. Quality measurements, library preparation and Illumina^Ⓒ^ sequencing (125 nt paired-end) were performed by GATC Biotech AG (Konstanz, Germany). All RNA-Sequencing data was made publicly available via bioproject PRJNA345404, see Table [2](#Tab2){ref-type="table"}. Table 2Summary of RNA-Seq libraries published along with this studySampleSRA IDRoleTreatment\# Reads*E*SRR4343845Both-8,432,479*E* ~+~SRR4343846BothRNase R7,690,777*N* ~+~SRR4343847NurseRNase R5,843,829*F* ~+~SRR4343848ForagerRNase R5,931,097Samples were taken from brains of nurse bees, forager bees or a blend of both

Identification of circular RNAs {#Sec16}
-------------------------------

We used two independent algorithmic approaches for the identification of circular RNAs. In one approach reads were mapped to the NCBI *A. mellifera* genome version 4.5 release 102 (RefSeq GCF_000002195.4) using segemehl (v0.2.0) with the split reads option \[[@CR65]\]. The alignment was subsequently screened for model-free splicing events using the accompanied testrealign tool. In the second approach we used BWA (v0.7.5a) as mapping tool and subsequently screened using CIRI2 (v2.0.6) with default parameters \[[@CR66], [@CR67]\]. Identified junctions were post-processed using custom scripts bundled in our Chiasm suite. Chiasm was also used to perform the statistical calculations later on (e.g. CpG-content, pairing-probability, see below). The full analysis pipeline is publicly available at <https://git.io/chiasm>. More precisely, junctions with almost identical start and end positions were merged if they differed by less than 6 nt. Junctions mapped +/-5 nt next to exon boundaries were corrected to exactly match the boundary. This accounts for small variations in sequencing and mapping, e.g. due to flanking intron sequence being potentially identical to the junctioning exon or indels in the genome. We assigned the respective gene and exon numbers to each hit and normalized the number of JSR to the host gene's total read number. Analogously to present studies in *Drosophila* \[[@CR35]\] we normalized BSJ read counts (*norm*(*n*~∘~)) by dividing the number or circular JSRs (*n*~∘~) by the number of mapped library reads *N* (in millions), divided by reads per kilobases in million reads (RPKM) of the host gene (*g*). The latter is defined as number or reads assigned to the host gene (*n*~*g*~) divided by the length of the gene (*l*~*g*~) in thousand bases and divided by library size of mapped reads, *N*, in millions. $$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document} $$ norm\left(n_{\circ}\right) = \frac{n_{\circ}}{\frac{N}{1,000,000} ~ RPKM_{g}}  $$ \end{document}$$

with $$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document} $$ RPKM_{g} = \frac{n_{g}}{\frac{l_{g}}{1,000} ~ \frac{N}{1,000,000}}  $$ \end{document}$$

We divided the identified circular RNAs into two sets limited by different stringency levels. The low stringency set contains all circRNAs picked up by both approaches (testrealign and CIRI2) with at least three JSRs. In the high stringency set, we only considered BSJs with more than ten JSRs across all libraries as suggested in literature \[[@CR35]\]. Thereby, the BSJ has to be found in library *E*~+~ and at least one other independent RNase R treated library. Moreover, a five-fold enrichment of JSRs in the RNase R treated library (*E*~+~ vs *E*~−~) is required.

Validation of circRNAs {#Sec17}
----------------------

Total RNA was extracted from ten worker bee brains and prepared as described for the RNA-Seq preparation (see above, without enrichment by RNase R). After DNA digestion, 1 *μ*g of RNA was transcribed into cDNA using RevertAid H minus reverse transcriptase (ThermoFisher Scientific) adhering to the manufacturer's specifications. For PCR amplification 15 *μ*mol of divergent primers were added to 10 ng of cDNA with 25 *μ*L of Phusion Polymerase master mix. PCR steps were 30 sec heating to 98 °C followed by 35 cycles of 10 sec denaturation at 98 °C, 10 sec annealing at 62 °C and 8 sec elongation at 72 °C. After a final extension period of 10 min at 72 °C, PCR products were either stored at -20 °C or subjected to agarose gel electrophoresis. Primer sequences are provided in Additional file [7](#MOESM7){ref-type="media"}. The results of PCR verification are provided in Additional file [2](#MOESM2){ref-type="media"}.

Quantification of circRNAs {#Sec18}
--------------------------

750 *μ*L of Isol-RNA lysis reagent (5PRIME, Hilden, Germany) was added to frozen brain samples and homogenized subsequently. After adding 150 *μ*L of chloroform and consequent phase separation the aqueous phase was transferred to 900 *μ*L ethanol (75%). RNA was purified using peqGOLD Total RNA Kit (Peqlab, Erlangen, Germany) following the standard protocol provided by the manufacturer including an optional DNase I digestion step. From each bee 1.5 *μ*g of total brain RNA was transcribed using qScriber cDNA Synthesis Kit (highQu, Kraichtal, Germany). Triplicates of each cDNA (5 *μ*L) were run in a quantitative Real-Time PCR on a Rotor-Gene Q (Qiagen, Hilden, Germany) in a total reaction volume of 25 *μ*L, containing each primer (0.25 *μ*M), TaqMan probe (0.1 *μ*M), Rotor-Gene Multiplex PCR 9Master Mix (Qiagen, Hilden, Germany). The number of biological replicates is specified in Fig. [7](#Fig7){ref-type="fig"}. TaqMan probe sequences are provided in Additional file [7](#MOESM7){ref-type="media"}. The following protocol was used: 60 °C for 1 min, 95 °C for 5 min and 45 cycles at 95 °C for 20 sec and 60 °C for 1 min. Afterward the relative expression to AmEF1 *α* \[[@CR68]\] with the *ΔΔC*~*t*~ method was determined using Rotor-Gene Q software (Qiagen, Chatsworth, CA). Expression of circRNA was compared only if respective groups did not differ in their AmEF1 *α* expression (*p*\>0.05, Student's t-test). For the circRNA candidates, circAm*rsmep2* and circAm*rad* the established TaqMan probe-based assays were designed using outward facing primers. PCR experiments for detection of circRNAs were designed analogously to \[[@CR30]\]. The TaqMan probe binds directly to the circular junction and thus signals can only derive from non-canonical spliced RNAs.

Homology screen and functional annotation {#Sec19}
-----------------------------------------

Predicted circRNAs were correlated to those previously reported for *D. melanogaster* \[[@CR32], [@CR35]\] and *B. mori* \[[@CR34]\]. We matched the loci based on the predicted homologs of the closest protein-coding gene with respect to OrthoDB v9 \[[@CR69]\]. CircRNAs from genes without homolog could thus not be accounted for. Homologous fruit fly genes were then submitted to the online PANTHER annotation platform for further over-representation analysis using Fisher's Exact test with false discovery rate (FDR) multiple testing correction. We included functional annotations with more than five-fold over-representation and FDR below 1%.

Sequence and structural analysis {#Sec20}
--------------------------------

Based on the genomic annotation and the largest spanning transcript of each circRNA that contained exon boundaries, we extracted whether the circRNA contained part of the 5'-UTR, 3'-UTR of a canonical protein-coding transcript or if it exclusively contained coding regions. The number of exons spanned by the transcript was noted for the 5' and 3' end of the BSJ. For comparison to potentially non-circular transcripts, a random control was generated by drawing genes with more than two exons proportionally from all chromosomes and picking exon boundary pairs that were neither from the start nor the end of the transcript. Genes harboring any JSRs found in this study were excluded from this control, see Additional file [3](#MOESM3){ref-type="media"}. A random control of 10,000 such junctions was generated for all following statistical tests. Flanking introns were determined by including the sequence outside of the BSJ until the next exon in the same transcript.

In order to screen for complementarity between flanking intron pairs, the 5' intron was matched to the 3' intron using BLAST \[[@CR47]\] with a word size of six to determine the highest scoring stretch of reverse complementarity. We repeated the procedure with 100 nt from the end of the upstream and 100 nt from the start of the downstream intron, to discern whether approximate regions showed increased complementarity. The same 100 nt portions were used for structural analysis utilizing RNAcofold \[[@CR48]\]. We applied soft constraints to ensure MFE scores solely based on base-pairing between both intronic regions. Both procedures were repeated with all combinations of starts and ends of the respective introns as educated control set (an interaction of the end of the upstream and the end of the downstream intron is probably not relevant). Surprisingly, the results for all combinations were similar. To rule out, that we bias for specific length effects at 100 nt, all calculations were also done with 50 and 200 nt without changing the outcome (data not shown). Introns were checked for GC-content ignoring undetermined residues in the genome sequence (N). Similarly, the mononucleotide frequency of cytosine and the relative frequency of CpG dinucleotides was calculated.

To asses, whether the observed increase of potential DNA-methylation sites is reflected in actual DNA-methylation, we used whole genome bisulfite sequencing data of worker bees that was publicly available. Precisely, we used all native worker libraries provided in BioProject PRJNA104931 \[[@CR52]\] and combined them for this analysis as no differences in average methylation was found between nurse and forager bee libraries for the genes relevant in this study (data not shown). Methylation patterns were analyzed using Bismark \[[@CR70]\] v0.19.1 with Bowtie2 \[[@CR71]\] v2.2.6 for bisulfite specific mapping and default parameters suggested by its authors. For each intron, we counted the average methylation per base on both strands. An average coverage of at least five reads for each intron was required. Calculations were done for 50, 100 and 200 nt as well as for the length of the complete intron where it exceeded 200 nt and numbers were normalized by the respective sequence length. A single-sided Wilcoxon-Mann-Whitney rank-sum test was used to determine the significance of the increase over the control.

miRNA interference analysis {#Sec21}
---------------------------

Predicted and experimentally verified miRNA sequences of *A. mellifera* were obtained from miRBase \[[@CR72]\] release 21. Potential target sites were screened in all exon sequences overlapping with the identified circRNAs using nucleotide two to seven of the mature miRNA sequence, see \[[@CR73]\]. The analysis pipeline is publicly available on git.io, see above. For each potential miRNA binding site, we determined conservation in further *Apis* species (*A. cerana*, *A. dorsata*, *A. florea*) and other eusocial insects (*E. dilemma*, *L. ventralis*, *M. quadrifasciata*, *B. impatiens*, *B. terrestris*) for the seed region with 100 nt up- and downstream using the best BLAST match \[[@CR47]\] in the respective genome. We considered a site conserved if the 6 nt seed region was perfectly conserved among three out of four *Apis* or four out of five eusocial insects, respectively. As a random control we used linear exons, see "Sequence and structural analysis". We split the control to sets of about equal size (42 sets) and applied the above procedure to each set. This results in 42 control datasets where each represents a subset of exons with similar length to avoid a bias due to an over-representation of certain length species. Identified target sites were normalized to sites per 1000 nt.

Additional files
================

 {#Sec22}

Additional file 1Splicing. PDF showing the splice site motif of circRNAs and details on exceptions. (PDF 547 kb)

Additional file 2Verification of circRNAs via PCR. Results of the circularity validation through PCR. (PDF 582 kb)

Additional file 3List of circRNAs in honeybee identified here. Excel table of all circRNAs identified here. The data is presented analogous to Table [1](#Tab1){ref-type="table"} but addressing additional information and details for all RNA-Seq libraries. A second sheet contains a similar list including all 3384 circRNAs identified by both algorithmic approaches based on JSRs in any library. (XLSX 411 kb)

Additional file 4GO term enrichment. Excel table showing the results of a GO term enrichment analysis for the host genes of circRNAs based on homologous fruit fly genes performed with the PANTHER annotation platform. Terms with at least five-fold over-representation and a false discovery rate (FDR) below 1% were considered. From these, we limited the interpretation with a relevant *p*-value threshold of 10^−4^ which is marked in the table. (XLSX 12 kb)

Additional file 5List of potential miRNA targets including conservations. Excel table of all miRNA target sites found on circRNA sequences identified here. It contains detailed data on the target circRNA, the potential position of interaction and its conservation in *Apis*, eusocial insects, *Drosophila* and *Bombyx*. (XLSX 231 kb)

Additional file 6RNase R enrichment control. Experimental control of circRNA enrichment over linear products. (PDF 86 kb)

Additional file 7List of primers and probes for PCR and TaqMan assay. Excel table with a list of all PCR primers used in this study. A second sheet lists the TaqMan probes. (XLSX 10 kb)

BSJ

:   Back-spliced junction

circRNA

:   Circular RNA

cDNA

:   Complementary DNA

C ~*t*~

:   Threshold cycle

FDR

:   False discovery rate

GO

:   Gene ontology

ID

:   Identifier

JSR

:   Junction-spanning read

MFE

:   Minimum free energy

miRNA

:   Micro RNA

PCR

:   Polymerase chain reaction

RNA-Seq

:   RNA-sequencing

SCC

:   Single cohort colony

SRA

:   Sequence read archive

UTR

:   Untranslated region

We thank Marietta Thüring for fruitful discussions.

Funding {#d29e2684}
=======

This work was partially funded by the PostDoc Plus Grant to MT of the Graduate School of Life Sciences, University of Würzburg.

Availability of data and materials {#d29e2689}
==================================

All RNA-Seq datasets used in this study are available via NCBI BioProject PRJNA345404 (<https://www.ncbi.nlm.nih.gov/bioproject/PRJNA345404>). For details on individual SRA IDs see Table [2](#Tab2){ref-type="table"}. The analysis pipeline is publicly available at <https://git.io/chiasm>.

ML and MT conceived the study. ML and CT advised on the sequencing procedure. CE grew bees and conducted the RNA-Seq preparation and experiments. CT carried out the bioinformatic analysis. CT and CE performed PCR verifications of circular transcripts. MT carried out the quantification assays. ML, MT and CT wrote the manuscript. All authors read and approved the final manuscript.

Consent for publication
=======================

Not applicable.

Competing interests
===================

The authors declare that they have no competing interests.

Publisher's Note
================

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
