RNA editing is abundant and correlates with task performance in a social bumblebee

Colonies of the bumblebee Bombus terrestris are characterized by wide phenotypic variability among genetically similar full-sister workers, suggesting a major role for epigenetic processes. Here, we report a high level of ADAR-mediated RNA editing in the bumblebee, despite the lack of an ADAR1-homolog. We identify 1.15 million unique genomic sites, and 164 recoding sites residing in 100 protein coding genes, including ion channels, transporters, and receptors predicted to affect brain function and behavior. Some edited sites are similarly edited in other insects, cephalopods and even mammals. The global editing level of protein coding and non-coding transcripts weakly correlates with task performance (brood care vs. foraging), but not affected by dominance rank or juvenile hormone known to influence physiology and behavior. Taken together, our findings show that brain editing levels are high in naturally behaving bees, and may be regulated by relatively short-term effects associated with brood care or foraging activities.


Primers used for qPCR amplification of ADAR transcripts in B. terrestris and A. mellifera
The primer sequences and details are summarized in Supplementary Table 1.

Selection of reference genes for qPCR analyses in different tissues of B. terrestris
We tested several control genes for experiments comparing mRNA levels among different B. terrestris tissues. We compared genes that were used as control genes in previous qPCR studies with bumblebees 1 that show little variation in an RNAseq study comparing brain expression in bees differing in task performance, JH levels, or dominance (Shpigler et al, unpublished). The candidate genes compared were: Ribosomal protein L13 variant 1 (RPL13), Phospholipase A2 (PLA2) ,Arginine kinase (Argk), 40S Ribosomal protein S16-like (S16), Elongation factor Tu (Ef-Tu), and Elongation factor 1a (EF1a). Primers were tested on a standard curve using a dilution series of cDNA. The Correlation coefficients (R 2 ) and the amplification efficacy (%E) were determined based on the standard curve. The amplification efficacy was determined based on the formula: E(%) = 100 * (1 − 10 − 1 Slope ). Candidate genes that showed amplification efficacy values below 95% or above 105% were excluded from further analyses. Second, we determined mRNA levels in various tissues of 1-day-old bees. Candidate reference genes showing an expression difference of more than 2 cycles among body parts were excluded from further analysis. Third, we tested the expression variance using 3 designated software programs: geNorm 2 , NormFinder 3 and BestKeeper 4 . Each program received the qRT-PCR results and calculated the variance between and within tissues based on different algorithms. All three programs identified S16 as the most stable (Supplementary Table 2). RPL13 and S16 were selected based on the lowest tvariance across the three programs. The primvers used for RPL13 are (the asterisk mark exon-exon boundary): F -GTTGCACCAAG*ACCTGTGAAG; R -

Experiment 3. The influence of juvenile hormone on brain RNA editing
We manipulated JH levels by allatectomy and treatment with JH-III. We used newly emerged (< 18 hrs post pupal emergence) worker bees. At this age the cuticle of adult bumble bees is relatively soft and easy to manipulate. The collected bees had free access to sugar syrup and pollen ad libitum both before and after the allatectomy operation. For the allatectomy treatment ('CA-'), we anesthetized the bees on ice for 5-20 min (the variation in chilling duration was due to individual differences in body size, and consequence of the dissections order as the bees were chilled in groups of four). The anesthetized bees were fixed under a stereoscopic microscope (Nikon SMZ645, X50) to an ice-chilled metal stage. The bees were fixed with molded modeling clay such that their dorsal side faces up, and the head bent down to expose the thin neck cuticle connecting the thorax and the head. We used fine scalpel to open a latitudinal incision in the posterior part of the head capsule, and moved the inner membrane and trachea to expose the corpora allata (CA) glands. Both corpus allatum were gently grasped with fine forceps and detached. The entire procedure took between 2-5 minutes. The cuticle resumed to its original shape and the incision appeared self-sealed within few hours after the operation. Sham-operated bees ('Sham') were handled and dissected in a similar way to the CA-bees, but their CA glands were only touched gently and not detached. Control bees ('Control') were anesthetized and handled similarly, but were not operated. At the end of the operation the bees were placed in a small wooden cage (12x5x8 cm), with other similarly manipulated bees, and were left to recover overnight in an incubator (32ºC ± 1ºC, 70%±5% RH). On the second day, the surviving bees from each treatment group were assigned to groups of three, each transferred to a fresh wooden cage (12x5x8 cm).
The bees in each group were colored marked for individual recognition and were kept in the incubator (32ºC ± 1ºC, 70%±5% RH) for four days with ad-libitum food supply, 70% sugar syrup and fresh pollen cake. For Most of the mortality occurred during the first day, before assigning the bees to groups. The average survival rate for the first day in the three experiments detailed below was 50% for the allatectomized ('CA-') bees, 80% for the sham, and 100% for the control bees. Survival of the bees after placed in groups was checked daily on days two to five. The recorded survival during these days was: 86% (45/52) for the CA-and CA-+JH, 94% (34/36) for the sham, and 100% (36/36) for the control bees. Only groups in which all three bees survived for the whole experiment were used for RNAseq analysis in order to keep the social state similar across groups. The bees were flesh frozen in liquid nitrogen on day five, and immediately transferred into marked 2 ml tubes that were immersed in dry ice. The tubes were later transferred to a deep freezer (-80ºC) until further analysis. This procedure assured that the bee samples are deeply frozen from the time of collection until the time of RNA extraction.

Assessing genetic variation within a colony
To estimate the genetic variability in our samples, we applied the GATK pipeline to the DNA-seq data of four worker bees that were used for the RNA-seq analyses. Average coverage depths in these four DNA-seq samples were 32, 32, 39 and 34X. More than 96% of the genome (236/245 Mbps) was covered by at least one sample. The average number of single nucleotide variants (SNVs) with respect to the reference genome was 1,286,561 (1,268,774-1,322,571), where most of these, 717,255 (56%), were common to all four bees, likely representing genetic variation between the Israeli B. terrestris population used in our study and the individual used to generate the official B. terrestris genome draft. Only 3-15% (36,223-193,007) of the variants were unique to a single bee. These results agree with our expectation that the four bee genomes are much closer to each other than to the reference genome. In total, 1,906,220 single nucleotide variants The four conserved editing sites within B. terrestris and D. melanogaster are denoted with arrows. Amino-acid replacements resulting from editing, are designated above the arrows, and the amino-acids that are prone to editing are denoted with squares. The red arrows point to editing sites that are also conserved in cephalopods and the black arrow points to the editing site that is also conserved in mammals. colonies. Values are presented as mean ± SE. Sample size is shown inside each column.