Phylogenetic analyses

Gene alignments were completed for COI and Cytb by hand using the reading frame as a guide with Bioedit v.7.1.3 [54]. There were no indels present in either gene and thus the alignments were unambiguous. A modified [52] secondary structure model [51] was used to align 28S. Regions of ambiguous alignment (RAAs), and regions of expansion and contraction (RECs) were not excluded from the data set as we assumed most informative characters for this gene would be contained within these regions that can be hyper variable across genera. For each gene, the best fitting model of DNA sequence evolution for nucleotide analyses were determined using jModelTest v.0.1.1 [55] under the Bayesian Information Criterion (BIC). The model with the lowest calculated BIC score was considered the best-fitting model for each gene. Depending on the gene, either one or two species from other braconid subfamilies were used as outgroups Afrocampsis sp. (Acampsohelconinae) and Eumacrocentrus americanus (Cresson 1873) (Helconinae) to ensure the ingroup was monophyletic.

For individual and concatenated data sets, Bayesian inference with two independent runs each with four chains and default priors was run in MrBayes v.3.2.6 [56]. An independent molecular model was applied to each partition in the concatenated data set (partitioned by gene) and different parameters of the model were unlinked to allow each partition to have its own set of estimations for parameters. The rate parameter was set to vary across different partitions to incorporate rate heterogeneity across partitions. Data sets for all gene alignments were deposited in Figshare (https://figshare.com/): 10.6084/m9.figshare.6149219.

All analyses were performed for 5,000,000 generations sampling every 1000^th generation and the results from the two independent runs were summarized in a majority rule consensus tree after discarding the initial 25% of the trees for burn-in. Stationarity and appropriate mixing of the two independent runs were determined when the average standard deviation of split frequencies approached 0.01, the Potential Scale Reduction Factor (PSRF) for each parameter of the model was close to 1, and the overlay plots of both runs showed the number of generations versus the log probability of the data were heteroscedastic. Average intraclade and interclade genetic distances were calculated using Kimura’s two-parameter model [57] using MEGA v.7.0.14 [58].