Skip to main content
Download PDF

Abundant small RNAs in the reproductive tissues and eggs of the honey bee, Apis mellifera

BMC Genomics volume 23, Article number: 257 (2022) Cite this article

Abstract

Background

Polyandrous social insects such as the honey bee are prime candidates for parental manipulation of gene expression in offspring. Although there is good evidence for parent-of-origin effects in honey bees the epigenetic mechanisms that underlie these effects remain a mystery. Small RNA molecules such as miRNAs, piRNAs and siRNAs play important roles in transgenerational epigenetic inheritance and in the regulation of gene expression during development.

Results

Here we present the first characterisation of small RNAs present in honey bee reproductive tissues: ovaries, spermatheca, semen, fertilised and unfertilised eggs, and testes. We show that semen contains fewer piRNAs relative to eggs and ovaries, and that piRNAs and miRNAs which map antisense to genes involved in DNA regulation and developmental processes are differentially expressed between tissues. tRNA fragments are highly abundant in semen and have a similar profile to those seen in the semen of other animals. Intriguingly we also find abundant piRNAs that target the sex determination locus, suggesting that piRNAs may play a role in honey bee sex determination.

Conclusions

We conclude that small RNAs may play a fundamental role in honey bee gametogenesis and reproduction and provide a plausible mechanism for parent-of-origin effects on gene expression and reproductive physiology.

Peer Review reports

Introduction

Parents can influence the phenotype of their offspring not only through the genetic contribution that they pass on but also through so-called ‘parental effects’. Maternal effects, the best-studied parental effect, occur when the mother influences the phenotype of the offspring, regardless of the offspring’s own genotype. For example, strong maternal effects can be caused by maternal deposition into oocytes of proteins or messenger RNAs (mRNAs) that are critical for early development of the zygote. If the mother is unable to provide these critical proteins/mRNAs, the resultant embryo dies regardless of its genotype. Maternal effects have been recognised by biologists for many decades whereas paternal effects are less well studied, partially due to the relatively small size of sperm and the commonly held belief that sperm contribute only their DNA to the fertilised egg. However, it is increasingly clear that paternal effects exist, although the mechanisms by which they are mediated is often unclear [1].

Honey bees are a species in which parental manipulation of gene expression (such as via parental effects) is likely to evolve because i) females are polyandrous (mate with many males), ii) there is large investment in offspring colonies, the costs of which are borne by the parent colony collectively, but the majority of benefits accrue to a tiny number of offspring queens and their parents, iii) the value of male and female offspring to parents differs strongly [2,3,4,5,6]. Parental manipulation of offspring gene expression is particularly relevant for fathers, who die after mating [7], and are therefore unable to directly enhance the reproductive success of their daughters. Indeed, a drone’s only opportunity to enhance the reproductive outcomes of his daughters is to use epigenetic manipulations that provide his daughters with an advantage over the daughters of other males [4, 8,9,10]. For example, a male might attempt to minimise additional matings by the queen that he has just mated with, potentially increasing the reproductive success of his own daughters [11]. Further, it has been repeatedly shown that larvae reared as queens comprise a non-random set of patrilines, suggesting that some males enhance the probability that their daughters will be reared as reproductive queens, potentially by epigenetic means [12,13,14].

Small RNAs refer to non-coding RNA molecules between 18 and 50 nt long. Multiple classes of small RNAs are abundant in the semen of many animals and can alter offspring phenotypes when injected into zygotes [15,16,17,18,19,20,21]. Small RNAs are deposited maternally in Drosophila melanogaster [22, 23] and recent evidence suggests that small RNAs are also deposited paternally and can influence gene expression in the next generation [24]. Two examples demonstrate that biologically active small RNAs can be transferred from hymenopteran parents to offspring or from workers to larvae, providing a plausible mechanism by which epigenetic information could be transmitted between generations in honey bees. First, biologically active dsRNAs that confer acquired immunity are shared between honey bee generations via the glandular secretions that workers feed to larvae [25]. Second, female jewel wasps (Nasonia vitripennis - of the same taxonomic order as honey bees) include RNA molecules in their eggs that determine the sex of their offspring [26].

Small RNAs are classified according to their size and function. The first-discovered small RNA molecules, collectively known as microRNAs (miRNAs), are short (c.a. 22 nt) non-coding RNAs that typically bind to the 3′ UTR region of messenger RNAs, causing translational inhibition and/or mRNA degradation [27, 28]. Maternally inherited miRNAs are involved in sex determination during embryogenesis in Caenorhabditis elegans [29] and in mouse sperm, miRNAs degrade maternal mRNA stores in early zygotes to reprogram gene expression in the offspring [20]. Some miRNAs are maternally deposited in D. melanogaster, where they play important roles in embryogenesis [30,31,32].

Piwi interacting RNAs (piRNAs) are a class of small RNA between 24 and 32 nt in size, that typically have a uridine at the 5′ end (5′U bias) [33]. piRNAs interact with PIWI proteins, a class of Argonaute protein, to repress transposable element (TE) activity in the germline of animals during meiosis [34], including in insects [35, 36]. Additionally, piRNAs are involved in regulating gene expression in developing sperm cells and in somatic cells [37,38,39]. In Drosophila embryonic somatic cells, piRNAs destabilise/cleave target mRNAs to regulate embryonic development [40]. In Drosophila germ cells piRNAs interact with Aubergine and a germline specific poly(A) polymerase to facilitate the localisation of essential germline mRNAs to the germ plasm, where they are protected and can be passed from generation to generation [40]. Aubergine-piRNA-mediated epigenetic silencing of protein coding genes is well characterised in D. melanogaster, and can occur through piRNA induced silencing complexes [41,42,43] and/or spreading of piRNA-mediated heterochromatin into neighbouring loci [44].

tRNA fragments (tRFs) are small RNA molecules that are derived from the cleavage of mature transfer RNAs (tRNAs) [45]. The function of tRFs is mostly an enigma: some tRFs act similarly to miRNAs [46], while others interfere with global protein translation at the ribosome [47, 48]. Like miRNAs, tRF expression in Drosophila is age dependent. tRFs contain a seed region that has complementarity to 3’UTR regions of messenger RNAs and they interact with Argonaute proteins, suggesting that they form post-transcriptional RNA-induced silencing complexes (RISC) [49] and/or are involved in influencing mRNA stability and transport [50]. In mammalian sperm there is a global loss of piRNAs and an increase in tRFs and miRNAs that occurs during transit through the epididymis, a process that is essential for sperm maturation. tRFs have also been implicated in paternal transgenerational epigenetic inheritance [51,52,53].

The central roles played by miRNAs, piRNAs and tRFs in spermatogenesis, nurturing of germ line cells, somatic gene regulation and epigenetic inheritance in other species make them prime mechanistic candidates for parental manipulation of offspring development in honey bees [38, 54]. Previous studies of small RNAs in A. mellifera have primarily focussed on ovary, embryo and thorax tissues to investigate caste determination, oviposition or phylogenetic similarities [55,56,57,58]. Here, we focus on small RNAs in reproductive tissues to investigate their potential roles in gametogenesis, parent-of-origin effects and epigenetic inheritance. We find that honey bee reproductive tissues have distinctive small RNA profiles, with ovaries and eggs having a higher proportion of piRNAs relative to other tissues, semen having a higher proportion of tRNA fragments, and spermatheca and testes having a higher proportion of miRNAs. Our miRNA and piRNA target prediction suggests that germ cells utilise small RNAs to regulate processes during development (post-fertilisation), and in gametogenesis and sex determination. Our findings place small RNAs front and centre as the mechanism that mediates the parent- of-origin effects that are phenotypically observed in honey bees, but which cannot be ascribed to other epigenetic mechanisms, notably DNA methylation [59,60,61,62,63,64].

Results

Length distribution and biotype analysis of small RNAs

We surveyed the small RNA populations in the reproductive and germline tissues of Australian commercial honey bees (mainly A. m. ligustica heritage). We extracted RNA and generated small RNA libraries from a range of male (testes and semen) and female (spermatheca and ovary) reproductive tissues, as well as eggs from mated and virgin queens. The spermatheca is the organ in which sperm are stored, which a queen uses to fertilise female-destined eggs throughout her life. Sequencing of these libraries generated 12 to 35 Mb of reads per sample. On average 65% of reads mapped to the A. mellifera 4.5 genome. We first classified the small RNAs by size and biotype. Strikingly, the RNA populations varied greatly in size and biotype depending on their tissue of origin (Fig. 1A, Supplemental Fig. 1).

Fig. 1
figure 1

Small RNA species in the reproductive tissues of the honey bee. A. Size distribution of small RNAs between 13 and 43 nt in length mapped to the Amel 4.5 honey bee genome. Pie charts show the proportion of mapped reads that align to annotations of each biotype (tRNA, pre-miRNA, mRNA, other ncRNA and unannotated) for each tissue type. B. Euler plot illustrating number of novel miRNAs discovered in each reproductive tissue type and the overlap between tissue types. C. The 10 most abundant tRFs in semen and their relative abundance in other tissue types. tRF isoforms that originate from the same tRNA are shown as different shades of the same colour. D. A representation of the abundance of tRNA fragments derived from the tRNAGlyGCC in semen. The origin of each tRF is shown by its position around the molecule and relative abundance is indicated by colour

Full size image

Mapped reads from spermatheca and testis samples had a prominent peak at 21–23 nt which is the expected size of miRNAs. 58% and 82% respectively mapped to known precursor miRNAs from miRbase (Fig. 1A). In addition to honey bee miRNAs contained in miRbase, we identified 301 novel miRNAs using miRDeep2 [65]. The majority of these novel miRNAs were identified in the testes and spermathecal libraries (Fig. 1B) (Supplemental Fig. 1B, Supplemental Table S1). StemLoop qPCR [66] was conducted to validate three of the novel miRNAs that were differentially expressed between tissue types, according to small RNA sequencing. In all three cases the StemLoop qPCR results mirrored the expression levels suggested by the small RNA sequencing (Supplemental Fig. 1C).

In contrast to all other samples, semen samples had an abundance of mapped reads in the 32–33 nt range. Biotype analysis revealed that 60% of these reads mapped to sub-regions of tRNAs and were therefore tRNA fragments (tRFs) (Fig. 1A). Following the classifications of Loher and colleagues [67], tRFs were classified as being either 5′ or 3′ tRNA halves (tRH), 5′ or 3′ tRFs or internal tRFs (itRFs). Figure 1C shows the top 10 tRFs in semen, and their proportional abundance in other tissues. A full list of the relative abundance of tRFs across all tissue types is provided in Supplemental Table S2. The levels of individual tRFs appear to be highly regulated because their abundance differed greatly among tissues. For example, although the most abundant isodecoder, tRNAGlyGCC (Fig. 1D), generated 30–60% of all tRFs per tissue type, the first and second most abundant tRNAGlyGCC itRFs in semen were present at much higher levels in semen and spermatheca than in all other tissues (Fig. 1C). This difference indicates that tRF production is tissue-specific, and not simply a consequence of tRNA degradation. Strikingly, the high abundance of tRFs in semen was not observed in the spermatheca (Fig. 1A). This suggests that tRFs present in semen are carried in the seminal fluid and not by the spermatozoa, and/or that female contribution of small RNAs to the spermathecal fluid is abundant.

Mapped reads from egg and ovary samples showed a peak at 29 nt. This peak was also evident in the spermathecal samples but was relatively smaller than the larger peak at 21–23 nt (Fig. 1A). The high abundance of 29 nt length reads was evident in both mapped and unmapped reads (Supplemental Fig. 1B), suggesting that many were from repetitive regions that were not placed in the Amel 4.5 genome assembly. Egg and ovary samples also had a large proportion of total reads (~ 85%) that did not map to any recognised category of small RNA. The 29 nt reads are within the size range of piRNAs and nucleotide frequency plots generated using 26–31 nt reads showed a strong bias towards 1 U in ovaries and eggs (Supplemental Fig. 2). We therefore hypothesised that these 29 nt small RNAs were piRNAs. To generate a list of putative piRNAs for each tissue, reads mapping to other known ncRNAs (i.e. tRNAs, rRNAs, miRNAs) and reads below 26 nt or above 31 nt were removed.

piRNA cluster analysis

piRNAs are often transcribed from pericentromeric and telomeric heterochromatic regions, termed piRNA clusters, that are characterized by an abundance of transposable element (TE) remnants [68]. Using three algorithms: a custom script (CREST), proTRAC v2.4.2, and piClust, we identified 199 putative honey bee piRNA clusters totaling 3113.3 kB (1.25% of the genome) (Supplemental Table S3). Of the 199 clusters, 113 (57%) were located on unmapped scaffolds, probably due to the association with repetitive sequences that are hard to map [69]. In contrast to D. melanogaster [37, 70], we found that most honey bee piRNA clusters (171 of 199) are uni-directional and equally present on both genomic strands (Supplemental Table S3). This supports the suggestion of Wang et al. that species other than Drosophilids are incapable of dual-stranded piRNA cluster activity [58].

Clusters ranged in size from 1kB (Cluster 28) to 152 kB (Cluster 98) (Supplemental Table S3). Of the 199 clusters, 69 were identified in only one tissue, whereas only eight clusters were identified in all five tissues (Fig. 2A), demonstrating strong tissue-specific regulation of piRNA expression. The vast majority of putative piRNA reads mapped to clusters for ovaries (99%), eggs (99%), and spermatheca (90%). In contrast, only 56% of putative piRNAs mapped to clusters for testes and 23% for semen, suggesting that many of the piRNA-sized reads in these two tissues were not bona fide piRNAs, or are produced non-canonically. Of the piRNA reads that mapped to clusters, over 50% mapped to cluster 83 in eggs, ovaries and spermatheca (Fig. 2B). Cluster 83 contains many transposable element remnants (Fig. 2C) (predominantly large retrotransposon derivatives (LARDs)) and is the main piRNA-generating cluster in female tissues. This shows that honey bee piRNA clusters, like those in Drosophila [36], resemble transposon graveyards. In semen, 75% of clustered piRNA reads map to cluster 8, which contains one large peak of reads but does not overlap any TE or other genomic feature. This again suggests that TEs are not associated with piRNA-sized reads in semen.

Fig. 2
figure 2

Analysis of piRNAs in the reproductive tissues of the honey bee. A. Euler diagram showing the overlap of clusters identified within each individual tissue. B. Table of the proportion of clustered total/uncollapsed piRNAs that align to the top three most active piRNA clusters for each tissue. C. Coverage of piRNAs at Cluster 83 for each tissue. Y-axis values are log10 RPM (reads per million) normalised values of total reads. D. Left: Nucleotide frequency plots for unique piRNAs in ovaries, eggs, testes, semen and spermatheca samples. Right: The 1 U and 10A bias of each piRNA cluster, detected within each tissue individually. Red box indicates clusters with a 1 U bias. Cluster 25 (blue arrow) and cluster 8 (green arrow) are examples of the same piRNA cluster with different 1 U and 10A frequencies between tissue types. E. Feature mapping of total piRNAs sense and antisense to retrotransposons and DNA transposons (left) as well as introns and exons (right). All distributions compared with ovaries (chi-squared test). ** indicates p < 0.01, *** indicates p < 0.001, ns = non-significant

Full size image

Previous analysis of whole tissue honey bee larvae suggested that honey bees utilise ping pong biogenesis in the production of piRNAs [58]. To determine if this holds true for reproductive tissues, we assessed two signatures of ping pong biogenesis [37]: the prevalence of a 1 U/10A nucleotide bias and complementary sequence overlaps spanning 10 bp. Nucleotide frequency analysis of clustered piRNAs (in contrast to all piRNA-length reads) showed a strong 1 U bias in all tissues but no strong 10A bias (Fig. 2D), indicating that primary piRNAs are the most common overall. However, when 1 U/10A ratios were plotted for each individual cluster it was apparent that some clusters had both a 1 U and 10A bias, and that the extent of this bias differed between tissues (Fig. 2D). Furthermore, a 10 bp overlap between reads was enriched in testes and spermatheca, present in ovaries and eggs, but conspicuously absent in semen (Supplemental Fig. 3A). This suggests that some ping pong cycling is present in tissues except semen. We also measured piRNA phasing by plotting the 3′-5′ end distances between adjacent piRNAs. In D. melanogaster primary piRNA biogenesis, PIWI proteins use piRNAs as guides to fragment a pre-piRNA into a string of tail-to-head phased pre-piRNAs that are further processed into mature piRNAs [71]. We detected a significant signature of phased piRNA in semen (Z = 3.38, p < 0.001), but not in other tissues (Supplemental Fig. 3B). This indicates that primary piRNAs dominate in semen.

Overall, these data support previous findings that honey bees have ping-pong cycling capability [58], but strongly suggest that this activity is highly locus-specific, and not a general feature of piRNAs biosynthesis. piRNAs that are present in semen are mostly primary piRNAs – unsurprising, given that spermatozoa are mostly transcriptionally inert (see Discussion).

piRNA targeting

Most piRNA-length reads in maternal tissues mapped to TE-dense piRNA clusters. We hypothesized that semen and testes would have fewer TE-associated piRNA-length reads, as fewer piRNA clusters were identified in these tissues. As expected, we found that less than half of piRNA-length reads mapped to TEs in semen and testes, compared to around 65% in maternal tissues (Fig. 2E). Furthermore, TE-mapping piRNAs in ovaries, eggs and spermatheca predominantly mapped to retrotransposons in the sense direction, whereas piRNAs in semen predominantly mapped to DNA transposons (Fig. 2E), revealing clear tissue specificity in the types of TE targeted by piRNAs. Of the TEs that have been confidently assigned to TE families (more recently active TEs) [72], the class II transposon Mariner/TC1 was highly targeted by piRNAs in all tissues except semen, in which the class I LINE retrotransposon R2 was the most targeted (Table 1). A chi-squared test comparing the proportion of class I vs class II repetitive elements showed that only semen had a significantly different proportion to ovaries (chi-squared test, ovaries vs semen, p < 0.05, ovaries vs other tissues p > 0.05). Of putative piRNA reads that map to genes, most align sense to introns for all tissues except semen, in which the majority align sense to exons (Fig. 2E). The relative lack of piRNAs in semen and the stark difference of targeting and abundance of piRNA-length reads in paternal tissues compared to all other tissues suggests that piRNAs are a mechanism by which maternal effects could be mediated.

Table 1 Percentage of total piRNAs mapping to DNA transposons and retrotransposons for each tissue
Full size table

Differential expression and target analysis

Principal component analysis (PCA) revealed that miRNA, piRNA and tRF profiles differed between tissue types (Fig. 3A). To identify small RNAs that were contributing to the differences between tissues we performed pairwise differential small RNA expression analysis for each tissue (Fig. 3B, Table 2).

Fig. 3
figure 3

Differential expression analysis A. Principal component analysis (PCA) plots on tRF (top), miRNA (bottom left) and piRNA (bottom right) profiles between tissue types. B. Differential tRF, miRNA and piRNA expression for pairwise tissue comparisons (FDR < 0.01, FC > 5). Counts of differentially expressed and/or uniquely expressed small RNAs are on the x-axes. Red and grey bars extending to the right refer to the number of small RNAs that are uniquely (red) expressed or upregulated but not unique (grey) in tissue 2 relative to tissue 1. Conversely, the yellow bar extending leftward from the axis refers to the count of uniquely expressed small RNAs expressed in tissue 1 relative to tissue 2, the blue bar refers to the count of upregulated but not uniquely expressed sRNAs

Full size image
Table 2 Counts of differentially and uniquely expressed small RNAs for each tissue pairwise comparison
Full size table

Ovaries and eggs

There is preliminary evidence that queens may use supplementary epigenetic mechanisms to influence the sex of their eggs [73]. However, we found no difference in the abundance of small RNAs between eggs from mated queens collected from worker cells (fertilized) and eggs collected from drone cells (unfertilized) (Fig. 3B, row 2), indicating that at this age (< 24 h old) there is no difference in small RNA profiles between fertilised and unfertilised eggs. We therefore combined the data sets from mated queen eggs obtained from drone cells and worker cells into a single data set (mated queen eggs). We then compared sRNA profiles between eggs laid by virgin queens in worker cells with eggs laid in drone cells (both unfertilised). While we did not detect any differentially expressed (DE) miRNAs or tRFs we did detect 107 DE piRNAs (out of ~ 2.5 million piRNAs) (Fig. 3B, row 1). Given the minimal differences between these two classes of virgin queen eggs we also pooled these datasets together (virgin queen eggs) for greater statistical power for comparisons between mated queen eggs and virgin queen eggs. We identified 39 tRFs, 30 miRNAs and 3586 piRNAs that were DE between eggs from mated queens and eggs from virgin queens. The majority of these were upregulated in eggs laid by virgin queens (Fig. 3B, Fig. 4A). Due to these differences, we kept the virgin and mated queen egg datasets separate for comparisons with other reproductive tissues.

Fig. 4
figure 4

A. Volcano plots indicate the patterns of upregulated tRFs (left), miRNA (middle) and piRNAs (right) between eggs laid from virgin queens and eggs laid from mated queens. Grey dots are all small RNAs, red dots are those that are significantly DE (FDR < 0.01, FC > 5, counts > 5). Arrows refer to upregulation in either virgin or mated queen eggs. B. C. and D. Volcano plots showing differentially expressed tRFs (B), piRNAs (C) and miRNAs (D) between testes from 24 h pupae and 72 h pupae (left), 72 and 144 h (middle) and 24 and 144 h (right). Grey dots are all small RNAs, while dots above the horizontal dashed lines have a FDR p-value < 0.01 in the analysed comparison. Red dots in all panels indicate the significantly upregulated small RNAs from the 24–72 h comparisons (FDR < 0.01, FC > 5, counts > 5 in one tissue). E. Zetaview measurement of particles present in honey bee semen indicating particle size (nm) and concentration (particles/mL). F. Representative transmission electron micrographs of potential EVs in honey bee semen

Full size image

The small RNA expression profiles in the eggs derived from both virgin and mated queen eggs were more similar to those seen in ovaries than to other tissues. This is expected, as the ovaries are the site of oogenesis and egg maturation, and the ovaries that we collected came from laying queens containing hundreds of developing eggs. Nonetheless, differentially expressed tRFs, miRNAs, and piRNAs were observed between ovary and egg samples (Fig. 3B, Table 2). Differentially expressed small RNAs, particularly miRNAs and piRNAs, tended to be upregulated in ovaries. Any differentially expressed small RNAs must originate in either the maternal ovarian tissue, arise from differential transport of miRNAs into mature eggs, or arise due to differential expression early in egg development. The small RNAs upregulated in the ovaries are likely to be derived from the maternal ovarian tissue and not from the developing oocytes, while the small RNAs upregulated in the eggs must come from either active transport into the developing oocyte, or early expression in the embryo. Of the seven miRNAs upregulated in ovaries (relative to both mated and virgin eggs) three (mir-989, mir-278 and mir-12) have been implicated in ovary development or oogenesis in insects [74,75,76].

Developing Testes

During spermatogenesis meiosis occurs as pupation commences and ends with the appearance of clustered spermatids in red-eyed pupae [77]. To investigate whether paternally deposited small RNAs are produced during this process, we compared the testes samples from three different stages of pupal development; pink-eyed (24 h from the commencement of pupation), red-eyed (72 h) and brown-eyed (144 h). There were 139 differentially expressed tRFs between the developing testes of pink-eyed and red-eyed pupae and 140 differentially expressed tRFs between red-eyed and brown-eyed testes (Fig. 4B, Table 2). Strikingly, all differentially expressed tRFs were upregulated in the testes of red-eyed pupae, for each comparison. However, no tRFs were differentially expressed between pink-eyed and brown eyed pupae (Fig. 4B). These results imply a peak of tRF expression at around 72 h of pupal development with subsequent downregulation. Like tRF expression, there is a drastic increase in piRNA expression within the testes between 24 and 72 h of pupal development (Fig. 4C). These piRNAs show reduced abundance between 72 and 144 h, but not as strongly as tRFs. In contrast, miRNAs that are upregulated at 72 h testes relative to 24 h testes tend to remain upregulated in 72 h testes and 144 h testes, representing a trend towards increased miRNA production as the testis develop (Fig. 4D). Of tRFs that were significantly upregulated at 72 h (red-eyed pupae) relative to other pupal stages, many multi-mapped to eleven genes, two of which are important gametogenesis genes in D. melanogaster; Hephaestus (heph) [78, 79] and CUGBP Elav-like family member 2 (bru-2) [80]. 72 h post-fertilisation (red-eyed pupae) coincides with the end of meiosis and the start of spermiogenesis [77]: the striking changes in small RNA expression during this window strongly suggest that small RNAs play an important role in spermatogenesis. Might they also influence gene expression in the next generation? For further comparisons with other tissues we pooled all testes samples together.

Spermatheca, semen and testes

In order to determine whether paternal small RNAs are present in tissues relevant to the production of the next generation we considered the three environments encountered by honey bee sperm prior to fertilisation: the pupal testes, the adult seminal vesicles and the period of storage in the spermatheca [81, 82]. Strikingly, although tRFs made up 60% of small RNAs in semen and only 8% of spermathecal small RNAs by abundance (Fig. 1B), there are fewer DE tRFs between semen and spermatheca samples than between other tissue types, suggesting high similarity between semen and the semen-storage organ (Fig. 3B). In contrast, the stark difference in tRF abundance between semen and testes suggests that a large component of the tRF content in semen does not come from the testes but is a component of the seminal fluid, which is produced in the accessory glands [83]. This observation accords with the previous observation that tRFs are upregulated in the testes of red-eyed pupae but subsequently downregulated in brown-eyed pupae. The similarity in tRF expression between semen and spermatheca suggests that tRFs are a small RNA molecule uniquely placed to provide a paternal influence on early embryonic development.

Small RNAs are often present in extracellular vesicles (EVs) which have been shown to be present in semen of various organisms [84,85,86]. We performed particle tracking and electron microscopy of freshly collected honey bee semen and detected particles of the correct size and shape to be EVs (Fig. 4E, F). If EVs are indeed present in honey bee semen they are a plausible mechanism by which small RNAs could be transmitted paternally.

A different trend was observed for piRNAs between semen, spermatheca and testes. There were few DE piRNAs between semen and testes. In contrast, the spermatheca was strongly different, with approximately 28-fold higher expression of unique piRNAs than in testes (Fig. 3B, Table 2). The abundance of piRNAs in the spermatheca is most likely due to piRNA expression in the maternal cells of the spermatheca itself, or in the maternal component of the spermathecal fluid, but probably not from spermatozoa.

The miRNA compositions between semen, spermatheca and testes are different again. There are few miRNAs that are unique to the testes, while there are many uniquely expressed miRNAs in the semen and spermatheca. Although miRNAs make up 82% of small RNAs in the testes, and 58% in the spermatheca, they comprise only 9% of small RNAs in semen. There are slightly fewer miRNAs upregulated in semen relative to spermatheca. This could suggest that there is a miRNA component to the seminal fluid (although substantially less than the tRF contribution shown above), and a maternal miRNA contribution to the spermathecal fluid. The seven most upregulated miRNAs in spermatheca relative to semen are also more abundant in spermatheca relative to all other tissues but are of conspicuously low abundance in eggs (Supplemental Fig. 4). This suggests that these miRNAs regulate genes important to early embryogenesis. Indeed mir-317, mir-277 and mir-278 all play roles in insect development [75, 87, 88].

Female tissues compared to male tissues

We sought to investigate how small RNAs differ between male and female reproductive tissues, to gain insight into how they may utilise different strategies to regulate gene expression and genomic stability or contribute small RNAs to the next generation. Semen and spermatheca differentially express many more tRFs than do testes, ovaries, and eggs (Fig. 3B, Table 2, t-test p < 0.001), and thus tRFs are prime candidates for epigenetic marks passed by fathers to their offspring. However, at our significance thresholds we did not detect any tRFs present in semen, spermatheca and fertilised eggs from mated queens that were not also in unfertilised eggs from mated queens.

Overall, paternal tissues and spermatheca have more abundant and uniquely-expressed miRNAs relative to eggs and ovaries (Fig. 3B, Table 2, t-test p < 0.001), although not to the same degree as tRFs. GO analysis of miRNA target genes showed that DNA binding activity involved in transcription regulation was often enriched for both tissues of a comparison, suggesting that miRNAs that are differentially expressed between tissues might act on/suppress different DNA-transcription machinery to regulate the activity of large cellular pathways (Supplemental Table S5). In contrast to tRFs and miRNAs, many more piRNA species are strongly upregulated in maternal tissues relative to testes and semen (Fig. 3B, Table 2, t-test p < 0.05). This is further evidence that piRNAs are expressed in the maternal reproductive tissues and suppressed in the paternal reproductive tissues. In addition to mapping to TEs, piRNAs also map to genes. GO-term analysis of genes targeted by piRNAs in ovaries and eggs relative to semen and testes (and presumably therefore downregulated in maternal tissues) were enriched for terms related to glutamate receptor signalling. Conversely, genes targeted by piRNAs in semen and testes relative to eggs and ovaries were enriched for terms related to transcription and development. Genes/clusters which have notably different antisense-piRNA targeting between tissue types are listed in Table 3, and the top 10 piRNA-targeted genes for each tissue pairwise comparison are listed in Supplemental Table S4. Full GSEA output is in Supplemental Table S5. The targeting of genes by piRNAs suggests that piRNAs also regulate gene expression and developmental pathways.

Table 3 Notable piRNA-targeted genes and clusters
Full size table

Intriguingly, piRNA cluster 20 is much more active in ovaries, eggs and spermatheca relative to testes and semen (Table 3). Cluster 20 is a unidirectional cluster that maps antisense to most of the genes located in the cluster. Reads map along the length of the cluster and although there are transposable elements present in the cluster, they are not exclusively targeted. Cluster 20 overlaps the complementary sex determiner gene (csd), feminizer gene (fem) and GB47023/LOC408733 (pinocchio) (Fig. 5) [94]. csd, situated within the sex determination region on chromosome 3, determines the sex of the bee and arose from duplication of the feminizer gene [95]. Bees heterozygous for csd develop into females, whereas bees hemizygous at this locus develop as haploid males [94]. The high number of piRNAs at this locus suggest a role for piRNAs in sex determination.

Fig. 5
figure 5

Normalised log10 read counts of piRNAs at the csd locus (Cluster 20). Blue reads align sense and yellow reads antisense to Chromosome 3 as shown

Full size image

Discussion

Our pairwise comparisons of reproductive tissues of the honey bee has revealed strong differences in small RNA abundance, species and the genomic features that they target. Strikingly, these differences are strongly indicative that small RNAs function during gametogenesis, in the production, maturation and storage of semen, and in sex determination. Unique small RNA profiles are seen in fertilized and unfertilized eggs, and these differ from those seen in the spermathecal fluid. While our study was by its nature an exploratory survey rather than a test of specific hypotheses, it has revealed that piRNAs and miRNAs are present in ovaries and semen and that these target genes involved in DNA regulation and developmental processes. tRNA fragments are abundant in semen and have a similar profile to those seen in the semen of other animals. We conclude that small RNAs are likely to play an integral role in honey bee gametogenesis and reproduction and provide a plausible mechanism for parent-of-origin effects on gene expression and reproductive physiology. Below we discuss the potential function of the various small RNA classes in different tissues.

Small RNAs as potential mediators of paternal effects

tRFs made up around 60% of total small RNAs in semen, and although less abundant in spermatheca, the proportional abundance of tRF species was similar, suggesting that a large proportion of the tRFs present in the queen’s spermatheca are transferred there via spermatozoa or seminal fluid. The low frequency of tRFs in honey bee testis suggests that honey bee semen acquires tRFs from epididymosomes during transit through the epididymis, a process that coincides with spermatozoa acquiring their fertilising ability and motility. This process also occurs in mammals [51, 96] where seminal tRFs, particularly tRFGlyGCC, are acquired during epididymal transit. In mammals, tRFs are involved in the repression of embryonically expressed genes in developing offspring [51, 53, 96,97,98]. The high abundance of tRFGlyGCC in honey bee semen highlights a striking evolutionary conservation and suggests that tRFs have a role in mediating paternal effects in honey bees and potentially in other insects. For example, Liberti and colleagues [11] suggest that factors in seminal fluid cause queens to cease mating flights, a potential point of sexual conflict. Queens should prefer many matings to increase offspring diversity [99] and to reduce worker-worker conflicts [100], whereas drones can maximise their probability of fathering offspring if their partner queen’s mating flights cease after his mating.

The structure of a honey bee society is predicated by the fact that a single diploid queen mates with 20 or more haploid males to produce a colony of half-sisters. It is well established that not all sub-families are equal; some are more likely to be reared as queens [14], and workers of some rare subfamilies are much more likely to become reproductively active than average [13, 101, 102]. Further, offspring gene expression and phenotype are influenced by paternity in reciprocal crosses between honey bee subspecies (reviewed in Oldroyd and Yagound [64]).

Our work implicates seminal tRFs as a potential mechanism by which these effects may be mediated by affecting gene expression during the early development of the embryo. Further, the presence of extracellular vesicles in semen represents a potential route for transmission for tRFs from father to daughter. Nonetheless, it is difficult to see how tRFs produced by males could be transferred to eggs at an abundance that is sufficient to reliably influence embryonic development. In this context it is important to note that queens use stored sperm from the spermatheca to fertilise eggs over a lifespan that can exceed four years. This period of storage provides the opportunity for queens to ameliorate any paternal effects mediated by semen. We conclude that the function of paternal-origin tRFs in semen and spermatheca is far from clear.

Most novel miRNAs were identified in pupal testes and spermatheca. miR-210 was highly upregulated in testes, semen and spermatheca relative to ovaries and eggs. miR-210 is one of the most well-studied miRNAs across a range of species and is a master regulator of gene networks controlling neuronal development, circadian rhythm and photoreception [103]. In D. melanogaster, miR-210 overexpression during development causes visual defects, however miR-210 expression is also required for photoreceptor maintenance and survival via regulating acetyl coenzyme A synthetase expression [104, 105]. As strictly regulated expression of miR-210 is required for vision in D. melanogaster, miR-210 of paternal origin may be required for appropriate vision development in honey bees. miR-210 is also a candidate molecule to signal queens to cease mating flights after a successful mating [11]. Similarly, miR-34 is upregulated in testes, semen and spermatheca relative to eggs. Interestingly, miR-34 is maternally inherited in Drosophila and regulates synaptogenesis [106] and neuronal differentiation during embryogenesis [30].

Small RNAs may be maternally deposited in honey bees

In D. melanogaster piRNAs are maternally deposited in eggs, and this contribution is required for TE repression in the early embryonic development of offspring [22, 23, 107]. The high abundance of piRNAs that mapped to TEs in both egg and ovary samples suggests that the same is true in honey bees. During oogenesis in the honey bee ovary, intercellular cytoplasmic bridges allow cytoplasmic contents to be directly transferred from ovarian nurse cells into the oocyte in a process known as “nurse cell dumping” [108]. This suggests that, as is the case in Drosophila, piRNAs are transferred from honey bee nurse cells into the oocyte, priming it for TE repression upon zygotic genome activation. In addition to their role in TE repression, we found antisense gene-associated piRNAs in all reproductive tissues, indicating that some piRNAs function in gene regulation, rather than TE suppression, as has been shown in other species [109, 110]. Interestingly, we found that glutamate receptor signalling genes are enriched as targets of piRNAs in maternal tissues, whereas anatomical development genes are enriched as targets of piRNAs in semen and testes. This is suggestive that piRNAs are involved in regulating developmental processes post-fertilisation as has been shown in D. melanogaster [24].

Eggs of mated queens, (both fertilised and non-fertilised) showed different small RNA profiles to virgin queens. Queens have control over egg fertilisation, and can therefore control the sex of their offspring [111]. It is possible that virgin queens add more spermathecal fluid to their eggs than mated queens, perhaps as an attempt to compensate for a lack of sperm in their spermatheca, and that this alters the small RNAs present in the eggs of mated and virgin queens. In order to study this rigorously, freshly laid eggs (< 2 h) would need to be collected to detect differences before development begins. We note that in order to induce oviposition in virgin queens it is necessary to subject them to two rounds of CO2-induced narcosis [112] which is a potential confounding factor [113]. However, even if the changes in small RNA profiles that we observed in the eggs of virgin queens are caused by the CO2 treatment of the queens, it still indicates a maternal effect and communication of environmental factors between generations via small RNA molecules.

Control of transposable elements

We have identified evidence for ping-pong amplification of piRNAs at specific piRNA clusters and TEs (see below). However, the vast majority of piRNAs in honeybee reproductive tissues arise via primary biogenesis. The overall paucity of ping-pong signatures and the high proportion of sense-TE mapping reads at many piRNA clusters, including the highly active piRNA cluster 83, indicates that most piRNAs are derived from now-inactive TE remnants, as is the case in Drosophila [114]. While most piRNAs show no ping-pong signature, piRNAs targeting annotated TEs in female reproductive tissues mapped mostly to the PiggyBac and Mariner TE families, as shown previously [58], and had much stronger ping-pong signatures than piRNAs that mapped to unclassified TEs (probably inactive remnants) [72]. This suggests that PiggyBac and Mariner TEs are transcribed in female reproductive tissues, with potential for transposition, therefore requiring robust piRNA-mediated silencing that involves ping-pong amplification. In contrast, in semen the retrotransposon LINE R2 is the only TE highly targeted by piRNAs, suggesting that it is transcriptionally active. Interestingly, LINE-1 retrotransposons produce functional reverse transcriptase in murine spermatozoa [115], and have even been suggested as a mechanism by which extrachromosomal RNA carried by sperm could influence gene expression in progeny, perhaps by negatively regulating miRNA biogenesis [116]. As such, paternally inherited piRNAs that silence LINE-1 reverse transcription may permit biogenesis of certain miRNAs. Our data suggests that the presence of active LINE retrotransposons, and thus their reverse transcriptase, is highly conserved in spermatozoa, suggesting that drones may be able to influence embryonic development of their daughters via inhibition of LINE-1 activity.

A role for piRNAs in sex determination?

The unidirectional piRNA cluster 20 encompasses two key genes involved in the sex determination pathway: csd and feminizer. This cluster includes abundant piRNAs that align antisense to csd. We speculate that piRNAs transcribed from this locus may regulate the expression or splicing of genes involved in sex determination as has been demonstrated in the silk moth Bombyx mori, where a single piRNA transcribed from feminizer silences a gene that controls masculinization in male embryos [117] Involvement of piRNAs in sex determination is also known in C. elegans, where an X-chromosome derived piRNA silences a regulator of X-chromosome dosage compensation and sex determination [118]. In fertilized honeybee embryos, csd heterozygosity leads to an active CSD protein that splices feminizer transcripts. The spliced Fem protein is active and in turn splices doublesex transcripts, leading to female development [119]. Hemizygosity for csd leads to normal male development whereas homozygosity leads to diploid males that are eaten by workers early in embryonic development [120]. RNAi knockdown of feminizer and/or csd results in female to male phenotypic sex reversion [94] as does knock out of feminizer and/or doublesex [121]. It is not clear how heterozygosity at the csd locus results in fertile diploid females: one suggestion is that csd alleles form an active heterodimer [122]. Our finding that many piRNAs present in female reproductive tissues map antisense to the csd locus introduces another possibility; that piRNAs transcribed from one csd allele may affect expression or splicing of the second csd allele. In such a case a high level of sequence similarity between two csd alleles may be compensated for by piRNA-mediated silencing that results in female development despite a level of genetic homozygosity that would normally produce a male. Small RNA target prediction can be prone to false positives and experimental validation is essential, but it is important to note that piRNA target prediction, which relies on full complementarity along the length of the piRNA, is generally less prone to false positives than miRNA target prediction, which relies only on a short seed region.

Conclusions

The germline transmission of small RNAs in honey bees has not previously been studied. This study is the first to conduct pairwise comparisons of small RNA expression between reproductive tissues of maternal and paternal origin. We find evidence that, as with mammals, nematodes, and Drosophila, small RNAs are intimately involved in the regulation of gametogenesis and embryogenesis in the honey bee and provide a plausible but tentative pathway for parental manipulation of gene expression in offspring. Future studies will seek to experimentally validate the role of small RNAs at key loci such as the sex-determination locus (csd and feminizer) and to uncover the role of seminal tRFs. Studies such as these are essential to illuminate our understanding of how epigenetic mechanisms evolve to promote individual fitness.

Methods

Sample Collection

Semen

In February 2015 mature drones were sampled randomly from 3 colonies. Semen from 10 drones per colony was collected into sterile glass insemination tips using standard procedures used during honey bee artificial insemination [123]. Semen from each sample was expelled directly into 0.5 ml Trizol and the tube vigorously flicked to disperse the semen in the reagent. The preparation was then frozen at − 70°C until RNA extraction.

Eggs from unmated queens

In January 2015 we established four nucleus colonies and introduced one queen pupa to each using standard methods [123]. The sister queens were prevented from leaving the nucleus colonies via a ‘queen excluder’ grid placed over the entrance which was sufficiently large to allow the passage of workers, but too small to allow the passage of queens. When the queens were 7 days old we narcotized them with carbon dioxide for 10 min, once a day for two days, to induce oviposition, even though they had not mated [112]. Once the queens hard started laying we provided drone comb and worker comb to the colonies. Over several weeks we collected eggs 0–24 hr old into 0.5 ml Trizol from drone comb and worker comb. Eggs disintegrated on contact with the Trizol. Each sample (two from drone comb and two from worker comb) consisted of 59–200 eggs. Samples were frozen at − 70 °C after collection.

Eggs from mated queens

Eggs (59–200) were collected in September 2015 from three colonies headed by naturally mated queens of unknown age. Eggs were collected separately from drone and worker cells into 0.5 ml Trizol as above.

Spermathecae

In March 2016, 10 naturally-mated laying queens were removed from their colonies and taken to the laboratory. Queens were frozen, and the ovaries and spermatheca were removed by dissection as described in Dade (1977) [124]. Tissue was placed separately in 0.5 ml Trizol and immediately frozen.

Testes

In August 2016, a section of brood comb containing newly capped drone larvae was collected from a single colony and taken into the laboratory. White eyed drone pupae were extracted from their cells and placed in petri dishes lined with filter paper soaked in 12% glycerol. Pupae were incubated at 34.5°C, 50% relative humidity, and staged to collect developing testes at three time points. On day 1 (24 h post-incubation), pink-eyed pupae were collected and testes removed from the abdomen by piercing the cuticle with forceps and aspirating testes tissue using a P1000 pipette tip that had been widened by removing the end of the tip with sterile scalpel. Red-eyed pupae were collected on day 3 (72 h post-incubation) and brown-eyed pupae on day 6 (144 h post-incubation) and testes removed as described. Samples were frozen immediately on dry ice and stored at − 80°C prior to RNA extraction. Two individual drone pupae samples were used to generate two replicate small RNA libraries for each timepoint.

RNA extractions

Frozen tubes containing the sample and 500 μl Trizol were thawed on ice. The sample was then homogenized with a RNase free pellet pestle and RNA extracted using Trizol according to the manufacturer’s instructions (Life Technologies).

Transmission electron microscopy

For Extracellular Vesicles’ imaged using transmission electron microscopy (TEM), samples were diluted 1:200 in PBS and applied to formvar-coated grids with heavy carbon coating after fixation in glutaraldehyde (2% v/w in H2O) for 30 mins, room temperature, or overnight at 40 °C. Samples were visualized by negative staining with uranyl acetate (2% w/v in H2O), with images captured using a Joel JEM-2100 electron microscope.

Nanoparticle tracking analysis of isolated exosomes

Extracellular Vesicles’ were analyzed for size and concentration via the Zetaview instrument model Basic-PMX120 installed with a 405 nm laser diode (Particle Metrix). Vesicles were prepared with a 1:20,000 dilution in dPBS to an average 149 particles per frame. For each measurement, “high” (60fps) number of frames were captured for 11 positions; Camera sensitivity on the Zetaview was 80; Shutter speed was 1/100 s; After capture, analysis was performed by ZetaView Software version 8.05.12 SP1.

Small RNA Library prep and mapping

Small RNA libraries were prepared according to the manufacturer’s instructions (Illumina TruSeq Small RNA Library Prep Kit, Part # 15004197 Rev. G) unless stated otherwise. Input amount for library preparation was 1.5 μg of total RNA, quantified using a NanoDrop Spectrophotometer. PCR amplification was performed for 15 instead of 11 cycles. The amplified cDNA construct was purified on 6% TBE polyacrylamide gels. Gel bands between 145 bp and 160 bp in size, which correspond to adapter ligated 22 nt and 30 nt RNA fragments, were excised and purified. The final library was concentrated by ethanol precipitation and the pellet was resuspended in 10 μl 10 mM Tris-HCI, pH 8.5. Libraries were sequenced at the Australian Genome Research Facility using the Illumina HiSeq- Single Read (50/100) chemistry.

Small RNA libraries were mapped to the Amel4.5 genome using CLC Genomics (Qiagen). Biotype annotations were assigned using the Annotate and Merge tool. Differential expression analysis was performed using the empirical analysis of DGE tool in CLC Genomics, with FDR correction for multiple testing. To be included in our final pairwise DE gene lists, a differentially expressed small RNA had an FDR-corrected p value of < 0.01, a fold change of > 5, and greater than or equal to 5 counts in at least one tissue.

Novel miRNA identification

Novel miRNA prediction tool mirDeep2 [65] was used to predict novel miRNAs in each tissue (biological replicates pooled). CLC Genomics (Qiagen) was used to manually confirm expression of each predicted miRNA: those for which expression could not be confirmed (expression < 5 in all tissues) were removed.

Validation of miRNAs by stem–loop RT–PCR

For the validation of miRNAs by stem–loop RT–PCR independent samples were collected in September 2017. Total RNA of pink, red and brown testis, semen, eggs, ovaries and spermathecae was extracted, procedure for tissue collection and RNA extraction as described above. The protocol for miRNA validation by stem-loop RT-PCR was adapted from Ashby et al. (2015) [125] and Varkonyi-Gasic et al. (2007) [126]. Primers for stem-loop RT-PCR were designed with a miRNA primer design tool (Astridresearch: http://genomics.dote.hu:8080/mirnadesigntool/processor), using the base stack melting temp primer [127] (Table 4).

Table 4 Stem Loop primers
Full size table

For Reverse Transcription in a final volume of 20 ul 50 ng of total RNA was combined with 50 nM of the StemLoop RT primer, 0.25 mM dNTP mix and nuclease free water to 13.65 ul. The mixture was heated to 65 °C for 5 min followed by a 2-min incubation on ice. 6.35 ul of the enzyme mix consisting of 1x First First-Strand buffer, 10 mM DTT, 4 U RNaseOUT, 50 U SuperScript III Reverse Transcriptase (Life Technologies, Australia) was added and incubated in a Biorad Mycyler thermalcycler under following conditions: 16 °C for 30 min, 42 °C for 30 min, heat inactivated at 85 °C for 5 min and cooled to 4 °C. RT-PCR was performed as described by Ashby et al. (2015). For a 15 ul reaction 1 ul of undiluted cDNA was combined with 1x FastStart Universal Probe Master Mix with Rox (Roche Diagnostics, Australia), 50 nm of the forward and universal reverse primer, 10 nm of Universal Probe #21 (Roche Diagnostics, Australia). Cycling conditions were 95 °C for 5 min, 40 cyles of 95 °C for 10 s, 56 °C for 30 s, 72 °C for 10 s on a Applied Biosystems FAST 7500 real-time PCR machine using the normal ramp reaction mode.

piRNA Analysis

Prior to cluster identification, unmapped reads and reads that were annotated as other known RNAs (such as ribosomal and tRNAs, but excluding protein coding RNA) were removed. A custom script (CREST) was generated to identify genomic regions containing a high density of piRNA reads. After clusters were identified, those which contained less than 50% of reads of the appropriate size (26-31 nt) were filtered out. proTRAC v2.4.2 [128] and piClust [129] were run using the same dataset but excluding reads less than 26 nt or greater than 31 nt in length to ensure consistency across packages. Default parameters were used for proTRAC v2.4.2 and piClust was run using parameters; eps of 10,000, minread of 35 and a cut score of 2. Tissue-specific clusters which were either overlapping or within 1 kb of each other were merged and clusters less than 200 bp in length were removed. Cluster strandedness was determined by calculating the proportion of collapsed reads mapping to either strand. If one strand mapped more than double the number of piRNAs relative to the other strand, across all tissues, it was assigned as either sense or antisense. If there were fewer than double the number of reads (collapsed) on either strand, across all tissues, the cluster was assigned as a dual strand/bidirectional cluster. We identified enrichment of 10 nucleotide overlaps using the tool PPmeter (v 0.4) [130] and phased piRNA signatures were identified by calculating the 3′-5′ distances between putative piRNAs. piRNAs were mapped to the Amel 4.5 TE dataset [72], TE annotations are divided into TE class and TE superfamily as many repetitive elements could not be assigned to a superfamily (ie. Mariner, Copia etc) but were assigned to a class (ie LARD, TRIM etc).

Target prediction and gene ontology analysis

3’UTR sequences were obtained from OGS v3.2 (Beebase). Two miRNA target prediction algorithms, miRanda [131] and RNAhybrid [132] were used to identify miRNA 3’UTR-gene targets. Only miRNA-target interactions predicted by both algorithms were used for GO analysis. The targets were matched against each list of differentially expressed miRNA, for each pairwise tissue comparison, and the genes searched using the graphical gene set enrichment tool ShinyGO v 0.61 to obtain enriched biological process, cellular component and molecular function terms, KEGG pathways and promoter motifs [133]. To determine the cellular processes that piRNAs may regulate in each tissue type we sought to identify genes and TEs that are differentially targeted for each pairwise tissue comparison. Genes and TEs which had significantly more putative piRNAs mapping antisense to them in one tissue relative to another were assigned as differentially targeted. The Featurecounts command, as part of the Rsubread package (v. 1.22.2), was used to counts reads mapping to TEs and genes (including introns, exons, 3’UTR and 5’UTR). Intergenic reads were calculated by subtracting reads mapping to TEs and/or gene elements from the putative piRNA count (after length filtering).

Availability of data and materials

All raw and processed sequencing datasets generated in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE182720.

References

  1. Crean AJ, Bonduriansky R. What is a paternal effect? Trends Ecol Evol. 2014;29:554–9. https://doi.org/10.1016/j.tree.2014.07.009.

    Article  PubMed  Google Scholar 

  2. Haig D. Multiple Paternity and Genomic Imprinting. Genetics. 1999;151:1229–31. https://doi.org/10.1093/genetics/151.3.1229.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Haig D. The Kinship Theory of Genomic Imprinting. Annu Rev Ecol Syst. 2000;31:9–32. https://doi.org/10.1146/annurev.ecolsys.31.1.9.

    Article  Google Scholar 

  4. Queller DC. Theory of genomic imprinting conflict in social insects. BMC Evol Biol. 2003;3:15. https://doi.org/10.1186/1471-2148-3-15.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Burt A, Trivers R. Genes in Conflict. Harvard University Press; 2006.

  6. Haig D. Colloquium papers: Transfers and transitions: parent-offspring conflict, genomic imprinting, and the evolution of human life history. Proc Natl Acad Sci U S A. 2010;107(Suppl):1731–5. https://doi.org/10.1073/pnas.0904111106.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Winston ML. The biology of the honey bee. Cambridge: Harvard University Press; 1987.

    Google Scholar 

  8. Haig D. Intragenomic conflict and the evolution of eusociality. J Theor Biol. 1992;156:401–3. https://doi.org/10.1016/s0022-5193(05)80683-6.

    Article  CAS  PubMed  Google Scholar 

  9. Drewell RA, Lo N, Oxley PR, Oldroyd BP. Kin conflict in insect societies: a new epigenetic perspective. Trends Ecol Evol. 2012;27:367–73. https://doi.org/10.1016/j.tree.2012.02.005.

    Article  PubMed  Google Scholar 

  10. Pegoraro M, Marshall H, Lonsdale ZN, Mallon EB. Do social insects support Haig’s kin theory for the evolution of genomic imprinting? Epigenetics. 2017;12:725–42. https://doi.org/10.1080/15592294.2017.1348445.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Liberti J, Görner J, Welch M, et al. Seminal fluid compromises visual perception in honeybee queens reducing their survival during additional mating flights. Elife. 2019;8:e45009. https://doi.org/10.7554/eLife.45009.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Osborne KE, Oldroyd BP. Possible causes of reproductive dominance during emergency queen rearing by honeybees. Anim Behav. 1999;58:267–72. https://doi.org/10.1006/ANBE.1999.1139.

    Article  CAS  PubMed  Google Scholar 

  13. Châline N, Ratnieks FLW, Burke T. Anarchy in the UK: Detailed genetic analysis of worker reproduction in a naturally occurring British anarchistic honeybee, Apis mellifera, colony using DNA microsatellites. Mol Ecol. 2002;11:1795–803. https://doi.org/10.1046/J.1365-294X.2000.01569.X.

    Article  PubMed  Google Scholar 

  14. Withrow JM, Tarpy DR. Cryptic “royal” subfamilies in honey bee (Apis mellifera) colonies. PLoS One. 2018;13:e0199124. https://doi.org/10.1371/JOURNAL.PONE.0199124.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Vogel JP, Garvin DF, Mockler TC, et al. Genome sequencing and analysis of the model grass Brachypodium distachyon. Nature. 2010;463:763–8. https://doi.org/10.1038/nature08747.

    Article  CAS  Google Scholar 

  16. Rassoulzadegan M, Grandjean V, Gounon P, et al. RNA-mediated non-mendelian inheritance of an epigenetic change in the mouse. Nature. 2006;441:469–74. https://doi.org/10.1038/nature04674.

    Article  CAS  PubMed  Google Scholar 

  17. Gapp K, Jawaid A, Sarkies P, et al. Implication of sperm RNAs in transgenerational inheritance of the effects of early trauma in mice. Nat Neurosci. 2014;17:667–9. https://doi.org/10.1038/nn.3695.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Grandjean V, Fourré S, De Abreu DAF, et al. RNA-mediated paternal heredity of diet-induced obesity and metabolic disorders. Sci Rep. 2015;5:18193. https://doi.org/10.1038/srep18193.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Benito E, Kerimoglu C, Ramachandran B, et al. RNA-Dependent Intergenerational Inheritance of Enhanced Synaptic Plasticity after Environmental Enrichment. Cell Rep. 2018;23:546–54. https://doi.org/10.1016/j.celrep.2018.03.059.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Rodgers AB, Morgan CP, Leu NA, Bale TL. Transgenerational epigenetic programming via sperm microRNA recapitulates effects of paternal stress. Proc Natl Acad Sci. 2015;112:201508347. https://doi.org/10.1073/pnas.1508347112.

    Article  CAS  Google Scholar 

  21. Chen Q, Yan M, Cao Z, et al. Sperm tsRNAs contribute to intergenerational inheritance of an acquired metabolic disorder. Science. 2016;351(80):397–400. https://doi.org/10.1126/science.aad7977.

    Article  CAS  PubMed  Google Scholar 

  22. Brennecke J, Malone CD, Aravin AA, et al. An epigenetic role for maternally inherited piRNAs in transposon silencing. Science. 2008;322:1387–92. https://doi.org/10.1126/science.1165171.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Akkouche A, Grentzinger T, Fablet M, et al. Maternally deposited germline piRNAs silence the tirant retrotransposon in somatic cells. EMBO Rep. 2013. https://doi.org/10.1038/embor.2013.38.

  24. Lempradl A, Kugelberg U, Iconomou M, et al. Intergenerational metabolic priming by sperm piRNAs. bioRxiv. 2021:2021.03.29.436592. https://doi.org/10.1101/2021.03.29.436592.

  25. Maori E, Garbian Y, Kunik V, et al. A Transmissible RNA Pathway in Honey Bees. Cell Rep. 2019;27:1949–1959.e6. https://doi.org/10.1016/j.celrep.2019.04.073.

    Article  CAS  PubMed  Google Scholar 

  26. Verhulst EC, Beukeboom LW, van de Zande L. Maternal Control of Haplodiploid Sex Determination in the Wasp Nasonia. Science. 2010;328(80):620–3. https://doi.org/10.1126/science.1185805.

    Article  CAS  PubMed  Google Scholar 

  27. Cannell IG, Kong YW, Bushell M. How do microRNAs regulate gene expression? Biochem Soc Trans. 2008;36:1224–31. https://doi.org/10.1042/bst0361224.

    Article  CAS  PubMed  Google Scholar 

  28. Gebert LFR, MacRae IJ. Regulation of microRNA function in animals. Nat Rev Mol Cell Biol. 2019;20:21–37. https://doi.org/10.1038/s41580-018-0045-7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. McJunkin K. Maternal effects of microRNAs in early embryogenesis. RNA Biol. 2018;15:165–9.

    Article  Google Scholar 

  30. Soni K, Choudhary A, Patowary A, et al. miR-34 is maternally inherited in Drosophila melanogaster and Danio rerio. Nucleic Acids Res. 2013;41:4470–80. https://doi.org/10.1093/nar/gkt139.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Lee M, Choi Y, Kim K, et al. Adenylation of Maternally Inherited MicroRNAs by Wispy. Mol Cell. 2014;56:696–707. https://doi.org/10.1016/j.molcel.2014.10.011.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Kugler J-M, Chen Y-W, Weng R, Cohen SM. Maternal Loss of miRNAs Leads to Increased Variance in Primordial Germ Cell Numbers in Drosophila melanogaster. G3 Genes|Genomes|Genetics. 2013;3:1573–6. https://doi.org/10.1534/g3.113.007591.

    Article  CAS  PubMed Central  Google Scholar 

  33. Gunawardane LS, Saito K, Nishida KM, et al. A slicer-mediated mechanism for repeat-associated siRNA 5′ end formation in Drosophila. Science. 2007;315(80):1587–90. https://doi.org/10.1126/science.1140494.

    Article  CAS  PubMed  Google Scholar 

  34. Girard A, Sachidanandam R, Hannon GJ, Carmell MA. A germline-specific class of small RNAs binds mammalian Piwi proteins. Nature. 2006;442:199–202. https://doi.org/10.1038/nature04917.

    Article  PubMed  Google Scholar 

  35. Anand A, Kai T. The tudor domain protein kumo is required to assemble the nuage and to generate germline piRNAs in Drosophila. EMBO J. 2012;31:870–82. https://doi.org/10.1038/emboj.2011.449.

    Article  CAS  PubMed  Google Scholar 

  36. Brennecke J, Aravin AA, Stark A, et al. Discrete Small RNA-Generating Loci as Master Regulators of Transposon Activity in Drosophila. Cell. 2007;128:1089–103. https://doi.org/10.1016/j.cell.2007.01.043.

    Article  CAS  PubMed  Google Scholar 

  37. Czech B, Hannon GJ. One Loop to Rule Them All: The Ping-Pong Cycle and piRNA-Guided Silencing. Trends Biochem Sci. 2016;41:324–37. https://doi.org/10.1016/j.tibs.2015.12.008.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Thomson T, Lin H. The biogenesis and function of PIWI proteins and piRNAs: progress and prospect. Annu Rev Cell Dev Biol. 2009;25:355–76. https://doi.org/10.1146/annurev.cellbio.24.110707.175327.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Weick E-M, Miska EA. piRNAs: from biogenesis to function. Development. 2014;141:3458–71. https://doi.org/10.1242/dev.094037.

    Article  CAS  PubMed  Google Scholar 

  40. Dufourt J, Bontonou G, Chartier A, et al. piRNAs and Aubergine cooperate with Wispy poly(A) polymerase to stabilize mRNAs in the germ plasm. Nat Commun. 2017;8:1305. https://doi.org/10.1038/s41467-017-01431-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Rouget C, Papin C, Boureux A, et al. Maternal mRNA deadenylation and decay by the piRNA pathway in the early Drosophila embryo. Nature. 2010;467:1128–32. https://doi.org/10.1038/nature09465.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Barckmann B, Pierson S, Dufourt J, et al. Aubergine iCLIP Reveals piRNA-Dependent Decay of mRNAs Involved in Germ Cell Development in the Early Embryo. Cell Rep. 2015;12:1205–16. https://doi.org/10.1016/j.celrep.2015.07.030.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Wang C, Lin H. Roles of piRNAs in transposon and pseudogene regulation of germline mRNAs and lncRNAs. Genome Biol. 2021;22:1–21.

    Article  Google Scholar 

  44. Lee YCG. The Role of piRNA-Mediated Epigenetic Silencing in the Population Dynamics of Transposable Elements in Drosophila melanogaster. PLoS Genet. 2015;11:e1005269. https://doi.org/10.1371/journal.pgen.1005269.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Keam S, Hutvagner G. tRNA-Derived Fragments (tRFs): Emerging New Roles for an Ancient RNA in the Regulation of Gene Expression. Life. 2015;5:1638–51. https://doi.org/10.3390/life5041638.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Haussecker D, Huang Y, Lau A, et al. Human tRNA-derived small RNAs in the global regulation of RNA silencing. RNA. 2010;16:673–95. https://doi.org/10.1261/rna.2000810.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Ivanov P, Emara MM, Villen J, et al. Angiogenin-induced tRNA fragments inhibit translation initiation. Mol Cell. 2011;43:613–23. https://doi.org/10.1016/j.molcel.2011.06.022.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Kim HK, Fuchs G, Wang S, et al. A transfer-RNA-derived small RNA regulates ribosome biogenesis. Nature. 2017;552:57–62. https://doi.org/10.1038/nature25005.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Karaiskos S, Naqvi AS, Swanson KE, Grigoriev A. Age-driven modulation of tRNA-derived fragments in Drosophila and their potential targets. Biol Direct. 2015;10:1–13. https://doi.org/10.1186/s13062-015-0081-6.

    Article  CAS  Google Scholar 

  50. Göktaş Ç, Yiğit H, Coşacak Mİ, Akgül B. Differentially expressed tRNA-derived small RNAs co-sediment primarily with non-polysomal fractions in Drosophila. Genes (Basel). 2017;8:333. https://doi.org/10.3390/genes8110333.

    Article  CAS  Google Scholar 

  51. Sharma U, Conine CC, Shea JM, et al. Biogenesis and function of tRNA fragments during sperm maturation and fertilization in mammals. Science. 2016;351:391–6. https://doi.org/10.1126/science.aad6780.

    Article  CAS  PubMed  Google Scholar 

  52. Chen Q, Yan W, Duan E. Epigenetic inheritance of acquired traits through sperm RNAs and sperm RNA modifications. Nat Rev Genet. 2016;17:733–43. https://doi.org/10.1038/nrg.2016.106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Cropley JE, Eaton SA, Aiken A, et al. Male-lineage transmission of an acquired metabolic phenotype induced by grand-paternal obesity. Mol Metab. 2016;5:699–708. https://doi.org/10.1016/j.molmet.2016.06.008.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. He Z, Kokkinaki M, Pant D, et al. Small RNA molecules in the regulation of spermatogenesis. Reproduction. 2009;137:901–11. https://doi.org/10.1530/REP-08-0494.

    Article  CAS  PubMed  Google Scholar 

  55. Lewis SH, Quarles KA, Yang Y, et al. Pan-arthropod analysis reveals somatic piRNAs as an ancestral defence against transposable elements. Nat Ecol Evol. 2018;2:174–81. https://doi.org/10.1038/s41559-017-0403-4.

    Article  PubMed  Google Scholar 

  56. Pires CV, De Paula Freitas FC, Cristino AS, et al. Transcriptome analysis of honeybee (Apis Mellifera) haploid and diploid embryos reveals early zygotic transcription during cleavage. PLoS One. 2016;11:e0146447. https://doi.org/10.1371/journal.pone.0146447.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Chen X, Ma C, Chen C, et al. Integration of lncRNA–miRNA–mRNA reveals novel insights into oviposition regulation in honey bees. PeerJ. 2017;5:e3881. https://doi.org/10.7717/peerj.3881.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Wang W, Ashby R, Ying H, et al. Contrasting Sex-and Caste-Dependent piRNA Profiles in the Transposon Depleted Haplodiploid Honeybee Apis mellifera. Genome Biol Evol. 2017;9:1341–56. https://doi.org/10.1093/gbe/evx087.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Remnant EJ, Ashe A, Young PE, et al. Parent-of-origin effects on genome-wide DNA methylation in the Cape honey bee (Apis mellifera capensis) may be confounded by allele-specific methylation. BMC Genom. 2016;17:1–14. https://doi.org/10.1186/s12864-016-2506-8.

    Article  CAS  Google Scholar 

  60. Guzman-Novoa E, Hunt GJ, Page RE, et al. Paternal effects on the defensive behavior of honeybees. J Hered. 2005;96:376–80. https://doi.org/10.1093/jhered/esi038.

    Article  CAS  PubMed  Google Scholar 

  61. Wu X, Galbraith DA, Chatterjee P, et al. Lineage and Parent-of-Origin Effects in DNA Methylation of Honey Bees (Apis mellifera) Revealed by Reciprocal Crosses and Whole-Genome Bisulfite Sequencing. Genome Biol Evol. 2020;12:1482–92. https://doi.org/10.1093/gbe/evaa133.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Galbraith DA, Kocher SD, Glenn T, et al. Testing the kinship theory of intragenomic conflict in honey bees (Apis mellifera). Proc Natl Acad Sci U S A. 2016;113:1020–5. https://doi.org/10.1073/pnas.1516636113.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Kocher SD, Tsuruda JM, Gibson JD, et al. A Search for Parent-of-Origin Effects on Honey Bee Gene Expression. G3 (Bethesda). 2015;5:1657–62. https://doi.org/10.1534/g3.115.017814.

  64. Oldroyd BP, Yagound B. Parent-of-origin effects, allele-specific expression, genomic imprinting and paternal manipulation in social insects. Philos Trans R Soc Lond Ser B Biol Sci. 2021;376:20200425. https://doi.org/10.1098/rstb.2020.0425.

    Article  CAS  Google Scholar 

  65. Friedländer MR, Mackowiak SD, Li N, et al. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40:37–52. https://doi.org/10.1093/nar/gkr688.

    Article  CAS  PubMed  Google Scholar 

  66. Hurley J, Roberts D, Bond A, et al. Stem-Loop RT-qPCR for MicroRNA Expression Profiling. Methods Mol Biol. 2011:33–52.

  67. Loher P, Telonis AG, Rigoutsos I. MINTmap: fast and exhaustive profiling of nuclear and mitochondrial tRNA fragments from short RNA-seq data. Sci Rep. 2017;7:41184. https://doi.org/10.1038/srep41184.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Andersen PR, Tirian L, Vunjak M, Brennecke J. A heterochromatin-dependent transcription machinery drives piRNA expression. Nature. 2017;549:54–9. https://doi.org/10.1038/nature23482.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Huang X, Fejes Tóth K, Aravin AA. piRNA Biogenesis in Drosophila melanogaster. Trends Genet. 2017;33:882–94. https://doi.org/10.1016/j.tig.2017.09.002.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Akkouche A, Mugat B, Barckmann B, et al. Piwi Is Required during Drosophila Embryogenesis to License Dual-Strand piRNA Clusters for Transposon Repression in Adult Ovaries. Mol Cell. 2017;66:411–419.e4. https://doi.org/10.1016/j.molcel.2017.03.017.

    Article  CAS  PubMed  Google Scholar 

  71. Gainetdinov I, Colpan C, Arif A, et al. A Single Mechanism of Biogenesis, Initiated and Directed by PIWI Proteins, Explains piRNA Production in Most Animals. Mol Cell. 2018;71:775–790.e5. https://doi.org/10.1016/j.molcel.2018.08.007.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Elsik CG, Worley KC, Bennett AK, et al. Finding the missing honey bee genes: lessons learned from a genome upgrade. BMC Genomics. 2014;15:86. https://doi.org/10.1186/1471-2164-15-86.

    Article  PubMed  PubMed Central  Google Scholar 

  73. Oldroyd BP, Allsopp MH, Gloag RS, et al. Thelytokous parthenogenesis in unmated queen honeybees (Apis mellifera capensis): central fusion and high recombination rates. Genetics. 2008;180:359–66. https://doi.org/10.1534/genetics.108.090415.

    Article  PubMed  PubMed Central  Google Scholar 

  74. Kugler J-M, Verma P, Chen Y-W, et al. miR-989 is required for border cell migration in the Drosophila ovary. PLoS One. 2013;8:e67075. https://doi.org/10.1371/journal.pone.0067075.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Song J, Li W, Zhao H, et al. MicroRNA let-7 and miR-278 regulate insect metamorphosis and oogenesis via targeting juvenile hormone early response gene Krüppel-homolog 1. Development 145:dev170670. 2018. https://doi.org/10.1242/dev.170670.

  76. Macedo LMF, Nunes FMF, Freitas FCP, et al. MicroRNA signatures characterizing caste-independent ovarian activity in queen and worker honeybees (Apis melliferaL.). Insect Mol Biol. 2016;25:216–26. https://doi.org/10.1111/imb.12214.

    Article  CAS  PubMed  Google Scholar 

  77. Lago DC, Martins JR, Dallacqua RP, et al. Testis development and spermatogenesis in drones of the honey bee, Apis mellifera L. Apidologie. 2020;51:935–55. https://doi.org/10.1007/s13592-020-00773-2.

    Article  CAS  Google Scholar 

  78. Robida M, Sridharan V, Morgan S, et al. Drosophila polypyrimidine tract-binding protein is necessary for spermatid individualization. Proc Natl Acad Sci U S A. 2010;107:12570–5. https://doi.org/10.1073/pnas.1007935107.

    Article  PubMed  PubMed Central  Google Scholar 

  79. Sridharan V, Heimiller J, Robida MD, Singh R. High throughput sequencing identifies misregulated genes in the drosophila polypyrimidine tract-binding protein (hephaestus) mutant defective in spermatogenesis. PLoS One. 2016;11:e0150768. https://doi.org/10.1371/journal.pone.0150768.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Dasgupta T, Ladd AN. The importance of CELF control: molecular and biological roles of the CUG-BP, Elav-like family of RNA binding proteins. Wiley Interdiscip Rev RNA. 2012;3:104. https://doi.org/10.1002/WRNA.107.

    Article  CAS  PubMed  Google Scholar 

  81. Bishop GH. Fertilization in the honey-bee. I. The male sexual organs: Their histological structure and physiological functioning. J Exp Zool. 1920;31:224–65. https://doi.org/10.1002/jez.1400310203.

    Article  Google Scholar 

  82. Bishop GH. Fertilization in the honey-bee. II. Disposal of the sexual fluids in the organs of the female. J Exp Zool. 1920;31:266–86. https://doi.org/10.1002/jez.1400310204.

    Article  Google Scholar 

  83. Snodgrass RE. The anatomy of the honey bee /: Government Printing Office; 1910.

    Book  Google Scholar 

  84. Vojtech L, Woo S, Hughes S, et al. Exosomes in human semen carry a distinctive repertoire of small non-coding RNAs with potential regulatory functions. Nucleic Acids Res. 2014;42:7290–304. https://doi.org/10.1093/nar/gku347.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Zhao F, Cheng L, Shao Q, et al. Characterization of serum small extracellular vesicles and their small RNA contents across humans, rats, and mice. Sci Rep. 2020;10:1–16. https://doi.org/10.1038/s41598-020-61098-9.

    Article  CAS  Google Scholar 

  86. Chan JC, Morgan CP, Adrian Leu N, et al. Reproductive tract extracellular vesicles are sufficient to transmit intergenerational stress and program neurodevelopment. Nat Commun. 2020;11:1–13. https://doi.org/10.1038/s41467-020-15305-w.

    Article  CAS  Google Scholar 

  87. Shen Z-J, Liu Y-J, Zhu F, et al. MicroRNA-277 regulates dopa decarboxylase to control larval-pupal and pupal-adult metamorphosis of Helicoverpa armigera. Insect Biochem Mol Biol. 2020;122:103391. https://doi.org/10.1016/j.ibmb.2020.103391.

    Article  CAS  PubMed  Google Scholar 

  88. Yang H, Li M, Hu X, et al. MicroRNA-dependent roles of Drosha and Pasha in the Drosophila larval ovary morphogenesis. Dev Biol. 2016;416:312–23. https://doi.org/10.1016/j.ydbio.2016.06.026.

    Article  CAS  PubMed  Google Scholar 

  89. McQuibban GA, Lee JR, Zheng L, et al. Normal Mitochondrial Dynamics Requires Rhomboid-7 and Affects Drosophila Lifespan and Neuronal Function. Curr Biol. 2006;16:982–9. https://doi.org/10.1016/j.cub.2006.03.062.

    Article  CAS  PubMed  Google Scholar 

  90. Bier E, Jan LY, Jan YN. rhomboid, a gene required for dorsoventral axis establishment and peripheral nervous system development in Drosophila melanogaster. Genes Dev. 1990;4:190–203. https://doi.org/10.1101/gad.4.2.190.

    Article  CAS  PubMed  Google Scholar 

  91. Moberg KH, Mukherjee A, Veraksa A, et al. The Drosophila F box protein archipelago regulates dMyc protein levels in vivo. Curr Biol. 2004;14:965–74. https://doi.org/10.1016/j.cub.2004.04.040.

    Article  CAS  PubMed  Google Scholar 

  92. Shcherbata HR, Althauser C, Findley SD, Ruohola-Baker H. The mitotic-to-endrocycle switch in Drosophila follicle cells is executed by Notch-dependent regulation of G1/S, G2/M and M/G1 cell-cycle transitions. Development. 2004;131:3169–81. https://doi.org/10.1242/dev.01172.

    Article  CAS  PubMed  Google Scholar 

  93. Kanamori A, Nakayama J, Fukuda MN, et al. Expression cloning and characterization of a cDNA encoding a novel membrane protein required for the formation of O-acetylated ganglioside: A putative acetyl-CoA transporter. Proc Natl Acad Sci. 1997;94:2897–902. https://doi.org/10.1073/PNAS.94.7.2897.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. Beye M, Hasselmann M, Fondrk MKK, et al. The Gene csd Is the Primary Signal for Sexual Development in the Honeybee and Encodes an SR-Type Protein. Cell. 2003;114:419–29. https://doi.org/10.1016/s0092-8674(03)00606-8.

    Article  CAS  PubMed  Google Scholar 

  95. Hasselmann M, Gempe T, Schiøtt M, et al. Evidence for the evolutionary nascence of a novel sex determination pathway in honeybees. Nature. 2008;454:519–22. https://doi.org/10.1038/nature07052.

    Article  CAS  PubMed  Google Scholar 

  96. Sharma U, Sun F, Conine CC, et al. Small RNAs Are Trafficked from the Epididymis to Developing Mammalian Sperm. Dev Cell. 2018;46:481–494.e6. https://doi.org/10.1016/j.devcel.2018.06.023.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  97. Schorn AJ, Gutbrod MJ, LeBlanc C, Martienssen R. LTR-Retrotransposon Control by tRNA-Derived Small RNAs. Cell. 2017;170:61–71.e11. https://doi.org/10.1016/j.cell.2017.06.013.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. Conine CC, Sun F, Song L, et al. Small RNAs Gained during Epididymal Transit of Sperm Are Essential for Embryonic Development in Mice. Dev Cell. 2018;46:470–480.e3. https://doi.org/10.1016/j.devcel.2018.06.024.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  99. Oldroyd BP, Fewell JH. Genetic diversity promotes homeostasis in insect colonies. Trends Ecol Evol. 2007;22:408–13. https://doi.org/10.1016/J.TREE.2007.06.001.

    Article  PubMed  Google Scholar 

  100. Ratnieks FLW. Reproductive Harmony via Mutual Policing by Workers in Eusocial Hymenoptera. 2015;132:217–36. https://doi.org/10.1086/284846.

  101. Oldroyd BP, Smolenski AJ, Cornuet J-M. Crozler RH (1994) Anarchy in the beehive. Nat. 1994;3716500(371):749. https://doi.org/10.1038/371749a0.

    Article  Google Scholar 

  102. Montague CE, Oldroyd BP. The evolution of worker sterility in honey bees: an investigation into a behavioral mutant causing a failure of worker policing. Evolution (N Y). 1998;52:1408–15. https://doi.org/10.1111/J.1558-5646.1998.TB02022.X.

    Article  Google Scholar 

  103. Cusumano P, Biscontin A, Sandrelli F, et al. Modulation of miR-210 alters phasing of circadian locomotor activity and impairs projections of PDF clock neurons in Drosophila melanogaster. PLoS Genet. 2018;14:e1007500. https://doi.org/10.1371/journal.pgen.1007500.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  104. Weigelt CM, Hahn O, Arlt K, et al. Loss of miR-210 leads to progressive retinal degeneration in Drosophila melanogaster. Life Sci Alliance. 2019;2:e201800149. https://doi.org/10.26508/lsa.201800149.

    Article  PubMed  PubMed Central  Google Scholar 

  105. Lyu J, Chen Y, Yang W, et al. The conserved microRNA miR-210 regulates lipid metabolism and photoreceptor maintenance in the Drosophila retina. Cell Death Differ. 2021;28:764–79. https://doi.org/10.1038/s41418-020-00622-w.

    Article  CAS  PubMed  Google Scholar 

  106. McNeill EM, Warinner C, Alkins S, et al. The conserved microRNA miR-34 regulates synaptogenesis via coordination of distinct mechanisms in presynaptic and postsynaptic cells. Nat Commun. 2020;11:1092. https://doi.org/10.1038/s41467-020-14761-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  107. Chambeyron S, Popkova A, Payen-Groschene G, et al. piRNA-mediated nuclear accumulation of retrotransposon transcripts in the Drosophila female germline. Proc Natl Acad Sci. 2008;105:14964–9. https://doi.org/10.1073/pnas.0805943105.

    Article  PubMed  PubMed Central  Google Scholar 

  108. Cavaliere V, Taddei C, Gargiulo G. Apoptosis of nurse cells at the late stages of oogenesis of Drosophila melanogaster. Dev Genes Evol. 1998;208:106–12. https://doi.org/10.1007/s004270050160.

    Article  CAS  PubMed  Google Scholar 

  109. Gou L-T, Dai P, Yang J-H, et al. Pachytene piRNAs instruct massive mRNA elimination during late spermiogenesis. Cell Res. 2014;24:680–700. https://doi.org/10.1038/cr.2014.41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  110. Zhang P, Kang JY, Gou LT, et al. MIWI and piRNA-mediated cleavage of messenger RNAs in mouse testes. Cell Res. 2015;25:193–207. https://doi.org/10.1038/cr.2015.4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  111. Ratnieks FLW, Keller L. Queen control of egg fertilization in the honey bee. Behav Ecol Sociobiol. 1998;44:57–61. https://doi.org/10.1007/s002650050514.

    Article  Google Scholar 

  112. Mackensen O. Effect of Carbon Dioxide on Initial Oviposition of Artificially Inseminated and Virgin Queen Bees. J Econ Entomol. 1947;40:344–9. https://doi.org/10.1093/jee/40.3.344.

    Article  CAS  PubMed  Google Scholar 

  113. Manfredini F, Brown MJF, Vergoz V, Oldroyd BP. RNA-sequencing elucidates the regulation of behavioural transitions associated with the mating process in honey bee queens. BMC Genomics. 2015;16:563. https://doi.org/10.1186/s12864-015-1750-7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  114. Senti K-A, Brennecke J. The piRNA pathway: a fly’s perspective on the guardian of the genome. Trends Genet. 2010;26:499–509. https://doi.org/10.1016/j.tig.2010.08.007.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  115. Giordano R, Magnano AR, Zaccagnini G, et al. Reverse Transcriptase Activity in Mature Spermatozoa of Mouse. J Cell Biol. 2000;148:1107–14. https://doi.org/10.1083/jcb.148.6.1107.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  116. Spadafora C. Sperm-Mediated Transgenerational Inheritance. Front Microbiol. 2017;8:2401. https://doi.org/10.3389/fmicb.2017.02401.

    Article  PubMed  PubMed Central  Google Scholar 

  117. Kiuchi T, Koga H, Kawamoto M, et al. A single female-specific piRNA is the primary determiner of sex in the silkworm. Nature. 2014;509:633–6. https://doi.org/10.1038/nature13315.

    Article  CAS  PubMed  Google Scholar 

  118. Tang W, Seth M, Tu S, et al. A Sex Chromosome piRNA Promotes Robust Dosage Compensation and Sex Determination in C. elegans. Dev Cell. 2018;44:762–770.e3. https://doi.org/10.1016/j.devcel.2018.01.025.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  119. Biewer M, Schlesinger F, Hasselmann M. The evolutionary dynamics of major regulators for sexual development among Hymenoptera species. Front Genet. 2015;6:124. https://doi.org/10.3389/fgene.2015.00124.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  120. Woyke J. What happens to diploid drone larvae in a honeybee colony. J Apic Res. 1963;2:73–5. https://doi.org/10.1080/00218839.1963.11100063.

    Article  Google Scholar 

  121. McAfee A, Pettis JS, Tarpy DR, Foster LJ. Feminizer and doublesex knock-outs cause honey bees to switch sexes. PLoS Biol. 2019;17:e3000256. https://doi.org/10.1371/journal.pbio.3000256.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  122. Beye M. The dice of fate: thecsd gene and how its allelic composition regulates sexual development in the honey bee, Apis mellifera. BioEssays. 2004;26:1131–9. https://doi.org/10.1002/bies.20098.

    Article  CAS  PubMed  Google Scholar 

  123. Harbo JR. Propagation and Instrumental Insemination. Bee Genet Breed. 1986:361–89.

  124. Dade HA. Anatomy and dissection of the honey bee. London: International Bee Research Association; 1977.

    Google Scholar 

  125. Ashby R, Forêt S, Searle I, Maleszka R. MicroRNAs in Honey Bee Caste Determination. Sci Rep. 2016;6:18794. https://doi.org/10.1038/srep18794.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  126. Varkonyi-Gasic E, Wu R, Wood M, et al. Protocol: a highly sensitive RT-PCR method for detection and quantification of microRNAs. Plant Methods. 2007;3:12. https://doi.org/10.1186/1746-4811-3-12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  127. Czimmerer Z, Hulvely J, Simandi Z, et al. A versatile method to design stem-loop primer-based quantitative PCR assays for detecting small regulatory RNA molecules. PLoS One. 2013;8:e55168. https://doi.org/10.1371/journal.pone.0055168.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  128. Rosenkranz D, Zischler H. proTRAC - a software for probabilistic piRNA cluster detection, visualization and analysis. BMC Bioinformatics. 2012;13:5. https://doi.org/10.1186/1471-2105-13-5.

    Article  PubMed  PubMed Central  Google Scholar 

  129. Jung I, Park JC, Kim S. piClust: A density based piRNA clustering algorithm. Comput Biol Chem. 2014;50:60–7. https://doi.org/10.1016/j.compbiolchem.2014.01.008.

    Article  CAS  PubMed  Google Scholar 

  130. Jehn J, Gebert D, Pipilescu F, et al. PIWI genes and piRNAs are ubiquitously expressed in mollusks and show patterns of lineage-specific adaptation. Commun Biol. 2018;1:1–11. https://doi.org/10.1038/s42003-018-0141-4.

    Article  CAS  Google Scholar 

  131. Betel D, Koppal A, Agius P, et al. Comprehensive modeling of microRNA targets predicts functional non-conserved and non-canonical sites. Genome Biol. 2010;11:R90. https://doi.org/10.1186/gb-2010-11-8-r90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  132. Kruger J, Rehmsmeier M. RNAhybrid: microRNA target prediction easy, fast and flexible. Nucleic Acids Res. 2006;34:W451–4. https://doi.org/10.1093/nar/gkl243.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  133. Ge SX, Jung D, Yao R. ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics. 2020;36:2628–9. https://doi.org/10.1093/bioinformatics/btz931.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

The authors thank Dr. Ros Gloag for helpful discussions about the csd locus.

Funding

This work was supported by Australian Research Council grant DP150100151 to BPO and AA. AA is supported by the Australian Research Council (DE140100199 and FT180100653).

Author information

Authors and Affiliations

  1. School of Life and Environmental Sciences, The University of Sydney, Sydney, NSW, 2006, Australia

    Owen T. Watson & Alyson Ashe

  2. BEE Laboratory, School of Life and Environmental Sciences, University of Sydney, Sydney, NSW, 2006, Australia

    Gabriele Buchmann, Emily J. Remnant, Boris Yagound & Benjamin P. Oldroyd

  3. Molecular Cardiology and Biophysics Division, Victor Chang Cardiac Research Institute NSW 2010, Darlinghurst, Australia

    Paul Young

  4. School of Mathematics and Statistics, The University of Sydney, Sydney, NSW, 2006, Australia

    Kitty Lo

  5. Department of Biochemistry and Chemistry, School of Agriculture, Biomedicine and Environment, La Trobe Institute for Molecular Science, La Trobe University, Bundoora, Victoria, 3086, Australia

    Mitch Shambrook & Andrew F. Hill

  6. Institute for Health and Sport, Victoria University, Footscray, VIC, Australia

    Andrew F. Hill

  7. Wissenschaftskolleg zu Berlin, Wallotstrasse 19, 14193, Berlin, Germany

    Benjamin P. Oldroyd

Authors
  1. Owen T. Watson
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Gabriele Buchmann
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Paul Young
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Kitty Lo
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Emily J. Remnant
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Boris Yagound
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Mitch Shambrook
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Andrew F. Hill
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Benjamin P. Oldroyd
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Alyson Ashe
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

A.A., O.W., and B.P. wrote the main manuscript text. G.B. carried out the small RNA library preparation and miRNA validation. A.A. performed the miRNA and tRF analysis. O.W. performed the piRNA target analysis and GSEA. A.A. performed the DESeq analysis. P.Y. wrote the piCREST script. K.L. assisted with data visualisation. M.S. and A.H. performed the EV analysis. E.R. and B.Y. provided critical comments on the manuscript.

Corresponding authors

Correspondence to Benjamin P. Oldroyd or Alyson Ashe.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

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

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Watson, O.T., Buchmann, G., Young, P. et al. Abundant small RNAs in the reproductive tissues and eggs of the honey bee, Apis mellifera. BMC Genomics 23, 257 (2022). https://doi.org/10.1186/s12864-022-08478-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08478-9

BMC Genomics

ISSN: 1471-2164

Contact us