A Nonstationary Markov Model Detects Directional Evolution in Hymenopteran Morphology

Implementation of a Nonstationary Model for Binary Characters

We implemented a nonstationary model of evolution for binary characters in MrBayes 3.2 (Ronquist et al. 2012b). The model assumes a single instantaneous rate matrix over the whole tree, and allows for the state frequencies at the root to differ from the stationary frequencies of the process (Boussau and Gouy 2006). The model thus accounts for nonstationarity, even though the process along the branches is homogeneous. Since the model is directional, probabilities have to be computed on rooted trees, but the position of the root does not have to be specified a priori. The additional parameters of the model, the root frequencies, are updated during the Markov chain Monte Carlo (MCMC) procedure in the same way as the stationary frequencies, using standard sliding window and Dirichlet proposal mechanisms. The model, which is currently only available for binary data (data type “restriction”), can be invoked in MrBayes by the command “lset statefreqmodel = directional” (instead of “stationary”, the default), and the prior on the root frequencies is set with the command “prset rootfreq = dirichlet/fixed” (see the “help” function in MrBayes for details).

To allow for direct model testing and to integrate over model uncertainty, we implemented a RJ algorithm, that is, a transdimensional MCMC (Green 1995, 2003) which jumps between the stationary and nonstationary models. We assume equal prior probabilities of the two models, and use a split-merge move as follows. In the split move, we obtain new stationary and root frequencies by two independent draws from a beta distribution centered on the stationary frequencies in the previous generation, with the variance of the distribution used as a tuning parameter. The merge move draws new stationary frequencies from a beta distribution with the mean equal to the average of the previous stationary and root frequencies, under the same tuning parameter as in the split move. This move has a Jacobian of 1, as calculated following the approach detailed in Green (2003) and confirmed by running the model without data, which returned the prior probabilities for the stationary and directional models. As the directional-evolution model requires a rooted tree, this is assumed throughout, even while the MCMC is sampling under the stationary model and there is thus no information about the placement of the root. When there is a low but significant probability for the directional model, the information about the position of the root is gradually lost while sampling from the stationary model, reducing the acceptance rate of swaps back to the directional model. Under such conditions, a rooting constraint can be introduced to improve mixing. Both the root frequencies and stationary frequencies are updated with sliding-window and Dirichlet proposals. The command in MrBayes to invoke the RJ approach is “lset statefreqmodel = mixed”. The implementations of the directional model and RJ algorithm were tested by running the MCMC without data, which recovered the prior distributions of state frequencies, model probabilities, and topology (results not shown).

Simulation Study

To study the performance of the model and test for circumstances under which directional evolution could be discerned based on data from only extant taxa, we performed simulations in R using the “ape” package, including the “matexpo” function (Paradis et al. 2004; R Development Core Team 2009), and a number of newly developed functions. All R scripts are available as Supplementary Material on Dryad at http://dx.doi.org/10.5061/dryad.0sn23. We use the term “branch” to specify a single edge of the phylogeny (in contrast to “lineage” which can denote a collection of edges) and “tree” to denote a topology with branch lengths.

We tested the influence of the shape of the tree (see below), the number of taxa sampled (8 or 16), the evolutionary rate in expected substitutions from the root to any of the tips (0.01 to 1), and the strength of directionality, that is, the difference between root and stationary frequencies (several combinations of the state frequencies 0.1, 0.3, 0.5, 0.7, and 0.9 were used, either the same at the root and at equilibrium to simulate stationarity, or different to mimic directional evolution). We used three tree shapes that represent extremes in terms of topology and branch lengths: (i) an entirely balanced tree with Grafen branch lengths (where the height of each node is equal to the number of subtended terminals minus one, divided by the root height; Grafen 1989) (tree A); (ii) an entirely unbalanced tree with Grafen branch lengths (tree B); and (iii) an entirely unbalanced tree with all branches of equal lengths (tree C) (Fig. 2). The trees were scaled so that trees A and B, whose branch lengths conform to a molecular clock, have a root-tip distance of one, while tree C was scaled so that the average root-tip distance was 1. On these trees, we simulated 1000 (in the case of trees with eight taxa) or 100 (for 16 taxa) data sets of 100 binary characters each. These data sets were analyzed in MrBayes under the stationary and directional models and using the RJ approach, running two independent runs with one cold and one heated chain for 100,000 generations, sampling every 100th generation. We assumed that only the variable characters could be observed (using the command “lset coding = variable”), thus effectively excluding the constant characters from the analyses. Such an ascertainment bias is well suited for analyzing morphological or restriction-site data which are the focus here, but is less appropriate for data such as DNA or amino acid sequences where constant characters are also observed. We employed R scripts to study the power to distinguish between directional and stationary evolution, the impact on topology inference and root recovery, and the accuracy of the inferred state frequencies at the root and at equilibrium. Simulations and analyses were run at the National Supercomputer Center at Linköping University in Sweden (NSC).

Empirical Example: Hymenopteran Morphology

We modified the morphological data set in Ronquist et al. (2012a) by recoding all wing vein, muscle, and sclerite characters so that they followed the same state labeling, with zeros for the state “absent” (or “fused” in muscles and sclerites) and ones for “present” (or “separate” for muscles and sclerites). All other characters were left unchanged. A list of characters is given in Supplementary Table S1, available on Dryad at http://dx.doi.org/10.5061/dryad.0sn23.

First, we analyzed the wing vein (23 characters), muscle (56 characters), and sclerite data (67 characters) in combination with the remaining morphological characters (208 characters), invoking a separate directional or stationary model for each of the binary partitions, under a “variable” coding bias and gamma-distributed among-character rate variation. We used a more restrictive and a more diffuse branch length prior, that is, an exponential distribution of mean 0.1 and of mean 1, respectively. The morphology-only analyses employed four independent runs of one cold and three heated chains each, sampling every 1000th of 10 million generations. Convergence was assessed through the potential scale reduction factor (PSRF, smaller than 1.01 unless mentioned otherwise) for scalar parameters and the average standard deviation of split frequencies (ASDSF, below 0.02) for the topology. To calculate Bayes factors for model comparison, we used the stepping-stone algorithm (Xie et al. 2011) with 49 steps of 1,000,000 generations each plus 1,000,000 generations as initial burn-in, and 250,000 generations burn-in within each step. The alpha parameter was set to 0.4, and sampling went from the posterior to the prior. All of the steps reached an ASDSF below 0.03, except for the last 10 steps (sampling close to the prior) in which no split reached the minimum frequency of 0.1.

We then added the molecular data from Ronquist et al. (2012a) to improve the estimates of topology and relative branch lengths, once more running the binary morphological partitions under the stationary and directional models using the RJ approach. The molecular data were partitioned into genes and first and second versus third codon positions in the protein-coding genes, and a separate RJ algorithm over the GTR-model subspace (Huelsenbeck et al. 2004) was used for each partition, combined with discrete gamma-distributed among-site rate variation and a proportion of invariant sites. We used an exponential distribution with rate 10 (mean 0.1) as the branch-length prior. Settings for the MCMC search and convergence diagnostics were the same as in the case of the morphology-only analyses. Model testing was performed using the RJ algorithm, but no stepping-stone approach, given the slow convergence close to the prior observed for the morphology-only analyses and the large computational effort associated with this algorithm.

Finally, we repeated the total-evidence dating analysis as detailed in Ronquist et al. (2012a), including fossil morphology and fossil ages in the data matrix, with the updated morphological data set and under the stationary and directional models as well as with the RJ approach. Topology convergence was lower than in the nonclock analyses (ASDSF of 0.047, 0.055, and 0.035, respectively, for the directional, stationary, and RJ analyses), and the clock-rate estimates showed some deviations between runs as well (PSRF of 1.07, 1.04, and 1.02, respectively). All analyses were run at the National Supercomputer Center at Linköping University in Sweden (NSC). Alignments and trees are deposited on TreeBase (www.treebase.org, study number S16060).

Recovery of the Directional Model in Simulated Data Sets

Correct recovery of the directional model depends on several factors. First of all, the evolutionary rate has a major impact on correct inference (Figs. 2 and 3, see also Supplementary Files S3, available on Dryad at http://dx.doi.org/10.5061/dryad.0sn23). While very low rates (0.01–0.05 expected substitutions between root and tip) result in too little information for distinguishing between the models, rates of 0.5 and especially 1 seem already too saturated in many cases. An evolutionary rate of 0.1 or 0.2, sometimes 0.5, retains most information about the process and the tree.

The shape of the tree also plays a major role. Detecting directional evolution is most difficult on a balanced tree (tree A; Fig. 2a,b), followed by the unbalanced tree with ultrametric branch lengths (tree B; Fig. 2c,d). The unbalanced tree with branches of equal length (tree C; Fig. 2e,f) is the most favorable for inferring directional evolution. In all cases, fixing the topology to the true tree for the analysis has a positive effect. Increasing the size of the tree from 8 to 16 taxa also increases the success in correctly recovering the directional model (Fig. 3a,b), especially if directionality is strong. The proportion of analyses that recovers the correct topology is lowest for tree B, medium for tree A, and highest for tree C (Fig. 3e,f, see also Supplementary Files S3, available on Dryad at http://dx.doi.org/10.5061/dryad.0sn23). In general, the chance of inferring the entire topology correctly is much higher for trees of 8 than for those of 16 taxa, as one might expect for a limited number of characters. Compared with a stationary model, topology inference improves noticeably under the directional model when the data have been simulated with a minimum level of directionality. Rooting of the tree under the directional model is easiest for trees A and C and hardest for the unbalanced tree B with ultrametric branch lengths (Fig. 2b,d,f). In all cases, stronger directionality (i.e., larger differences between the root and stationary frequencies) improves the correct placement of the root.

The success of inferring the root frequencies and stationary frequencies of the directional model varies with the evolutionary rate measured from root to tip (Fig. 3c,d). At low rates, the process is still far from reaching its equilibrium state, and the estimated stationary state frequencies are strongly biased toward the starting frequencies. At rates above 0.5, saturation diminishes the signal about the true state frequencies at the root. The most accurate inference of both root and stationary state frequencies occurs at intermediate rates (0.05 to 0.2), but even in those cases, the difference between root and stationary frequencies is often underestimated.

When the trees are not fixed to the true topology and very high rates are assumed (two or more expected substitutions between root and tips), we observe false positives in the recovery of the directional model, and the same is the case with entirely random data (Supplementary File S4, available on Dryad at http://dx.doi.org/10.5061/dryad.0sn23). Further analyses show that under such circumstances, that is without reliable information about topology and relative branch lengths, and in combination with an inadequate branch length prior (we used an exponential distribution with rate 10 throughout), the analyses recover a balanced tree with extreme directionality because this parameter combination increases the probability of observing a large number of substitutions on short branches.

Nonclock Analyses of the Hymenopteran Data Set

The analysis of the morphological data alone results in highest posterior probabilities for the directional model in the case of wing veins and muscles, and just below significant support for the sclerites (Table 1 and Fig. 4a–c). The root and stationary state frequencies differ most strongly in the wing veins (Fig. 4a), followed by the muscles (Fig. 4b), while they are rather similar in the case of the sclerites (Fig. 4c). The Bayes factor estimated via stepping-stone sampling also strongly favors the directional model (71.6 on the log scale). While MCMC sampling produces very similar topology samples across runs under all models, convergence on tree height and partition-specific rate multipliers is slower under the directional model (PSRF values around 1.04), which indicates a problem with estimating relative branch lengths under this model when the amount of data is limited. The estimated tree height and rate multipliers are also quite sensitive to the branch-length prior. The tree height estimate increases from a median of 6.59 under an exponential branch-length prior with mean 0.1 to 117 under a prior with mean 1. Convergence, especially of the partition-specific rate multipliers, is much worse under the more diffuse branch-length prior, with PSRF values of up to 1.86.

Adding molecular data resolves the problem of inferring branch lengths and rate multipliers. The posterior probability of the directional model remains maximal for wing veins and muscles, but drops to 0.8 for the sclerite partition (Table 1), with partly overlapping posterior densities recovered for the root and stationary state frequencies (Fig. 4f).

Adding Fossil and Time Information

Fossils and the corresponding time information promise a much better chance for inferring the actual evolutionary process than when relying on data from extant taxa alone. The fossil data added here are very unevenly distributed among the character partitions, with 99% of the wing veins coded for the fossils, 7.2% of the sclerites, but none of the muscles. Adding the fossil information in a total-evidence dating approach does not change the main conclusions, with wing veins and muscles still achieving highest support for the directional model (Table 1). The difference between estimated root and stationary frequencies increases even further and estimates are more precise than those from the analyses with extant taxa alone (Fig. 4g,h). The posterior densities of the state frequencies of the sclerite partition, however, now overlap broadly (Fig. 4i), indicating no support for the directional model for this partition. The extent of directionality in certain character systems becomes apparent when plotting the wing veins of selected fossil and extant taxa on a random tree drawn from the posterior distribution (the topology of which is concordant with the majority-rule consensus tree from the RJ analysis, Fig. 5).

The stationary and directional models, as well as the RJ approach sampling under both models, all retrieve similar median ages for Hymenoptera of approximately 305 Ma (Table 2). However, the precision of the estimates increases from the stationary to the directional model, and even further in the RJ approach, with the 95% credibility intervals decreasing from 63 to 38 million years (Table 2).

A Nonstationary Markov Model

While the standard model for analyzing discrete morphological characters assumes stationarity of the evolutionary process (Lewis 2001), the nonstationary Markov model for binary characters suggested here allows the modeling of biological processes that have not yet reached their equilibrium state. The only relaxation of the Lewis model is that the state frequencies at the root of the phylogenetic tree are allowed to differ from the equilibrium frequencies of the homogeneous process running along the rest of the tree. The rate matrix is potentially asymmetric, allowing different forward and backward rates between the two states. The model can also accommodate heterogeneous distributions of state frequencies among terminals, if these are caused by differences in evolutionary rates across lineages.

In terms of macroevolutionary inference, the model proposed here might be useful to distinguish between stasis and directional change. While the state frequencies at the root coincide with the stationary frequencies in the former, the difference between them drives directionality in the latter. This model is especially useful in contexts where the evolutionary process within a group is believed to differ from the conditions in its sister groups and thus at the root of the tree, for instance after the appearance of a key innovation, the colonization of a new geographical area, or the occupation of a new ecological niche. While we apply the model to morphological data coded with nonarbitrary (“presence–absence”) state labels (Lewis 2001), it is also suited for other binary data types, for example presence or absence of introns, restriction sites, gene duplications, or lateral gene transfers, or even behavioral or ecological traits (e.g., Nylander et al. 2004; Alekseyenko et al. 2008; Klopfstein et al. 2010; Stern et al. 2010; Klopfstein and Ronquist 2013). While currently only implemented for binary data, an extension of the model to more states is possible, for instance to model the evolution of nucleotide or amino acid sequences.

Other variations of the finite-state Markov model that also relax the assumptions of stationarity, time-reversibility, and homogeneity have been proposed (for an overview, see Jermiin et al. 2008; Jayaswal et al. 2014). Barry and Hartigan (1987) detailed what is now referred to as the GMM, which assumes separate, arbitrary substitution matrices along each branch of the phylogeny and is thus both very flexible and very parameter rich. It has been demonstrated to suffer from identifiability problems (Zou et al. 2011). For amino acids, Blanquart and Lartillot (2008) suggested a model that is both site and time heterogeneous; even though heterogeneous, it can be seen as stationary in the broad sense because they assume that the state frequencies at the root of the tree correspond to the stationary frequencies of the mixture of their heterogeneous processes. The same applies to a three-state Markov model employed by Stern et al. (2010) to study lateral gene transfer, and to two more general models proposed by Jayaswal et al. (2011b).

Truly nonstationary models were suggested repeatedly (Yang and Roberts 1995; Galtier and Gouy 1998; Pagel and Meade 2004; Jayaswal et al. 2005, 2007; Blanquart and Lartillot 2006; Zou et al. 2012) to address phylogenetic inference from nucleotide or amino acid sequences that vary in base composition over time. Most of them introduced restricted versions of the GMM which allow for some heterogeneity of the process, but use far less parameters. In their models, only the state frequencies are allowed to vary among branches, while a single, simplified vector of exchangeability rates is assumed for the entire tree. To circumvent the problem of overparameterization that comes with branch-wise models, alternatives were suggested that either use a fixed number of different compositions (Foster 2004), use clustering methods and/or likelihood-based model testing to determine optimal numbers of parameters (Jayaswal et al. 2011a, 2014; Dutheil et al. 2012; Groussin et al. 2013), or employ a Poisson process to model changes in the stationary state frequencies on the tree (Blanquart and Lartillot 2006).

All these models relax the assumption of stationarity while at the same time allowing for heterogeneity of the process along the tree. However, assuming time-heterogeneity beyond the origin of the studied group might not always be necessary to achieve an adequate fit of the model to the data at hand. Our model can be seen as the simplest possible nonstationary model; it will thus be particularly useful for situations where there is a limited amount of data and in contexts where the largest deviations from stationarity are expected at the root of the tree.

Performance in Simulations

Our simulations suggest that the best conditions for detecting directional evolution are unbalanced trees with nonclock branch lengths, a fixed topology, many taxa, and a strong directional drive (the latter being determined by the difference between the state frequencies at the root and at equilibrium). But even under such favorable conditions, the signal left by the directional process can only be detected for a range of evolutionary rates that typically span 0.1 to 0.5 expected substitutions between the root and the tips of the tree. Interestingly, similar rates have been shown to also be optimal for phylogenetic inference and branch length estimation (Goldman 1998; Felsenstein 2004; Susko and Roger 2012).

The shape of the tree also has a profound effect on the detection of directional evolution, with unbalanced trees being beneficial, especially if they have short terminal branches attaching close to the root. Such trees might result in contexts where a key innovation speeds up speciation and evolutionary rates at the same time (Mooers and Heard 1997). The shape of the hymenopteran tree, with the very slowly evolving family Xyelidae as the sister group to the remaining taxa (Ronquist et al. 2012a), is certainly reminiscent of such a pattern.

Even under favorable conditions and at suitable evolutionary rates, the differences between stationary and root frequencies are often underestimated in our simulations, as stationary frequencies were not yet reached in the terminals before the signal of the root frequencies began to erode. This relates to the difficulties associated with inferring evolutionary processes in the distant past from extant taxa only. Even when using adequate models that allow for directionality, inference based on extant taxa only can be misleading, both positively and negatively, while the incorporation of fossil information leads to a dramatic increase in accuracy and precision both in simulated and empirical data (Slater et al. 2012; Finarelli and Goswami 2013). Nevertheless, especially under certain tree shapes, the signal left by directional evolution might still be detectable in extant taxa using parametric statistical inference under appropriate models.

Allowing the state frequencies at the root to differ from the stationary state frequencies requires the computation of rooted instead of unrooted trees. Our simulations show that adequately modeled, directional evolution can even be used to infer the root of a tree, at least for a certain period after the onset of the process (Fig. 3e,f, Huelsenbeck et al. 2002). Moreover, not accounting for directionality cannot only mislead estimation of the ancestral states in a phylogeny (Susko and Roger 2013), but also of the phylogeny itself, as is evident from our simulation results and has been shown previously for model misspecifications in general (e.g., Felsenstein 1978; Huelsenbeck and Rannala 2004).

Our simulations were intended to mimic morphological data sets and were thus performed assuming an ascertainment bias which excludes constant characters. Especially at low evolutionary rates, which result in numerous constant characters, this will probably have an impact on the inference of directionality, as constant characters can retain a lot of information about the conditions at the root. For studies based on different data types, for instance nucleotide or amino acid sequences which contain constant characters as well, recovery rates of the directional model might thus be higher. Further studies on simulated and empirical data sets are needed to evaluate the interplay between ascertainment biases and the inference of directional evolution.

At high evolutionary rates and also with random data, we found erroneous preference for the directional model under inappropriate branch-length priors. Support for the directional model was often very high under such conditions, but largely disappeared when additional information was provided about the topology and relative branch lengths, or if a more adequate branch-length prior was used. To improve confidence in the inference of directionality, it is thus advisable to examine the sensitivity to and adequacy of the branch-length prior, and to add complementary information about the topology and relative branch lengths, for instance by adding molecular data. Including information from fossil taxa, if available, can also help address the problem of false positives.

The Evolution of Hymenopteran Morphology

Our results strongly support directional evolution in hymenopteran wing venation and musculature, while the picture is more complex in the case of sclerites. Even though the posterior probability of directionality in the evolution of sclerites was just shy of the 95% mark in the pure-morphology analyses, adding molecular or fossil data weakened the support considerably. Therefore, it seems likely that the signal in the pure-morphology analyses is spurious. Why is there no directionality in the evolution of sclerites? One possibility is that the hymenopteran exoskeleton evolves under conflicting selection pressures. On the one hand, there may be advantages to the structural simplifications that can result from the loss of sclerites. On the other hand, the exoskeleton might compensate for the loss of muscles by increasing in complexity, using more parts and more complex interactions between parts to allow better control of movement by fewer muscles. Several examples could be mentioned to support this idea. However, until we achieve a more profound understanding of the functional morphology of the hymenopteran body, this largely remains the realm of speculation.

With respect to wing venation and musculature, our results suggest that hymenopteran evolution has essentially been a story of losses. In the case of wing venation, the systematic trend toward simplification across multiple lineages is beautifully illustrated by the fossil record (Fig. 5). In fact, the frequencies of presence and absence of wing veins among extant terminals included in our analysis are quite different from the estimated stationary frequencies, suggesting that this character system is still far from reaching equilibrium. The same holds true for muscle characters. It is interesting to note that we were able to document such strong directional trends despite the fact that our analysis included so few representatives of the Apocrita, the explosive terminal radiation of originally insect-parasitic hymenopterans. It is especially evident within this group that reduction of morphological features is correlated with reduction in body size, with the most simplified morphologies found among the smallest hymenopteran taxa, for example Ceraphronoidea (Miko et al. 2013), Chalcidoidea (Krogmann and Vilhelmsen 2006), Mymarommatidae (Vilhelmsen and Krogmann 2006), and Platygastridae (Austin et al. 2005). Whether the trend toward simplification in hymenopteran morphology is due to directional selection, for instance through lower costs of producing fewer wing veins or muscles, or due to developmental constraints favoring losses over gains, is still unclear. However, the fact that the trends are so strong even among basal hymenopteran lineages, mostly including insects of moderate to large body size, suggests that reduction in body size is at least not the only causal factor. Dollo-like evolution (Gould 1970; Dollo 1893), that is, the irreversibility of losses of complex traits, probably also plays an important role. As evolution toward smaller body sizes has happened several times independently and was probably reversed multiple times during hymenopteran evolution, and as the loss of many muscles and wing veins seems to be a consequence of small body size (Vilhelmsen et al. 2010), one could even envision an interaction of these two causes, combining in a ratchet-like mechanism leading to an accumulation of losses in different body-size-dependent character systems (c.f. Lukes et al. 2011).

Few other studies have found empirical support for directional evolution. Hunt (2007) compiled a comprehensive set of over 250 fossil sequences monitoring trait evolution in single lineages over mostly 0.5 to 6 million years. He found that only 5% of these sequences showed statistical support for directional change as opposed to stasis or unbiased random walks, even though these fossil sequences had often been recorded with such a change in mind. Hunt's study differs strongly from our analysis, both in temporal or taxonomical scope and in methodology. Using data from several lineages in a phylogenetic framework, as we exemplify here, certainly increases the power to detect directional patterns dramatically. Furthermore, instead of examining single characters, our data set included 23–67 characters in each of the morphological partitions, and our simulations even assumed 100 characters, which certainly increased the statistical power considerably. However, a researcher might in many cases only have a few characters or even a single character for which she would like to test directionality; it is reasonable to assume that our model is useful also in such a context. The reduced power can be compensated for by increasing the number of terminals sampled; additional simulations are needed to evaluate the required number. In the context of character-dependent diversification (the BiSSE model, Maddison et al. 2007), the authors of a recent simulation study suggest using at least 300 terminals (Davis et al. 2013), but our model is less complex than theirs and might thus be sufficiently powerful even with a smaller sample. Sanderson (1993) studied rate asymmetry and found that given a sufficient evolutionary rate, as few as 32 taxa can be enough to distinguish between evolutionary models.

Even though our model was clearly better than the stationary model in explaining the evolution of hymenopteran morphology, it might be possible to improve the fit to this particular data set even further. In particular, it would be interesting to explore the possibility of allowing a limited number of changes between evolutionary regimes across the tree, instead of assuming a single regime. For instance, the poorly sampled Apocrita represent a group that has undergone rather dramatic changes at its base, with the shift to a parasitic lifestyle, the appearance of a wasp waist, and strongly increased rates of diversification (Heraty et al. 2011; Aguiar et al. 2013; Klopfstein et al. 2013). One might suspect that these changes are associated with detectable modifications in the overall process of morphological evolution. A first step toward modeling process heterogeneity would thus be to allow for a separate model in this part of the tree. However, similar shifts might of course have occurred in other parts of the tree as well. Ideally, one would like to infer both the number of different evolutionary regimes and the points of change between them. In the extreme case, such a heterogeneous model of morphological evolution might estimate different regimes on each branch of the tree (cf. Yang and Roberts 1995; Galtier and Gouy 1998; Jayaswal et al. 2005), which might lead to overparameterization. Some model component that puts a fairly strong constraint on the number of regimes would be advisable for typically rather scarce morphological characters, such as a compound Poisson process or a Dirichlet process prior (cf. Blanquart and Lartillot 2006). Unfortunately, models like these are often difficult to sample from efficiently using standard Bayesian MCMC algorithms, and the statistical power is likely to decrease rapidly with the number of postulated evolutionary regimes. How to model the process that changes evolutionary regimes is a critical issue in this context.

Another interesting possibility would be to consider dependence among discrete morphological characters. We treated the characters as independent of each other, which could potentially be problematic if their evolution was in fact tightly linked, even though the effect of such linkage among characters on the inference of directionality is not entirely clear. If character change is fundamentally directional, and the linkage is a simple matter of the state of some characters affecting the rate of change of others, then the net effect might simply be overconfidence in the directional nature of evolution. One could also imagine a scenario in which atomic changes are essentially nondirectional, but the dependence among characters creates the impression of directionality. However, to be able to distinguish such character-dependence scenarios from independent directional and stationary models, we will probably need detailed models of the nature of the potential linkage among morphological characters. In formulating such models, both developmental and functional links need to be considered.

In hymenopteran wing veins, we generally observe a reduction from the posterior to the anterior margin of the wing and from its apex toward the base, which could form the basis of character-dependence models. In muscles and sclerites, there are a few cases of potential developmental or functional constraints, for instance in serially homologous structures such as the meso- and metafurca, meso- and metatibial spurs, and certain sets of muscles in the meso- and metathorax, or in functionally linked muscle groups. Even though it would be interesting to try and model such character-linkage phenomena, we probably need larger data sets (more terminals and more characters) to convincingly demonstrate their presence. For instance, it seems likely that the underlying developmental constraints in wing venation (e.g., the presence of more basal veins being necessary for the formation of a more apical vein) will only start playing a major role once a certain number of wing veins have been lost, which is not the case in the relative modest range of venational complexity seen in the taxa sampled for this study. Similarly, the observed distributions of skeletal and musculature features among the taxa examined here suggest that it would be difficult to demonstrate any potential character linkages in these systems (Vilhelmsen et al. 2010). Nevertheless, modeling character dependence certainly has a lot of potential, as exemplified for molecular evolution (e.g., accounting for secondary or tertiary structure in rRNA and proteins; Schoeniger and von Haeseler 1994; Robinson et al. 2003). While not yet feasible for most morphological character systems, the formulation of more complex evolutionary models will in the future greatly increase the rigor of phylogenetic analyses involving character sets that violate the standard models; further research in this area is sorely needed. Another possible extension to our approach would be to model the correlation between the evolution of morphological characters and the evolutionary process itself. Using this approach, one might for instance be able to account for the influence of body size on the extent and speed of evolutionary reductions in morphological character systems. Such a framework has been presented for quantitative characters based on multivariate Brownian motion models (Lartillot and Poujol 2011). This type of model remains to be extended to discrete characters, but the Brownian framework would actually be adequate for the study of the correlation between body size and evolutionary transition rates, both of which are quantitative traits. For Hymenoptera, the problem here lies more with the data than the analysis framework: It appears unlikely that a sample of less than 1% of the described taxa would be sufficient to infer body size evolution with any degree of precision.

These are just a few examples of possible extensions to our model. Clearly, there is still much to learn from more sophisticated probabilistic modeling and statistical inference of morphological evolution.