1 Introduction

Division of labor is the hallmark of social insects, and in the highly eusocial bees, ants, and wasps, it comes in two modes, reproductive division of labor and division of labor concerning colony maintenance tasks. Reproductive division of labor means the monopolization of reproduction by a long-lived and highly reproductive queen caste, while the short-lived workers are subfertile, or even completely sterile. Division of colony tasks, in turn, is restricted to the workers, which, at the beginning of their adult life cycle, predominantly perform the less risky within-nest task, especially brood care, and, later in life, the higher risk foraging tasks.

The mechanisms underlying the two modes of division of labor are best understood in the Western honey bee, Apis mellifera, where queen pheromones, especially the mandibular gland compounds, as well as brood pheromones, inhibit the ovarian activity in workers (Plettner et al. 1993; Le Conte et al. 2001; Cardoso-Junior et al. 2020). The age-related behavioral maturation of the workers, marked by the nurse-to-forager transition, is regulated by the hemolymph titers of vitellogenin (Vg) and juvenile hormone (JH), which function as a mutual repressor system (Huang et al. 1991; Amdam et al. 2003; Piulachs et al. 2003) that can adaptively adjust the worker tasks to the colony needs (Amdam et al. 2005).

The mutual repressor configuration between Vg and JH in A. mellifera workers represents a major change with regard to their ancestral association in insects, where JH functions as a gonadotropin, stimulating the Vg production in the abdominal fat body of the females (Santos et al. 2019). This change in regulatory circuitry is also relevant for the queen caste, where JH no longer drives the continuous production of Vg (Hartfelder and Engels 1998). In the primitively eusocial bumble bee, Bombus terrestris, however, JH still acts as a gonadotropin (Bloch et al. 2000; Amsalem et al. 2014; Shpigler et al. 2014). Hence, the question arises as to whether the reconfiguration of the JH/Vg circuitry is a result of the transition in the level of eusociality from the primitively eusocial bumble bees to the highly eusocial honey bees, or whether this is an idiosyncratic change that occurred in the evolutionary history of the honey bees (Apini).

To address this question, we decided to look at the stingless bees (Meliponini), which, together with the honey bees (Apini), are the only two tribes of highly eusocial bees within the monophyletic corbiculate bees (Hedtke et al. 2013). The Meliponini are ideally suited for comparative studies because, different from the tribe Apini with only a single genus (Apis) and less than a dozen species, the stingless bees comprise over 60 genera and more than 500 described species (Michener 2007; Rasmussen and Cameron 2010). In the Neotropical region, the Meliponini are the native highly eusocial bees, and here, they are also the most diverse, not only in terms of species number (Rasmussen and Cameron 2010), but also in their biological characteristics. Besides variation in nesting sites and colony size, they are remarkably diverse with respect to worker reproduction, i.e., the workers’ contribution to the males that are produced in a colony (Hartfelder et al. 2006; Vollet-Neto et al. 2018).

In this respect, the species Frieseomelitta varia and Melipona quadrifasciata are of particular interest, as they represent extremes in the spectrum of worker reproduction. In F. varia, the workers’ ovaries undergo complete atrophy during pupal development, while in M. quadrifasciata, the workers produce trophic as well as reproductive eggs, particularly when they are middle-aged nurses (Hartfelder et al. 2006). The two species are also the currently only two Neotropical meliponids for which genomic information is available (Kapheim et al. 2015; de Paula Freitas et al. 2020), which considerably facilitates gene expression analyses.

With this in mind, we investigated in the workers of these two species of stingless bees the transcript levels of a set of candidate genes that are well-known for being related to the behavioral phenotypes of honey bee workers. We hypothesized that these genes would show variation in their expression levels related to the specificities of the life history and reproductive capacity in the stingless bees. Stingless bees exhibit an age polyethism similar to the honey bee, and life table data have been collected for several species. However, as mainly unpublished student theses, they became accessible only once they were compiled in the recently published book “Stingless Bees” (Grüter 2020).

We analyzed the expression of the following genes: vitellogenin (Vg), related to insect reproduction and division of labor in A. mellifera (Amdam et al. 2003; Piulachs et al. 2003); methyl farnesoate epoxidase (mfe), involved with juvenile hormone synthesis (Bomtorin et al. 2014); juvenile hormone esterase (jhe), related to juvenile hormone degradation (Bomtorin et al. 2014); malvolio (mvl) and foraging (for), both related to the nurse-to-forager transition (Ben-Shahar et al. 2004; Ben-Shahar, 2005; Thamm and Scheiner 2014); and, finally, the gene target of rapamycin (tor) encoding the serine/threonine kinase mTOR protein, which is related to nutrient sensing and metabolism in honey bee caste development (Patel et al. 2007), and to a lesser extent and with seasonal differences, also with the progression to foraging (Ament et al. 2008). Quantitative expression analyses were performed on the tissues of the head and abdomen in the three major life cycle stages: newly emerged workers, nurse-stage workers, and foragers. The considerable difference that we found between the two stingless bee species, and compared to A. mellifera and bumble bees, can be interpreted as molecular signatures of the life histories of these eusocial bees.

2 Materials and methods

2.1 Bees

Workers of Frieseomelitta varia and Melipona quadrifasciata were collected from colonies kept in the apiary of the Department of Genetics of the University of São Paulo in Ribeirão Preto, Brazil. Newly emerged workers were collected as soon as they hatched from brood cells; nurses were captured while building or supplying cells with larval food; and foragers were captured when returning to the nest with pollen on their hind legs. The bees were briefly anesthetized on ice before separating the body compartments. The heads and abdomens were then stored in TRIzol Reagent (Thermo Fisher Scientific, Waltham, MA, USA) and kept at − 80 °C until RNA extraction.

2.2 RNA extraction and cDNA synthesis

Total RNA was extracted using the standard protocol of the TRIzol manufacturer. RNA quantity and quality were assessed spectrophotometrically (NanoVue, GE Healthcare, Chicago, IL, USA). After treatment with RNase-free DNase (Thermo Fisher Scientific) to remove residual DNA, cDNA was synthesized from 1 µg of total RNA using the UltraScript 2.0 cDNA Synthesis Kit (PCR Biosystems, London, UK).

2.3 Quantitative analysis of relative gene expression

For the quantitative analysis of gene expression by the RT-qPCR methodology, we used species-specific primers to amplify the Vg, mfe, jhe, mvl, for, and tor genes. The primers were designed using Primer3Plus software (http://www.bioinformatics.nl/cgibin/primer3plus/primer3plus.cgi) based on the respective predictions of these genes in the F. varia (de Paula Freitas et al. 2020) and M. quadrifasciata (Kapheim et al. 2015) genomes. The primers were synthesized by Integrated DNA Technologies (IDT, Coralville, IA, USA). All genes, accession numbers, primer sequences, amplicon sizes, and primer efficiencies are listed in the Supplementary Table S1. Amplifications were performed using the 2 × qPCR SyGreen Mix (PCR Biosystems). The reactions consisted of 5 µL of SyGreen, 1.25 pM of each primer, 1 µL of 10 × diluted cDNA, and ultrapure water for a final volume of 10 µL. The detection was done in a StepOnePlus™ Real-Time PCR system (Applied Biosystems, Waltham, MA, USA) under the following conditions: 95 °C for 2 min followed by 40 cycles at 95 °C for 5 s and 60 °C for 10 s. Primer specificity was determined by dissociation curve analysis, and the standard amplification efficiencies of the primers (E) were estimated based on the slope values of the standard curve (E = 10(−1/slope)) generated from serial 1:10 dilutions of 1 µg of cDNA, consisting of a mixture of cDNAs from all F. varia or M. quadrifasciata samples.

The assays were carried out on seven biological replicates for F. varia (pool of three heads or abdomens of the same bees per biological replicate) and, respectively, on eight biological replicates (individual heads and abdomens) for M. quadrifasciata for each age class and body compartment. Each sample was analyzed as a technical triplicate. The ribosomal protein l32 (rpl32) and ribosomal protein s18 (rps18) genes were used for normalizing the relative expression of the target transcripts. These genes have been validated as the most appropriate endogenous control genes for RT-qPCR studies with stingless bees (Freitas et al. 2019). For each of the two species, always the same sample of the head from the newly emerged worker group was used as calibrator. The threshold cycle (Ct) values were used to calculate the relative expression of the genes of interest using the 2−ΔΔct formula (Livak and Schmittgen 2001).

2.4 Statistical analysis

The statistical analyses were performed using GraphPad software (GraphPad Software, Inc. Prism 7 for Windows, San Diego, CA, USA). First, the conformity of the data with normality was checked using the Shapiro–Wilk’s test, and since all the data were in conformity, the levels of relative gene expression for the three life cycle stages and the two body compartments were compared using two-way ANOVA tests followed by Bonferroni post hoc tests. The significance level was set at p < 0.05. Pearson’s correlation analyses were performed separately for the heads and abdomens of the two species, but jointly for the three developmental stages.

3 Results

3.1 Vitellogenin gene expression

The relative expression of the Vg gene showed significant variation in the head of the life cycle stages of adult F. varia workers, but not in the abdomen (Figure 1A; Supplementary Table S2). This contrasts strongly with the results observed for M. quadrifasciata, where the Vg mRNA levels were constantly low in the head, but significantly elevated in the abdomen of the nurse bees (Figure 1B; Supplementary Table S3). Remarkably, the dynamics of Vg expression in the head of F. varia, with the highest transcript levels in foragers, is exactly opposite to what is seen for the abdomen of M. quadrifasciata, where Vg expression is highest in the nurse bees.

Figure 1.
figure 1

Transcript levels of the vitellogenin (Vg) gene. Relative expression levels in the head (black dots) and abdomen (gray squares) of F. varia (A) and M. quadrifasciata (B). Data for newly emerged workers (NE), nurses (N), and foragers (F) are presented as means ± SEM, as well as individual data points (n = 7 for F. varia and n = 8 for M. quadrifasciata). Different letters indicate statistical significance (p < 0.05) analyzed by two-way ANOVA and Bonferroni post hoc tests.

Full size image

The elevated Vg mRNA levels in the head of newly emerged F. varia workers, and even more so in the foragers, strongly contrast with the low levels in nurses. This indicates that in the sterile workers of F. varia, Vg may have shifted from its ancestral reproductive function to an exclusively behavioral one. This contrast seen for the two stingless bee species, and the possible shift in Vg function, implicitly leads to the question as to how JH signaling may have been configured in the two species.

3.2 Methyl farnesoate epoxidase and juvenile hormone esterase gene expression

In both species, we found detectable levels of mfe transcripts only in the head (Figure 2A, C) and not in the abdomen (Figure 2B, D). This is consistent with the role of this gene in JH synthesis in the corpora allata, as specifically shown for the honey bee (Bomtorin et al. 2014).

Figure 2.
figure 2

Transcript levels of the methyl farnesoate epoxidase (mfe) and juvenile hormone esterase (jhe) genes. Relative expression levels in the head (black dots) and abdomen (gray squares) of F. varia (A, C) and M. quadrifasciata (B, D). Data for newly emerged workers (NE), nurses (N), and foragers (F) are presented as means ± SEM, as well as individual data points (n = 7 for F. varia and n = 8 for M. quadrifasciata). Different letters indicate statistical significance (p < 0.05), analyzed by two-way ANOVA and Bonferroni post hoc tests.

Full size image

For the heads of the two species, we found that in F. varia, the mfe mRNA levels were significantly lower in the head tissue of nurses compared to those of newly emerged workers and foragers (Figure 2A; Supplementary Table S2). For M. quadrifasciata, the analysis of the levels of mfe transcripts in the head showed no statistical difference among the three life cycle stages (Figure 2B; Supplementary Table S3).

Regarding the jhe gene, F. varia did not show significant expression differences in the head across the three life cycle stages, but in the abdomen, it was found higher expressed in nurses compared to newly emerged workers and foragers (Figure 2C). In M. quadrifasciata, the jhe mRNA levels did not differ in the heads of the three life cycle stages, but were significantly lower in the abdomen of the nurses (Figure 2D).

Hence, for M. quadrifasciata, the results indicate that the circulating JH levels in hemolymph are primarily regulated by the JH degrading activity of the JH esterase (Figure 2D) and not by JH synthesis (Figure 2B). Furthermore, the lower jhe expression in the abdomen of M. quadrifasciata nurses (Figure 2D) would be consistent with a gonadotropic function of JH, as the predicted higher JH levels in the nurses would coincide with their higher Vg expression (Figure 1B). In contrast, the higher jhe mRNA levels in the foragers would mean that their predicted JH titer would be low and, thus, inconsistent with a possible role of this hormone in driving the nurse-to-forager transition.

For F. varia, the statistically significant variation in mfe and jhe expression (Figure 2A, C) indicates that JH synthesis in the corpora allata and JH degradation in the abdomen both contribute to the regulation of the hemolymph JH levels. In terms of predicted JH titers, these are, thus, expected to be low in nurse bees and high in the foragers, i.e., consistent with a role for JH as a driver for the nurse-to-forager transition. Furthermore, the low JH levels in the nurses would explain their diminished Vg expression in the head (Figure 1A). Hence, taken together, the results for the JH-related gene expression in the two stingless bees indicates that in M. quadrifasciata, JH has retained its ancestral gonadotropic role, while in F. varia, it is likely involved in the nurse-to-forager transition.

3.3 mTOR serine/threonine protein kinase gene expression

The tor mRNA levels in the head of F. varia nurses were significantly lower when compared to those in the newly emerged workers and foragers. They also were continuously lower in the abdominal than in the head tissues (Figure 3A). In contrast, for M. quadrifasciata, no variation was seen, neither for the head nor the abdominal tor transcript levels (Figure 3B).

Figure 3.
figure 3

Transcript levels of the mTOR serine/threonine protein kinase (tor) gene. Relative expression levels in the head (black dots) and abdomen (gray squares) of F. varia (A) and M. quadrifasciata (B). Data for newly emerged workers (NE), nurses (N), and foragers (F) are presented as means ± SEM, as well as individual data points (n = 7 for F. varia and n = 8 for M. quadrifasciata). Different letters indicate statistical significance (p < 0.05) analyzed by two-way ANOVA and Bonferroni post hoc tests.

Full size image

As far as we know, we are the first to report on the expression of the tor gene for stingless bees, and to our surprise, we found that in both species, the expression was low in the abdomen, despite the fact that it contains the major portion of the fat body, which is the insects’ metabolic center. This indicates that the tor gene is probably only marginally related to reproduction in these stingless bees, despite its importance as a major nutrient sensing system in developmental processes and reproduction throughout the animal kingdom (Oldham 2011). The tor gene expression in F. varia heads, however, presented an interesting coincidence with the Vg and mfe expression dynamics, indicating a potentially positive regulatory link between TOR, JH, and Vg in the behavioral maturation of the workers of this species.

3.4 Malvolio and foraging gene expression

The mRNA levels of the mvl gene did not show statistical differences in the heads of the life cycle stages in F. varia, while in the abdomen, there was a statistically significant increase in the foragers (Figure 4A). For the for gene, we found an increased expression in foragers, both in the head and abdomen of F. varia (Figure 4C), which would be consistent with its putative function in the foraging behavior.

Figure 4.
figure 4

Transcript levels of the malvolio (mvl) and foraging (for) genes. Relative expression levels in the head (black dots) and abdomen (gray squares) of F. varia (A, C) and M. quadrifasciata (B, D). Data for newly emerged workers (NE), nurses (N), and foragers (F) are presented as means ± SEM, as well as individual data points (n = 7 for F. varia and n = 8 for M. quadrifasciata). Different letters indicate statistical significance (p < 0.05) analyzed by one-way ANOVA and Tukey’s post hoc test.

Full size image

In M. quadrifasciata, the expression dynamics for mvl were similar to those of F. varia, exhibiting a tendency to increase as the bees become foragers (Figure 4B). The for gene, however, showed exactly opposite expression patterns for the head vs. the abdomen, with statistically significant higher transcript levels in the head of the workers performing brood care functions (Figure 4D).

Taken together, it is only in the sterile workers of F. varia where the expression of the mvl and for genes is consistent with a function in behavioral maturation, while the elevated expression of for in the heads of M. quadrifasciata coincides with the nursing stage, i.e., just when the workers are reproductively most active. Also worthy of note is the finding that the dynamic range for the expression of the mvl gene is much broader than that of the for gene.

3.5 Pearson’s correlation analysis

For an in-depth statistical analysis on possible correlations in the transcript levels of the gene set, we performed Pearson’s correlation analyses. In the two species, these were done separately for the two body compartments, but jointly for the three life cycle stages. The results are shown in Figure 5 in the form of heat maps, with numerical values indicating the positive or negative correlation coefficients among pairs of genes. The statistical significance for the respective correlations is presented in the Supplementary Tables S4 and S5.

Figure 5.
figure 5

Heat map of pairwise Pearson’s product-moment correlations. Correlation coefficients for the six candidate genes, Vg, mfe, jhe, tor, mvl, and for, in the head and abdomen of Frieseomelitta varia (A, B) and for Melipona quadrifasciata workers (C, D).

Full size image

For F. varia, all statistically significant correlations were in the positive direction, both in the head and the abdomen (Figure 5A, B). In both body compartments, a strong correlation was seen between Vg, mfe, tor, and for expression, and in the abdomen, Vg transcript levels were also correlated with those of mvl. In contrast, for the jhe expression, which encodes the negative regulator for the hemolymph JH titer, we saw only a weak, non-significant negative correlation with the other genes. This indicates that the control, over JH synthesis, but not JH degradation, is likely a major driver in the life history of the workers of F. varia.

For M. quadrifasciata, we found a strong positive correlation between Vg and jhe expression in the head, while in the abdomen, the sign of the correlation is the inverse (Figure 5C, D); Supplementary Table S5. This is in agreement with the postulated positive regulation of Vg expression in the abdomen of the nurse-stage workers. A similar opposite directionality was also noted in the correlations between the Vg and mvl transcript levels. The mvl and for genes, as modulators of food search activity, show a particularly interesting pattern of correlations. In the head, they are both strongly positively correlated with the expression of Vg and jhe and also with each other, while in the abdomen, there is only a significant negative correlation between mvl and Vg expression.

The results of the correlation analysis, thus, further underline the conclusion drawn above that the reproductive potential of the M. quadrifasciata workers differentially affects the expression patterns and the apparent connectivity among these life cycle-relevant genes in their head versus abdominal tissues.

Although better than using whole bodies as sources for gene expression analyses, we are aware that whole heads and whole abdomens are also not ideal, because the head, besides the brain, also contains glands and a small amount of fat body and muscles, and the abdomen, in addition to the fat body and ovaries as metabolic and reproductive tissues, also contains the gut, abdominal ganglia, and the dorsal heart tube.

4 Discussion

Our comparative analysis on two species of stingless bees revealed that the reproductive potential of their workers has an important impact on the expression pattern of a set of candidate genes known to play important roles in the life cycles and life histories of social insects. First of all, and quite surprisingly, we found that the vitellogenin-encoding gene is strongly expressed in the heads of the sterile F. varia workers, especially in the forager stage (Figure 1A). Hence, in this species, Vg could be related to behavioral maturation, whereas in the reproductive workers of M. quadrifasciata, it is strongly expressed in the abdomen of nurse bees (Figure 1B), contributing to egg production, as expected.

A strength in our study is that we analyzed the mRNA levels of these genes for both the head and the abdomen of the same samples of bees, while in most studies on honey bees, Vg expression was either investigated for the whole body, or, when done separately for the head or abdomen, these were not the same sample sets. Though unexpected, Vg expression in the head of honey bees has been documented in several studies and by different methods. One of the first studies (Corona et al. 2007) showed that Vg expression in the head occurs predominantly in fat body cells located in the head of queens, while the expression in worker heads is much lower. At the same time, Seehuus et al. (2007) used an immunogold protocol and detected Vg protein in head fat body cells and also in the hypopharyngeal glands of workers, lending support to the previously proposed hypothesis of the exploitation of Vg for the organization of social life (Amdam et al. 2003). Subsequently, a vitellogenin receptor encoding gene was found expressed in head tissue, primarily in the hypopharyngeal glands of nurse bees (Guidugli-Lazzarini 2008). While none of these studies could propose a behavioral, i.e., a neural function for Vg, the study by Münch et al. (2015) was the first to detect Vg immunoreactivity in a set of glial cells. However, they could not confirm the presence of Vg mRNA by in situ hybridization. Only recently, a single-cell transcriptomics analysis (Zhang et al. 2022) did provide the lacking evidence, showing the expression of Vg in a glial-cell subtype.

The results that we now present for the workers of two stingless bee species expand the facets of the relationship of Vg with the social evolution of the bees. The elevated Vg expression in the abdominal tissue in M. quadrifasciata nurses is consistent with the fact that in this life cycle stage, they produce trophic eggs as nutrition for the queen, as well as reproductive eggs that contribute to the production of males of the colony (Velthuis et al. 2005). Similarly, the continuously low levels of Vg expression in the abdomens of F. varia workers also come to no surprise, because their ovaries are functionally completely degraded during pupal development (Boleli et al. 1999), hence leaving no room for a reproductive function.

Our results for M. quadrifasciata, and those of Dallacqua et al. (2007) for the closely related species M. scutellaris, are also in line with the ancestral reproductive function of this protein. This is also the case for the bumble bees, which are the sister group to the stingless bees in the monophyletic corbiculate bees (Porto and Almeida 2021). For Bombus terrestris, Amsalem et al. (2014) found that the head Vg transcript levels did not differ significantly between nurses and foragers. For the abdomen, Jedlicka et al. (2016) found that 9-day-old workers had Vg mRNA levels similar to those of egg-laying queens, and in another bumble bee species, B. hypocrita, the Vg transcript levels in the abdominal fat body showed a temporal pattern that was very similar to that of honey bee workers, with a peak around 10 to 20 days, followed by a progressive decline (Li et al. 2010). Hence, the temporal dynamics of Vg expression in the abdominal tissues of bumble bees and the Melipona stingless bees are highly similar, and like in stingless bees, the bumble bee workers can become potential egg layers, especially at the end of the annual colony cycle (Duchateau and Velthuis 1989; Amsalem et al. 2015).

In the genus Frieseomelitta, however, Vg has apparently shifted from its ancestral reproductive function to a potentially purely behavioral one in the worker caste. Yet, this shift occurred in an apparently inverse directionality compared to the honey bees, where Vg transcript levels in the abdominal fat body, and consequently the Vg hemolymph titer, drop during the transition from the nurse to the forager stage (Hartfelder and Engels 1998; Amdam et al. 2003; Corona et al. 2007). In the honey bee workers, this transition defines the terminal phase of their life span (Guzmán-Novoa et al. 1994; Prado et al. 2021), in consequence of the onset of immunosenescence caused by the diminished Vg levels (Amdam et al. 2005). As mentioned, this stands in contrast with the elevated Vg expression that we saw in the heads of F. varia foragers, compared to the nurse bees (Figure 1A). This raises the question as to how worker life span may be regulated in these bees, including a possibly altered configuration of the Vg/JH regulatory module, or via a JH-independent regulation of Vg via a DNA methyltransferase (Cardoso-Júnior et al. 2017a).

The JH hemolymph titer in insects is determined by the rate of JH biosynthesis in the corpora allata and the rate of JH degradation in the hemolymph and target tissues. A critical step in JH-III biosynthesis is the epoxidation of farnesoic acid or JH acid, which is catalyzed by the enzyme methyl farnesoate epoxidase (MFE) (Goodman and Cusson 2012). In bees, as in most insects, the gene encoding this enzyme is exclusively expressed in the corpora allata of the retrocerebral complex (Bomtorin et al. 2014). Circulating JH, in turn, is degraded either by the removal of the methyl ester group by a juvenile hormone esterase (JHE), or the opening of the epoxy ring by a juvenile hormone epoxy hydrolase (JHEH). In the honey bee, it is the JHE that is the primary JH degrading enzyme (Mackert et al. 2008, 2010). Hence, to understand the potential regulatory dynamics for the JH levels, we analyzed the expression of the mfe and the jhe gene in the two stingless bee species (Figure 2).

The results indicated that in M. quadrifasciata, JH is likely to have a gonadotropic function, stimulating Vg expression (Figures 1B and 2B, D). This claim is supported by the levels of mfe expression and JH titer measurements in queens and workers of the closely related species M. scutellaris (Cardoso-Júnior et al. 2017b) and also by the data on mfe and jhe expression in M. interrupta workers (Brito et al. 2021). The maintenance of the ancestral gonadotropic JH function in the genus Melipona finds its explanation in the fact that the workers make a major contribution to the production of males in their colonies (Velthuis et al. 2005; Vollet-Neto et al. 2018). In contrast, in the sterile workers of F. varia, JH is likely to drive the behavioral nurse-to-forager transition, but in a way that would involve the regulation of Vg expression in the brain (Figures 1A, 2A, C) and not, like in A. mellifera, as an inhibitor of Vg expression in the abdominal fat body (Piulachs et al. 2003; Amdam et al. 2003). Furthermore, and taking into account the gonadotropic function of JH in the primitively eusocial bumble bees (Bloch et al. 2000; Amsalem et al. 2014), the comparison with the highly eusocial A. mellifera and the two stingless species suggests that the likely factor that shaped the regulatory circuitry between JH and Vg was not the evolutionary transition from a primitively to a highly eusocial lifestyle, but the status of worker sterility in the respective species.

We figured that the nutrient sensing insulin/insulin-like signaling (IIS) and TOR pathway could be a mechanism underlying this connection between Vg and JH seen in the stingless bees, as has already shown for the honey queens and workers (Corona et al. 2007) and also for the age polyethism and immunosenescence in the workers (Ament et al. 2008; Cardoen et al. 2011; Aurori et al. 2014; Ihle et al. 2019; Lourenço et al. 2019). This connection also involves changes in the gustatory responsiveness of the bees (Wang et al. 2012; Mengoni Goñalons et al. 2016) that is related to the switch in diet preference at the nurse-to-forager transition, accompanied by a major change in the metabolic profile (Cardoen et al. 2011). We found that for the completely sterile Frieseomelitta workers, the expression profiles for Vg, mfe, and tor (Figures 1A, 2A, 3A) are exactly aligned for both the head and the abdomen, while for the reproductive Melipona workers, they either go in opposite direction (Figure 1B) or are not variable (Figures 2B, 3B). This finding is of interest in view of the recent evidence from B. terrestris, showing that the brain-reproduction metabolic tradeoff is regulated by JH in the bumble bee, but not in A. mellifera (Shpigler et al. 2020). When putting the data on tor gene expression in the two stingless bees in this context, our interpretation is that this disruption in the proposed brain-reproduction tradeoff may, again, be more related to the degree of worker sterility than to the actual level of sociality.

The two candidate genes related to foraging behavior in the honey bee, malvolio (mvl) and foraging (for), also showed unexpected and contrasting expression patterns in the two stingless bee species (Figure 4A, B). The malvolio protein and its human ortholog Nramp1 belong to a group of divalent ion transporters that are associated with food choice behavior in Drosophila (Orgad et al. 1998; Søvik et al. 2017) and also in mammals (Evans et al. 2001). In A. mellifera (Ben-Shahar et al. 2004) and A. cerana (Ma et al. 2021), mvl expression is related to the behavioral transition from nurse to forager. Our results on mvl expression are only partially consistent with a role in the division of labor of the workers, as statistically significant higher mvl transcript levels were only seen in the abdomen of both species, but not in the head Figure 4A, B. Hence, we consider that the connection between mvl expression in the abdominal tissues and foraging could be related to its transport function for divalent ions, especially iron, since in Drosophila (Wu et al. 2022), malvolio and its mammalian homologs Nramp1 and Nramp2/DMRT1 are important players in the iron metabolism and in infectious diseases (Ganz 2018; Yanatori and Kishi 2019). In the honey bee ectoparasitic mite Varroa destructor, mvl transcripts were found in very different organs, from synganglia to Malpighian tubules, and the levels were higher in reproductive mites than in phoretic ones (Cabrera et al. 2014).

In contrast to mvl, the for gene showed a clearer association with foraging behavior in F. varia workers, both in the head and in the abdomen (Figure 4C), while, surprisingly, in M. quadrifasciata, it is significantly higher expressed in the head of the nurse bees (Figure 4D). The foraging protein and its human homolog PKG1 are cyclic guanosine monophosphate (cGMP)-dependent protein kinases (PKG) related to food choice and foraging behavior in fruit flies (Osborne et al. 1997; Allen and Sokolowski 2021) and also in humans (Struk et al. 2019). In brains of A. mellifera and A. cerana, the for gene is significantly higher expressed in the foragers than in the nurse bees (Ben-Shahar et al. 2002; Heylen et al. 2008; Ma et al. 2018). For bumble bees, there are still conflicting results, as Tobback et al. (2011) found higher for transcript levels in B. terrestris foragers, but this was not the case in the study by Lockett et al. (2016), and in B. ignitus, Kodaira et al. (2009) found higher for expression in nurses, just as we did in M. quadrifasciata. A lack of a consistent link between for expression and foraging behavior was also noted in other social Hymenoptera (see Table 1 in Lockett et al. 2016), and especially so in ants, where it is higher in foragers of Solenopsis (Chen et al. 2022), but not of Pogonomyrmex (Ingram et al. 2005), or related only to age but not behavior, as in Cardiocondyla (Oettler et al. 2015). This apparent inconsistency in the neuroethological role of the foraging gene in the division of labor of social Hymenoptera may be related to the respective sensory modalities associated with the age-related polyethism (Lucas and Ben-Shahar 2021), such as sucrose responsiveness in the honey bee (Thamm and Scheiner 2014), or photoperiod in harvester ants (Ingram et al. 2016) and fire ants (Lei et al. 2019). Furthermore, Thamm et al. (2018) found that the higher for mRNA and PKG protein levels in the honey bee brain are not necessarily reflected in a higher PKG enzymatic activity. For the two stingless bee species, the divergent dynamics of for transcript levels in the brains of M. quadrifasciata and F. varia workers also indicate that there is no simple rule that explains the expression of this gene. This indicates that, in addition to the behavioral state, the degree of worker sterility, as well as age, may also play roles in the expression pattern of these genes.

5 Conclusions

The comparative analysis of candidate genes related to the worker life cycle in two species of stingless bees revealed striking differences among the two species and also in comparison with the honey bee and the bumble bees. Notably, we chose for this study two species occupying the extreme ends of the workers’ reproductive potential. For M. quadrifasciata, the expression patterns for the Vg gene and the JH-related genes mfe and jhe were consistent with their ancestral reproductive function in female insects, similar to what is seen in the bumble bees. In contrast, in F. varia, the expression of these genes and of the nutrient sensing pathway gene tor showed a close association with the behavioral nurse-to-forager transition. But this association apparently involves a configuration of their expression patterns that is just opposite to those known for the honey bee. For the foraging gene, which is commonly related to the neuroethology of foraging behavior in social insects, we found such an association for F. varia, but not for M. quadrifasciata. These findings place a caveat on the notion of a major change in the expression pattern of “universal toolkit genes” that accompanied the transition from a solitary or primitively eusocial to a highly eusocial lifestyle in bees. So far, most studies compare the highly eusocial honey bees with the primitively eusocial bumble bees, or with other social Hymenoptera, hence crossing gaps with respect to levels of sociality, as well as in the hymenopteran phylogenetic tree. Taking advantage of the variability seen in the reproductive biology of the stingless bees (Meliponini), we conclude that, in addition to age, the molecular signatures associated with the division of labor in the workers of these highly eusocial bees are strongly influenced by the respective degree of worker sterility. This should make this species-rich and understudied clade of highly eusocial corbiculate bees an interesting group for insights into life history evolution and associated tradeoffs in a phylogenetically well-defined context.