Genetic Diversity of Tamarixia radiata Populations and Their Associated Endosymbiont Wolbachia Species from China

: Tamarixia radiata is one of the established biocontrol pests against the major Asian citrus psyllid, Diaphorina citri , a vector of Candidatus Liberibacter that is a causal agent of citrus Huanglongbing (HLB) disease. Updated information and regional exploration on biocontrol pests are important elements for effective disease management strategies. In this study, the diversity and parasitism rate of T. radiata populations were evaluated. Due to the importance of the host–parasitoid relationship, the presence of Wolbachia as an endosymbiont was also investigated. The parasitism rate of various T. radiata populations from Ecuador and China ranged between 57.27% and 66.32%, respectively, with a non-signiﬁcant emergence rate and a statistically similar sex ratio. Sequence analysis of ITS and COI from T. radiata populations was consistent with the morphological hypothesis that the collections represent a single species, whereas phylogeny of the wsp gene conﬁrmed the presence of Wolbachia pipientis as an endosymbiont within T. radiata populations. Based on partial COI sequences, the maximum genetic diversity such as total haplotype diversity (Hd = 0.788), nucleotide, diversity ( π = 0.2439), and average nucleotide difference (k = 171.844) was also estimated for different T. radiata populations. Furthermore, neutrality tests based on COI sequences indicated an overall contraction in T. radiata populations, whereas an expansion trend was observed in associated W. pipientis strains. This study clearly demonstrated the presence of genetically diverse T. radiata populations that were able to parasitize D. citri effectively, and these can be further explored as promising biocontrol candidates in integrated pest management strategies to solve citriculture economic loss caused by D. citri .


Introduction
Parasitoids constitute a species-rich group containing more than 20% of the world's insect species [1], divided into six families (Aphelinidae, Encyrtidae, Eulophidae, Eupelmidae, Pteromalidae, and Signiphoridae) within the order Hymenoptera. According to the parasitoids' taxonomic position, Eulophidae contains many species associated with insect A total of twelve T. radiata, populations, four each from Fujian (Fuzhou), Guangdong (Zhaoqing), and Jiangxi (Ganzhou) Figure S1, were collected from May to July in 2020. We also used one population from Guayaquil city, Guayas, Ecuador. In the previous study, the life table of Ecuador population of T. radiata has already been reported without any comparative molecular data [44]. Therefore, we decided to include Ecuador population of T. radiata (which is already maintained in our lab) for a preliminary comparative molecular analysis with other T. radiata populations collected from China. The different populations of T. radiata were maintained on D. citri nymphs and adults collected from citrus orchards in Fujian in May 2020. The stock populations were reared on orange jasmine (Murraya paniculata) (Sapindales: Rutaceae) and kept in mesh cages (0.60 m wide; 0.50 m deep; 0.50 m high) under control conditions (25 ± 2 • C; 70 ± 5% RH; a photoperiod of 14 h light and 10 h dark) in the laboratory of Insect Ecology and Biological Control, Fujian Agricultural and Forestry University, Fuzhou, China [44]. This study was conducted in three phases: biological parameters, molecular analysis, and endosymbiont characterization.

Biological Aspects (Parasitism, Emergence, and Sex Ratio) of Different Populations of Tamarixia radiata
Each population of T. radiata was reared and maintained in 5 plastic rearing bottles (0.2 m in height and 0.08 m in diameter). Each bottle contained 20 nymphs of D. citri (thirdfifth instar) placed on a 15 cm tall M. paniculata shoot. To keep the shoot fresh for several days, it was implanted in a water-soaked cotton plug in a polypropylene deli container (25 mL) ( Figure S2). Previously matted females T. radiata (24 h old) were introduced into each bottle containing nymphs for 48 h and then removed. Subsequently, potentially parasitized nymphs in rearing bottles were kept at 27 ± 1 • C, RH 70 ± 5% with L: D = 14:10 photoperiod. The experimental setup was observed for each population and replicated five times, which means there was a total of 100 D. citri nymphs investigated per population. (All experiments were run at the same time.)

DNA Extraction and Isolation
Before DNA extraction, all the T. radiata populations were surface-sterilized in 75% ethanol for 90 s, followed by three rinses in sterilized distilled water. For DNA extraction, each population's specimens were crushed separately using a DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer's guidelines. The quality of DNA was assessed by electrophoresis on 1% agarose gel and quantified using Nanodrop 1000 spectrophotometer (Thermo-Scientific, Waltham, MA, USA). DNA samples were stored at −20 • C for further use.

Amplification of DNA and Sequencing
The internal transcribed spacer (ITS1 and ITS2) regions and cytochrome oxidase subunit I (COI) was amplified using specific primers for the identification of T. radiata populations, while endosymbiont Wolbachia was confirmed using wsp gene-specific primers (Table 1). Each PCR reaction was carried out in a 25 uL reaction mixture containing 2 uL DNA template (100 ng/uL), 2 uL (10 mM) of each forward and reverse primers, and 12.5 uL of 2XTaq PCR Master mix (Tiangen Biotechnology, Beijing, China). The amplification protocol, annealing temperatures, and expected product size (bp) are summarized in Table 1 according to the literature [45][46][47]. The PCR products were visualized by 1% agarose gel electrophoresis and purified using purification kit QIAQuick (QIAGEN, Hilden, Germany). Sequencing of amplified DNA fragments was performed at BioSune Biotechnology Co., Ltd. (Shanghai, China).

Sequence Composition and Phylogenetic Analysis
Obtained sequences of ITS1 and ITS2 regions, COI, and wsp genes were first assembled with DNAMAN v5 (Lynnon BioSoft, Vaudreuil-Dorion, Quebec, QC, Canada). Assembled sequences were subjected to BLAST analysis (https://blast.ncbi.nlm.nih.gov/Blast.cgi, accessed on 1 May 2021) to determine their identity with those deposited in GenBank [46]. The sequences showing the highest similarity index were retrieved and aligned using the Clustal W program BioEdit v7 [48]. For COI and wsp genes, nucleotide sequences were translated to amino acid sequences with ExPASY tools [49]. Our ITS1, ITS2, COI, and wsp sequences of all T. radiata populations were deposited in the NCBI GenBank database and their accession numbers are shown in Table S1.
The phylogenetic analysis of ITS1, ITS2, COI, and wsp sequences of different populations of T. radiata was conducted using MEGA-X [50]. For additional support and diversity, already reported ITS1 sequences of T. radiata populations from Bhutan, Mexico, Vietnam, and USA (Florida and Texas regions) had been used [25,27]. Meanwhile, sequences from Bhutan, Finland, India, Puerto Rico, Vietnam, and USA (Florida and Texas regions) were employed for ITS2 analysis of T. radiata populations [25,27]. For phylogenetic tree constructions, the General Time Reversible under the initial tree option of maximum parsimony method was used for the maximum likelihood phylogeny [51], while the Maximum Composite Likelihood method was considered for Neighbor-Joining topology [52]. The heuristic search trees were obtained with 1000 bootstrap support [50]. Furthermore, MEGA-X software implemented the Maximum Likelihood model to calculate the genetic distances of the aligned sequences for province-wise and overall populations of T. radiata from China [50].

Genetic Diversity and Distribution Analysis
The nucleotide polymorphism evaluation parameters, including haplotype number (h), haplotype diversity (Hd), nucleotide diversity (π), and the average number of nucleotide differences (k) were calculated using DnaSP 5 within and overall populations of T. radiata collected from China [53]. The neutrality test indices, i.e., Tajima's D, Fu, and Li's D* statistical test values for deviating neutral evolution pattern, were calculated within and overall T. radiata populations. Further, population analysis was performed with the Raggedness indices model [53].

Data Analysis
The parasitism percentage was determined on the fifth day based on the number of nymphs parasitized. The percentage of parasitism was calculated as the ratio of the number of emerged T. radiata to the combined total number of emerged D. citri and T. radiata. The number of emerged T. radiata evaluated the emergence percentage divided by the number of parasitized T. radiata. The sex ratio is based on the number of females (T. radiata) divided by the total number of emerged T. radiata for each population. The differences in the parasitism, emergence, and sex ratio among different populations of T. radiata were analyzed using one-way analysis of variance (ANOVA) followed by Tukey's multiple comparison test with a significant difference at p < 0.05.

Results
In this study, nymphs of ACP were found to be parasitized with T. radiata on M. paniculata plants were collected in Fujian, Guangdong, and Jiangxi provinces in China, while a previously collected population of T. radiata from Guayas, Ecuador, maintained in insectary since 2018, was also included in analysis. All the T. radiata populations were kept in separate cages under controlled conditions and reared successfully for the present investigation. The study of T. radiata populations determines biological aspects, including parasitism rate, emergence rate, and sex ratio from each population by parasitizing third-fifth instar nymphs of D. citri. The obtained sequences of ITS1, ITS2, COI, and wsp were blasted in the NCBI database to identify populations of T. radiata and their associated Wolbachia strains. The blast results showed more than 97-99% sequence identity with their respective species or strains deposited in the NCBI database and represented in the phylogenetic trees.

Efficacy of Parasitism, Emergence, and Sex Ratio on a Different Population of Tamarixia radiata
The parasitism rate of different T. radiata populations after parasitization of D. citri third-fifth instar nymphs are provided in Table 2. The results showed a statistical parasitism percentage difference between China and Ecuador populations as well as within T. radiata populations from China (F 3.39 = 15.57, p < 0.001). The parasitism percentage of the Fujian population was significantly higher (68.20%) than the Jiangxi population (63.54%), whereas there was no significant difference observed between populations from Guangdong and Fujian provinces. In addition, the minimum parasitism (57.28%) was observed by the Guayas population of T. radiata. After the potential parasitism occurred post-T. radiata-release, the highest emergence rate (86.61%) was observed in the Fujian population. However, no statistical dissimilarities were observed in the emergence percentage among and in between China and Ecuador populations. The trends in the sex ratio of different T. radiata populations were substantial and similar to the emergence rate. Furthermore, the female proportion results showed that the F1 generation of Guayas' population was lower (58.40%), whereas there was a statistical similarity in the overall population of China and Ecuador (Table 2).

Sequence and Phylogenetic Analysis of the ITS1 Region
The length of the ITS1 sequences from all populations of T. radiata ranged from 503 to 515 bp. The alignment results showed maximum sequence similarity among Fujian-1, Guangdong-2, Guayas-1, and Jiangxi-1 populations with a difference of one nucleotide at their respective positions of 18, 22, 45, and 47. The phylogenetic analysis showed that T. radiata individuals collected from Fujian, Guangdong, Jiangxi, and Guayas formed a separate clade with the Bhutan and Vietnam populations of T. radiata ( Figure 1A,B). However, the Fujian-1 population (MW537748) is grouped with the population reported from Vietnam, while Tamarixia triozae Burks is grouped into a distinct clade. The pairwise distance divergence ranged between 0.0020% and 0.0048% among the populations obtained from Fujian, Guangdong, and Jiangxi ( Table 3) 86% similarity among all population sequences ( Figure 1A). In the Neighbor-Joining tree, all populations showed a resemblance of 96% with all T. radiata populations ( Figure 1B).
Agronomy 2021, 11, x FOR PEER REVIEW 6 of 18 divergence ranged between 0.0020% and 0.0048% among the populations obtained from Fujian, Guangdong, and Jiangxi ( Table 3). The 34 nucleotide sequences involved 563 positions in the final dataset at pairwise deletion and gamma distribution options, including 1000 bootstrap values. Furthermore, the Maximum Likelihood method presented 86% similarity among all population sequences ( Figure 1A). In the Neighbor-Joining tree, all populations showed a resemblance of 96% with all T. radiata populations ( Figure 1B).

Sequence and Phylogenetic Analysis of the ITS2 Region
Obtained sequences of the ITS2 region of parasitoid T. radiata ranged from 489 to 490 bp. A difference of nine nucleotides was observed among nucleotide sequences of Fujian-1, Guangdong-2, and Jiangxi-1. However, minimum nucleotide difference was observed within populations of Fujian, Guangdong, Jiangxi, and Guayas. All the collected T. radiata populations (Fujian, Guangdong, Jiangxi, and Guayas) are grouped into a single clade with other T. radiata populations collected from D. citri host insect insects from several geographical locations. Similarly, the other species of Tamarixia (Tamarixia drukyulensis Yefremova and Yegorenkova) clustered into a separate clade concerning their host insect (Diaphorina communis Mathur) (Figure 2A,B). According to the sequenced database, the maximum distance divergence ranged between 0.0043% and 0.0243% in Fujian, Guangdong, and Jiangxi populations, while minimum variation was observed within each group (Table 3). This analysis involved 44 nucleotide sequences with 561 positions in the final dataset at pairwise deletion and gamma distribution option, including 1000 bootstrap values. Moreover, according to the Maximum Likelihood method, all populations showed 100% resemblance with the entire group ( Figure 2A). The Neighbor-Joining tree presented up to 94% similarity among all group sequences ( Figure 2B).

Sequence and Phylogenetic Analysis of the ITS2 Region
Obtained sequences of the ITS2 region of parasitoid T. radiata ranged from 489 to 490 bp. A difference of nine nucleotides was observed among nucleotide sequences of Fujian-1, Guangdong-2, and Jiangxi-1. However, minimum nucleotide difference was observed within populations of Fujian, Guangdong, Jiangxi, and Guayas. All the collected T. radiata populations (Fujian, Guangdong, Jiangxi, and Guayas) are grouped into a single clade with other T. radiata populations collected from D. citri host insect insects from several geographical locations. Similarly, the other species of Tamarixia (Tamarixia drukyulensis Yefremova and Yegorenkova) clustered into a separate clade concerning their host insect (Diaphorina communis Mathur) (Figure 2A,B). According to the sequenced database, the maximum distance divergence ranged between 0.0043% and 0.0243% in Fujian, Guangdong, and Jiangxi populations, while minimum variation was observed within each group (Table 3). This analysis involved 44 nucleotide sequences with 561 positions in the final dataset at pairwise deletion and gamma distribution option, including 1000 bootstrap values. Moreover, according to the Maximum Likelihood method, all populations showed 100% resemblance with the entire group ( Figure 2A). The Neighbor-Joining tree presented up to 94% similarity among all group sequences ( Figure 2B).

Sequence and Phylogenetic Analysis of COI
The obtained sequence length of COI was 488-565 bp from each of T. radiata populations collected in Fujian, Guangdong, Jiangxi, and Guayas. Phylogenetic analysis showed that two variable sites corresponding to transitions and transversions were identified in the COI sequences of T. radiata populations clustered into two distinct clades based on the haplotype diversity (H1 and H2). In the tree topology, the first clade showed H-2 from the Fujian population grouped with the Mexican population of T. radiata. The second clade represents H-1, which is frequently encountered from China (Guangdong and Jiangxi) to Ecuador (Guayas) grouped with the Mexican and USA (Florida and Taxes) popula- tions of T. radiata ( Figure 3A,B). However, the distance divergence was calculated using the Maximum Composite Likelihood model, and variation ranged between 0.0021% and 3.05% (Table 3). The 32 nucleotide sequences with 628 positions were involved in the final dataset at pairwise deletion and gamma distribution options, including 1000 bootstrap values. Moreover, according to the Maximum Likelihood and Neighbor-Joining methods, all populations showed 100% resemblance with the entire group of T. radiata populations ( Figure 3A,B).
that two variable sites corresponding to transitions and transversions were identified in the COI sequences of T. radiata populations clustered into two distinct clades based on the haplotype diversity (H1 and H2). In the tree topology, the first clade showed H-2 from the Fujian population grouped with the Mexican population of T. radiata. The second clade represents H-1, which is frequently encountered from China (Guangdong and Jiangxi) to Ecuador (Guayas) grouped with the Mexican and USA (Florida and Taxes) populations of T. radiata ( Figure 3A,B). However, the distance divergence was calculated using the Maximum Composite Likelihood model, and variation ranged between 0.0021% and 3.05% (Table 3). The 32 nucleotide sequences with 628 positions were involved in the final dataset at pairwise deletion and gamma distribution options, including 1000 bootstrap values. Moreover, according to the Maximum Likelihood and Neighbor-Joining methods, all populations showed 100% resemblance with the entire group of T. radiata populations ( Figure 3A,B).

Phylogenetic Analysis of Wolbachia Based on the wsp Gene
The presence of Wolbachia as endosymbionts is well-known in the previous host-parasitoid relationship of D. citri and its endoparasitoid (Diaphorencyrtus algerihnses Shafee, Alam and Agarwal) [54]. However, the presence of Wolbachia in T. radiata has been confirmed using the wsp gene in the present study. Phylogenetic analysis showed that Wolbachia sequences obtained in this study clustered in a single sub-clade with already reported Wolbachia pipientis from Drosophila melanogaster and Aedes albopictus with maximum bootstrap support. Based on sequence alignment and phylogenetic analysis, we

Phylogenetic Analysis of Wolbachia Based on the wsp Gene
The presence of Wolbachia as endosymbionts is well-known in the previous hostparasitoid relationship of D. citri and its endoparasitoid (Diaphorencyrtus algerihnses Shafee, Alam and Agarwal) [54]. However, the presence of Wolbachia in T. radiata has been confirmed using the wsp gene in the present study. Phylogenetic analysis showed that Wolbachia sequences obtained in this study clustered in a single sub-clade with already reported Wolbachia pipientis from Drosophila melanogaster and Aedes albopictus with maximum bootstrap support. Based on sequence alignment and phylogenetic analysis, we could identify Wolbachia endosymbiont as W. pipientis from different populations of T. radiata. Previous research describes at least two Wolbachia strains (A and B groups) based on host-parasitoids relationships [55]. To evaluate the phylogenetic position and diversity of Wolbachia within T. radiata populations, we performed a detailed analysis using the wsp sequence data. Our study contained four 16s RNA used as an outgroup for Wolbachia sequenced from different hosts (D. algerihnses (EF433794), D. citri (EF433793, AB038370), Culex pipiens (U23709) according to Meyer and Hoy (2008)). Based on the Maximum Likelihood and Neighbor-Joining method, the phylogenetic tree clustered the entire Wolbachia sequences obtained Agronomy 2021, 11, 2018 9 of 17 from T. radiata populations into Subdivision-B [55], which showed a strong relationship with previously reported Wolbachia strains. Analysis results also showed that previous Wolbachia sequences (data from the literature) in D. citri also belong to Subdivision-B, as observed in T. radiata populations. Wolbachia sequences found in T. radiata and D. citri are closely related but not identical and therefore positioned into two separate subclades. Accordingly, the horizontal transfer might have occurred between host and parasitoid, but if so, it was not a recent one. Pairwise sequence divergence of 0.0035%-0.0117% was observed within the Wolbachia strains (Table 3). This analysis involved 77 nucleotide sequences with 670 positions in the final dataset at pairwise deletion and gamma distribution options, including 1000 bootstrap values. Moreover, in the Maximum Likelihood and Neighbor-Joining method, all populations assembled in both trees showed a resemblance of 97% and 99% within the group, respectively ( Figure 4A,B).

Population Genetic Diversity Analysis
All the populations of T. radiata were collected from different geographical regions and showed some genetic variation within and overall populations. The genetic diversity was calculated based on haplotype diversity (Hd), nucleotide diversity (π), and the average number of nucleotide differences (k) ranges from (0.500-0.833), (0.00193-0.00252), and (1.0000-1.50000) for the COI gene, respectively (Table 4). Likewise, the COI gene showed maximum nucleotide variations up to 243 mutations overall among the Chinese populations based on the Eta value and revealed a high genetic diversity (Table 4).

Demographic Analysis
The neutrality tests were performed to confirm the significance of genetic diversity within and overall populations of T. radiata and its associated endosymbiont W. pipientis. The results of the COI gene from both neutrality indices showed a significant positive p-value of Tajima's D (1.6893) and Fu and Li's D* test and (2.1857), respectively (Table 4). A positive value of these tests is evidence for a deficiency of alleles and an expected decrease in population size. The neutrality tests showed a non-significant negative value of neutrality indices analysis for the wsp gene (Table 4). These negative p-values indicate an excess of rare mutations that favor population expansion. These results may suggest an excess number of low-frequency alleles typical of positive selection or recent population expansion. DNA sequences from different genes were analyzed for population size changes to enrich the genetic diversity results among different T. radiata populations and their associated Wolbachia strains to observe pairwise nucleotide differences (mismatch distribution). The pairwise nucleotide differences graph showed a bimodal curve in the COI gene under a sudden contraction model that expected the population bottleneck ( Figure 5A). In contrast, the pairwise mismatch distributions plot, which comprises the W. pipientis strain of the wsp gene, was smooth and unimodal, which may indicate that a strong population subdivision confers a stable population size ( Figure 5B).

Discussion
The ACP is the primary vector of the gravest disease of the citrus industry, known as a citrus greening disease [56], and T. radiata is one of the recommended biological controls for sustainable management of D. citri in different geographic areas worldwide. It has been reported that T. radiata has a very high reproductive rate with a short life span and is frequently reared in control conditions [57]; nevertheless, little work has been completed in insectary to assess the parasitism rates and genetic diversity of different geographical populations of T. radiata from China. The objective of this study was first to examine the biological parameters of different populations of T. radiata, their genetic di-

Discussion
The ACP is the primary vector of the gravest disease of the citrus industry, known as a citrus greening disease [56], and T. radiata is one of the recommended biological controls for sustainable management of D. citri in different geographic areas worldwide. It has been reported that T. radiata has a very high reproductive rate with a short life span and is frequently reared in control conditions [57]; nevertheless, little work has been completed in insectary to assess the parasitism rates and genetic diversity of different geographical populations of T. radiata from China. The objective of this study was first to examine the biological parameters of different populations of T. radiata, their genetic diversity, and the presence of their associated Wolbachia species and its diversity in the various populations. T. radiata is a solitary, arrhenotokous ectoparasitoid of third, fourth, or fifth instars of D. citri nymphs as it can kill approximately 500 nymphs by oviposition and feeding [6,58,59]. The environmental factors may affect the parasitism percentage of T. radiata and directly influence the foraging behavior of parasitic wasps. Among these factors, temperature, hostage, and geographical locations may directly affect the parasitism and emergence of T. radiata [23]. The previous studies revealed that the T. radiata preferred older nymphs of D. citri (fourth to fifth) at a temperature range from 25 to 30 • C and reported the highest parasitism percentage (85.50%) and emergence rate (88.30%) [22]. The parasitism of older instars psyllid nymphs in orange jasmine ranging from 48% to 70% was recorded during spring and summer at Río Piedras and San Juan, respectively [21]. According to Wang et al. 1999, temperatures ranging from 20 to 30 • C were favorable for parasitism and reproduction of T. radiata. However, T. radiata may feed and kill all instars of D. citri nymphs but mostly preferred fifth instar host because the progeny sex ratio of parasitoids mainly influences the host instar stage and oviposition [60]. In the present study, the different populations of parasitoids performance (parasitism rate, emergence rate, and sex ratio) may provide insights into the behavior and interactions of the parasitoid concerning its geographical location. The current findings show that T. radiata preferred parasitized (third-fifth instar) nymphs of D. citri on M. paniculate plants under control conditions (27.5 ± 1 • C, 70 ± 5º% RH) and the maximum percentage of parasitism was observed from T. radiata populations from Fujian (68.20%) and minimum from Guayas (57.28%). Similarly, the rate of emergence may range from 86.61% to 82.95% for third-fifth instar psyllid nymphs. Likewise, results from the life table of T. radiata have reported a higher parasitism rate (77.24%) at a temperature of 26.3 • C. [23]. Previous research revealed that a series of complex environmental factors had significant effects on survival, parasitism, and reproduction of T. radiata. Accordingly, the highest parasitism percentage was observed in summer and the lowest in winter in Sao Paulo [61].
In the molecular study, the ITS (ITS1, ITS2) regions and COI were used to determine nucleotide diversity between and within populations of T. radiata sampled across China (Fujian, Guangdong, and Jiangxi) and one population from Ecuador (Guayas). We also amplified the wsp gene to confirm the identification and variability of Wolbachia strains as an associated endosymbiont in T. radiata populations. Molecular analysis of T. radiata is primarily related to genetic diversity within the species using conserved ITS regions and the COI gene [24,27,62]. Among the conserved regions, the ITS loci are useful diagnostic markers due to their high divergence levels between species but low levels of variation within a species [33,63]. The ITS regions have been used to study phylogenetic relationships and resolve taxonomic controversies in various insect species, including T. triozae, T. mercet, and T. drukyulensis [25,27]. In the current study, the ITS regions have been used to construct the phylogenetic tree with twenty-one (ITS1) and thirty-one (ITS2) reference sequences of Tamarixia reported from the different agro-ecological zones, including T. radiata sequence from China populations and one from Ecuador (Guayas). Sequence analysis of ITS1 and ITS2 showed no variability among the populations of T. radiata from China and Ecuador. The phylogenetic tree's topology based on ITS regions represented strong intermingled relationships between and within all the populations of T. radiata ( Figure 1A,B and Figure 2A,B). However, sequence analysis of ITS regions showed minor distance divergence for province and overall T. radiata populations from Fujian, Guangdong, and Jiangxi, ranging from 0.0028% to 0.0048% (ITS1) and 0.0043% to 0.0243% (ITS2), respectively ( Table 3). The minimum divergence in conserved regions in the T. radiata populations supports the that these different populations represent a single species.
In previous studies, mitochondrial genes have been used to analyze the genetic diversity and genetic differentiation of the Tamarixia species [24,25,27]. Proper assimilation and understanding of the arthropods' genetic variability found it essential to mitigate and improve its monitoring, which further facilitates the implementation of need-based management strategies. In this study, T. radiata populations' maximum nucleotide diversity was observed based on Eta mutations in the total COI sequences, as observed in the Chinese population ( Table 4). The previous studies also revealed two haplotypes (Hap-1 and Hap-2) for T. radiata in the Sonora region of Mexico and America (Florida, Texas) [24,27]. The previous analysis also revealed that H1 and H2 have been present in reported populations from Mexico, Texas, Florida, but the distribution and frequency bias by H2 was greater than that of H1. Moreover, haplotype H2 was found in all samples obtained from 11 states of Mexico, and H1 was only found in Yucatan [64]. COI gene analyses have also been widely used for the phylogeny of population and genetic inferences. Based on the COI phylogeographic analysis, it was observed that T. radiata populations from China showed two types of haplotype diversity (H1, H2). However, T. radiata populations from Fujian fall into separates clades based on the H1 type of haplotype ( Figure 3A,B). In contrast, the population of T. radiata from Jiangxi Guangdong and Guayas showed H2 haplotype diversity and is similar to the previously reported populations of T. radiata from Mexico and Florida. Furthermore, the previous work showed a high number of haplotypes H2 samples obtained from the northeast states of Nuevo León, and Tamaulipas and the western states of Colima and Michoacán, and H1 was only obtained from the Tamaulipas state of Mexico. These findings suggest gene flow due to the release and introduction of T. radiata between countries [61]. The high haplotype diversity and nucleotide diversity (Hd > 0.78, π > 0.24) observed in overall populations from China represent stable populations with a long evolutionary history [65], and these populations are the contributors to the nucleotide diversity of T. radiata. Likewise, the averaged maximum genetic distance deviation among all T. radiata populations was observed at 0.0053% and 3.0552% (Table 3). The neutrality indices of Tajima's D tests' positive values represent a lack of rare alleles, which may be due to balancing selection or population contraction. In contrast, the negative values indicate an excess of low-frequency alleles and represent a recent selective sweep or population expansion [66,67]. However, p-values of neutrality tests indicated a positive significant genetic diversity for total T. radiata populations, which may support population contraction, and the mismatch distributions could represent the bottleneck growth. The significant positive values of Tajima's D (2.1857) indicated that the population was undergoing balanced selection and had no further rapid expansion in the future. Additionally, Fu and Li's D* (1.6893) positive significance indicated that background selection had a significant effect on haplotypes of the total populations from China.
The facultative endosymbiont Wolbachia was previously reported in bacteriocytes and multiple somatic and reproductive tissues in various insect hosts and parasitoids [68]. Wolbachia has been found in multiple tissues, bacteriocytes, and somatic and other reproductive tissues in various insect pests [33,69]. The association of Wolbachia can affect its host insect's behavior and biology by manipulating their resident strains [70]. Therefore, the characterization of Wolbachia in T. radiata is of particular interest due to its impact on host biology and its potential to control D. citri. By analyzing Wolbachia infection status using the wsp gene within T. radiata populations and its distribution across the three provinces from China and one province from Ecuador, we found the presence of W. pipientis in association with T. radiata. Our results showed that W. pipientis obtained from T. radiata formed a clade with different insects, including different species of Drosophila, Aedes, Armadillidium, and Culex. These results revealed that the previous literature also showed that the W. pipientis host range is broader than initially thought [71]. However, we did not find any evidence of an association between Wolbachia and its host's geography. The phylogenetic analysis showed that the strains of W. pipientis found in T. radiata related to supergroup B ( Figure 4A,B) but distinct from D. citri positioned into a separate subclade. None of these results is surprising since the phylogeny of Wolbachia is known to be incongruent with that of its arthropod hosts, precisely because horizontal transfers between species are not uncommon at a large phylogenetic scale. However, it has been proposed that supergroup G be decommissioned, as it is based primarily on recombinant wsp sequences and clusters with A supergroup based on five multilocus sequence typing genes [72,73], and eight supergroups (A-H) are still widely used in the research community. An MLST system based on five house-keeping genes (coxA, gatB, hcpA, ftsZ, and fbpA) has been developed for Wolbachia [72] and is widely used for strain typing and to characterize strain variation within Wolbachia [74]. With the differentiation of Wolbachia from T. radiata and D. citri within the B subdivision, the idea of recent and rapid expansion of the B-clade Wolbachia could result from more frequent transfers of Wolbachia between host and parasitoid. In previous studies, rapid expansion has also been observed in Wolbachia strains associated with Drosophila species and their associated parasitoids [55,75]. However, the genetic diversity was calculated through the particular mutations in the wsp gene. The neutrality tests (Tajima's D, −1.1671 and Fu and Li's D*−1.1084) and sequence mismatch distributions analysis of total populations of Wolbachia strains were negative and not significant, indicating that the expansion of these populations was limited [66,76].

Conclusions
This study investigates the population genetic diversity of twelve populations of T. radiata from different provinces in China and one population from Ecuador. The results showed that different populations of T. radiata were distinguished into two clades of haplotypes (H1-H2). The populations of T. radiata from Fujian fall into H-1, and Jiangxi, Guangdong, and Guayas were H-2, the most frequent population clade. The total genetic diversity of T. radiata was maximum as indicated by the total haplotype diversity (Hd = 0.788), nucleotide diversity (π = 0.243), and average nucleotide difference (k = 1177). Tajima's D and Fu and Li's D* analyses indicated that these populations of T. radiata had population contraction and were limited by local area. The presence of Wolbachia strains from all T. radiata populations had a strong interaction with previously reported wsp gene strains. The negative values of the neutrality test may indicate long-term balancing selection or population expansion. The perceptive study of different T. radiata populations can be used to design any independent strategy to carry out integrated pest management of D. citri.   Table S1.