Low genetic diversity but strong population structure reflects multiple introductions of western flower thrips (Thysanoptera: Thripidae) into China followed by human‐mediated spread

Abstract Historical invasion scenarios based on observational records are usually incomplete and biased, but these can be supplemented by population genetic data. The western flower thrips (WFT), Frankliniella occidentalis, invaded China in the last 13 years and has rapidly become one of the most serious pests in the country. To assess whether this invasion involved a single event or multiple events, we examined patterns of genetic diversity and population structure of WFT across 12 Chinese populations and a native US population based on mitochondrial DNA and/or 18 microsatellite loci. The average allelic richness and haplotype diversity in Chinese populations were significantly lower than in a population from its native range. The distribution of mitochondrial haplotypes suggested multiple independent invasions of WFT into China, including two invasions into the Beijing region. Based on microsatellite data, two distinct clusters were identified, with both of them splitting further into two clusters; in the Beijing region, the microsatellite data also provided evidence for two introductions. Both the absence of isolation by distance and the fact that distant populations were similar genetically suggest patterns of WFT movement linked to human activities. Our study therefore suggests multiple introductions of WFT into China and human‐assisted spread.

has rapidly become one of the most serious pests in the country. To assess whether this invasion involved a single event or multiple events, we examined patterns of genetic diversity and population structure of WFT across 12 Chinese populations and a native US population based on mitochondrial DNA and/or 18 microsatellite loci. The average allelic richness and haplotype diversity in Chinese populations were significantly lower than in a population from its native range. The distribution of mitochondrial haplotypes suggested multiple independent invasions of WFT into China, including two invasions into the Beijing region. Based on microsatellite data, two distinct clusters were identified, with both of them splitting further into two clusters; in the Beijing region, the microsatellite data also provided evidence for two introductions. Both the absence of isolation by distance and the fact that distant populations were similar genetically suggest patterns of WFT movement linked to human activities. Our study therefore suggests multiple introductions of WFT into China and human-assisted spread.

K E Y W O R D S
Frankliniella occidentalis, invasive species, microsatellite, mitochondrial gene, multiple introductions, population genetics

| INTRODUCTION
Biological invasions, as by-products of travel and trade, are becoming a worldwide environmental, health and agricultural issue (Pimentel et al., 2001;Ricciardi, 2007). Although many invasion events recorded in the literature are based on sample collections, introductions and dispersal routes are often undetectable by direct observations, particularly when introductions have taken place following failure in quarantine (Estoup & Guillemaud, 2010). Studies on the genetic diversity and population structure of an invasive species in both its introduced and native area can help to elucidate the source, path and times of introduction (Boissin et al., 2012;Lombaert et al., 2014;Rollins, Woolnough, Wilton, Sinclair, & Sherwin, 2009) and provide critical information on management and prevention of future invasion (Signorile et al., 2014;Veale et al., 2013;Zhang, Edwards, Kang, & Fuller, 2014). Multiple introductions, secondary expansion of introduced populations and control efforts can all lead to structured populations in introduced regions (Berthouly-Salazar et al., 2013;Cao, Wei, Hoffmann, Wen, & Chen, 2016;Konecny et al., 2013), while gene flow between introduced populations eventually decreases levels of differentiation (Tsuchida, Kudo, & Ishiguro, 2014) except at loci under selection.
The western flower thrips (WFT), Frankliniella occidentalis (Pergande) (Thysanoptera: Thripidae), is a global invasive pest causing serious damage to many crops (Childers, 1997;Rotenberg, Jacobson, Schneweis, & Whiffleld, 2015). It originated in southwestern USA (Morse & Hoddle, 2006) and has invaded other areas around the world since the late 1970s (Kirk & Terry, 2003). WFT was introduced into the eastern states of USA and Canada before 1995, most states in Australia before 1996, nearly all countries across Europe before 1998, Central and South America and four countries in Africa before 2001 and four countries in Asia before 1994 (Kirk & Terry, 2003). In Asia, it was first reported in Israel in 1987, peninsular Malaysia in 1989, Japan in 1990and South Korea in 1994(Kirk & Terry, 2003. The invasion success of WFT is likely to reflect high levels of polyphagy, resistance to most classes of insecticides and rapid population growth due to a short generation time and parthenogenetic reproduction (Gao, Lei, & Reitz, 2012). Establishment and expansion of WFT are thought to be enhanced by rapid cropping rotations and the presence of glasshouses, which provide continuous food and habitats (Wan & Yang, 2016).
The invasion of WFT into China has occurred relatively recently.
The first report of WFT establishment was in 2003 on glasshouse pepper in a suburban area of Beijing in northern China (Zhang, Wu, Xu, & Zhu, 2003). Prior to this time, WFT was first found in 2000 at an international exhibition on floriculture in Kunming, Yunnan province, in southwestern China (Jiang, Bai, Xiao, & Yang, 2001). While these locations are separated by 2,200 km, WFT began to cause heavy damage in these two regions almost simultaneously before spreading to other areas (Ren, Lei, Zhang, Hong-Ling, & Hua, 2006;Wu, Zhang, & Li, 2009;Zhang et al., 2003;Zheng, Liu, Zhang, & Zhao, 2007). WFT has rapidly become a common pest of vegetable crops and flowers in nearly all provinces across China in 8 years since the first report of its establishment (Chen, Yuan, Du, Zhang, & Wang, 2011;Zhang et al., 2003).
While patterns of establishment of WFT in China may point to two invasions in Beijing and Yunnan, it is also possible that southwestern populations in China served as a single point of origin. Previous genetic studies based on microsatellite markers and mitochondrial variation suggested low levels of genetic differentiation of WFT among populations in China in the early phase of establishment (Yang, Sun, Xue, Li, & Hong, 2012;Yang et al., 2015). However, there is variation among Chinese WFT populations in their level of resistance to insecticides, even within the Beijing region, which points to some level of genetic differentiation among populations (Wang, Hou, 2011;Wang, Liu, et al., 2011;Wang et al., 2016).
To clarify pathways of invasion in China, we characterize the genetic diversity and population structure of WFT using new genetic markers, and we also consider variation in a population from the native range of WFT. We use 18 newly developed microsatellite loci selected by stringent criteria (Cao, Li, et al., 2016) as well as a segment of a mitochondrial gene to characterize patterns of variation. The results point to high levels of differentiation among some populations suggestive of multiple origins and long distance gene flow. (558-3,081 km) as well as six sites within the vicinity of Beijing (separated by 16-94 km). Given limited dispersal of WFT (Rhainds & Shipp, 2004) and the fact that semirural areas in Beijing are typically separated by residential housing developments, and topographical variation within the Beijing region, there is potential for a degree of isolation among WFT from the six sites sampled in the Beijing area. The thrips were collected from flowers by beating the flowers over a white enamel tray, and each thrips used in this study came from a different host plant.

| Sample collection and genotyping
All specimens were preserved in absolute ethanol and stored at −80°C prior to DNA extraction. In total, 411 female WFT were used.
Genomic DNA was extracted from individuals with the DNeasy Blood and Tissue Kit (Qiagen, Germany). Eleven populations were genotyped for both nuclear and mitochondrial markers, while the other two populations from Beijing (BJDX and BJCP) were genotyped for mitochondrial variation (insufficient DNA was available from these populations for microsatellite analyses). For nuclear markers, 18 microsatellite loci designed by Cao, Li, et al. (2016) were validated in our study following the methods and criteria they proposed (   Nei, 1978). Furthermore, we compared the number of alleles (A S ) and

| Genetic diversity analysis
heterozygosity (H ES ) among samples with different sample sizes using a rarefaction method in GENCLONE. Allelic richness (A R ) and allelic richness of private alleles (P AR ) were calculated with a rarefaction approach in HP-RARE version 1.1 (Kalinowski, 2005) on a minimum sample size of 11 diploid individuals.
FreeNA was used to estimate the frequency of null alleles and calculate the adjusted fixation index (F ST ) using the excluding null alleles correction method with 10,000 bootstraps (Chapuis & Estoup, 2007).
We also constructed split networks to reveal relationships among mitochondrial haplotypes using SPLITSTREE version 4.13.1 (Huson & Bryant, 2006) with the neighbor-net method under a distance model of K2P after 1,000 bootstraps. A maximum-likelihood phylogenetic tree was con-

| Population genetic structure
Two Bayesian cluster methods (STRUCTURE and GENELAND) and a discriminant analysis of principal components (DAPC) were used to investigate genetic structure across the populations. Population structure was first characterized with the Bayesian model-based clustering method implemented in STRUCTURE version 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). An admixture model with correlated allele frequencies was chosen. Thirty replicates for each K (from 1 to 10) were run with 200,000 Markov chain Monte Carlo (MCMC) iterations after a burn-in of 100,000 iterations. The outputs of STRUCTURE were submitted to STRUCTURE HARVESTER WEB version 0.6.94 (Earl & Vonholdt, 2012) to estimate the optimal value of K using the Delta (K) method (Evanno, Regnaut, & Goudet, 2005). The membership coefficient matrices (Q-matrices) of replicated runs for each K were combined using CLUMPP version 1.1.2 (Jakobsson & Rosenberg, 2007) with the Greedy algorithm, and then visualized using DISTRUCT version 1.1 (Rosenberg, 2004). To investigate more subtle substructuring, additional analyses were carried out on each cluster separately.
We used GENELAND to estimate the number of clusters and their spatial patterns. The analysis was conducted over 10 replicates for each K value (from 1 to 15) with the null allele model selected to infer the number of clusters with 1,000,000 MCMC iterations-of which every 1,000th was retained. The replicates with the highest mean logarithm of posterior probability were used to compute posterior probabilities of population membership for each pixel of the spatial domain with a burn-in of 500.
To identify the number of different genetic clusters in the microsatellite data, DAPC was performed using adegenet version 2.0.1 (Jombart, 2008) in the R environment.

| Gene flow analysis
Recent migration rates among Chinese populations were estimated using BAYESASS 3.0.4 (Wilson & Rannala, 2003) based on microsatellite data. We performed preliminary runs (10,000,000 steps) to adjust mixing parameters for allele frequencies and inbreeding coefficients.
We then conducted ten longer runs of 100,000,000 steps using different start seeds with sampling every 1,000 steps. Trace files of the ten longer runs were combined to calculate mean migration with a burn-in of 10,000,000 using TRACER 1.6 (Rambaut, Suchard, & Drummond, 2013).

| Genetic diversity based on mitochondrial gene
We observed 18 mitochondrial haplotypes from the WFT populations

| Phylogenetic relationships
In an unrooted NJ tree constructed using the microsatellite loci

| Isolation by distance
The Mantel tests based on microsatellite data sets did not reveal a significant association between genetic and geographical distances

| Gene flow
The estimates of gene flow obtained from the BAYESASS analysis were consistent among 10 independent runs, and mean migration rates of each population are provided as a heat map ( Figure 5). A low level of gene flow was found between C1 and SB (from C1 to SB, 0.0084 on average; from SB to C1, 0.0097 on average), while a relatively higher rate of gene flow was detected within groups (within C1, 0.0255 on average; within SB, 0.0561 on average). Gene flow from or into JSYZ was estimated as consistently low.

| Population genetic structure of WFT after invasion
Invasive species may show genetic structure in their introduced regions, which might be caused by multiple introductions (Konecny et al., 2013), dispersal patterns following invasion (Berthouly-Salazar et al., 2013) and/or control efforts . In this study, we have used microsatellite markers and mitochondrial DNA sequences to examine the genetic structure of WFT in China. We found structured populations of WFT in the invasive area, in contrast to previous studies based on mitochondrial DNA and microsatellites isolated from an SSRenriched library or expressed sequence tags (Yang et al., 2012(Yang et al., , 2015. The reason for these inconsistent results is not entirely clear but may reflect the nature of the markers used and sampling variation. The microsatellite markers in this study were isolated from a draft WFT genome with stringent criteria, including the highest amplification rate during PCR, avoidance of single nucleotide insertions, and unambiguous peaks during genotyping, as well as the use of markers that did not deviate strongly from HWE and that were not in linkage disequilibrium. These criteria were adopted to ensure that the markers could minimize genotyping errors and reliably describe population structure, resulting in higher levels of population differentiation (F ST ) and improving assignment of individuals to populations (Cao, Li, et al., 2016).
In terms of sampling sites, we included a population from eastern China (JSYZ) as well as more populations from the Beijing area where establishment of WFT was first reported, in contrast to a more limited set of locations considered in two previous studies (Yang et al., 2012(Yang et al., , 2015. This allowed us to highlight the fact that the population from eastern China was distinct and also that the four populations from Beijing had a nested hierarchical structure and that three of these populations separated from the six other Chinese populations.
Genetic differentiation was low and gene flow was high among the six geographically separate populations (YNHH, GZGY, LNCY, XJWS, XZLS, and BJYQ). This was largely consistent with previous studies that found low genetic differentiation and strong evidence for gene flow among geographically separate populations (Yang et al., 2012(Yang et al., , 2015 that included locations near those considered here. The different conclusions of our study compared with the earlier work may therefore reflect the additional samples included in the current study.

| Mitochondrial haplotypes
The population of WFT that has spread around the world is considered to be an insecticide resistant strain or strains from California established in glasshouses (glasshouse strain, WFTG) (Kirk & Terry, 2003). There is also evidence for a more susceptible "lupin strain" (WFTL) that spread to New Zealand prior to the spread of WFTG (Martin, Workman, & Butler, 2005;Rugman-Jones et al., 2010). Based on mitochondrial DNA sequences and complete linkage disequilibrium between nuclear ribosomal DNA and the mitochondrial DNA, these two distinct lineages found both within and outside the native range of WFT were regarded as sympatric cryptic species by Rugman-Jones et al. (2010). However, the two lineages have also been described as two allopatric ecotypes with distinct mitochondrial DNA and nuclear DNA sequences-a hot/dry (HD) ecotype and a cool/moist (CM) ecotype (Brunner & Frey, 2010).
The two phylogenetic lineages in this study correspond to WFTG/HD (Hap1-Hap2, Hap4-Hap15) and WFTL/CM (Hap3 and Hap16) based on mitochondrial DNA sequence variation, but we were unable to link mitochondrial DNA variation with nuclear marker variation ( Figure 3). This may suggest hybridization of the two forms in China (Yang et al., 2012) in contrast to the situation in California (Brunner & Frey, 2010;Rugman-Jones et al., 2010). The common haplotypes found in this study were similar to those found in previous studies (Yang et al., 2015), except for haplotype 8, which was common in five populations from Beijing. This haplotype was absent from other Chinese populations (except for one individual from LNCY) or from populations in the native range of WFT that have been characterized (Brunner & Frey, 2010

| Microsatellite data
The genetic structure identified from microsatellite data, which was largely in accordance with that identified from the mitochondrial DNA  (Ren et al., 2006;Wu et al., 2009;Zhang et al., 2003 (Yang et al., 2012(Yang et al., , 2015. The southern Beijing populations and JSYZ formed two independent clusters in the DAPC and the second batch of STRUCTURE analyses. However, the POPTREE analyses based on mitochondrial DNA showed that JSYZ was distant to the southern Beijing populations but closest to YNHH. The relatively high frequency of rare alleles for some microsatellite loci in the JSYZ population may reflect an independent introduction into this region. However, founder effects can also lead to genetic differentiation between invading and/ or expanding front populations (Pierce et al., 2014;Schulte, Veith, Mingo, Modica, & Hochkirch, 2013). Founder events are particularly likely to contribute to the genetic distinctness of JSYZ given that this population was sampled as soon as it was detected and therefore prior to ongoing gene flow reducing differentiation among populations.

| Dispersal of WFT after invasion in China
Western flower thrips across China showed no IBD even when populations with likely different origins were separately analyzed. The

| Genetic diversity of WFT
Although we were unable to pinpoint the origins of Chinese populations accurately because only one sample was available from its native range, diversity levels in USCA provide a reference point for comparison.
Average allelic richness and mitochondrial haplotype diversity of WFT populations in China were lower than in USCA and this may be a consequence of founder events during invasion. The average allele richness of WFT in Chinese populations was reduced to about half (60% in the case of LNCY) that found in USCA, while the average number of mitochondrial haplotypes was also reduced to 45.8% (up to 75% in GZGY and BJMT).
In total, 14 mitochondrial haplotypes have been found in China (12 haplotypes in this study), accounting for 24.6% of the haplotypes found so far in WFT (Brunner & Frey, 2010;Rugman-Jones et al., 2010;Yang et al., 2015).
Nevertheless, genetic diversity in WFT based on microsatellite markers and mitochondrial DNA was higher than in other invasive species introduced to China from America, such as the fall webworm Hyphantria cunea , red turpentine beetle Dendroctonus valens (Cai, Cheng, Xu, Duan, & Kirkendall, 2008) and black locust gall midge Obolodiplosis robiniae (Shang, Yao, Huai, & Zhao, 2015). For instance, only one mitochondrial haplotype and an allelic richness of 2.92 were observed in 24 populations of fall webworm in China, in contrast to 12 haplotypes and an allelic richness of 6.29 in two populations from North America. The relatively high level of genetic diversity in invasive WFT populations may explain the relatively rapid development of resistance to insecticides in Chinese WFT populations  and other forms of local adaptation.

| Implications for WFT management
Understanding the sources and likely dispersal pathways of invasive species is crucial for accurate risk assessment and for successful management (Stepien, Brown, Neilson, & Tumeo, 2005).

| CONCLUSION
In conclusion, we found that introduced WFT populations in China were genetically differentiated in contrast to findings from previous work.
Both the mitochondrial DNA and microsatellite markers point to multiple introductions of WFT. Populations in the Beijing region arose from at least two independent introduction events, while populations in other areas likely resulted from secondary expansions out of Yunnan Province.