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.


| INTRODUC TI ON
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, 1994;Tammaru, Esperk, Ivanov, & Teder, 2010;Teder, 2014). That is, selection acts for high fecundity in females, while sexual selection to monopolize females drives body size increase in males (Blanckenhorn, 2000;Webb & Freckleton, 2007). 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, 2010;Teder, 2014;Teder & Tammaru, 2005), 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., 2010).
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, 1974;Wilson, 1971). 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, 2008). Although studies of highly eusocial taxa are limited, it is expected that queens should exhibit consistently larger sizes than males (Boomsma, Baer, & Heinze, 2005). 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, 2014;Paxton, 2005).
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, 1997;Rensch, 1950). These trends are explained by greater evolutionary divergence and plasticity in male size, compared with female size, a pattern known as Rensch's rule Rensch, 1950).
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, 2007;Stillwell et al., 2010). High levels of competition among males are consistent with Rensch's rule (Dale et al., 2007). 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 Foellmer & Moya-Larano, 2007;Webb & Freckleton, 2007) or on small males that perform aerial displays (Dale et al., 2007).
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 Fairbairn et al., 2007). 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, 2012). Therefore, fecundity selection seemed stronger compared with sexual selection in corbiculate eusocial bees (Stubblefield & Seger, 1994;Boomsma et al., 2005). Surprisingly, in a recent intraspecific study, Medina et al. (2016) 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., 2016). 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, 1982). 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, 2010), in contrast with the only other group of highly eusocial bees, the Apini, in which only ten species have been recognized (Oldroyd & Wongsiri, 2006). 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 (2012) 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, 2010;Linksvayer & Wade, 2009). 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.

| ME THODS
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, 1999;Greenleaf, Williams, Winfree, & Kremen, 2007).
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 &  we first homologated the taxonomic information found in the literature search with the updated revision by Camargo and Pedro (2007) and Rasmussen and Cameron (2010). 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, 2004). 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, 2015) 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 conducted a concatenated Bayesian inference (BI) analysis in MrBayes 3.2.6 (Ronquist et al., 2012) 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, 2018). 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., 2010;Rasmussen & Cameron, 2010).

| 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, 2007). 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, 2013).
MCMC searches were run for 100,000,000 generations, sampling every 10,000 generations. We performed three independent analysis in CIPRES (Miller, Pfeiffer, & Schwartz, 2010). We assessed parameter convergence and proper mixing of independent analysis using TRACER 1.6 (Rambaut, Drummond, Xie, Baele, & Suchard, 2018). 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, 2015) package "ape" (Paradis, Claude, & Strimmer, 2004). 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 (1992)  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, 2013Revell, , 2014. For the reconstruction and visualization of ancestral state reconstruction of SDI, we used the R package "Phytools" (Revell, 2012).

| 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, 1985), as implemented by the package "Caper" (Orme et al., 2013) in R (ver. 3.0.1; R Development Core Team 2013) to control for the phylogenetic nonindependence of species (Harvey & Pagel, 1991). 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, 2012). 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, 2012). 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.

| Ancestral reconstruction of SSD
We found a continuum from a moderate SSD bias to males to a SSD bias to females (Figure 1). In 23 species (57%), there is a female-biased SSD, in six species, males and females showed similar body sizes (Austroplebeia  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 1).
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 1), but the mean thorax size of the last species is larger (X ≈ 3.63 mm) than M. beecheii (X ≈ 2.36 mm) (Table 1)

| 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 2 = 0.82; df = 42, p < 0.0001), between queens and workers (r 2 = 0.86; df = 39, p < 0.0001), and between males and workers (r 2 = 0.93; df = 39, p < 0.0001; Figure 2). 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 3) 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.

| D ISCUSS I ON
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, 1997;Blanckenhorn et al., 2007;Cueva del Castillo & Fairbairn, 2012;Webb & Freckleton, 2007). Perhaps due to the independent evolution of differences in size and SSD bias to females and males (Figures 2 and 3), 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 F I G U R E 1 Maximum likelihood ancestral reconstruction of SSDi for 37 social and solitary bee taxa performed in R package "phytools" (Revell, 2012). 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. fecundity also can explain the size divergence between queens and the sterile workers in non-Melipona species .
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, 2000;Kratochvíl & Frynta, 2007 (Berg, Koeniger, Koeniger, & Fuchs, 1997;Boomsma & Ratnieks, 1996), there is no reason to discard its importance in the evolution of male traits (Baer, 2005;Boomsma et al., 2005) 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, 2005). 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., 2016). Cuckoo bees, for instance, show a tendency toward less female-biased SSD compared to nest-building species (Shreeves & Field, 2008). Likewise, in M. beecheii, Medina et al. (2016) 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, 2014). 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, 2018). In Melipona, female larvae become queens not only by food but also by a possible genetic mechanism (Jarau et al., 2010;Kerr, 1948 Schwarz (1932). b Schwarz (1948). c Darchen (1971). d Halcroft (unpublished data). e Nates-Parra (unpublished data). f Sung, Yamane, Ho, Wu, and Chen (2004). g Camargo and Pedro (2007). h Camargo and Pedro (2008). i Smith (unpublished data). j de Araujo-Alves (2010). k Cueva del Castillo, Sanabria-Urbán, and Serrano-Meneses (2015). l Medina et al. (2016). m Quezada-Euán (unpublished data). n Cueva del Castillo (unpublished data).
TA B L E 1 (Continued) 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.

ACK N OWLED G M ENTS
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 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 F I G U R E 2 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) F I G U R E 3 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) México (Hymenoptera: Meliponini)" and 291333 "Manejo sustentable de polinizadores."

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
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 ACCE SS I B I LIT Y
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