Patterns of sexual size dimorphism in stingless bees: Testing Rensch’s rule and potential causes in highly eusocial bees (Hymenoptera: Apidae, Meliponini)

Abstract

Eusocial insects offer a unique opportunity to analyze the evolution of body size differences between sexes in relation to social environment. The workers, being sterile females, are not subject to selection for reproductive function providing a natural control for parsing the effects of selection on reproductive function (i.e., sexual and fecundity selection) from other kinds of natural selection. Patterns of sexual size dimorphism (SSD) and testing of Rensch's rule controlling for phylogenetic effects were analyzed in the Meliponini or stingless bees. Theory predicts that queens may exhibit higher selection for fecundity in eusocial taxa, but contrary to this, we found mixed patterns of SSD in Meliponini. Non‐Melipona species generally have a female‐biased SSD, while all analyzed species of Melipona showed a male‐biased SSD, indicating that the direction and magnitude of the selective pressures do not operate in the same way for all members of this taxon. The phylogenetic regressions revealed that the rate of divergence has not differed between the two castes of females and the males, that is, stingless bees do not seem to follow Rensch's rule (a slope >1), adding this highly eusocial taxon to the various solitary insect taxa not conforming with it. Noteworthy, when Melipona was removed from the analysis, the phylogenetic regressions for the thorax width of males on queens had a slope significantly smaller than 1, suggesting that the evolutionary divergence has been larger in queens than males, and could be explained by stronger selection on female fecundity only in non‐Melipona species. Our results in the stingless bees question the classical explanation of female‐biased SSD via fecundity and provide a first evidence of a more complex determination of SSD in highly eusocial species. We suggest that in highly eusocial taxa, additional selection mechanisms, possibly related to individual and colonial interests, could influence the evolution of environmentally determined traits such as body size.

INTRODUCTION

Sexual size dimorphism (SSD) is a widespread phenomenon across the animal and plant kingdoms. It is generally accepted that patterns of SSD result from the selective forces acting on male versus female body size in relation to their different reproductive roles (Andersson, ; Tammaru, Esperk, Ivanov, & Teder, ; Teder, ). That is, selection acts for high fecundity in females, while sexual selection to monopolize females drives body size increase in males (Blanckenhorn, ; Webb & Freckleton, ). However, the ultimate causes of SSD have remained elusive and it has been difficult to determine the relative impact of the interplay between natural and sexual selection in the evolution of SSD. In most invertebrates, females are the larger sex (Stillwell, Blanckenhorn, Teder, Davidowitz, & Fox, ; Teder, ; Teder & Tammaru, ), and in insects, it has been estimated that between 72% and 95% of species within each order exhibit female‐biased SSD suggesting that fecundity selection may be the major driving force acting on SSD in most orders (Stillwell et al., ).

The evolution of SSD has remained understudied in social insects. However, highly eusocial insects are an interesting model to analyze the relative importance of fecundity and sexual selection, because in this group of organisms the females are divided into reproductive queen and (mostly) sterile worker castes (Michener, ; Wilson, ). As workers perform nest duties, queens are presumably released from the constraints of foraging and nest building, and might more readily respond to fecundity selection (Shreeves & Field, ). Although studies of highly eusocial taxa are limited, it is expected that queens should exhibit consistently larger sizes than males (Boomsma, Baer, & Heinze, ). This may be particularly reinforced in highly eusocial bees in which males do not fight to monopolize groups of females and engage in scramble competition for sexual partners (Beani, Dessì‐Fulgheri, Cappa, & Toth, ; Paxton, ).

In many animal taxa, the allometry of body size among sexes varies in accordance with the pattern of SSD, it increases when males are the larger sex, but decreases with body size when females are larger than males (Abouheif & Fairbairn, ; Fairbairn, ; Rensch, ). These trends are explained by greater evolutionary divergence and plasticity in male size, compared with female size, a pattern known as Rensch's rule (Fairbairn, ; Rensch, ). Although there are several explanations for Rensch's rule, this allometric trend is usually attributed to sexual selection acting on male body size (Fairbairn, Blanckenhorn, & Székely, ; Stillwell et al., ). High levels of competition among males are consistent with Rensch's rule (Dale et al., ). The converse trend, where female size varies more than male size, is less common, but seems to be the result of strong fecundity selection acting on females (Blanckenhorn et al., ; Foellmer & Moya‐Larano, ; Webb & Freckleton, ) or on small males that perform aerial displays (Dale et al., ).

Interestingly, in insects, allometric patterns have not always conformed to Rensch's rule. Notably, in solitary Hymenoptera a consistent opposite trend has been found, with variance in female size tending to be greater than variance in male size (Blanckenhorn et al., ; Fairbairn et al., ). Similar to solitary species, a first macroevolutionary study of SSD and Rensch's rule on the primitively social corbiculate bumblebees (Tribe Bombini) revealed a predominantly female‐biased SSD but allometric patterns not conforming to Rensch's rule (Cueva del Castillo & Fairbairn, ). Therefore, fecundity selection seemed stronger compared with sexual selection in corbiculate eusocial bees (Stubblefield & Seger, ; Boomsma et al., ). Surprisingly, in a recent intraspecific study, Medina et al. () found moderately male‐biased SSD in the highly eusocial stingless bee Melipona beecheii (Tribe Meliponini) and the honeybee Apis mellifera (Tribe Apini), and only little difference between males and females in Euglossa mellifera, an eusocial primitive species of the Tribe Euglossini. Thus, contrasting differences in the patterns of SSD have been revealed among the four Tribes comprising the corbiculate bees, possibly related to the level of sociality in this group (Medina et al., ). However, an extensive macroevolutionary analysis in the highly eusocial bees with which to identify proximate and ultimate explanations is lacking.

We used the stingless bees (Tribe Meliponini) as a model to test hypotheses related to body size and allometry in the highly eusocial bees. The Meliponini are a clade of exclusive highly eusocial species living in colonies with morphologically different queen and workers (Sakagami, ). The stingless bees offer the opportunity to better analyze body size patterns because of the high diversity of species (ca. 500 with pantropical distribution; Rasmussen & Cameron, ), in contrast with the only other group of highly eusocial bees, the Apini, in which only ten species have been recognized (Oldroyd & Wongsiri, ). However, in contrast with the Apini, the biology of most species of stingless bees remains unstudied.

We analyzed the evolutionary divergence in body size and sexual size dimorphism among stingless bees compared with primitively social and solitary closer relatives using a series of allometric predictions tested by Cueva del Castillo and Fairbairn () in bumblebees. After controlling by phylogenetic effects, we expected that our results would reflect that selection for female fertility would act equally in all species of stingless bees and that female‐biased SSD would be the norm in this clade. We also expected that selection to act more strongly on queens and males than on workers because workers, being sterile, experience selection only indirectly through their effects on colony success (Kovacs, Hoffman, Marriner, & Goodisman, ; Linksvayer & Wade, ). Thus, a comparison of the evolutionary divergence of queens and males to that of workers should reveal the effects of selection on reproductive function (i.e., fecundity and sexual selection).

If the divergence caused by fecundity selection on queens’ size has been greater than sexual selection on males, the regression of male size on queen size should have a slope below 1, not following Rensch's rule. Otherwise, a slope >1 would indicate that evolutionary divergence caused by sexual selection on males has exceeded fecundity selection on queens. In addition, if the effect of sexual selection on males exceeds that of natural selection on workers, we predict that the regression of male against worker size will have a slope >1. Because it is expected that queens are under strong fecundity selection, we predict that the regression of queen size on worker size selection for foraging and brood care should have a slope >1.

METHODS

For the analyses, we selected an estimator of body size at emergence of the different individuals. For this, we used intertegular distance, a measure currently accepted as standard of body size in different bee taxa, that does not seem to change with age (Bullock, ; Greenleaf, Williams, Winfree, & Kremen, ).

The data on queens, males, and workers used in this study were directly measured on specimens from collections at Universidad Autónoma de Yucatán (J Quezada‐Euán, 12 species), the American Museum of Natural History (Corey Smith, six species), Universidad Nacional de Colombia (Guiomar Nates‐Parra, three species), and the University of Western Sydney (Megan Halcroft, five species). Additionally, we searched published data on size of different species of stingless bees by executing a search on Google Scholar using the terms “body size,” “intertegular width,” and “thorax width.” Google Scholar was used as the search engine because it usually provides full‐text versions of published papers. Data on fifteen more species were obtained from this search. In addition, data from females and males of solitary corbiculate species Euglossa imperialis and Exaerete smaragdina collected in Los Tuxtlas, Veracruz, (Cueva del Castillo, unpublished data) were included in the study.

Phylogenetic reconstruction

To remove the phylogenetic effects from the analysis of data, we built a phylogeny of Meliponini species mainly based on the genetic information available on the GenBank. This genetic information involved partial sequences of two mitochondrial loci (cytochrome oxidase subunit I, COI; and 16S ribosomal RNA, 16S) and four nuclear loci (arginine kinase, ArgK; elongation factor 1 alpha, EF‐1α; long‐wavelength rhodopsin, LWR; and 28S ribosomal RNA, 28S) that were used in previous phylogenetic studies in stingless bees (Françoso & Arias, ; Halcroft et al., ; May‐Itza, Quezada‐Euán, Medina, Enríquez, & De la Rúa, ; Ramírez et al., ; Rasmussen & Cameron, ; Ruiz, May‐Itza, Quezada‐Euán, & De la Rúa, ). Particularly, COI sequences for three Meliponini species (Frieseomelitta nigra, Lestrimelitta niitkib, and Partamona bilineata) were provided by the Canadian Barcode of Life Network. In order to match the species with the data in the molecular sources, we first homologated the taxonomic information found in the literature search with the updated revision by Camargo and Pedro () and Rasmussen and Cameron (). We also retrieved from GenBank genetic information of five corbiculate bee species (Apis dorsata, Bombus terrestris, Euglosa imperialis, Eulaema boliviensis, and Exaerete smaragdina) that were used as outgroups for the phylogenetic analysis as in previous studies. We provide the accession numbers of the GenBank's used sequences on the Supporting Information Table S1.

We aligned the nucleotide sequences of each locus using the algorithm implemented in MUSCLE (Edgar, ). We removed the ending sequence fragments that did not overlap with at least 80% of the species in each aligned dataset. The final number of aligned positions for each dataset was as follows: 1200 for COI, 606 for 16S, 1003 for ArgK, 874 for EF‐1α, 605 for LWR, and 865 for 28S. We then used MESQUITE 3.51 (Maddison & Maddison, ) to concatenate the sequences of each gene and to make the final total evidence matrix for all genes, which comprised 5034 aligned positions and 258 terminals.

We initially partitioned our concatenated matrix in nine partitions considering the coding (COI) and noncoding mitochondrial loci (16S and 28S), as well as the intron and exon regions of the nuclear coding loci (ArgK‐introns, ArgK‐exons, EF‐1α‐introns, EF‐1α‐exons, LWR‐introns, and LWR‐exons). We then estimated the best partition scheme and nucleotide substitution models for each partition using the greedy algorithm implemented in PartitionFinder 1.1.1 (Lanfear, Calcott, Kainer, Mayer, & Stamatakis, ). According to this analysis, we subdivided the final concatenated matrix in seven partitions (P1–P7) that involved the following loci: P1 (16S), P2 (28S and LWR‐exons), P3 (COI), P4 (ArgK‐exons), P5 (ArgK‐introns), P6 (EF‐1α‐exons), and P7 (EF‐1α‐introns and LWR‐introns). The best substitution models estimated for these partitions where as follows: GTR+I+G for P1, P3, and P7; K80+I+G for P2 and P4; and HKY for P5 and P6.

We conducted a concatenated Bayesian inference (BI) analysis in MrBayes 3.2.6 (Ronquist et al., ) with the total evidence matrix by applying the specific substitution model estimated for each partition individually. The BI analysis consisted of four independent runs, each with 20,000,000 generations and four chains, sampling every 2,000 generations. We used default priors for other parameters in the analysis. We assessed parameter convergence and proper mixing of independent runs using TRACER 1.6 (Rambaut, Drummond, Xie, Baele, & Suchard, ). All parameter values sampled during the MCMC of the analysis resulted in ESS values >200. We discarded 25% of the samples obtained prior to stability as burn‐in to obtain a final consensus phylogeny. We used this consensus phylogeny for our posterior comparative analysis considering it was largely congruent with the most robust phylogenies estimated in the group at present (Ramírez et al., ; Rasmussen & Cameron, ).

Comparative analyses

To convert branch lengths of the consensus phylogeny to units of time, we conducted a divergence time analysis in BEAST v1.8.4 (Drummond & Rambaut, ). For this analysis, we used a relaxed clock uncorrelated lognormal method and we assumed a birth–death speciation process for the tree prior to using the topology of our consensus phylogeny. We also applied the same partition scheme and substitution models as in the MrBayes analysis. For the age parameters of the root node (Corbiculates) and the crown group (Meliponini), we assumed normal prior distributions with means and standard deviations as follows: 78 ± 5 Mya and 54 ± 3 Mya, respectively. These two age calibration points were obtained from the most robust divergence study of bees’ lineages at present (Cardinal & Danforth, ). MCMC searches were run for 100,000,000 generations, sampling every 10,000 generations. We performed three independent analysis in CIPRES (Miller, Pfeiffer, & Schwartz, ). We assessed parameter convergence and proper mixing of independent analysis using TRACER 1.6 (Rambaut, Drummond, Xie, Baele, & Suchard, ). All parameter values sampled during the MCMC of the analyses resulted in ESS values >200. We discarded 25% of the samples obtained prior to stability as burn‐in to obtain the chronogram. We used this resulting ultrametric phylogeny to perform all comparative analysis, pruning the species for which we did not obtain morphologic data with the R (version 3.1.3; R Core Team, ) package “ape” (Paradis, Claude, & Strimmer, ). A data file comprising the total evidence matrix of DNA alignments, the consensus BI phylogeny and the chronogram that resulted from the above described analyses is available on the TreeBASE (accession number: 23779).

Ancestral reconstruction of SSD

In order to explore the evolutionary trends in body size and SSD, we used the thorax width of males and queens to build a SSDi index following the Lovich and Gibbons () criteria. This index expresses SSD as [(length of larger sex/length of smaller sex) − 1]. For convention, the SSDi is arbitrarily changed to negative when males are the larger sex and positive when females are the larger sex (Cox, Butler, & John‐Alder, ). Then, we performed an ancestral character reconstruction following Revell () and considering the SSDi values estimated from Meliponini species, as well as from other corbiculate bees species (A. dorsata, B. terrestris, E. imperialis, and E. smaragdina). This method estimates the maximum likelihood value for internal nodes and then interpolates the states along the branches of the phylogenetic tree (see for details Revell, ). For the reconstruction and visualization of ancestral state reconstruction of SDI, we used the R package “Phytools” (Revell, ).

Rensch's rule

Rensch's rule predicts that the slope of a regression of male body size on queen and workers body size will be steeper than 1. To test these predictions in Meliponini species, we used the phylogenetic independent contrasts method (PIC method; Felsenstein, ), as implemented by the package “Caper” (Orme et al., ) in R (ver. 3.0.1; R Development Core Team ) to control for the phylogenetic nonindependence of species (Harvey & Pagel, ). We examined the studentized residuals for outliers >|±3|, but found none in our dataset. Ultimately, we tested the allometric relationship between independent contrasts of log10 male thorax width on log10 queen and worker thorax width, and the allometric relationship between independent contrasts of log10 queen thorax width on the log10 worker thorax width by fitting major axis regression using the R package “smatr” (Warton, Duursma, Falster, & Taskinen, ). Females and males of the other social and solitary corbiculate species were included in the regression of males on queens. Major axis regression offers an accurate approach to test the null hypothesis of isometry (h0: β = 1), because both variables were measured on the same scale and residual variance is minimized in both variables (Warton, Duursma, Falster, & Taskinen, ). Because species of the genus Melipona showed a different trend in the evolution of SSD than the other bees (see below), we performed another phylogenetic independent contrasts analyses excluding this genus in order to evaluate their relative weight in the evolutionary divergence of female and male size in the Tribe.

RESULTS

Ancestral reconstruction of SSD

We found a continuum from a moderate SSD bias to males to a SSD bias to females (Figure ). In 23 species (57%), there is a female‐biased SSD, in six species, males and females showed similar body sizes (Austroplebeia magna, Austroplebeia cassia, Austroplebeia essingtoni, F. nigra, Melipona fallax, and Plebeia frontalis), and 32% (12 species) showed some degree of male‐biased SSD, particularly those of Melipona (Table , Supporting Information Figure S1). Considering that solitary species are predominantly female‐biased SSD, a moderate female‐biased dimorphism or no dimorphism could be the ancestral condition in the corbiculate bees. In our sample, we found that male‐biased SSD possibly evolved independently two times, in Melipona and Trichotrigona extranea (Figure ). Interestingly, the level of SSD seems not related to the body size of the species. Thus, M. beecheii and Melipona fuliginosa show similar levels of SSD bias to males (≈−0.15; Table ), but the mean thorax size of the last species is larger ($\overline{X}$ ≈ 3.63 mm) than M. beecheii ($\overline{X}$ ≈ 2.36 mm) (Table ). Similarly, L. niitkib and Nannotrigona perilampoides have similar levels of SSD bias to females (≈0.19). However, the first species is around 30% larger than the latter. Moreover, species like Trigona amalthea, Trigona chanchamayoensis, and Cephalotrigona capitata with a relatively large body size ($\overline{X}$ ≈ 2.35, 2.38, and 3.19 mm, respectively) show high levels of SSD bias to females (>0.35), than the average for the other non‐Melipona species (Thorax width $\overline{X}$ ≈ 0.18), and the smallest species (Austroplebeia group, P. frontalis, Scaura longula, and Tetragonula iridipennis) show small levels of SSD bias to both females and males (Table ).

Maximum likelihood ancestral reconstruction of SSDi for 37 social and solitary bee taxa performed in R package “phytools” (Revell, ). For the analysis, we used the ultrametric phylogeny and the values of SDI estimated for each species. The values in the color ramp represent the range of SSDi registered for the study species. Negative values indicate male‐biased SSD (blue to gray) and positive values female‐biased SSD (green to violet). SSD: sexual size dimorphism.

Intertegular width values for workers, males, and queens from 37 social and solitary bee species, SSD_i: Lovich and Gibbons Sexual Size Dimorphism index estimated for queen and males (negative values indicate when males are the larger sex), and average thorax width per species (considering just reproductive individuals)

Species	Thorax width (mm)	Species mean thorax size (mm)
Worker	Male	Queen	SSDi
Austroplebeia australis 4	1.31	1.31	1.43	0.091	1.37
Austroplebeia cassiae 4	1.36	1.37	1.45	0.058	1.41
Austroplebeia cincta 4	1.11	1.15	1.23	0.070	1.19
Austroplebeia essingtoni 4	1.10	1.13	1.17	0.035	1.15
Austroplebeia magna 4	1.33	1.35	1.39	0.030	1.37
Axestotrigona ferruginea 3	1.925	2.00	2.32	0.160	2.16
Bombus terrestris 11	4.77	7.76	7.86	0.013	7.81
Cephalotrigona capitata 2	2.88	2.63	3.75	0.425	3.19
Cephalotrigona zexmeniae 13	2.081	1.925	2.468	0.282	2.196
Euglossa imperialis 14		5.38	5.61	0.043	5.50
Exaerete smaragdina 14		7.05	8.11	0.150	7.58
Frieseomelitta nigra 13	1.227	1.307	1.352	0.034	1.329
Lepidotrigona ventralis 6	1.44	1.48	2.05	0.385	1.765
Lestrimelitta limao 2	2.00	2.00	2.75	0.375	2.375
Lestrimelitta niitkib 13	1.493	1.584	1.89	0.191	1.737
Melipona beecheii 12	2.528	2.556	2.158	−0.155	2.357
Melipona bicolor 2	3.875	4.00	3.25	−0.188	3.625
Melipona colimana 13	2.894	2.985	2.699	−0.095	2.842
Melipona eburnea 5	4.89	4.708	4.409	−0.063	4.558
Melipona fallax 8	3.43	3.53	3.48	−0.014	3.51
Melipona favosa 9	2.633	2.675	1.968	−0.264	2.321
Melipona fuliginosa 8	3.68	3.92	3.33	−0.151	3.625
Melipona marginata 9	2.299	2.307	1.86	−0.193	2.083
Melipona quadrifasciata 9	2.888	2.884	2.61	−0.095	2.747
Melipona rufiventris 1	3.75	4.00	3.25	−0.188	3.63
Melipona scutellaris 10	3.01	2.87	2.71	−0.056	2.79
Melipona yucatanica 13	2.235	2.245	1.97	−0.11	2.107
Nannotrigona perilampoides 13	1.309	1.224	1.453	0.187	1.338
Oxytrigona tataira 9	1.526	1.522	2.16	0.424	1.841
Partamona bilineata 13	1.589	1.561	1.655	0.060	1.608
Plebeia frontalis 13	1.02	1.018	0.993	−0.024	1.00
Scaptotrigona mexicana 13	1.453	1.524	1.82	0.194	1.672
Scaptotrigona pectoralis 13	1.613	1.582	1.831	0.157	1.706
Scaura latitarsis 2	1.25	1.25	1.50	0.200	1.375
Scaura longula 2	1.50	1.50	1.66	0.107	1.58
Schwarziana quadripunctata 9	1.748	1.865	2.039	0.093	1.952
Tetragonisca angustula 5	1.019	1.061	1.886	0.777	1.473
Trichotrigona extranea 7	NA	1.48	1.25	−0.155	1.365
Trigona amalthea 5	2.218	1.982	2.729	0.376	2.355
Trigona chanchamayoensis 2	1.88	2.00	2.75	0.375	2.375
Trigona corvina 2	2.13	2.25	2.75	0.222	2.50
Trigona fulviventris 13	1.429	1.606	1.871	0.165	1.738
Trigona nigerrima 2	2.38	2.50	2.75	0.100	2.625
Trigona pallens 9	1.416	1.722	2.047	0.188	1.884

¹Schwarz ().

²Schwarz ().

³Darchen ().

⁴Halcroft (unpublished data).

⁵Nates‐Parra (unpublished data).

⁶Sung, Yamane, Ho, Wu, and Chen ().

⁷Camargo and Pedro ().

⁸Camargo and Pedro ().

⁹Smith (unpublished data).

¹⁰de Araujo‐Alves ().

¹¹Cueva del Castillo, Sanabria‐Urbán, and Serrano‐Meneses ().

¹²Medina et al. ().

¹³Quezada‐Euán (unpublished data).

¹⁴Cueva del Castillo (unpublished data).

Rensch's rule

After controlling for phylogenetic nonindependence among of the species studied, the results of the major axis regression of independent contrasts indicated strong coevolution between queens and males (r² = 0.82; df = 42, p < 0.0001), between queens and workers (r² = 0.86; df = 39, p < 0.0001), and between males and workers (r² = 0.93; df = 39, p < 0.0001; Figure ). However, none of the regressions showed a slope significantly different from 1 (males on workers: β = 1.01, p = 0.89; queens on workers: β = 1.12, p = 0.09; and males on queens: β = 0.89, p = 0.12). Thus, it seems that evolutionary divergence has been similar in males and the female sexual and sterile castes. Nonetheless, after removing Melipona species from the analyses, we found that the slope of males on queens was significant and smaller than 1 (β = 0.84, p = 0.04; Figure ) and the slope of the regressions of queen on workers (β = 1.19, p = 0.03) was significant and stepper than 1, suggesting that queen size has diverged more than male and worker size in non‐Melipona species, even though males on workers (β = 1.01, p = 0.91) remain no significant.

Allometric major axis regressions of independent contrasts (IC) of (a) males’ on queens’, (b) males’ on workers’, and (c) queens’ on workers’ thorax width of stingless bees. Dashed lines indicate isometry (β = 1)

Allometric major axis regressions of independent contrasts (IC) after removing Melipona species of (a) males’ on queens’, (b) males’ on workers’, and (c) queens’ on workers’ thorax width of stingless bees. Dashed lines indicate isometry (β = 1)

DISCUSSION

Contrary to expectations, our results quantitatively confirm the existence of two different patterns of SSD in stingless bees, revealing the Meliponini as a taxon with mixed SSD. Nevertheless, compared to solitary species moderate levels of SSD were found. In addition, no evidence of hyperallometry was evident when the size of males was regressed on queens, making Meliponini a first clade of highly eusocial insects in which it is shown that Rensch's rule does not hold (Abouheif & Fairbairn, ; Blanckenhorn et al., ; Cueva del Castillo & Fairbairn, ; Fairbairn, ; Webb & Freckleton, ). Perhaps due to the independent evolution of differences in size and SSD bias to females and males (Figures and ), we did not detect significant departures from isometry for our phylogenetically controlled regressions for the thorax width of males on queens, males on workers, and queens on workers. Nonetheless, when we removed Melipona from the analysis, the phylogenetically controlled regressions for the thorax width of males on queens have a slope significantly smaller than 1, whereas the regression of queens on workers has a slope significantly stepper than 1, suggesting that similar to other eusocial species the magnitude of evolutionary divergence has been larger in queen than males, and could be explained by stronger selection on female fecundity. This strong selection of queen's fecundity also can explain the size divergence between queens and the sterile workers in non‐Melipona species (Fairbairn et al., ).

Apparently, independent evolution of differences in size and SSD bias to queens and males suggests that the interplay between natural and sexual selection has operated in different ways in Meliponini, resulting in two patterns of SSD, male‐ and female‐biased, which also suggest that as in other animals, the direction and magnitude of the selective pressures do not operate in the same way for all members of this taxon (Colwell, ; Kratochvíl & Frynta, ; Serrano‐Meneses, Cordoba‐Aguilar, Azpilicueta‐Amorín, González‐Soriano, & Székely, ). The reproductive behavior of stingless bees remains largely unknown (Quezada‐Euán, ). The lack of sufficient empirical data prevents an accurate explanation of the evolutionary meaning and reproductive implications of the opposite patterns of SSD found in this taxon. One possible explanation for the male‐biased SSD found in Melipona species and T. extranea could be sexual selection acting on males’ body size. Even though there are few studies of sexual selection in social insects (Berg, Koeniger, Koeniger, & Fuchs, ; Boomsma & Ratnieks, ), there is no reason to discard its importance in the evolution of male traits (Baer, ; Boomsma et al., ) and SSD bias to males in social insects, especially when it is well known that males from eusocial bees usually engage in scramble competition during mating events (Paxton, ). Nevertheless, differences in SSD in stingless bees could also be related to selection on the size of queens rather than males. Indeed, it has been suggested that reproductive females when released from the burden imposed by the collection of food and construction of nests can experience more flexibility in their development (Medina et al., ). Cuckoo bees, for instance, show a tendency toward less female‐biased SSD compared to nest‐building species (Shreeves & Field, ). Likewise, in M. beecheii, Medina et al. () found that weight at emergence varied widely among individual queens (up to fourfold) compared with males (less than twofold).

A proximate cause to explain size differences among sexes in solitary insects is the difference in development time of males and females (Teder, ). Differences in development time of queens relative to males could also extend to explain differences among genera in Meliponini. Empirical data indicate that in Melipona males, workers and queens receive similar amounts of food as larvae and pass through the same stages of development, but the pupae stage of queens is shorter than in males and workers (Jarau et al., ; Macías‐Macías & Quezada‐Euán, ; Moo‐Valle, Quezada‐Euán, Canto‐Martín … González‐Acereto, ; Quezada‐Euán, ; Quezada‐Euán et al., ). Such accelerated rate of development can explain why the queens in Melipona species are comparatively smaller than the queens of non‐Melipona, in which queens take longer to develop than males and workers, probably as a consequence of the excess food they receive as larvae. In support for this argument, in non‐Melipona species, female larvae receiving less food can also become queens, but of minute size which also have reduced development times (Morales, ). In Melipona, female larvae become queens not only by food but also by a possible genetic mechanism (Jarau et al., ; Kerr, ).

Interestingly, Melipona colonies produce excess queens that are eliminated by the workers (Wenseleers, Hart, Ratnieks, & Quezada‐Euán, ). Such apparently unnecessary production of surplus queens in Melipona has been explained as result from the selfish behavior of female larvae attempting to evade a worker fate with less reproductive advantages (Ratnieks, ; Ratnieks & Wenseleers, ). Perhaps, a fast rate of development (and the resulting small size) is a result of the competition among the many queens that are produced at one time in Melipona, but not in the other genera of stingless bees, which produce substantially less queens. Queens relieved from the burden of nest construction, and provision can trade traits requiring large body size for faster development that would result in smaller size. However, this would be adaptive only if queen fitness may not suffer as a consequence of reduced size up to one point. Indirect empirical evidence from non‐Melipona species suggests that minute queens (presumably experiencing fast development rates) are capable of heading colonies, which seems to confirm this hypothesis (Imperatriz‐Fonseca & Zucchi, ; Ribeiro, Wenseleers, Santos Filho, & Alves, ).

Our results in the stingless bees provide a first evidence of additional evolutionary forces that may influence SSD in highly eusocial species. We suggest that in highly eusocial taxa, a trade‐off between individual and colonial interests, in addition to sexual and fecundity selection, could influence the evolution of environmentally determined traits such as body size. The collection of more empirical data on the relative size of males and females in highly eusocial insects is necessary to test such hypotheses. Studies comparing body size variation of males and females within and across regions are also needed to analyze environmental effects that may determine male and female fitness in highly eusocial insects.

ACKNOWLEDGMENTS

The authors wish to thank Victor Hugo Jimenez Arcos for his valuable help in performing comparative analyses. We are most thankful to colleagues who kindly measured specimens in their collections and shared the information with us: Guiomar Nates‐Parra (Universidad Nacional de Colombia), Klaus Hartfelder (Universidade de São Paulo), Megan Halcroft (University of Western Sydney), Octavio Macías‐Macías (Universidad de Guadalajara), and Miguel Guzman (Ecosur‐Tapachula). We thank Laurence Packer, Ricardo Ayala, and Eunice Enriquez, who made available DNA sequences form the the Canadian Barcode of Life Network (Genome Canada, the Natural Sciences and Engineering Research Council of Canada, and sponsors listed at www.BOLNET.ca). We also thank Salvador Medina Peralta (UADY) and Ismael Hinojosa (UNAM) for their help in data acquisition and organization. This study was partially supported with funds from the project CONACyT 103341 “Conservación de las abejas sin aguijón de México (Hymenoptera: Meliponini)” and 291333 “Manejo sustentable de polinizadores.”

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTIONS

J. Javier G. Quezada‐Euán collected and provided part of the data required for this work and contributed in the process to write the paper. Corey Smith collected and provided part of the data required for this paper. Salomón Sanabria‐Urbán performed the phylogenetic reconstructions, performed part of the comparative analyses, and contributed in the process to write the paper. Raul Cueva del Castillo performed part of the comparative analyses and contributed to writing of the paper.

DATA ACCESSIBILITY

The genetic and morphological information used in this study was mainly obtained from previously published studies which are cited throughout the article. The GenBank accession numbers and source references of all the used gene sequences are provided in the Supplementary Information Table S1. The total evidence datasets used for our analyses, including the obtained DNA alignments, consensus BI phylogeny and Chronogram, are also available on the TreeBASE with the following accession number: 23779.