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 over-represented in the data matrix, we added a representative of an additional family, theMegalodontesidae (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α inHymenoptera 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 inSupplementary 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% inHymenoptera), POL for 76% and CAD for 80% of the taxa. Genbank accession numbers are given inTable 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 ofChalcidoidea [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 withSymphyta, then added the remainingApocrita 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 andFigure 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 protein-coding 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 datasets (43 ingroup taxa), using an R script [32] that was based on the Ape package [33].

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 reversible-jump 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 four-category 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 withinHymenoptera. 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 athttp://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 re-calculated. 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 total-evidence 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 inFigure 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 ofOrussoidea. 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.

Phylogeny ofHymenoptera as inferred from protein-coding genes

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 withOrussoidea placed firmly as the sister group ofApocrita. WithinApocrita, 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 paraphyleticXyeloidea,Evanioidea,Chrysidoidea,Vespoidea, andPlatygastroidea. TheXyeloidea 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 enigmaticMymarommatoidea, 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 withinIchneumonoidea, but with low support, and sitting on a very long branch.

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 ofIchneumonoidea andMymarommatoidea, and the monophyly of Diaprioidea. A notable difference is the sister group of Aculeata, which is theTrigonaloidea +Megalyroidea clade according to the rRNA tree andEvaniidae orStephanoidea according to the analysis of the protein-coding genes, depending on whether third codon positions were excluded or included.

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 inFigures 4 and5, 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 protein-coding 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 ofStephanoidea andCeraphronoidea, the three large groups withinApocrita — theIchneumonoidea, 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 superfamilyXyeloidea is paraphyletic, withMacroxyela more closely related to the remainder ofHymenoptera. Second, the recently proposed Diaprioidea — includingDiapriidae, Maamingidae andMonomachidae — are not supported, although the evidence against its monophyly is weak. Finally, relationships within Aculeata are unstable, and neitherChrysidoidea norVespoidea 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 andCephoidea) leading to the Vespina (Apocrita +Orussoidea) is fully reversed in the two analyses, with the sequenceCephoidea –Siricoidea –Xiphydrioidea – Vespina supported by the former, andXiphydrioidea –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 ofStephanoidea andCeraphronoidea within theApocrita but involves relationships that are less well supported.

Maximum likelihood estimates based on both the combined molecular and the total-evidence datasets are given inFigures S1 andS2, 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 betweenTenthredinoidea 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 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 inHymenoptera [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 total-evidence 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) allHymenoptera, 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 ofIchneumonoidea 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 ofXyeloidea andEvanioidea, and the placement ofOrussoidea 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 ofHymenoptera 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 withoutXyeloidea,Tenthredinoidea, and Pamphilioidea) [42], the grouping ofOrussoidea withApocrita, 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 includeTrigonaloidea +Megalyroidea, core Proctotrupomorpha (Proctotrupomorpha excludingCynipoidea andPlatygastroidea), coreProctotrupoidea (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 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 deep-level relationships. The performance of CAD [53] was especially outstanding. With less than 1,000 bp, this marker recovered a largely resolved phylogeny ofHymenoptera 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 parsimony-informative 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–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 forHymenoptera 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 sister-group 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 forHymenoptera 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–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.

• 1.Huber JT (2009) Biodiversity of Hymenoptera. In: Foottit RG, Adler P.Insect Biodiversity: Science and Society. 1st edition ed. Oxford, UK: Wiley-Blackwell. [Google Scholar]
• 2.Rasnitsyn AP (2006) Ontology of evolution and methodology of taxonomy. Palaeontol J40: S679-S737. doi:10.1134/S003103010612001X. [Google Scholar]
• 3.Ronquist F, Klopfstein S, Vilhelmsen L, Schulmeister S, Murray DLet al. (2012) A total-evidence approach to dating with fossils, applied to the early radiation of the Hymenoptera. Syst Biol61: 973–999. doi:10.1093/sysbio/sys058. PubMed:22723471. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 4.Davis RB, Baldauf SL, Mayhew PJ (2010) Many hexapod groups originated earlier and withstood extinction events better than previously realized: inferences from supertrees. Proc R Soc Lond B Biol Sci277: 1597-1606. doi:10.1098/rspb.2009.2299. PubMed:20129983. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 5.Vilhelmsen L, Turrisi GF (2011) Per arborem ad astra: Morphological adaptations to exploiting the woody habitat in the early evolution of Hymenoptera. Arthropod Struct Dev40: 2-20. doi:10.1016/j.asd.2010.10.001. PubMed:20951828. [DOI] [PubMed] [Google Scholar]
• 6.Cardinal S, Danforth BN (2011) The antiquity and evolutionary history of social behavior in bees. PLOS ONE6: e21086. doi:10.1371/journal.pone.0021086. PubMed:21695157. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 7.Pilgrim EM, von Dohlen CD, Pitts JP (2008) Molecular phylogenetics of Vespoidea indicated paraphily of the superfamily and novel relationships of its component families and subfamilies. Zool Scripta37: 539-560. doi:10.1111/j.1463-6409.2008.00340.x. [Google Scholar]
• 8.Bacher S (2012) Still not enough taxonomists: reply to Joppaet al. Trends Ecol Evol27: 65-66. doi:10.1016/j.tree.2011.11.003. PubMed:22138045. [DOI] [PubMed] [Google Scholar]
• 9.Heraty JM (2009) Parasitoid biodiversity and insect pest management. In: Foottit RG, Adler P.Insect biodiversity: science and society. Springer-Verlag Press; pp. 445-462. [Google Scholar]
• 10.Grissell EE (1999) Hymenopteran diversity: Some alien notions. Am Entomol45: 235-244. [Google Scholar]
• 11.Stork NE (1996) Measuring global biodiversity and its decline. In: Reaka-Kudla ML, Wilson DE, Wilson EO.BiodiversityII. Washington, D.C.: Joseph Henry Press; pp. 41-68. [Google Scholar]
• 12.Sharkey MJ (2007). hylogeny Classifications Hymenoptera Zootaxa1668: 521-548 [Google Scholar]
• 13.Rasnitsyn AP (1988) An outline of evolution of the hymenopterous insects (Order Vespidae). Oriental Insects22: 115-145. [Google Scholar]
• 14.Ronquist F, Rasnitsyn AP, Roy A, Eriksson K, Lindgren M (1999) Phylogeny of the Hymenoptera: a cladistic reanalysis of Rasnitsyn’s (1988) data. Zool Scripta28: 13-50. doi:10.1046/j.1463-6409.1999.00023.x. [Google Scholar]
• 15.Vilhelmsen L, Mikó I, Krogmann L (2010) Beyond the wasp-waist: structural diversity and phylogenetic significance of the mesosoma in apocritan wasps (Insecta: Hymenoptera). Zool J Linn Soc159: 22-194. doi:10.1111/j.1096-3642.2009.00576.x. [Google Scholar]
• 16.Heraty JM, Ronquist F, Carpenter JM, Hawks D, Schulmeister Set al. (2011) Evolution of the hymenopteran megaradiation. Mol Phylogenet Evol60: 73-88. doi:10.1016/j.ympev.2011.04.003. PubMed:21540117. [DOI] [PubMed] [Google Scholar]
• 17.Sharkey MJ, Carpenter JM, Vilhelmsen L, Heraty JM, Liljeblad Jet al. (2012) Phylogenetic relationships among superfamilies of Hymenoptera. Cladistics28: 80-112. doi:10.1111/j.1096-0031.2011.00366.x. [DOI] [PubMed] [Google Scholar]
• 18.Lunter GA, Miklós I, Drummond AJ, Jensen JL, Hein J (2005) Bayesian coestimation of phylogeny and sequence alignment. BMC Bioinformatics6: 83. doi:10.1186/1471-2105-6-83. PubMed:15804354. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 19.Redelings BD, Suchard MA (2005) Joint Bayesian estimation of alignment and phylogeny. Syst Biol54: 401-418. doi:10.1080/10635150590947041. PubMed:16012107. [DOI] [PubMed] [Google Scholar]
• 20.Katoh K, Asimenos G, Toh H (2009) Multiple alignment of DNA sequences with MAFFT. In: Posada D.Bioinformatics for DNA Sequence, 537Analysis, 537: 39-64. [DOI] [PubMed] [Google Scholar]
• 21.Suchard MA, Redelings BD (2006) BAli-Phy: simultaneous Bayesian inference of alignment and phylogeny. Bioinformatics22: 2047-2048. doi:10.1093/bioinformatics/btl175. PubMed:16679334. [DOI] [PubMed] [Google Scholar]
• 22.Klopfstein S, Ronquist F (2013) Convergent intron gains in hymenopteran elongation factor-1α. Mol Phylogenet Evol67: 266-276. doi:10.1016/j.ympev.2013.01.015. PubMed:23396205. [DOI] [PubMed] [Google Scholar]
• 23.Tamura K, Peterson D, Peterson N, Stecher G, Nei Met al. (2011) MEGA5: Molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol28: 2731-2739. doi:10.1093/molbev/msr121. PubMed:21546353. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 24.Morrison DA (2009) Why would phylogeneticists ignore computerized sequence alignment?Syst Biol58: 150-158. doi:10.1093/sysbio/syp009. PubMed:20525575. [DOI] [PubMed] [Google Scholar]
• 25.Notredame C (2007) Recent evolutions of multiple sequence alignment algorithms. PLOS Comput Biol3: e123. doi:10.1371/journal.pcbi.0030123. PubMed:17784778. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 26.Munro JB, Heraty JM, Burks RA, Hawks D, Mottern JLet al. (2011) A molecular phylogeny of the Chalcidoidea (Hymenoptera). PLOS ONE6: e27023. doi:10.1371/journal.pone.0027023. PubMed:22087244. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 27.Redelings BD, Suchard MA (2007) Incorporating indel information into phylogeny estimation for rapidly emerging pathogens. BMC Evol Biol7: 40. doi:10.1186/1471-2148-7-40. PubMed:17359539. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 28.Wheeler TJ, Kececioglu JD (2007) Multiple alignment by aligning alignments. Bioinformatics23: i559-i568. doi:10.1093/bioinformatics/btm226. PubMed:17646343. [DOI] [PubMed] [Google Scholar]
• 29.Katoh K, Frith MC (2012) Adding unaligned sequences into an existing alignment using MAFFT and LAST. Bioinformatics, 28: 3144–6. doi:10.1093/bioinformatics/bts578. PubMed:23023983. PubMed:23023983 [DOI] [PMC free article] [PubMed] [Google Scholar]
• 30.Swofford DL (2002) PAUP*. Phylogenetic analysis using parsimony (* and other methods), version 4. Sunderland, MA: Sinauer Associates. [Google Scholar]
• 31.Ronquist F, Teslenko M, Van der Mark P, Ayres DL, Darling Aet al. (2012) MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol61: 539-542. doi:10.1093/sysbio/sys029. PubMed:22357727. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 32.R Development Core Team (2009) R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. [Google Scholar]
• 33.Paradis E, Claude J, Strimmer K (2004) APE: analyses of phylogenetics and evolution in R language. Bioinformatics20: 289-290. doi:10.1093/bioinformatics/btg412. PubMed:14734327. [DOI] [PubMed] [Google Scholar]
• 34.Lewis PO (2001) A likelihood approach to estimating phylogeny from discrete morphological character data. Syst Biol50: 913-925. doi:10.1080/106351501753462876. PubMed:12116640. [DOI] [PubMed] [Google Scholar]
• 35.Stamatakis A, Hoover P, Rougemont J (2008) A rapid bootstrap algorithm for the RAxML web servers. Syst Biol57: 758-771. doi:10.1080/10635150802429642. PubMed:18853362. [DOI] [PubMed] [Google Scholar]
• 36.Wilkinson M (1996) Majority-rule reduced consensus trees and their use in bootstrapping. Mol Biol Evol13: 437-444. doi:10.1093/oxfordjournals.molbev.a025604. PubMed:8742632. [DOI] [PubMed] [Google Scholar]
• 37.Aberer AJ, Krompass D, Stamatakis A(2013) Pruning rogue taxa improves phylogenetic accuracy: an efficient algorithm and webservice. Systematic Biology 62: doi: 10.1093/sysbio/sys1078 (epub ahead of print). [DOI] [PMC free article] [PubMed]
• 38.Dowton M, Austin AD (1995) Increased genetic diversity in mitochondrial genes is correlated with the evolution of parasitism in the Hymenoptera. J Mol Evol41: 958-965. PubMed:8587141. [DOI] [PubMed] [Google Scholar]
• 39.Varón A, Vinh LS, Wheeler WC (2010) POY version 4: phylogenetic analysis using dynamic homologies. Cladistics26: 72-85. doi:10.1111/j.1096-0031.2009.00282.x. [DOI] [PubMed] [Google Scholar]
• 40.Fleissner R, Metzler D, Von Haeseler A (2005) Simultaneous statistical multiple alignment and phylogeny reconstruction. Syst Biol54: 548-561. doi:10.1080/10635150590950371. PubMed:16085574. [DOI] [PubMed] [Google Scholar]
• 41.Saurabh K, Holland BR, Gibb GC, Penny D (2013) Gaps: an elusive source of phylogenetic information. Syst Biol. doi:10.1093/sysbio/sys043. [DOI] [PubMed] [Google Scholar]
• 42.Schulmeister S, Wheeler WC, Carpenter JM (2002) Simultaneous analysis of the basal lineages of Hymenoptera (Insecta) using sensitivity analysis. Cladistics18: 455-484. doi:10.1016/S0748-3007(02)00100-7. [DOI] [PubMed] [Google Scholar]
• 43.Gatesy J, DeSalle R, Wahlberg N (2007) How many genes should a systematist sample? Conflicting insights from a phylogenomic matrix characterized by replicated incongruence. Syst Biol56: 355-363. doi:10.1080/10635150701294733. PubMed:17464890. [DOI] [PubMed] [Google Scholar]
• 44.Goremykin VV, Nikiforova SV, Biggs PJ, Zhong B, Delange Pet al. (2013) The evolutionary root of flowering plants. Syst Biol (epub ahead of print). PubMed:22851550. [DOI] [PubMed] [Google Scholar]
• 45.Wiegmann BM, Trautwein MD, Kim J-W, Cassel BK, Bertone MAet al. (2009) Single-copy nuclear genes resolve the phylogeny of the holometabolous insects. BMC Evol Biol7: 34-. PubMed:19552814. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 46.Meusemann K, von Reumont BM, Simon S, Roeding F, Strauss Set al. (2010) A phylogenomic approach to resolve the arthropod tree of life. Mol Biol Evol27: 2451-2464. doi:10.1093/molbev/msq130. PubMed:20534705. [DOI] [PubMed] [Google Scholar]
• 47.Beutel RG, Vilhelmsen L (2007) Head anatomy of Xyelidae (Hexapoda: Hymenoptera) and phylogenetic implications. Organ Divers Evol7: 207-230. doi:10.1016/j.ode.2006.06.003. [Google Scholar]
• 48.Gibson GAP (1986) Evidence for monophyly and relationships of Chalcidoidea, Mymaridae, and Mymarommatidae (Hymenoptera: Terebrantes). Can Entomol118: 205-240. doi:10.4039/Ent118205-3. [Google Scholar]
• 49.Rasnitsyn AP, Zhang H (2010) Early evolution of Apocrita (Insecta, Hymenoptera) as indicated by new findings in the middle Jurassic of Daohuou, northeast China. Acta Geologica Sinica84: 834-873. doi:10.1111/j.1755-6724.2010.00254.x. [Google Scholar]
• 50.Oeser R (1961) Vergleichend-morphologische Untersuchungen über den Ovipositor der Hymenopteren. Mitteilungen Zool Museum Berl37: 3-119. [Google Scholar]
• 51.Debevec AH, Cardinal S, Danforth BN (2012) Identifying the sister group to the bees: a molecular phylogeny of Aculeata with an emphasis on the superfamily Apoidea. Zool Scripta41: 527-535. doi:10.1111/j.1463-6409.2012.00549.x. [Google Scholar]
• 52.Castro LR, Dowton M (2006) Molecular analyses of the Apocrita (Insecta: Hymenoptera) suggest that the Chalcidoidea are sister to the diaprioid complex. Invert Syst20: 603-614. doi:10.1071/IS06002. [Google Scholar]
• 53.Moulton JK, Wiegmann BM (2004) Evolution and phylogenetic utility of CAD (rudimentary) among Mesozoic-aged Eremoneuran Diptera (Insecta). Mol Phylogenet Evol31: 363-378. doi:10.1016/S1055-7903(03)00284-7. PubMed:15019631. [DOI] [PubMed] [Google Scholar]
• 54.Regier JC, Zwick A (2011) Sources of signal in 62 protein-coding nuclear genes for higher-level phylogenetics of arthropods. PLOS ONE6: e23408. doi:10.1371/journal.pone.0023408. PubMed:21829732. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 55.Philippe H, Brinkmann H, Lavrov DV, Littlewood DTJ, Manuel Met al. (2011) Resolving difficult phylogenetic questions: why more sequences are not enough. PLOS Biol9: e1000602 PubMed:21423652. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 56.Yang Z (1998) On the best evolutionary rate for phylogenetic analysis. Syst Biol47: 125-133. doi:10.1080/106351598261067. PubMed:12064232. [DOI] [PubMed] [Google Scholar]
• 57.Townsend JP (2007) Profiling phylogenetic informativeness. Syst Biol56: 222-231. doi:10.1080/10635150701311362. PubMed:17464879. [DOI] [PubMed] [Google Scholar]
• 58.Klopfstein S, Kropf C, Quicke DLJ (2010) An evaluation of phylogenetic informativeness profiles and the molecular phylogeny of Diplazontinae (Hymenoptera, Ichneumonidae). Syst Biol59: 226-241. doi:10.1093/sysbio/syp105. PubMed:20525632. [DOI] [PubMed] [Google Scholar]
• 59.Rokas A, Williams BL, King N, Carroll SB (2003) Genome-scale approaches to resolving incongruence in molecular phylogenies. Nature425: 798-804. doi:10.1038/nature02053. PubMed:14574403. [DOI] [PubMed] [Google Scholar]
• 60.Collins TM, Fedrigo O, Naylor GJP (2005) Choosing the best genes for the job: the case for stationary genes in genome-scale phylogenetics. Syst Biol54: 493-500. doi:10.1080/10635150590947339. PubMed:16012114. [DOI] [PubMed] [Google Scholar]
• 61.Sharanowski BJ, Robbertse B, Walker JA, Voss SR, Yoder Ret al. (2010) Expressed sequence tags reveal Proctotrupomorpha (minus Chalcidoidea) as sister to Aculeata (Hymenoptera: Insecta). Mol Phylogenet Evol57: 101-112. doi:10.1016/j.ympev.2010.07.006. PubMed:20637293. [DOI] [PubMed] [Google Scholar]
• 62.Peters RS, Meyer B, Krogmann L, Borner J, Meusemann Ket al. (2011) The taming of an impossible child: a standardized all-in approach to the phylogeny of Hymenoptera using public database sequences. BMC Biol9: 55. doi:10.1186/1741-7007-9-55. PubMed:21851592. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 63.Graybeal A (1998) Is it better to add taxa or characters to a difficult phylogenetic problem?Syst Biol47: 9-17. doi:10.1080/106351598260996. PubMed:12064243. [DOI] [PubMed] [Google Scholar]
• 64.Heath TA, Hedtke SM, Hillis DM (2008) Taxon sampling and the accuracy of phylogenetic analyses. J Syst Evolution46: 239-257. [Google Scholar]
• 65.Venditti C, Meade A, Pagel M (2006) Detecting the node-density artifact in phylogeny reconstruction. Syst Biol55: 637-643. doi:10.1080/10635150600865567. PubMed:16969939. [DOI] [PubMed] [Google Scholar]
• 66.Phillips MJ, Delsuc F, Penny D (2004) Genome-scale phylogeny and the detection of systematic biases. Mol Biol Evol21: 1455-1458. doi:10.1093/molbev/msh137. PubMed:15084674. [DOI] [PubMed] [Google Scholar]

Supplementary Materials