Nine wild fruit fly populations were collected from four geographical areas: Indonesia (6), Malaysia (1), Thailand (1), and Suriname (1) (Table 1 and Figure 1). These populations were from hosts and locations within the known ranges of Bactrocera carambolae (http://www.cabi.org/isc/datasheet/8700). In Indonesia, six populations were collected from three main islands; two populations (i.e., North Sumatra-NS and Pekanbaru-PK) were sampled from Sumatra; three populations (i.e., Depok-DP, Jakarta-JK, and Bandung-BD) were sampled from Java; and another (West Kalimantan-WK) was sampled from Borneo. In Thailand, one population was collected from the southern region (Nakhon Sri Thammarat-NT). All fruit fly samples were firstly characterized as B. carambolae based on Drew and Hancock (1994). In addition, the male pheromone profile (Wee and Tan 2005b) and/or ITS1 marker (Armstrong and Cameron 2000, Boykin et al. 2014) were also used to confirm the characterization (Table 1). Only populations BD and WK showed conflicting identifications and were classified as unidentified populations.

Figure 1. 

Sampling collections of Bactrocera carambolae and B. dorsalis in this study. Seven populations of B. carambolae (blue dots) were collected from Southeast Asia and Suriname. Three populations of B. dorsalis (red dots) were sampled from East and Southeast Asia. Two other unidentified populations (purple dots) were included. Information for each population is described in Table 1.

Three other populations of B. dorsalis were included in this study as outgroup samples for the investigation of genetic relationship between two cryptic species. These populations were collected from the known distributions of B. dorsalis (http://www.cabi.org/isc/datasheet/17685) and characterized as B. dorsalis using the same methods described before (Table 1). One population is from Ratchaburi, Thailand (RB), which is a representative source of B. dorsalis in Southeast Asia (Aketarawong et al. 2007, 2014a). Another is from the northern part of Thailand, Chang Mai (CM). The other was collected from Kaohsiung, Taiwan (KS). This sample is supposed to have been introduced from mainland China and has become an isolated population (Aketarawong et al. 2007, 2014a).

To record the genetic relationship between the genetic sexing Salaya5 strain (Isasawin et al. 2014) and wild populations, a sample of the Salaya5 colony (SY5) was included. The Salaya5 strain was created by hybridization and introgression of the genetic sexing B. dorsalis strain, named Salaya1 (Isasawin et al. 2012), and B. carambolae from Jakarta, Indonesia. This strain has a genetic background close to B. carambolae (99.9%) and a part of the Y-pseudo linked autosome carrying a dominant allele of white pupae alleles of Salaya1 strain (Isasawin et al. 2014). The strain was confirmed to be B. carambolae described by Isasawin et al. (2014). The Salaya5 colony has been maintained in the conditions presented in Isasawin et al. (2014).

All samples were preserved in 95% ethanol and kept at -20 °C until use. The genomic DNA of each fly was extracted using the method of Baruffi et al. (1995).

Twelve microsatellite loci (Bd1, Bd9, Bd15, Bd19, Bd39, Bd42, and Bd85B derived from B. dorsalis s.s. (Aketarawong et al. 2006), and Bp58, Bp73, Bp125, Bp173, and Bp181 derived from B. papayae (Shearman et al. 2006)), previously established for study of members in the same complex (Aketarawong et al. 2014b), were analyzed in B. carambolae. Amplifications were set up in a 15- μl volume reaction containing 1 × buffer, 2.5 mM MgCl₂, 25 μM dNTPs, 0.5 U Taq polymerase (Vivantis), 5 μM of each primer, and 100 ng of genomic DNA. PCRs were performed using the thermal cycler Flexcycler (AnalytikJena, Germany) using the conditions described by Aketarawong et al. (2014b). Amplicons of sizes consistent with fragments of B. dorsalis were cloned and sequenced using an ABI PRISM 310 genetic analyser (Macrogen, Korea). All sequences, except Bp173, showed homology to the original sequences in Genbank (Table 2).

To improve null allele problems due to mutations on primer-binding sites and/or unsuitable PCR conditions, a new set of primers were designed and renamed for the 11 loci using OLIGO version 4.0-s (Rychlik and Rhoads 1989). Moreover, the new annealing temperature for each primer pair is shown in Table 2.

Fifteen flies from Jakarta and North Sumatra, Indonesia were initially screened with the 11 sets of new primers using the PCR conditions mentioned above. Electrophoresis and allele scoring were determined as in Aketarawong et al. (2014b). Only eight loci were selected based on sharpness, specificity, and polymorphism of amplified products for further testing in all samples (Table 2).

The descriptive parameters of population genetics were estimated using GENALEX v.6.5 (Peakall and Smouse 2012). These parameters include the mean number of alleles (n_a), mean effective number of alleles (n_e), mean number and frequency of private alleles (n_p and A_p, respectively), and mean observed and expected heterozygosity (H_O and H_E, respectively). Null allele frequency (A_n) was estimated following Brookfield (1996). Departure from Hardy-Weinberg equilibrium and the effect of linkage disequilibrium were tested using GENEPOP v.4 (Rosset 2008), with their critical levels set according to the sequential Bonferroni test (Rice 1989).

Genetic differentiation (F_ST) among 13 populations was measured using MICROSATELLITE ANALYSER (MSA) (Dieringer and Schlötterrer 2003). In addition, genetically distinct groups (or clusters) were determined using the Bayesian approach implemented in STRUCTURE v.2.3.1 (Pritchard et al. 2000, Falush et al. 2003). The admixture model (the F model), assuming correlated allele frequencies, was run. The program computed the number of possible clusters (K) from one to 13, with the condition of the burn-in period being 100,000 steps, followed by 500,000 MCMC repetitions. For each K value, five iterations were performed. The other parameters were set at default values: a standard deviation of 0.05, prior F_ST mean of 0.01, and different values of F_ST for different subpopulations. The optimal number of hypothetical clusters was indicated by the Delta K method (Evanno et al. 2005). To identify the potential admixed individuals between B. carambolae and B. dorsalis, an individual-based genetic cluster was plotted. STRUCTURE analysis under different assumptions was also run to verify the consistency. The repeated analyses considering (1) uncorrelated allele frequency and (2) missing data as recessive homozygotes for the null alleles were set to avoid the shared descent of samples and to verify possible bias from null alleles, respectively.

Principle Coordinate Analysis (PCoA) was used to display genetic divergence among fruit fly populations in multidimensional space. This analysis was based on allele frequency data and performed on genetic distance using GENALEX v.6.5 (Peakall and Smouse 2012). The subprogram MOD3D in NTSYS-pc v.2.1 (Rohlf 2005) was subsequently used for plotting the first three principal coordinates.

To test genetic homogeneity in different hierarchical population structures, Analysis of Molecular Variance (AMOVA) in ARLEQUIN v.3.1.1 was used (Excoffier and Lischer 2010). After 1,000 permutations, populations were grouped according to nine criteria described in Table 5.

The correlation analysis between genetic and geographic distance was performed using the subprogram ISOLDE in the GENEPOP package (Rousset 2008).

Analyses of the genetic networks were performed using EDENetworks (Kivelä et al. 2015). This program is advantageous in that it can provide graphical representations of the structure and dynamics of a system of interaction between populations in multidimensional space, without a priori assumptions of the clustering of populations and some of the constraints (e.g., binary branching) compulsory in phylogenetic trees. Data types of ‘genotype matrix, diploid, sampling site based’ were used as input. The genetic distance metrics or network files were calculated using the F_ST-based distance of Reynolds (Reynolds et al. 1983); networks were constructed.

Networks consist of nodes (or vertices), corresponding to populations, connected by links (or edges), representing their relationships or interactions. Connectivity degree (or Degree) is the number of edges connected to a node summarizing how strongly a population associated with the other populations in the system and whether or not it is a source population. Betweenness-centrality (BC) determines the relative importance of a node within the network as an intermediary in the flow of information. Each network was weighted demonstrating genetic similarity associated with each link.

The network can be analyzed at various meaningful thresholds (thr). thr is the maximum distance considered as generating a connection in the network. One meaningful distance is the one corresponding to the percolation threshold (D_p), edges with weights below the threshold were removed from the weighted network, and only the most important links were retained. Above the D_p level, there is a giant component containing almost all the nodes in the networks while below the D_p level, the network is fragmented into small disconnected components and the system therefore loses its ability to transport information across the whole system. Therefore, scanning at different thresholds was performed to analyze possible sub-structured systems to observe the sequential forms of clusters (Kivelä et al. 2015). The D_p and thr levels were determined by automatic and manual thresholding options, respectively.

All eight microsatellite loci tested within 13 populations have different levels of polymorphism in terms of number of alleles (ranging from moderately polymorphic, at five (Bcar181) to highly polymorphic at 16 (Bcar9)) and allele size range, as presented in Table 2. At locus Bcar73, allele 113 appeared to be fixed in male individuals of the SY5 strain, as reported by Isasawin et al. (2014). Linkage disequilibrium (LD) and Hardy-Weinberg equilibrium tests (HWE) were therefore performed for seven microsatellite loci for all populations. No significant evidence of LD among all loci was detected. After Bonferroni correction (Rice 1989), 28 out of 91 comparisons of loci and populations significantly deviated from the HWE results.

Overall genetic variations detected in each population is summarized in Table 3. Bactrocera carambolae samples collected from Southeast Asia showed relatively higher genetic variation than the introduced population (PR) for all parameters (i.e., n_a, n_e, n_r, A_r, n_p, A_p, H_O and H_E). In addition, the values of genetic variation of three B. dorsalis populations, two unidentified populations, and the SY5 strain were in the same range as B. carambolae. The n_p values were observed in all populations, except for the NS, DK and PR populations, varying from 0.125 (DP, JK, KS and SY5) to 0.625 (RB) per locus. Likewise, rare alleles (allele frequency less than 0.05) were also detected in all populations (n_r: ranging from 0.375 to 1.875 per locus), except KS. The average H_E values varied from 0.185 (PR) to 0.668 (RB). However, a deficiency in the average H_O values was found in all populations. Inbreeding (F_IS) was detected in all populations, ranging from 0.095 to 0.628. The average for null alleles was 0.12, varying from low (0.01) to high frequency (0.25), which may contribute to the deficiency of heterozygosity observed in the given populations.

Genetic differentiation among 13 populations was measured by the fixation index (F_ST), as shown in Table 4. Among seven populations, the pairwise F_ST values range from 0.134 (between PK and JK) to 0.631 (between DK and PR). Genetic differentiation was relatively high between native and introduced populations, ranging from 0.444 to 0.631. The introduced population (PR) was most similar to population PK (F_ST = 0.444) and JK (F_ST = 0.448). Within native areas of B. carambolae, the pairwise F_ST values ranged from 0.134 (PK and JK) to 0.344 (PK and DK). On the other hand, the pairwise F_ST values among samples of B. dorsalis varied from 0.210 (RB and CM) to 0.357 (CM and KS). Without PR, the degree of genetic differentiation between B. carambolae and B. dorsalis ranged from 0.181 (DP and RB) to 0.524 (NS and KS). The SY5 strain was revealed to be genetically closest to JK (F_ST = 0.278). The degree of pairwise F_ST among pairs of the SY5 strain and others varied from 0.372 (SY5 and DP) to 0.507 (SY5 and PR) (Table 4).

STRUCTURE analysis demonstrated the proportion of co-ancestry (Q) distributed in hypothetical clusters (K) whereas PCoA illustrated the genetic divergence of fruit fly populations in multidimensional space, as shown in Figure 2. Among seven B. carambolae populations, the Delta K value (Evanno et al. 2005) was determined to be K equals two (K = 2) as the optimal number. At K = 2, genetic clusters were divided into two groups: native and introduced B. carambolae. Cluster 1 contained all native populations (NS (Q = 0.979), PK (Q = 0.969), DP (Q = 0.984), JK (Q = 0.953), DK (Q = 0.991), and NT (Q = 0.995)) whereas the introduce population (PR) was distinguished, forming its own cluster (Q = 0.988) (Figure 2A). This subdivision corresponded to the first axis (44% of total variation) of the principal coordinate.

Figure 2. 

Three-dimensional plot of Principal Coordinate Analysis (PCoA) and STRUCTURE analysis. A the planes of the first three principal coordinates explain 43.65%, 20.13%, and 16.91% of total genetic variation, respectively, for seven B. carambolae populations using eight SSRs B the planes of the first three principal coordinates explain 33.05%, 23.17%, and 15.87%, respectively, for B. carambolae and B. dorsalis groups using eight SSRs C the planes of the first three principal coordinates explain 30.50%, 22.14%, and 18.53%, respectively, for the SY5 strain and wild populations using seven SSRs. Pie graphs, consisting of different colored sections, represent co-ancestor distribution of 185, 289, and 321 individuals in A two, B three, and C two hypothetical clusters, respectively.

When three additional populations of B. dorsalis (RB, CM and KS) and two unidentified populations (BD and WK) were included in the STRUCTURE analysis, the optimal number for K was three. Genetic clusters were separated into two groups: B. dorsalis belonged to cluster 1 while B. carambolae belonged to clusters 2 and 3 (Figure 2B). Six native populations of B. carambolae shared genetic memberships in cluster 2 (NS (Q = 0.976), PK (Q = 0.954), DP (Q = 0.519), JK (Q = 0.910), DK (Q = 0.970), and NT (Q = 0.758)). However, PR formed its own cluster, cluster 3 (Q = 0.986). Three B. dorsalis samples (RB (Q = 0.927), CM (Q = 0.955), and KS (Q = 0.952)) were grouped into the same genetic memberships (cluster 1) with two unidentified populations BD (Q = 0.927) and WK (Q = 0.953). Co-ancestor distribution between the B. carambolae and B. dorsalis clusters (Q ≥ 0.001) was observed in populations DP (Q = 0.457) and NT (Q = 0.236). Using PCoA, the first plane (33% of total variation) of multidimensional space also separated PR from the rest of the populations. Although the other two axes (23% and 16%, respectively) did not clearly divide samples into two groups in accordance with the STRUCTURE results, the group of B. dorsalis appeared to be plotted separately from the B. carambolae group.

The individual-admixture plot for K = 3 is presented in Figure 3. Individuals contained in the proportion of genetic cluster (Q) between 0.100 to 0.900 (0.100 ≤ Q ≤ 0.900) were identified as admixed individuals (or potential hybrids). Among the B. carambolae group, 18 of 185 individuals (9.73%) were admixed individuals. On the other hand, 13 of 104 individuals (12.5%) of the B. dorsalis group were admixed individuals. Pure individuals (Q > 0.900) were identified between the two groups. Eight individuals from the B. carambolae group were labelled as B. dorsalis whereas one individual from the B. dorsalis group was indicated to be B. carambolae.

Figure 3. 

The individual admixture plot for K = 3. Each bar reveals a single individual. Each color of bars represents each genetic cluster. Samples of B. carambolae belong to clusters 2 and 3 (green and blue, respectively) while samples of B. dorsalis belong to cluster 1 (red). Potential hybrids have a proportion of genetic cluster (Q) between 0.100 to 0.900 (0.100 ≤ Q ≤ 0.900) as identified with asterisk (*).

Adding the SY5 strain to the genetic cluster analysis, two was the optimal number for K. At K = 2, genetic clusters formed two groups likely corresponding to B. carambolae and B. dorsalis groups. Six populations of B. carambolae (NS (Q = 0.981), PK (Q = 0.913), JK (Q = 0.948), DK (Q = 0.964), and PR (Q = 0.993) belonged to cluster 1 (Figure 2C). The SY5 strain was also a member of the B. carambolae cluster (Q = 0.988). Cluster 2 consisted of the B. dorsalis group (CM (Q = 0.909), RB (Q = 0.960), KS (Q = 0.964), BD (Q = 0.918), and WK (Q = 0.933)). However, populations DP (Q = 0.596) and NT (Q = 0.677) also shared membership in this cluster and became an admixture structure. The principal coordinates did not clearly delineate the two clusters using the three-dimensional plot. The PR population was still isolated from the others while the B. dorsalis group seemed to form a peripheral group.

Analysis of Molecular Variance (AMOVA) was used to study the hierarchical structure of populations for different scenarios (Table 5). Overall, up to 23% of variation was attributed to the differences among groups whereas approximately 52% to 62% of variation was attributable to differences within populations. In scenarios 1 and 2, populations of B. carambolae were divided into subpopulations according to micro- and macro- geographies. At the micro-geographical level, seven populations were divided into five groups: Sumatra, Indonesia (NS and PK); Java, Indonesia (DP and JK); Malaysia (DK); Thailand (NT); and Suriname (PR). Non- significant differences among those groups were detected (scenario 1: P = 0.1877). Likewise, when all seven populations were divided based on macro-geography (Southeast Asia vs. South America), a non-significant difference between groups was still observed (scenario 2: P = 0.1953). A test of genetic homogeneity between B. carambolae and B. dorsalis is presented in scenario 3, illustrating a significant difference (P = 0.0479). Even though BD and WK were included in the B. dorsalis group, a significant difference was still detected (scenario 4: P = 0.0156). In scenarios 5 to 7, the SY5 strain was included to compare the genetic variation among other samples. The results revealed a non-significant difference between B. carambolae and the SY5 strain (scenario 5: P = 0.2483) and among B. carambolae, B. dorsalis and the SY5 strain (scenario 6: P = 0.0694). However, a significant difference was indicated when BD and WK were included in the B. dorsalis group (scenario 7: P = 0.0342). The last two scenarios was set following the STRUCTURE results (Figures 2B and 2C, respectively). Significant differences were detected in both scenarios (P = 0.0010 and P = 0.0039, respectively).

The correlation between genetic and geographic distance was analyzed using only wild samples consisting of B. carambolae and B. dorsalis populations. The correlation between genetic and geographical distance became non-significant [R² = 0.394, P = 0.106, F_ST/(1-F_ST) = 0.146 Ln (geographical distance) - 0.572] when B. carambolae samples were analysed. This fact indicates that there is no limitation of gene flow among B. carambolae across the region. Among B. carambolae and B. dorsalis populations, analysis showed significant correlation between genetic and geographical distances [R² = 0.449, P = 0.002, F_ST/(1-F_ST) = 0.180 Ln (geographical distance) + (0.868)], even though only the PR sample was excluded (R² = 0.119, P = 0.021).

Simplified networks were constructed for three different scenarios: (1) among seven B. carambolae populations, (2) among 12 populations belonging to B. carambolae and B. dorsalis clusters, and (3) among 13 populations, including the SY5 strain (Figures 4 to 6, respectively). Genetic distance metrics for the first and second data sets were estimated using eight microsatellite loci, but the third data set was analyzed using seven loci because the omitted locus is Y-pseudo linked in the SY5 strain (Isasawin et al. 2014). The networks were scanned for decreasing thresholds from the fully connected network to the percolation threshold (D_p) and lower threshold chosen (thr) to reveal the sequential substructure at decreasing thresholds. The first scenario was constructed based on data among seven B. carambolae populations (Figure 4). The percolation threshold D_p = 0.52 showed the emergence of a connection between native and introduced populations (Figure 4B). The node of JK had a larger degree (Degree = 6.0) and betweenness-centrality (BC = 5.0) than others and played a role connecting between native and introduced populations (Suppl. material 1: Table 1). The scanning at decreasing thresholds illustrated sub-cluster of native and introduced populations (thr = 0.40; Figure 4C). Moreover, JK and PK populations were the first to be jointed in the network (thr = 0.15; Figure 4D).

Figure 4. 

Simplified network of seven Bactrocera carambolae populations, and the sequential forms of cluster. The network was constructed using eight SSRs. Scanning was done for decreasing thresholds A is the fully connected network B is the percolation threshold (D_p = 0.52, with all links corresponding to distances superior to D_p excluded). JK plays an important role connecting between native and introduced populations C–D are the lower thresholds chosen (thr = 0.40 and 0.15, respectively) to reveal sub-structured network.

Figure 5. 

Simplified network of Bactrocera carambolae and B. dorsalis groups, and the sequential disconnection of the network. The network was constructed using eight SSRs. Scanning was done for decreasing thresholds A is the fully connected network B is the percolation threshold (D_p = 0.20, with all links corresponding to distances superior to D_p excluded). DP, JK, and NT are connecting between B. carambolae and B. dorsalis groups. Red dashed lines with number are corresponded to the threshold values, revealing serial disconnection of the network C is the lowest threshold (thr = 0.15).

Figure 6. 

Simplified network of the SY5 strain and wild populations, and the sequential disconnection of the network. The network was constructed using seven SSRs. Scanning was done for decreasing thresholds A is the fully connected network B is the percolation threshold (D_p = 0.23, with all links corresponding to distances superior to D_p excluded). DP, JK, and NT are connecting between B. carambolae and B. dorsalis groups C is the lowest threshold (thr = 0.15). Red dashed lines with number are corresponded to the threshold values, revealing serial disconnection of the network.

For the second scenario, five more populations belonging to the B. dorsalis cluster were included in the network analyses (Figure 5). At the percolation threshold D_p = 0.20, non-substructured system was observed (Figure 5B). Genetic connections between the B. carambolae and B. dorsalis groups were identified through the nodes of DP (Degree = 4; BC = 11.00), JK (Degree = 3; BC = 5.00) and NT (Degree = 2; BC = 5.00) (Suppl. material 2: Table 2). The network scanned below D_p demonstrated a serial disconnection of the system (Figure 5B). The first relationship of the network was detected between JK and PK populations (thr = 0.15; Figure 5C). To the contrary, the genetic links between native (JK and PK) and introduced (PR) populations of B. carambolae were revealed when increasing the threshold to thr = 0.55 (Degree = 11 and BC = 4.61) (Suppl. material 2: Table 2).

For the third scenario, the SY5 strain was added into the network to evaluate the genetic similarity between the SY5 strain and wild populations (Figure 6). The links between the SY5 strain and wild populations were detected above the percolation threshold (thr = 0.34 to 0.75) (Figure 6A). The SY5 strain shared greatest genetic similarity to JK (thr = 0.34). At the percolation threshold D_p = 0.23, non-substructured network was recognized (Figure 6B). Scanning below D_p showed disconnection of the system (Suppl. material 3: Table 3). Again, the first relationship of the network was detected between JK and PK populations (thr = 0.15; Figure 6C).

At the macro-geographical level, comparing among seven populations of B. carambolae from Southeast Asia and South America, the populations across the species’ native range possessed higher genetic variation than the introduced population, which is generally expected for invasive species. The first genetic connections between native and introduced populations were identified as Jakarta (JK) and Pekanbaru (PK), Indonesia. However, the genetic structure of the Suriname population (PR) (based on F_ST, STRUCTURE, PCoA, and genetic network analysis) was differentiated from the Carambola fruit fly in Southeast Asia. These results were congruent with multilocus phylogenetic analysis established by Boykin et al. (2014). They deduced that factors and processes shaping the genetic variation and population structure of B. carambolae in the introduced area potentially include genetic drift during the colonization process and local adaptation. Based on the genetic data in this study, JK and PK are possible sources of the PR population. Bactrocera carambolae accidentally invaded South America, likely by tourists and trade from Indonesia to Suriname. Even though this species was detected in 1975, it took up to 21 years (1975–1996) to establish its new populations in other areas, expanding 500 km eastward and westward from the original point of introduction in Suriname (Malavasi et al. 2000, 2013). In the meantime, according to the report of van Sauers-Muller (2005), host plants of B. carambolae (e.g., guava, Malay apple, carambola, West Indian cherry, and mango) in Suriname were limited to backyards; agricultural production had not yet developed. Hosts for the fly were occasionally recorded to be different between native and introduced areas (van Sauers-Muller 2005), as shown in Table 6. Regarding that evidence, the conditions of habitats, including sufficient time for genetic drift to take effect, may be natural factors causing high genetic differentiation between native and introduced populations. Likewise, the same factors, including limitations of human activity, among four countries near Suriname may slow down species distribution in South America. This pattern is similar to the case of the related species B. dorsalis in several introduced areas such as Hawaii, Myanmar, and Bangladesh, in that their genetic variation and population structure were shaped by genetic drift and different local adaptations (Aketarawong et al. 2007). However, the current study included only one population from South America. Research using more samples, if available, from these regions should provide a better understanding of the demographic dynamics of the Carambola fruit fly within South America as well as between the two continents.

Within Southeast Asia, B. carambolae demonstrates high genetic variation and homogeneous population structure. West Java, in particular JK, is also a potential source of B. carambolae populations in Southeast Asia. JK showed relatively higher genetic variation and greater values of Degree and betweenness-centrality than other populations in the genetic network. Moreover, this area is important for the cultivation of star fruit and is a center for transportation to other cities and countries. The homogeneous genetic pattern of B. carambolae in native areas is similar to B. dorsalis collected from neighboring countries, including Thailand, Laos, and Cambodia (Aketarawong et al. 2007, 2014a, Schutze et al. 2012). Both not having limitations on human migration and the intensive cultivation of similar host plants could enhance gene flow, shaping genetic homogeneity among flies from nearby countries. Although the geography of Indonesia, Malaysia and Thailand is not entirely contiguous, increasing trade can promote the migration of insect pests within the country, as well as among other countries (Suputa et al. 2010). Therefore, effective quarantine measures are important to reduce the pest’s distribution in Southeast Asia.

Pairwise F_ST between native and introduced populations was significantly higher than zero. Approximately 13.4% to 32.9% and 44.4% to 63.1% of genetic diversity were results of genetic difference among populations within native areas and among populations between native and introduced areas, respectively. From the comparison of the genetic diversity of other closely related species using eight similar microsatellite loci with different primer sets as the current study, lower levels of genetic diversity (approximately 2% to 18%) were estimated from B. dorsalis (Aketarawong et al. 2014b). The most likely explanation of the situation involves a high level of gene flow and/or recent population divergence. In case of B. carambolae, it implies that the level of colonization of this invasive fly may not be as high as B. dorsalis. Alternatively, it can be deduced that populations, in particular between native and introduced areas, became diverted, congruent with studies of multilocus phylogeny using sequence analyses (Boykin et al. 2014). However, the F_ST value can be influenced by the geographical difference of sampling locations (Neigel 2002). Samples of B. dorsalis were collected at no more than 200 km intervals whereas in this study, B. carambolae populations were collected from locations varying from 20 (Depok) to 18,022 km (Paramaribo) from Jakarta, Indonesia. Research using more samples from different locations on a fine scale and/or more genetic markers may provide more details.

Departures from HWE were quite high and null alleles might be responsible for the departures. The average null allele frequency was moderate (0.12), varying from low (0.01) to high frequencies (0.25). The high departure from HWE was also observed in other study using a single set of microsatellite markers for different species such as the Ceratitis FAR complex (52.4%, Virgilio et al. 2013). We verified possible bias produced by the presence of null alleles in our data set. The STRUCTURE analysis was also repeated considering missing data as recessive homozygotes for the null alleles. We found that the effect of null alleles was negligible (Suppl. material 4: Figure 1).

We found 44 alleles shared between B. carambolae and B. dorsalis while 14 and 27 alleles were detected in only B. carambolae and B. dorsalis, respectively. Coincidence of similar/different allele profiles between them at microsatellite loci may be due to several phenomena including retention of ancestral alleles in both sister species; substantial gene flow between the two species; size homoplasy (Estoup et al. 2002). For the latter case, homoplasy is possibly ruled out. It does not represent a significant problem for many types of population genetics (i.e., only 1–2% underestimation of genetic differentiation) considering only when microsatellite with high mutation rate and large population size together with strong allele size constraints were involved. The choice of more realistic mutation models than a common strict-Stepwise Mutation Model (SMM) should alleviate the homoplasy effect (review in Estoup et al. 2002, review in Selkoe and Toonen 2006).

Using the genetic cluster approach, assuming correlated allele frequencies in different clusters were likely to be similar due to migration or shared ancestral, species’ genetic structures were determined. We found admixed individuals (potential hybrids) in both clusters with relatively similar ratio (9.73% in B. carambolae cluster and 12.5% in B. dorsalis cluster). To avoid the shared descent, a stricter model using uncorrelated allele frequencies was tested. The results still presented species’ genetic cluster and admixed individuals (Suppl. material 4: Figure 1). Genetic connectivity also revealed that populations DP, JK, and NT in the B. carambolae cluster and population WK in the B. dorsalis cluster consisted of several admixed individuals, playing the role of linker between B. carambolae and B. dorsalis in the genetic network. Comparing to B. dorsalis and the junior synonym, B. papayae, they share better genetic profiles (Schutze et al. 2012, Krosch et al. 2013, Aketarawong et al. 2014b). Weak and no genetic structure were presented using both similar (but different primer sets) (Aketarawong et al. 2014b) and different microsatellite loci (Krosch et al. 2013). It was observed that a single cluster dominated the ancestral of all samples, when uncorrected allele frequencies were assumed in the analysis (Krosch et al. 2013). Using a different assumption of allowing similar allele frequencies between populations, the population’s genetic structure was observed rather than species and admixed individuals were found among population clusters (Aketarawong et al. 2014b).

The current study therefore provided additional evidence to support an incomplete reproductive isolation between B. carambolae and B. dorsalis (Wee and Tan 2005b, Augustinos et al. 2014, Isasawin et al. 2014). Boundaries between B. carambolae and B. dorsalis may be semipermeable, varied as a function of genome region. Alleles at microsatellite loci used in this study could be introgressed between two species rather than other nuclear and mitochondrial DNA sequences (Boykin et al. 2014). To achieve a clearer picture of species boundary, genome-wide comparisons (ranging from modest number of microsatellite loci to full genome sequences, transcriptome, or SNPs analysis) between recently diverged forms or species should be performed. This may not only offer patterns of differentiation across the genome, but also define the dynamics of species boundary (Harrison and Larson 2014).

Isasawin et al. (2014) reported a proof of concept to characterize and use the new genetic sexing Salaya5 strain (SY5) for controlling B. carambolae. At that time, two wild populations of B. carambolae were then included in the experiment. However, this study was focusing on how it would be possible to use the SY5 strain for pest control on a broader level, not only locally. Therefore, more samples of B. carambolae from native and introduced areas (i.e., Indonesia, Malaysia, Thailand, and Suriname) were included. More analyses on genetic variation, population structure and genetic network were performed between the SY5 strain and wild populations. We found that the SY5 strain had genetic variation, population structure, and genetic similarity comparable to B. carambolae, rather than B. dorsalis, in Southeast Asia. The results strongly confirmed that the Salaya5 strain had not diverted away from its original genetic makeup. Under laboratory condition, at least 12 generations apart, selection did not strongly affect genetic compatibility between the strain and wild populations. Therefore, the SY5 strain could be included in the pest control programs using male-only SIT for control B. carambolae at local and regional levels. However, an actual mating test is still required between the strain and samples from introduced populations.

SIT is a species-specific control method that can deliver environmental benefit. However, it may be restricted where at least two major target pests coexist and need to be controlled. Releasing sterile males of only one target may not ensure a reduction of all problems (Alphey et al. 2010). Bactrocera carambolae and B. dorsalis were reported to be sympatric species in some areas (e.g., Kalimantan). Although several lines of evidence suggested that both species showed mating compatibility to some degree (McInnis et al. 1999, Schutze et al. 2013, Isasawin et al. 2014), release of either sterile genetic sexing Salaya1 or Salaya5 strain may be not enough to control both target pests. The release both of the sterile genetic sexing Salaya5 and of Salaya1 strains for controlling B. carambolae and B. dorsalis, respectively, should maximize the success of SIT programs. A simulation of mating competitiveness tests and field operation, when releasing two species together, will be further required to gain knowledge for confirmation of the proposed idea.

In order to identify the fruit fly samples, although microsatellite data showed significantly different population structure of the two species, eight of 185 individuals (4.32%) and one of 104 individuals (0.96%) belonging to the B. carambolae and B. dorsalis clusters were identified as opposite to their original assumed identity. At the individual level, microsatellite data in this study may not provide definitive data for studying systematic questions of incipient species. However, at the population level, microsatellite data can be used to distinguish species. This is similar to the case of the Ceratitis FAR complex in that genetic clustering can solve three species’ statuses whereas other data (i.e., morphology, phylogenetics based on DNA sequence analyses, and niche) cannot (Virgilio et al. 2013). In this study, the two unidentified populations should be good examples to support this particular advantage of using microsatellite markers. Therefore, a combination of microsatellite data with other approaches should be better than a stand-alone approach to define and confirm individual and/or population.