Evidence of weak genetic structure and recent gene flow between Bactrocera dorsalis s.s. and B. papayae, across Southern Thailand and West Malaysia, supporting a single target pest for SIT applications

Background Bactrocera dorsalis s.s. (Hendel) and B. papayae Drew & Hancock, are invasive pests belonging to the B. dorsalis complex. Their species status, based on morphology, is sometimes arguable. Consequently, the existence of cryptic species and/or population isolation may decrease the effectiveness of the sterile insect technique (SIT) due to an unknown degree of sexual isolation between released sterile flies and wild counterparts. To evaluate the genetic relationship and current demography in wild populations for guiding the application of area-wide integrated pest management using SIT, seven microsatellite-derived markers from B. dorsalis s.s. and another five from B. papayae were used for surveying intra- and inter-specific variation, population structure, and recent migration among sympatric and allopatric populations of the two morphological forms across Southern Thailand and West Malaysia. Results Basic genetic variations were not significantly different among forms, populations, and geographical areas (P > 0.05). Nonetheless, two sets of microsatellite markers showed significantly different levels of polymorphisms. Genetic differentiation between intra- and inter-specific differences was significant, but low. Seventeen populations revealed three hypothetical genetic clusters (K = 3) regardless of forms and geographical areas. The genetic structure of sympatric populations slightly changed during the different years of collection. Recent gene flow (m ≥ 0.10) was frequently detected whether samples were sympatric or allopatric. Ninety-five of 379 individuals distributed across the given area were designated as recent migrants or of admixed ancestry. As a consequence of substantial migration, no significant correlation between genetic and geographic distances was detected (R2 = 0.056, P = 0.650). Conclusions According to the 12 microsatellite variations, weak population structure and recent gene flow suggest that there is no status for cryptic species between B. dorsalis s.s. and B. papayae forms in Southern Thailand and West Malaysia. Both forms can be treated as a single target pest for the SIT program in an area-wide sense. Additionally, the result of species identification based on molecular data and morphological character are not congruent. The use of independent, multiple approaches in the characterization of the target population may ensure the effectiveness and feasibility of SIT-based control in the target area.


Background
The sterile insect technique (SIT) is a powerful biological control method. It relies on mating between sterile insects (biological control agent) and their wild counterparts in order to reduce reproductive potential. SIT is effective because it is a species-specific approach and one of the most environmentally friendly solutions to insect pest management. However, there are several steps: mass-rearing, sexseparation for male-only release, sterilization, marking, and mass-releasing [1]. SIT is widely included in areawide pest control programs against high-profile key pests of economic and biomedical importance. Examples of successful insect pest control programs using SIT include the eradication of the New World screwworm fly from the United States, Mexico, and Libya, the Mediterranean fruit fly from the northern part of Chile and the southern part of Peru, and the melon fly from Japan [2]. To maximize SIT effectiveness, the mating success between the sterile males and wild females should be enhanced. Consequently, these females will lay sterile eggs, leading to population disruption, which reduces overall pest damage.
Problems such as potential sexual isolation due to the existence of cryptic species and/or population isolation must be investigated before SIT application [3]. Cryptic species comprise two or more nominal species, which are sometimes morphologically indistinguishable. They may be recently or deeply diverged in sympatric or allopatric locales [4]. Resolution of cryptic species is needed for the management of many types of field operations related to pest control, conservation, and infectious diseases [4,5]. Ignoring cryptic species potentially undermines the SIT program due to ambiguous mating compatibility [4]. Mating incompatibility may contribute to the pattern of population isolation. An unsuccessful case of eradication using the SIT-based approach occurred for the New World screwworm [3] in Jamaica. Population isolation was present after several years of application, indicating the existence of a sexually cryptic species or a mating incompatibility between released sterile and wild insects [6].
The Bactrocera dorsalis species complex, a large group of nominal tephritid fruit flies, is of interest to scientists in several fields of study. Regarding morphological characters and host preferences, the B. dorsalis complex consists of almost 100 similar species [7,8] [7][8][9][10][11]. B. dorsalis s.s. and B. papayae were reported to be parapatric species. B. dorsalis s.s. is distributed from the far north to the northern part of the Malay Peninsula while the species range of the other extends throughout the northern Malay Peninsula all across to the southern end, overlapping the range of B. dorsalis s.s around the Isthmus of Kra [7][8][9][10][11].
Allopatrically, B. carambolae and B. philippinensis are separately found in the Indonesian archipelago and the Philippines, respectively. Focusing on only Southern Thailand and the Malay Peninsula areas, the delimited distribution of B. dorsalis s.s. in the Malay Peninsula is uncertain [7][8][9][10]. According to only morphological forms, B. dorsalis s.s. and B. papayae are still classified as distinct species within the complex, although species limit has often been ambiguous [12,13]. On the other hand, the classification of B. papayae has recently been revised to be synonymized with B. philippinensis [8]. Nowadays, research into the species complex status has arrived at the question regarding the actual number of true economic species in such a complex, which is important for pest quarantine and management, as well as international trade [11].
Non-morphological characters (molecular markers, mating behaviour, and male pheromones), combined with morphological characters, have provided better resolution in the analysis of species limitation within several cryptic species [4]. Similarly, the biological status of members of B. dorsalis complex -B. dorsalis s.s., B. papayae, and B. phillipinensis -has recently been analyzed using independent multiple approaches such as studying several genes of mitochondrial DNA, microsatellite markers, and mating competitiveness [13][14][15][16]. This evidence suggests that they are the same entity. Among several molecular traits, the genetic diversity of microsatellite DNA could be a powerful tool for the study of species complex [17]. Strong inferences with less bias can be made when the microsatellite markers derived from both of the cryptic members are used. In addition, population samples are frequently collected throughout their geographical distribution in order to observe cline.
The current research aims to use microsatellite markers for generating high-resolution population genetic data from morphological forms -B. dorsalis s.s. and B. papayaespanning from the top of Southern Thailand to the tip of the Malay Peninsula (West Malaysia). Population sampling was carried out using intervals of no more than 200 km. Sympatric populations of the two forms were also collected.
To avoid the negative effect from an ascertainment bias of microsatellite DNA variation, we used microsatellite DNA markers derived from both B. dorsalis s.s. [31] and B. papayae [18] for surveying the intra-and inter-specific variation, genetic structure, and population dynamics. Furthermore, these data were subsequently used for evaluating the feasibility of SIT application in given areas. The degree of interbreeding between the two morphological forms can also be inferred.

Fruit fly samples and genomic DNA extraction
Males of Bactrocera dorsalis s.s. and B. papayae were collected using methyl eugenol traps. Seventeen fruit fly populations corresponding to 12 areas -spanning from the top of Southern Thailand to the tip of the Malay Peninsula (West Malaysia) -were collected at no more than 200 km intervals (Table 1, Figure 1). The identification of fruit fly forms was carried out based on published descriptions by Drew and Hancock [7]. Three-hundred and seventy-nine individuals, with conforming morphological characters, were subsequently genotyped and analyzed.
Fruit fly samples were preserved in 95% ethanol and kept at −20°C until use. Total genomic DNA was individually extracted from each fly [31].

Microsatellite amplification and genotyping
Sets of seven (Bd1, Bd9, Bd15, Bd19, Bd39, Bd42, and Bd85B) and five (Bp58, Bp73, Bp125, Bp173, and Bp181) microsatellite loci -previously isolated and characterized from B. dorsalis s.s. [31] and B. papayae [18], respectivelywere used (Additional file 1: Table S1). Microsatellite DNA amplifications were set up in a 15-μl volume reaction containing 100 ng of genomic DNA, 1 × buffer, 2.5 mM MgCl 2 , 25 μM dNTPs, 0.5 U Taq polymerase (Vivantis) and 5 μM of each primer. PCRs were performed using a thermal cycler Flexcycler (AnalytikJena, Germany) using the following conditions: 5 min at 94°C, 29 cycles of 30 s at 94°C, 90 s at T a of each primer pair (Additional file 1: Table S1) and 90 s at 72°C, and an additional 5 min of elongation at 72°C at the end of the process. Electrophoresis and allele scoring were determined in 6% or 12% (w/v) non-denaturing polyacrylamide gel in 1xTBE buffer at 800 V for 6 or 10 hrs, respectively, stained with 0.5 mg/l ethidium bromide [32] and photographed under UV light. PCR product size was measured with a 25 bp DNA ladder (Promega, USA). For each locus analysis, samples with no PCR product were declared null allele (missing data), if two PCR attempts had been carried out.

Data analysis
Test of the ascertainment bias hypothesis: Basic genetic variability (i.e., variance in PCR product size range (V m ), number of allele (n a ), and expected heterozygosity (H E )) was estimated for each population. However, in this study, we used the effective number of allele (n e ), instead of n a because it is less sensitive to sample sizes and rare alleles [33].
The n e was calculated as n e = 1/(1-H E ). H E was calculated as H E = 1-(Σq i 2 ), where q i is the frequency of the i th allele in the population. V m and n e were in an approximately normal distribution. Each of the two measurements was analyzed by ANOVA using these factors: microsatellite DNA loci derived-species, B. dorsalis s.s. and B. papayae, populations, and countries. Conversely, H E was transformed to 1/H E and subsequently tested with the non-parametric Mann-Whitney U test with the following factors: microsatellite DNA loci derived-species, B. dorsalis s.s. and B. papayae, and countries. Only the population factor was tested with a Kruskal-Wallis test. All statistical analyzes were performed using PASW statistics software v18.0 (ó SPSS). The significant level is below 0.05 (P < 0.05).
Genetic diversity: The parameters for population genetic analyzes, i.e., n a , n e , number and frequency of private alleles (n p and A p , respectively), observed heterozygosity (H O ), H E , and inbreeding coefficient (F IS ), were estimated using GENALEX v.6.5 [34]. In addition, rarefaction allele and private allele richness (R s and R p , respectively) were estimated using HP_Rare [35]. The frequency of null alleles (A n ) was estimated following Brookfield [36]. Deviation from Hardy-Weinberg equilibrium and linkage disequilibrium was determined using GENEPOP v.4 software [37] together with their critical levels after the sequential Bonferroni test [38].
Population genetic structure: We evaluated genetic structure using four different approaches: (i) measuring genetic differentiation (F ST ) among populations, (ii) Bayesian model-based clustering, (iii) the Principle Coordinate Analysis (PCoA), and (iv) a hierarchical AMOVA. The first approach quantified genetic differentiation, pairwise F ST , among 17 populations using MICROSATELLITE ANA-LYSER (MSA) [39]. Estimation of F ST values and their statistical significance was done for 10,000 permutations.
The Bayesian approach was used to determine genetically distinct groups (or clusters) using the program STRUC-TURE v.2.3.1 [40,41]. To identify the number of genetic clusters within a location, in particular with low sample sizes from some populations, the no admixture model was initially studied [5]. This model hypothesizes that each individual belongs to one cluster. STRUCTURE was run with the Infer Lambda option set at the number of clusters (K) equal to one and five independent iterations. Subsequently, the mean of five lambda values was set with the option of correlated allele frequency for all additional runs. The other parameters were set at default: different values of F ST for different subpopulations, prior F ST mean of 0.01, and a standard deviation of 0.05. We ran STRUCTURE using K = 1 to K = 17. For each K value, five iterations were performed with the condition of a burn-in period of 100,000 iterations followed by a run of 500,000 Markov Chain Monte Carlo (MCMC) repetitions. The optimal K value was estimated by examining the Ln P(X/K) output from STRUC-TURE [40] and calculating the ΔK statistic [42].
Also, to recognize the potential admixture between genetic clusters, the admixture model was additionally run with the same parameters in the no admixture analysis. . Information for each sample is described in Table 1. The Isthmus of Kra (connecting zone between yellow and green areas) is the putative transition zone between Bactrocera dorsalis s.s. (yellow area) and B. papayae (green area) distribution [7,11]. Intra-and inter-specific differences (F ST ) among nearby populations are reported. The red line represents intra-specific difference (pairwise F ST between populations with the same morphological form) while the green line represents inter-specific difference (pairwise F ST between populations with different morphological forms). An asterisk (*) indicates a non-significant F ST value. Rectangular boxes detail the designated locations (by arrow) where either sympatric or temporal populations were collected.
The location prior option (or sample assignment prior in this case) was included in order to allow for revealing weak population structure, but does not find structure when none is present. Moreover, it is able to ignore prior assignment when a correlation between the clusters and sampling locations is not observed [43]. All data were summarized using CLUMPP v.1.1.2b [44] and visualized through DISTRUCT v.1.1 [45].
The PCoA using GENALEX v.6.5 [34], performed on genetic distance, was used to display genetic divergence among the individual fruit flies in multidimensional space using allele frequency data. A plot of the first three principal coordinates was constructed using the subprogram MOD3D in NTSYS-pc v.2.1 [46].
A hierarchical spatial AMOVA was performed using ARLEQUIN v.3.1.1., with 1,000 permutations [47]. Populations were grouped corresponding to three major criteria, i.e., morphological form, geographical area, and population genetic structure, to test genetic homogeneity in different hierarchies.
Analysis of genetic distances: The MSA program [39] was also used for estimating genetic distance based on Nei's genetic distance [48] and the proportion of shared alleles [49]. The programs NEIGHBOR and CONSENSE in the PHYLIP package [50] were used to reconstruct the neighbor-joining trees after 1,000 bootstraps of the original data. Subsequently, the TreeView program was used for phylogenetic tree visualization [51].
Demographic inferences: Information regarding demography was investigated through three analyzes: (i) assignment test, (ii) population contraction/expansion, and (iii) isolation by distance (IBD). GENECLASS v.2.0 [52] was run to estimate the probability of each individual belonging to its own population, the probability of being an immigrant from each of the other populations, and the probability of being a migrant to other populations. A standard criterion, the Bayesian method [53], with enabled probability computation and Monte-Carlo resampling, following the simulation algorithm for population assignment was used [54]. An arbitrary threshold probability value of 0.100 was determined after simulating 10,000 genotypes for each population. BOTTLENECK v.1.2.02 [55] was run to detect the signal of demographic expansion/contraction in each population. The heterozygosity excess was tested under two proposed models of microsatellite mutation: the stepwise mutation model (SMM) and the two-phased mutation model (TPM). The latter model comprised 90% single-step and 10% multiplestep mutation. To avoid the effect of too few individuals and loci tested per population, Wilcoxon signed-rank test was performed [55]. Finally, the IBD analysis was performed using the ISOLDE option in the GENEPOP package to find the correlation between genetic and geographic distances [39].

Genetic variability
Hardy-Weinberg exact tests were performed for 12 microsatellite loci. After sequential Bonferroni correction [38], 55 out of the 204 population by locus comparisons significantly departed from Hardy-Weinberg expectations. However, all deviations were not concentrated in any population or at any locus. No significant linkage disequilibrium was detected between genotypes at the 12 loci.
All microsatellite loci reveal different levels of polymorphism among populations as summarized in Table 2. Within each population, the mean of R s values is 3.23 ± 0.64, with the highest value observed in RN population (3.51). The n p values are detected in all samples, except for three populations (ST, NSa, and PK1b). They range from one (KD and SL) to seven (SK2) with a low average frequency (0.02 and 0.06, respectively). JH shows the highest A p value (0.13) at a low number (n p = 3). Likewise, using the rarefaction approach, the mean R p value is 0.15 ± 0.19, with the highest value detected in SK2 (0. To identify the number of hypothetical genetic clusters, STRUCTURE analyzes were run using both no admixture and admixture models. The ΔK value indicates that K equals three (K = 3) is the optimal number of hypothetical genetic clusters in both no admixture and admixture models ( Figure 2). However, members in each genetic cluster are not related to any morphological forms. At K = 3, nine populations (i.e., RB, PK1a, PK2a, PK2b, RN, NSb, SK1, TR, and JH) reveal an admixed structure (Figure 3). Populations ST, TR, SL, and PH are separated from the rest with a major proportion of co-ancestry (Q) ranging from 0.789 (TR) to 0.997 (ST). Within the Prachub Kirikan location (PK code), most individuals are admixed. Such samples collected from PK1a, PK1b, and PK2a generally share their proportion of ancestry in the same clusters whereas PK2b samples appear to be admixed to a different cluster (Figure 3, Additional file 2: Table S2).
Principle Coordinate Analysis (PCoA) demonstrates the genetic divergence of fruit fly populations in multidimensional space (Figure 4). This result is also consistent with STRUCTURE analyzes. The first axis accounts for 30.25% of total variation, which separates populations ST, TR, SL, and PH from the remaining populations. The second (23.08% of total variation) mainly distinguishes SK1 and SK2 from the others. The third (16.40% of total variation) distinguishes KD and JH from the rest of the populations.
The main feature of neighbor-joining trees, based on Nei's genetic distance [48] and the proportion of shared alleles [49], is a clear-cut separation of four populations (ST, TR, SL, and PH) from the others. Additional neighborhoods include SK1 & SK2 and NSa & NSb ( Figure 5).
The analysis of molecular variance (AMOVA) shows the extent of genetic variation in different hierarchies. The populations were grouped by three criteria: morphological forms, geographical areas, and genetic structure (Table 4). In the overall analyzes, less than 5.0% of variation is attributed to the differences among groups while approximately 90% of variation is attributable to differences within populations. Non-significant difference among groups is observed with regard to morphological forms (scenario 1: B. dorsalis s.s. vs. B. papayae, P = 0.555). Likewise, geographical considerations -scenario 2 (Southern Thailand vs. West Malaysia) and scenario 3 (above vs. below the Isthmus of Kra) -reveal non-significant differences (P = 0.058 and P = 0.143, respectively). However, significant differences are observed (P < 0.050) when all samples were grouped by the genetic coancestry clusters at K = 2 or 3.

Demographic analyzes
Recent migration: The individual assignment analysis (m) was performed using GENECLASS 2.0 [52] as shown in Table 5

Figure 2
Log-likelihood probability (Ln P(X/K)) and the Delta K values of data. Three is indicated to be the most likely number of hypothetical genetic cluster (K) using admixture model.
In the case of NSa and NSb, the other sympatric populations, an asymmetric migration rate is detected from NSb to NSa (m = 0.611). We combined the results of two programs (GENE-CLASS and STRUCTURE) in order to interpret congruency between genetic data and morphological forms. According to the STRUCTURE analysis, 95 individuals (approximately 25% of the total) are categorized as admixed based on their proportion of Q that ranges from 0.200 to 0.800 (Additional file 2: Table S2). Five populations, including ST, NSa, SK2, SL, and KL, do not contain any of these admixed individuals (Figure 3). Based on the results from the GENECLASS program, 45 of 95 individuals appear to be nonimmigrant considering the consistency between the original sampling site and the most probable population. The remaining 50 individuals are significantly classified as migrants (m ≥ 0.100) from at least one population. Individual numbers 10, 14, 16, 33, and 43 are potentially admixed and are also migrants from elsewhere (Additional file 2: Table S2).
Bottleneck: Under the Stepwise Mutation Model (SMM), a significant heterozygote deficit (population expansion) was observed in three populations (PK1a (P = 0.043), PK2b (P = 0.001), and SK2 (P = 0.003)) and no significant heterozygote excess (population bottleneck) was detected based on a two-tailed Wilcoxon signed-rank test. Using a more stringent model, the two-phase model (TPM), only two populations (PK2b (P = 0.0034) and SK2 (P = 0.0342)), still showed a recent population expansion.  . papayae assigned to two and three genetic clusters (K = 2 and K = 3, respectively). Each horizontal stripe represents an individual. Each color represents the proportion of membership with regard to the each hypothetical genetic cluster. Five replicates were combined into one figure using CLUMPP [44] and DISTRUCT [45]. Arrows indicate admixed individuals with a mean proportion of genetic cluster (Q) between 0.200 to 0.800, most present in samples from Southern Thailand.

Discussion
A fine-scale population genetics study of two morphological forms -B. dorsalis s.s. and B. papayae -was carried out in areas spanning from the top of Southern Thailand to the tip of the Malay Peninsula using high-resolution microsatellite DNA markers derived from each form. It can be inferred that these two morphological forms comprise a panmictic unit. Genetic variation, population structure, and recent population demography suggest a high feasibility for area-wide integrated pest management using sterile insect technique (AW-IPM-SIT) control programs. A combination of population genetic tools and morphological characterization may be necessary to better understand target pest populations.

Comparable intra-and inter-specific variations
The factors affecting microsatellite DNA variability include locus-and genome-specific mutation rates [56]. Microsatellite DNA markers are always isolated and characterized from a single species. In such, specific loci with high genetic variability are chosen. However, this genetic variability may not be inherently transferred when the marker system is applied to another species. If the recipient genome is very different, the level of genetic variability of the orthologues will be biased (ascertainment bias) [28][29][30]. In this study, we did not directly analyze the mutation rates at each microsatellite locus, but we statistically compared genetic variation using common parameters (i.e., effective number of alleles, variance in allele size range, and heterozygosity) among 12 loci for 17 populations. We found that the given microsatellite loci were highly polymorphic for all populations, although the B. papayae-derived loci were less variable than B. dorsalis s.s. for all parameters. The explanation involves several factors that influenced microsatellite mutation rate, such as repeat number, repeat type, flanking sequence (GC content), chromosome location, and base substitutions in the microsatellite arrays [57]. In our case, all microsatellite DNA loci used have similar repeat numbers, the same type of GT/CA motif (Additional file 1: Table S1) and lack data regarding chromosome location. The different level of variability between two sets of microsatellite loci may be influenced by flanking sequences and base substitution in the microsatellite motifs. The low GC content of flanking sequences may have effect on high Figure 5 Neighbor-joining trees based on the genetic distance derived from (a) Nei's genetic distance [48] and (b) the proportion of shared alleles [49]. The number at each node indicates the bootstrap percentile values after 1,000 replications. microsatellite variability [58]. In this case, primer sequences In addition, no significant genetic differences were observed among all populations, the two forms, the two countries, and also the interaction between marker derivation and species factors. Consequently, it can be deduced that the genetic polymorphisms of all tested loci may be shaped by the issue of population histories, not by ascertainment bias. Therefore, using two different species-derived sets of microsatellite loci, if available, could help us to clarify the source of genetic variability as well as to avoid biasing the data. This may help us to not undermine the determination of sympatric species that have no interbreeding between groups [4]. Studies using high-resolution microsatellite markers had been successfully used to delineate the limits of other species complexes (e.g., [5,19,59]); however, our study does not support the hypothesis that B. dorsalis s.s. and B. papayae forms are cryptic species in the investigated areas. Independent studies using different samples and sets of microsatellite markers (i.e., microsatellite-derived markers from only B. dorsalis s.s. [14] and a combination of microsatellite-derived markers from B. dorsalis s.s. & B. papayae (this study)) present a comparable genetic variation between two morphological forms. These facts imply that F ST of the given species cannot be resolved. A plausible explanation is that both forms comprise a panmictic unit or recently diverged, and/or they are connected through a high level of hybridization [4]. If both species have recently diverged, intra-specific microsatellite variation should correlate with inter-specific variation as a result of incomplete lineage sorting or still maintain a low level of gene flow. We found that an observed intra-specific difference was generally comparable to an observed inter-specific divergence. Moreover, other data derived from various DNA barcodes such as COI [13], cox1, nad4-3′, CAD, period, ITS1, ITS2 [16], EF-1α, PER [60], and other non-morphological characteristics such as pheromone profile [61], and mating competition between both species [15], along with morphometric analyzes [13,14], provide congruent results, supporting the recent divergence hypothesis. Nonetheless, microsatellite DNA markers may provide additional insights into population genetics.
Other evidence of F ST between two morphological forms and among populations was significant, but quite low. Approximately 2% to 18% of genetic diversity is the result of genetic differentiation among B. dorsalis s.s. and B. papayae populations. A high level of gene flow between populations or recent species divergence is the most likely explanation. Considering the genetic diversity of the single species B. dorsalis s.s., using at least the same seven microsatellite loci as the current study, we found that approximately one to  six percent (1-6%) and up to 30% of genetic diversity are the consequence of genetic differentiation among the populations at the level of micro-geography (populations collected from Thailand) and macro-geography (covering populations collected from the hypothetical area of origin and other areas), respectively [62]. Comparing the two morphological forms and the single species B. dorsalis s.s., the F ST of former falls between micro-and macro-geographical values. However, this interpretation must be viewed with caution because the F ST value is independent of the particular characteristics of individual loci or alleles and is influenced by the geographical differences of sampling locations [63]. Finally, the AMOVA result indicated that there is no genetic heterogeneity when all samples were grouped based on morphological and geographical criteria. We can infer that the genetic variation detected from all samples belong to the population level but not the cryptic species level. For that reason, we have no strong evidence to disprove the hypothesis that both morphological forms are actually a panmictic unit [13][14][15][16]60] or very recently divergent. More sample collections from other geographical areas or different ecological niches may be required in order to generalize the status of the population genetic divergence, although other diagnosable characters have not yet evolved or been established.

Population dynamics in allopatric and sympatric areas
We inferred population dynamics in allopatric and sympatric areas with regard to recent migration and population expansion/contraction. The narrowest part of the Thai-Malay Peninsula -Kra Isthmus -is located between Ranong and Chumphon Provinces, Thailand. Despite the fact that this area was proposed to be a biogeographic transition zone for B. dorsalis s.s. and B. papayae [8,11], our data do not indicate any genetic barrier or genetic heterogeneity of both forms across this zone (congruent with [14]). In addition, the RN population appears to be a genetic source for the morphological forms of B. dorsalis s.s. and B. papayae in our studied areas. The data shows a high level of genetic variation, a low level of genetic differentiation, and recent asymmetric migration to almost all populations. Similarly, Ratchaburi was previously proposed as a genetic source of B. dorsalis s.s. in Thailand when samples from Southern Thailand and West Malaysia were not taken [62]. However, Ratchaburi is such a legitimate area because there are plenty of commercial fruit orchards and distribution centers. The current study illustrates that RB has as high a genetic variation as RN but serves as a recent genetic source for a few populations in Southern Thailand. The different types of preferred hosts and climates between the south of Southern Thailand and northward [10] may be the underlying factors affecting the different population dynamics. On the other hand, ST and TR may indicate recently introduced populations. These populations have a low level of genetic variation, clearly indicated by the low number and frequency of private alleles. These populations have recently received genetic information from other populations, but not vice versa. The genetic ancestry of a few fruit fly samples (1.32%) cannot be traced back to any investigated populations, although nearby populations were no more than 200 km away. Adequate representative populations could be sampled in order to expand on this type of population dynamics detail, which may help the implementation of SIT programs. Population dynamics is also illustrated using the data from two different collections at two different time points, but in the same area of Prachub Kirikan Province. Populations PK1a and PK2a, the same morphological form B. dorsalis s.s., seem to maintain their stable populations, which is inferred from the comparable population size and genetic variability. On the other hand, populations of the other morphological form appear to be expanding. The first collection, PK1b, has a small population size (N < 10); however, two years later, PK2b was observed to have increased in population size as well as in genetic variability. According to inferences made from allopatric and sympatric populations, the fruit fly population dynamics may not be at equilibrium in Southern Thailand and West Malaysia. This demonstrates that population genetic study using multiple-time-point sample collection is encouraged to infer delicate population dynamics processes.
Feasibility of an AW-IPM-SIT program for B. dorsalis s.s. and B. papayae in Southern Thailand and West Malaysia Fruit fly surveillance using a combination of tools (e.g., monitoring traps, species identification, and fruit sampling) is a key concept for planning and implementation of an AW-IPM-SIT program. Utilization of the traditional approach based on morphology and type of fruit fly hosts for characterization of pest species may result in misidentification. Sample collection from significantly isolated geographical areas without intermediate sites could provide samples with more discrete morphological characters than geographically closer sites [14]. Although our genetic study does not agree with the traditional morphological-based taxonomic species status of B. dorsalis s.s. and B. papayae forms, it is consistent with other independent identification approaches using non-morphological characteristics and morphometric analyzes (as mentioned before). This implies that morphological variation is not a standalone indicator of species boundaries. Extreme environments can impose stabilizing selection on morphological characters which have nothing to do with the species differentiation process. In contrast, selection of non-morphological traits such as behavior, molecular markers, or reproductive character, can accompany speciation [4]. Incorporation of data from non-morphological-based approaches (such as population genetics) may supply strong evidence that can describe species status.
The release of sterile male only is a very crucial step to improve SIT, in terms of the cost and biological efficiency in the field [2]. Development and evaluation of genetic sexing strains (GSSs), separating males from females, is still a challenge for researcher. In the tephritid fruit fly, several GSSs have been developed using classical and transgenic technology [64]. Over the last two decades, few of them have been available due to genetic instability, poor mating performance, and delayed regulatory approval, for both mass-rearing and field application. There are Ceraitis capitata (Vienna 8) [65], Bactrocera dorsalis (Salaya1) [66], and Anastrepha ludens (Tapachula-7) [67]. Therefore, it is very importance that one strain can be used for controlling other members of the same complex. Whether or not both of the current morphological forms are definitely defined as distinct taxa, releasing the sterile Salaya1 strain to control B. dorsalis s.s. and B. papayae is possible, in an area-wide sense, in Southern Thailand and West Malaysia. Several reliable lines of evidence (e.g., mating competitiveness field cage tests [15] and estimation of gene flow in this study) support that B. dorsalis s.s. and B. papayae forms may be a single target pest. However, mating competitiveness field cage tests need to be done for additional confirmation before actual implementation in the studied areas. A successful showcase study of released sterile B. dorsalis s.s. (originally derived from Hawaii) against B. carambolae in Suriname had been reported [68].
However, we detect ongoing delicate population dynamics processes such as population reduction/expansion and migration within studied areas. The phenomena could produce subpopulations by various scenarios and subsequently impact the effectiveness of AW-IPM-SIT. For example, new immigrants may be introduced into a new microhabitat which may form a pocket population. When the SIT activity is less intense or no longer practiced, they may be founders for an incursion scenario [69]. This reminds us that before and during the implementation of an AW-IPM-SIT program, a population genetic survey is highly recommended, especially when cryptic species and/ or population isolation issues are involved.

Conclusions
In summary, there is no status for cryptic species between two morphological forms -B. dorsalis s.s. and B. papayaein Southern Thailand and West Malaysia based on the variations of microsatellite DNA markers derived from both species. Hence, both forms may be treated as a single target pest for an SIT control program. However, resolution of genetic isolation and morphology are not congruent in species identification. The characterization of a pest population using multiple approaches may ensure the effectiveness and feasibility of the SIT-based method.