The evolution of parental care strategies in subsocial wasps

Insect parental care strategies are particularly diverse, and prolonged association between parents and offspring may be a key precursor to the evolution of complex social traits. Macroevolutionary patterns remain obscure, however, due to the few rigorous phylogenetic analyses. The subsocial sphecid wasps are a useful group in which to study parental care because of the diverse range of strategies they exhibit. These strategies range from placing a single prey item in a pre-existing cavity to mass provisioning a pre-built nest, through to complex progressive provisioning where a female feeds larvae in different nests simultaneously as they grow. We show that this diversity stems from multiple independent transitions between states. The strategies we focus on were previously thought of in terms of a stepping-stone model in which complexity increases during evolution, ending with progressive provisioning which is a likely precursor to eusociality. We find that evolution has not always followed this model: reverse transitions are common, and the ancestral state is the most flexible rather than the simplest strategy. Progressive provisioning has evolved several times independently, but transitions away from it appear rare. We discuss the possibility that ancestral plasticity has played a role in the evolution of extended parental care. Parental care behaviour leads to prolonged associations between parents and offspring, which is thought to drive the evolution of social living. Despite the importance of insect parental care for shaping the evolution of sociality, relatively few studies have attempted to reconstruct how different strategies evolve in the insects. In this study, we use phylogenetic methods to reconstruct the evolution of the diverse parental care strategies exhibited by the subsocial digger wasps (Sphecidae). Contrary to expectations, we show that parental care in this group has not increased in complexity over evolutionary time. We find that the ancestral state is not the simplest, but may be the most flexible strategy. We suggest that this flexible ancestral strategy may have allowed rapid response to changing environmental conditions which might explain the diversity in parental care strategies that we see in the digger wasps today.


Introduction
Parents often invest resources to increase the fitness of their offspring, at a direct cost to their own future reproduction. When this investment continues after birth or hatching, a prolonged period of interaction between parents and offspring can lead to the evolution of social traits, such as traits associated with parental care, group living and eusociality. In insects, the evolution of extended parental care has been attributed to a variety of ecological pressures, including harsh environments, defence against natural enemies and ephemeral food resources when offspring are immobile (Field and Brace 2004;Gardner and Smiseth 2010;Wong et al. 2013). However, there have been relatively few phylogenetically based analyses of insect parental care (Wong et al. 2013).

Communicated by L. Keller
Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00265-020-02853-w) contains supplementary material, which is available to authorized users. Such analyses can be revealing. First, they can pinpoint when parental care transitions most probably occurred and in what direction. For example, phylogenetic work has revealed frequent losses of eusociality in bees, and that supposedly 'ancestral' strategies are in fact 'derived' (e.g. Wcislo and Danforth 1997;Schwarz et al. 2003). Second, phylogenetically based comparative analyses can shed light on the factors associated with evolutionary origins and losses of parental care traits (Gilbert and Manica 2010;Gittleman 1981;Székely and Reynolds 1995;Gonzalez-Voyer et al. 2008). For instance, in frogs, increased investment in parental care is associated with a switch to terrestrial living (Vági et al. 2019), and in Old world babblers, ground-nesting species have transitioned to building domed nests which offer greater protection (Hall et al. 2015). While such studies do not always allow inferences to be made about the direction of causality (Wong et al. 2013), they can help to guide more targeted empirical work.
The nesting aculeate (stinging) bees and wasps (Hymenoptera) are a promising group in which to study the evolution of parental care. They exhibit a diverse range of often elaborate parental care strategies, some of which are thought to act as precursors to the evolution of eusociality. In subsocial aculeates, mothers provision each offspring with arthropod prey (wasps) or pollen and nectar (bees). Evans and West-Eberhard (1970) identified a series of distinct strategies in terms of the sequence of care events and the degree of association between mothers and immature offspring that occurs during nesting. Four critical steps in the series are as follows. The simplest strategy, seen only in wasps, we refer to as 'Prey First': prey capture precedes nest preparation (state 1; Fig. 1b). A single large prey item is captured and either placed in a pre-existing cavity which serves as the nest, or cached temporarily while a nest is constructed. After laying a single egg, the nest containing egg plus prey is sealed permanently. Prey First strategies probably require the least maternal investment. The vast majority of wasps and all bees, however, use 'Nest First' strategies: they construct a nest (or offspring cell within a nest) before prey capture (Fig. 1b). Following construction, a single large prey (state 2: Nest First/ Single Prey) or multiple smaller prey captured one at a time (state 3: Nest First/Multiple Prey) are transported to the nest, and an egg laid on one of the prey items. Nest First sequences are likely to be energetically more costly because provisions may have to be transported long distances if the nest is built first instead of building it close to the eventual site of prey capture. In addition, the nest's location must be recalled via mechanisms such as orientation flights, presumably requiring costly sensory and neural machinery (Niven and Laughlin 2008).
Species that adopt the three strategies mentioned so far are termed mass provisioners: before the egg hatches, the nest (or cell) is permanently sealed containing the full mass of provisions required for offspring maturation. Contact between mothers and immature offspring is therefore minimized. Under the fourth strategy, progressive provisioning (state 4; Fig. 1b), however, the mother feeds each larva gradually as it grows, fostering larva-adult communication and resulting in a longer period of offspring dependency that can favour group living (Field and Brace 2004;Field 2005). Under progressive provisioning, the mother interacts repeatedly with her offspring and must gauge its needs as it develops. Neural costs associated with progressive provisioning are therefore likely to be greater than under mass provisioning. Progressive provisioning can be further complicated when mothers divide their attention between multiple offspring developing simultaneously, instead of provisioning successive offspring one at a time (Field 2005). In some of these species, different offspring are in spatially separate nests, whose locations must therefore be held in memory at the same time.
The sequence of evolutionary transitions between these four strategies is traditionally thought of in terms of a stepping-stone model as follows: Prey First ➔ Nest First/ Single Prey ➔ Nest First/Multiple Prey ➔ Progressive Provisioning (see Evans and West-Eberhard 1970). Prey First is thought of as ancestral because it is most similar to the ancestral parasitoid lifestyle, where a female lays eggs on a single prey (host) item which she then leaves where she finds it in the environment (Evans and West-Eberhard 1970). The order of subsequent transitions is based on logical interdependence. It seems likely that the evolution of Nest First/Single Prey is a preadaptation necessary for the evolution of Nest First/Multiple Prey: provisioning multiple prey would probably not be viable if the nest was constructed only after prey capture, since all but the last prey item would have to be left undefended in the environment while others were captured. In turn, Nest First/Multiple Prey seems likely to be a preadaptation necessary for the evolution of progressive provisioning, since progressive provisioners by definition feed their offspring gradually using multiple food items.
Although successive strategies in the above evolutionary sequence may be associated with increasing energetic costs, they also provide increasing potential benefits through protection against natural enemies (Evans and West-Eberhard 1970;Field 2005). Unlike the Prey First sequence, Nest First does not involve mothers leaving prey undefended while they construct the nest, reducing prey losses to scavengers, and reducing exposure to parasites that will later compromise offspring development (either directly through consumption of eggs or larvae, or indirectly through competition over food; Rosenheim 1989Rosenheim , 1990Field 1993). Progressive provisioning can provide further advantages, for example by enabling mothers to abandon offspring that fail part-way through development, rather than wasting a full quota of investment on them (Field and Brace 2004;Field 2005). However, although previous discussions have emphasized these ecological advantages, present day exemplars demonstrate that all four strategies are stable endpoints (Evans and West-Eberhard 1970). Furthermore, the neural costs likely to accompany more derived strategies make it conceivable that transitions in the opposite direction ('reversals') could occur (Evans 1966;O'Neill 2001;Shreeves and Field 2008). While the stepping-stone model is logical and intuitive, phylogenetic reconstructions of other behavioural traits associated with  (unshaded) in the context of relationships with other major clades (grey shaded) recovered by Branstetter et al. (2017) and Peters et al. (2017). Colours represent the parental care strategies reported for the majority of species in each group, which are illustrated in b. In b, dashed lines ending in arrows indicate repeated prey capture for strategies that involve providing multiple prey per offspring. parental care have reveal interesting and unexpected results, such as the frequent losses as well as gains reconstructed for eusociality (Wcislo and Danforth 1997;Schwarz et al. 2003). Combining inferences from phylogenetic work with empirical studies also suggests that the the diversity of social systems seen in the Hymenoptera may stem from plasticity in facultatively eusocial ancestors (Jones et al. 2017).
In this study, we test the stepping-stone model in subsocial wasps by viewing parental care strategies through a phylogenetic lens. We use a recent phylogeny to investigate parental care transitions in sphecid wasps, focussing particularly on the subfamily Ammophilinae. We use ancestral state reconstructions and a reversible-jump approach to infer the sequence in which parental care strategies have evolved, with a focus on three questions: (1) have transitions followed the stepping-stone model and proceeded in the 'forwards' direction (without any reversals) implied by previous discussions?
(2) what is the ancestral state for Sphecidae as a whole, and is it the 'simplest' state (Prey First)? (3) how frequently have different transitions occurred, and where have they occurred on the phylogeny? In the light of our results, we consider the evolutionary pathways and selection pressures that may have contributed to the strategies observed, and discuss the possible importance of behavioural plasticity in the evolution of parental care.

Study system
Sphecid and especially ammophiline wasps are an excellent group in which to investigate evolutionary transitions in parental care. Parental care has been characterised across a range of species, for which a well-resolved molecular phylogeny is available (Field et al. 2011). There are no eusocial ammophilines, but they otherwise represent a microcosm of strategies seen in aculeates more broadly, encompassing all four of the behavioural sequences and degrees of motheroffspring association described above (Evans 1959;Field et al. 2011;see Fig. 1a). The progressively provisioning ammophilines that have been studied all maintain more than one nest simultaneously (Baerends 1941;Evans 1965;Kazenas 1970;Weaving 1989a, b;Hager and Kurczewski 1986;Field et al. 2018).
Ammophilinae is a monophyletic subfamily of apoid wasps in the family Sphecidae s.l. (Melo 1999;Field et al. 2011). The subfamily comprises approximately 300 described species in 6 genera: Ammophila, Eremnophila, Eremochares, Hoplammophila, Parapsammophila and Podalonia. Ammophilines are relatively homogeneous in terms of gross morphology and basic nesting biology. They are relatively large wasps with long, thin abdomens. Nests of nearly all species are short burrows, usually dug in sandy soil, with each nest containing a single offspring. Mothers typically provisioning offspring with lepidopteran caterpillars which are paralysed using the sting, and most species appear to be prey generalists within Lepidoptera (e.g. Evans 1959;Weaving 1989a, b;Field 1992b). The wasp larva feeds on the prey provided by its mother, then pupates in the nest. Females are relatively large compared with males in species that capture relatively larger prey, probably reflecting the energetic demands of prey carriage (Field et al. 2015). The small genus Hoplammophila (4 species) is exceptional in that females make multicellular nests in pre-existing cavities such as old beetle borings in dead wood. In common with most aculeates, ammophilines are attacked by various heterospecific cuckoo parasites and parasitoids (e.g. Rosenheim 1989;Field and Brace 2004). Facultative intraspecific parasitism is also widespread, with prey commonly being stolen or eggs replaced by conspecifics (Field 1989a, b;Field et al. 2018).
Three recent phylogenies for the wider Apoidea based on different molecular datasets nest a monophyletic Sphecidae within other clades containing largely species that use the Nest First/Multiple Prey strategy (Peters et al. 2017;Branstetter et al. 2017;Sann et al. 2018). However, Sphecidae and especially ammophilines themselves exhibit a diverse range of strategies ( Fig. 1a, b). We focus primarily on ammophilines in this study, but we include representatives from 9 of the 13 non-ammophiline genera of sphecids, and 2 non-sphecid apoids from the large clade identified as the sister-group to Sphecidae (Tachysphex and Anacrabro: Fig.  1a; Peters et al. 2017). Recent systematic work means that the nomenclature of suprageneric apoid groupings is somewhat in flux (e.g. Ohl 1996; Field et al. 2011;Peters et al. 2017;Branstetter et al. 2017;Sann et al. 2018). To avoid confusion, we use the subfamily Ammophilinae (equivalent to Bohart and Menke's (1976) tribe Ammophilini) to encompass ammophiline wasps as a subfamily within the monophyletic Sphecidae (equivalent to Bohart and Menke's Sphecinae). Sphecidae contains two other principle subfamilies, Sceliphrinae and Sphecinae.

Phylogeny
We used a published molecular phylogeny for 40 ammophilines, which additionally included 9 nonammophiline genera within Sphecidae s.l. (Field et al. 2011). The maximum credibility (MaxC) Bayesian tree was based on a 2229 base pair DNA dataset combining three gene sequences: Mitochondrial COI, EF-1α and Opsin (Field et al. 2011). We carried out ancestral state reconstructions on this tree, using Maximum likelihood and Bayesian methods (see below). While the MaxC phylogeny is well-resolved, with the majority of nodes present in most trees (posterior probabilities generally between 0.8 and 0.9), we conducted additional Bayesian ancestral state reconstructions using reversible-jump MCMC to simultaneously take into account uncertainty in the phylogeny as well as model uncertainty.

Parental care strategies
Information on parental care strategies was extracted from the literature. Bohart and Menke (1976) provide a thorough and well-balanced description of the literature prior to 1976 for each genus, as well as information for many individual species. We supplemented this using additional papers known to the first author through more than 30 years of experience working on sphecids. In addition, we searched for publications concerning each genus in the Web of Science using each of the 19 sphecid genus names as a search term. Work published since 1976 provided new information about individual ammophilines species (Table S1), while adding little to Bohart and Menke (1976)'s genus-level summaries for other sphecids. We were able to assign states to 37 species in the phylogeny: 26 of the 40 ammophilines and 11 nonammophilines, including representatives of 9 of the 13 genera from the other sub-families within the Sphecidae and 2 outgroup taxa from the sister lineage of 'Crabronidae' (Fig.  1a, Table S1). All four parental care strategies described in the "Introduction" were reported (see Fig. 1b). Parental care in each non-ammophiline sphecid genus was categorised according to the majority strategy observed among its species (which was the only strategy in 6 out of 9 genera; Bohart and Menke 1976; Table S1). As with most phylogenetic comparative studies, sampling was necessarily limited to species where both phylogenetic position and focal parental care traits are known. Nevertheless, our sample is generally a good reflection of the known diversity of species and their parental care strategies. The number of ammophiline species per genus in our dataset is approximately in proportion to the number of described species (Bohart and Menke 1976), with representation of the genus Podalonia, and therefore of strategy 1 (Prey-First) which is typical of Podalonia species, being slightly low. In the genus Ammophila, 7 species groups are discussed by Bohart and Menke (1976), and our sample includes representatives of 6 of these, the exception being the small induta group. The genus Parapsammophila is not included in our sample because its biology has not yet been described.

Ancestral state reconstructions
We ran ancestral character reconstructions using maximum likelihood and Bayesian methods in the packages ace and Simmap respectively in R (R core team 2019; Paradis and Schliep 2018;Revell 2012), using the maximum credibility tree from Field et al. (2011). We repeated reconstructions using reversible-jump MCMC to estimate ancestral states across a posterior sample of 3000 phylogenies which comprised the Bayesian estimated maximum credibility tree in Field et al. (2011). We used rjMCMC in addition to other methods because it searches amongst all possible models over a sample of possible trees, visiting each model in proportion to its posterior probability, thus taking both phylogenetic and model uncertainty into account (Pagel and Meade 2006). We used the package Multistate in Bayestraits (Pagel and Meade 2007) to implement rjMCMC. We ran MCMC chains for 100,000,000 iterations with a burn-in of 50,000 generations and sampling every 50,000 generations. Parameter means and ranges from maximum likelihood models helped to inform our choice of priors. After checking the effects of different priors, we selected an exponential prior with a uniform hyperprior ranging between 0 and 1 for all parameters.
The two non-sphecid outgroups, Tachysphex and Anacrabro, are effectively placeholders for the speciose 'crabronid' clade that is the sister group to Sphecidae (> 3000 species: Fig. 1a). Tachysphex and Anacrabro were included to more reliably estimate the ancestral state in the Sphecidae. While the vast majority of species in the 'crabronid' clade are reported as state 3 (Nest First/Multiple Prey) in the literature (e.g. Bohart and Menke 1976), there are sporadic instances of other parental care strategies, in particular state 2 (Nest First/Single Prey). We therefore ran models twice, once where Tachysphex was scored as Nest First/Single Prey (state 2) and once where it was scored as Nest First/ Multiple Prey (state 3). Altering the coding did not change the overall results, in particular, regarding the ancestral state for Sphecidae, and we therefore present the results of the more conservative analysis where Tachysphex was scored as state 2.
In order to check the robustness of our findings to phylogenetic incompleteness, we attempted to allow for sphecid taxa missing from the maximum credibility tree of Field et al. (2011). To achieve this, we added parental care data for an additional 27 notional species in 8 genera (see Table S1). In the resulting 'MaxC+' dataset, parental care strategies and numbers of notional species were added approximately in proportion to their actual occurrence in each sphecid genus (Bohart and Menke 1976). We repeated ML and Bayesian ancestral state reconstructions on a tree that included these additional species (MaxC+ dataset). Additional species were included with congeners as polytomies in the tree, with missing genera assigned positions based on morphology-based taxonomy (Bohart and Menke 1976). Branch length for polytomies was minimized, and set as 10 -20 based on minimum branch sizes throughout the tree. We additionally ran Simmap on 1000 of these trees for which the polytomies were resolved at random to account for uncertainty in the topology associated with this method.

Phylogenetic signal
We calculated measures of phylogenetic signal: the tendency for more closely related species to resemble each other in terms of their trait values. Revell et al. (2008) caution against using phylogenetic signal alone to estimate evolutionary processes, but it does provide information about the extent to which related species are similar (Blomberg et al. 2003). As such, we use estimates of phylogenetic signal to help interpret the robustness of results for ancestral state reconstructions repeated on the three different datasets (Table S1): the maximum credibility tree of Field et al. (2011) (MaxC); the tree with extra species treated as polytomies (MaxC+); and the tree with extra species with polytomies randomly resolved over 1000 trees (MaxC+ RRP). We calculated two measures of phylogenetic signal. The first was Blomberg's K (Blomberg et al. 2003), which is a quantitative estimate of phylogenetic signal with an expected value of 1 if the trait has evolved under Brownian motion. As Blomberg's K was developed for continuous traits, we also calculated δ, a measure of phylogenetic signal for categorical traits (Borges et al. 2019). We used the R package 'phytools' (Revell 2012)

Estimating rates of state change
If transitions can occur in either direction, there are 12 possible transitions between the four parental care states represented in our dataset. However, because of the logical interdependence between states discussed in the "Introduction", we focussed a priori on transitions in the stepping-stone sequence 1 ↔ 2 ↔ 3 ↔ 4. We refer to transitions from lower to higher states in this sequence as 'forward' transitions, and from higher to lower states as 'reverse' transitions. In order to test specific hypotheses about rates of change between states, we ran multiple models that allowed different transition rate parameters to vary (see SI Figure S1 and Table S2). There were insufficient species in the phylogeny to reliably estimate all rate parameters simultaneously, and so we used models that restricted certain transition rates either to zero (i.e. the transition cannot occur) or to be equal to other rates.
The number of rate parameters that each model estimated varied from one to three. Using ace and Simmap, we initially tested whether forward and reverse transitions as a whole occur at different rates, by testing whether a one-rate model that restricted all transitions to occur at the same rate was a worse fit than a two-rate model that allowed reverse and forwards transitions to differ. Model comparisons were based on the likelihood ratio test (LRT) and on Bayes Factor (BF) for ML and Bayesian models respectively. We also tested whether models that allowed 'jumps' as well as single step transitions (jumps are transitions that skip intermediate states, such as from state 1 to 3, or 2 to 4) fit the data better than models that allow only single steps (such as from 1 to 2). We went on to use increasingly complex models to test more specific hypotheses (see Figure S1, Table S2). For instance, by allowing or restricting particular rate parameters, we tested whether reverse and forward jumps occur at different rates; whether different sized jumps (1➔3 or 1➔4) occur at different rates; and whether any specific transitions do not occur (see Figure S1 for a visual summary of all models). In order to assess how robust the results were to missing species and uncertainty in the tree topology, we compared best fitting models using the three different datasets: including or excluding the additional 27 species (MaxC/MaxC+ datasets) and with polytomies in MaxC+ resolved randomly (MaxC+ RRP).
We carried out additional ancestral state reconstructions using rjMCMC by restricting certain groups of transitions to be equal (SI : Table S3). In this case, we did not restrict any rates to zero (as we did using Simmap and ace, above), since the rjMCMC function in multistate itself estimates the probability that certain transitions/groups of transitions are zero. This allows the model to determine, for example, whether forward rates as a group occur at the same rate as reverse rates as a group, and whether either rate is zero.

Ancestral state reconstructions
State at the root Over all models, with (MaxC+) or without (MaxC) the additional 27 species, Nest First/Multiple Provisioning (state 3) was consistently reconstructed as the most likely state at the root ( Fig. 2; see SI Tables S3 and S4 for likelihoods and posterior probabilities). The posterior probability at the root tended to be lower when phylogenetic and model uncertainty was taken into account in the rjMCMC models (SI Tables S3 and S4), particularly when Tachysphex was coded as 2 (Nest First/Single Prey). In general however, altering the coding of Tachysphex from a multiple to single prey did not have a substantial impact on our findings. The proportion of time spent in state 3 (Nest First/Multiple Prey) across the tree was reconstructed as longer than the time spent in any other state across all Simmap models, suggesting that this is the most evolutionarily stable state (SI Table S4).

Estimating rates of state change
Full details are given in the following two paragraphs but overall, our results suggest that reversals as well as forward transitions between logically adjacent parental care states have occurred during evolution, and that if anything, the reverse rate is higher than the forward rate. Our results also suggest that forward 'jumps' that miss out intermediate states may occur, but that reverse 'jumps' are unlikely. The estimated rate for forward jumps was similar to the rate for forward single steps using ace and Simmap, but the estimated rate for single steps was higher using rjMCMC.
A robust finding was that models that allowed reversals and forward single steps to occur at different rates performed consistently better than models which forced the two rates to be the same. Additionally, models allowing forward jumps performed better than models that restricted them to zero, but the opposite was true for reverse jumps: models that permitted reverse jumps performed poorly compared to when they were restricted to zero. Maximum likelihood (ace) and Bayesian (Simmap) analyses on the single maximum credibility trees (MaxC and MaxC+) gave reasonably consistent results (SI Table S2). Among models that specified different rates for different kinds of transitions, the 2-rate model where reverse jumps were restricted to zero, forward steps and jumps were equal and occurred at a lower rate than reverse steps performed consistently well based on the LRT and BF scores (SI Table S2; Fig. 3a). This model was the best fit to the data when only species included in the original phylogeny (MaxC: SI Table S2) were considered, but also performed well when extra species were included as polytomies (MaxC+) in ML models (Table S2). For Bayesian analyses (Simmap), this model was the best fit when only species present in the original phylogeny (MaxC; BF compared to null 'all rates equal model' = 5.02) were included, but not when the 27 additional species were added as polytomies (MaxC+ ; Table S2). For MaxC+, the best fitting model was slightly different: forward but not reverse jumps were allowed as before, but jumps occurred at a lower rate than forward single steps. However, this model was not a substantial improvement on the null model where all rates were equal (BF = 1.76).
Only in the extreme case where polytomies were randomly resolved (MaxC+ RRP) were results less consistent: most models then performed poorly compared to the null (all rates equal) model (Table S2). Measures of phylogenetic signal suggest a possible explanation. Where only species in the original phylogeny were included (MaxC), and where additional species were included (MaxC+), trait distribution across the phylogeny was non-random (MaxC: K = 0.56, p = 0.04; δ = 2.51; MaxC+: K = 0.7, p = 0.01, δ = 2.59): species that were more closely related had similar trait va lue s. H o w ev er, w hen w e ra ndo mly r eso l ve d polytomies, the phylogenetic signal was broken up (MaxC +RRP: K = 0.31, p = 0.20, δ = 1.07 ± 0.21 SD). Given that our sampling was generally representative of species and trait diversity (see "Methods"), it seems most likely that trait distribution is actually non-random. Resolving polytomies randomly is likely to be overly conservative, because by artificially creating random, finescale patterns, it dilutes the signal from the observed broad-scale clustering in the trait distribution.  Figure S1) where reverse jumps were restricted to zero, forward steps and jumps occurred at equal rates and backward steps occurred at a different rate.
The rjMCMC models produced similar results to the ace and Simmap analyses. rjMCMC visits chains of rate parameters in proportion to their posterior probability and estimates rate parameters based on the frequency that each chain is landed on. In the current study, the chain that received the greatest support had reversals and forward single steps and forward jumps occurring at similar rates, but reverse jumps occurring at a modal rate of zero (Table S3; Fig. 3b). This differs only slightly from the results with Simmap and ace, where the rate for backward single steps was estimated as higher than that for forward single steps (Fig. 3a), but the general pattern of no reverse jumps, and other transitions occurring at similar rates, remained. The results were broadly consistent whether Tachysphex was scored as states 2 or 3, in terms of both the state reconstructed at the root and the transition rates (Tables S2, S3 and S4).

Estimating where transitions occur in the phylogeny
Based on the likelihoods (ML) and posterior probabilities (Bayesian) of each state at the nodes in the phylogeny (Fig.  2), all models were consistent in reconstructing four independent transitions to progressive provisioning (state 4; in extant taxa Eremochares dives, Ammophila pubescens, A. rubiginosa and A. azteca). These transitions were forward, arising from state 3 (Nest First/Multiple Provisioning), except that in some simulations, a jump from (2) to (4) was reconstructed for A. pubescens and E. dives (Fig. 2, SI Figure S2a and b).
Prey First (state 1) was consistently reconstructed as a derived trait requiring backward transitions and arising at the tips of the phylogeny. Transitions to Prey First occurred exclusively from state 2 (Nest First/Single Prey). Ancestral state reconstructions suggest that the Prey First strategy evolved between Fig. 3 Transition rates estimated from a Simmap, based on 1000 simulations across the MC tree of Field et al. (2011), with Tachysphex scored as state 2. Rates are proportional to line thickness, and are estimated for the best model, where reverse jumps are restricted to 0, forward single steps and jumps occur at equal rates (1.15) but backwards single steps occur at a higher rate (2.51). b Rates estimated using rjMCMC when Tachysphex was scored as 2. Rates are proportional to line thickness and are also shown as numbers adjacent to lines. Grey dashed lines represent rates with a mode = 0. For each transition, the median rate X̃± SD is shown, alongside a frequency distribution of rates estimated from 2000 runs over 3000 trees.
3 and 5 times independently, once in Prionyx kirbii, once or twice in the genus Podalonia (P. affinis/P. hirsuta) and once or twice in Ammophila marshii/A. wrightii. Transitions away from Prey First were reconstructed as less likely, reflecting its 'tippy' distribution (Fig. 2).
Transitions between states 2 and 3 (Nest First/Single Prey and Nest First/Multiple Prey) appeared to be the most frequent. They were the most common transitions in ML models (Table S5) and in both Simmap and ace, transitions to state 2 were exclusively from state 3. The transition from 3 to 2 almost certainly occurred once prior to the divergence of Eremochares and Podalonia. The reverse transition back to 3 then occurs after Ammophila pictipennis and A. nigricans diverge.

Discussion
Sphecid wasps exhibit four major parental care strategies. We found that phylogenetic signal for these strategies was moderate but significant: more closely related species exhibit similar strategies, reflecting phylogenetic relatedness and/or shared environments during evolution. In particular, groups of Nest First/Single Prey and Nest First/Multiple Prey taxa are clumped together (Fig. 2), with Prey First and Progressive Provisioning sprinkled among them. The signal is likely to be even stronger once more species in the genus Podalonia can be included, since most Podalonia probably use the Prey First strategy (Bohart and Menke 1976). The scattered distribution for Progressive Provisioning is not confined to Sphecidae: a similar pattern occurs in nyssonine wasps (Evans 1966).
The four different strategies form a logically connected stepping-stone sequence, from the parasitoid-like Prey First strategy through to Nest First provisioning of a single prey, then multiple prey and finally progressive provisioning (see Evans and West-Eberhard 1970). Our phylogenetically based analysis, however, suggests that parental care evolution has in fact often not progressed 'forward' through this sequence. Instead, transitions in both directions seem to have occurred, with some reconstructions suggesting that reverse transitions are the more frequent. For example, Nest First/Multiple Prey (state 3) is the likely ancestral state for Sphecidae as a whole, 'reversing' to Nest First/ Single Prey (state 2) at the point where Prionyx and the Ammophilinae diverge (Fig. 2). Similarly, rather than being the starting point from which other strategies evolved, the Prey First strategy (state 1) is a derived state, with transitions to it occurring exclusively on the terminal branches of the phylogeny from Nest First/Single Prey (four times).
While our sampling is generally representative of the known trait and species diversity in the group (see "Methods"), it is necessary to be cautious in interpreting our results, given that many taxa remain unstudied. We therefore discuss our key findings in the light of qualitative patterns seen across the group more broadly, in order to assess their robustness and significance. For instance, while we could not include 4 of the 19 sphecid genera, Trigonopsis, Dynatus, Palmodes and Chilosphex, current phylogenetic hypotheses suggest that their future inclusion is unlikely to alter our main conclusions (Bohart and Menke 1976;Sann et al. 2018). Dynatus and Trigonopsis are in the same clade as Podium and Penepodium (Fig. 2), and all species are reported as Nest First/ Multiple Prey (state 3), typical of the relatively basal subfamily Sceliphrinae (Figs. 1 and 2), except that one species of Trigonopsis may be a progressive provisioner (Eberhard 1972;Bohart and Menke 1976;Kimsey 1978). Palmodes (Nest First/Single Prey) and Chilosphex (Nest First/Multiple Prey) are closely related phylogenetically to Prionyx (see Fig. 2; Bohart and Menke 1976;Sann et al. 2018). Inclusion of these four genera would thus bolster our conclusion that progressive provisioning evolves sporadically ( Fig. 2 and see below), and is unlikely to alter broad conclusions such as the reconstructed ancestral state for Sphecidae and the occurrence of reverse transitions. On the other hand, the forward jumps that our reconstructions imply, omitting intermediate states, should be interpreted with caution, since it is possible that denser taxon sampling will reveal intermediate states.
A noteworthy result is that provisioning of a pre-built nest with multiple prey (Nest First/Multiple Prey) is ancestral for sphecids. We are confident in the robustness of this finding, given that the vast majority of species in the Crabronidae sister clade (Fig. 2) are Nest First/Multiple Prey (Bohart and Menke 1976). Indeed, preliminary analyses using a broader scale phylogeny of the Apoidea (Sann et al. 2018) are consistent with this (work in progress). Nest First/Multiple Prey is probably the most flexible parental care strategy exhibited by aculeates, allowing females to exploit a wider range of prey sizes (since each individual item need not be enough to fully provision an offspring), and to tailor their provisioning to short-term changes in environmental conditions. For example, mothers can provide nests with more or fewer prey, depending on prey size availability, with some Ammophila species known to switch routinely between use of multiple or single prey (Table S1; e.g. Field 1992a). Mothers can thus more easily adjust total provision mass according to offspring sex (females are always larger than males ;Field 1992a;Field et al. 2015). It may be that this ancestral flexibility contributes to the evolutionary pathways that follow and to the diversity of parental care strategies seen in ammophilines. The 'flexible stem hypothesis' suggests that plasticity in an ancestral lineage can result in replicate descendant lineages in which one of the once environmentally induced phenotypes has become fixed ('genetic assimilation';West-Eberhard 2003;Wund et al. 2008). A similar process has been suggested to occur during evolutionary pathways to eusociality, where phenotypic plasticity exhibited by facultatively eusocial species is the substrate for selection (Jones et al. 2017). When selection is consistent, traits are thought to become canalised, environmental sensitivity is lost and obligate eusociality and caste dimorphism arise (Jones et al. 2017).
Another robust result is that a forwards-only stepping-stone model is not a good fit to the evolution of parental care in the Sphecidae. Such a model seems intuitive but anthropomorphism, emphasizing the evolution of cognitive complexity, may be part of its appeal (Dacey 2017). Parental care evolution in the Sphecidae, however, probably relates more to ecological complexity and ecological variability than to selection for cognitive complexity per se. While we expect evolution to progress 'forwards' in the sense that sphecid hosts co-evolve with their inter-and intraspecific parasites, the benefits of more complex parental care are likely to be highly variable in both time and space, reflecting variation in prey availability and parasite pressure. Our phylogenetic reconstructions suggest that building a nest before prey capture arose prior to the divergence of Sphecidae and Crabronidae. Ancestral Nest First provisioning may have provided resilience to changes in parasitism risk, with flexibility inherent in the Nest First/ Multiple Prey strategy also being particularly robust to changes in prey type and availability. Although our results do not support a forwards-only stepping-stone model, logical interdependence between states is probably still important. For example, our reconstructions suggest that Prey First (state 1) arose only through 'reversals' from the state adjacent in the sequence, Nest First/Single Prey (state 2), perhaps when parasite pressure was relaxed.
The transition from a Nest First to a Prey First nesting sequence could represent heterochrony, i.e. a simple rearrangement of the timing of events during nesting. In contrast, the transition to progressive provisioning is more profound. Progressive provisioners exhibit distinctive behaviours that are unique among ammophilines (Baerends 1941;Evans 1965;Kazenas 1970;Hager and Kurczewski 1986;Weaving 1989a, b;Field et al. 2018). They simultaneously maintain multiple spatially separated nests, rather than provisioning nests sequentially one after the other, necessitating the scheduling of parental care among different offsprings. Females also wait for 2 or more days after egg-laying, before inspecting inside each nest, adding more prey only if the inspection reveals that the offspring is alive and well (Field and Brace 2004). This provides a major advantage under high levels of parasitism, because it enables inspecting mothers to terminate investment in nests containing parasites that have destroyed their part-grown offspring, and because it reduces exposure to parasites until wasp larvae have become immune (Field and Brace 2004;Field 2005). Progressive provisioning and the Prey-First strategy nevertheless exhibit similar phylogenetic patterns. Both appear to have evolved multiple times independently, but transitions away from either are unlikely. This may indicate that these strategies tend to become more canalised than ancestral Nest First strategies. In the early stages of transitions to Prey First or progressive provisioning, flexible expression of the ancestral state may occur. Indeed, occasional mass provisioning occurs in the progressive provisioner A. pubescens (J. Field and W.A. Foster, unpublished data), and both Nest First and Prey First provisioning occur within a few species of Podalonia (Field 1993 and references therein). Under consistent selection, however, for instance when rates of parasitism are extremely high or low, the costs of plasticity may become prohibitive, such that progressive provisioning or Prey First strategies becomes canalised (Snell-Rood 2013;Niven and Laughlin 2008;West-Eberhard 2003). If plasticity facilitates evolutionary transitions, this loss of plasticity could produce the 'tippy' phylogenetic distribution we observe, where progressive provisioning and prey-first are derived states, with transitions away from them being rare.
The diversity of parental care strategies seen in sphecid wasps stems from multiple independent transitions between different states. We have shown that transitions do not always follow a traditional 'forwards' stepping-stone model that correlates with increasing cognitive complexity. The sequence of evolutionary transitions suggests that more cognitively complex strategies, with associated costs of neural architecture, are not always derived. Our analysis instead reveals backwards transitions, and perhaps even jumps between non-adjacent states. When viewed together with the ancestral state for sphecids being mass provisioning with multiple prey, our results are consistent with a 'flexible stem' model of evolution of extended parental care in the digger wasps. We suggest that like facultative eusociality, mass provisioning with multiple prey is a flexible ancestral state that may lead to the same evolved solutions to similar problems arising repeatedly through genetic assimilation (Pfennig et al. 2010). We hope that our findings will stimulate further theoretical, empirical and comparative study into factors that drive the evolution of parental care in the Sphecidae. Particularly welcome will be work assessing the importance of phenotypic plasticity as a factor that acts alongside life history preadaptations and ecology in facilitating the evolution of extended parental care and perhaps also eusociality (Field 2005). We also encourage studies that incorporate phylogeography and spatial/ecological predictors, for example to assess the validity of our hypothesis that variation in parasitism rates and prey availability have driven the evolution (and perhaps canalisation) of Prey First and progressive provisioning.
Data accessibility All data will be archived on Dryad digital repository on acceptance.

Compliance with ethical standards
Conflict of interest/competing interests The authors declare that they have no conflict or competing interests.
Ethics approval and consent Data analysed for this study were compiled from the available literature, ethical approval and consent to participate and publish is not applicable.
Code accessibility Code used for this study will be archived on Dryad digital repository on acceptance.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.