Introduction

Polyandry is the term applied to mating systems in which an individual female mates with more than one male (Pizzari and Wedell 2013). Once considered a rarity, polyandry is now recognized as a feature common across sexually reproducing organisms ranging from angiosperms (Pannell and Labouche 2013) to insects (Arnqvist and Nilsson 2000) and to primates (Muller and Thompson 2012). Early polyandry research focused on its implications for sexual selection in individuals, with an emphasis on sperm competition (Parker 1970; Simmons 2002), cryptic female choice (Childress and Hartl 1972; Firman et al. 2017), and diploid male load (Page 1980). Today polyandry is appreciated as a dynamic in processes as fundamental as population extinctions, social networks, altruism and cooperation (Pizzari and Wedell 2013), and group formation and complexification across evolutionary transitions (Bourke 2011).

Polyandry played a central role in the evolution of complex sociality in insects. Evidence for this is suggested in the derived loss of the spermatheca in sterile castes with its simultaneous retention as a fully functional and capacious organ in the queens (Treanore et al. 2020). Thus equipped to store sperm from her multiple mates, nourish it and use it across her lifetime, a polyandrous queen is able to produce a colony of genetically diverse daughters, a condition often (Bourke 2011) though not necessarily (Linksvayer and Johnson 2019) associated with caste polymorphisms and increasing colony size and ecologic success.

Although a queen’s promiscuity reduces average relatedness among nestmates, appearing to threaten the bases of eusocial cooperation (Ratnieks and Helanterä 2009), the increase in within-nest genetic diversity increases colony fitness, as shown from examples across the ants, wasps, and bees (Keller and Reeve 1994; Schmid-Hempel 1998; Hughes and Boomsma 2004; Goodisman et al. 2007; Seeley and Tarpy 2007). Polyandry is thus integral to the transition from simple to complex eusociality (Hughes et al. 2008; Bourke 2011) and by extension to social group formation once response to between-group selection exceeds response to within-group selection (West et al. 2007).

Among the mechanistic explanations for the evolution of polyandry is a group known as the genetic variance (GV) hypotheses (Keller and Reeve 1994). Each is a variant of the idea that within-colony genetic diversity improves colony fitness through better social integration, group task efficiency, and ecologic competitiveness (Gove et al. 2009). The GV hypotheses are reviewed and appraised by Crozier and Fjerdingstad (2001).

The GV hypotheses are consistent, but only up to a point, with computed rates of in-nest genetic diversity as a function of effective mating number (me = the queen’s mating number adjusted for relative contribution of each of her mates to actual progeny (Starr 1985)). The relationship between me and a colony’s genetic diversity is not linear, so that over 90% of the increase in a colony’s diversity potential is achieved by its mother’s 10th mating (Fig. 1) (Palmer and Oldroyd 2000). In other words, any hypothesis for polyandry based on genetic diversity per se, runs out of explanatory power after the tenth mating. This is a problem when we consider that some females mate at rates far in excess of the genetic diversity (GD) asymptote, rates hereafter called hyperpolyandry, as shown in the African army ant Dorylus (Anomma) molestus me = 17.5 (Kronauer et al. 2004), neotropical army ants me ≤ 25 (Barth et al. 2014), and honey bees Apini me ≤ 77 (Neumann and Moritz 2000; Tarpy et al. 2004; Withrow and Tarpy 2018). How does natural selection produce such high rates of promiscuity when marginal increases in nestmate heterogeneity—the presumed effect driver—seem so small?

Fig. 1
figure 1

Average hymenopteran within-colony genetic diversity as a function of the queen’s effective mating number, me (Palmer and Oldroyd 2000)

Full size image

One answer could be that hyperpolyandry is an adaptation by queens for avoiding sperm depletion (Boomsma and Ratnieks 1996; Fjerdingstad and Boomsma 1998; Kraus et al. 2004).

Another answer could be that fitness-improving synergisms among workers of different patrilines continue increasing at mating frequencies beyond the GD asymptote. Honey bees are well-known for their age-based polyethism schedule by which an individual progresses from interior brood- and nest-oriented tasks toward exterior guarding and foraging tasks (see Table 6.1 in Winston (1987)). Subfamily-based proclivities are known for several of these life-critical common tasks including guarding, undertaking, pollen collecting, nectar foraging, scouting, and nursing larvae (Calderone and Page 1991; Chapman et al. 2007). Optimizing ratios of such subfamily-based proclivities may be one adaptive outcome of hyperpolyandry according to a model by Fuchs and Moritz (1999). And while average within-nest diversity as a function of mating number is capped by rules of hymenopteran inheritance (Fig. 1), the number of genotypic compositions has no such ceiling. Optimal compositions may continue increasing at mating frequencies in excess of me = 10, and indeed it is in this zone that the model of Fuchs and Moritz predicts the highest fitness returns. These authors envisioned these optimal compositions as ratios of specialists to non-specialists, but we propose that the mechanism can apply more generally to the simple variance that exists in task proclivities among patrilines with the result being more efficient execution of common tasks.

A second aspect of the model of Fuchs and Moritz (1999) applies to the inheritance of phenotypes that are rare and beneficial. Rarity is selected for because a colony dominated by specialists will be short on workers who perform life-critical common tasks. In the case of uncommon tasks such as behavioral parasite resistance, rarity can also be an outcome of selection that balances costs of maintaining the character against benefits realized by sustaining it (Boots and Begon 1993; van Baalen 1998). As a queen enlarges her sampling of local drones, she more closely matches her colony’s allele frequencies to those of the breeding population. This averaging effect reduces the likelihood of environmental failures from phenotypic anomalies (Page et al. 1995a; Rueppell et al. 2008) and applies notably to rare beneficial alleles vulnerable to cliff edge effects (Mountford 1968; Boyce and Perrins 1987). A colony “falls off the cliff” when resistance specialisms are so rare that the colony fails to procure any specialists at all (Fuchs and Moritz 1999). Hyperpolyandry can be seen as a bid for avoiding such catastrophes.

It seems to us, therefore, that two testable hypotheses emerge from this discussion. (1) For colony phenotypes associated with common tasks, colony fitness will continue increasing with mating levels in excess of the GD asymptote, consistent with the idea of hyperpolyandry increasing optimum genotype compositions. (2) For colony phenotypes associated with uncommon tasks, benefits of the rare specialism will be particular to conditions of hyperpolyandry, consistent with the idea of hyperpolyandry capturing rare allele combinations, resisting cliff-edge effects.

In this study, we used instrumental insemination to create a 2 × 2 factorial experiment looking at two levels of polyandry crossing with two genetic lines of drone semen. The two levels of polyandry bracket the ranges of interest, with one near the GD asymptote (observed mating number = mo = 9) and the other well into the range of hyperpolyandry (mo = 54). By crossing the experiment with two genetic lines, we insert allelic diversity into the design, letting us observe outcomes of interactions between polyandry and common or rare traits. Drones of the first genetic line were sons of unselected wild type queens from whom we expected the full range of common colony phenotypes. Drones of the second line were sons of queens whose mothers were instrumentally inseminated breeders from a program selecting for Varroa Sensitive Hygiene (VSH), a rare behavioral resistance trait, and propagated by Jeffrey Harris at Mississippi State University. Genetic proclivities at the subfamily or individual level are known for the common tasks brood nursing (Chapman et al. 2007) as well as for VSH (Spivak and Danka 2020).

VSH is a subset of the more general state of hygienic behavior, a heritable defensive response of honey bees to pathogens and parasites in brood. Worker bees expressing this character detect and uncap brood cells containing pathogens or parasites and remove the contents, interrupting the parasite’s reproductive cycle. Originally a natural adaptation associated with the bacterial brood pathogen Paenibacillus larvae (Woodrow and Holst 1942; Rothenbuhler 1964), hygienic behavior has since been associated with brood invaded by the parasitic brood mite Varroa destructor (Harbo and Harris 2005; Ibrahim and Spivak 2006).

It was Rothenbuhler (1964) who originally ascribed hygienic behavior to two independent loci—brood cell uncapping and removing—which when occurring at high rates of recessive homozygosity result in the hygienic colony phenotype. Hygienic behavior directed at V. destructor (VSH) is known to be uncommon, with 768 of 22,000 (3.49%) (Spötter et al. 2016) to 270 of 1000 (27%) (Scannapieco et al. 2017) of phenotyped worker bees expressing some components of the syndrome. The genes controlling the behavior number at least 4–5 and are generally associated with associative learning, sensitivity to external stimuli, memory, and olfaction (Spötter et al. 2016; Scannapieco et al. 2017). Although VSH is not the only resistance character in A. mellifera against Varroa mites, colonies of A. mellifera lacking any resistance phenotypes invariably die without artificial interventions (Büchler et al. 2010; Mondet et al. 2020). The combination of control by few loci of moderate to strong effect, high fitness returns, and vulnerability to cliff edge catastrophes makes VSH a good candidate for the rare specialist model for hyperpolyandry evolution.

There is, however, one criterion under which hygienic behavior does not meet expectations for a driver of polyandry. Uncommon resistance tasks are predicted to exert positive selection on polyandry (and by extension hyperpolyandry) when they are expressed by behavioral dominance, i.e., when only a small number of actors is needed to confer resistance to the whole colony (Sherman et al. 1988; Pamilo et al. 1994). This does not describe the hygienic behavior directed at brood disorders. Using freeze-killed brood as a stimulus in colonies composed of 0, 25, 50, or 100% hygienic workers, Arathi and Spivak (2001) showed that removal of freeze-killed brood increased linearly with hygienic compositions, achieving 100% only in colonies with 100% hygienic workers. Ecologically speaking, therefore, hygienic behavior operates not by “few actors for big effect,” but rather by “many actors for any effect.” To the extent these dynamics describe VSH, this suggests a non-mutually excluding variant of our hypothesis (2): hyperpolyandry avoids cliff edge catastrophes when the environmental stressor demands high colony expression frequencies of a rare task genotype.

Methods

Queens and inseminations

We performed an experiment near Athens, Georgia consisting of four treatments in a 2 × 2 factorial arrangement with two levels of genetic Varroa resistance (semen from VSH-selected drones or non-selected wild type drones) and two levels of polyandry (high or low). There was one dedicated research apiary in which all treatments were represented. Each colony consisted of a 5-frame Langstroth nucleus hive and was started with one 2-lb (0.9-kg) package of worker bees from one supplier practicing a common V. destructor treatment regimen, thus normalizing incipient parasite levels across the experiment. Colonies were arranged in one large circle at the apiary to minimize effects of bee and parasite drift (Dynes et al. 2019a). The experiment was set up in April 2018 and ran until November 2018. The year before, we made up two sets of dedicated drone-source colonies for future inseminations, one set headed by wild type unselected queens and the other headed by queens from the VSH line.

In spring of 2018, we reared one wild type virgin queen and inseminated her with one wild type male. This single-drone inseminated queen was used to rear supersister virgin daughters who were subsequently inseminated and used in the experiment. This procedure minimizes maternal variation in experimental colonies (Harbo 1985). Supersister virgins were each randomly assigned one of four treatments in the 2 × 2 design: (1) VSH semen/low polyandry, (2) VSH semen/high polyandry, (3) wild type semen/low polyandry, or (4) wild type semen/high polyandry.

We had nine strong VSH drone-source colonies come out of winter, a number that dictated the multiples for our polyandry ranges. For our low polyandry treatment, we inseminated each supersister virgin with one male from each of the VSH drone source colonies (mo = 9), and for our high polyandry treatment, each queen was inseminated with 6 males from each VSH drone source colony (mo = 54). Supersister virgins in the wild type treatments similarly received one drone from each of 9 wild type drone source colonies (mo = 9) or 6 from each wild type source colony (mo = 54).

Each virgin queen was emerged and housed in her 5-frame nucleus colony before and after the artificial insemination (AI) procedure, following methods of Cobey et al. (2013). Semen from 9 or 54 males, as per protocol, was collected into one common capillary tube. A different tube was used for each prescribed polyandry/genetic line; no contamination across treatments was possible. However, more than one queen was inseminated from any given tube. The volume of semen in each tube was increased by 15% with physiological saline (recipe 2.2.1 in Cobey et al. (2013)), thoroughly mixed in an Eppendorf tube to encourage a homogeneous drone representation per batch, re-drawn into the capillary tube, then used to inseminate 2–20 queens at an average dose of 4-μL mixed semen each, thus removing effects of semen volume (Niño et al. 2012).

Queens were narcotized three times with CO2—once on the day of insemination and once on each of two successive days thereafter; CO2 narcosis stimulates egg laying (Mackensen and Roberts 1948). Experimental queens were housed in their home 5-frame nuclei and observed until all competent egg layers were identified. We waited 6 weeks thereafter to allow worker populations to turn over to progeny of the experimental queens.

Dependent variables

Colonies were sampled at each of four time points (July, August, September, November). We sampled colonies to determine brood production, brood survival, Varroa mite levels, and adult bee populations.

Brood production and worker bee population were derived by visually summing proportions of whole deep frames covered by brood or workers (Delaplane et al. 2013), converting frames of brood to cm2 by the observation that one deep Langstroth comb (both sides) = 1760 cm2, and converting frames of adult bees to bee populations with the regression model of Burgett and Burikam (1985). When necessary, we converted cm2 brood to cells of brood with the conversion of 3.8 cells per cm2 (Delaplane et al. 2013).

Brood survival was measured by removing a comb of open brood, laying over it a sheet of transparent acetate, using a felt-tip permanent marker to mark the location of 100 brood cells each containing a 1st or 2nd instar larva, then returning the comb to the hive; 3 days later, the comb was retrieved; the same acetate laid on top of it, and the number of surviving brood cells recorded.

Relative Varroa mite numbers are derived by inserting sticky sampling sheets into bottom board hive inserts and recording the number of Varroa mites trapped after 24 h (Dietemann et al. 2013).

Patriline determination

As a check on our success at creating two discrete classes of polyandry with AI, we were able to genotype workers from 18 colonies to determine effective realized paternity (me) in both polyandry classes.

A 50-worker sample was taken from each of 8 colonies in the mo = 9 polyandry class and 10 colonies in the mo = 54 polyandry class. DNA was extracted from the right hind leg for all honey bee specimens using a Qiagen DNeasy Blood & Tissue Kit. Hind legs were thoroughly pulverized with micro-scissors and digested with proteinase K in a 56 °C water bath for 3 h with sample agitation every hour. Extracted DNA was stored at − 80 °C until microsatellite amplification. Ten variable microsatellite loci were screened for all samples using two multiplexes: Plex 1 and Plex 2. Plex 1 included loci A107, A113, AP043, A024, and A006; Plex 2 included loci A28, A88, AP66, AP81, and B124 (Shaibi et al. 2008; Delaney et al. 2009). Both plexes were amplified for one cycle at 95 °C for 7 min, 30 cycles of 95 °C for 30 s, 54 °C for 30 s, 72 °C for 30 s, and a final extension at 72 °C for 60 min. A 10-μL final reaction volume containing 5 μL of PCR Master Mix (Promega, Madison, WI), 1.0–2.5 μL of fluorescent dye-labeled primer, 0.9 μL of nuclease-free water, and 2 μL of DNA extract per sample was analyzed at the University of Delaware Sequencing and Genotyping Center at the Delaware Biotechnology Institute with an Applied Biosystems 3130 XL Genetic Analyzer via capillary electrophoresis.

Microsatellite repeat sizes were scored using Geneious Prime® 2020.0.5 software (Biomatters Ltd). The number (N0) and frequency (pi) of full sibling patrilines were determined using raw microsatellite data in the software program COLONY 2.0.6.5 (Jones and Wang 2010). The observed patriline number (N0), proportion of each patriline within each sample (pi), and number of worker bees in each respective sample (n) were used to calculate me after (Nielsen et al. 2003; Tarpy et al. 2015):

$${m}_{e}=\frac{({n-1)}^{2}}{\sum_{i=1}^{{N}_{0}}{p}_{i}^{2}\left(n+1\right)\left(n-2\right)+3-n}$$

Statistical analyses

Data were analyzed with the JMP Pro platform 15.0 and SAS Studio software, Version 3.8 for Linux (© 2018 SAS Institute). To improve compliance with ANOVA assumptions for normality, data for cm2 brood and bee population were square root transformed, and data for mite counts were transformed to the natural log. Mixed model ANOVAs were used, recognizing as fixed effects polyandry level (high or low), genetics (VSH or wild), sampling month, and all interactions thereof. Surviving colonies were sampled at each of four time points: July, August, September, and November. Colonies were thus treated as random repeating measures; their auto-correlations across time controlled with the compound symmetry covariance structure method. The exception was for brood survival which was measured only in July; for this, we ran a Wilcoxon non-parametric rank sums analysis, comparing brood survival by polyandry and genetics. Values in figures are presented as natural numbers.

Colony strength data such as ours are often autocorrelated, especially in repeated measures studies (Dynes et al. 2019b). Of special concern to us was the possibility that colony mite levels were independently affected by bee colony demographics. Increasing ratios of brood/adult bee are associated with increasing mite reproduction and population growth (Boot et al. 1994; Martin and Kemp 1997). So, for mite levels we included cells of brood per bee as a fixed effect covariate.

Results

Effective realized paternity (me) was significantly different (F = 34.7; df = 1,16; P = 0.0001; one-way ANOVA) between polyandry mo = 9 (9.1 ± 2.2; lsmeans ± SE; range = 6.2–15.1) and mo = 54 (26.5 ± 2.0; range = 12.9–37.2), confirming that our AI protocols produced two discontinuous classes of polyandry.

Table 1 gives ANOVA terms for three dependent variables. For cm2 brood, there were significant effects of month with a predictable seasonal decline in November (Fig. 2).

Table 1 ANOVAs of main effects genetics (VSH or wild), level of polyandry (observed mating number = mo = 9 or 54), sampling month, and interactions on whole-colony strength measures. Differences were accepted at the α ≤ 0.05 level
Full size table
Fig. 2
figure 2

Experiment-wide effect of sampling month on cm2 brood produced. Analyses were run on square root-transformed data, but natural means ± SEs are shown here. Month lsmeans separated by Tukey–Kramer’s test at α ≤ 0.05

Full size image

For mite drop in 24 h, there were significant effects of month and an interaction between polyandry*genetics. Mite numbers increased from incipient lows in July to peaks in August and September, then dropped in November (Fig. 3). The interaction between polyandry*genetics for mite drop (Fig. 4) shows divergent responses to genetics between polyandry classes. Post hoc t tests showed that within mo = 9, there was no difference in mite counts between genetic lines (t = 1.19, P = 0.2438), whereas within mo = 54, mite counts were significantly lower in VSH colonies (t =  − 2.69, P = 0.0117). There was no effect of the covariate cells of brood per bee.

Fig. 3
figure 3

Experiment-wide effect of sampling month on 24-h Varroa mite counts on bottom board sticky sheets. Analyses were run on log-transformed data, but natural means ± SEs are shown here. Month lsmeans separated by Tukey–Kramer’s test at α ≤ 0.05

Full size image
Fig. 4
figure 4

Interaction between polyandry*genetics for mite counts. Analyses were run on log-transformed data, but natural means ± SEs are shown here. Post hoc t tests showed that within mo = 9, there was no difference in mite counts between genetic lines, whereas within mo = 54 mite counts were significantly lower in VSH colonies

Full size image

For colony bee population, there were significant effects of month with a predictable decline in November. There was a significant interaction between month and genetics. Post hoc t tests showed no differences between genetic groups in July, August, and September, but in November, VSH colonies had significantly more bees than wild colonies (t = 2.7, P = 0.0085) (Fig. 5).

Fig. 5
figure 5

Interaction between sampling month and genetics for colony adult bee population. Analyses were run on square root-transformed data, but natural means ± SEs are shown here. Post hoc t tests showed no differences between genetic groups in July, August, and September, but in November VSH colonies had significantly more bees than wild colonies

Full size image

The Wilcoxon rank sums analysis for brood survival over 3 days found differences by polyandry (χ2 = 7.18, P = 0.0074), with significantly higher brood survival in mo = 54 colonies than in mo = 9 colonies (Fig. 6). There were no differences by genetics (χ2 = 0, P = 0.9841). There were marginal (χ2 = 7.77, P = 0.0510) differences when all four P*G combinations were compared, with a positive change vector in both genetic groups with increasing polyandry (Fig. 7).

Fig. 6
figure 6

Experiment-wide effect of polyandry on survival rate (3-day) of uncapped larvae. Wilcoxon ranked sum values are separated with Kruskal–Wallis test

Full size image
Fig. 7
figure 7

Brood survival by polyandry*genetics. Wilcoxon ranked sum analysis showed only marginal (P = 0.0510) differences when all four P*G combinations were compared; however, natural means ± SEs are shown here to illustrate the change vectors from mo = 9 to mo = 54 across both genetic groups

Full size image

Discussion

By limiting our genotyping to only 50 workers per colony, we most certainly underestimated me (Tarpy and Nielsen 2002). This is especially true of the mo = 54 treatments. Additionally, by constraining the number of drone source colonies, we limited the divergence of genetic diversity possible between polyandry levels. For instance, any six brothers pooled into a mo = 54 insemination were related by r = 0.5. Each such sextet constitutes a comparatively exhaustive sampling of their mother’s alleles, but it remains that all inseminations were made up of sons of no more than nine mothers. Our intention for this constraint was to limit background variation and to treat polyandry in a narrow sense, in distinction from genetic diversity per se. But this protocol also means that effect sizes detected for polyandry in this study are conservative and underestimated.

We used a line of VSH drones as a model source of uncommon specialist alleles, but there were no direct effects of VSH on any dependent variable in our study, including colony levels of the parasite V. destructor against which the VSH stock is specifically selected. VSH did, however, interact with month to deliver the highest colony bee populations by November (Fig. 5). Hyperpolyandry, on the other hand, directly improved brood survival (Fig. 6) and interacted favorably with VSH alleles to decrease overall colony parasite loads (Fig. 4).

Brood production and survival

The results for cm2 brood production did not support our hypothesis (1). There were no significant terms in the model except for the predictable seasonal effect of sampling month (Fig. 2). It is worth noting, however, that Seeley and Tarpy (2007) showed increases in brood as one increased below the GD asymptote from mo = 1 to = 10, and Delaplane et al. (2015) showed an increase in brood rearing efficiency as one increased within the hyperpolyandry zone from mo = 15 to ≥ 30, a pattern consistent with the change vector for cm2 brood in the present study: mo = 9 (2026 ± 1367; means ± SD; range = 0–4532) vs. mo = 54 (2226 ± 1248; range = 0–4136).

However, support for hypothesis (1) is found with brood survival over 3 days for which there was a significant increase from mo = 9 to = 54 colonies (Fig. 6). Brood rearing efficiency is a global outcome of a cluster of common worker phenotypes including ones involved in foraging (Page et al. 1995b), detoxification (Gregorc et al. 2012), nursing (Chapman et al. 2007), thermal regulation, and cell cleaning (Schmickl and Crailsheim 2004), to name a few. Although only marginally significant, results for brood survival when examined by all P*G combinations show that the change vector was positive in both genetic classes when moving from mo = 9 to = 54 (Fig. 7). The sum of evidence from this experiment and others supports that hyperpolyandry facilitates inheritance of gene families in a manner that improves brood rearing efficiency even beyond limits predicted by a simple GV model.

Colony mite levels

When it came to colony levels of V. destructor, there was a strong effect of sampling month, with mite levels tracking the seasonal increase and decrease of available hosts (Fig. 3). There were no direct effects of either polyandry or genetics; however, there was a significant interaction of polyandry*genetics such that the parasite suppressing effect of the VSH genotype expressed only at mo = 54 (Fig. 4), supporting our hypothesis (2).

An increase in the ratio of brood cells per bee is expected to vary positively with mite population growth (Boot et al. 1994; Martin and Kemp 1997), thus presenting an independent source of variation. In our experiment, this positive relationship was non-significant and barely detectable (mite drop = 16.3 + 3.7 (brood per bee), r2 = 0.006, P = 0.3394) and, accordingly, explained no variation in the whole model when included as a covariate. It appears, therefore, that the favorable interaction between hyperpolyandry and specialist mite resistance alleles operates independently of bee population demographics.

Bee population

Bee population can be seen as a sum of outcomes of all tasks, common and rare. Even though there were no differences by genetics in colony bee populations from July through September, by November with the approach of winter colony bee populations were significantly higher in VSH colonies (Fig. 5). As there were no direct or indirect effects of genetics on cm2 brood or brood survival, the only candidate for an early predictor of stronger adult autumn populations in VSH colonies is the favorable genetics*polyandry interaction (Fig. 4) that delivered the best overall parasite suppression to VSH mo = 54 colonies. Given the large negative effect size of Varroa on bee populations (Guzmán-Novoa et al. 2010), it appears that the benefit of hyperpolyandry interacting with the VSH genotype was sufficient to tip overall autumn bee populations in favor of VSH colonies.

Implications for the evolution of hyperpolyandry

By bracketing our queen inseminations within (mo = 9) and above (mo = 54) the GD asymptote (Fig. 1), our results directly show that phenotypic improvements from hyperpolyandry are not constrained by the ceiling on population allele capture imposed by rules of hymenopteran inheritance. Brood survival was significantly higher in mo = 54 than in mo = 9 colonies (Fig. 6), and the VSH genotype reduced colony parasite numbers only under conditions of hyperpolyandry, mo = 54 (Fig. 4). To our knowledge, these are the first direct phenotypic comparisons across this qualitative divide. Our results are consistent with natural hyperpolyandry observed in army ants (Kronauer et al. 2004; Barth et al. 2014) and the Apini (Neumann and Moritz 2000; Tarpy et al. 2004; Withrow and Tarpy 2018). And while they do not refute alternate explanations such as sperm limitation (Simmons 2002) or chance (Tarpy and Page 2000), our results do support colony level adaptive benefits as a direct driver of hyperpolyandry.

It appears that hyperpolyandry increases acquisition of favorable genotype combinations for common tasks (hypothesis 1) in a process that is both independent of within-colony genetic diversity and possibly even continuous. Indeed, with its emphasis on optimal genotypic compositions over average diversity, the model of Fuchs and Moritz (1999) admits no upper limit to mating numbers. In the case of a common task outcome like brood survival, this benefit holds across both genetic types (Fig. 7).

In contrast, effects of the specialist VSH phenotype on mite numbers were detectable only in mo = 54 colonies, consistent with the ideas that hyperpolyandry is adaptive when allele acquisition for rare tasks (hypothesis 2) is more stochastic, i.e., vulnerable to “either/or” cliff edge catastrophes or when the environmental stressor demands high expression frequencies of a rare task.

Hypothesis (1) is consistent with social heterosis—the idea that intracolonial genetic heterogeneity promotes advantageous synergisms of interacting genotypes (Nonacs and Kapheim 2007), but hypothesis (2) calls for evidence of a positive association between queen mating number and environmental constraints in natural populations. A variance-based model by Rueppell et al. (2008) predicts directional selection toward monandry under environmental conditions favoring exceptional colonies, and selection toward polyandry under conditions favoring average colonies, while a model by Gove et al. (2009) suggests that polyandry is favored when conditions demand an increasing to intermediate number of colony tasks, as would be expected in mild environments that support an accumulation of inter-species interactions. Together, these models have been interpreted to mean that monandry is favored in adverse, simple environments and that polyandry is favored in mild, complex environments, and indeed climate-based comparisons of populations of the ant Lasius niger seem to support this (Corley and Fjerdingstad 2011) (see, however, El-Niweiri and Moritz (2011)). A parasite-based context such as ours would seem to align with these modeled expectations to the extent that the parasite is a pioneer (monandry favoring exceptional colonies) or ubiquitous (polyandry favoring average colonies). In any case, our study signals that selection for hyperpolyandry will be sensitive to the kinds of tasks demanded by the environment and the scarcity of their controlling alleles.

Conclusions

Colony fitness increases as a function of the queen’s mating number, but the mechanisms of this genome to phenome translation shift along that mating continuum. Up to levels of me ~ 10, colony fitness improvements are consistent with the genetic variance (GV) hypotheses that see improvements as a function of increasing within-colony average genetic diversity. Mating levels in excess of me = 10 are common in nature, and colony fitness gains in that range are shown here experimentally with brood survival; however, fitness gains from a rare task—behavioral resistance to a dangerous parasite—only express in the presence of the controlling alleles and at the highest polyandry level tested mo = 54. Our results are consistent with at least two mechanistic hypotheses for the evolution of hyperpolyandry in social Hymenoptera: hyperpolyandry improves colony fitness by (1) increasing optimum genotypic compositions for common tasks, and (2) by capturing rare specialist allele combinations, resisting cliff-edge ecological catastrophes.

Ultimately, a formal treatment of our hypotheses will require describing changes to subfamily-based task distributions under the influence of hyperpolyandry and showing direct effects of these changes on colony fitness.