Classical biological control requires parasitoids to efficiently locate their hosts and exhibit a host range that is narrow enough to eliminate or reduce the chances of unwanted non-target effects. Phylogenies can reveal the mechanisms that contribute to these traits by determining if they are phylogenetically constrained or are highly variable within the genus. The primary candidate as a biological control agent for H. halys is Trissolcus japonicus (Ashmead), a species for which adventive populations are now in USA, Canada, Switzerland, and Italy (Talamas et al. 2015a, Abram et al. 2019, Stahl et al. 2018, Sabbatini-Peverieri et al. 2018). However, it is not the only species of Trissolcus that can parasitize H. halys, and there are many species that attempt to parasitize of H. halys eggs with limited or no success. Our analysis examines how these traits of host acceptance and host competence are distributed within Trissolcus. Taxon sampling, while focused on species associated with H. halys, includes additional representatives for the basalis, flavipes, and thyantae species groups from the Holarctic. We consider this phylogeny to provide a backbone for future molecular studies on Trissolcus wasps that will undoubtedly occur as interest in the group continues to grow and specimens from a broader geographic sampling become available.

The present study is not the first phylogenetic effort for Trissolcus. Johnson (1987) provided a phylogenetic hypothesis for Trissolcus species which was cladistic in its argumentation, but not computationally optimized (i.e. characters were optimized on a tree, but tree space was not searched). The result was a partially resolved phylogeny and only concerned Nearctic species. Later, Johnson (1991) improved the resolution of Trissolcus by employing a character matrix and analysis in PAUP. However, computer limitations prevented a complete data matrix, as a great deal of homoplasy occupied the tree memory storage of PAUP at that time. As a result, Johnson (1991) reduced the size of the matrix and relied on ground-plan coding for some taxa; the resulting data matrix recovered a single tree.

The first molecular sequence data for Trissolcus were provided by Murphy et al. (2007), who investigated higher level relationships in Platygastroidea using three gene fragments. The results of that study confirmed the placement of Trissolcus in Telenominae with very high bootstrap support. Guz et al. (2013) were the first to investigate the relationships within Trissolcus using molecular data. Here the focus was on Trissolcus species that were natural enemies of the sunn pest (Eurygaster integriceps Puton (Hemiptera: Scutelleridae) of wheat and barley. While the study has limited utility with respect to relationships within the genus, Guz et al. (2013) demonstrated the usefulness of the COI marker for species-level questions, and reported that due to insertions, ITS2 was difficult to align, and that 28S, 18S, and 5.8S were too conserved to be informative.

Taekul et al. (2014) included 12 species of Trissolcus in an analysis that redefined the limits of Telenominae. Their phylogeny was based on four molecular markers (18S, 28S, CO1 and Ef-1α) and focused on shifts in host selection in Telenominae, Gryon Haliday, and the Psix-cluster of genera. Importantly, it demonstrated the utility of 18S and 28S sequence data for phylogenetic analysis of Trissolcus and its relatives.

Species determinations were made with the identification tools provided in Talamas et al. (2015a), Talamas et al. (2017), and Tortorici et al. (2019). The data associated with these specimens, including host associations, are deposited in Hymenoptera Online (hol.osu.edu) and can be accessed via the Collectin Unit Identifiers listed in Suppl. material 1. Voucher specimens from this study are deposited in the National Museum of Natural History (Washington, DC) and the Florida State Collection of Arthropods (Gainesville, FL).

Most specimens were collected alive and fixed in 95% or absolute ethanol and some were gleaned from material stored in ethanol in entomological collections. These specimens were used for nondestructive DNA extraction using the Qiagen DNeasy kit (Hilden, Germany) following the protocol published in Taekul et al. (2014), but with minor modifications specified in Sabbatini Peverieri et al. (2018). Individual specimens were bathed three times at room temperature in molecular grade water for five minutes prior to overnight incubation in lysis buffer at 55 °C. In step 7 of the Qiagen protocol, the elution buffer was warmed to 55 °C and allowed to rest on the membrane for 15 minutes before centrifugation. The collected flow-through was reloaded onto the spin column to increase the DNA yield. When we started this study, the nondestructive method was not employed, and therefore, some specimens were entirely ground using the Qiagen DNeasy kit (Hilden, Germany) following the manufacturer’s recommendations. These specimens thus have no corresponding voucher specimen. A negative control (no insect tissue) was included in each extraction to detect potential contamination. The genomic DNA was stored at -24 °C for further use. All voucher specimens are deposited in the Florida State Collection of Arthropods (Gainesville, FL), and the National Insect Collection, National Museum of Natural History (Washington DC, USA).

Five molecular markers were sequenced. These included the mitochondrial 5’ end of the cytochrome c oxydase subunit I gene (CO1) also named the barcode region (~660bp), the nuclear ribosomal gene 18S rRNA (variable region V3-V5, ~780bp) , the 28S rRNA (D2-D3 expansion regions, ~800bp), the internal transcribed spacer 2 (ITS2), (~550bp to 650bp), and the nuclear gene Wingless (exon, ~450bp). The choice of these markers was partly guided as a compromise between ‘‘top down” and “bottom up’’ approaches (Wiens et al. 2005). We apply the bottom up approach to resolve higher level relationships (the bottom of the tree) using relatively slowly evolving markers (18S rRNA, 28S rRNA, Wingless) and then apply the top down approach to resolve species level relationships (the top of the tree) using faster evolving markers (CO1, ITS2). Primers and PCR conditions used in this study are described in Tables 1, 2, respectively. All PCRs were performed in a 30 µl total volume with 2 µl of DNA template, 0.2 mM of each dNTP, 0.3 µM of each primer, 1× CoralLoad PCR Buffer (including 1.5mM of MgCl₂) and 1 Unit of Taq DNA Polymerase (Qiagen). PCR amplifications were run on a 9700 thermocycler (Applied Biosystem). The PCR products were purified and sequenced in both directions using the same sets of PCR primers, by Genoscreen, Lille, France, whereas others were cloned prior to sequencing (especially ITS2). Both strands for each overlapping fragment were assembled using the sequence editing software Bioedit, version 7 (Hall 1999). All sequences have been deposited in GenBank and accession numbers are provided in Suppl. material 1. All residual DNAs are archived (-24°C) at the European Biological Control Laboratory (EBCL, USDA/ARS), Montpellier, France.

The protein coding genes CO1 and Wingless were aligned using ClustalW with default gap opening, extension, and substitution costs as implemented in Mega 6 (Tamura et al. 2013). These sequences were checked for stop codons and frame shifts, and sequences were translated to amino acids using the invertebrate mitochondrial code and the standard code respectively as implemented in MEGA 6 (Tamura et al. 2013). Secondary structural alignments were implemented for ribosomal RNA sequences of 18S, 28S and ITS2. The ClustalW alignment conventions followed Kjer (1995) with slight modifications (Gillespie 2004). Ambiguous regions in ITS2 were excluded from the final analyses using GBlock as implemented in PhyML 3.0 (Guindon et al. 2010). The aligned, partitioned sequence data is provided in Suppl. material 2.

Bayesian inference. The resulting concatenated matrix was exported from Mesquite for Mr. Bayes 3.2 applying the GTR+I+G rate matrix for each data partition (COI divided into three partitions, one for each position) and running 15 million generations with a burn-in of 25%; explanation and justification of these protocols are in Buffington et al. (2007).

Parsimony. The parsimony searches were conducted using PAUP* (Swofford 2002), employing an initial 10000 replicate searches of TBR under equal weights with branches of maximum length zero collapsed and steepest descent set to ‘off’. For bootstrap analyses (Felsenstein 1985), a simple addition sequence was employed with Telenomus (Te. californicus complex sensu Johnson (1984b)) set as the reference taxon, followed by 1000 bootstrap replicates, with each employing 100 TBR swapping replications. As many equally parsimonious trees were found in the initial tree search, successive approximations (Farris 1969) were used to converge on a topology favored by the characters with the best tree score. A separate analysis was run in TNT (Goloboff et al. 2008) employing sectorial searches, parsimony ratchet, and tree fusing.

Maximum likelihood. These analyses were run using RAxML version 8.2.10. The model used was GTRGAMMA+I. Automatic bootstopping criterion was selected as the appropriate number of bootstraps; 300 replicates were run. Six partitions were identified using PartitionFinder 2. The proportion of gaps/undetermined sites in the alignment was 11.47%. All resulting trees were visualized in FigTree 1.3.1, and the out-group (Telenomus) was assigned; the final tree figure was generated using Adobe Illustrator. The commands used to perform each analysis are listed in Suppl. material 3.

flavipes group

The flavipes group sensu Talamas et al. (2017) was retrieved as a monophyletic clade in the parsimony analysis, but with T. mitsukurii, T. latisulcus, and T. thyantae included (Figure 3). The Bayesian and ML analyses both retrieved the flavipes group as two separate clades in a polytomy with the basalis group (Figs 1–2), with T. mitsukurii sister to a flavipes clade comprised of primarily Asian species. Talamas et al. (2017) treated T. mitsukurii and T. latisulcus as part of the basalis group based on the number of clypeal setae (6), absence of a hyperoccipital carina on the medial vertex, and glabrous metapleuron. However, each of the analyses indicate that T. mitsukurii is better accommodated in the flavipes group. A morphological character supports this new hypothesis: in T. mitsukurii the orbital furrow is expanded at its intersection with the malar sulcus, which is found only in the flavipes group, at least among the species in this phylogeny.

Figure 1. 

Phylogenetic tree, RaxML analysis.

Figure 2. 

Phylogenetic tree, Bayesian analysis.

Figure 3. 

Phylogenetic tree, PAUP* analysis.

We thus retain much of the previous concept of the flavipes group, but the inclusion of T. mitsukurii means that the number of clypeal setae can be 2, 4, or 6. The number of clypeal setae remains a useful character because having 4 or fewer clypeal setae is limited to this group. We therefore redefine the flavipes group based on the following characters: clypeus with 2–6 clypeal setae; hyperoccipital carina usually complete, sometimes weakened or absent between lateral ocelli; orbital furrow often expanded at intersection with malar sulcus; metapleuron glabrous. This approach ignores the ambiguity of the polytomy in the Bayesian and ML analyses, and the presence of T. latisulcus and T. thyantae retrieved within the flavipes group by the parsimony analysis. Because the results do not fully agree, we prefer an approach that minimizes changes to the infrageneric organization until consensus and better supported resolution is achieved through increased sampling of species and molecular markers. To further examine the degree of homoplasy in morphological characters and their utility for delimiting species groups, we recommend that future efforts include species that do not fit into the current species groups (e.g. T. atys (Nixon), T. tersus Lê, T. levicaudus Talamas) and species from Asia and Africa that are morphologically similar to T. mitsukurii, of which there are many.

thyantae group

The thyantae group was represented by a single species, T. thyantae. The Bayesian and RaxML analyses retrieved it as the most basal lineage of Trissolcus (Figs 1–2), whereas the parsimony analysis placed it within the flavipes group as sister to T. latisulcus (Figure 3). However, Trissolcus thyantae and T. latisulcus are not morphologically similar to each other or to other species in the flavipes group. Increased taxon sampling is needed to address the ambiguity in the placement of the thyantae group and assess relationships between the morphologically similar species that constitute it.

basalis group

The basalis group remains largely unchanged regarding its constituent species and the characters that delimit it: clypeus with 6 or more setae; hyperoccipital carina absent between lateral ocelli; metapleuron glabrous; orbital furrow not expanded near intersection with malar sulcus. In both the Bayesian and RaxML analyses, Trissolcus latisulcus and T. erugatus were retrieved as a paraphyletic group sister to the other members of the basalis group (Figures 1–2). Excluding the aberrant placement of T. latisulcus in the parsimony analysis, the basalis group was consistently retrieved as monophyletic, but with varying topologies among its species.

Numerous species were treated by Talamas et al. (2017) as junior synonyms of T. semistriatus (Nees von Esenbeck). Tortorici et al. (2019) reexamined characters previously treated as variable within T. semistriatus and further updated the classification of Palearctic Trissolcus, resurrecting T. colemani (Crawford) and T. manteroi (Kieffer) as valid species and establishing name usage for T. belenus (Walker). Although T. belenus was described in 1836, this species name was largely ignored in literature on Palearctic Trissolcus because it had not been reliably characterized. Tortorici et al. (2019) examined the lectotype of this species, established a means of separating it from other members of the T. semistriatus complex and provided records of it parasitizing H. halys eggs in Europe. Although our analysis did not include T. manteroi, it confirms the conclusion of Tortorici et al. (2019) that T. belenus, T. colemani, and T. semistriatus are distinct species.

The ability to develop in H. halys eggs is not constrained phylogenetically, but the species with high rates of successful parasitism are all found in the flavipes group (T. mitsukurii now included). The closest relative of T. japonicus in our analysis, T. plautiae (Watanabe), has been reared from H. halys eggs in Asia, but accounted for only 2% of parasitism in a study by Zhang et al. (2017). Trissolcus cultratus (Mayr) and T. mitsukurii have appreciable rates of parasitism on H. halys eggs (Zhang et al. 2017; Sabbatini-Peverieri et al. 2018), leading to host range testing for these species. Trissolcus euschisti (Ashmead) and T. edessae Fouts (flavipes group) have been reared from H. halys eggs in North America, but at very low rates if the eggs are viable, indicating that they recognize H. halys as a potential host but are largely unable to complete development (Abram et al. 2017). Outside of the flavipes group, T. basalis (Wollaston) and T. solocis Johnson (basalis group) have been reared from live, sentinel H. halys eggs (Balusu et al. 2019a, Balusu et al. 2019b), but these records are considered to be rare events.

A phenomenon that deserves future attention is the geographic division in the ability of Trissolcus cultratus to successfully develop in live H. halys eggs. Our analysis retrieved a European specimen of T. cultratus (TFLA4) nested well within a clade of Asian specimens, supporting the conclusion that this is a single widespread species. However, European populations of T. cultratus fully develop and emerge from H. halys eggs only if they were previously frozen or had defenses compromised by parasitism from another species (Haye et al. 2015, Konopka et al. 2016). Given the rate at which adventive populations of Asian parasitoids follow the movement of H. halys, it is likely that an Asian population of T. cultratus will eventually appear in Europe. If this occurs, identification of the exotic population will require molecular diagnostics because they cannot be differentiated based on external morphology. In any biocontrol project such as this one, molecular characterizations of different populations should be done as soon as possible before population mixing has a chance to occur.

In recent years, DNA barcode sequences have increasingly been used to confirm morphology-based identification of Trissolcus species (Ganjisaffar et al. 2018, Balusu et al. 2019b, Talamas et al. 2015b, Abram et al. 2019, Stahl et al. 2018, Sabbatini Peverieri et al. 2018). In some cases, this is primarily a supplement to morphological diagnosis, and in others it is an invaluable means of confirmation. For example, in Ganjisaffar et al. (2018) the initial detection of Trissolcus hyalinipennis Rajmohana & Narendran in California was based on single male specimen that lacked some of the diagnostic female characters, and specimens of Trissolcus basalis had reduced morphology because of the diminutive size the bagrada bug eggs from which they emerged. For both species, the use of DNA barcoding provided an additional level of confidence in their identification. More recently, Gariepy et al. (2019) published a method that enables identification of parasitoids via residual DNA in empty, parasitized eggs. This method has since been used to detect the first population of Trissolcus japonicus in eastern Canada (Gariepy and Talamas 2019) and confirm that a species of Idris Förster (Scelionidae) parasitized bagrada bug eggs, the first non-spider host for the genus (Lomeli-Flores et al. 2019). Each of these examples relied on a pre-existing library of CO1 sequences that were reliably matched to species names. In this study, we provide CO1 sequences for 20 species of Trissolcus. The names associated with these sequences are provided with the highest level of confidence possible, given that the specimens were identified in the context of the most recent and thorough taxonomic treatments and with direct comparison to primary types.