Figure 1. Presence of Ingested dsRNA in the Honey Bee Hemolymph and Secretion into Worker and Royal Jellies
(A) Presence of ingested dsRNA in the hemolymph. Probe-free northern blot analysis performed on pooled raw hemolymph extracts (10 μl per well). Raw hemolymph extracts were collected from bees that were fed on 50% sucrose solution (w/w) containing DIG-labeled dsRNA-GFP or sucrose solution only. The 430-bp band represents free full-length dsRNA. Multiple lanes with a single label represent independent biological replicates.
(B) Association of ingested dsRNA with a protein complex in the hemolymph. Probe-free northern blot analysis performed on untreated or protease K-treated pooled raw hemolymph extracts (10 μl per well). The hemolymph in both treatments was derived from the same hemolymph sample.
(C) Occurrence of ingested dsRNA in the larval diets and newly laid eggs from dsRNA-GFP treated and untreated mini hive. Northern blot analysis of 1 μg total RNA extracted from worker jelly (WJ), royal jelly (RJ), and eggs. Eggs were laid in untreated mini hives by queens that were transferred from dsRNA-treated or untreated mini hives. Purified dsRNA-GFP was used as a positive control and a size marker for full-length dsRNA.
Interestingly, two additional dsRNA∗ bands of higher molecular weight (corresponding approximately to 2.5 and 4 Kbp) were also observed in raw hemolymph extracts (Figure 1A). Two scenarios can be drawn to explain the appearance of high molecular weight labeled RNAs: (1) recombination between dsRNA∗ and other RNAs or (2) association of the dsRNA∗ with other components of the hemolymph. Following Proteinase K digestion of the hemolymph extract, the size of the higher bands shifted back to the expected 430 bp, indicating that the higher compounds were complexes of circa full-length dsRNA and hemolymph protein(s) (Figure 1B). Consistent with a previous report (Garbian et al., 2012), no processed dsRNA forms could be detected in the hemolymph extracts tested.
Presence of Ingested dsRNA in the Worker and Royal Jelly
The finding of ingested dsRNA-protein complexes in the hemolymph indicates a possible active mechanism for uptake and translocation of dsRNA through the honey bee’s circulatory system and possible spread to the jelly-producing glands. If this happens, mobile dsRNA may be present in the jelly. Therefore, we next tested whether ingested dsRNA is secreted in the diet of worker- and queen-destined larvae.
Mini hives with circa 250 worker bees and reproductive queens were established and fed on sucrose only (control) or sucrose solution mixed with dsRNA carrying a foreign GFP sequence (dsRNA-GFP). Worker jelly was collected from brood cells containing 5th instar worker larvae and royal jelly was harvested from queen brood cells with 3rd–4th instar larvae (see STAR Methods). We could detect the presence of dsRNA-GFP by northern blot analysis performed on total RNA extracted from worker and royal jellies. Notably, while the full-length dsRNA could be detected, additional degraded or processed GFP-RNA forms appeared in both jellies (Figure 1C).
RNA Is Horizontally Transferred among Bee Generations
The presence of ingested dsRNA in the circulatory system indicates systemic spread and the possibility of dsRNA transport to the food glands (Figures 1A and 1B). This was further supported by dsRNA presence in both worker and royal jellies (Figure 1C). Hence, dsRNA could presumably be transmitted from nurse bees to the young larvae through jelly consumption. In order to test whether RNA can be transferred horizontally among bee generations down the line, we established reproductive mini hives that were fed on sucrose solution (control hives) or sucrose solution containing dsRNA-GFP. During the experiment, adult workers, larvae, pupae, and newly emerged bees were collected and analyzed for the presence of dsRNA-GFP (Figure 2A). An RNA slot-blot assay with a GFP-specific probe was conducted and showed, as expected, the presence of dsRNA-GFP in adult workers that had been consuming directly dsRNA in their diet. It also demonstrated the occurrence of dsRNA-GFP in larvae that had been consuming jelly secreted by treated bees (Figure 2B). As previously reported, dsRNA is stable and persists in treated adult bees for a few days post feeding (Maori et al., 2009). Here, we show that dsRNA, which is transmitted from nurse bees to the larvae, persists in subsequent developmental stages including pupae and newly emerging bees. It appears that ingested foreign dsRNA diminished with time but could be detected at least 14 days after the last dsRNA feed (Figure 2B).
Figure 2. Biologically Active dsRNA Is Horizontally Transferred among Bee Generations
(A) Experimental design for detecting RNA transmission within the hive. Reproductive mini hives contained circa 250 adult bees and brood at different developmental stages and an active queen. Treated colonies were provided with 300 μg dsRNA-GFP in 50% sucrose solution (w/w) per feeding. Control mini hives were fed on 50% sucrose solution (w/w) only. Adult bees, 2nd–3rd instar larvae and pupae were sampled from the available developmental stages at each time point.
(B) Occurrence of dsRNA in adult bees and its transfer to the next generation. RNA slot-blot analysis of 1.5 μg total RNA extracted from individual larvae, pupae, and adult bees from dsRNA-GFP treated or untreated colonies (Figure 2A).
(C) Experimental design to test horizontal RNA transfer among bees. Reproductive mini hives contained circa 250 bees and an active queen. Treated colonies were provided with 200 μg dsRNA-GFP in 50% sucrose solution (w/w) per feeding. Control mini hives were fed on 50% sucrose solution (w/w) only.
(D) Transfer of dsRNA from worker bees to nourished larvae. 1st instar larvae were transferred from untreated hive and nourished by workers from dsRNA treated or untreated mini hives for four days. Northern blot of 5 μg total RNA extracted from individual 5th instar larvae.
(E) Experimental design to test vertical RNA transfer among bees. Reproductive mini hives contained circa 250 bees and an active queen. Treated colonies were provided with 200 μg dsRNA-GFP in 50% sucrose solution (w/w) per feeding. Control mini hives were fed on 50% sucrose solution (w/w) only. Queens from dsRNA treated or untreated mini hives were removed to untreated mini hives. The queens were then isolated for two days to acclimatize. Next, the queens were released and their newly laid eggs were collected.
(F) Transmissible RNA is biologically active. Vitellogenin (Vg) knockdown in ten-day-old workers that were nourished as larvae by dsRNA-Vg treated bees. Vg expression was quantified by RT-qPCR. Individual workers were tested in the untreated (N = 6), dsRNA-GFP (N = 6), and dsRNA-Vg (N = 7) treatments. Left y axis corresponds to ΔCt values; individual bee values are shown as dark circles; mean value as dark rectangle; and error bars representing standard error of the mean. Right y axis corresponds to Vg-RNA quantification relative to the untreated, and mean values are shown as colored bars. Different letters above the bars indicate statistically significant difference between treatments according to the Tukey-Kramer (HSD) test (p < 0.05). Higher ΔCt values indicate lower RNA quantification.
While the presence of dsRNA-GFP in the jellies supports horizontal transfer to the progeny, it could also be explained by vertical dsRNA transmission from queen to egg. To distinguish between the two possible transmission routes, we designed an experiment in which only horizontal nurse-bee-to-larvae transmission could occur (Figure 2C). A similar mini hive set up was established, and sucrose solution containing dsRNA-GFP or sucrose solution only (control hives) was provided for five days. On day six, combs containing 1st instar larvae from control hives were transferred either to a dsRNA-treated hive or to another control hive. The transferred larvae were then allowed to develop four additional days in the new untreated or dsRNA-treated host colonies. On day ten, the 5th instar larvae were collected from the transferred combs, washed rigorously, and analyzed for the presence of dsRNA-GFP. Northern blot analysis performed on total RNA detected a GFP-RNA signal in larvae that were nourished by dsRNA-treated bees. Consistent with the dsRNA-GFP signal in jellies (Figure 1C), degraded or processed GFP-RNA forms are detected in recipient larvae (Figure 2D). All together, these data demonstrate an environmentally mediated horizontal RNA transfer route among honey bee generations.
RNA transfer from nurse bee to larvae does not rule out maternal queen deposition of RNA to eggs and its persistence throughout the progeny development. Therefore, we attempted exploring whether vertical dsRNA transmission also occurs in honey bees (Figure 2E). To this end, the reproductive mini hive system was applied and dsRNA-GFP-containing sucrose solution or sucrose-only solution (control hive) was provided for six days. On day six, simultaneous queen swaps among treatments were performed as follows: (1) the queen from dsRNA-GFP hive replaced the queen from a control hive, and (2) the substituted control queen replaced a queen from a different control hive. On day eight, after two days of acclimatization, the newly introduced queens were released and allowed to lay eggs in their new colony. On day ten, eggs were collected, pooled, and analyzed for the presence of dsRNA-GFP. We could not detect GFP-RNA signal in eggs laid by a queen that was previously nourished with royal jelly provided by dsRNA-treated bees (Figure 1C). We acknowledge that such negative detection could be explained by the sensitivity limitation of the northern blot assay. Thus, we conclude that, while vertical transmission from queen to larvae might occur to some extent, horizontal nurse-to-larvae transfer is the main route that RNA flows among generations and members of the hive.
Transmissible RNA is Biologically Active in Recipient Bees
We next asked whether horizontally transferred RNA (i.e., transmissible RNA) is biologically active within recipient individuals. It has been previously demonstrated that supplementing dsRNA into the natural larval diet induces potent RNAi against endogenous RNA (Guo et al., 2013, Nunes and Simões, 2009). Moreover, similar application with dsRNA, corresponding to deformed wing virus (DWV) and sacbrood virus (SBV), effectively reduced viral RNA titer as well as disease symptoms in infected worker larvae and adults (Desai et al., 2012, Liu et al., 2010). Notably, these studies showed that ingestion of jelly-containing dsRNA results in sustainable gene silencing that lasts until adulthood (Desai et al., 2012, Guo et al., 2013, Nunes and Simões, 2009).
Unlike previous reports, we examined whether dsRNA, which is originated from nurse bees and secreted into the jelly, could elicit RNAi response in the progeny. To answer this question, we established a minimal hive system in plastic boxes containing circa 150 workers and a comb with eggs and young larvae. The adult bees were fed for eight days on sucrose solution (untreated), sucrose solution mixed with dsRNA-GFP (non-specific dsRNA control), or dsRNA that matched the vitellogenin mRNA sequence (dsRNA-Vg). When the brood cells were sealed, the adult bees were removed and the combs were kept until new workers emerged. We then collected ten-day-old workers and analyzed the expression levels of vitellogenin by RT-qPCR. Consistent with persistence of jelly-transmitted dsRNA (Figure 2B) and in agreement with the aforementioned reports (Desai et al., 2012, Guo et al., 2013, Nunes and Simões, 2009), we observed vitellogenin knockdown in adult workers that were nourished as larvae by dsRNA-fed nurse bees (ANOVA: F 2,16 = 5.91, p = 0.012; Figure 2F). Five out of seven bees from the dsRNA-Vg treatment had the highest delta CT (ΔCT) values, higher than any of the 12 bees in the controls (p = 0.023 by the conservative non-parametric Wilcoxon test). An ANOVA comparing the dsRNA-Vg treatment with the two controls combined yields F1,17 = 12.47, p = 0.0026. Therefore, we concluded that transmissible RNA, at least in a dsRNA form, is biologically active in recipient bees.
Naturally Occurring RNA in Worker and Royal Jellies
The previous experiments described a mechanism that enables, through jelly secretion and consumption, transmission of biologically active RNA among individuals within and between generations in the hive. These findings suggest a naturally occurring transmissible RNA in honey bees. In line with this hypothesis, it has been reported that both worker and royal jellies contain small honey bee RNA populations, demonstrating endogenous bee RNA secretion into the larval diet (Guo et al., 2013, Zhu et al., 2017). To further explore the RNA repertoire in the jellies, including the natural occurrence of exogenous and pathogen-related RNA, we adapted a small RNA-seq protocol to sequence full-length RNA up to 200 nt (see STAR Methods). Samples of royal jelly were collected from 3rd instar queen larvae brood cells, and worker jelly was collected from 5th instar worker larvae brood cells. The jellies were harvested from untreated healthy-looking hives.
Size distribution analysis of sequenced full-length RNA indicated that worker and royal jellies have different profiles, with RNAs corresponding to 39 and 72 nt mainly differentiating between the two jellies (Figures 3A and S1A). We next applied a metagenomics analysis to identify the origin of jelly RNA and, again, found different profiles in the worker and royal jellies (Figures 3B and S1B; Table S1). Surprisingly, bee RNA represents only a minor fraction in both jellies, representing 0.58% and 3.55% of worker and royal jelly, respectively. Instead, large proportions of plant, fungi, and bacteria were identified alongside sequences originating from unknown sources. Remarkably, RNA fragments corresponding to various exogenous bee-affecting viruses could also be detected in both jellies.
Figure 3. Differential Naturally Occurring RNA Populations in Worker and Royal Jellies
(A) Size distribution of naturally occurring royal and worker jelly RNA. RNA size was determined through sequencing full-length jelly RNA by a small RNA-seq protocol that was adapted to sequence broad full RNA length spectrum (i.e., 15–200 nt). Data represent a merged analysis of three biological repeats per jelly and are presented as the normalized number of reads per million (RPM). Common peak sizes are marked in black font, and differential sizes in red font.
(B) RNA-based metagenomic analysis of royal and worker jellies. Data represent a merge analysis of three biological repeats per jelly.
(C) Classification of honey bee (Apis mellifera) RNA in royal and worker jellies. RNA types classification derived from annotated sequences in worker and royal jellies (54.02% and 86.96%, respectively). Non-coding RNA classification derived from the RNA types analysis of worker and royal jellies (3.14% and 0.19%, respectively).
(D) Occurrence of honey bee transposable elements RNA in royal and worker jellies. TE classification derived from RNA types analysis of worker and royal jellies (0.2% and 0.01%, respectively)
(E) Detection and proportion of putative endogenous (Apis mellifera) dsRNA in worker and royal jelly samples. DsRNA is detected when two distinct RNA molecules had at least 25-nt base pairs, and the overhang on either side did not exceed 100 nt.
(F) Size distribution of putative worker and royal jelly dsRNAs that are mapped to the Apis mellifera genome.
(G) Worker and royal jelly contain diverse putative endogenous dsRNAs derived from un-annotated, protein-coding and non-coding honey bee genes. Data in (A), (B), and (C) represent a merge analysis of three biological repeats per jelly.
See also Figure S1 and Tables S1, S2 and S3.
We next further characterized jelly RNA that corresponds to the honey bee genome (Figures 3C and S1C–S1E; Table S2). Like the previous jelly RNA analyses, different honey bee RNA profiles were detected in worker and royal jellies. In both jellies, the large proportions represent bee RNA derived from protein coding genes followed by tRNAs. However, worker jelly is relatively enriched in ribosomal, transposable elements (TE) and non-coding RNA. Interestingly, differential TE RNA occurrence could be detected among the jellies, which is mainly associated with LTR-retrotransposons and TIR transposons (Figure 3D).
Previous experiments and other reports demonstrated that ingestion of jelly-containing dsRNA downregulates specific gene expression in recipient larva and adult honey bees. Therefore, we next screened for duplexed honey bee RNA in worker and royal jellies. We found substantial proportions of putative endogenous dsRNA, ranging from 12% to 37% of total honey bee jelly RNA (Figure 3E). The putative dsRNA fragments have a broad size distribution and are mainly derived from protein-coding genes but also from tRNAs, rRNA, TEs, and un-annotated DNA (Figures 3F and 3G; Table S3). Notably, royal jelly contained higher proportion of bi-directional RNA fragments derived from tRNA genes.
We hypothesized that bees treated with an IAPV-specific dsRNA in previous field trials (Hunter et al., 2010) may have transmitted the antiviral RNA to other bees, resulting in transmissible protection against the viral infection. Therefore, we also characterized the viral RNA in worker and royal jellies.
Overall, RNA corresponding to ten and four bee-affecting viruses could be detected in royal and worker jelly, respectively (Figures 4A and 4B; Table S4). The most abundant viruses in both jellies were DWV and Varroa destructor virus 1 (VDV-1). Both sense and antisense viral RNA strands are detected for most viruses. The presence of replicative forms (anti-sense viral genome) suggests an intracellular origin of the viral RNA and its secretion into the jelly rather than RNA derived exclusively from environmental capsids. Additionally, we analyzed the size distribution of viral sequences and identified diverse sense and anti-sense viral RNA fragments in both jellies (Figures 4C and 4D). Interestingly, while both worker and royal jellies contain large populations of small RNAs (Figures 3A and S1A), almost no small viral RNAs (20–25 nt) were identified (Figures 4C and 4D). To assess fragment diversity as well as potential occurrence of base-paired viral RNA, reads were mapped against corresponding viral genomes (Figures 4E, 4F, and S2). Multiple mappings were observed in most viruses, especially in the abundant VDV-1 and DWV. Remarkably, presence of long (>25 nt) overlapping viral sense and antisense RNA fragments is somewhat common, suggesting naturally occurring viral dsRNA in worker and royal jellies (Figures 4E, 4F, and S2).
Figure 4. Diverse Viral RNA Fragments Naturally Occur in Royal and Worker Jellies
(A) Occurrence of diverse viral RNAs in royal jelly.
(B) Occurrence of diverse viral RNAs in worker jelly.
(C) The viral RNA in royal jelly is fragmented.
(D) The viral RNA in worker jelly is fragmented.
(E) Genome distribution of VDV-1 RNA fragments from royal jelly.
(F) Genome distribution of VDV-1 RNA fragments from worker jelly.
Three biological samples were individually sequenced per jelly; sequencing outcomes were merged and analyzed. Data are presented as the number of reads normalized to log transcripts per kilobase million (TPM).
See also Figure S2 and Table S4.
Discussion
Employing dsRNA as a model, this study reports on an environmentally mediated RNA cycle among honey bees. The cycle is engaged by consumption of RNA-containing diet by an individual bee. Then, the ingested RNA is spread from the digestive system through the gut cells to the hemolymph, where it is associated with a protein complex. A systemic RNA signal reaches the food secretion glands of nurse bees and is transmitted to the progeny, again, through RNA-containing jelly consumption (Figure 5). This phenomenon is driven by horizontal RNA transfer among individual bees and across generations. Hence, it demonstrates an inherent non-organism autonomous RNA—a transmissible RNA route in honey bees. Such a route could involve transmission of diverse exogenous and endogenous RNA types, including double- and single-stranded RNA corresponding to protein-coding and non-coding genes.
Figure 5. Working Model for Transmissible RNA Pathway in Honey Bees
Bees are able to take up RNA from the environment through ingestion. The ingested RNA is taken up from the digestive system through the gut cells to the circulatory system, the hemolymph. In the hemolymph, ingested extracellular RNA is associated with a protein complex and systemically spread, including to the jelly-producing hypopharyngeal and mandibular glands. Then, ingested and other RNAs are secreted in the royal and worker jellies. A new environmental RNA cycle is initiated through ingestion of jelly-containing RNA. A few potential RNA sources could participate in the transmissible RNA pathway including systemic antiviral RNAi and endogenous mobile RNA, jelly-secreting gland transcription as well as environmental hive RNA (e.g., plant, fungi, and bacteria).
Larva and adult honey bees can ingest biologically active dsRNA (Maori et al., 2009, Nunes and Simões, 2009). However, the ability of bees to efficiently take up dietary small RNAs (e.g., miRNAs and siRNAs) is currently debatable (Chen et al., 2014, Guo et al., 2013, Masood et al., 2016, Zhu et al., 2017). While exchange of dietary RNA among adult bees can reasonably occur via trophallaxis (LeBoeuf et al., 2016), in our experiments, secretion of jelly that contained dsRNA-GFP required systemic distribution of the environmentally consumed dsRNA within the nurse bee. A few RNA uptake mechanisms have been reported including internalization by extracellular vesicles, dsRNA channels, and receptor-mediated endocytosis (Feinberg and Hunter, 2003, McEwan et al., 2012, Saleh et al., 2006, Valadi et al., 2007). Although these mechanisms have demonstrated cellular RNA import, little is known about the systemic spread of RNA through the insect’s circulatory system. Here, we show that ingested-dsRNA is exported into the bee’s hemolymph where at least part of it is not naked, but associated with proteins, forming extracellular ribonucleoprotein complexes. In agreement with previous studies, which determined import preference of long dsRNA (Bolognesi et al., 2012, McEwan et al., 2012, Saleh et al., 2006) in C. elegans and Drosophila, we found that the hemolymph dsRNA-protein complex is comprised of long and mostly non-processed dsRNA (Figure 1B). These findings led us to propose potential parallel complementary functions of the extracellular dsRNA-protein complex involving stabilization, translocation, and introduction of the circulatory RNAi signal to recipient cells and tissues in a specific and/or non-specific manner. Yet, investigating such potential roles would require the identification and characterization of the RNA-binding hemolymph proteins in future studies. Northern blot analysis of dsRNA-GFP in worker and royal jellies indicated that RNA secreted by the feeding gland cells had undergone partial processing or degradation (Figure 1C), therefore, suggesting that once systemic hemolymph dsRNA is taken up by cells, it is engaged with RNA processing (e.g., Dicer) or degrading factors.
Suspected prolonged viral disease resistance in field hives fed on dsRNA homologous to a bee virus (IAPV) (Hunter et al., 2010) suggested long-term effect of RNAi in treated colonies 3–4 months after the last dsRNA treatment. RNAi maintenance via dsRNA amplification, driven by the viral RdRp and/or endogenous expression (Goic et al., 2013, Maori et al., 2007, Tassetto et al., 2017), can explain silencing persistence in an individual bee. Nevertheless, it is unlikely to explain long-term protection at the colony level since the honey bee’s lifespan during the summer is circa six weeks. Worker larvae are fed exclusively jelly for three days and then predominantly jelly, honey, and pollen mix. Therefore, the presence of dsRNA-GFP in jelly demonstrated jelly-secretion mediated RNA transfer to next generations (Figure 1C). Our experiments show that horizontal RNA transfer is the main route to share and spread RNA, at least in a dsRNA form, among the bee population in the hive. However, the data cannot rule out additional vertical transmission of RNA from queen to eggs. Ingestion of dsRNA-supplemented jelly under natural or in vitro conditions has been reported to confer efficient endogenous and exogenous gene silencing in honey bee larvae and newly emerged adults (Desai et al., 2012, Guo et al., 2013, Nunes and Simões, 2009, Zhu et al., 2017). Here, we show that dsRNA that is secreted into the jelly and consumed by larvae is also biologically potent and can induce a long-lasting RNAi that persists until adulthood (Figure 2F). Therefore, an interpretation of our results leads us to conclude that RNA transfer to larvae could potentially prime anti-viral RNAi and explain the suspected long term protection against viral disease in infected hives fed on IAPV-dsRNA (Hunter et al., 2010).
RNAi has been established in insects, including honey bees, as a key immune response against viruses (Gammon and Mello, 2015, Maori et al., 2009, Tassetto et al., 2017). During infection, local RNAi develops into a systemic signal to control viral spread and propagation in distant cells and tissues. Our data indicate that systemic RNAi signal is not limited to the infected bee, but spreads beyond to other individuals in the hive. Diverse fragments of bacterial, fungal, and viral RNA naturally occur in both jellies, representing 25% and 16.75% of total worker and royal jelly RNA, respectively (Figures 3B, 4A, and 4B). The presence of both sense and anti-sense viral sequences in the jellies suggests secretion of viral RNA originated from cells. Additionally, potential Dicer substrates seem to be common among natural viral jelly RNA (Figures 4E, 4F, and S2), supporting transfer of naturally occurring RNAi triggers to the larvae. We previously demonstrated cross species bi-directional RNAi transfer between the honey bee and Varroa mite, which has been applied for Varroa control (Garbian et al., 2012, Sela et al., 2015). Thus, in addition to individual defense, transmissible RNA could elicit a colony-level protective outcome. Relative to other insects, the honey bee’s genome encodes a reduced number and variety of immune gene families. It has been therefore suggested that the bees’ behavioral group defense provides a complementary level of immunity, compensating for the reduction of immune genes (Evans et al., 2006). We hypothesize another form of collective defense: social immunity that is engendered through transfer of pathogen-related RNA among members in the hive. Natural RNA-driven social immunity requires further evidence in future research to elucidate the physiological and biochemical responses of recipient larvae to ingested jelly RNA, including engagement with RNAi factors and processing into small RNAs.
It is generally agreed that RNAi evolved as a defense mechanism against selfish nucleic acids and further diversified to regulate endogenous gene expression. The presence of differential naturally occurring RNA among worker and royal jellies points toward a potential effect of transmissible RNA on genome function in recipient bees. Indeed, supplementing jelly with endogenous or exogenous miRNAs that are naturally enriched in worker jelly affected gene expression as well as developmental and morphological characters of newly emerged workers and queens (Guo et al., 2013, Zhu et al., 2017). Here, we showed that worker and royal jellies contain substantial proportions of putative honey bee dsRNA (Figures 3E–G), presumably derived from bi-directional transcription (Pelechano and Steinmetz, 2013). We speculate that bee-to-larva RNA transfer could also play a role in epigenetic dynamics among honey bees. A general involvement of epigenetics in the phenotypic plasticity of female bees has been demonstrated (Kucharski et al., 2008, Lyko et al., 2010, Spannhoff et al., 2011). Nonetheless, it is still not understood how caste-specific epigenetic marking is directed. Continuous uptake of regulatory jelly RNA could contribute toward a specific gene-expression profile of a genome with a potential to differentiate into two castes; hence, shifting toward an emergence of worker or queen. Making of a queen is a multilevel process and presumably there are numerous factors interplaying one with the other. Further research should molecularly determine the impact of jelly RNA and its relation to other identified players, such as Royalactin (Kamakura, 2011).
The presented mechanistic experiments with artificial RNA and the occurrence of natural RNA populations in the jellies indicated not only an RNA share among bees but also its ability to persist in an external non-sterile environment. It has been reported that royal jelly contains bactericide components (Romanelli et al., 2011). Although inhibition of microbial growth could contribute to the stability of RNAs in the jelly, it cannot solely protect against environmentally distributed nucleases and physical degradation. Moreover, after ingestion, jelly RNAs need to be further stabilized in the digestive system of the larva and adult queen, which harbors a diverse microbial population (Kwong and Moran, 2016). This raises the question of how does the jelly support environmental persistence and activity of RNA (Maori et al., 2019)?
STAR★MethodsKey Resources Table
REAGENT or RESOURCESOURCEIDENTIFIERBiological SamplesChemicals, Peptides, and Recombinant ProteinsCritical Commercial AssaysDeposited DataExperimental Models: Organisms/StrainsOligonucleotidesSoftware and Algorithms
Royal jelly | B. Triwaks Bee Research Center, Hebrew University of Jerusalem; Springwell Apiaries, UK | n/a |
Worker jelly | B. Triwaks Bee Research Center, Hebrew University of Jerusalem; Springwell Apiaries, UK | n/a |
Bee hemolymph | B. Triwaks Bee Research Center, Hebrew University of Jerusalem | n/a |
Polen supplement patties | ManLake | FD-374 |
DIG RNA Labeling Kit (SP6/T7) | Roche Life Science | 11175025910 |
PCR DIG Probe Synthesis Kit | Roche Life Science | 11636090910 |
Cap-Clip Tobacco Acid Pyrophosphatase | CellScript | C-CC15011H |
5x HOT FIREPol® EvaGreen®qPCR Mix Plus | Solis Biodyne | 08-24-00001 |
RNA 6000 Pico Chip and reagents | Agilent | 5067-1513 |
NEXTflex Small RNA-Seq Kit v3 | Bioo Scientific | NOVA-5132-05 |
RNA-seq | This paper; ArrayExpress data | https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-6557/ |
Raw data | This paper; Mendeley data | https://data.mendeley.com/datasets/wv4s233nx7/draft?a=21d6454c-ead9-4a13-bcd8-f6395652424c |
Honey bee (Apis mellifera) | B. Triwaks Bee Research Center, Hebrew University of Jerusalem | n/a |
See Table S5 of this paper | n/a | |
In-house scripts | Github | https://doi.org/10.5281/zenodo.1302437 |
JMP Statistical Software, Version 13 | SAS Institute | https://www.jmp.com/en_gb/home.html |
TREP Database, Release 16 | Wicker et al., 2007 | http://botserv2.uzh.ch/kelldata/trep-db/ |
RepBase, Version 22.04 | Jurka, 1998 | https://www.girinst.org/repbase/ |
gProfiler | Reimand et al., 2016 | https://biit.cs.ut.ee/gprofiler_archive2/r1760_e93_eg40/web/ |
pairedBamToBed12 | Quinlan and Hall, 2010 | https://github.com/Population-Transcriptomics/pairedBamToBed12 |
intersectBed, bedtools version 2.17.0 | Quinlan and Hall, 2010 | https://github.com/Population-Transcriptomics/pairedBamToBed12 |
Contact for Reagent and Resource Sharing
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Dr. Eyal Maori (eyalmm@gmail.com).
Experimental Model and Subject DetailsReproductive mini hive system
Caged fertile queen bees, together with approximately 250 worker bees, were placed in mini hives (26x17x15.5 cm polystyrene) fitted with four mini-combs each. The mini hives were sealed and placed in a temperature-controlled room (28°C) for three days in which the combs were constructed and queen-workers recognition had been established. During the first three days, the bees were fed on a mixture of 33% honey and 67% sucrose powder (candy). Next, the mini hives were transferred into two net-houses separating between dsRNA treated, and untreated mini hives. The bees were free to fly within the net-houses and to forage for water from buckets. The first 14 days were an adaptation period, during which the colonies were fed on demand with candy, and pollen supplement patties (5 g each) were placed on top of the combs and replaced once a week. An established mini-colony was determined by at-least two constructed combs and egg-laying activity of the queen; only these hives were included in the experiment. During all the experiments, established colonies (two per treatment) were fed on pollen supplement patties (5 g each), and had an unlimited water supply.
Method DetailsdsRNA synthesis
The GFP sequence and primers as well as dsRNA synthesis procedure were described in Maori et al., 2009. Apis mellifera sequence that corresponds to the Vitellogenin mRNA (bases 4648-5084; NCBI accession number: NM_001011578.1) served as a template for dsRNA-Vg transcription, and was amplified by the following primers: 5′ TAATACGACTCACTATAGGGCGAAAAGCTTATCAGAAGGTGGAAGAAAA 3′; 5′ TAATACGACTCACTATAGGGCGACAATGTTTGTTAACGTTATGGTGGTA 3′ (T7 promoter in bold). Labeling of GFP-dsRNA was performed using a DIG RNA Labeling Kit (Roche) with a DIG-11-UTP concentration of 70 μm per reaction.
Molecular procedures
Cloning, transcription, RNA preparations, cDNA synthesis, RNA slot blot, Northern-blot, PCR and Proteinase K digestions were carried out according to published protocols or to the manufacturers’ instructions. Northern-blot analyses for the detection of labeled-dsRNA were conducted without using a probe; a standard DIG Northern-blot protocol was modified by omitting the probe-hybridization step. Briefly, 10 μl of pooled whole hemolymph extracts (detailed below) were electrophorased in the Northern-blot analyses of labeled dsRNA, using native 1.2% agarose gel. Slot- and Northern-blot analyses for detection of non-labeled dsRNA were probed with a DIG-labeled PCR probe (Roche Diagnostics Indianapolis, IN, USA) of a sequence corresponding to GFP-sequence used as template for the dsRNA-GFP synthesis.
Hemolymph extraction from bees
Young worker bees were collected from a single hive and immobilized within plastic straws. Individual bees were fed on 10 μl 50% sucrose solution (w/w) containing DIG-labeled dsRNA-GFP in a final concentration of 50 ng/μl. The control group was fed on 50% sucrose solution (w/w) only. To control dsRNA spillage and cuticle contamination, the immobilized bees were fed directly to their glossa and complete uptake was monitored. Hemolymph was collected from a small incision at the level of the 3rd dorsal tergite, using a microcapillary. The hemolymph of 10 workers was pooled per sample and stored at −80°C for later analysis.
Detection of dsRNA-GFP in royal jelly
Treated colonies were fed daily on 10 mL of 50% sucrose solution (w/w) containing dsRNA-GFP in a final concentration of 20 ng/μl (200 μg dsRNA-GFP per day), and untreated mini hives were fed daily on 10 mL of 50% sucrose solution (w/w) only. Feeding dsRNA-GFP was continued for five days and subsequently the colonies were fed only on 50% sucrose solution until the experiment’s end. At the beginning of day four, the queens were removed in order for the bees to rear new queens. We waited two hours, and placed artificial queen cells containing 1st-2nd instar larvae grafted from a different untreated hive. On day 6, we carefully removed intact 3rd-4th instar larvae with a fine paintbrush and harvested royal jelly. Royal jelly was harvested from five artificial queen cells, pooled and stored at −80°C.
Detection of dsRNA-GFP in worker jelly
Treated colonies were fed daily on 10 mL of 50% sucrose solution (w/w) containing dsRNA-GFP in a final concentration of 20 ng/μl (200 μg dsRNA-GFP per day), and untreated mini hives were fed daily on 10 mL of 50% sucrose solution (w/w) only. Feeding dsRNA-GFP was continued for eight days. On day five and eight, 4th-5th instar larvae were carefully removed from worker brood cells and checked for any physical damage. Worker jelly was collected by washing cells with nuclease free water to resuspend the low jelly quantity available. Samples were stored at −80°C.
Horizontal RNA transfer in the hive
Two treated colonies were fed daily on 15 mL of 50% sucrose solution (w/w) containing dsRNA-GFP in a final concentration of 20 ng/μl (300 μg dsRNA-GFP per day), and two untreated mini hives were fed daily on 15 mL of 50% sucrose solution (w/w) only. Feeding dsRNA-GFP was continued for seven days and subsequently all the colonies were fed only on 50% sucrose solution until the experiment’s end. The day in which dsRNA was first introduced represents ‘day-1’ (Figure 2A). Samples were collected in different time points, immediately frozen in liquid nitrogen and kept at −80°C for further analysis. Prior to RNA extraction, samples were rigorously washed with nuclease-free water.
Combs transfer experiment
For five days, treated colonies were fed daily on 10 mL of 50% sucrose solution (w/w) containing dsRNA-GFP in a final concentration of 20 ng/μl (200 μg dsRNA-GFP per day), and untreated mini hives were daily fed on 10 mL of 50% sucrose solution (w/w) only. On day six, combs containing 1st instar larvae were removed from an untreated hive and transferred either to another untreated or dsRNA-treated colony. On day ten, 5th instar worker larvae were collected, immediately frozen in liquid nitrogen and kept at −80°C for further analysis (Figure 2C). Prior to RNA extraction, samples were rigorously washed with nuclease-free water.
Queens swap experiment
For six days, one treated colony was fed daily on 10 mL of 50% sucrose solution (w/w) containing dsRNA-GFP in a final concentration of 20 ng/μl (200 μg dsRNA-GFP per day), and three untreated mini hives were fed daily on 10 mL of 50% sucrose solution (w/w) only. On day six, the queens from both untreated and dsRNA-treated mini hives were caged, and on day seven swapped with other queens as follows: the queen taken from the dsRNA-treated colony replaced the queen of an untreated mini hive, and a queen from the untreated colony replaced a queen from another untreated mini hive (Figure 2E). To allow acclimation, the queens remained caged for a day and then they were released to lay eggs for three days. 100 eggs were collected per treatment on day ten and stored at −80°C for further analysis.
Vitellogenin silencing by horizontal transfer of dsRNA-Vg
Nine plastic boxes were populated each with ca. 150 workers (collected nearby open brood cells) and a piece of comb containing eggs and young larvae. The plastic boxes were placed in an incubator at 34°C with 45%–55% humidity. The boxes were divided into three groups (three boxes per treatment) and were fed daily for eight days with: i) 100 μg dsRNA-Vg in 2ml 35% sucrose solution; ii) 100 μg dsRNA-GFP in 2ml 35% sucrose solution; and iii) 2ml 35% sucrose solution. Bees were fed with additional 4 mL 35% sucrose solution per day, to avoid hunger. In addition, the bees were routinely fed pollen-sugar supplement to encourage larvae rearing (70% pollen and 30% sugar powder). On day eleven, all cells were sealed and the adult bees were removed. The newly emerged bees were paired for ten days. The pairs were prepared by placing together two bees from different treatments that emerged at the same day and originated from the same hive. Paired bees were fed with 35% sucrose solution and pollen-sugar supplement. On day ten, we collected workers samples from each treatment for Vitellogenin expression analysis.
Real-time RT-PCR
qRT-PCR reactions were performed using the EvaGreen® qPCR Plus Mix kit (SOLIS BIODYNE) according to the manufacturer’s instructions. The Vitellogenin mRNA was amplified using primers: 5′-CCAAACTGGAACGGGACCTGC-3′ and 5′-TGTAGCTGTCAGTCGGCGTGC-3′. Calmodulin was used as a control for normalization using primers: 5′-CGAGAGAGAACGGTGGACTC-3′ and 5′-ATACGACACAGCCGACGAG −3′.
Sequencing of full length RNA from worker and royal jelly
Royal jelly was collected from queen brood cells containing 3rd instar larvae, and worker jelly was collected from 5th instar worker larvae brood cells; all brood belonged to untreated healthy looking hives. The jelly samples were immediately frozen in liquid nitrogen and kept at −80°C for further analysis. Total RNA extracted from worker and royal jellies was subjected to Tobacco Acid Pyrophosphatase (TAP) treatment using the Cap-Clip enzyme (CellScript). Modified RNA was purified again by standard Phenol/Chloroform extraction followed by Ethanol precipitation in the presence of Glycogen. The RNA pellet was taken up in 12 ul nuclease-free water. RNA quality and quantity were verified using Agilent RNA 6000 pico chip on a Bioanalyzer 2100 instrument. Libraries construction and sequencing were provided by Cambridge Genomic Services (Cambridge University). Briefly, full-length RNA libraries were prepared using the NEXTflex small RNA-seq kit v3 (Bioo Scientific) according to the manufacturer’s instructions with the following modifications; Adaptor ligation was performed at 16°C overnight (step A). The bead cleanup (step F) was performed following an amended one-sided purification protocol to retain also long fragments (no size selection protocol) as provided by the manufacturer. The final purification of the PCR product (step H1) followed also the amended protocol without size selection as provided by the manufacturer. Sequencing was performed in a NextSeq 500 instrument, in a 150 bp paired-end read run using the high output kit (300 cycle). RNA-seq data have been deposited in the ArrayExpress database at EMBL-EBI (www.ebi.ac.uk/arrayexpress) under accession number: E-MTAB-6557.
Quantification and Statistical AnalysisReal-time RT-PCR statistics
Statistical analyses were conducted with JMP statistical software version 13 (SAS Institute, Cary, NC, USA). Statistical significance was set at p < 0.05. We present both ΔCt values and the more intuitively understood measure of relative quantification ( = 2ˆ(-ΔΔCt)); to test for significant differences in relative expression of Vg in workers, one-way ANOVA was conducted on ΔCt values as previously described (Yuan et al., 2006). Significant differences between treatments were tested by the Tukey-Kramer (HSD) test.
Reads trimming and QC
An in-house script utilizing cutadapt (version 1.11) (Martin, 2011), fastx_trimmer (FASTX Toolkit 0.0.13) (Gordon and Hannon, 2010) and FastQC (version v0.10.1) (Bioinformatics, 2011) was used to trim raw fastq files. Briefly, Illumina NEXTflex small RNA 5′ and 3′ adaptor sequences were trimmed from paired-end fastq sample files, while retaining sequences that were at least 23 bp long. Due to the two color method of sequencing and other technical sequencing considerations, reads that were shorter than the total number of sequencing cycles had a poly-A tail followed by a poly-G tail that were both trimmed. Then, the four random index bases were trimmed from both ends of the sequences. Next, 5′ and 3′ Illumina NEXTflex small RNA adaptor sequences were trimmed anew and only sequences for which the length of both paired reads was at least 15bp long were retained. QC, performed with FastQC, revealed low quality bases at both ends of the reads. These low quality bases were trimmed using fastx_trimmer and an additional QC run indicated that all the samples are properly trimmed. Total number of reads that passed trimming and QC per library is summarized in Table S1. All in-house scripts have been deposited in Github and can be downloaded: https://doi.org/10.5281/zenodo.1302437.
Alignment of royal and jelly samples to the genome of A. mellifera
The genome of A. mellifera was downloaded from ensembl (release 32). Using in-house scripts (https://doi.org/10.5281/zenodo.1302437), the fastq files (R1 & R2) of each royal and worker jelly sample were converted into fasta files, which were then aligned to the A. mellifera genome using blat (Kent, 2002). The best hit for each R1 mapped read was matched with the best hit for its R2 mapped read. Therefore determining the final mapping of the read as well as the size of its matching RNA.
RNA size analysis
For each royal and worker jelly sample, we computed the length of the sequenced RNA’s independent of any genome alignment. That is, we used blat to align an R1 read and the reverse complement of its matching R2 read. We then calculated the length of the sequence the two reads span from the 5′ of R1 to the 3′ of R2. The size of reads that did not overlap is indicated as being greater than the length of the longer between the R1 and R2 reads (e.g., > 38 where 38 is the length in bp of the longer read). Using the described method, we managed to map the size of approximately 98.5% of all the reads in each of the samples.
Comparative RNA size analysis
To compare the size distributions of RNA in royal jelly and worker jelly samples, we computed a histogram of all RNA sizes present in the jellies and calculated the log of Reads Per Million, hereof denoted as RPM, using the following formula:(1)log(RPMi)=log((|sizei|∑i=15300|sizei|)∗106),where|sizei|=#ofRNA'sofsizei
To compute the log(RPM) of royal and worker jelly samples in general, we calculated the average number of RNA’s in each size group for all royal / worker jelly samples, that is |sizei¯¯|, and used it to calculate the general royal / worker log(RPM) value using the following formula:(2)log(RPMi)=log((|sizei¯|∑i=15300|sizei¯|)∗106)where|sizei¯|=average#ofRNA'sofsizeiinroyalorworkerjellysamples
Double-stranded honey bee RNA screen
Reads that mapped to the Apis mellifera genome were analyzed for the occurrence of putative honey bee dsRNA. Mapped paired-end reads were extracted from bam files using pairedBamToBed12 (https://github.com/Population-Transcriptomics/pairedBamToBed12). Reads were separated into plus and minus strands. Subsequently, we tested for pairwise overlaps of reads on the forward and reverse strand via bedtools intersect 2.27.1 (Quinlan and Hall, 2010), where the overlap was at least 25 nt, and the overhang on either side did not exceed 100 nt. We quantified the relative count of dsRNA candidates falling on each annotated gene (normalized by their library size). We further classified unique dsRNA candidates (characterized here by unique start and end coordinates) by noting their length distribution as well as classifying the annotation of their loci. The fraction of dsRNA reads in a sample was calculated by dividing the number of mapped reads in dsRNA regions to the total number of reads.
Comparative metagenomics of sequenced RNA populations
To identify the origin of royal and worker jelly RNAs we used blat (Kent, 2002) to map the RNA sequences against the genomes of A. mellifera, 13 viral genomes, 5243 bacterial genomes, 1038 fungi genomes and two plant genomes. The names and the sizes of the genomes are depicted in Table S1. Each of the royal and worker jelly RNA samples was mapped to the above mentioned genomes and the percentage of RNA mapped to each species was recorded as the number of RNA sequences mapped to each specie out of the total number of RNA sequences in the sample. The RNA in each sample was first searched against the bee genome and then against viral genomes, bacterial genomes, fungi genomes and finally against plant genomes. Each blat search was performed with the RNAs that did not match the previous genome they were searched against. That is, all the RNAs in each sample were searched against the bee genome but only RNA reads that did not match the bee genome were searched against viral genomes, only RNAs that did not match any of the viral genomes were searched against bacterial genomes etc. In order to calculate the number of RNA sequences mapped to each genome for each type of jelly (i.e., royal / worker), the data from the three royal jelly samples were treated as a single sample and the data from the three worker jelly samples were treated as a single sample. We then recorded the percentage of RNAs that were mapped to each species for the merged royal jelly samples and for the merged worker jelly samples. For individual sample analysis, we simply recorded the percentage of RNAs that were mapped to each species for each royal or worker jelly sample.
Characterization of jelly RNA corresponding to the honey bee genome
We used the ensembl annotation file for the genome of A. mellifera (release 32) to characterize the jelly RNA corresponding to the honey bee genome. Using an in-house script, we converted the ensembl annotation file into a bed format file and intersected it with a bed file version of the RNAs corresponding to the honey bee genome using intersectBed (bedtools version 2.17.0) (Quinlan and Hall, 2010). The number of RNAs in each annotation category was then recorded, excluding matches in which the overlap between the RNA and the annotation element was shorter than 7bp. To further characterize the jellies corresponding to the honey bee genome, we performed a blast search against two transposable elements (TE) databases, the TREP database (release 16) (Wicker et al., 2007) and RepBase (version 22.04) (Jurka, 1998). The number of RNAs that mapped to each of the TE categories was then recorded from which the percentage of each TE category was calculated.
Enrichment analysis of protein-coding jelly RNAs
Functional enrichment profiling of protein-coding RNAs in worker and royal jelly was conducted using gProfiler (https://biit.cs.ut.ee/gprofiler_archive2/r1760_e93_eg40/web/) with default g:GOSt parameters.
Detection of viral RNAs
To detect the presence of pathogen-related RNA, we performed a blast search (Kent, 2002) against 13 honey bee viral genomes (Table S1). That is, the RNAs of each of the six jelly samples were searched against each of the viral genomes and their genomic position, length (i.e., insert length), orientation (i.e., forward/reverse) and abundance, relative to each viral genome, were obtained (Table S4).
Abundance comparison of viral RNAs between royal and worker jellies
To enable a comparison between the abundance of virus-related RNAs in each of the 13 viral genomes, we calculated for each jelly (i.e., royal or worker) its log(TPM), where TPM = Transcripts Per Million, using the following formula:(3)TPMvirusi=Xvirusilvirusi(1∑virusjXvirusjlvirusj)∗106,whereXvirusj=rawRNAcountsforvirusjinroyalorworkerjellyandlvirusj=lengthofviralgenomej(inKb)
To analyze the size distribution of viral RNAs for each jelly type (i.e., royal / worker) or for each sample separately, we grouped the various insert lengths into nine distinct size groups: 15-19nt, 20-25nt, 26-31nt, 32-50nt, 51-75nt, 76-100nt, 101-125nt, 126-200nt and > 200nt and calculated the log(TPM) of each viral genome in each size group using the following formula:(4)TPMvirusi,size_groupj=Xvirusi,size_groupjlvirusi,size_groupj(1∑virusj,size_groupkXvirusj,size_groupklvirusj,size_groupk)∗106,whereXvirusj,size_groupk=rawRNAcountsforvirusjinsize_groupkandlvirusj'size_groupk=lengthofviralgenomej−middleofsize_groupk(inKb)
To evaluate the abundance of viral RNAs along the genome of each virus, we calculated the log(TPM) at each position along the viral genome of each virus using the following formula:(5)TPMvirusj,positioni=Xpositionilvirusj(1∑positionkXpositionklvirusj)∗106,whereXpositioni=rawRNAcountsatpositionialongthegenomeofvirusjandlvirusj=lengthofviralgenomej(inKb)
Data and Software Availability
The RNA-seq data have been deposited in the ArrayExpress database at EMBL-EBI under accession number: E-MTAB-6557. All in-house scripts have been deposited in Github and can be downloaded: https://doi.org/10.5281/zenodo.1302437. Other Software used in this work are all publicly available, with the links to them in the Key Resources Table. The raw imaging data, including images of gels and blots, have been deposited in Mendeley Data and can be accessed: https://data.mendeley.com/datasets/wv4s233nx7/draft?a=21d6454c-ead9-4a13-bcd8-f6395652424c. All the rest of the data are available in the manuscript or the supplementary materials.
Acknowledgments
Professor Ilan Sela, who died aged 81, was Professor of Virology and Molecular Biology at the Hebrew University of Jerusalem and was regarded as one of the pioneer virologists in Israel. Ilan was a creative mind and a passionate scientist. He found joy in understanding fundamental principles in biology and reveled in their translation into real world solutions. For more than five decades, he nourished generations of students with curiosity and admiration toward science. Ilan is greatly missed. The authors thank Dr. Ed Farnell and Dr Emily Clemente from Cambridge Genomic Services for their contribution in the RNA-seq protocol development. We thank Mr. John Rayner (Springwell Apiaries, UK) for his support on bee and jelly samples. We are grateful to Prof. Eric Miska, Prof. Jonathan Heeney, and their groups, as well as to Professor Sir Peter Lachmann, Mr. David Seilly, Dr. Laurence Tiley, and Dr. Shahar Avin (Cambridge University) for inspiring discussions and advice. This work was supported by internal B. Triwaks Bee Research Center funds (HUJI no. 0356043), the Orion Foundation (HUJI no. 0368598), and Israel Science Foundation (grant no. 1456/10) to S.S. Research was also supported by Marie Curie Intra-European Fellowship For Career Development (PIEF-GA-2010-274406), Leo Baeck Scholarship, and Herchel Smith Postdoctoral Fellowship (XXACC_AFGLTRB2) to E.M.
Author Contributions
E.M. conceptualized and supervised the study; E.M., Y.G., O.M., I.S., and S.S. designed experiments; E.M., Y.G., and H.K. performed the beehive and hemolymph experiments; O.M. performed the dsRNA-Vg experiment; E.M., Y.G., and R.M.K. performed the molecular analyses; V.K. and N.S. analyzed the sequencing data; and E.M. wrote the paper with input from all authors.
Declaration of Interests
The authors declare no competing interests. E.M. is a founder, Director, and Chief Science Officer of Tropic Biosciences Ltd. Tropic Biosciences had no role in the design of the study, in collection, analysis, and interpretation of data, nor in writing the manuscript.
Supplemental Information