Astrid Cruaud

Gérard Delvare

Sabine Nidelet

Laure Sauné

Sujeevan Ratnasingham

Marguerite Chartois

Bonnie B Blaimer

Michael Gates

Seán G Brady

Sariana Faure

Simon Van Noort

Jean-Pierre Rossi

Jean-Yves Rasplus
email: jean-yves.rasplus@inra.fr

Ultra-Conserved Elements and morphology reciprocally illuminate conflicting phylogenetic hypotheses in Chalcididae (Hymenoptera, Chalcidoidea)

Keywords: morphology, systematic bias, topological conflict, Ultra Conserved Elements

Recent technical advances combined with novel computational approaches promised the acceleration of our understanding of the tree of life. However, when it comes to hyperdiverse and poorly known groups of invertebrates, studies are still scarce. As published phylogenies will be rarely challenged by future taxonomists, careful attention must be paid to potential analytical bias. We present the first molecular phylogenetic hypothesis for the family Chalcididae, an emblematic group of parasitoid wasps, with a representative sampling (144 ingroups and 7 outgroups) that covers all described subfamilies and tribes and 82% of the known genera. Analyses of 538 Ultra-Conserved Elements (UCEs) with supermatrix (RAxML and IQTREE) and gene-tree reconciliation approaches (ASTRAL, ASTRID) resulted in highly supported topologies in overall agreement with morphology but reveal conflicting topologies for some of the deepest nodes. To resolve these conflicts, we explored the phylogenetic tree space with clustering and gene genealogy interrogation methods, analyzed marker and taxon properties that could bias inferences and performed a thorough morphological analysis (130 characters encoded for 40 taxa representative of the diversity).

This joint analysis reveals that UCEs enable attainment of resolution between ancestry and convergent /divergent evolution when morphology is not informative enough, but also shows that a systematic exploration of bias with different analytical methods and a careful analysis of morphological features is required to prevent publication of artefactual results. We highlight a GC-content bias for ML approaches, an artefactual mid-point rooting of the ASTRAL tree and a deleterious effect of high percentage of missing data on gene tree reconciliation methods. Based on the results we propose a new classification of the family into eight subfamilies and 10 tribes that lay the foundation for future studies on the evolutionary history of Chalcididae.

.

Chalcididae [START_REF] Delvare | Order Hymenoptera, family Chalcididae[END_REF], the type family of Chalcidoidea that comprises 1,548 described species and 83 genera classified into six subfamilies and seven tribes [START_REF] Noyes | Universal Chalcidoidea Database[END_REF] (Table 1).

The family is found on all continents except the polar regions but has its greatest diversity in the tropics. Species diversity within genera varies greatly, and four genera (Antrocephalus, Brachymeria, Conura, Hockeria) represent more than half (54%) of the species diversity while 65 genera (78%) comprise less than five described species. The evolutionary history of the Chalcididae has been the focus of a single study based on a limited sampling (22 taxa) and a few (34) morphological characters (Wijesekara 1997b). Another morphological study has addressed relationships within two tribes (Wijesekara 1997a). Currently, our knowledge of the infra-familial relationships comes mainly from the three analyses that have focused on the higher level classification of Chalcidoidea using 18S and 28S ribosomal gene regions [START_REF] Campbell | Molecular systematic of the Chalcidoidea using 28S-rDNA[END_REF]; 11 taxa and [START_REF] Munro | A molecular phylogeny of the Chalcidoidea (Hymenoptera)[END_REF]; 41 taxa) or morphology plus rRNA [START_REF] Heraty | A phylogenetic analysis of the megadiverse Chalcidoidea (Hymenoptera)[END_REF]; 25 taxa). In two of the three analyses, the family was retrieved polyphyletic and infrafamilial relationships were never resolved. Although Chalcididae are understudied, they are among the best documented chalcid wasps. Their medium size (1.5 to 15mm) enables the observation of multiple morphological features. Furthermore, characters used in previous studies can serve as a starting point for a thorough morphological analysis.

We chose to infer the molecular phylogeny of the family using Ultra-Conserved Elements (UCEs) and their flanking regions [START_REF] Faircloth | Ultraconserved elements anchor thousands of genetic markers spanning multiple evolutionary timescales[END_REF], McCormack et al. 2012) that are increasingly used to solve ancient and recent divergences in insects [START_REF] Blaimer | Phylogenomic methods outperform traditional multi-locus approaches in resolving deep evolutionary history: a case study of formicine ants[END_REF][START_REF] Faircloth | Target enrichment of ultraconserved elements from arthropods provides a genomic perspective on relationships among Hymenoptera[END_REF], Blaimer et al. 2016a[START_REF] Bossert | The impact of GC bias on phylogenetic accuracy using targeted enrichment phylogenomic data[END_REF], Branstetter et al. 2017a, Branstetter et al. 2017b[START_REF] Jesovnik | Phylogenomic species delimitation and host-symbiont coevolution in the fungus-farming ant genus Sericomyrmex Mayr (Hymenoptera: Formicidae): ultraconserved elements (UCEs) resolve a recent radiation[END_REF][START_REF] Prebus | Insights into the evolution, biogeography and natural history of the acorn ants, genus Temnothorax Mayr (Hymenoptera: Formicidae)[END_REF][START_REF] Van Dam | Ultraconserved elements (UCEs) resolve the phylogeny of Australasian smurfweevils[END_REF], Ward and Branstetter 2017[START_REF] Bossert | Combining transcriptomes and ultraconserved elements to illuminate the phylogeny of Apidae[END_REF][START_REF] Cruaud | Optimised DNA extraction and library preparation for minute arthropods: application to target enrichment in chalcid wasps used for biocontrol[END_REF][START_REF] Kieran | Insight from an ultraconserved element bait set designed for hemipteran phylogenetics integrated with genomic resources[END_REF].

We sequenced UCEs from 144 ingroup and 7 outgroup taxa, analyzed UCE /taxa properties, performed exploration of the phylogenetic tree space and used different analytical approaches to detect possible systematic bias and conflicts among gene regions. We interpreted results in the light of a thorough morphological analysis (130 characters encoded on 40 taxa representing all major lineages) to propose a new classification for the higher relationships within the family.

MATERIALS AND METHODS

Sampling for the molecular study

The data set contained 144 ingroup taxa (Table S1). All subfamilies and tribes as well as 82% of the world genera were included in the analysis. At the beginning of our study, 83 genera of Chalcididae were considered as valid and two as incertae sedis (Antrochalcis Kieffer, 1910 andChalcitiscus Ghesquière, 1946, a fossil genus). The within tribe classification was not the purpose of this study and will be reviewed elsewhere. However, the examination of numerous specimens in several museums and personal collections, together with a review of the literature, suggested that we need to remove six genera from synonymy (Eniacomorpha Girault, Hontalia Cameron, Invreia Masi, Pareniaca Crawford, Parinvreia Steffan, Peltochalcidia Steffan). We also discovered seven new genera awaiting description. Thus, we now consider 94 genera of Chalcididae as valid, 77 of which were included in our study. Most missing genera are extremely rare, known mostly from the type series only. Relationships within the Chalcidoidea are unclear [START_REF] Heraty | A phylogenetic analysis of the megadiverse Chalcidoidea (Hymenoptera)[END_REF]) but morphological features [START_REF] Lotfalizadeh | Phylogenetic analysis of Eurytominae based on morphological characters (Chalcidoidea: Eurytomidae)[END_REF]) as well as preliminary results obtained with anchored hybrid enrichment and UCEs (work in prep) suggest that the Eurytomidae may be the sister family of the Chalcididae. As a consequence, five species of Eurytomidae representing three subfamilies, as well as two species of Cerocephalinae (that form a clade with Eurytomidae + Chalcididae in the chalcidoid tree from the anchored/UCE approach) were used as outgroups.

We used dried specimens (35% of the samples), specimens kept in 75-96% EtOH (50%), as well as DNA extracts (15%) kept at -20°C for about the last 10 years. The oldest specimen was collected in 1951 though most specimens were collected in the last 20 years. The paratypes of two species: Kopinata partirubra Bouček, 1988 andChalcis vera Bouček, 1974 housed in NHMUK, London were also included.

DNA extraction, library preparation and sequencing

DNA was extracted non-destructively and vouchers were subsequently remounted on cards.

DNA was extracted using the Qiagen DNeasy Blood and Tissue kit following manufacturer's protocol with a few modifications detailed in [START_REF] Cruaud | Optimised DNA extraction and library preparation for minute arthropods: application to target enrichment in chalcid wasps used for biocontrol[END_REF]. Library preparation followed [START_REF] Cruaud | Optimised DNA extraction and library preparation for minute arthropods: application to target enrichment in chalcid wasps used for biocontrol[END_REF]. Briefly, input DNA was sheared to a size of ca 400 bp using the Bioruptor® Pico (Diagenode). End repair, 3'-end adenylation, adapters ligation and PCR enrichment were then performed with the NEBNext Ultra II DNA Library prep kit for Illumina (NEB). Adapters that contained amplification and Illumina sequencing primer sites, as well as a nucleotide barcode of 5 or 6 bp long for sample identification were used to tag samples. Pools of 16 samples were made at equimolar ratio. Each pool was enriched using the 2749 probes designed by [START_REF] Faircloth | Target enrichment of ultraconserved elements from arthropods provides a genomic perspective on relationships among Hymenoptera[END_REF] using a MYbaits kit (Arbor Biosciences) and following manufacturer's protocol. The hybridization reaction was run for 24h at 65°C. Post enrichment amplification was performed on beads with the KAPA Hifi HotStart ReadyMix.

The enriched libraries were quantified with Qubit, an Agilent Bioanalizer and qPCR with the Library Quantification Kit -Illumina/Universal from KAPA (KK4824). They were then pooled at equimolar ratio. Paired-end sequencing (2*300bp) was performed on an Illumina Miseq platform at UMR AGAP (Montpellier, France).

UCE data analysis (from raw reads to UCEs)

Quality control checks were performed on raw sequence data with FastQC v.0.11.2 [START_REF] Andrews | FastQC -A Quality Control application for FastQ files[END_REF]. Quality filtering and adapter trimming were performed with Trimmomatic-0.36 [START_REF] Bolger | Trimmomatic: a flexible trimmer for Illumina sequence data[END_REF]. Overlapping reads were merged using FLASH-1.2.11 [START_REF] Magoc | FLASH: fast length adjustment of short reads to improve genome assemblies[END_REF]. Demultiplexing was performed using a bash custom script [START_REF] Cruaud | Optimised DNA extraction and library preparation for minute arthropods: application to target enrichment in chalcid wasps used for biocontrol[END_REF]. Assembly was performed with Trinity [START_REF] Haas | De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[END_REF]. UCE loci were identified with PhylUCE [START_REF] Faircloth | PHYLUCE is a software package for the analysis of conserved genomic loci[END_REF]) (phyluce_assembly_match_contigs_to_probes, phyluce_assembly_get_match_counts --incomplete matrix, phyluce_assembly_get_fastas_from_match_counts, all scripts were used with default parameters).

UCE data analysis (from UCEs to phylogenetic trees)

UCEs present in more than 70% of the taxa were retained for analysis. UCEs were aligned with MAFFT v7.245 (-linsi option) [START_REF] Katoh | MAFFT multiple sequence alignment software version 7: improvements in performance and usability[END_REF]. Sites with more than 50% gaps were removed from the alignments using the program seqtools implemented in the package PASTA (Mirarab et al. 2014b). Sequences exhibiting up to 500 gaps and 10 substitutions in pairwise alignments between all members of a UCE set were flagged by a custom script (available from https://github.com/DNAdiversity/UCE-Cross-Contamination-Check) as potential contaminations and reviewed before exclusion in subsequent analysis.

Individual gene trees were inferred with raxmlHPC-PTHREADS-AVX version 8.2.4 [START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF]. As α and the proportion of invariable sites cannot be optimized independently from each other [START_REF] Gu | Maximum likelihood estimation of the heterogeneity of substitution rate among nucleotide sites[END_REF] and following Stamatakis' personal recommendation (RAxML manual), the proportion of invariant sites was not included in the model. A rapid bootstrap search (100 replicates) followed by a thorough ML search (-m GTRGAMMA) was performed. TreeShrink [START_REF] Mai | TreeShrink: fast and accurate detection of outlier long branches in collections of phylogenetic trees[END_REF] was used to detect and remove abnormally long branches in individual gene trees (e.g. due to misalignment). The per-species mode was used and preliminary analyses together with gene tree visualization showed that the optimal value of b (the percentage of tree diameter increasing from which a species should be removed) was 20. Indeed, as reported on TreeShrink tutorial, lower value of b led to the removal of species even when they were not on particularly long branches. All other parameters were set to default values. Once outliers were removed, UCEs were realigned with MAFFT-linsi and alignments were cleaned using seqtools. TreeShrink was used a second time to clean gene trees from long branches that might have been missed due to the presence of extra-long branches in the original gene trees. The per-species mode was used with b set to 20 but only the longest outliers were removed (k was set to 1). All other parameters were set to default values. Once outliers were removed, UCEs were re-aligned with MAFFT-linsi and alignments were cleaned using seqtools.

The final data set was analyzed using supermatrix and coalescent-based summary methods.

Phylogenetic trees were estimated from the concatenated data set using Maximum Likelihood (ML) as implemented in raxmlHPC-PTHREADS-AVX version 8.2.4 [START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF]IQTREE v1.5.3 (Nguyen et al. 2015). For the RAxML analysis, a rapid bootstrap search (100 replicates) followed by a thorough ML search (-m GTRGAMMA) was performed. IQTREE analysis employed an ML search with the GTR+G model with branch supports assessed with ultrafast bootstrap [START_REF] Minh | Ultrafast approximation for phylogenetic bootstrap[END_REF]) and SH-aLRT test [START_REF] Guindon | New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0[END_REF]) (1000 replicates). For both RAxML and IQTREE, two analyses were conducted: i) on the unpartitioned data set, ii) on the dataset partitioned according to the best partitioning scheme inferred by PartitionFinder 2.1.1 [START_REF] Lanfear | PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses[END_REF]) using the Sliding-Window Site Characteristics (SWSC) method. This method has been recently developed to account for within-UCE heterogeneity (conserved core versus flanking, variable regions) [START_REF] Tagliacollo | Estimating improved partitioning schemes for Ultraconserved Elements[END_REF]). To fit with RAxML models and as α and the proportion of invariable sites cannot be optimized independently from each other, only the GTR+G model was evaluated. Branch lengths were considered as linked, model selection and partitioning scheme comparison were performed with the corrected Akaike Information Criterion (AICc) and the rclusterf algorithm. Finally, we used the GHOST model implemented in IQTREE to account for heterotachous evolution as it does not require a priori data partitioning, a possible source of model misspecification [START_REF] Crotty | GHOST: Recovering historical signal from heterotachously-evolved sequence alignments[END_REF]. ASTRAL-III v5.6.1 (Zhang et al. 2018) and ASTRID [START_REF] Vachaspati | ASTRID: Accurate Species TRees from Internode Distances[END_REF] were used to infer a species tree from the individual UCE trees inferred by RAxML. To improve accuracy (Zhang et al. 2018) nodes with BP support < 10 were collapsed in individual gene trees with the perl script AfterPhylo.pl (Zhu 2014). For the analysis with ASTRID and following recommendations for incomplete distance matrices, BioNJ was used to compute the phylogeny. Node supports were evaluated with local posterior probabilities (local PP) for the ASTRAL tree and 100 multi-locus bootstrapping [START_REF] Seo | Calculating bootstrap probabilities of phylogeny using multilocus sequence data[END_REF] for the ASTRID tree. Summary statistics were calculated using AMAS [START_REF] Borowiec | AMAS: a fast tool for alignment manipulation and computing of summary statistics[END_REF]. Tree annotation was performed with TreeGraph 2.13 [START_REF] Stöver | TreeGraph 2: Combining and visualizing evidence from different phylogenetic analyses[END_REF]. Correlation analysis between properties of the UCEs was performed with the R package Performance Analytics [START_REF] Peterson | PerformanceAnalytics: econometric tools for performance and risk analysis[END_REF].

Exploration of topological conflicts

New approaches have been recently developed to identify markers /sites supporting conflicting topologies that either make explicit assumptions about the biological basis of conflict (e.g. horizontal transfer, incomplete lineage sorting ILS, recombination, gene duplication) or not. Early approaches (e.g. [START_REF] Abby | Detecting lateral gene transfers by statistical reconciliation of phylogenetic forests[END_REF][START_REF] Heled | Bayesian inference of species trees from multilocus data[END_REF][START_REF] Szöllősi | Modeling gene family evolution and reconciling phylogenetic discord[END_REF] were computationally too intense to be implemented on large data sets. Furthermore, they are limited in their focus, constrained to one or two sources of incongruence and may not be robust to additional sources [START_REF] Gori | Clustering genes of common evolutionary history[END_REF]). Here we used more recent approaches, that do not rely on any assumptions about the biological basis of conflicts, and are computationally tractable on large data sets. This involves: partitioning of the data into coherent groups by clustering of tree to tree distances [START_REF] Gori | Clustering genes of common evolutionary history[END_REF][START_REF] Duchêne | Analysis of phylogenomic tree space resolves relationships among marsupial families[END_REF] or statistical tests of incongruence (Gene Genealogy Interrogation GGI; [START_REF] Arcila | Genome-wide interrogation advances resolution of recalcitrant groups in the tree of life[END_REF], Zhong and Betancur-R 2017[START_REF] Betancur-R R | Phylogenomic incongruence, hypothesis testing, and taxonomic sampling: The monophyly of characiform fishes[END_REF]. For the clustering approach, geodesic distances between all pairs of trees were calculated with TreeCl [START_REF] Gori | Clustering genes of common evolutionary history[END_REF] that requires overlapping set of samples but allows missing samples. Then, the optimal number of clusters obtained with the Partitioning Around Medoids (PAM) algorithm [START_REF] Kaufman | Partitioning around medoids (Program PAM)[END_REF] was estimated with the gap statistics as implemented in the R package cluster [START_REF] Maechler | cluster: Cluster Analysis Basics and Extensions[END_REF][START_REF] Core | R version 3.5.1 (Feather Spray): A language and environment for statistical computing[END_REF] (Kmax was set to 10 and number of bootstrap samples was set to 500 to keep computation time reasonable). The R package factoextra [START_REF] Kassambara | factoextra: extract and visualize the results of multivariate data analyses[END_REF] was used to visualize clusters (ggplot2based scatter plot [START_REF] Wickham | ggplot2: Elegant graphics for data analysis[END_REF]. Possible effects of missing data or gap content were evaluated by testing whether these variables were phylogenetically clustered on any of the conflicting topology. Tests were conducted using the K statistic [START_REF] Blomberg | Testing for phylogenetic signal in comparative data: behavioral traits are more labile[END_REF] with the R package Phytools 0.6-00 [START_REF] Revell | Phytools: An R package for phylogenetic comparative biology (and other things)[END_REF]. The null expectation of K under no phylogenetic signal was generated by randomly shuffling the tips of the phylogeny 1000 times.

We also used Gene Genealogy Interrogation (GGI) to compute the relative support of the UCEs for each competing topology following [START_REF] Arcila | Genome-wide interrogation advances resolution of recalcitrant groups in the tree of life[END_REF]. RAxML was used to infer trees from UCEs using each of the three competing topologies (Figure 1) as multi-furcating constraints (the structure of the backbone was fixed, but taxa within clades were free to move around). Thus, three constrained trees were inferred from each UCE. Per-site log likelihood scores for all constrained trees were calculated with RAxML and used to perform AU tests [START_REF] Shimodaira | An approximately unbiased test of phylogenetic tree selection[END_REF] in the package CONSEL [START_REF] Shimodaira | CONSEL: for assessing the confidence of phylogenetic tree selection[END_REF]. The program makermt was used to generate K = 10 sets of bootstrap replicates with each set consisted of 100 000 replicates of the row sums. For each UCE, constrained trees were ranked by p-values of the AU test. The constrained tree with the highest p-value was considered as the best explanation of the data.

Finally, we examined whether compositional heterogeneity among loci and nucleotide positions as well as evolutionary rate heterogeneity among taxa could explain topological conflict. GC content and long branch (LB) score heterogeneity were calculated for each UCE and each taxon in all UCEs. GC content was calculated with AMAS [START_REF] Borowiec | AMAS: a fast tool for alignment manipulation and computing of summary statistics[END_REF]) and LB heterogeneity scores were calculated with TreSpEx [START_REF] Struck | TreSpEx -Detection of misleading signal in phylogenetic reconstructions based on tree information[END_REF]. In a given tree, the per sample LB score measures the percentage deviation of the patristics distance (PD) of a sample to all others, from the average PD across all samples [START_REF] Struck | TreSpEx -Detection of misleading signal in phylogenetic reconstructions based on tree information[END_REF]. For a given tree, the LB score heterogeneity is the standard deviation of the LB scores of the samples present in the tree. Thus, the LB score heterogeneity reflects the branch length heterogeneity of a given tree and is independent of the root of the tree. Consequently, calculation of LB score heterogeneities can help to prevent possible Long Branch Attraction (LBA) artefacts [START_REF] Phillips | Genome-scale phylogeny and the detection of systematic biases[END_REF][START_REF] Bergsten | A review of long-branch attraction[END_REF]. Hierarchical clustering of taxa based on GC content and LB scores was performed with the R package cluster. The LBA artefact was also tested by removing outgroups from the analysis [START_REF] Bergsten | A review of long-branch attraction[END_REF].

Morphological analyses

A matrix of 130 characters (Appendix S1) was assembled for 40 species that covered all major lineages of the Chalcididae (Appendix S2). The matrix was analyzed with parsimony (PAUP* version 4.0a; Swofford 2003), maximum likelihood (RAxML) and Bayesian (MrBayes-3.2.6; [START_REF] Ronquist | MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space[END_REF] approaches. PAUP analyses were performed with unordered, equally weighted and non-additive characters. A traditional heuristic search was conducted using 1000 random addition sequences (RAS) to obtain an initial tree and "tree bisection and reconnection (TBR)" as branch swapping option. Ten trees were retained at each step. Robustness of the topology was assessed by bootstrap procedures (100 replicates; TBR RAS 10; one tree retained at each step). RAxML and MrBayes analyses were conducted with the Mk model [START_REF] Lewis | A likelihood approach to estimating phylogeny from discrete morphological character data[END_REF] with only variable characters scored, equal state frequencies and assuming a gamma-distributed rate variation across characters. Robustness of the ML topology was assessed by bootstrap (100 replicates). Two independent runs of 1 million generations were run for the MrBayes analysis, with each run including a cold chain and three incrementally heated chains. The heating parameter was set to 0.02 in order to allow swap frequencies from 20% to 70%. Parameter values were sampled every 1000 generations.

Convergence was examined in Tracer 1.6 [START_REF] Rambaut | Tracer v1.6[END_REF]) and the results were based on the pooled samples from the stationary phases of the 2 independent runs.

Computational resources

Analyses were performed on a Dell PowerEdge T630 server with two 10-core Intel(R) Xeon(R) CPUs E5-2687W v3 @ 3.10GHz and on the Genotoul Cluster (INRA, Toulouse).

RESULTS

UCE data set

The final data set included 151 taxa (Table S1). No cross-contaminations were detected. In the first round, TreeShrink detected outlier long branches in 155 gene trees. Between 1 to 12 samples were flagged and removed per gene tree (average = 2 samples; Table S2). In the second round, TreeShrink detected outlier long branches in 28 gene trees. The final matrix (70% complete) included 538 UCEs and taxa were represented by 19 to 528 UCEs (median 478, Table S1). Eight taxa had more than 80% missing UCEs. The alignment contained 283,634 bp, 70.6% of which were parsimony informative. The percentage of missing data was 20.1%, the percentage of gaps was 16.1% and the GC content was 40.9%. Gap content for taxa (that can either result from the alignment of full length UCEs or capture of partial UCEs) ranged from to 2.2 to 37.0% (median 13.3%).

Phylogenetic inference from the UCE data set

PartitionFinder2 used in combination with the SWSC method split the data set into 890 partitions. ML (Figures S1-S2, Appendix S3), ASTRAL (Figure S3) and ASTRID (Figure S4) trees were globally well resolved. The normalized quartet score of the ASTRAL tree was 0.91, which indicates a high degree of congruence between the species-tree and the input gene Notably, these conflicts were still observed when nodes with BP < 50 were collapsed in individual gene trees prior to species tree inference by ASTRAL (Figure S3; normalized quartet score = 0.98). For brevity, we will hereafter refer to the Brachymeriini + Dirhininae + Epitraninae + Phasgonophorini + Smicromorphinae clade (Figure 1) as the "BDEPS" clade.

Exploration of topological conflicts

Analysis of bootstrap ML trees showed that topology A was recovered in ca 80% of the replicates and topology C was recovered in the remaining 20%. Topology B was recovered when ML trees were mid-point rooted (Figure S1C). Visualization of individual UCE trees revealed that ingroups were monophyletic in only 51 of the 538 trees (i.e 9.48%, Appendix S3). Interestingly, when the set of input trees for the ASTRAL analysis was reduced to these 51 trees, ASTRAL inferred topology C (Figure 3C). Neither missing data nor gap content appeared phylogenetically clustered on the ML, ASTRAL or ASTRID topologies (p-values > 0.05). The optimal number of clusters of loci as estimated by the gap statistics on the matrix of geodesic distances among individual gene trees was 1 (Figure S5) and the GGI approach showed that none of the topology was significantly preferred over the other (Figure 2).

Topology A had a significantly best fit for 32 UCEs according to the GGI approach (p-value of the AU-test < 0.05).

Spearman's rank correlation tests showed a significant negative correlation between GC content of the UCEs (Table S3) and the average support of individual gene trees (Figure S6).

There was a higher difference between observed base composition and that predicted under the substitution model for GC-rich UCEs. The alpha parameter of the Gamma distribution was positively correlated with the number of sites informative for the parsimony (i.e. UCEs with more homogeneous rates among sites are more informative) and the average support of individual gene trees. LB score heterogeneity (Table S3) was negatively correlated with the number of sites informative for the parsimony and the average support of individual gene trees.

Hierarchical clustering of taxa based on the heterogeneity of LB scores (Table S4, Figure S7A) suggested that, with a few exceptions (Zavoyini, Tropimeridini, Belaspidia), the diversification dynamics of the Haltichellinae was somehow different compared to other subfamilies. Visual observation of the branch lengths of the ML trees confirmed this pattern (Figures S1-S2). Brachymeriini, Chalcidini, Cratocentrinae, Dirhininae, Epitraninae, Phasgonophorini and Smicromorphinae are supported by more elongated external branches.

However, the dendrogram did not show obvious evidence to suspect that the position of Cratocentrinae as sister to the other Chalcididae in the topology A could result from an LBA artefact. In the dendrogram built from GC content (Table S4, Figure S7B), outgroups were not monophyletic and the Haltichellinae split into two groups. One of the subgroups shared similar properties as a few species of Chalcidini and Phasgonophorini.

Outgroup removal did not result in a shift of position for Cratocentrinae and the GHOST tree that accounts for heterotachous evolution was identical to other ML trees (Appendix S3).

Joint analyses (GGI+UCE properties) revealed that heterogeneities of LB scores were not significantly different among the UCEs that supported either topology (Figure 2B-C). On the contrary, UCEs supporting topology C (with p-values significant or not) had a significantly lower GC content (Figure 2B). Furthermore, the 32 UCEs that significantly supported topology A had a significantly higher GC content (Figure 2C).

These results suggest that GC content may bias the results towards topology A. To test this hypothesis, data subsets were constructed by incrementally removing the most biased UCEs and nucleotide sites (Tables 2&3). To keep computation time reasonable, ML analyses were performed only with RAxML without partitioning. Indeed, all other approaches /models /softwares infer the same topology from the complete data set and are certainly subject to the same bias, if any. While ASTRAL and ASTRID inferred the same topology whatever the data subset analyzed, RAxML inferred topology C instead of topology A when UCEs with GC content > 0.48 (i.e. 28% of the UCEs) were removed from the data set (Table 2,Appendix S3). However, as already observed with the complete data set, bootstrap support for position of Chalcidini never reached 100%. The removal of nucleotide sites with GC > 0.57 (i.e 18.9% of the sites) also induced a shift from topology A to topology C for the ML analyses (Table 3, Appendix S3), though here again, the position of Chalcidini was not supported by 100% bootstrap (BP=80).

Analysis of morphological data.

This study is the first in the whole superfamily Chalcidoidea to investigate, among other characters, the inner skeleton of the cephalic capsule and its external landmarks on the back of the head. Detailed results about character coding in the different taxa are provided in Appendix S4. Parsimony analysis of the 130 morphological characters (Appendix S1, S2) produced 1439 equally parsimonious trees (length = 330, consistency index (CI) = 0.545, retention index (RI) = 0.840, rescaled consistency index (RC) = 0.458; Figure 3). ML and Bayesian analyses recovered Cratocentrinae as sister to Phasgonophorini with moderate support (Figure S8). Otherwise, the ML, Bayesian and parsimony topologies were similar.

Unsurprisingly, trees inferred from the morphological matrix were poorly supported as compared to the UCE trees. Just as UCE trees, morphological data support the monophyly of the family, the monophyly of all subfamilies except Chalcidinae, and the polyphyly of the tribes Haltichellini and Hybothoracini. In contrast with UCE trees, Zavoya was recovered sister to Notaspidium with high support and Dirhininae + Epitraninae formed a strongly supported clade. Tree backbones were too poorly supported to draw any firm conclusions. An in-depth analysis of the character states that support each clade is provided in Appendix S5.

Mapping of character state transformation revealed that no one morphological character supported a sister taxa relationship between Chalcidini and Haltichellinae as observed in the UCE ML trees. Five characters unambiguously supported Cratocentrinae as sister to all other Chalcididae as observed in the UCE ML+ASTRID trees: mandibular base exposed, condyles elongate and visible externally, mouth margin not incised for reception of mandible (lateral to clypeus); subforaminal bridge depressed compared to postgena; lateral lamella on anterior tentorial arm narrow but with broad apical lobe; single metafurcal pit medially; two parallel and short submedian carinae present between metacoxae. Only two characters unambiguously supported a sister taxa relationship between Cratocentrinae, Chalcidini and the BDEPS clade as observed in the UCE ASTRAL tree: no raised carina on the inner margin of axillula and apex of metatibia diagonally truncate (Appendix S5). Consequently, morphological data provided no support for topology B and provides more support for topology C over topology A.

DISCUSSION

The power of UCEs to resolve the tree of life of poorly known groups So far, UCEs have been successfully used to infer the phylogenies of a few groups of large to medium-sized insects (bees, wasps, ants and weevils: [START_REF] Blaimer | Phylogenomic methods outperform traditional multi-locus approaches in resolving deep evolutionary history: a case study of formicine ants[END_REF], Blaimer et al. 2016a, Branstetter et al. 2017a[START_REF] Jesovnik | Phylogenomic species delimitation and host-symbiont coevolution in the fungus-farming ant genus Sericomyrmex Mayr (Hymenoptera: Formicidae): ultraconserved elements (UCEs) resolve a recent radiation[END_REF][START_REF] Prebus | Insights into the evolution, biogeography and natural history of the acorn ants, genus Temnothorax Mayr (Hymenoptera: Formicidae)[END_REF][START_REF] Van Dam | Ultraconserved elements (UCEs) resolve the phylogeny of Australasian smurfweevils[END_REF][START_REF] Bossert | Combining transcriptomes and ultraconserved elements to illuminate the phylogeny of Apidae[END_REF]). Here, we highlight the power of UCEs for the exploration of hyperdiverse and poorly known groups of small to medium-sized chalcid wasps.

Interestingly, we were able to successfully use the universal Hymenoptera probe set designed by [START_REF] Faircloth | Target enrichment of ultraconserved elements from arthropods provides a genomic perspective on relationships among Hymenoptera[END_REF] to capture UCEs from chalcid wasps without any optimisation. This confirms the genericness of UCEs [START_REF] Bossert | On the universality of target-enrichment baits for phylogenomic research[END_REF]. Besides Hymenoptera, universal probe sets have been designed for Arachnida, Coleoptera, Diptera, Hemiptera and Lepidoptera [START_REF] Faircloth | Identifying conserved genomic elements and designing universal bait sets to enrich them[END_REF]) that could likely be used to infer relationships within many other poorly known insect groups. Easy to set-up and affordable protocols are freely available (http://ultraconserved.org; [START_REF] Faircloth | Target enrichment of ultraconserved elements from arthropods provides a genomic perspective on relationships among Hymenoptera[END_REF][START_REF] Cruaud | Optimised DNA extraction and library preparation for minute arthropods: application to target enrichment in chalcid wasps used for biocontrol[END_REF], that do not require advanced skills or the use of an expensive robotic molecular biology platform. A software with a detailed tutorial [START_REF] Faircloth | PHYLUCE is a software package for the analysis of conserved genomic loci[END_REF]) exists for the easy analysis of raw data.

Three situations can be distinguished in our study that may be encountered in other groups.

First, when morphology is informative enough to circumscribe taxa and establish relationships, UCEs and morphology converge to the same results. Second, when morphology is not informative enough or misleading, UCEs are helpful to circumscribe taxa (e.g. genera or tribes) and clarify relationships within taxa. Additional studies of the morphology can then be performed and help discern characters that support the taxa. In summary, UCEs help to enable attainment of resolution between ancestry, convergent evolution or divergent evolution. We confirm that UCEs can be captured from museum specimens (Blaimer et al. 2016b[START_REF] Mccormack | Sequence capture of ultraconserved elements from bird museum specimens[END_REF], which is key requirement when working on poorly known/rare groups, as sequencing of type specimens is often required to fix species names.

Lastly, when neither morphology nor molecules are informative enough to resolve relationships, a systematic exploration of bias as well as the use of different analytical methods and a careful feedback assessment with morphological features is required.

Resolving conflicts among methodological approaches

While there is a global agreement between morphology and molecules on the one hand, and between the different UCE analytical approaches on the other hand, there are a few conflicts that we discuss below.

First, there are unsupported conflicts between concatenation and gene tree reconciliation approaches for the position of taxa with a high level of missing UCEs (>90% e.g Hastius Schmitz, Solenochalcidia Steffan, Pseudeniaca Masi). The impact of missing data on phylogenetic inference has been widely discussed but no consensus has been reached.

Missing data are either considered as deleterious [START_REF] Lemmon | The effect of ambiguous data on phylogenetic estimates obtained by maximum likelihood and Bayesian inference[END_REF] or not problematic if enough informative characters are available to infer relationships [START_REF] Wiens | Missing data, incomplete taxa, and phylogenetic accuracy[END_REF][START_REF] Wiens | Missing data in phylogenetic analysis: reconciling results from simulations and empirical data[END_REF][START_REF] Hosner | Avoiding missing data biases in phylogenomic inference: an empirical study in the landfowl (Aves: Galliformes)[END_REF][START_REF] Streicher | How should genes and taxa be sampled for phylogenomic analyses with missing data? an empirical study in iguanian lizards[END_REF]. Gene tree reconciliation approaches have been proven robust enough to a high, global level, of missing data [START_REF] Nute | The performance of coalescentbasedspecies tree estimation methods under modelsof missing data[END_REF].

However, according to our knowledge of their morphology, the position of taxa with a high level of missing data is more accurate in the concatenation approach than in the gene tree reconciliation approaches.

Second, while shallow and intermediate relationships are similar among analytical approaches, two strongly supported conflicts are highlighted for the deepest nodes of the UCE phylogeny (Figure 1): 1) Cratocentrinae is recovered either as sister to all other Chalcididae (ML + ASTRID) or nested within Chalcididae (ASTRAL).

2) Chalcidini is either sister to Haltichellinae (ML), sister to the BDEPS clade (ASTRID) or sister to the BDEPS + Cratocentrinae clade (ASTRAL).

Spearman correlation tests suggest that evolutionary rate heterogeneity among taxa and nucleotide sites as well as compositional heterogeneity among UCEs could bias the analyses (Figure S6). Thus, the first hypothesis we explored was that the position of Cratocentrinae as sister to all other chalcidids in the ML tree resulted from an LBA artefact. Indeed, supermatrix approaches are more sensitive to LBA which, in addition, tend to be reinforced as more and more markers are considered [START_REF] Boussau | Strepsiptera, Phylogenomics and the Long Branch Attraction Problem[END_REF]). However, hierarchical clustering of taxa properties (Figure S7), analysis with complex models that considers heterotachous evolution (GHOST, Appendix S3), outgroup removal analysis (Appendix S3) and the non-significant difference of scores of LB heterogeneity among UCEs that support topology A and other UCEs (Figure 2), show that LBA should be excluded. Sampling used in this study could certainly be improved but is representative of the group. More importantly, the outgroups used for this study are the closest relatives to the ingroups (work in prep., see the method section). Finally, morphological analysis reveals that Chalcididae minus Cratocentrinae possess apomorphic characters (Figure 3, Appendix S5) for which a scenario of loss and reacquisition implied by the ASTRAL tree seems unlikely.

On the other hand, our results suggest that the position of Chalcidini as sister to Haltichellinae in the ML trees, could result from a GC content bias. Indeed, in the morphological analysis, Chalcidini was never recovered sister to Haltichellinae. Besides, Chalcidini does not share any uniquely derived characters with Haltichellinae (Figure 3, Appendix S5). The largest difference between the observed GC content and that predicted under the substitution model is obtained for GC-rich UCEs and GGI analyses show that topology A is preferred by UCEs with a significantly higher GC content (Figure 2). Furthermore, hierarchical clustering of taxa based on their GC content show that some Haltichellinae share more similar GC content with Chalcidini or other subfamilies than with members of their own subfamily. GC-rich markers are more subject to recombination (Lartillot 2013, Romiguier and[START_REF] Romiguier | Analytical biases associated with GC-content in molecular evolution[END_REF] and, as a consequence, to ILS [START_REF] Pease | More accurate phylogenies inferred from low recombination regions in the presence of incomplete lineage sorting[END_REF]. Thus, it is not surprising that only methods that are statistically consistent under the multi-species coalescent model (i.e. ILS-aware methods ASTRAL, ASTRID) do not recover Chalcidini sister to Haltichellinae from the analysis of the complete data set. ILS-aware methods have indeed been developed primarily to solve the deepest relationships for rapid radiations (e.g. birds (Mirarab et al. 2014a, Mirarab et al. 2014c). For the concatenation approach, when the most GC-biased UCEs (GC content >0.48; ca 28% of the total UCEs which corresponds to 23.7% of the total sites) or nucleotide sites (GC content > 0.57; 18.9% of the total sites) were removed, a sister taxa relationship between Chalcidini and the BDEPS clade was inferred, though with moderate bootstrap support (Tables 23Figure 4,Appendix S3). This result agrees with other studies on reduced data sets, which revealed that UCEs can be GC-biased and support conflicting topologies [START_REF] Sun | The evolution of peafowl and other taxa with ocelli (eyespots): a phylogenomic approach[END_REF][START_REF] Bossert | The impact of GC bias on phylogenetic accuracy using targeted enrichment phylogenomic data[END_REF]). However, we encourage the removal of the most GC biased nucleotide sites instead of the full UCEs to preserve phylogenetic signal in the analyzed subset.

That said, one question still remains unanswered: how can we explain that Cratocentrinae are nested within Chalcididae in the ASTRAL tree? It is difficult to draw firm conclusions regarding this point. Our analyses suggest that individual UCE trees are not sufficiently resolved to allow a proper inference on the position of the outgroups from all bipartitions present in the input gene trees. Indeed, a majority of trees with intermixed outgroups and ingroups can complicate species tree inferences [START_REF] Mai | Minimum variance rooting of phylogenetic trees and implications for species tree reconstruction[END_REF]. The consequence of this lack of information is that outgroups are placed in an intermediate position. Indeed, the position of the outgroups in the ASTRAL tree corresponds to a mid-point rooting of the ML trees (Figure S1C) and topology C is recovered by ASTRAL when the set of input trees is reduced to those for which the ingroup is monophyletic. Outgroup choice is a difficult decision. Ideally, outgroups should be as close as possible to the ingroups to reduce LBA artefact while remaining sufficiently distantly-related to reduce impacts of ILS. However, an informed choice is often difficult if not impossible for hyperdiverse and poorly known groups.

Bringing together all species for phylogenetic inference is also impossible and such large data sets would be impossible to analyse with current methods [START_REF] De Vienne | Pitfalls in supermatrix phylogenomics[END_REF]. This is why exploration of phylogenetic incongruence and systematic bias is particularly needed. The results of the ASTRAL analysis on the complete set of gene trees could also be explained by the nature of the marker used. UCEs are short which leads to poorly resolved gene trees with unsupported, incorrect bipartitions. Alignment cleaning (here indel removal) contributes to further reduce UCE length and, to some extent, to a loss of signal that should be quantified.

Further research is needed to better understand possible drawback of tree reconciliation methods based on gene tree topology when they are used on UCEs and what could be the best strategy for alignment cleaning to preserve signal contained in gaps [START_REF] Donath | Split-inducing indels in phylogenomic analysis[END_REF].

When it comes to the analysis of ancient groups that have undergone an explosive radiation like the Chalcidoidea [START_REF] Heraty | A phylogenetic analysis of the megadiverse Chalcidoidea (Hymenoptera)[END_REF], another issue is the presence of monotypic or speciespoor groups that are the only extant representatives of a long line of ancestors. These lineages are characterized by long external branches and insidious LBA artefacts may occur (not necessarily with the outgroups). On the morphological side, these taxa can be highly transformed and homologies between their features and those observed in the remaining species difficult to assess. Here, the position of two genera, Zavoya and Smicromorpha, remains ambiguous. Contrary to the results from the UCEs, Zavoya (3 species known) is strongly supported as sister to Notaspidium in the morphological tree, with which it shares many characters (Figures 4-S8, Appendix S4 & S5), though only one is an unambiguous synapomorphy. Additionally, this synapomorphy is a character loss (the ventral carinae of the metatibia are absent), that may confound interpretation of homology and relationships [START_REF] Bleidorn | The role of character loss in phylogenetic reconstruction as exemplified for the Annelida[END_REF]. For all other characters shared by Zavoya and Notaspidium, the same character state is observed only in species that do not belong to the Haltichellinae. These characters could thus be considered as local synapomorphies and, together with the absence of the ventral carinae on the metatibia, could reveal an undetected LBA artefact in our UCE analysis that may be reduced with an increasing sampling of Zavoya, Notaspidium, and Haltichellinae species.

Further studies are nevertheless required to assess whether morphological convergence or systematic bias drive the position of Zavoya.

The position of Smicromorpha is also doubtful both in the morphological and the molecular trees. Even after UCE removal based on TreeShrink results (Smicromorpha is the most flagged taxon, Table S2), the long branch is still obvious (Appendix S3, Figure 4). Smicromorphinae are highly transformed parasitoids of weaver-ant larvae (Oecophylla, Formicinae) and character homologies are difficult to assess (Appendix S4-S5). During the day or at dusk, females Smicromorpha deposit their eggs on silk-spinning larva of weaverants held by workers to seal the leaves that are being pulled together by other workers when building their nest. Of the seven known species of Smicromorpha, two were reared from and three were collected flying around nests of Oecophylla smaragdina (Fabricius) in the Oriental and Australasian regions [START_REF] Naumann | A revision of the Indo-Australian Smicromorphinae (Hymenoptera: Chalcididae)[END_REF][START_REF] Darling | A new species of Smicromorpha (Hymenoptera, Chalcididae) from Vietnam, with notes on the host association of the genus[END_REF]. The biology of other species is unknown. While there is no fossil record for Smicromorpha, multiple fossils of Oecophylla are known from Europe [START_REF] Barden | Fossil ants (Hymenoptera: Formicidae) : ancient diversity and the rise of modern lineages[END_REF] suggesting that numerous species have become extinct since the Eocene Epoch (the estimate age of the first fossil is about 56 Ma). Interestingly, Oecophylla is also on a long branch and its putative relationships might therefore be artefactual [START_REF] Ward | A revised phylogenetic classification of the ant subfamily Formicinae (Hymenoptera: Formicidae), with resurrection of the genera Colobopsis and Dinomyrmex[END_REF]). The highly specialized interaction Smicromorpha-Oecophylla and the numerous extinct lineages may explain the long branches and the difficulties encountered in correctly placing these taxa in the phylogenies. Morphological data support a sister taxa relationship between Smicromorpha, Dirhininae and Epitraninae but this clade is not recovered in the UCE tree. A close relationship between Epitraninae and Smicromorpha could be hypothesized based on UCEs and morphology /behavioral data (Appendix S5, Figure 4). Further studies are required to clarify the position of Smicromorpha in the chalcidid phylogenetic tree. The addition of more species of Smicromorpha and the inclusion of an undescribed genus, the probable sister taxon of Smicromorpha in the Afrotropical region, that we unfortunately failed to sequence due to poor specimen preservation, may improve the results. However, all these species are extremely rare and difficult to collect.

Finally, a conflict is also observed between morphological and molecular data for Dirhininae and Epitraninae. These subfamilies are recovered as sister taxa in the morphological tree while they are shown to be distantly related in the molecular trees. The sister taxa relationship is mostly supported by a closely related structure of the tentorium and other characters that are homoplastic especially those of the mesoscutellum, the fore wings and the hind legs. Therefore, their close relationships in our morphological analysis may reflect convergence.

inferred from genome-scale data are usually highly supported which is often wrongly confused with accuracy, and falsely reinforces confidence in molecular results. Several authors have strongly advocated a systematic exploration of biases with different analytical methods. Indeed, this exploration may highlight better, alternative topologies that would not be revealed by a point and click approach. However, such studies are scarce, especially when data sets are composed of hundreds of taxa and genes, which makes computation time prohibitive and not compatible with the current publish or perish system. A thorough exploration of phylogenetic tree space may nevertheless reveal alternative hypotheses that are difficult to rank without independent sources of evidence. Our study highlights the power of a systematic exploration of biases to sort among conflicting phylogenomic hypotheses and the usefulness of a careful analysis of morphological features by expert taxonomist to corroborate (or not) the most likely topology. It may provide guidelines to build the tree of life of other hyperdiverse groups of animals on which little phylogenetic knowledge has been acquired, which is the rule rather than the exception in non-vertebrate taxa. Topology B is observed when the complete data set is analyzed with ASTRAL. All nodes are highly supported (localPP=1.0). Topology B is also inferred when mid-point rooting instead of outgroup rooting is used to root the trees obtained with the concatenation approaches. Topology C is observed when the complete data set is analyzed with ASTRID. However, the node grouping Chalcidini with Dirhininae-Smicromorphinae-Epitraninae-Brachymeriini-Phasgonophorini is poorly supported (BP=32). Topology C is also observed in about 20% of the bootstrap trees obtained with the concatenation approaches. Topology C is inferred with concatenation approaches when UCEs (resp. nucleotides) with GC content strictly superior to 0.48 (resp. strictly superior to 0.57) are removed from the data set, though the node grouping Chalcidini with the Dirhininae-Smicromorphinae-Epitraninae-Brachymeriini-Phasgonophorini group is only moderately supported (BP=67 (resp. 80)). Finally, Topology C is observed when the set of input gene trees for the ASTRAL analysis is reduced to trees for which ingroups are monophyletic Nota: Current higher classification is used to annotate the trees. For brevity, we will refer to the clade grouping Brachymeriini + Dirhininae + Epitraninae + Phasgonophorini + Smicromorphinae as the "BDEPS" clade throughout text. The full list of morphological character and character states is provided in Appendix S1. A detailed discussion of the character states that support the different clades can be found in

ACKOWLEDGEMENTS

FIGURE AND TABLE CAPTIONS

  trees. Regardless of the method used, Chalcididae was always recovered as monophyletic with strong support. Of the six recognized subfamilies, only Chalcidinae was not monophyletic. With the exception of Hybothoracini and Haltichellini that were recovered as polyphyletic, all tribes were monophyletic with high support values. Of the 29 non-monotypic genera in the data set only 16 were recovered monophyletic. Shallow and intermediate relationships were similar among analytical methods. About 10 unsupported topological conflicts were highlighted. Most of the time, conflicts involved taxa with a high level of missing UCEs (> 70%). Highly supported conflicts were observed in the deepest nodes of the phylogeny. While four moderately to highly supported clades were inferred by ML (Figure 1, topology A), Chalcididae clustered into two highly supported clades in the ASTRAL tree (Haltichellinae versus all other Chalcididae; Figure 1, topology B) and three poorly to highly supported clades were inferred by ASTRID (Figure 1, topology C).

  We are grateful to Cedric Mariac and Leila Zekraoui (IRD, DIADE, France) for providing access to the Bioruptor; Audrey Weber (INRA, AGAP, France) for sequencing of the libraries and the Genotoul bioinformatics platform Toulouse Midi-Pyrenees, France for providing computing resources. We thank Natalie Dale-Skey Papilloud (NHM, London) and Nicole Fisher (ANIC, Canberra) for the loan of specimens as well as the Queensland government for collecting permits (WITK18248017-WITK18278817). This work is part of a large NSF project (Award#1555808) led by John Heraty (UC Riverside USA) that attempts to solve the phylogeny of the Chalcidoidea with NGS approaches and was funded by the INRA SPE department (recurrent funding to JYR and AC). USDA is an equal opportunity employer and provider. Trade names mentioned herein are for informational purposes only and do not imply endorsement by USDA. Zhang C., Rabiee M., Sayyari E., Mirarab S. 2018. ASTRAL-III: polynomial time species tree reconstruction from partially resolved gene trees. BMC Bioinformatics, 19:153. Zhong B., Betancur-R R. 2017. Expanded taxonomic sampling coupled with gene genealogy interrogation provides unambiguous resolution for the evolutionary root of angiosperms. Genome Biology and Evolution, 9:3154-3161. Zhu Q. 2014. AfterPhylo. A Perl script for manipulating trees after phylogenetic reconstruction. Available from https://github.com/qiyunzhu/AfterPhylo/.

Figure 1 .

 1 Figure 1. Summary of the topologies recovered from the UCE data set.

Figure 2 .

 2 Figure 2. Results of the Gene Genealogy Interrogation approach (GGI) and correlation

Figure 3 .

 3 Figure 3. Results of the morphological analysis.

  Figure 2

Table 2 . GC content threshold to remove UCEs

 2

	Complete
	data set

Table 3 . GC content threshold to remove nucleotide sites

 3

	Complete
	data set

DATA ACCESSIBILITY

Demultiplexed cleaned reads are available as a NCBI Sequence Read Archive (ID#XXXX).

Custom script to detect cross contamination is available from https://github.com/DNAdiversity/UCE-Cross-Contamination-Check.

A new higher classification for the Chalcididae.

The monophyly of the family is confirmed by our study. Chalcididae is supported by multiple synapomorphies: the presence of a genal carina and a postoral bridge, a sclerotized labrum with an exposed ventral plate, the parascutal and axillar carina ∩ -shaped, a metanotal scutellar arm reduced to a thin carina, the cubital vein present as a non-pigmented fold, and a petiole fused ventrally without suture (Appendix S5). Based on our results, we propose a revised higher classification (subfamilies and tribes) for the family (Table 1, Figures 34, Appendix S5). The within tribe classification was not the purpose of this study and will be reviewed elsewhere. To be conservative, we propose to keep subfamilies in their historical taxonomic rank, but to raise the tribes Phasgonophorini and Brachymeriini to subfamily rank.

Within Haltichellinae, we recognise six major monophyletic groups that should be considered The RAxML tree inferred from the less biased UCEs (GC content ≦ 0.48) is used as a template and bootstrap support values less than 100 are reported at nodes as follows: RAxML on less biased UCEs /RAxML on less biased nucleotides (i.e. with GC content ≦ 0.57). This topology received the highest support from our morphological analysis (Appendix S5). The new classification proposed in this paper is used to annotate the tree. Color coding is similar to Figure 3. Within tribe classification including synonymization of invalid genera was not the purpose of this study and will be reviewed elsewhere. Single quotes indicate new genera awaiting description, the genus name used in the annotation is the one obtained when using current identification keys.

Table 1. Current and new higher classification of the Chalcididae.

Table 2. Properties of the data subsets analyzed to detect possible systematic bias due to compositional heterogeneity among loci (GC content) and impact on the branching patterns of the deepest nodes. Table 3. Properties of the data subsets analyzed to detect possible systematic bias due to compositional heterogeneity among nucleotide sites (GC content) and impact on the branching patterns of the deepest nodes.

AUTHOR CONTRIBUTIONS

Designed the study: JYR, AC, GD; obtained funding: JYR, AC; contributed samples: GD, JYR, SF, SvN; identified samples: GD, JYR; performed lab work: SN, LS, SF; contribute sequences (UNSM samples): BBB, MG, SB; wrote script to detect cross-contaminations: SR; assembled the morphological matrix: GD; analyzed morphological data: JYR, GD, AC; analyzed molecular data: AC, JYR; provided help with R: MC, JPR; drafted the manuscript: AC, JYR, GD. All authors commented on the manuscript.