Population structure of a global agricultural invasive pest, Bactrocera dorsalis (Diptera: Tephritidae)

Abstract Bactrocera dorsalis, the Oriental fruit fly, is one of the world's most destructive agricultural insect pests and a major impediment to international fresh commodity trade. The genetic structuring of the species across its entire geographic range has never been undertaken, because under a former taxonomy B. dorsalis was divided into four distinct taxonomic entities, each with their own, largely non‐overlapping, distributions. Based on the extensive sampling of six a priori groups from 63 locations, genetic and geometric morphometric datasets were generated to detect macrogeographic population structure, and to determine prior and current invasion pathways of this species. Weak population structure and high genetic diversity were detected among Asian populations. Invasive populations in Africa and Hawaii are inferred to be the result of separate, single invasions from South Asia, while South Asia is also the likely source of other Asian populations. The current northward invasion of B. dorsalis into Central China is the result of multiple, repeated dispersal events, most likely related to fruit trade. Results are discussed in the context of global quarantine, trade, and management of this pest. The recent expansion of the fly into temperate China, with very few associated genetic changes, clearly demonstrates the threat posed by this pest to ecologically similar areas in Europe and North America.


| INTRODUC TI ON
The horticultural sector is one of the largest within the global agricultural economy. Global fruit and vegetable production was estimated at 676.9 and 879.2 million tonnes, respectively, in 2013, with global fruit exports in the same year valued at USD 97.02 billion (Gyan Research and Analytics, 2014). While developed nations such as the United States have long-established horticultural sectors, significant growth in commercial horticulture is also occurring in the developing world. An expanded horticulture sector is seen not only as a mechanism to increase human health in developing nations, but also as a way of increasing general living standards through cash generation from the sale of fresh commodities for export (Virchow & Jaenicke, 2016;Weinberger & Lumpkin, 2007).
One of the major biological impediments to horticultural production and export in tropical and subtropical regions of the world is the frugivorous tephritid fruit flies (Diptera: Tephritidae) (Hendrichs, Vera, De Meyer, & Clarke, 2015). The tephritid or "true" fruit flies (not to be confused with drosophilid fruit flies) lay their eggs into sound, near-ripe fruit on plant, where the resultant larvae feed.
Depending on commodity, and in the absence of controls, fruit fly damage can easily lead to 80% to 100% crop loss (White & Elson-Harris, 1992). The global fruit fly problem is exacerbated by a small group of highly polyphagous, highly invasive pest species which competitively dominate local fauna if they enter and establish in a region (Duyck, David, & Quilici, 2004;Duyck et al., 2006), and which can subsequently stop fresh commodity trade because of the quarantine risk they pose (Dohino et al., 2016). The two best known of these invasive tephritids are the Mediterranean fruit fly, Ceratitis capitata (Wiedemann), and the focus of this paper, the Oriental fruit fly, Bactrocera dorsalis (Hendel).
Oriental fruit fly is one of the world's most invasive and polyphagous pests of agriculture, with a recorded host range of over 250 fruits and vegetables (Clarke et al., 2005). Endemic to the Indo-Asian region, the fly first established outside this native range in Hawaii in 1945, where it remains a major pest (Vargas, Piñero, & Leblanc, 2015). The fly has subsequently invaded the continental United States on numerous occasions and, while the formal regulatory position is that it is currently absent from the continental United States, debate exists in the scientific literature as to whether it is permanently established in California (Papadopoulos, Plant, & Carey, 2013), or is a repeat invader (Barr et al., 2014). Regardless of the position in the United States, B. dorsalis is invasive and permanently established in several South Pacific countries (Vargas et al., 2015), has invaded and been eradicated twice in Australia (Cantrell, Chadwick, & Cahill, 2002), is currently actively invading Central China (Chen, Zhang, Ji, Yang, & Zheng, 2014), and is an "A1" quarantine pest for the European Union (EPPO, 2015). However, it is its invasion, spread, and establishment in sub-Saharan Africa that has received most attention in recent time. The fly was first detected in Kenya in 2003 (Lux, Copeland, White, Manrakhan, & Billah, 2003), and within a span of 14 years has spread across all of sub-Saharan Africa and only small parts of South Africa remain free of the pest (Manrakhan, Venter, & Hattingh, 2015). The cost of lost export markets to Africa due to the invasion has been estimated at $2 billion (Ekesi, De Meyer, Mohamed, Massimiliano, & Borgemeister, 2016  , these latter three species are now recognized as junior synonyms of B. dorsalis (Drew & Romig, 2013;Schutze, Mahmood et al., 2015;Schutze, Aketarawong et al., 2015).
While the synonymization clarifies taxonomic identity and helps some aspects of pre-and post-harvest control and market access (Dohino et al., 2016;Hendrichs et al., 2015), they also create new challenges. An organism, whose geographic range extends from Africa, across Asia to the Pacific, might be predicted to exhibit macrogeographic population structuring (Ascunce et al., 2011;Gloria-Soria et al., 2016;Virgilio, Delatte, Backeljau, & De Meyer, 2010;Zhang, Edwards, Kang, & Fuller, 2014). As the International Plant Protection Convention (FAO, 2011) recognizes "Pest" as "any species, strain or biotype of plant, animal or pathogenic agent injurious to plants or plant products," synonymization of taxa does not negate the issue that different geographic populations may still show high levels of population structuring and so be of potential quarantine and trade concern at the "strain" level.
To address this issue, this paper presents the most comprehensive global assessment of B. dorsalis population structuring yet undertaken. We sampled from across the entire range of B. dorsalis occurrence, including invasive locations, and used morphological data (geometric morphometric analysis of wing shape) and molecular temperate China, with very few associated genetic changes, clearly demonstrates the threat posed by this pest to ecologically similar areas in Europe and North America.

K E Y W O R D S
Bactrocera dorsalis, geometric morphometrics, microsatellites, mitochondrial genes, population structure markers (cox1 and nad6 genes and microsatellite loci) to determine global population structuring in this species. Three independent markers were used in an integrative framework (Schlick-Steiner et al., 2010), and in alignment with previous studies of B. dorsalis that have shown these markers to be informative at different temporal and spatial scales of population structuring and invasion biology (Boontop, Schutze, Anthony, Cameron, & Krosch, 2017;Schutze, Krosch et al., 2012;Shi, Kerdelhue, & Ye, 2012). Using DIYABC analysis, our global population data also allow us to make informed comment on the likely origin of B. dorsalis within the Indo/Asian region (an issue under debate (Choudhary, Naaz, Prabhakar, & Lemtur, 2016)) and global invasion pathways.
A second major component of our study is to document morphological and genetic changes associated with the current, ongoing northward invasion of B. dorsalis into Central China. Although well-documented in tropical and sub-tropical China (Wan, Nardi, Zhang, & Liu, 2011), B. dorsalis was historically absent from Central China because of climatic unsuitability (specifically cold stress (Stephens, Kriticos, & Leriche, 2007;De Villiers et al., 2016)).
Nevertheless, B. dorsalis is now able to successfully overwinter in central Chinese provinces, such as Hubei Province (Han et al., 2011). This poses a great concern not only for China, but must also to temperate Europe and North America. Climate models predict these regions to be "unsuitable" for B. dorsalis (De Villiers et al., 2016) but, given the Chinese situation, must now be con-

| Sample collection
Two-thousand eight-hundred and sixty-seven B. dorsalis adults were collected from 63 locations and assigned into six a priori groups  Table S1, and the localities are shown in Figure 1. All samples were identified using available taxonomic keys prior to conducting molecular analyses (Liang, Yang, Liang, Situ, & Liang, 1996;White & Elson-Harris, 1992). The sampling involved the use of maleonly lures, so females were not collected. The legs or a portion of the body were removed for genetic analysis and one wing (usually the right) for geometric morphometric shape analysis. The rest of the body and the DNA were stored at −20°C with voucher references for morphological and molecular verification at the Plant Quarantine and Invasion Biology Lab in China Agricultural University.

| Geometric morphometric analysis
Usually, the right wing was dissected from each fly for slide mounting, image capture and analysis, the left was used instead if the right wing was damaged. Wings were slide mounted using DPX mounting agent and air-dried prior to image capture using an AnMo Dino-Eye microscope eyepiece camera (model # AM423B) mounted on a Leica MZ6 stereomicroscope. Fifteen wing landmarks were selected following Schutze, Jessup, & Clarke, 2012 and digitization using tps-DIG2 v2.16 (Rohlf, 2010).
Raw landmark coordinate data were imported into the computer program MORPHOJ v1.04a (Klingenberg, 2011) for shape analysis. Data were first subjected to Procrustes superimposition to remove all but shape variation (Rohlf, 1999). Multivariate regression of the dependent wing shape variable against centroid size (independent variable) was conducted to assess the effect of wing size on wing shape (i.e., allometry) (Drake & Klingenberg, 2008;Schutze, Jessup et al., 2012). The statistical significance of this regression was tested by permutation tests (10,000 replicates) against the null hypothesis of independence. Subsequent analyses used the residual components as determined from the regression of shape on centroid size to correct for allometric effect.
The size of each wing (centroid size) was calculated in MORPHOJ v1.04a. Centroid size is an isometric estimator of size calculated as the square root of the summed distances of each landmark from the center of the landmark configuration. ANOVA with post hoc Tukey's test was used to assess for significant differences among sample sites.
Samples were a priori assigned to the above six groups (as for centroid size analysis) and 16 provinces within China, from which subsequent canonical variates analysis (CVA) was applied to determine relative differences in wing shape among groups (Krosch et al., 2013). Significant differences were determined via permutation tests (1000 permutation rounds) for Mahalanobis distances among groups. We regressed geographic distance (km) between each pair of sampling sites determined by Google Earth 7.1.7.2606 against Mahalanobis distances calculated from CVA to test for "isolation-bydistance" (IbD) effects (Wright, 1943). The strength of the association was determined by linear regression analysis using the program SPSS v17.0.

| Genetic analysis
Total genomic DNA was extracted from each individual fly following the manufacturer's protocol from the TIANamp Genomic DNA kit (DP304, TIANGEN, China) for animal tissue, and slight modifications were made to increase DNA concentration (Jiang et al., 2013).
Gene amplification and sequencing methods were reported previously (Jiang et al., 2013). The primers designed for this study are shown in Supporting Information Table S2. Both directions of the cox1 (divided into two fragments at first) and nad6 sequences from each individual were reviewed using Chromas (version 2.33) and assembled using DNAMAN 5.
For the sequences, the nucleotide composition and variable positions were visualized using MEGA 7 (Kumar et al., 2016). The nucleotide diversity (π), haplotype diversity (Hd), and the number of haplotypes were estimated using DnaSP 6 (Rozas et al., 2017).

| Population genetic structure
Pairwise F ST was calculated for both types of markers using Arlequin 3.5 (Excoffier & Lischer, 2010) to measure the degree of genetic  Table S1. Note: Insert figure: China, Hawaii. The map was created in ArcGIS 10.2 software (ESRI Inc., Redlands, CA, USA). URL http://www.esri.com/sofware/arcgis/arcgis-for-desktop differentiation between pairs of populations. We grouped populations into the six prior groups to explore differences between specific regions.
Isolation by distance (IBD) was examined by testing the correlation between F ST and geographic distances using the program SPSS v17.0.
Bayesian clustering of individuals based on microsatellite genotypes was performed in STRUCTURE 2.0 (Pritchard, Stephens, & Donnelly, 2000) to infer genetic structure among the six a priori defined groups and 63 populations of B. dorsalis. We set the number of clusters (K) from 1 to 10 and conducted 10 independent runs for each value of K. Each run consisted of a burn-in period of 50,000 steps, followed by 100,000 Markov chain Monte Carlo (MCMC) repetitions with a model allowing admixture. ΔK values (Evanno, Regnaut, & Goudet, 2005) were computed to select the most likely number of K using the online resource Structure Harvester (Earl, 2011) that explained the structure in data. We then conducted model to summarize cluster membership coefficient matrices for each value of K with CLUMPP 1.1.2 (Jakobsson & Rosenberg, 2007), and plotted using DISTRUCT 1.1 (Rosenberg, 2004).

| Geometric morphometric analysis
One-thousand two-hundred and sixteen wings from flies collected from 63 locations covering six a priori defined geographic groupings ( Figure 1 and Supporting Information Table S1) were used for geometric morphometric analysis. Wing centroid size varied significantly across the six population groups (F 5, 1210 = 6.128; p < 0.05), with significant variation found between populations from Central China  (Table 1). For the full dataset, a significant isolation-by-distance relationship (Pearson correlation = 0.514; p < 0.01) was detected between geographic distances and Mahalanobis distances (Supporting Information Figure   S1).

| Population genetic structure
Significant genetic differentiation based on among-site F ST indices was observed between all regions and across all molecular markers (  (Table 3). Isolation by distance (IBD) was detected from the two sets of molecular markers (Supporting Information Figure S1). Within China, genetic differentiation was low among all sample sites (Supporting Information Tables S4 and S5).  including many that were central to the network (Figure 4a,b).
Further, the four Asian groups showed high haplotype diversity relative to Africa and Hawaii. There was no evidence for macrogeographic population structure among regions. Similarly, there was no apparent structure observed within China, with many haplotypes shared among provinces/cities, and with all locations showing similarly high haplotype diversity (Figure 4c, d).

| Demographic history
Nine scenarios were assessed in DIYABC analysis (Supporting Information Figure S2) Information Table   S6). Furthermore, the unimodal mismatch distribution (Supporting Information Figure S4) supported a model of population expansion (p SSD > 0.05).

| D ISCUSS I ON
This study represents the largest-ever sampling of B. dorsalis populations and encompasses both the accepted native range and recent (<70 years) invasive locations. This extensive sampling provides higher resolution of population relationships and more rigorous tests regarding invasive history and region of origin of this highly invasive and pestiferous fruit fly than previously available.

| Global population structure
Microsatellite markers divided B. dorsalis into two genetic units, an Asian group and a non-Asian group, while the network of mitochondrial haplotypes did not suggest any subdivision in the data.
The non-Asian group inferred in the microsatellite data corresponded to known recent invasive populations in Africa and Hawaii.
Interpretation of this pattern would normally imply recent common TA B L E 2 Genetic diversity indices in six groups of Bactrocera dorsalis π: nucleotide diversity. ancestry between these locations; however, it is much more likely in this case that the pattern shows the effect of separate founder events in each location (Nardi, Carapelli, Dallai, Roderick, & Frati, 2005;Sakai et al., 2001). Both Africa and Hawaii have much lower molecular diversity than Asia overall, and the two locations are supported by ABC analysis as arising from a South Asian source in separate colonization events. Further, genetic differentiation indices support these populations as significantly different to most Asian locations, and also to each other. Overall, we argue that these data strongly suggest that Hawaii and Africa are the result of single separate invasions, with little or no subsequent gene flow with source populations in South Asia. This supports previous analyses of African populations of B. dorsalis (Khamis et al., 2009;, but suggests also that founder effects in separate invasive populations that arose from the same source can cause these populations to appear similar under some analytical scenarios. Within Asia, weak genetic structure and/or isolation-by-distance trends have been recorded in all previous studies (Aketarawong et al., 2007;Schutze, Krosch et al., 2012;Shi et al., 2012;Wan, Liu, & Zhang, 2012) and were explained by repeated long-distance migration events, facilitated by the polyphagy of the fruit fly. Our genetic and morphological results within Asia agree with this hypothesis and infer that there is no macrogeographic sub-structuring across the Asian region. Microsatellite data suggested some individual  Figure 1 and Supporting Information Table S1,  Southern China had been inferred as the potential origin by previous authors (Aketarawong et al., 2007;Li, Wu, Chen, Wu, & Li, 2012;Schutze, Krosch et al., 2012;Shi et al., 2012;Wan et al., 2012); however, these studies lacked critical Indian/South Asian populations.
Recent work within India has shown high genetic diversity within Indian populations (cox 1 Hd 0.833-1.00), but there was very limited integration of that data with non-Indian populations (Choudhary et al., 2016). Our global population data allow us to test multiple global invasive pathway scenarios under an ABC framework, and these tests continue to support the hypothesis of South Asia (= India + Bangladesh) as the original location, even with extensive sampling through the rest of Asia.

| Implications for management
Despite our extensive sampling across B. dorsalis's entire geographic range, no macrogeographic population structuring was observed: a result fully consistent with prior studies which have covered individual components of the range (Choudhary et al., 2016;Khamis et al., 2009;Schutze, Krosch et al., 2012;Shi et al., 2012). The genetic uniqueness of Hawaii and Africa is linked to the recent invasions (60-15 years, respectively, assuming the invasive populations were detected soon after their establishment) of those locations, not because of long-term population structuring. Given such consistent results from multiple studies, we consider it highly unlikely that different geographic populations of the fly will show marked biological differences, such as differences in host use or thermal tolerances, which might impact on quarantine, trade, or pest management.  Walsh et al., 2011).
Yet, as damaging as spotted-wing drosophila is, its impact is still largely restricted to berries and soft fruits. In contrast, B. dorsalis attacks fruits from well over 20 plant families and is generally regarded as one of global agriculture's most damaging insects (Clarke et al., 2005). Should it permanently establish in North Asia, America or Europe, we anticipate a far more northerly spread than currently predicted by climate matching models (De Villiers et al., 2016).

ACK N OWLED G EM ENTS
The authors would like to thank Abdeljelil Bakri, Yubing Huang,