The Hymenopteran Tree of Life: Evidence from Protein-Coding Genes and Objectively Aligned Ribosomal Data

Previous molecular analyses of higher hymenopteran relationships have largely been based on subjectively aligned ribosomal sequences (18S and 28S). Here, we reanalyze the 18S and 28S data (unaligned about 4.4 kb) using an objective and a semi-objective alignment approach, based on MAFFT and BAli-Phy, respectively. Furthermore, we present the first analyses of a substantial protein-coding data set (4.6 kb from one mitochondrial and four nuclear genes). Our results indicate that previous studies may have suffered from inflated support values due to subjective alignment of the ribosomal sequences, but apparently not from significant biases. The protein data provide independent confirmation of several earlier results, including the monophyly of non-xyelid hymenopterans, Pamphilioidea + Unicalcarida, Unicalcarida, Vespina, Apocrita, Proctotrupomorpha and core Proctotrupomorpha. The protein data confirm that Aculeata are nested within a paraphyletic Evaniomorpha, but cast doubt on the monophyly of Evanioidea. Combining the available morphological, ribosomal and protein-coding data, we examine the total-evidence signal as well as congruence and conflict among the three data sources. Despite an emerging consensus on many higher-level hymenopteran relationships, several problems remain unresolved or contentious, including rooting of the hymenopteran tree, relationships of the woodwasps, placement of Stephanoidea and Ceraphronoidea, and the sister group of Aculeata.


Introduction
The Hymenoptera (sawflies, wasps, bees and ants) are one of the four largest insect orders, with more than 146,000 described species [1] (J.T. Huber, personal communication). The oldest fossils belong to the family Xyelidae and date back to the middle Triassic (about 235 Ma) [2], but recent age estimates based on molecular data suggest a much earlier origin in the late Carboniferous (about 311 Ma) [3,4]. Hymenoptera assume a wide range of different life styles, from phytophagous to parasitic and predatory [1,5], occupy a wide range of ecological niches, and have undergone several transitions to eusociality [6,7]. Most species live as parasitoids of other insect larvae and thus fulfill a vital role in most terrestrial ecosystems, and many aculeates are economically important pollinators or predators. Despite their ecological and economic importance, especially the parasitic hymenopterans are one of the most severely understudied insect groups, with large regions of the world virtually unexplored, and undescribed species discovered at a regular pace even in well-studied faunas in the Western Palearctic and Nearctic [8]. Conservative estimates suggest that over 600,000 species of Hymenoptera may exist [9], although much higher numbers of 1-2.5 million species have been proposed [10,11].
The history of hymenopteran phylogenetic research dates back to pre-cladistic times, when the traditional division into the Symphyta (sawflies and woodwasps, without a wasp waist) and Apocrita (hymenopterans with a wasp waist) was established, as was the paraphily of the former with respect to the latter [12]. Apocritans were further divided into the Parasitica (parasitoid wasps) and Aculeata (stinging wasps), with the latter believed to be nested within the former. Rasnitsyn, in a series of seminal papers examining the morphology of both recent and fossil taxa [2,13] (and references there-in), proposed a very influential phylogenetic hypothesis. One of the most innovative aspects of this hypothesis was the division of the Apocrita into four clades, the Evaniomorpha, Proctotrupomorpha, Ichneumonoidea ('Ichneumonomorpha') and Aculeata ('Vespomorpha'), only the last two of which had been recognized previously. Rasnitsyn was also the first to provide convincing evidence for the monophyly of 'Vespina', consisting of the sawfly family Orussidae and the Apocrita. However, Rasnitsyn never provided an explicit quantitative analysis, and a later attempt to specify the character observations and subject them to cladistics analysis [14] indicated that there was little objective support for the proposed groupings in the Apocrita.
Since then, several morphological and early molecular studies have improved our understanding of hymenopteran relationships while leaving many questions open. Sharkey [12] summarized earlier attempts to reconstruct the hymenopteran tree of life, setting the stage for a concerted effort of many international specialists collaborating under the Hymenoptera Tree of Life project (HymAToL). Three papers on higher-level hymenopteran relationships stemming from this project have recently been published, relying on morphology [15], molecular data [16], and both [17]. Vilhelmsen et al. [15] described 273 morphological characters from mesosomal anatomy, scored them for 89 species across the hymenopteran tree, and assessed their phylogenetic information content. Heraty et al. [16] analyzed approx. 6.2 kb of molecular sequences from four markers: the ribosomal 18S and 28S, the mitochondrial cytochrome oxidase 1 (CO1) and one copy of the nuclear elongation factor 1-α. They used both parsimony and statistical approaches. Sharkey et al. [17] combined the molecular dataset, Vilhelmsen et al's [15] mesosomal characters, and 115 additional morphological characters from other parts of the body into a total-evidence dataset which they analyzed under the parsimony criterion.
Briefly, these studies show that morphological data resolve part of the basal sawfly grade but contain little information about relationships above the superfamily level in the Apocrita. The molecular data, in contrast, shed considerable light on apocritan relationships. For instance, they support the monophyly of Proctotrupomorpha, while showing that the Aculeata are nested within a paraphyletic Evaniomorpha. They also corroborate the monophyly of the much discussed Evanioidea (including Gasteruptiidae, Aulacidae and Evaniidae), and identify several novel groupings, such as the 'core Proctotrupomorpha' (Proctotrupomorpha without Cynipoidea and Platygastroidea), the Diaprioidea (Diapriidae, Monomachidae and Maamingidae), the 'core Proctotrupoidea' (Proctotrupoidea without Diaprioidea), and a clade consisting of Trigonaloidea + Megalyroidea. At the same time, the molecular data leave many parts of the apocritan tree unresolved, in particular relationships within Aculeata and Evaniomorpha. More disturbingly, they also suggest groupings that conflict strongly with morphology-based conclusions on sawfly relationships. In particular, they fail to support the established consensus view on woodwasp relationships and, depending on alignment, even fail to support the monophyly of Apocrita itself, placing the Orussoidea among Evaniomorpha lineages.
One of the major problems in interpreting the molecular results is that they are based to a large extent on ribosomal data. The ribosomal sequences (18S and 28S) comprise almost three quarters of the HymAToL data, and an even larger fraction of the phylogenetically informative sites. Ribosomal sequences are challenging to align correctly, especially at the evolutionary distances involved in higher hymenopteran phylogeny, and all currently available methods involve some compromises. Heraty et al. [16] employed two approaches, a by-eye alignment and an alignment based on predicted secondary structure; Sharkey et al. [17] chose to use the former. Both methods rely on human judgment and hence the results may have been influenced by preconceived notions of phylogenetic relationships. As evidenced by the differences between the results based on the by-eye and secondarystructure alignments [16], the alignment method can strongly affect phylogenetic inference.
One way to remove potential alignment bias from the equation is to align the ribosomal sequences using methods that do not involve subjective human input. Another possibility is to infer relationships based entirely on easily aligned proteincoding sequences, but until now there have not been enough protein-coding data available. In this paper, we explore both tactics. First, we explore objective alignment of the ribosomal data. Ideally, the alignment should be based on models including nucleotide substitutions as well as insertion and deletion events, and phylogenetic inference should accommodate alignment uncertainty. In principle, such methods are available in a Bayesian framework [18,19], but they are still too computationally expensive to be applied to the HymAToL data. Instead, we use a two-step approach in which we obtain a ribosomal alignment without or with very little subjective human input first and then subject it to analysis using standard methods. Specifically, we use two methods for obtaining the ribosomal alignments: i) a fully objective, iterative approach using MAFFT [20]; and ii) a semi-objective Bayesian approach based on an explicit model of indel evolution, as implemented in the program BAli-Phy [21], applied to subalignments that are then pieced together. Second, we add three nuclear protein-coding genes to the HymAToL dataset: RNA polymerase II, the carbamoyl phosphate synthase domain of CAD, and the F1 copy of elongation factor 1-α. This allows us for the first time to infer higher relationships across the Hymenoptera based entirely on protein-coding data (4.6 kb from five markers). Finally, in order to identify the origins of different, sometimes conflicting, phylogenetic signals in the resulting data, we conducted in-depth analyses of the different data partitions separately and combined in a fully stochastic, Bayesian framework.

Taxon sampling and molecular methods
Our taxon sampling is largely based on the HymAToL sampling as described in Heraty et al. [16] and Sharkey et al. [17], with minor modifications. While excluding some of the aculeate taxa with low gene coverage that were overrepresented in the data matrix, we added a representative of an additional family, the Megalodontesidae (Pamphilioidea). In total, we included 110 hymenopteran species covering 66 families and all 22 superfamilies [12], and 27 outgroup taxa ( Table 1).
The previous data matrix from the HymAToL project encompassed, unaligned, about 1,400 bp of 18S rRNA (by-eye alignment: 2,014 bp, secondary structure alignment: 1,860 bp), about 3,000 bp of 28S rRNA (by-eye alignment: 4,681 bp, secondary structure alignment: 3,252 bp; both after exclusion of unreliably aligned portions), 770 bp of CO1 mtDNA, and 1,040 bp of the coding region of the F2 copy of elongation factor 1-α (EF1α-F2). To these four markers, we added sequences from three nuclear, protein-coding genes: 990 bp of the carbamoylphosphate synthetase domain of the Conserved ATPase Domain (CAD), 800 bp of RNA polymerase II (POL), and 1,040 bp of the F1 copy of the elongation factor 1-α (EF1α-F1). The F1 and F2 copies of EF1α in Hymenoptera originate from a duplication event that took place before the radiation of the order and the two copies evolved independently since [22].
Laboratory protocols followed Heraty et al. [16] and Klopfstein and Ronquist [22]. The wide taxonomic scope of this study necessitated the use of a range of primer pairs for different taxonomic groups. While primers and protocols used for the EF1α-F1 sequences are given elsewhere [22], primers for CAD and POL are listed in Supplementary Table S1. Gene coverage was 84%, so on average six of the seven genes were sequenced per taxon. The 18S and 28S genes were sequenced for all taxa, CO1 for 92%, EF1α-F2 for 85%, EF1α-F1 for 61% (75% in Hymenoptera), POL for 76% and CAD for 80% of the taxa. Genbank accession numbers are given in Table 1.

Multiple-sequence alignment
Protein-coding genes were aligned in Mega5 [23] after translation into amino acids. Few gaps were detected, and alignment was straightforward. Introns were identified by alignment against known coding regions from Genbank (Table  1) and their exact position conditioned on the presence of GT-AG splicing sites. Introns were not objectively alignable and were removed from all further analyses.
For the MAFFT alignment of the ribosomal sequences, we used the E-INS-i algorithm as available on the web server at http://mafft.cbrc.jp/alignment/server/ with all parameters at their default values [20]. This algorithm has been shown to be more accurate for difficult alignments than other iterative alignment procedures on a wide range of benchmarks, in several simulation studies [24,25], and also was the preferred alignment algorithm for ribosomal stem regions in analyses of Chalcidoidea [26].
As an alternative approach, we used the program BAli-Phy [19,21] with a model of indel evolution that takes branch lengths into account [27] to obtain MAP (maximum posterior probability) subalignments of subsets of taxa that were later pieced together into a complete alignment. We split our data in four different taxon sets. The first set included all outgroup taxa and one representative of each hymenopteran superfamily, the second set contained all remaining symphytan taxa, the third the species of Proctotrupomorpha, and the fourth the rest of the apocritan taxa. Each of these four taxon sets was then aligned separately in BAli-Phy under a GTR + Γ + I substitution model. In order to speed up convergence, we introduced multiple alignment constraints. To do so, we examined the secondary-structure alignment from Heraty et al. [16] for length-constant stem regions of at least length 10 bp, and fixed the alignment at a conserved base in the middle of each such stem. A total of 48 and 85 alignment constraints were invoked for 18S and 28S, respectively. Because the 28S alignment used too much memory to be run in a single analysis, we cut the alignment into two parts at one of the constraint points around the middle of the sequence, and ran it in two separate analyses. For all four taxon sets, which included from 16 to 33 taxa each, we ran four independent runs for seven days (the maximum period) at the National Supercomputer Center in Linköping, Sweden (NSC). Most of the runs did not reach the aspired topology convergence (the average standard deviation of split frequencies (ASDSF) between runs for the different taxon sets was 0.003-0.09 for 18S and 0.04-0.17 for 28S), but the sample of other parameters had reached convergence as judged from effective sample sizes > 100. The MAP alignments obtained from these runs were combined using OPAL [28]. First, we merged the outgroup and backbone taxa with Symphyta, then added the remaining Apocrita without Proctotrupomorpha, and finally merged all of these with Proctotrupomorpha. Alignment and polishing methods were set to "exact", the distance type to "normalized alignment costs", and the polishing approach to "random three-cut". Nineteen of the 18S and 36 of the 28S sequences had missing parts, which were not sequenced. Because BAli-Phy relies on an explicit indel model of evolution and gaps thus become informative, these sequences had to be removed from the BAli-Phy analyses. We added these fragmentary sequences to the final BAli-Phy alignment using the "add" option in MAFFT [29].

Data properties
The variation present in the different genes and gene partitions was examined using the "cstatus" command, and a basic test of non-stationarity of nucleotide composition was performed with the command "basefreqs" in PAUP* [30]. Saturation plots for each gene and for the third codon positions of protein-coding genes were produced by retrieving pairwise uncorrected p-distances in Mega 5 [23], and plotting them against inferred branch-length distances on the tree with the highest likelihood found during the Bayesian tree search based on the single genes (R script available from the first author on request). The third codon positions of all genes showed clear signs of saturation and non-stationarity (Table 2 and Figure 1), so we also analyzed our data after excluding them.
In order to get a rough estimate of the performance of the different genes (or of their contribution to the final phylogenetic inference) and to assess the quality of the two alignment approaches for the rRNA partition, we compared the Bayesian tree samples obtained from the single-gene analyses and from an analysis of morphology alone (see below) to the proteincoding and total-evidence tree samples. As a measure of topological distance, we used ASDSF values as obtained with the 'sumt' command in MrBayes 3.2 [31]. We compared 10,000 trees from each set after reducing the trees to taxa shared in all Parnips nigripes

Phylogenetic analyses
We performed a number of different Bayesian analyses on parts of the dataset in order to discern the sources of different signals and conflict (Table 3). These analyses include two different alignment options for the ribosomal RNA genes (18S and 28S), molecular-only and total-evidence analyses which included the morphological partition, analyses of the ribosomal and protein-coding genes separately, in the latter case including or excluding third codon positions (third codon positions of CO1 were always excluded), and finally single-gene analyses. The protein-coding genes were also analyzed after translation into amino acids and applying a reversiblejump algorithm to integrate over the fixed-rate amino-acid models implemented in MrBayes. The data matrices and associated consensus trees of all analyses are deposited on TreeBase (URL for reviewers: http://purl.org/phylo/treebase/ phylows/study/TB2:S13902?x-access-code=44421680b40bc7867da8bbe7cece2e9c&format=html).
All analyses were performed in MrBayes 3.2 [31]. Where applicable, data were partitioned into genes and into first and second versus third codon positions, with substitution models unlinked across partitions. We used model jumping to integrate over the GTR model subspace ("nst=mixed" option in MrBayes) and modeled among-site rate variation with a four-category gamma distribution and a proportion of invariable sites. The morphology partition was modeled using the standard discrete model [34], a "variable" ascertainment bias, and a fourcategory gamma distribution to model among-character rate variation. For each analysis, we ran four independent runs of four chains each until they had reached topological convergence (ASDSF < 0.05, preferably lower, with 25% of samples discarded as burn-in). In the case of the single-gene analyses of the EF1-α copies, we ran 100 million generations, but ASDSF values remained at about 0.095. In order to capture the uncertainty that might arise through a lack of convergence of the MCMC in these and also in all the other analyses, we scanned the MrBayes output for bipartition frequencies with a standard deviation larger than 0.1 between runs. The corresponding support values are preceded in each tree figure by a question mark, as they might not have been estimated accurately. Samples of all substitution model parameters were adequate in all runs, as judged from the PSRF values being close to 1.0 and effective sample sizes of (usually much) more than 200. In the single-gene analysis of CAD, the outgroup taxa were recovered within Hymenoptera. In order to obtain meaningful signal from this data partition, we repeated the analysis with all outgroups removed, which strongly improved topology convergence. Although we focus on the Bayesian analyses, we also performed maximum likelihood (ML) analyses for comparison. These analyses were conducted on the combined molecular data and the total-evidence dataset, each under both alignment strategies for the ribosomal partitions. We obtained an estimate for the maximum-likelihood tree from RAxML [35] under a partitioned GTR model for the molecular and the Mk model [34] for the morphological partitions, respectively. To assess support, we performed 1000 bootstrap replicates.

Rogue taxa identification
We used a new algorithm to search for rogue taxa, i.e., taxa that are highly inconsistent in their phylogenetic placement [36] in our set of Bayesian trees. The algorithm aims to optimize the relative improvement in clade support achieved by removing single or groups of taxa [37]. As input, we used 1,000 evenly spaced trees from the post-burnin phase of the MrBayes tree sets. The program was accessed via the webserver at http:// exelixis-lab.org/roguenarok.html under the majority-rule threshold, optimizing overall support, and using maximum dropset sizes of two, three and ten taxa. In all cases, these three dropset sizes led to the same rogue taxa being identified. Rogue taxa associated with a raw improvement (sum of increase in support values) of at least 0.5 (Table 3) were excluded and support values of the consensus tree recalculated. On the tree graphs, we indicate these new values for all nodes except those directly below the rogue taxon, which show the original value. Rogues (or groups of rogues) are indicated by dashed branches.

Alignment and analysis of ribosomal RNA
The MAFFT runs resulted in the shortest alignments, 2027 bp and 5,418 bp for 18S and 28S, respectively. The BAli-Phy alignments are much longer, i.e. 2310 bp and 10,557 bp. The harmonic means of the likelihoods of the Bayesian tree samples retrieved from these alignments (treating gaps as missing data) reflect the alignment lengths, with the longer BAli-Phy alignment reaching a much higher likelihood than the shorter MAFFT alignments (lnL values of -119,599 and -106,394 for the MAFFT and BAli-Phy alignments, respectively). Congruence with the trees retrieved from the protein-coding genes, from the total-evidence analysis that included the BAli-Phy alignment, and even from the totalevidence analysis based on the MAFFT alignment is higher for the BAli-Phy than for the MAFFT alignment ( Table 4).
The consensus tree retrieved from the rRNA data based on the MAFFT alignment is provided in Figure 2, together with support values from the BAli-Phy alignment. Despite the very different alignment approaches and resulting alignment lengths, the consensus trees do not differ much, but the support values for the MAFFT alignment are usually lower. Interestingly, differences between alignment approaches concern some of the relationships which also differed between the by-eye and secondary structure alignment in the Heraty et al. study [16], i.e. the rooting of the hymenopteran tree and the placement of Orussoidea. Independent of alignment strategy, the rRNA tree is only poorly resolved around the deeper nodes, in contrast to the results from a similar number of base pairs of protein-coding data.  Figure 3 shows the tree retrieved from first and second codon positions of the protein-coding genes, along with support values obtained when including third codon positions of the nuclear genes (but not of CO1). The symphytan grade is well resolved, with maximal support on most of the nodes, and with Orussoidea placed firmly as the sister group of Apocrita. Within Apocrita, the Proctotrupomorpha, Ichneumonoidea and (Evaniomorpha + Aculeata) clades are recovered, although only the former two have high support. The relationships among these three are unresolved. In general, superfamilies are recovered, with the exception of paraphyletic Xyeloidea, Evanioidea, Chrysidoidea, Vespoidea, and Platygastroidea. The Xyeloidea are however monophyletic both when including the third codon positions and when analyzing the data as amino acids. As with the rRNA data, resolution is rather low among the evaniomorph superfamilies and within Aculeata. Mymaromma, the only representative of the enigmatic Mymarommatoidea, has an incomplete coverage in terms of gene sampling (Table 1). It was identified as a rogue taxon, appearing in different places in the Bayesian tree sample. In the consensus tree, it ended up within Ichneumonoidea, but with low support, and sitting on a very long branch.

Phylogeny of Hymenoptera as inferred from proteincoding genes
Most conflicts with the rRNA tree are weakly supported and/or in areas of the tree which are poorly resolved in both analyses, e.g. the relationships within Evaniomorpha, the placement of Ichneumonoidea and Mymarommatoidea, and the monophyly of Diaprioidea. A notable difference is the sister group of Aculeata, which is the Trigonaloidea + Megalyroidea clade according to the rRNA tree and Evaniidae or Stephanoidea according to the analysis of the protein-coding genes, depending on whether third codon positions were excluded or included. The Bayesian tree samples obtained from single data partitions are compared to the total-evidence and protein-coding trees using the average standard deviation of split frequencies as a measure of topological distance. 1 Resolution of the respective consensus tree after reduction to 43 ingroup taxa present in each dataset, given as the percentage of nodes that were resolved. 2 Total-evidence trees

Combined molecular results and total-evidence results
The Bayesian total-evidence tree (molecular and morphological data combined) based on the BAli-Phy alignment is given in Figures 4 and 5, including support values from the total-evidence analysis based on the MAFFT-aligned rRNA sequences, and from analogous analyses of the molecular data partition only. The tree also includes symbols summarizing the results from the rRNA data and the proteincoding genes when analyzed separately. Most of the deeper nodes and well-established groupings like the Holometabola, Apocrita, and Aculeata are well supported. When ignoring the uncertain positions of Stephanoidea and Ceraphronoidea, the three large groups within Apocrita -the Ichneumonoidea, Proctotrupomorpha and, with less support, (Evaniomorpha + Aculeata) -are also corroborated. Although most of the proposed superfamilies are recovered as monophyletic, usually with high support, there are several exceptions. First, the most basal superfamily Xyeloidea is paraphyletic, with Macroxyela more closely related to the remainder of Hymenoptera. Second, the recently proposed Diaprioidea -including Diapriidae, Maamingidae and Monomachidae -are not supported, although the evidence against its monophyly is weak. Finally, relationships within Aculeata are unstable, and neither Chrysidoidea nor Vespoidea are recovered.
Comparing the total-evidence topology, which included morphological data, to the phylogeny obtained from the molecular data alone, there is considerable congruence, but also two areas where the morphological data have the power to change the molecular results ( Figure 6). First, the grade of woodwasps (Siricoidea, Xiphydrioidea and Cephoidea) leading to the Vespina (Apocrita + Orussoidea) is fully reversed in the two analyses, with the sequence Cephoidea -Siricoidea -Xiphydrioidea -Vespina supported by the former, and Xiphydrioidea -Siricoidea -Cephoidea -Vespina by the latter. The molecular signal is fairly strong in the BAli-Phy alignment but weaker in the MAFFT alignment, showing that there is some alignment-dependent signal from the ribosomal sequences. The second example concerns the positions of Stephanoidea and Ceraphronoidea within the Apocrita but involves relationships that are less well supported.
Maximum likelihood estimates based on both the combined molecular and the total-evidence datasets are given in Figures  S1 and S2, with bootstrap support values obtained under both the MAFFT and the BAli-Phy alignment approaches for rRNA. The ML trees are similar to those obtained from the Bayesian analyses, but differ with respect to the placement of the hymenopteran root, which is between Tenthredinoidea and the remaining hymenopterans in the total-evidence and between a monophyletic Xyeloidea + Tenthredinoidea + Pamphilioidea and Unicalcarida in the molecular analysis. Furthermore, the total-evidence analyses did not recover a monophyletic Evaniomorpha, but the conflicting nodes were associated with very low bootstrap support.

Phylogenetic signal in different data partitions
In order to assess the contribution of the different genes and of morphology, we investigate patterns of variation, resolution of the single-gene or single-partition consensus trees, and their The Hymenopteran Tree of Life PLOS ONE | www.plosone.org    congruence with trees derived from other data partitions. Table  2 summarizes some basic properties of the molecular data by gene and by gene partitions. The ribosomal genes and the third codon positions of the two EF1-α copies showed no to moderate GC-biases (up to 61.5% in the third codon positions of EF1-α F1), whereas CO1 had moderate to strong AT bias (62% and 90% for first plus second and third codon positions, respectively), as is the rule for mitochondrial genes in Hymenoptera [38]. All third codon positions of the protein coding genes are heavily saturated, while there appears to be a favorable signal-to-noise ratio in the ribosomal genes and at first and second codon positions of CAD and CO1 (Figure 1). The first and second codon positions of the other three genes (POL, EF1-α F1, EF1-α F2) show comparatively little variation.
Resolution of the single-partition consensus trees varies strongly (employing the majority rule criterion). Table 4 shows the percentage of resolved nodes after reducing each tree to the 43 ingroup taxa common to all datasets. The rRNA data resolved 95% of nodes irrespective of alignment approach, and morphology and the CAD gene each reached 91%. The other single genes lag behind at 67% to 79%. The ranking of partitions is very similar when not only the 43 completely sampled ingroup taxa, but all taxa available per partition are included, with the difference that CAD now outperforms the MAFFT-aligned rRNA data. A similar picture appears when comparing the topological distances between trees obtained from the single-gene analyses to the protein-coding and totalevidence trees ( Table 4). The rRNA data and the CAD gene consistently rank highest, followed by morphology and CO1, while POL and the two EF1-α copies result in more conflicting topologies [22].

Objective and semi-objective alignments of ribosomal DNA sequences
Arguably the best approach to phylogenetic inference based on ribosomal sequences is to analyze unaligned sequences directly. There are several methods that simultaneously estimate alignment and phylogeny: POY in the parsimony framework [39], and ALIFRITZ [40], BAli-Phy [19] and Luntner et al. [18] in a Bayesian setting. These approaches make use of the information present in gaps when reconstructing the phylogeny, which can greatly improve phylogenetic inference [41]. The Bayesian approaches are particularly compelling in that they integrate over alignment uncertainty when reconstructing phylogenetic relationships. Unfortunately, they are still too computationally complex and converge too slowly to be applicable to most empirical datasets.
In addition to the entirely objective alignment approach using MAFFT, we attempted Bayesian analysis of our unaligned data. However, we had to split our dataset both by taxa and sites in order to reach even marginally acceptable convergence on topology and the parameters of the indel model. Merging the MAP sub-alignments and analyzing the composite matrix in a traditional two-step procedure meant we had to forego the possibility of using the information in the indels and integrating over alignment uncertainty. Nonetheless, our partitioned Bayesian method compared favorably to the MAFFT alignment method, as indicated by higher support values and higher congruence with the protein-coding tree.
Although largely objective, our partitioned Bayesian method does involve human decisions on how to decompose the alignment by taxa and by sites. Splitting by taxa has the greatest potential to bias the results. Difficult sequence positions might be aligned in a different manner in the different sub-problems, and the merging of the sub-alignments by OPAL [28] might not be able to resolve such conflict and thus lead to an exaggerated alignment similarity between taxa that were aligned in the same batch. We minimized such problems by basing the decomposition on results from the analysis of the protein-coding genes. We also searched the final results for potential alignment-induced biases. Nodes that might be involved were at the bases of i) all Hymenoptera, ii) Apocrita, iii) Proctotrupomorpha, and iv) all non-proctotrupomorph apocritans. The first three nodes are present in both the trees derived from the BAli-Phy aligned sequences and those resulting from the MAFFT alignment. The fourth group was not recovered in either analysis. In fact, the grouping of Ichneumonoidea with Proctotrupomorpha instead of with Evaniomorpha plus Aculeata, across alignment decomposition lines, was even retrieved with higher support in the BAli-Phy than in the MAFFT analyses ( Figure 2). The apparent absence of decomposition-induced artifacts, the fact that clade support values were almost always higher in the BAli-Phy than in the MAFFT analysis, and the higher congruence of the tree sample obtained from the BAli-Phy alignment with the trees from the protein-coding genes indicate that splitting the alignment problem based on a few explicit and well-grounded assumptions about relationships may be a good general strategy for improving alignment quality.
Several candidate alignment artifacts were identified based on a comparison of the by-eye and secondary-structure alignments of Heraty et al. [16], and by comparison with the results from the protein-coding sequences. These include the monophyly of Xyeloidea and Evanioidea, and the placement of Orussoidea among Evaniomorpha. If they were artifacts of a subjective alignment in the previous analyses of ribosomal data, they should disappear in our analyses of objective alignments. However, all these signals were clearly present in our re-analyses of the ribosomal data ( Figure 2), even though the support values are generally lower, indicating that subjective bias has possibly augmented these signals.

Implications for the hymenopteran tree of life
Our analyses recover a large part of the higher-level phylogeny of Hymenoptera with high support and strong corroboration from independent data sources. Many of these relationships have been uncontroversial at least in the more recent past, e.g. the monophyly of the Unicalcarida (Hymenoptera without Xyeloidea, Tenthredinoidea, and Pamphilioidea) [42], the grouping of Orussoidea with Apocrita, the monophyly of Aculeata and Proctotrupomorpha, and of most of the superfamilies as outlined in Sharkey [12]. More recent suggestions that we could corroborate here with independent protein-coding data include Trigonaloidea + Megalyroidea, core Proctotrupomorpha (Proctotrupomorpha  excluding  Cynipoidea  and  Platygastroidea), core Proctotrupoidea (Proctotrupoidea without Diaprioidea), and finally the placement of Aculeata within a paraphyletic Evaniomorpha. These results appear to be robust and will probably pass the test of time. As they were discussed at some length in a previous study [17], we will not go into further detail here, but only discuss equivocal relationships.
Several parts of the hymenopteran tree remain unresolved and most of these unstable areas include taxa that were also identified as rogue taxa in one or more of the analyses. Rogues can arise due to several reasons, e.g., insufficient gene coverage or particularly long branches. While on average, the twenty taxa identified as rogues did not have a lower number of genes sampled (one missing gene being the average both of the whole dataset and among the rogue taxa), missing data might still be behind the formation of some of the rogues (e.g., Mymarommatoidea, see below). Most of the controversial relationships were also ambiguous in earlier analyses, and might represent difficult phylogenetic histories like rapid radiations (e.g., Aculeata, see below). Figure 7 summarizes the areas of conflict or uncertainty, and we here give a short summary of the evidence for conflicting hypotheses.
The three at the root. It has been recognized early on that Xyeloidea, Tenthredinoidea and Pamphilioidea are the three superfamilies closest to the root of Hymenoptera [13,14,42]. However, the relationships among these superfamilies are not resolved. The three competing hypotheses that result from the current and recent analyses are shown in Figure 7a. Most of the controversy boils down to the uncertain placement of the root of the order Hymenoptera -either within Xyeloidea, between Xyeloidea and the remaining Hymenoptera, or between (Xyeloidea + Tenthredinoidea) and the rest. Undoubtedly, the problem is caused to a large extent by the deep roots of the order that probably date back to the Carboniferous [3], in combination with the long branches connecting it to the outgroups [43,44]. As Hymenoptera are probably the sister group to all other holometabolous insects [45,46], only a much denser taxon sampling of both holometabolan and hemimetabolan outgroups could help improve the reconstruction of ancestral sequences, and hence help resolve these deep relationships. Unfortunately, such an approach is limited by the fact that extant outgroups are placed in isolated crown groups of their own and cannot break down the long branch leading to Hymenoptera.
In addition to the rooting problem, it is somewhat unclear whether Tenthredinoidea or Pamphilioidea are more closely related to the Unicalcarida. The former hypothesis was found in the Sharkey et al. [17] total-evidence analysis, but with low support, and in the CAD single-gene analysis, but again with a posterior probability of only 0.53. In contrast, the combined protein-coding genes show maximum support for Pamphilioidea as the sister to Unicalcarida. Morphological evidence is also somewhat equivocal about the placement of the hymenopteran root. Putative synapomorphies that could The Hymenopteran Tree of Life support a monophyletic Xyeloidea can be found among the mouthparts, e.g. the labral brush, asymmetric mandibles and elongate maxillary palpi [47]. These features are associated with pollen feeding in the adults and are unique within Hymenoptera. In contrast, the long, compound third segment of the antenna which results from the fusion of several flagellomeres might represent a symplesiomorphy, as it is also found in many early fossil hymenopterans and in the tenthredinoid families Blasticotomidae and Argidae [2].
The woodwasp grade.
The remaining symphytan superfamilies in most analyses form a grade towards Vespina (Orussoidea + Apocrita). The sequence in which they branch off is strongly dependent on the dataset and constitutes one of the two strong conflicts between the morphological and molecular data partitions (Figures 6, 7b). Morphological evidence supporting Xiphydrioidea as sister to Vespina is rather strong; the most convincing proposed synapomorphies for this relationship include a number of characters in the dorsal part of the thorax, e.g., the presence of a transscutal articulation, the reduction of the posterodorsal part of the metapleuron (possibly an incipient step in the formation of the wasp waist in Apocrita), and the loss of a number of thoracic muscles [15]. However, none of the single-gene or various combined molecular datasets supported this relationship, and the combined molecular, protein-coding and CAD single-gene analysis are strongly against. Nevertheless, the signal in the morphological partition is strong enough to resolve this conflict in favor of morphology in the total-evidence analyses.
Placement of Stephanoidea and Ceraphronoidea. These two groups are notoriously difficult to place. In the totalevidence analyses, Stephanoidea is placed as the sister-group of all remaining apocritans, a placement that is supported by several morphological, in particular mesosomal, characters [15]. Again, this conflicts with the protein-coding genes, which place stephanids within Evaniomorpha and potentially as the sister clade to Aculeata. The rRNA data do not provide stable resolution around the nodes in question. The situation is complicated by Ceraphronoidea, which assume very differing positions in different analyses, grouping alternatively with Stephanoidea, with Ichneumonoidea, or as sister to Ichneumonoidea plus Proctotrupomorpha. A sister-group relationship between Ceraphronoidea and Megalyroidea, as recovered in Sharkey et al. [17], was never observed here. Morphology does not provide many reliable characters due to the small size of these wasps. Confidence about the placement of Stephanoidea and Ceraphronoidea will depend on additional data, and will help to refine the status of the highly contested Evaniomorpha.
Placement of Mymarommatoidea. The placement of this family is complicated by their small size, associated reduction of many otherwise informative morphological character systems, and risk of homoplasy in other character states associated with size. The gene sampling for this taxon was not complete in our analysis, and they came out as a rogue taxon on a very long branch in the protein-coding tree. The rRNA data recover them as the sister group of Diaprioidea plus Chalcidoidea, but the support for this placement disappears in the combined molecular analysis. Including morphology added support for the common interpretation of mymarommatids as the sister group of Chalcidoidea [14,15,17,48], but this was sensitive to the alignment approach. More molecular data is needed to resolve this conflict, especially because of the limitations inherent in the morphological data for these tiny wasps.
The sister group of Aculeata. Aculeata are firmly placed within a paraphyletic Evaniomorpha (see next paragraph) in all our analyses. A similar placement was recovered in previous analyses [16,17], and contradicts early hypotheses of a sistergroup relationship between Aculeata and Ichnemonoidea [13,14]. Within Evaniomorpha, however, the relationships are highly unstable, and the sister-group of aculeates remains unclear. Although there is some indication that the strongly supported Trigonaloidea + Megalyroidea clade is sister to aculeates, support is weak, alignment-dependent, and contradicted by the analysis of the concatenated protein-coding genes, which favored either Stephanoidea or Evaniidae as the sister group. Given the low resolution both among evaniomorph superfamilies and within Aculeata, a denser taxon sampling within these groups is probably needed to clarify this question.
Evaniomorpha. The concept of Evaniomorpha, as originally proposed by Rasnitsyn [13], included the superfamilies Stephanoidea, Ceraphronoidea, Megalyroidea, Trigonaloidea and Evanioidea, while excluding Aculeata. The morphological and fossil evidence supporting this somewhat heterogeneous assemblage has always been weak, and Rasnitsyn himself recently proposed that the Evaniomorpha be restricted to the Evanioidea [49]. The circumscription of Evaniomorpha remains unclear even after our analyses, especially with respect to Stephanoidea and Ceraphronoidea, but it should definitely be revised to include Aculeata if it is retained as a concept defining a major apocritan lineage.
Non-monophyletic superfamilies. The superfamilies not recovered as monophyletic in the total-evidence analyses are the following: Xyeloidea, Chrysidoidea, Vespoidea and Diaprioidea. While the Xyeloidea are discussed above, the remaining superfamilies deserve further attention. There are several rather convincing morphological synapomorphies for Chrysidoidea, e.g., the subdivision of the second valvifer of the ovipositor into two articulating parts [50], and their nonmonophyly was in fact not strongly supported; rather, the relationships among aculeate families are very poorly resolved in all our analyses, and a much denser taxon and gene sampling is obviously required to address these relationships. The same is true for Vespoidea, although it has been hypothesized previously that they are paraphyletic with respect to Apoidea [7,51].
The superfamily Diaprioidea was suggested by Sharkey [12] to include Diapriidae, Maamingidae and Monomachidae, based on an earlier molecular analysis [52]. While not retrieved in the total-evidence analysis, which instead suggested paraphily with respect to Mymarommatoidea and Chalcidoidea (but with very low support), Diaprioidea are recovered in the CAD singlegene, the protein-coding, and the combined molecular analyses.
Although the Evanioidea were recovered as monophyletic in the total-evidence and combined molecular tree, they were split into Gasteruptiidae + Aulacidae versus Evaniidae in the protein-coding and CAD single-gene analyses. This superfamily may thus deserve more attention, especially given the weak support from morphology, the most striking putative synapomorphy being the attachment of the metasoma high above the hind coxal cavities [e.g. 13,15].

The future of hymenopteran phylogenetics
Although we present here the most comprehensive study of higher-level hymenopteran relationships to date, many questions of great taxonomic and evolutionary interest remain unresolved; the search for more and better data must thus continue. In the light of the large differences in information content in the genes studied here, it becomes clear that data quality can strongly influence the outcome of studies of deeplevel relationships. The performance of CAD [53] was especially outstanding. With less than 1,000 bp, this marker recovered a largely resolved phylogeny of Hymenoptera that was in close agreement with the total-evidence tree. Overall, data partitions that did not show signs of saturation and at the same time included a relatively large number of parsimonyinformative sites consistently achieved higher congruence with trees derived from independent and total-evidence partitions. This is in line with recent theoretical and phylogenomic studies, which found a connection between evolutionary rate, saturation, and phylogenetic utility of different markers [54][55][56][57][58]. Data quality might thus play a very important role when it comes to utility for phylogenetic inference, and could render it unnecessary to accumulate huge quantities of data even (or maybe especially) for difficult phylogenetic problems.
On the other hand, the lack of resolution in vital parts of the hymenopteran tree as inferred here from seven genes might simply demonstrate the limits of few-gene approaches. Estimates of the numbers of genes necessary for reliable phylogenetic inference depend strongly on the phylogenetic context and the inference method, but have been suggested to lie around 20 [43,59,60]. Gene sampling for Hymenoptera phylogenetics has until now relied mostly on very few genes, with two exceptions. A study of 24 expressed sequence tags (ESTs) in 10 disparate hymenopteran taxa [61] recovered deep-level relationships which were almost invariably controversial and in conflict with any previous study, e.g. Chalcidoidea placed outside Proctotrupomorpha and a sistergroup relationship between the latter and Aculeata. These relationships are likely due to the extremely low taxonomic coverage and potentially also to limited phylogenetic signal in the different markers. Another analysis of phylogenomic proportions made use of all sequence data for Hymenoptera present in Genbank [62]. By developing a bioinformatics pipeline that filtered the vast amount of data for genes with compositional stationarity and defined levels of density and taxonomic overlap, they retrieved about 80,000 sites for 1,100 taxa. The main problem with this dataset was the amount of missing data (more than 98%). The resulting tree had very low resolution, recovered many of the included families as para-or polyphyletic, and placed some taxa in obviously erroneous positions. Nevertheless, some of the superfamilies and undisputed higher-level relationships were recovered in this analysis, which demonstrates the potential of such an approach.
The future of hymenopteran phylogenetics lies in datasets that combine the advantages of each of the afore-mentioned studies, i.e., a dense and balanced taxon sampling [63][64][65], sufficiently large amounts of molecular data, a careful assessment of the quality of this data [55,66], and appropriate analysis methodology. Only the combination of these is likely to resolve the remaining uncertainties in the evolutionary history of a group that originated hundreds of million years ago and diversified into hundreds of thousands of species. Table S1. Commented table of primers used in this study. (PDF) Figure S1. Maximum likelihood tree recovered from the analysis of the combined molecular data (rRNA MAFFT aligned). Support values next to the nodes (or after species pairs) are bootstrap supports obtained from 1000 replicates based on both the MAFFT and the BAli-Phy alignments, respectively. Asterisks stand for maximal support. (TIF) Figure S2. Maximum likelihood tree recovered from the analysis of the combined molecular and morphological data (rRNA BAli-Phy aligned). Support values next to the nodes (or after species pairs) are bootstrap supports obtained from 1000 replicates based on both the BAli-Phy and the MAFFT alignments, respectively. Asterisks stand for maximal support. (TIF)