Phylogenomic Evidence Overturns Current Conceptions of Social Evolution in Wasps (Vespidae).

The hypothesis that eusociality originated once in Vespidae has shaped interpretation of social evolution for decades and has driven the supposition that preimaginal morphophysiological differences between castes were absent at the outset of eusociality. Many researchers also consider casteless nest-sharing an antecedent to eusociality. Together, these ideas endorse a stepwise progression of social evolution in wasps (solitary → casteless nest-sharing → eusociality with rudimentary behavioral castes → eusociality with preimaginal caste-biasing (PCB) → morphologically differentiated castes). Here, we infer the phylogeny of Vespidae using sequence data generated via anchored hybrid enrichment from 378 loci across 136 vespid species and perform ancestral state reconstructions to test whether rudimentary and monomorphic castes characterized the initial stages of eusocial evolution. Our results reject the single origin of eusociality hypothesis, contest the supposition that eusociality emerged from a casteless nest-sharing ancestor, and suggest that eusociality in Polistinae + Vespinae began with castes having morphological differences. An abrupt appearance of castes with ontogenetically established morphophysiological differences conflicts with the current conception of stepwise social evolution and suggests that the climb up the ladder of sociality does not occur through sequential mutation. Phenotypic plasticity and standing genetic variation could explain how cooperative brood care evolved in concert with nest-sharing and how morphologically dissimilar castes arose without a rudimentary intermediate. Furthermore, PCB at the outset of eusociality implicates a subsocial route to eusociality in Polistinae + Vespinae, emphasizing the role of mother-daughter interactions and subfertility (i.e. the cost component of kin selection) in the origin of workers.

Introduction

In most eusocial insects, an individual isn’t destined to become a worker or queen based on their genetics (Schwander et al. 2010). Rather, environmental cues they experience during their lifetime fate them to their caste. Polyphenism, in which two or more distinct phenotypes are produced by the same genotype, is prevalent in nature (West-Eberhard 2003; Simpson et al. 2011). Examples include the swarming and solitary phase of locusts, the castes of eusocial insects, and the differentiation of stem cells (Simpson et al. 2011). Despite polyphenism being ubiquitous and fundamental to the success of multicellular life, our understanding of how it originates is limited (West-Eberhard 2003; Suzuki and Nijhout 2006). In addition to insights on polyphenisms, the study of social insects is critical to elucidating how complex adaptive systems with emergent properties (e.g. multicellular life, societies, ecosystems) arise from the interaction of simpler entities. To fully unravel how complex units of biological organization emerged, such as the colonies of eusocial insects, we require a bottom-up approach that first considers the ancestor, and its properties, from which it evolved (van Gestel and Tarnita 2017).

Vespidae, which includes potter wasps, paper wasps, hornets and yellowjackets, is a socially diverse family comprised of solitary, subsocial and eusocial species (Hunt 2007). Furthermore, much variation in the degree of eusociality exists among Vespidae, including facultative, primitive and advanced eusociality (Hunt 2012; Turillazzi 2012; fig. 1). Consequently, Vespidae are recognized as one of the few groups that could elucidate the transitional stages of social evolution, including stages that culminate into crossing the threshold of eusociality and those involved in the advancement of eusociality (Hines et al. 2007; Hunt 2012; Ferreira et al. 2013; Rehan and Toth 2015). Thus, the phylogeny of vespid wasps provides a framework for testing how eusocial groups originate and evolve.

. 1.

The traditional phylogeny of Vespidae (Carpenter 1982; Pickett and Carpenter 2010) supporting a single origin of eusociality. The ground plan condition of all eusocial wasps is inferred to be rudimentary castes with no preimaginal caste-biasing (PCB). That is, in the earliest eusocial ancestor, discrete physiological differences between workers and queens were established only during adulthood. In solitary vespids, mothers nest alone and do not interact with their offspring. Facultative eusociality is a term exclusive to describing stenogastrine social organization, where all females are morphologically similar and able to produce their own offspring; helper status is often temporal, wherein young females are helpers until their ovaries have developed; colony size is exceptionally small with a maximum of six females on a nest in most species (Hunt 2012; Turillazzi 2012). Primitive eusociality refers to independent-founding paper wasps: worker status is typically lifelong even though workers are reproductively plastic, and in some cases queens and workers have morphological differences. Highly eusocial species (i.e. advanced eusociality) have preimaginal determination of a life-time unmated worker caste that is typically morphologically distinct. In swarm-founding species, swarms of worker wasps are responsible for initiating new nests rather than queens; in some cases, they have morphologically differentiated castes. (*) Not all Ropalidiini are primitively eusocial; some are swarm-founders.

The first phylogenetic treatment of Vespidae recovered Stenogastrinae as most closely related to the social Polistinae and Vespinae (Carpenter 1982; fig. 1). This traditional view is consistent with nearly the entire taxonomic history of vespids, from de Saussure’s (1855) pre-evolutionary “Disposition naturelle” diagram with stenogastrines (Ischnogaster) grouped with the social wasps, to Richards’ (1962) classification and tree with stenogastrines as sister to Polistinae + Vespinae. Only Richards (1971), who speculated stenogastrines evolved from a Eumenes-like solitary ancestor, Spradbery (1975), who suggested evolution from an early vespoid ancestor, and Van der Vecht (1977), who argued they were closely related to zethines, had challenged the affinity of stenogastrines to the polistines and vespines. Despite stenogastrines being recognized as sister to polistines and vespines (Carpenter 1982), the single origin of eusociality hypothesis dates to Carpenter (1988). As knowledge of stenogastrine ethology accumulated it became clear that eusociality was the ground plan of stenogastrines and therefore evolved once in Vespidae (Carpenter 1988).

The first molecular study to suggest Stenogastrinae was distantly related from the other eusocial wasps (and thus eusociality evolved twice) leveraged evidence from nuclear 28S rDNA and mitochondrial 16S rDNA (Schmitz and Moritz 1998). This study was dismissed due to small taxon sampling (only 13 vespids), an anomalous result (honey bees placed within Vespidae), and parsimonious realignment of their data set supporting a single origin of eusociality (Carpenter 2003). The second molecular study to suggest two origins of eusociality utilized sequences from four genetic loci (28S, 18S, abdominal-A, and RNA polymerase II) across 27 vespid species (Hines et al. 2007). However, simultaneous analysis of that genetic data set with morphological and behavioral data favored the single origin hypothesis (Pickett and Carpenter 2010). Furthermore, simultaneous analysis of an independent and more comprehensive data set, with data from four loci (CO1, 28S, 12S, 16S) and 333 phenotypic characters for 130 vespid species, also supported the traditional (one origin) phylogeny (Pickett and Carpenter 2010). A more recent study reanalyzed the Pickett and Carpenter (2010) data set and confirmed that simultaneous analysis of molecular and phenotypic evidence does support a single origin of eusociality, but that analysis of the molecular data alone supports a dual origin (Piekarski et al. 2014). It was clear that molecular and phenotypic evidence were at odds, but that the single origin hypothesis is favorable if considering all the available evidence (Piekarski et al. 2014). A single origin of eusociality has received further support as morphological data accumulates (da Silva et al. 2014), and is also the result following from the most recent paleontological analysis of Vespidae (Perrard et al. 2017) and a phylogenetic analysis of tribal relationships within Eumeninae (Hermes et al. 2014).

More recently, a phylogenomic study exploring relationships among Hymenoptera contained six vespid exemplars and recovered relationships in support of the dual origin hypothesis (Peters et al. 2017). Likewise, another recent phylogenomic study (Bank et al. 2017) suggested two origins of eusociality, with taxonomic sampling of 49 exemplars (but only five being eusocial). Only one representative of Stenogastrinae was included in that study, and it is the placement of this subfamily that has been the most questioned and important. As limited taxon sampling can cause substantial systematic error and strong support for erroneous relationships (Nabhan and Sarkar 2012), especially through long branch attraction (Graybeal 1998; Hillis et al. 2003), this result must be taken cautiously. Nonetheless, Bank et al. (2017) provided a substantial body of new evidence in favor of the dual origin of eusociality hypothesis.

The polygynous family hypothesis (West-Eberhard 1978a) and subsocial hypothesis (Wheeler 1923) are competing ideas regarding the origin of reproductive castes in Vespidae. The polygynous family hypothesis has two central tenets: interactions amongst nest-sharing females of multiple generations were involved in the incipience of a worker caste, and casteless nest-sharing preceded eusociality. In this scenario, reproductive castes evolved gradually in a polygynous (multiple egg-layers on a shared nest) context. Some have argued that the presence of polygyny in a few species of potter wasps (Eumeninae) is suggestive of castes in eusocial vespids having originated in a casteless nest-sharing context (West-Eberhard 1978a). In contrast, the subsocial hypothesis (Wheeler 1923) posits that mother–daughter interactions are salient to the incipience of a worker caste, and that reproductive division of labor (i.e. castes) evolved simultaneously with nest-sharing.

Espoused by many researchers is the idea that social evolution is incremental, with intermediate stages between solitary life and advanced eusociality (Wilson 1971; West-Eberhard 1978a; Pardi and Marino Piccioli 1981; Wheeler 1986; Nowak et al. 2010; Ferreira et al. 2013; Kapheim, Pan, et al. 2015; Rehan and Toth 2015; Boomsma and Gawne 2018). It is posited that social evolution in vespids progressed in a stepwise manner, such that: 1) the transition from solitary to eusocial behavior involved a casteless nest-sharing intermediary stage (West-Eberhard 1978a; Gadagkar 1991); and 2) the beginnings of vespid eusociality were characterized by a lack of preimaginal caste-biasing (PCB), where all females eclose undifferentiated and have relatively equal reproductive physiology (Pardi 1948; Carpenter 1991; Field and Foster 1999; Noll et al. 2004). Workers of primitively eusocial wasps (e.g. Polistes) retain the ability to reproduce (i.e. reproductively plastic) and can serve as replacement queens, and so it is often posited that such wasps lack PCB (Pardi 1948; Reeve 1991; Sumner et al. 2010). However, the presence of reproductive plasticity in workers does not rule out PCB, and so absence of PCB as the ground plan for primitively eusocial wasps has been challenged before (Gadagkar 1991; Hunt 1991; O’Donnell 1998; Hunt 2006, 2007; Hunt et al. 2007). These challenges have gone largely unheeded because the traditional phylogeny with a single origin of eusociality implies rudimentary castes (i.e. no PCB) as the ground plan condition for all eusocial vespids (fig. 1). However, if the single origin of eusociality hypothesis is demonstrably false, then the possibility of PCB at the ancestor of Polistinae + Vespinae needs to be re-evaluated.

We present a phylogenomic data set for Vespidae with much broader taxon sampling compared with recent phylogenomic studies (Bank et al. 2017; Peters et al. 2017). Our study utilizes different taxa, loci, and phylogenetic methods than Bank et al. (2017). In addition to corroborating whether eusociality evolved twice, our study aims to produce a comprehensive phylogeny that will enable phylogenetically grounded comparative approaches using burgeoning -omics data (Ferreira et al. 2013; Berens et al. 2015; Patalano et al. 2015; Standage et al. 2016) to elucidate the molecular basis of social evolution in wasps. We also present an extensive probe set that captures DNA from a wide range of hymenopteran taxa, which will be of utility to many research groups that study ants, bees and even parasitoid wasps. Further, we provide a resolution to the competing polygynous family (West-Eberhard 1978a) and subsocial (Wheeler 1923) hypotheses regarding the origin of reproductive castes in Vespidae. Lastly, we test whether rudimentary and monomorphic castes characterized the early stages of eusocial evolution, as is the current supposition (West-Eberhard 1978a; Carpenter 1991).

Results and Discussion

Origins of Eusociality and PCB

For a detailed discussion on the recovered phylogenetic relationships and taxonomic changes [necessitated by the results, see supplementary results, Supplementary Material online]. Counter to the traditional phylogeny (Carpenter 1982, 1988, 2003; Pickett and Carpenter 2010; Hermes et al. 2014; Piekarski et al. 2014; Perrard et al. 2017), our robust molecular phylogeny inferred from novel genetic data supports two independent origins of eusociality in Vespidae (fig. 2). In corroboration of previous molecular studies (Schmitz and Moritz 1998; Hines et al. 2007; Bank et al. 2017; Peters et al. 2017), Stenogastrinae is sister to all other vespids, and thus its social habits are independently evolved from those of Polistinae and Vespinae. With this well-supported result across studies, Stenogastrinae should no longer serve as a phylogenetic intermediate between a solitary ancestor and the eusocial ancestor of polistines and vespines. Prior to Bank et al. (2017) and this study, the lack of resolution to the single versus dual origin debate has for decades prevented a unified effort among researchers to study vespid eusociality through the lens of dual origins. A resolution of dual origins prompts us to re-evaluate the supposition that PCB is derived and thus had no part in the origins of reproductive castes.

. 2.

Maximum-Likelihood tree of Vespidae inferred from 235 selected loci sequenced across 163 taxa (showing only 138 ingroup taxa). A six-celled box summarizes branch support values (bootstrap or posterior probabilities), where each cell corresponds to a separate phylogenetic analysis (ML = Maximum Likelihood; AST = ASTRAL; MP = Maximum Parsimony; 228, 235, and 378 = # of loci used) and its color represents support obtained in that analysis (white = 0–69; pink = 70–89; purple = 90–100). Branches that had >90 support in all six analyses are unmarked. Previous taxonomic classifications on the right; new classifications following the color code. For trees recovered in individual analyses, including outgroup relationships and support values, see supplementary figures S6–S11, Supplementary Material online.

It is important to first clarify the different usages of the term PCB in the literature, and which usage we adopt here. In its first usage, PCB refers to continuous variation in traits related to reproductive potential (e.g. fecundity, size, aggression) stemming from ontogenetic influences (e.g. nutrition), which results in some variants being more primed for reproduction and others for allomaternal care (Kapheim et al. 2012). In this usage, PCB is universal in eusocial insects given that environmental factors experienced during ontogeny affect the reproductive output of adults even in solitary insects (Cowan 1981; Tibbetts et al. 2013). Such phenotypic variation across individuals was undoubtedly conducive to the origin of reproductive castes (Sakagami and Maeta 1987; Kapheim et al. 2012). In its second usage (the usage we adopt here), PCB refers to cases where females, prior to becoming adults, exhibit discrete phenotypes (i.e. a polyphenism) that are at least physiologically decoupled (Hunt et al. 2007; Hunt et al. 2010). Morphologically discrete workers and queens are clear evidence of PCB, but PCB need not always produce discrete castes at a morphological level (O’Donnell 1998; Hunt et al. 2007). In contrast to PCB, preimaginal caste-determination necessitates a loss of reproductive plasticity in workers (i.e. life-time unmated workers).

From the ancestral state reconstructions, PCB was reconstructed as present in the ancestor of Polistinae + Vespinae (table 1, supplementary fig. S1, Supplementary Material online) despite assuming absence of PCB as the default state for paper wasps (proportional likelihood for presence: Assym.2 = 0.9760). A likelihood ratio test identified the Assym.2 model as superior for this trait (χ2 = 29.73, df = 1, P = 0.001). When ambiguous cases were conservatively coded as unknown, presence of PCB was again the most likely ancestral state (proportional likelihood for presence: Assym.2 = 0.9783; χ2 = 22.21, df = 1, P = 0.001; table 1, supplementary fig. S2, Supplementary Material online). Morphological castes, where workers are on average different from queens, were also reconstructed as present in the ancestor of Polistinae + Vespinae (proportional likelihood for presence: Assym.2 = 0.9655; χ2 = 24.67, df = 1, P = 0.001; table 1, supplementary fig. S3, Supplementary Material online). Presence of morphological castes was the most likely ancestral state for Polistinae + Vespinae in all taxon sampling schemes, but the result was not significant in all analyses (table 1). When forward and reverse transitions are constrained to be equal, both maximum likelihood reconstructions and stochastic character mapping (Revell 2013) retrieve monomorphic castes as most likely in the ancestor of Polistinae + Vespinae (table 1; supplementary fig. S4, Supplementary Material online). However, when forward and reverse transitions are assigned separate estimated rates, presence of morphological castes is most likely (supplementary figs. S3 and S5, Supplementary Material online). Although incomplete taxon sampling affects ancestral state reconstructions, this is the most comprehensive test to date on the supposition that castes were rudimentary at the outset of eusociality. Despite conservative coding that biases the result to favor an absence of PCB and morphological castes at the ancestor of Polistinae + Vespinae, presence of both is retrieved as the most likely ancestral state (table 1; supplementary figs. S1–S3, Supplementary Material online). Consequently, the idea that castes with ontogenetically established morphophysiological differences surfaced alongside nest-sharing in the shared ancestor of Polistinae + Vespinae must be considered. If true, a conception of stepwise social evolution in vespids can be rejected.

Table 1.

Likelihood of Presence for Morphological Castes (MC) and Preimaginal Caste-biasing (PCB) at the Common Ancestor of Polistinae and Vespinae.

Trait	Taxa Used	Model	Likelihood of Presence	χ2
PCB1	All vespids
		Mk1	0.6134
		Assym2	0.9760*	29.7280
		Parsimony	Ambiguous
	Z + P + V
		Mk1	0.7507
		Assym2	0.8019	1.6618
		Parsimony	Ambiguous
	P + V
		Mk1	0.7012
		Assym2	0.7414	0.3867
		Parsimony	Present
PCB2	All vespids
		Mk1	0.6916
		Assym2	0.9783*	22.2113
		Parsimony	Ambiguous
	Z + P + V
		Mk1	0.7891
		Assym2	0.8217	0.5152
		Parsimony	Ambiguous
	P + V
		Mk1	0.7377
		Assym2	0.8280	2.1951
		Parsimony	Ambiguous
MC	All vespids
		Mk1	0.1939
		Assym2	0.9655*	24.6741
		Parsimony	Absent
	Z + P + V
		Mk1	0.6777
		Assym2	0.6490	1.4050
		Parsimony	Absent
	P + V
		Mk1	0.6608
		Assym2	0.6824	0.0479
		Parsimony	Present

Note.—Two ancestral state reconstructions were performed for PCB. In the first (PCB1), ambiguous cases were coded as absent, and in the second (PCB2) ambiguous cases were coded as unknown. Reconstructions were done on three taxon sets: all vespids; only Zethinae, Polistinae, and Vespinae (Z + P + V); or only Polistinae and Vespinae (P + V). Asterisk (*) indicates support as the ancestral state by a likelihood decision threshold of 2.0. χ2 values correspond to a likelihood ratio test comparing Mk1 and Assym2 models. Significantly better (P < 0.05) models are italicized. See also supplementary figures S1–S3, Supplementary Material online.

A Developmental Bifurcation and Subfertile Females at the Outset

Counter to the view that PCB is characteristic of derived species (Pardi and Marino Piccioli 1981; Noll et al. 2004; Sumner et al. 2010), our results support PCB as the ancestral condition of paper wasps, yellowjackets and hornets (table 1; supplementary figs. S1 and S2, Supplementary Material online). In the presence of PCB, physiological and sometimes morphological differences between castes are established before adulthood (O’Donnell 1998; Hunt et al. 2007, 2010). This can include differences in number of Malpighian tubules, Van der Vecht organ shape, levels of hexameric storage proteins, maturation time, timing of ovarian activation, body size, fat content, and responsiveness to juvenile hormone (O’Donnell 1998; Hunt et al. 2007; Tibbetts and Sheehan 2012; Judd et al. 2015; de Souza et al. 2016). The latter three differences are known to affect reproductive potential (Cowan 1981; Markiewicz and O’Donnell 2001; Arrese and Soulages 2010; Tibbetts et al. 2011; Shukla et al. 2013; Tibbetts and Sheehan 2012). Furthermore, body size and fat content are influenced by nutritional input (Wheeler 1986; Karsai and Hunt 2002; Judd et al. 2015). Given a solitary ancestor where size, fat content, or both are associated with fertility (Cowan 1981; Honěk and Honek 1993; Tibbetts et al. 2013), mothers could manipulate daughters to be workers by feeding them less (Alexander 1974)—or, more accurately, mothers could exploit poorly nourished daughters (Hunt 1991; Karsai and Hunt 2002). If PCB was the ancestral condition of Polistinae + Vespinae, then a subset of daughters expressed a distinct, likely subfertile, phenotype at the outset of eusociality.

It has previously been argued that preimaginal differentiation of workers and queens is an underlying, unifying feature of polistines and vespines (Gadagkar 1991; Hunt 1991; O’Donnell 1998; Hunt 2006, 2007; Hunt et al. 2007). The results (fig. 2; supplementary figs. S1 and S2, Supplementary Material online) support this proposed framework of eusocial evolution, whereby two distinct developmental trajectories exist for both polistine and vespine female offspring during larval development (Hunt 2007): the putative worker (nondiapause) and the gyne (diapause) trajectory. Since gynes that undergo diapause are the future foundresses, they are more likely to become queens than nondiapausing females. Thus, the putative worker developmental trajectory biases a female to become an “ontogenetic worker” (Hunt 1991, 2006), but she remains reproductively plastic and can become a replacement queen (Reeve 1991). The gyne trajectory primes females to be the future foundresses, but a subset of foundresses will become subordinate cofoundresses that behave as workers (Reeve 1991). PCB is evident in temperate and subtropical independent-founding species, where gynes eclose with a distinct physiology that primes them for diapause (Hunt and Amdam 2005; Gobbi et al. 2006; Hunt et al. 2007; Judd et al. 2015). However, it may be that overwintering only serves to sharpen the distinctiveness of the gyne phenotype, and that rudiments of the developmental bifurcation between putative workers and gynes are also present in primitively eusocial paper wasps living in nonseasonal tropical regions (Hunt and Amdam 2005; Hunt 2007). For an overview of previous observations suggesting a developmental bifurcation of female larvae (i.e. PCB) in tropical primitively eusocial wasps, see supplementary discussion,Supplementary Material online.

There are many examples of polistines with morphologically distinct workers that are reproductively plastic, including various independent-founding Ropalidia (Wenzel 1992; Fukuda, Kojima, Tsuchida, et al. 2003), Belonogaster (Pardi and Marino Piccioli 1981; O’Donnell 1998; Keeping 2002), and Polistes (West-Eberhard 1969; Gobbi et al. 2006; Tibbetts 2006; Hunt 2007; de Souza et al. 2016), as well as some swarm-founding epiponines (Noll et al. 2004; Noll and Wenzel 2008), Ropalidia (Yamane et al. 1983; Kojima and Tsuchida 2000; Fukuda, Kojima, Jeanne 2003), and probably Polybioides (Turillazzi et al. 1994). Thus, morphological differences between castes can evolve despite females being reproductively plastic. Despite the presence of morphophysiological differences between castes at the origin of eusociality in Polstinae + Vespinae, we do not claim superorganismality (sensuWheeler 1911) at this origin. Preimaginal determination of life-time unmated workers demarcates the major transition to superorganismality, because permanent unmatedness of workers is what makes them analogous to somatic cells (Wheeler 1911; Boomsma and Gawne 2018). We also do not claim there is evidence for reversals of superorganismality in Vespidae, but of morphological differences in caste. When looking at specific morphological traits that differ between workers and queens, it is uncontroversial to say that they come and go throughout the history of a given lineage (examples within Epiponini shown in Noll and Wenzel 2008). A size dimorphism is apparently absent in Vespa analis and Vespa tropica (Matsuura 1991), which are not closely related to each other (Perrard et al. 2013). However, V. tropica does have a caste dimorphism with respect to wing venation (Perrard et al. 2012). If most or all morphological trait differences between castes are lost, there would be an apparent reversal of morphologically differentiated castes. As an example, our results show a loss of morphological castes within the swarm-founding Epiponini (node B in supplementary fig. S3, Supplementary Material online).

A possible explanation for the lack of morphophysiological queen/worker dimorphism in all of Stenogastrinae may be the longer developmental time of brood and reduced number of brood that are simultaneously reared by a single foundress compared with polistines and vespines (Turillazzi 2012). These aspects of stenogastrine biology raise the worker-to-larva ratio in the early stages of the colony cycle, meaning brood raised by a single foundress will generally not be malnourished. This would also explain why the subfertility hypothesis (West-Eberhard 1975) is inapplicable for explaining worker behavior in Stenogastrinae (Field and Foster 1999).

Nest-Sharing and Reproductive Castes Evolved in Concert

Since our results suggest PCB was present at the outset of eusociality (table 1), this implies that nest sharing and reproductive castes originated simultaneously and the worker caste of Polistinae + Vespinae evolved as a result of interactions between adults (mothers) and larvae (their daughters) in a context of solitary nesting (Hunt and Amdam 2005; Hunt 2007); the worker caste corresponds to the subfertile daughters that are reared early in the colony cycle. This contrasts with the polygynous family hypothesis (West-Eberhard 1978a), in which the worker caste emerged from interactions among nest-sharing adult females of multiple generations in a relatively casteless setting, wherein losers of social contests would become more worker-like (i.e. lay less eggs, and forage more). The observation of casteless nest sharing in at least three species of eumenines is interpreted by some researchers as evidence supporting the polygynous family hypothesis, in which the ancestor of eusocial vespids passed through an intermediary stage of communal nesting (West-Eberhard 1978a, 2005a). A previous molecular study recovering Zethinae as sister to Polistinae and Vespinae interpreted the close relationship to mean that the social system exhibited by Zethus miniatus could have been reminiscent of the transitional stage between solitary life and polistine eusociality (Hines et al. 2007). However, Z. miniatus is the only zethine known to be communal and can be extrapolated as a derived lineage within Zethinae (Lopes and Noll 2018; also fig. 2: Z. miniatus belongs to the subgenus Zethoides, represented in our phylogeny by Z. binodis and Z. toltecus). Given that only one species of Zethinae is known to exhibit communal nesting, the ancestral state of Zethinae was unlikely to be communal nest-sharing. Thus, the social system of Z. miniatus evolved independently from polistine/vespine eusociality. Consistent with Wheeler’s (1923) subsocial model, the inferred phylogeny (fig. 2) and ancestral state reconstructions (table 1; supplementary figs. S1–S3, Supplementary Material online) place mother–daughter interactions and daughter subfertility (daughters are malnourished due to mothers simultaneously provisioning multiple offspring) as the foundation from which eusociality in Polistinae + Vespinae emerged.

The evolution of a worker caste via the polygynous family hypothesis is also not applicable to Stenogastrinae. The predominant mode of nest initiation in hover wasps is by a single foundress (i.e. haplometrosis), and thus interactions among same generation females likely had no part in the origin of reproductive castes (contraRoubaud 1916; West-Eberhard 1978a). Also, most eggs in a shared nest are produced by a single female (i.e. monogyny) despite the presence of other fertilized females (Sumner et al. 2002; Bolton et al. 2006). In hover wasps, there are clear behavioral castes and females are reproductively plastic (Turillazzi 2012). Therefore, the ground plan of hover wasps is haplometrosis, short-term monogyny, and rudimentary castes; not casteless nest sharing.

Mechanisms for Morphologically Distinct Castes at the Outset

Cooption of a developmental switch regulating diapause in a bivoltine ancestor (Hunt and Amdam 2005; Hunt 2006; Hunt et al. 2007, 2010) is one means by which castes would necessarily start off preimaginally biased and physiologically decoupled. If a pre-existing seasonal polyphenism (i.e. bivoltinism) was coopted for the evolution of castes (Hunt and Amdam 2005; Quinones and Pen 2017), then the underlying morphophysiological differences between the diapause (gyne) and nondiapause (putative worker) phenotype would have preceded any behavioral specialization of castes. It is also possible that a diapause/nondiapause switch, known to be deeply conserved among insects (Denlinger 1986; Hahn and Denlinger 2011), served as a gateway to castes even if the preceding ancestor was not bivoltine. Regardless of whether the diapause ground plan was coopted in the advent of reproductive castes, the accumulation of trait differences between castes was not necessarily a gradual process based on sequential mutation.

Among evolutionary biologists, the “plasticity-first” hypothesis posits that phenotypic plasticity precedes and facilitates most adaptation, rather than mutation (Levis and Pfennig 2016). That is, novel traits may first appear as developmental variants in response to a changing environment—environment includes developmental contexts caused by other traits which may or may not have originated directly from a mutation. Phenotypic accommodation, adaptive adjustment without genetic change, is potentially a significant source of phenotypic innovation (West-Eberhard 2003, 2005b; Levis and Pfennig 2016), and has been postulated to be of chief importance in the emergence of vespid eusociality (West-Eberhard 1987; Wilson 2008; Hunt 2012). Through phenotypic plasticity, a novel environment or trait can trigger the expression of a new (previously cryptic) phenotypic variant. This process could lead to the simultaneous incipience of many novel traits and caste differences in the absence of a mutation event specific to each newly emerged trait/difference (West-Eberhard 2005b). That is, phenotypic plasticity could explain the synchronized appearance of simultaneous progressive provisioning (mothers provide provisions in increments to multiple offspring at once), nest-sharing among adults, cooperative (allomaternal) brood care, and reproductive division of labor, as well as the morphophysiological differences between castes at the outset of eusociality recovered in our results (see supplementary fig. S3, Supplementary Material online).

Genetic accommodation follows phenotypic accommodation, whereby selection refines the novel phenotypic variants (i.e. traits) that were a consequence of plasticity by acting on previously cryptic alleles (Levis and Pfennig 2016). Since the cost of altruism is minimized for subfertile females (West-Eberhard 1975), reduced size and fertility among first cohort daughters, an inherent feature of polistines and vespines, would bolster kin selection (Hamilton 1964a, 1964b) to favor genetic variants that express allomaternal care under malnourished conditions. Even if altruism is disadvantageous to the subfertile daughter, selection could still theoretically favor altruism if the benefits to mothers rearing altruistic subfertile daughters are high (West-Eberhard 1975; Kapheim, Nonacs, et al. 2015). We hypothesize that pre-existing natural variation of females in variables such as size and fertility, which in Polistinae + Vespinae was heightened by the appearance of simultaneous progressive provisioning, decreased the cost of altruism for a subset of daughters and thus promoted the spread of facultative (i.e. condition dependent) altruism via kin selection. In the presence of facultatively expressed allomaternal care, the developmental mechanisms producing a pre-imaginal bifurcation among females would be reinforced by genetic accommodation.

Conclusions

Our results demonstrate that eusociality emerged twice within Vespidae. For both eusocial lineages there is no indication of casteless nest-sharing (i.e. polygyny) preceding eusociality. Further, the results suggest that at the outset of eusociality in the shared ancestor of Polistinae + Vespinae, a class of subfertile daughters (i.e. putative workers), likely with morphophysiological differences (including size differences that can arise via plasticity), were produced early in the colony cycle. Thus, a conception of stepwise social evolution in vespid wasps, which holds that casteless nest-sharing was a precursor to eusociality and PCB was absent in the early stages of polistine/vespine eusociality, is unfounded. We hypothesize that nest-sharing, altruism, and reproductive division of labor, whether castes were physically differentiated or not, appeared in concert via “plasticity-first,” likely in response to simultaneous progressive provisioning. We suspect that future studies exploring the role of phenotypic plasticity and standing genetic variation in the origins of eusociality will further corroborate the saltatory nature of social evolution and caste divergence.

Materials and Methods

Taxon Sampling and Probe Design

As a starting point, we leveraged the 941-locus Coleoptera alignments developed by Haddad et al. (2018) for anchored hybrid enrichment (AHE). Methods for AHE locus selection and probe design followed Hamilton et al. (2016). Using the red flour beetle (Tribolium castaneum) sequence as a reference, we scanned the assembled genomes of two divergent hymenopteran species: the bumblebee Bombus impatiens and the parasitoid wasp Nasonia vitripennis. For each of the 941 Coleoptera loci, the best matching region in each of the two Hymenoptera genomes was identified and used to extract up to a 4,000 bp region that contained that locus. The extracted regions were then aligned using MAFFT (v7.023b, with –genafpair and –maxiterate 1000 flags, Katoh and Standley 2013).

To increase representation of lineages across Hymenoptera, we scanned 12 assembled and 11 unassembled genomes (supplementary material data set S1, Supplementary Material online) for anchor regions using the N. vitripennis sequences as a reference from the pairwise alignments generated above. For the 12 assembled genomes, we isolated up to 4,000 bp surrounding the region that best matched the reference. For the unassembled genomes, we first prepared indexed libraries with ∼250 bp inserts following Lemmon et al. (2012), then sequenced to ∼15× genomic coverage on two HiSeq2500 lanes (79 Gb total output) with 8 bp indexing and paired 100 bp reads. After reads were merged and trimmed following Rokyta et al. (2012), we mapped the reads to the N. vitripennis sequences. The consensus sequences from the resulting assemblies were then extended into the flanks, producing up to a 4,000 bp sequence for each species at each locus. Using MAFFT, we generated an alignment for each locus that contained N. vitripennis plus consensus sequences from each of 23 target species. After visual inspection in Geneious (Kearse et al. 2012), we trimmed and masked the alignments to exclude poorly aligned sequences. To ensure sufficient enrichment efficiency, we removed loci with <50% of the taxa. This process resulted in 528 anchor loci and 13 functional loci (total target size is 212,392 bp). Following Hamilton et al. (2016), we analyzed Kmer profiles to identify and mask repetitive alignment regions. Lastly, we tiled 120 bp probes at 2× density across each of the 24 sequences for each of the loci to produce 57,066 probes. Probe sequences are available at https://dx.doi.org/10.6084/m9.figshare.c.4135511, Last accessed June 22, 2018.

We obtained sequence data for 137 species of wasp (136 ingroup, 1 outgroup) out of 142 that were processed through the AHE workflow (supplementary material data set S2, Supplementary Material online). Publicly available hymenopteran genomes and genomic sequences from 11 ichneumonoid wasps that were collected for probe design were also included for analysis, which added two ingroup and 24 outgroup representatives (supplementary material data set S1, Supplementary Material online). Thus, a total of 163 taxa are incorporated into phylogenetic analyses, 138 of which are ingroup taxa. All recognized subfamilies and tribes of Vespidae were sampled (except Raphiglossinae, which at the time of collecting samples was not considered a distinct tribe or subfamily). More taxa are sampled here than any other previous molecular phylogenetics study of Vespidae, with samples collected from a total of 26 countries (supplementary material data set S2, Supplementary Material online). Our sampling is also broad at shallow timescales, and includes representatives from all vespine genera, 24 of 25 polistine genera, four of seven stenogastrine genera, and even all Polistes subgenera.

Genomic DNA was extracted from specimens following Qiagen protocols in conjunction with the DNeasy Tissue Kit (Qiagen, Valencia, CA, U.S.A.). DNA was successfully extracted from dry and ethanol preserved specimens with collection dates ranging from 1987 to 2015. Voucher specimens are deposited in a variety of institutions (supplementary material data set S2, Supplementary Material online).

Library Preparation, Enrichment, and Sequencing

Genetic data were obtained at the Center for Anchored Phylogenomics at Florida State University (http://www.anchoredphylogeny.com; Last accessed June 22, 2018) using the general methods outlined in Lemmon et al. (2012). First, each genomic DNA sample was sonicated to a fragment size of 150–350 bp using a Covaris E220 focused-ultrasonicator. Subsequently, library preparation and indexing were performed following a protocol outlined in (Meyer and Kircher 2010). Indexed samples were then pooled at equal quantities determined using a Qubit fluorometer (16 samples per pool), and samples in each pool were enriched using an Agilent Custom SureSelect kit (Agilent Technologies) with custom designed probes. The general enrichment process uses streptavidin-coated magnetic beads to isolate targeted genomic fragments that hybridize with the oligonucleotide probes (Gnirke et al. 2009). Enriched fragments were pooled in groups of three (three pools total with 48 samples in each) before sequencing on the PE150 Illumina HiSeq2500 (three lanes, 48 samples per lane). Sequencing was performed in the Translational Science Laboratory in the College of Medicine at Florida State University.

Data Processing

Data were processed following methods described in Prum et al. (2015) and Hamilton et al. (2016). Briefly, a bioinformatics pipeline completes a workflow that: merges overlapping paired-end reads (details given in Rokyta et al. 2012), assembles reads into contigs and generates a consensus sequence for each locus per sample (details given in Hamilton et al. 2016), filters out consensus sequences derived from assemblies with low coverage (i.e. <25 reads mapped), determines orthologues (with a minimum of 80 sequences required; details in Hamilton et al. 2016), aligns sequences within each orthologous set using MAFFTv7.023b (Katoh and Standley 2013), and trims/masks the alignment for each locus (with mingoodsites = 11, minpropsame = 0.8 and 75 missing taxa allowed; details in Hamilton et al. 2016). For a summary of the sequencing and assembly results, see supplementary material data set S1, Supplementary Material online.

Phylogenetic Analyses

Gene Tree Generation

The best evolutionary model for each locus was determined using jModelTest 2.1.7 (Darriba et al. 2012). The loci were not partitioned by codon position, and the single model with the highest BIC score was assigned to each locus. Gene trees were then inferred using MrBayes3.2.6 (Ronquist et al. 2012). A total of 378 gene trees were estimated. See supplementary materials and methods, Supplementary Material online for specific settings applied in MrBayes to generate gene tree estimates.

Locus Exclusion

Gene trees that did not reach stationarity, perhaps due to model inadequacy or homoplasy, were excluded in one of the coalescent-based analyses (see below). A locus was considered to have failed the stationarity test if at least one criterion was true: 1) Average standard deviation of split frequencies above 0.05; 2) Potential scale reduction factor for any one estimated parameter above 1.025; 3) Average estimated sample size for any one estimated parameter below 100. Only 12 of 378 loci failed to reach stationarity (not shown). Ambiguously aligned regions were removed using Gblocks 0.91b (Castresana 2000), with details described in supplementary materials and methods,Supplementary Material online. Sixty of 378 loci were entirely trimmed by Gblocks. For the remaining 318 loci a base composition heterogeneity test was performed on the original alignments within PAUP4.0b10 (Swofford 2002). Loci with base composition heterogeneity introduce systematic error (Van Den Bussche et al. 1998; Romiguier et al. 2016), and so 83 loci that failed the heterogeneity test were excluded.

Coalescent-based Analyses

Species tree estimations were done using ASTRAL4.10.6 (Mirarab et al. 2014; Mirarab and Warnow 2015; Sayyari and Mirarab 2016). Two ASTRAL trees were estimated. The first is based on all 378 gene trees, and the second on 228 gene trees (83 loci had base heterogeneity, 60 dropped post gblocks and seven failed the stationarity test). Both analyses used gene trees that were generated from original alignments, prior to trimming with Gblocks. Local posterior probabilities calculated within ASTRAL were favored over multi-locus bootstrapping (MLBS) support values (Sayyari and Mirarab 2016). These posterior probabilities have higher precision (percentage of supported branches that are correct) and recall (the percentage of all true branches that are supported) compared with MLBS, and are not prone to false positives even with high levels of gene tree estimation errors (Sayyari and Mirarab 2016).

Concatenation Based Analyses

Parsimony analyses were performed within TNT v1.5 (Goloboff and Catalano 2016). Analyses were done on two primary data sets consisting of 378 loci (all loci) and 235 loci (excludes 83 loci with base heterogeneity and 60 loci that dropped post gblocks). The second data set was composed of Gblock refined loci, not original alignments. The concatenations are of the following length: 378 loci = 162,919 bp; 235 loci = 63,712 bp (supplementary table S1, Supplementary Material online). For details on the tree search method and TNT settings used, see supplementary materials and methods,Supplementary Material online. The best tree length scores recovered were 364,354 for the 235-gene tree (CI = 0.161, RI = 0.653), and 963,153 for the 378-gene tree (CI = 0.181, RI = 0.655). After the tree search, 1,000 standard bootstraps were performed using traditional searches (defaults), with support values output as frequency differences (Goloboff et al. 2003).

Maximum Likelihood

PartitionFinderV2.0.0 (Lanfear et al. 2016) was used to find the best-fit models and partitioning scheme for the two concatenations. Each locus was treated as a subset, for a total of either 378 or 235 subsets. PartitionFinder was run with and without a user-specified fixed topology (the TNT parsimony tree inferred from 378 loci). Also, the kmeans algorithm (Frandsen et al. 2015) was implemented, with the entire concatenation treated as one subset, and the 378 gene parsimony tree used as the starting tree. Subsequently, PartitionFinder was run using the subsets retrieved by the kmeans algorithm without a starting tree specified. For each treatment an rcluster search (Lanfear et al. 2014) was performed, with branch lengths set as linked and model selection as BIC. The treatment with the lowest BIC score was favored (supplementary table S2, Supplementary Material online). Two maximum likelihood (ML) analyses were performed using RAxML v8.2.4 (Stamatakis 2014), each using the best partitioning scheme found for the individual concatenations (supplementary table S2, Supplementary Material online). A GTRΓ model was assigned to each partition following the recommendations of the creator of RAxML (Stamatakis 2014). Each RAxML analyses consisted of three stages: 1) executing 100 ML inferences (50 GTRΓ, 50 GTRΓX) using 100 distinct randomized MP trees; 2) executing 1000 rapid bootstrap replicates (GTRΓ); and 3) mapping bootstrap values onto the best ML tree.

Ancestral State Reconstruction

Reconstructions were performed on two traits (morphologically discrete castes; PCB), using the 235 loci ML phylogeny. A parsimony unordered model and two ML models (Markov k-state one parameter, Mk1, and asymmetrical two-parameter Markov k-state, Assym.2) were implemented for ancestral state reconstruction within Mesquite version 3.2 (Maddison and Maddison 2017). The Mk1 model assumes a similar rate of evolution for gains (0 → 1) and reversals (1 → 0). The Assym.2 model estimates a separate rate for gains and reversals. For both traits, the Assym.2 model favored reversals as more likely, especially when more solitary vespid taxa were included. This is because as more solitary vespids are included, the more observed instances of no gains. However, we also performed analyses after excluding solitary taxa to check the robustness of the reconstructions to variation in taxon sampling. The sampling schemes implemented were: all vespids; only Zethinae, Polistinae and Vespinae; only Polistinae + Vespinae. Regarding the reconstruction of morphological castes, we also provide a summary of 10,000 stochastic character maps obtained by using the functions “make.simmap” and “densityMap” of the R package “phytools” (Revell 2013; supplementary figs. S4 and S5, Supplementary Material online), which mirror the results of the maximum likelihood reconstructions (table 1; supplementary fig. S3, Supplementary Material online).

All vespines were coded as having morphological castes, and thus PCB (Greene 1991; O’Donnell 1998). The only independent-founding polistines with documented evidence of morphological caste differences sampled in this study are Polistes olivaceous (Hunt 2007) and Polistes dominulus (Tibbetts 2006). Coding of morphological castes for Epiponini follows other studies (Noll et al. 2004; Noll and Wenzel 2008) and swarm-founding Ropalidia and Polybioides were coded as having morphological castes based on the literature (Yamane et al. 1983; Turillazzi et al. 1994; Fukuda, Kojima, Jeanne 2003). Coding species for state of PCB carries some ambiguity; there are limited studies directly testing for PCB in tropical independent-founding polistines. Temperate species have diapausing gynes that have ontogenetically established physiological differences from workers (Hunt et al. 2007). Diapausing females (i.e. future foundresses) have an increased likelihood of becoming a queen relative to nondiapausing females. Thus, diapause is a caste-biasing factor. All independent-founding polistines with seasonal nesting cycles (in temperate or subtropical regions) have PCB present. However, coding is complicated by the fact that many paper wasps have distributions that span tropical and subtropical/temperate regions, resulting in some populations having seasonal nesting cycles and others with aseasonal nesting cycles (Gadagkar 1991). Polistes canadensis is distributed in both temperate and tropical regions (West-Eberhard 1969), and Polistes carnifex is too (Carpenter 1996). Belonogaster is primarily distributed in subSaharan Africa and most species presumably have aseasonal nesting cycles. However, the exemplars of Belonogaster (and Ropalidia flavoviridis) in this study are from Madagascar, which has distinct dry/wet seasons. A seasonal nesting cycle and queen/worker dimorphism has been shown in Ropalidia from Madagascar (Wenzel 1992), as well as Belonogaster in South Africa (Keeping 2002), which suggests Belonogaster and Ropalidia in Madagascar likely have seasonal nesting cycles. However, all these ambiguous cases were coded as having PCB absent (state = 0) following the prevailing idea that tropical independent-founding polistines lack PCB, which ultimately biases the ancestral state reconstruction to support the current view that PCB was absent at the shared ancestor of Polistinae + Vespinae. Another analysis was done with these ambiguous cases coded as unknown (state =?) for comparative purposes. Swarm-founding epiponines that lack morphological castes were coded as having an absence of PCB (West-Eberhard 1978b; Felippotti et al. 2010). Polistes nimpha inhabits temperate regions and Parapolybia indica are known to have seasonal nesting cycles with hibernating gynes (Gadagkar 1991; Carpenter 1996), and so were coded as having PCB. For further justifications and references related to how we coded character states, see supplementary material data set S3, Supplementary Material online.

Author Contributions

P.K.P., J.M.C., and B.J.S. conceived and designed the study. P.K.P. and J.M.C. identified and selected the taxa studied. A.R.L. designed the probes. A.R.L., E.M.L., and P.K.P. generated the data. P.K.P. performed the phylogenetic analyses. P.K.P. and B.J.S. performed ancestral state reconstructions. P.K.P. wrote the paper with contributions from all other authors. The authors declare no competing financial interests.

Click here for additional data file.