Highly variable COI haplotype diversity between three species of invasive pest fruit fly reflects remarkably incongruent demographic histories

Distance decay principles predict that species with larger geographic ranges would have greater intraspecific genetic diversity than more restricted species. However, invasive pest species may not follow this prediction, with confounding implications for tracking phenomena including original ranges, invasion pathways and source populations. We sequenced an 815 base-pair section of the COI gene for 441 specimens of Bactrocera correcta, 214 B. zonata and 372 Zeugodacus cucurbitae; three invasive pest fruit fly species with overlapping hostplants. For each species, we explored how many individuals would need to be included in a study to sample the majority of their haplotype diversity. We also tested for phylogeographic signal and used demographic estimators as a proxy for invasion potency. We find contrasting patterns of haplotype diversity amongst the species, where B. zonata has the highest diversity but most haplotypes were represented by singletons; B. correcta has ~7 dominant haplotypes more evenly distributed; Z. cucurbitae has a single dominant haplotype with closely related singletons in a ‘star-shape’ surrounding it. We discuss how these differing patterns relate to their invasion histories. None of the species showed meaningful phylogeographic patterns, possibly due to gene-flow between areas across their distributions, obscuring or eliminating substructure.

A fundamental assumption of species-level research is that sampling specimens from across their range accurately represents intraspecific genetic diversity [1][2][3][4] . This tenet is based on population-level processes such as genetic drift and selection, which follow the distance decay principle 5 resulting in regional differences that increase with geographic distance as the probability of gene-flow between more distant populations decreases. These divergences may be augmented in heterogenous landscapes or under strong environmental gradients 6 . Indeed, studies using various, predominantly native, taxa over large geographic ranges typically find that species with broader distributions, on average, have more intraspecific diversity than species with more restricted distributions (e.g. 3,4,7 ). However, depending on the life-history characters and evolutionary history, intraspecific genetic diversity can be regionally heterogeneous and may vary widely between taxa (e.g. 6,8 ).
We tested distance decay assumptions with three of the worst agricultural pests in Southeast Asia, since this has implications for their control, quarantine and the reconstruction of invasion patterns based on demographic history. When intraspecific variation is correlated with geography it becomes possible to trace the geographic origins of invasive species and provides a framework to prevent future incursions and better manage invasive populations [9][10][11] . Moreover, populations can be adapted to their local environments, which influences the chance of a successful future invasion 12 or can even reinvigorate invasions by enabling access to new hosts or through biocontrol resistance 13,14 . If, on the other hand, no intraspecific variation is detected in the invasive range, it suggests a single small initial invasive population or multiple invasions from the same source. In this situation,

Material & Methods
Taxon sampling. Based on morphology and a seven-gene molecular phylogeny with 178 Dacini species 40 , B. correcta and B. zonata are sister species, and Z. cucurbitae is closely related to Z. tau Walker, Z. synnephes Hendel and Z. choristus May. Zeugodacus cucurbitae is widespread and is documented to be expanding its distribution, which currently includes much of subtropical Africa, South-East Asia and Pacific Islands, including Guam and Hawaii 24,37 (Fig. 1C). Bactrocera zonata is present in northern Africa and Saudi Arabia, the Indian subcontinent and has been reported from Thailand and Vietnam (Fig. 1A). Partly overlapping with this distribution, B. correcta can be found throughout most of mainland South-East Asia (Fig. 1B) 41 . Samples were collected between 2005 and 2017 ( Fig. 1, BOLD, https://doi.org/10.5883/DS-ZOCOCU). Bactrocera correcta and B. zonata were collected with the male attractant methyl eugenol, and Z. cucurbitae with cue lure 42 . Bucket traps with one of these lures and a 1 × 1 cm dichlorvos strip as a killing agent were maintained for three to five days along trails in natural areas or near agricultural land. Specimens were preserved in 95% ethanol and stored at −20 °C. We selected 441 B. correcta, 214 B. zonata and 372 specimens of Z. cucurbitae for DNA extraction. All identifications were done based on morphology.
DNA extraction, amplification, and sequencing. One to three legs of each specimen were typically used for total genomic DNA extraction, voucher specimens are deposited in the University of Hawaii Insect Museum (UHIM). Some specimens were additionally used for genomic studies (e.g. 43 ), such flies were fully ground up to obtain higher concentration DNA extracts. DNA was extracted using the Qiagen (Valencia, CA) DNeasy animal blood and tissue kit, following manufacturer's protocols. A 780 base-pair (bp) section of the cytochrome c oxidase subunit I 3′P region (COI-3P) of the mitochondrial DNA was amplified using the forward primer HCO-2198rc (5′-3′ GCT CAA CAA ATC ATA AAG ATA TTG G) 44 and reverse primer Pat-k508 (5′-3′ TCC AAT GCA CTA ATC TGC CAT ATT A) 45 . The PCR thermal conditions were 2 min. at 94 °C, 40 cycles of (94 °C for 30 sec., 53 °C for 30 sec. and 70 °C for 60 sec.) with a final extension at 70 °C for 10 min. Bidirectional Sanger sequencing was outsourced to Eurofins (Louisville, Kentucky, USA). We combined chromatograms of the forward and reverse reads in Geneious R10, trimmed primer regions, and created alignments for each species. The alignments contained no insertions or deletions, and all sequences were checked for stop-codons, which can indicate the presence of nuclear pseudogenes. The sequences are deposited both in BOLD (https://doi. org/10.5883/DS-ZOCOCU) and NCBI Genbank (accession numbers MT257264 -MT258375).
B. correcta, 756 bp aligned, 214 B. zonata, 752 bp aligned and 372 Z. cucurbitae, 764 bp aligned (https://doi. org/10.5883/DS-ZOCOCU). Some papers have studied the COI diversity of B. correcta regionally 46,47 , but because of only partially matching COI target regions (i.e. COI-5P versus COI-3P) we have decided to not include their data in our study. Also, Kunprom et al. 46 reported a strongly divergent COI lineage within B. correcta, matching earlier results based on the COI-3P region 47 . They attributed this anomaly to Wolbachia bacterial infection, but upon comparison of the published COI-3P sequences with our COI reference dataset (data not shown) we conclude this lineage to match a cluster that includes B. nigrotibialis, B. nigrifacia and B. nigrofemoralis, which share some of the same hosts as B. correcta and can be morphologically similar 36,41 . For maximum likelihood tree inference of the relationships of the target species and their closest relatives (based on 40 ), we used unique haplotypes only and added B. nigrotibialis (Perkins), Zeugodacus tau (Walker), Z. synnephes (Hendel) and Z. choristus (May) to the alignments. Maximum likelihood tree inference was done using RaxML v8 with 50 searches for the best tree and subsequent multiparametric bootstrapping with the extended majority rule parameter to estimate statistical support values. We calculated haplotype networks using the TCS statistical parsimony algorithm 48 , employed in the software "TCS" v.1.21 49 . The raw output of TCS was visualized in the web implementation of tcsBU 50 and optimized for publication in Adobe Illustrator CC. Haplotype rarefaction curves were estimated in R using the HaploAccum() function from the package SpideR with 1,000 permutations, and we used the chaoHaplo() function from the same package to obtain the non-parametric Chao 1 estimator of total haplotype diversity. The "tools and recipes for evolutionary genetics and genomics" (EggLib 51 ); python library was used to calculate the descriptive statistics for nucleotide diversity π, Tajima's D, Watterson's estimator of θ, and Fu's Fs.

COI Phylogeny.
We successfully obtained COI sequences of 428 Bactrocera correcta, 214 B. zonata and 372 Zeugodacus cucurbitae. The COI maximum likelihood trees (Fig. 2) based on unique haplotypes, with added outgroup sister taxa, show that each of the study species can reliably be separated from its closest sister species using a monophyly criterion, and that the intraspecific variation does not overlap with the variation between Green areas indicate countries for which the species has been reported and we have sampled, brown areas indicate countries for which the species has been reported but we were unable to sample. Red circles indicate sampling localities. The map insert in C shows the Hawaiian Islands. Note that the actual distributions of the species do not always follow political boundaries, the flies do not occur in northern China, or the northern Islands of Japan.
species. All haplotypes are unique to each respective species. The intraspecific variation ranges from 1.32% in Z. cucurbitae to 3.06% in B. zonata (Table 1), and this variation is distributed across the COI region (Supplemental Information Table 1).
Haplotype networks. We find contrasting patterns of haplotype diversity between the three species (Fig. 3), where B. zonata has the highest diversity and most haplotypes are represented by singletons, B. correcta has ~7 dominantly represented haplotypes and Z. cucurbitae has a single dominant haplotype with closely related singeltons in a 'star-shape' surrounding it. None of the species displays a biogeographic signal, where certain haplotypes or parts of the haplotype network can be correlated with a biogeographic region. For B. zonata, many of the haplotypes are unique to a region, but also unique in the network and further sampling will likely show that these are shared between regions.
Haplotype rarefaction. The haplotype rarefaction curves with randomly sampled diversity do not reach an asymptote for any of the species (Fig. 4), indicating that we have only sampled part of the actual diversity. Despite their strongly different haplotype networks (Fig. 3), Bactrocera correcta and Zeugodacus cucurbitae have largely overlapping haplotype rarefaction curves and confidence intervals, although the Chao 1 estimated total diversity for the former is 1,110 haplotypes and only 270 for the latter. The angle on the curve for Bactrocera zonata is close to a 1:1 ratio where every new sample represents a new haplotype, indicating a much higher diversity for this species, which is underlined by a Chao 1 estimated total haplotype diversity of 1,890.
Population genetic statistics. From the lack of geographic partitioning in the haplotype networks (Fig. 3), we assumed the different sampling points for each species to be part of a panmictic population and estimated descriptors of their genetics ( Table 1). The number of observed haplotypes ranged from 70 to 184. With Chao 1 total estimated haplotypes per species ranging from 270 to 1,890 (Fig. 4.), we sampled 7.2% (in B. correcta) -25.9% (in Z. cucurbitae) of the total estimated diversity. Watterson's estimator of θ and nucleotide diversity π estimators show that B. zonata is the most diverse species, followed by B. correcta, and Z. cucurbitae shows little diversity, matching the results from the haplotype network (Fig. 3). Wattersons's estimator of θ, calculated through the formula θ w = 4N e µ can also be used as an estimate of effective population size (N e ), and shows that Z. cucurbitae has the smallest effective population size, followed by B. correcta and B. zonata. Negative values of Fu's Fs are evidence for an excess number of alleles relative to the expected. Fu's Fs is negative for all three species, but the largest negative is in B. zonata. An excess of low frequency polymorphisms relative to expectation is expressed by negative values of Tajima's D, and are commonly interpreted as demographic estimators for population expansion or selection. We found negative values of Tajima's D for all three species, but most strongly in Z. cucurbitae, followed by B. zonata and B. correcta.

Discussion
Distance decay COI diversity. We examined the role of the distance-decay model in the mitochondrial gene COI as a tool to understand the degree of spatially heterogeneity in genetic diversity and whether important pest species have experienced recent range expansions. In general, intraspecific genetic diversity is a result of demographic history, mutation and selection 52 . The different genetic patterns observed can reveal much about the forces, such as rapid expansion or significant bottlenecks, directly affecting species of interest. Considering the prevalence of horticulture in human history, it is difficult to assess what should be considered the native range of the three species. The native range can sometimes be assessed a posteriori, because it will commonly have the highest genetic diversity (e.g. 53 ), but a priori assumptions rely on historical records of observations and these may  be inaccurate due to, e.g., faulty taxonomy 54 . For our study, we therefore attempted to include as many samples from throughout the current ranges of the species as possible, and we did not detect variation in haplotype diversity between regions that would enable the distinction between native and non-native ranges. Despite the fact that we did not detect regional variation, all three species in our study show markedly different COI diversity patterns. This might be expected to influence their invasion success if genetic diversity is a boon for invasive species. The differences may be a result of the differing topology of the landscape of their range, ecological interactions and/ or species-specific variation in the mutation rate (μ) of COI. Zeugodacus cucurbitae has the widest distribution in our study, yet, conversely, it shows a lack of haplotype diversity. This could be due to very recent expansion in Z. cucurbitae from a very genetically limited source population that has adapted to the agricultural environments and since spread across the globe. The pattern of a single adapted population spreading rapidly due to humans has also been observed in other organisms (e.g. 13,55,56 ). Regardless of the underlying reasons for the contrasting patterns found in this study, the diverging haplotype patterns suggest that distance decay based principles likely do not apply to these pest species, at least not anymore. As a consequence, the management strategies undertaken for each pest should be based on the absence or presence of intraspecific diversity in regional populations observed. However, it should be kept in mind that different markers may reveal different intraspecific patterns, depending on their resolution and evolutionary history 10 . These results also pertain to species of conservation concern, since understanding current and past population retractions as measured by genetic diversity across their range has important implications for the evaluation and identification of declining taxa.  58,59 . In stark contrast with these results, we found a new haplotype for almost every additional sample of B. zonata that we included. The demographic estimators for the latter do not suggest that there has been rapid population expansion, at least not from a small source population, and in the regions covered by our sampling. A study of the COI variation of 1,600 specimens of two other prominent pests of the tribe, B. dorsalis and B. carambolae, also showed high levels of haplotype variation and a large global distribution with multiple recently invaded continents 60 . Overall, the results of Tajima's D estimator seem congruent with the historical record of invasions. Zeugodacus cucurbitae has the largest negative value for Tajima's D and is the most invasive and widespread of the three pests, B. zonata is the second and B. correcta third. However, it is worth noting that the melon, gourd and squash hosts of Z. cucurbitae are also the most widely grown and transported. The demographic estimators can reveal if there have been recent changes in population sizes, but the values depend on the original population structure in the 'native' range and only partly reflect the species' pest 'potential' . The pest potential, which can be regarded as the combined factors of the likelihood of reaching new areas and the possibility of becoming established 24 , is likely mostly determined extrinsically, by human actions -such as the growing of host plants and transporting of fruit. It appears that depending on these circumstances, any of the species in the tribe can develop into pests, with serious implications for quarantine which currently focuses on the handful of already known pest species. Undetected(?) biogeographic patterns. Although COI has been shown to reveal geographic patterning in some species (e.g. 61-63 ), we did not detect any biogeographic structure in our three pest species. The effects of historical climatic change were not as severe in Southeast Asia as in temperate regions, with more continuous possibilities for gene-flow throughout the region which may have resulted in less strongly pronounced biogeographic patterns 64,65 . However, a high resolution (1,097 SNPs) study of the melon fly Z. cucurbitae revealed geographic structuring with six to ten clusters that were only partly detected using microsatellite data 43 , and are not detected in our current study using COI. Despite being a widespread species occurring from Africa across Asia into the Pacific, we only found a single dominant haplotype. Our sampling is spread over the full range of Z. cucurbitae, including different invasive populations in Africa. The results from the demographic estimators and the star-shape of the haplotype network therefore all suggest that Z. cucurbitae was restricted to a small geographic area and began to expand rapidly relatively recently-an expansion that is still ongoing as more of Africa is invaded 37 . There are, however, some alternative explanations that we cannot exclude based on our present data, such as the presence of Wolbachia bacteria that could have influenced the mitochondrial diversity 66 . Nonetheless, the SNP data supports biogeographic patterning. The differences between mtDNA versus nuclear DNA are partly due to the differences in mutation rates and inheritance, but also to the fact that the method to derive nuclear DNA SNP's results in a much larger data set of independent sites with a greater ability to resolve small differences between populations. As more SNP data becomes available for more species, it will be interesting to test at which time-scales, or under which evolutionary circumstances, the resolution between mtDNA and SNP data are similar.
In contrast to the data available for Z. cucurbitae, the biogeographic structures of B. zonata and B. correcta have not been studied with high-resolution (i.e. genomic) methods, which may similarly reveal patterning that we did not detect using COI. The peach fruit fly B. zonata is mostly found on the Indian subcontinent and is expanding www.nature.com/scientificreports www.nature.com/scientificreports/ westward into Africa. Of the three species included in our study, our sampling of B. zonata was most limited, with high density sampling in Bangladesh but lacking specimens from the distribution west of Bangladesh. Nonetheless, it had the highest haplotype diversity. Despite the peach reference in its common name, it was likely first introduced into Africa with guava in 1993 25 . The area sampled by us likely is part of the original native range. The hosts of B. zonata and B. correcta largely overlap and similar biogeographical and population patterns might have been expected. Based on their current distribution, overlapping host range, phylogenetic relation and similar morphology, the most parsimonious conclusion is that B. correcta and B. zonata speciated through vicariance. Bactrocera correcta is now found mostly in eastern Asia, and B. zonata is mostly found in the Indian subcontinent with recent westward expansion 37 . The two species only overlap in Sri Lanka, India, Nepal, Bangladesh and Thailand. These vicariant distributions, despite the likely possibilities for human-mediated dispersal (fruit trade) into either range, may indicate strong competition or selection pressures that limit their avenues for invading regions already occupied by the other species. Fragments of the COI barcode region (COI-5P) of B. zonata have been studied in India 67 , where they also concluded that B. zonata has high intraspecific variability, but showed no biogeographic patterning. The third species in our study, the guava fruit fly Bactrocera correcta is mostly found in mainland Southeast Asia. Our sampling covers a large part of the distribution, only lacking Malaysia and India, and we likely included most of its native distribution and COI haplotype diversity. A regional study of the COI-5P diversity of B. correcta in Thailand concluded that there was much genetic homogeneity, likely because host plants are spatially and temporally continuously available 46 . However, we found samples from Thailand to belong to at least nine different haplotypes, mostly those dominantly represented in our dataset. All three species pose large invasion risks and biogeographical analysis of intercepted specimens can help to diagnose pathways and block them (e.g. 60,68 ). However, from the lack of biogeographic patterning we conclude that COI data cannot be used to trace the origins of introductions or intercepted flies of Z. cucurbitae, B. correcta and B. zonata.

COI-based identification.
Percentage-wise, the ranges for intra-and interspecific COI variation we find in our study species are similar to those found in other groups of life, such as moths, beetles, or birds, which commonly show <3% intraspecific and >2% interspecific variation 4,16,20,69,70 . Even though the estimated total number of haplotypes between the three species varies seven-fold, the maximum intraspecific diversity is at most ~3% and all species can be identified reliably using COI. The high mutation rates in COI, although predominantly in the third-codon position and synonyms are generally known across all animal taxa 71 , and have made the marker a useful part of a broader dataset for species delimitation and recognition, and to some extent population-level studies. Contrasting patterns of COI evolution between higher groups has been observed before 71,72 , which has been linked to the functional properties of the gene, but here we show that large differences can exist between closely related (sister) species as well, or species within the same tribe. Although we only studied three species, the patterns of haplotype diversity are starkly contrasting with very low haplotype diversity in Z. cucurbitae to very high in B. zonata. The estimates of the proportion of the diversity that we sampled range from 7.2-25.9%. All three of these species are pests associated with crops and would be expected to have high effective population sizes (N e ), given the potentially high population density in agricultural systems. The statistic Watterson's θ is dependent on N e and genetic diversity µ, and it is possible that only µ varies significantly between these three pests. Mitochondrial genetic diversity, in contrast to nuclear genetic diversity, may not be linked to N e but only to µ and selection 52,73 . In Bactrocera dorsalis and B. carambolae, high mutation rates were postulated to cause hyper-diversity in their mitochondrial genomes 60 . However, another explanation for this pattern may have to do with the demographic and invasion history of each species. Despite the common application of COI sequence data for large groups of life across the planet 2,19,20,74 , there are many open-ended questions on how this molecule functions, and to what extent there is active selection 71,75 . There have been few studies thus far on the intraspecific variation of COI using large sample sizes, but as DNA sequencing continues to become cheaper (e.g. 76 ) there will hopefully be larger sample sizes for more groups and we can evaluate to what extent general rules apply and how many samples are required to cover the haplotype diversity, to ensure reliable species identification.

conclusion
Our results indicate that COI was effective for species level identification across the ranges of these three invasive pest species, even in the presence of considerable intraspecific heterogeneity among species. The intraspecific heterogeneity compared between two closely and one more distantly related species with a non-overlapping host range suggests strongly different demographic histories between all three. Generalized distance decay predictions may not be applicable for pest species, because recently expanded distributions due to globalization may obscure historic biogeographic patterns. As a result, we could not use COI intraspecific variation to identify source populations or invasion pathways to augment pest management. Future studies based on high-resolution genomic data are likely to uncover additional patterns, which may be suitable for identifying invasion pathways. We found that simple demographic estimators based on COI haplotype diversity can be used to infer the historical expansion speed and complement otherwise documented distribution data, but that it cannot predict the virulence of a pest because its spread is more likely dependent on extrinsic factors, such as globalization.

Data availability
All data used in this paper is available through BOLD, https://doi.org/10.5883/DS-ZOCOCU as well as Genbank accession numbers MT257264 -MT258375.