Introduction
Three invasive leafminer species, Liriomyza huidobrensis (Blanchard), Liriomyza sativae Blanchard, and Liriomyza trifolii (Burgess), are damaging too many ornamental and vegetable crops (Spencer, Reference Spencer1973, Reference Spencer1990; Reitz et al., Reference Reitz, Kund, Carson, Phillips and Trumble1999; Shiao, Reference Shiao2004). In China, L. sativae was first found in Hainan in 1993 (Kang, Reference Kang1996), and then dispersed across most of provinces in mainland China within 5 years (Wang & Zhao, Reference Wang and Zhao1998). The L. sativae is especially difficult to control due to its polyphagous nature, high reproductive rate, and short development period. In recent decades, many studies have focused on the biology (Petitt & Wietlisbach, Reference Petitt and Wietlisbach1994; Araujo et al., Reference Araujo, Nogueira, Menezes Netto and Bezerra2013), ecology (Yildirim et al., Reference Yildirim, Civelek, Dursun and Eskin2012; Tian et al., Reference Tian, Zhang, Li and Cheng2013), and control (Saberfar et al., Reference Saberfar, Garjan, Naseri and Rashid2012; Saryazdi et al., Reference Saryazdi, Hejazi, Rashidi and Ferguson2014) of L. sativae whether in China or other countries. As for genetic aspects of L. sativae, there are several researches on this topic (Scheffer, Reference Scheffer, Wiegmann and Yeates2005; Scheffer & Lewis, Reference Scheffer and Lewis2005; Scheffer et al., Reference Scheffer, Lewis and Joshi2006; Amin et al., Reference Amin, Scheffer, Lewis, Pasha and Bhuiya2014), and we previously reported its complete mitochondrial genome (Yang et al., Reference Yang, Du, Wang, Cao and Yu2011, Reference Yang, Du, Cao and Huang2013), developed microsatellite markers (Ji & Du, Reference Ji and Du2013) and preliminarily analyzed the genetic differentiation of the host- and geo-populations of L. sativae based on mitochondrial DNA (Du et al., Reference Du, Wang, Lu, Zheng and Lu2008, Reference Du, Tang, Wang, Shen and Chang2014; Wang et al., Reference Wang, Du, He, Zheng and Lu2008). Given its economic and ecological importance for agriculture, more comprehensive investigations of genetic structure and variation of L. sativae populations throughout its range in China, as well as the mechanism of its dispersion and invasion, are needed.
In this study, we investigated the genetic diversity, population structure, and invasion mechanism of this species, based on 288 individuals from 12 populations of L. sativae in China, using polymorphic microsatellite markers we developed previously (Ji & Du, Reference Ji and Du2013).
Materials and methods
Sampling and DNA extraction
A total of 288 L. sativae adults were collected from common vegetables (cowpea, kidney bean and towel gourd) in 12 locations across China (table 1). These individuals were preserved in 100% ethanol at −20°C until DNA extractions were performed. Genomic DNA was extracted from L. sativae adults using an AxyPrep Multisource Genomic DNA Miniprep Kit (Axygen, Suzhou, China) as recommended by the manufacturer.
Table 1. Collection site, code, number of studied individuals, and host of Liriomyza sativae.
Microsatellite analysis
Variable microsatellite loci previously identified for L. sativae by Ji & Du (Reference Ji and Du2013) were examined and seven loci with strong, unambiguous banding patterns were selected for use in this study (Supplementary Table S1). These primers were attached to FAM, HEX and TAMRA fluorophores at the 5′ ends for genotyping. All PCR amplifications were performed in 25 µl reaction volumes containing 5 µl of 10 × PCR buffer, 50 ng of template DNA, 1.5 mM of MgCl2, 2.5 µl of dNTPs (2.5 mM each), 1 µl each of the primer (10 µM), and 0.75 U of rTaq DNA polymerase (Takara, Dalian, China). PCR amplifications were performed, which included an initial denaturation at 94°C for 4 min, followed by 42 cycles of 50 s at 94°C, 50 s at 54–63°C depending on the primer pair (Supplementary Table S1), 1 min at 72°C, and a final extension for 10 min at 72°C. PCR products were analyzed using an ABI 3730XL DNA sequencer. Electropherograms were derived using Gene Scan 4.0 and used to deduce DNA fragment sizes using Gene Mapper 4.0 (Sangon Biotech, Shanghai).
Data analysis
Genetic polymorphism for each population was assessed by calculating the number of alleles (A), the effective number of alleles (Ae), and observed (HO) and expected heterozygosities (HE), as well as Nei's genetic distance (1978) using POPGENE version 1.31 (Yeh et al., Reference Yeh, Yang and Boyle1999).
In addition, the number of private alleles (Ap) was calculated by CONVERT 1.31 (Glaubitz, Reference Glaubitz2004). F ST (Differentiation index) was performed by ARLEQUIN 2.0 (Excoffier et al., Reference Excoffier, Laval and Schneider2005). For each population-locus combination, departure from Hardy–Weinberg expectation was assessed using exact tests (Guo & Thompson, Reference Guo and Thompson1992), with unbiased P-values estimated through a Markov-chain method (Guo & Thompson, Reference Guo and Thompson1992). The null allele frequency was calculated using Micro-Checker (Van Oosterhout et al., Reference Van oosterhout, Hutchinson, Wills and Shipley2004). Nei's genetic distance (1978) was then imported into PHYLIP computer package version 3.66 (Felsenstein, Reference Felsenstein1995) from which a neighbor joining phenogram was generated using the program NEIGHBOR using the unweighted pair-group with the arithmetic mean (UPGMA) method. Bootstrap values were calculated using 1000 replicates.
A Mantel test for isolation by distance was tested using the IBDWS web service (Jensen et al., Reference Jensen, Bohonak and Kelley2005) with 10,000 randomizations. In addition, STRUCTURE 2.3.3 (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000; Falush et al., Reference Falush, Stephens and Pritchard2003) was used to determine a reasonable number of partitions (K) for the studied populations. Clustering results were then visualized by the program CLUMPP v1.1.2 (Jakobsson & Rosenberg, Reference Jakobsson and Rosenberg2007). In this clustering analysis, we specified an initial range of potential genotype clusters (K) from 1 to 10 under the admixed model and the assumption of correlated allele frequencies among populations. For each value of K, ten runs were performed with 1,000,000 iterations discarded as burn-in, followed by an additional 1 million iterations. The most probable number of K values in the data was detected both by comparing the log probability of the data lnP (D) for each value of K across all ten runs of Structure and by examining the standardized second-order change of lnP (D), ΔK (Evanno et al., Reference Evanno, Regnaut and Goudet2005). We also carried out an independent analysis of spatial structure using the R package GENELAND 3.1.4 (Guillot et al., Reference Guillot, Estoup, Mortier and Cosson2005a , Reference Guillot, Mortier and Estoup b ), which has been explicitly tested for robustness in the presence of null alleles. Like STRUCTURE, the software uses a markov chain monte carlo (MCMC) strategy to determine the most likely number of populations (K), and assigns individuals to the most appropriate population based on individual multi-locus genotypes. We carried out ten independent MCMC simulations (10,000,000 iterations, thinnings of 10,000 iterations) using the spatial model, null allele model and correlated allele frequency model to analyze, during which we allowed K to vary between 1 and 12. Analysis of molecular variance (AMOVA) analyses were performed using ARLEQUIN 2.0 (Excoffier et al., Reference Excoffier, Laval and Schneider2005). We used MIGRATE version 3.6 (Beerli, Reference Beerli2006) to estimate the effective number of migrants (M = m µ−1) entering and leaving each population per generation to verify if there was asymmetrical gene flow between populations. We relied on Bayesian search strategy and used ten short chains, three long chains with 10,000 trees discarded as initial ‘burn-in’, replicates = YES: 5, randomtree = YES, heating = ADAPTIVE: 1{1 1.2 1.5 3.0}, and ran migrate four times, verifying consistency in our results. The estimates from the final run are reported here. The method described Wilcoxon's signed-rank test (Piry et al., Reference Piry, Luikart and Cornuet1999), which were performed using the Bottleneck version 1.2.02 (Piry et al., Reference Piry, Luikart and Cornuet1999), to evaluate whether the natural populations examined here had experienced recent bottlenecks under three mutational models: an infinite allele model (IAM), a stepwise mutation model (SMM) and a two-phase model (TPM).
Results
Genetic diversity and Hardy–Weinberg equilibrium (HWE)
All of the seven microsatellite markers (Supplementary Table S1) proved to be polymorphic and informative. In detail, the number of alleles per locus ranged from 12 (JY60) to 28 (JY4) with an average of 16.3. The expected heterozygosities (HE) ranged from 0.7366 (JY14) to 0.9033 (JY4), while the observed heterozygosities (HO) ranged from 0.2639 (JY14) to 0.7203 (JY60) (table 2). In addition, there were significant deviations from HWE at multiple loci from all sampling locations. For example, BJ and SXTY showed significant departure from HWE for only one locus after sequential Bonferroni's correction (P < 0.01), but five loci for SH, GZRJ and GXNN populations (P < 0.01). Of these loci, JY60 showed no significant departure from HWE (table 2), while JY4 showed significant departure from HWE in nine populations (P < 0.01), which may be the result of heterozygote deficiencies in all cases. What's more, null allele frequencies varied from −0.1767 (JY60) to 0.8507 (JY14), and three of the seven loci (JY15, JY60, and JY70B) showed low levels of estimated null alleles (table 2).
Table 2. Characterization of seven polymorphic microsatellite loci of Liriomyza sativae.
A, mean number of alleles; Ae, effective number of alleles; HO, observed heterozygosities; HE, expected heterozygosities; P-HW, test for Hardy–Weinberg equilibrium;
n.s., denotes no significant deviation from Hardy–Weinberg equilibrium.
N. f., null allele frequence.
*denotes a significant deviation from Hardy–Weinberg equilibrium (P < 0.05).
**denotes a significant deviation from Hardy–Weinberg equilibrium (P < 0.01).
The average number of alleles (A) per population varied from 5.429 (ZJHZ) to 8.714 (SH and GZRJ), and SH population also possessed the most effective number of alleles (Ae), whereas SXTY had the least. Private alleles were distributed in most of the populations except SXTY, ZJHZ, and HNDZ with the most (8) in the SH population (table 3). HE and HO values were moderate, ranging from 0.677 to 0.789 and 0.440–0.645, respectively. All values of HO were less than that of HE.
Table 3. Genetic diversity of 12 Liriomyza sativae populations based on seven microsatellite loci.
A, mean number of alleles; Ap, private alleles; Ae, effective number of alleles; HO, observed heterozygosities; HE, expected heterozygosities.
Mantel test for matrices of different loci
To explore if the dispersal of L. sativae was limited by distance, a Mantel test for isolation by distance was constructed. Our results showed that the isolation by distance correlation was insignificant and slightly negative (Z = 9.1746, r = −0.0814, P = 0.7610; fig. 1), suggesting a possible anthropogenic influence on the spread of the L. sativae in China, possibly the transportation of vegetables.
Fig. 1. Scatter plots of genetic distance versus geographical distance for pairwise population comparisons.
Population genetic structure
The genetic differentiation among 12 L. sativae populations in China was estimated by pairwise F ST values between each pair. All population pairwise comparisons were no more than 0.15, and most of them ranged from 0.05 to 0.15 (table 4), revealing only moderate genetic differentiation, since pairwise F ST is usually regarded as a standard measurement of population differentiation, and Wright (Reference Wright1938) assumed that 0.05 < F ST < 0.15 exhibited moderate level of genetic diversity. We also found no genetic differentiation between GZRJ and GDGZ populations (F ST = −0.00644). Conversely, the F ST value of SXTY was significantly higher (F ST = 0.16366) than that of ZJHZ (table 4). In addition, the molecular diversity between most of the populations was significant (P < 0.05) (table 4).
Table 4. Pairwise FST (below the diagonal) and P value (above the diagonal) of Liriomyza sativae.
Bold indicates F ST > 0.05 (below the diagonal) and P > 0.05 (above the diagonal).
The UPGMA tree generated from the genetic distance matrix (Nei) grouped the 12 L. sativae populations into two major groups. Unexpected, SH and GXNN populations formed one cluster and the remaining ten populations were assigned to the other cluster, which we named Group 1 (fig. 2).
Fig. 2. UPGMA dendrogram between 12 Liriomyza sativae populations.
Further STRUCTURE analysis of these 12 populations showed the same major patterns as those detected by the analysis of the UPGMA dendrogram. The highest ΔK value was obtained for K = 2 (ΔK = 469.25; fig. 3A). Fig. 3B shows the proportion of each population that contributed to each of the two clusters. Likewise, most individuals from the SH and GXNN populations were grouped into one group.
Fig. 3. Clustering analysis by STRUCTURE for full-loci dataset. (A) Inference of the number of genetic cluster (K) for Liriomyza sativae populations. (B) Proportion of the genome of each individual assigned to each of the two clusters. Each individual is represented by a vertical bar. (C) Proportion of the genome of each individual assigned to each of the nine clusters. Each individual is represented by a vertical bar.
The most likely number of inferred populations was K = 9 based on analysis in GENELAND. Although the results from GENELAND indicated a greater degree of population structure compared to the STRUCTURE results (K = 2), there were no qualitative disagreements between the two analyses. For example, SH and GXNN clustered together based on the GENELAND analysis (fig. 4) also clustered together in the STRUCTURE analysis (fig. 3), whereas, there were no obvious clusters when K = 9 in STRUCTURE (fig. 3C).
Fig. 4. (A) Posterior distribution of the estimated number of populations using GENELAND. (B) Population structure inferred in GENELAND at K = 9. The black dots indicate the sampling locations. The abscissa and ordinate show the coordinates of sampling locations.
The AMOVA analysis revealed that more than 89% of variation was among individuals within populations when K = 2. However, only a small portion of the variation attributed among populations was within groups (4.27%) and among groups (6.73). In addition, our results revealed significant genetic differentiation between L. sativae both among populations (F SC = 0.04579, P < 0.0001) and among groups (F CT = 0.06726, P < 0.05).
Patterns of gene flow
Estimates of gene flow calculated by MIGRATE indicated that the levels of gene flow between populations of L. sativae are generally high. Unidirectional estimates of M ranged from 40.1 (HNDZ→BJ) to 817.7 (SH→ZJHZ) (Supplementary Table S2). Of the 66 pairwise comparisons, 14 had asymmetrical gene flow, as indicated by non-overlapping 95% confidence intervals (CI) around the estimate of M into each population (indicated in bold in Supplementary Table S2). We also marked asymmetrical gene flow on the map (fig. 5). Overall, the trend of asymmetrical gene flow was from south to north, with southern populations having an asymmetric migration toward northern populations, especially from HNDZ to GZRJ and HUBJS, then to SXTY and HBBD, respectively (fig. 5).
Fig. 5. Asymmetric gene flow between Liriomyza sativae populations. Values of gene flow (M) were showed in Supplementary Table S2. Maps were created using Esri's ArcGIS platform (http://www.esri.com/software/arcgis). Arrows represent asymmetric gene flow.
Bottleneck test
Bottleneck analysis with 12 populations of L. sativae across China showed that none of these populations exhibited excess heterozygosity under the SMM. Only two populations (HBBD and ZJHZ) had a statistically significant excess of heterozygotes under the TPM. However, heterozygosity excess was observed in all populations except five (HNDZ, GDGZ, BJ, GZRJ, and SDJN) under the IAM (table 5), suggesting that most L. sativae populations that we studied had not undergone a genetic bottleneck.
Table 5. Bottleneck test for Liriomyza sativae populations.
Hde, heterozygote deficiency; Hex, heterozygote excess.
Discussion
Low genetic diversity and moderate genetic structure
All of the microsatellite markers exhibited high polymorphism (table 2), whereas A, Ap, and Ae values (table 3) suggests that the genetic diversity within each population of L. sativae in China is relatively low. The loss of genetic diversity is consistent with another research which showed that mitochondrial variation of L. sativae across sampled New World populations was higher than that from Old World populations (Scheffer & Lewis, Reference Scheffer and Lewis2005). Indeed, most successful invasive insect species show a reduction in genetic diversity from the native to invaded areas (Ahern et al., Reference Ahern, Hawthorne and Raupp2009; Lozier et al., Reference Lozier, Roderick and Mills2009; Chu et al., Reference Chu, Gao, De Barro, Wan and Zhang2011; Yang et al., Reference Yang, Sun, Xue, Li and Hong2012). The most likely reason for reduced variation of L. sativae is that initial introductions to other places such as South and Southeast Asia involved bottlenecks and L. sativae has then spread with reduced variation from those locations. What's more, genetic drift (Baker et al., Reference Baker, Loxdale and Edwards2003; Schmitt et al., Reference Schmitt, Cizek and Konvicka2005) and selection pressure (Suarez & Tsutsui, Reference Suarez and Tsutsui2008) may not be the reasons for that as this pest typically reached outbreak proportions shortly after they arrive (Kang, Reference Kang1996; Wang & Zhao, Reference Wang and Zhao1998). In addition, the reduced variation character is distinct from native species such as Sesamia inferens (Walker) (Tang et al., Reference Tang, Tao, Wang and Du2014a ) or another exceptional invasive species, Bactrocera dorsalis (Hendel) (Wan et al., Reference Wan, Nardi, Zhang and Liu2011), which exhibited fairly high levels of genetic diversity.
We noted that three of the 12 populations (SH, GZRJ, and GXNN) in this study deviated significantly from HWE (table 2), which is obviously associated with significant heterozygote deficiency. The HO of these populations was much lower than HE at less than 0.500 (table 3). This may arise from the Wahlund effect (http://www.dorak.info/genetics/popgen.html) related to recurrent inbreeding and subpopulation structure. It is reasonable to conclude that inbreeding has a major impact on heterozygote deficiency since sampling mainly from greenhouses restricts the randomness of mating, and since this species produces multiple generations in a year, unlike some other species with special reproductive strategies such as Aleurocanthus spiniferus (Quaintance) (Tang et al., Reference Tang, Xu, Sun, Xie and Du2014b , Reference Tang, Tao and Du2015).
Mutation/drift equilibrium of populations usually leads to an equal rate of heterozygosity excess or deficiency (Maruyama & Fuerst, Reference Maruyama and Fuerst1985). However, natural populations might not remain in a steady state, and when a population has contracted in size there is a transient deficiency in the number of alleles present in the population compared with that expected in a population in equilibrium that has an equivalent heterozygosity (Maruyama & Fuerst, Reference Maruyama and Fuerst1985). Thus, we could explore population dynamics by assessing heterozygosity excess or deficiency (Cornuet & Luikart, Reference Cornuet and Luikart1996). It is well known that the most useful markers for detecting bottlenecks are those evolving under IAM, TPM, and SMM (Di Rienzo et al., Reference Di Rienzo, Peterson, Garza, Valdes, Slatkin and Freimer1994; Cornuet & Luikart, Reference Cornuet and Luikart1996; Primmer et al., Reference Primmer, Saino, Møller and Ellegren1998; Estoup & Cornuet, Reference Estoup, Cornuet, Goldstein and Schlotterer2000). Of these, the strict SMM is obviously the most conservative model. Our results showed that heterozygosity excess was not obvious under the SMM for all the populations, although two and seven populations had a statistically significant excess of heterozygotes under the TPM and IAM, respectively. Considering the TPM is thought to more closely simulate microsatellite mutation (Primmer et al., Reference Primmer, Saino, Møller and Ellegren1998; Estoup & Cornuet, Reference Estoup, Cornuet, Goldstein and Schlotterer2000), coupled with the low genetic diversity, we may conclude that only a few populations experienced a population bottleneck (HBBD and ZJHZ) (table 5). These bottlenecks could be due to two potential factors. The first is interspecific competition among these three Liriomyza species. Recent surveys have found that L. trifolii has successfully replaced L. sativae and become the dominant population in most regions in China (Gao et al., Reference Gao, Reitz, Wei, Yu and Lei2012; Xiang et al., Reference Xiang, Lei, Wang and Gao2012; Wang, Reference Wang2013; Yi, Reference Yi2014). The second factor is the effects of mortality such as natural enemies or pesticides. Previous studies have shown that natural enemies or pesticides are effective in regulating Liriomyza populations whether in native or invaded regions (Johnson, Reference Johnson1993; Murphy & LaSalle, Reference Murphy and LaSalle1999; Rauf et al., Reference Rauf, Shepard and Johnson2000; Chen et al., Reference Chen, Lang, XU, He and Ma2003; Bai et al., Reference Bai, Gu, Xu, Hu, Hao and Liu2009), which may reduce the population size of L. sativae even though resistance to pesticides has become more prevalent.
Analysis of L. sativae population structure in China based on F ST and AMOVA analysis revealed moderate (0.05 < F ST < 0.15) but significant differentiation (global F ST = 0.10998, P < 0.001; F ST = 0.08008, P < 0.001) (tables 4 and 6). Such genetic diversity may be due to the relatively high number of generations of this species each year (for example, 14–17 generations per year in Guangdong Province), which may increase its genetic diversity to some extent. Moderate population subdivision exists despite the species’ capacity for long-distance dispersal, which mainly occurs through wind transport of adults or through transfer of infested host plants, especially Solanaceae, Cucurbitaceae, and Leguminosae. The SH population, with more private alleles (Ap = 8) seems to be an isolated population with a significantly F ST value than the other populations except GXNN, which is consistent with phylogenetic analysis (figs 2–4) but at odds with the great geographic distance between the two populations. This suggests that anthropogenic transport is the most likely explanation for the large-scale dispersion of L. sativae, as there was no isolation-by-distance effect found in the populations we studied. The moderate differentiation and the absence of isolation by distance may be due to the relatively high gene flow, especially over long distance. Another explanation for such moderate differentiation is the effect of host. On the basis of our previous studies, the trend of genetic differentiation in the host populations was consistent with the preference of L. sativae to the plant hosts based on mitochondrial cytochrome oxidase subunit I (mtDNA-COI) gene, the ribosomal internal transcribed spacer 1 (rDNA-ITS1) gene (Wang et al., Reference Wang, Du, He, Zheng and Lu2008) and β-tubulin gene (Du et al., Reference Du, Wang, Lu, Zheng and Lu2008).
Table 6. Results of AMOVA test of Liriomyza sativae based on seven microsatellite loci.
Migration patterns
The direction of gene flow was asymmetrical, generally from southern to northern populations (fig. 5), which potentially reveals the invasion route of L. sativae. This hypothesis is supported by the fact that L. sativae was first found in Hainan in 1993 (Kang, Reference Kang1996), and then dispersed across most of provinces in mainland China (Wang & Zhao, Reference Wang and Zhao1998). More importantly, we noticed that the Hainan population displayed two asymmetric migration routes north to the GZRJ and HUBJS populations, which then reached the SXTY and HBBD populations, respectively (fig. 5). Furthermore, we also found that southern populations of L. sativae did not show higher genetic diversity (table 3). There appears to be no Bridgehead effect for the L. sativae invasion, as the leafminer expanded quickly in mainland China and dispersed to 21 provinces in only 5 years (Wang & Zhao, Reference Wang and Zhao1998). This is supported by the high levels of gene flow (40.1–817.7, Supplementary Table S2). Above all, unlike other invasive species such as F rankliniella occidentalis (Yang et al., Reference Yang, Sun, Xue, Li and Hong2012) and Solenopsis invicta Buren (Ascunce et al., Reference Ascunce, Yang, Oakey, Calcaterra, Wu, Shih, Goudet, Ross and Shoemaker2011), the leafminer did not show an established a bridgehead population in Hainan Province, possibly because one of the main means of spreading of L. sativae is anthropogenic transport through the movement of horticultural products. Of course, we cannot rule out wind as a factor as well. Our study is the first to examine the migration patterns of L. sativae in China, which might help to improve management strategies for L. sativae, preventing contaminated plants from being transported to other regions.
Supplementary Material
The supplementary material for this article can be found at http://dx.doi.org/10.1017/S0007485315000905.
Acknowledgements
We sincerely thank Dr. Li-Ping Wang for collecting samples. English editing by Van Driesche Scientific Editing. We also thank two anonymous reviewers for providing useful comments on the manuscript. We express our deep gratitude to the Testing Center of Yangzhou University. This research was funded by National Science and Technology Support Program (grant no. 2012BAD19B06), Wuxi Science and Technology Support Program (grant no. CLE01N2012), and Jiangsu Science & Technology Support Program (grant no. BE2014410).
