We examined more than 2,000 parasitoid wasps associated with psyllids sourced from both historical museum collections and extensive fieldwork that covered 22 provincial-level administrative divisions of China. To rear parasitoids, plant materials (twigs, leaves, and bark) infested with late-stage nymphs of psyllids were collected and maintained in glass vials or fine-mesh nylon bags. These rearing efforts yielded approximately 1,000 encyrtid wasps from 400 psyllid-rearing batches. Of all specimens examined, 263 representative individuals were selected for sequencing (Suppl. material 1). When multiple psyllid species infested the same host plant, each species was isolated and reared separately for at least two months to ensure complete parasitoid emergence. The identified host psyllids included gall-forming and non-galling, foliage-dwelling taxa representing four families–Aphalaridae, Carsidaridae, Psyllidae, and Triozidae; detailed host identities and feeding habits are provided in Suppl. material 1. Upon emergence, parasitoids were preserved in 95% ethanol for molecular analysis.

For morphological identification of species, we used the following keys: Trjapitzin (1989), Noyes and Hanson (1996), Hayat (2006), Zhang and Huang (2004) and Zou et al. (2023). All specimens used for molecular analyses and listed in Suppl. material 1 are deposited in the Institute of Zoology, Chinese Academy of Sciences.

Approximately 67% of the encyrtid parasitoids examined were used in DNA barcoding study. Total genomic DNA was extracted from each individual using the DNeasy Blood & Tissue Kit (Qiagen, Germany) following the manufacturer’s protocols. The COI gene was amplified and sequenced using the universal DNA barcoding primers LCO1490 and HCO2198 (Folmer et al. 1994). For specimens that proved difficult to amplify, the primer pair FWPTF (Li et al. 2010) and LepR1 (Hebert et al. 2004) was used. Sequences of the 28S-D2 region (D2 expansion segment of 28S rDNA) were generated using the primer pairs D2-3549 and D2-4068 (Campbell et al. 1993), or D2-3566 and D2-4057 (Gillespie et al. 2005). Reactions (30 μl) consisted of 3.3 μl of template DNA (not quantified), 9 μl of reaction buffer, 4 μl of dNTP mixture (2.5 mM each), 0.5 μM of each primer and 1.25 units of ExTaq polymerase. Following Deng et al. (2012) and Zhang et al. (2011), amplification comprised an initial denaturation at 95 °C for 3 min; 35 cycles of denaturation at 94 °C for 30 s, annealing at 52 °C for 50 s, and extension at 72 °C for 1 min; and a final extension at 72 °C for 10 min. All of the sequences were queried in the National Center for Biotechnology Information (NCBI) nucleotide database (http://blast.ncbi.nlm.nih.gov/Blast.cgi) with default parameters.

DNA sequences were checked and assembled with BioEdit (Hall 1999). In MEGA X (Kumar et al. 2018), consensus sequences were translated with the mitochondrial invertebrate genetic code to check for stop codons, and genetic interspecific and intraspecific distances were calculated with the Kimura two-parameter (K2P) model (Kimura 1980). Both distance-based (ASAP and ABGD) and tree-based (GMYC and bPTP) methods were used to evaluate barcoding success rates in the COI data set.

To accurately assess host utilization by parasitoids, the food web was reconstructed using species delimitation results combining the morphological and DNA barcoding analyses. Network-level metrics, including connectance and the specialization index H2' (Blüthgen et al. 2006), were calculated using the bipartite package (v2.18) (Dormann et al. 2009) in R (v4.4.2). Species-level specialization indices were also calculated: d’ (degree of host monopolization; Blüthgen et al. 2006) and Paired Difference Index (PDI; reliance on primary host; Poisot et al. 2010), and Species Specificity Index (SSI; evenness of interaction frequencies across hosts; Julliard et al. 2006). All analyses were performed using the bipartite (v2.18) (Dormann et al. 2009) and ESM (Poisot et al. 2012) packages in R (v4.4.2).

Morphological examination placed these specimens within 4 genera: Diaphorencyrtus Hayat, 1981, Psyllaephagus Ashmead, 1900 (the most abundant genus, more than 85% of specimens), TrechnitesThomson 1876, and Prionomitus Mayr, 1876. Preliminary sorting resulted in the recognition of 29 distinct morphospecies, of which Psyllaephagus contributed 24 species, while the remaining four were distributed among Prionomitus (2 species), Trechnites (2 species), and Diaphorencyrtus (1 species).

A total of 212 COI sequences representing 28 morphospecies in 3 genera were used for barcode analysis, with 209 newly generated from collected individual psyllids and three downloaded from NCBI (GenBank: MZ438311, MZ436006 and OP787025). We were not able to obtain COI sequences from Psyllaephagus calligonicola or Diaphorencyrtus aligarhensis. The length of all sequenced fragments was 540 bp after edge trimming, with 288 conserved sites, 246 variable sites, and 250 parsimony-informative sites. No insertions, deletions, or stop codons were noted in the sequences. All COI sequences had a bias toward low GC content (A = 26.3%, T = 45.6%, C = 13.5%, and G = 14.7%).

The mean intraspecific divergence was 0.012 (ranging from 0.00 to 0.0879), the mean interspecific divergence was 0.1492 (ranging from 0.0537 to 0.2276). The interspecies divergence within genera varied from 4.88% to 13.58%: Trechnites (mean = 7.55%), Prionomitus (mean = 4.88%) and Psyllaephagus (mean = 13.58%) (detailed values of each genus are provided in Suppl. material 2). For intraspecific differentiation, 1) the lowest divergence occurred between Prionomitus tiliaris, Trechnites sp 1, Psyllaephagus femoralis, Psyllaephagus sp 2, Psyllaephagus sp 4 (0.00) and Psyllaephagus obliquus (0.0004); 2) the highest diversity was observed in Psyllaephagus caillardiae (0.0879), followed by Psyllaephagus arenarius (0.0471), suggesting higher population-level variability; 3) most groups (e.g., Psyllaephagus nartshukae, Psyllaephagus colposceniae) exhibited moderate divergence (0.0012–0.0318). For interspecific differentiation, 1) the lowest divergence occurred between Psyllaephagus sp 2 and Psyllaephagus sp 3 (0.0537), indicating potential close relatedness; 2) the highest divergence was observed between Psyllaephagus sp 2 and Trechnites sp 1 (0.2276), supporting their distinct species status (see Fig. 1 and Suppl. material 2). The length of the 255 nuclear 28S sequences ranged from 519 bp in Psyllaephagus stenopsyllae to 481 bp in Prionomitus tiliaris. The mean inter- and intraspecific divergences were 4.44% and 0.406%, respectively.

In the ABGD analysis, K2P and JC69 models yielded identical results: initial partitioning suggested 32 groups at P = 0.02976, while recursive partitioning increased group numbers to 33–45 at intermediate P values, stabilizing at 71 groups for P ≤ 0.002069. The simple distance model produced fewer groups initially (32–40) but stabilized at 17 groups under recursive partitioning for P ≤ 0.00336. At a priori genetic distance thresholds of 0.02967 and 0.04833, ABGD strongly supported the presence of 32 genetically distinct partitions in the COI data set (see Suppl. material 4: table S4-1). The number of partitions produced by ABGD was in accord with morphological identifications. ABGD strongly supported the presence of 32 genetically distinct partitions in the COI data set. Among these morphological species, P. colposceniae, P. chinensis, P. caillardiae, and P. arenarius were each distinguished into two groups (Fig. 2). The ASAP analysis produced 30 MOTUs under JC69, the simple models (ASAP scores 6.00, 6.50) and K2P model (ASAP score 7.00) (Suppl. material 4: table S4-2). In contrast to the ABGD analysis, the ASAP analysis with the K2P model revealed more MOTUs within three morphospecies, splitting P. colposceniae into three and P. longiventris into two. Additionally, the JC69 and simple models grouped P. desertus and P. caillardiae (from Changji, Xinjiang) as a single MOTU (Fig. 2).

A total of 113 COI haplotypes and 81 nuclear 28S haplotypes were used to construct the NJ trees. Based on tree topology and branch support, the COI haplotypes were resolved into 32 distinct clades corresponding to the 32 consensus MOTUs, whereas the 28S haplotypes were resolved into 34 clades (Suppl. material 5, 6).

The GMYC single-threshold analysis identified 38 species (95% CI: 32–43) and 39 entities (95% CI: 33–44). The likelihood ratio test was highly significant (LR = 77.67, P < 0.001), rejecting the null coalescent model. Threshold time was -0.0092 relative to tree root. Consistent with distance-based methods, GMYC supported delimitation of P. chinensis, P. arenarius, P. colposceniae, and P. caillardiae as distinct MOTUs, while further splitting P. desertus, P. densiciliatus, and P. stenopsyllae (Fig. 2). The bPTP analysis showed an acceptance rate of 34.7%, with 249,749 merge events and 250,251 split events. The estimated number of species ranged broadly between 42 and 68, with a mean of 53.6.

Maximum Likelihood analyses of our dataset distinguished 49 MOTUs; and Bayesian analyses distinguished 47 MOTU (Fig. 2). Although the tree-based bPTP method yielded a significantly higher count of MOTUs (47–49), likely due to over-splitting of intraspecific population structure, the distance-based ABGD and tree-based GMYC analyses showed strong congruence in revealing cryptic diversity. Both methods consistently supported the splitting of four morphospecies (P. colposceniae, P. chinensis, P. caillardiae, and P. arenarius) into distinct lineages. Consequently, adopting a conservative consensus approach across multiple delimitation methods, we recognized a total of 32 MOTUs from the 29 initial morphospecies in this study (Fig. 2).

To assess host use pattern, a tanglegram showing host–parasitoid interactions was established for those species with host records (Fig. 3). The network comprised 21 host species and 23 parasitoid species. Nine of the 32 MOTUs were excluded from the network because host records were unavailable. The resulting bipartite network was sparse but highly structured, with 19 parasitoid species utilizing only a single host species; two parasitoid species utilizing two host species, and two parasitoid species utilizing three host species. At the network level, connectance was low (0.06) and the specialization index was exceptionally high (H2’ = 0.9205), indicating strong overall specialization. Species-level metrics reinforced this finding: over half of the parasitoid species showed extreme host specificity (d’ = 1.0). The few apparent generalists—P. stenopsyllae and P. guangxiensis—each attacked three host species and had relatively low SSI values (<0.6). However, their high PDI values (>0.9) revealed that even these species strongly preferred a single primary host. Additionally, we found seven psyllid species were parasitized by two or three parasitoid species.

DNA barcoding analyses demonstrated that COI sequences effectively discriminate all encyrtid morphospecies. The presence of a clear barcoding gap and high bootstrap support for species-level clades in phylogenetic trees (Suppl. material 5) indicate that COI possesses sufficient variation to distinguish most species in this group. These findings are consistent with successful COI barcoding applications in other encyrtid parasitoid wasps (Zhang et al. 2011; Chesters et al. 2012; Wang et al. 2016; Qin et al. 2018) and parasitoid wasp families more broadly (e.g., Smith et al. 2008, 2013; Xu et al. 2023). Notably, even between closely related species (e.g., Psyllaephagus sp 2 and Psyllaephagus sp 3), COI exhibited an average interspecific genetic distance of 5.37%, comparable to values reported for the Encyrtus sasakii complex (Chesters et al. 2012).

In the present study, the minimum interspecific threshold observed for COI (approximately 5.37%) may be attributable to its slower evolutionary rate compared to other mitochondrial markers. This interpretation is supported by recent mitogenomic analyses of Encyrtidae, which revealed that COI exhibits the lowest substitution rate among all 13 mitochondrial protein-coding genes (Zhang et al. 2024). Using mitochondrial Cytochrome b (Cytb), Hall et al. (2017) found the three Psyllaephagus morphospecies reared from eight Eucalyptus-feeding psyllid species might contain nine putative cryptic species, suggesting that Cytb, with its relatively higher interspecific genetic distance threshold (approximately 8.4%), may resolve additional cryptic species. The available evidence suggests that COI may evolve relatively slowly in encyrtid parasitoids associated with psyllids, although comparisons with a broader range of mitochondrial genes are needed to confirm this inference.

Although the 28S-D2 region is known to be more evolutionarily conserved than COI (Heraty et al. 2004), it remains a valuable complementary tool for species confirmation in Encyrtidae (Zhang et al. 2011; Wang et al. 2016; Qin et al. 2018). Our results support the utility of this marker for distinguishing encyrtid parasitoid species associated with psyllids, as it exhibited a clear barcode gap in interspecific genetic distances (Suppl. material 3: table S3-1, S3-2). Notably, 28S-D2 distinguished between closely related species (e.g., P. chinensis and P. obliquus) based on fixed diagnostic nucleotide substitutions, demonstrating that even a conserved nuclear marker can provide reliable species-level resolution when diagnostic sites are present. The presence of such diagnostic characters enhances the reliability of molecular identification and provides a stable reference for future taxonomic and biological control studies.

The advent of molecular taxonomy has greatly facilitated the discovery of cryptic diversity (Bickford et al. 2007). Numerous studies have demonstrated that cryptic species are prevalent in parasitoid wasps (Heraty et al. 2007; Smith et al. 2008; Nicholls et al. 2010), with some morphospecies in Encyrtidae comprising multiple MOTUs (Wang et al. 2016). However, the extent of cryptic diversity and the contribution of molecular approaches vary among taxonomic groups. Evidence from large-scale taxonomic surveys increasingly supports high concordance between morphology-based and molecular-based species delimitation in lineages with unambiguous diagnostic characters (Jablonski and Finarelli 2009; Chen et al. 2011; Pentinsaari et al. 2014; Toro-Delgado et al. 2022). In such morphologically challenging groups associated with psyllids, we found some characters, such as the shape of the oblique portion of clava and the coloration of coxae, provide complementary evidence to molecular data, helping to resolve species boundaries and improve taxonomic resolution.

Among Chinese encyrtid parasitoids associated with psyllids, we observed a high degree of concordance between morphological and molecular boundaries, with cryptic lineages detected in only several morphospecies (e.g. P. chinensis and P. arenarius). The congruence between molecular delimitations and morphological classifications demonstrates that morphological characters provide reliable and practical means for species identification in this group, reducing the likelihood that morphologically cryptic lineages have been overlooked.

The degree of host specificity of encyrtid parasitoids of psyllids is not known with any certainty, but it appears that host specificity is presumed to be high (Noyes and Hanson 1996); however, this critical biological trait has been rigorously tested in only a few cases. By integrating molecular delimitation with extensive geographic sampling and network analysis, we provide a framework for host specificity of these parasitoids. Our result revealed that 19 species with known host records are restricted to a single host species, and two species are associated with three host species. Our findings indicated a high level host specificity exists in encyrtid parasitoids associated with psyllids. Even a few species expand beyond a single host species, but in these cases the hosts are usually closely related species. This result is consistent with other studies of encyrtid parasitoids associated with scale insects (Zhang et al. 2011; Chesters et al. 2012; Wang et al. 2016; Qin et al. 2018). As endoparasitoids, they must overcome specific host defenses (Salt 1968), whereby many physical and chemical factors may become increasingly refined in the interactions with the hosts (Vinson 1976).

The introduction of molecular data has refined the understanding of parasitoid host specialization (Smith et al. 2006, 2008, 2013; Kaartinen et al. 2010). Patterns of host use in Psyllaephagus often reveal high specificity, even within species that exploit multiple host species (Fromont et al. 2017; Hall et al. 2017). Our findings for these parasitoids align with this trend: while monophagy is the predominant pattern (H2'=0.9205), several species, such as P. stenopsyllae and P. guangxiensis, were found to utilize multiple hosts. However, the high Paired Difference Index (PDI ≥ 0.9) for these species indicates a strong reliance on a single primary host. This outcome is consistent with patterns observed in encyrtid parasitoids of scale insects (Zhang et al. 2011; Qin et al. 2018).

Collectively, these insights suggest that host-parasitoid associations may be tighter than morphological assessments imply. This has significant implications for biological control: while high specialization minimizes non-target risks, it necessitates precise host-parasitoid matching to avoid control failure in dynamic environments.