Global mtDNA genetic structure and hypothesized invasion history of a major pest of citrus, Diaphorina citri (Hemiptera: Liviidae)

Abstract The Asian citrus psyllid Diaphorina citri Kuwayama is a key pest of citrus as the vector of the bacterium causing the “huanglongbing” disease (HLB). To assess the global mtDNA population genetic structure, and possible dispersal history of the pest, we investigated genetic variation at the COI gene collating newly collected samples with all previously published data. Our dataset consists of 356 colonies from 106 geographic sites worldwide. High haplotype diversity (H‐mean = 0.702 ± 0.017), low nucleotide diversity (π‐mean = 0.003), and significant positive selection (Ka/Ks = 32.92) were observed. Forty‐four haplotypes (Hap) were identified, clustered into two matrilines: Both occur in southeastern and southern Asia, North and South America, and Africa; lineages A and B also occur in eastern and western Asia, respectively. The most abundant haplotypes were Hap4 in lineage A (35.67%), and Hap9 in lineage B (41.29%). The haplotype network identified them as the ancestral haplotypes within their respective lineages. Analysis of molecular variance showed significant genetic structure (FST = 0.62, p < .0001) between the lineages, and population genetic analysis suggests geographic structuring. We hypothesize a southern and/or southeastern Asia origin, three dispersal routes, and parallel expansions of two lineages. The hypothesized first route involved the expansion of lineage B from southern Asia into North America via West Asia. The second, the expansion of some lineage A individuals from Southeast Asia into East Asia, and the third involved both lineages from Southeast Asia spreading westward into Africa and subsequently into South America. To test these hypotheses and gain a deeper understanding of the global history of D. citri, more data‐rich approaches will be necessary from the ample toolkit of next‐generation sequencing (NGS). However, this study may serve to guide such sampling and in the development of biological control programs against the global pest D. citri.


| INTRODUCTION
Genetic structure is the nonrandom distribution of alleles or genotypes in space or time (Mahy, Vekemans, & Jacquemart, 1999). Studies on population genetic structure can estimate genetic diversity and gene flow among populations, infer ancestral populations, and population level phylogeography. Such information is essential for informed management of pest species (Porretta, Canestrelli, Bellini, Celli, & Urbanelli, 2007). In recent years, understanding of biological invasion mechanism has increasingly relied on knowledge of the genetic structure of invasive species (Lee, 2002). Such studies can help to illuminate the basic biology and ecology of invasive species, the relationships between the intrinsic genetic characteristics of biological invaders, and their successful invasions (Xu, Zhang, Lu, & Chen, 2001). Furthermore, they can aid in the identification of effective biological control agents (Hoy, 2013).
The Asian citrus psyllid Diaphorina citri Kuwayama (Hemiptera: Liviidae) ( Figure 1) is a world-wild economic pest because it transmits the bacterial pathogen, Candidatus Liberibacter, that causes citrus greening disease (huanglongbing, HLB), considered one of the most serious diseases of citrus (da Graça, 2008). The psyllid was first described from Taiwan in 1907 (Halbert & Manjunath, 2004). However, the available evidence suggests that D. citri originates from the Indian subcontinent (Hollis, 1987) and has subsequently expanded to other subtropical and tropical regions of Asia and also to the Indian Ocean islands of Réunion and Mauritius (Waterhouse, 1998), and South America (Halbert & Núñez, 2004). Recently, it has reached Central American countries and southern USA (Halbert & Manjunath, 2004).
Environmental factors, such as temperature, humidity, host plants, and human activities, influence the survival and reproduction of D. citri (Wang, 2002;Waterhouse, 1998;Xu, Xia, & Ke, 1994;Yang, 1989;Zhao, 1987). However, D. citri is characterized by high reproductive output and rapid generation turnover (2-7 weeks) and may thus rapidly evolve in response to new environments after an invasion. If semi-isolated subpopulations acquire different ecological traits, such as in host preference, development time, pathogenicity, vector capacity, susceptibility to pesticides, and tolerance to heat, and other abiotic stressors (Carmichael et al., 2007;Jourdie et al., 2010;Pinho, Harris, & Ferrand, 2007;Remais, Xiao, Akullian, Qiu, & Blair, 2011), then genetically divergence among these subpopulations may rapidly emerge. Genetic studies have the potential to trace the D. citri demographic history that determines population genetic diversity of the invasive species and thus to illuminate geographic origin and secondary translocation events. However, to date, the global genetic differentiation, phylogeographic relationships, and invasion history were not resolved mainly owing limited geographic scope and relatively small sample sizes of prior studies (Boykin et al., 2012;de León et al., 2011;Guidolin, Fresia, & Cônsoli, 2014;Lashkari et al., 2014).
To assess the mtDNA population genetic structure, and potential dispersal routes of D. citri, we study global mtDNA genetic differentiation and phylogeography. We sequenced partial fragments of mitochondrial COI gene of recently collected individuals invasive to China and augmented our sampling with sequences from GenBank covering nearly the entire geographic range of D. citri. The current dataset covers 106 distinct localities around the globe. It thus offers a more holistic insight into D. citri matriline history than have prior studies.

| Insect collection and sequencing
We sampled D. citri in China between 2015 and 2016. A total of 412 specimens were taken from 11 sampling sites (Table S1 in Appendix S1). Collected insects were stored in 95% alcohol at −20°C. In total, 48 samples representing all our sampling sites were sequenced for genetic analyses. Additional sequences of D. citri from Iran, Pakistan, Saudi Arabia, India, Vietnam, Thailand, Indonesia, China, Reunion, Mauritius, Guadeloupe, Puerto Rico, Mexico, USA, and Brazil were taken from GenBank. We analyzed the COI sequences of the pest from a total of 356 colonies from 106 geographic sites in Asia, Africa, and America (Figures 2 and 3, and Table S1 in Appendix S1).
Amplification products were purified and sequenced with primers from the original PCR reactions on an ABI 3730 automatic sequencer (Applied Biosystems, Foster City, CA, USA), using BigDye technology. The chromatographs were interpreted with Phred and Phrap (Green, 2009;Green & Ewing, 2002) using the Chromaseq module (Maddison & Maddison, 2011a) implemented in the Mesquite (Maddison & Maddison, 2011b). Quality threshold was set high, for base trimming to 49 and for calls of ambiguity minimum secondary peak fraction for ambiguity at 0.3. Sequences were checked for F I G U R E 1 A photograph of Diaphorina citri Kuwayama female adult errors and edited manually with the program BioEdit (Hall, 1999) and then aligned in both CLUSTALX (Jeanmougin, Thompson, Gouy, Higgins, & Gibson, 1998)

| Genetic diversity and differentiation
In total, 356 COI sequences of the Asian citrus psyllids from Asia, Africa, and America were included in the analyses. Pairwise genetic distances (p-distances) were calculated using the software MEGA F I G U R E 2 Haplotype network derived from partial sequences of the COI gene from the global Diaphorina citri populations, built using PopART program. Circle size is proportional to haplotype frequency. SA, South Asia; WA, West Asia; SEA, Southeast Asia; EA, East Asia; NM, North America; SM, South America; AF, Africa F I G U R E 3 (a) Sampling regions of the Asian citrus psyllid. Detailed sampling information is presented in Table S1 in Appendix S1. A proposed scenario of early dispersal routes is illustrated with arrows. (b) Assignment of 356 individuals to K = 2 genetic clusters inferred from STRUCTURE simulations using the mitochondrial COI gene sequences. Green represents the lineage A; Red represents the lineage B v.5.05 (Tamura et al., 2011). Polymorphic sites, parsimony informative sites, and haplotype frequencies were obtained, and nucleotide (π) and haplotype (H) diversities were estimated as defined by Nei, Muruyama, and Chakraborty (1975) using the DnaSP ver. 5.10 software (Librado & Rozas, 2009). The genetic variability within the whole sample, and among lineages, was assessed by analysis of molecular variance (AMOVA) with the program Arlequin ver. 3.5.1.2 (Excoffier & Lischer, 2010). The estimation of nonsynonymous (Ka) and synonymous (Ks) substitution rates has been used as a way of understanding the evolutionary dynamics of protein-coding genes across populations (Fay & Wu, 2003;Ohta, 1995). Ks and Ka values were calculated in MEGA V.5 (Tamura et al., 2011) using the Kimura 2P model.

| Haplotype network construction and STRUCTURE simulation
A haplotype network was constructed to compare genetic relationships between geographic populations based on the genetic data. The The population substructure among the samples from the analyzed seven regions was inferred using the software STRUCTURE ver. 2.3.3 (Pritchard, Stephens, & Donnelly, 2000). We performed the analysis based on the single-nucleotide polymorphism (SNP) of COI datasets of all individuals sampled (K = 1-10) using the admixture model. Ten independent runs were evaluated for each K. The total length of the Markov Chain was set to be 100 000, and the length of burn-in period was 20 000. We adopted the parsimony criterion, in which the lowest K capturing the most differentiation among populations was selected as the most likely K value (DiLeo, Row, & Lougheed, 2010).

| Demographic analysis
To test signatures of demographic changes in D. citri, we first examined Tajima's D (Posada, 2008) and Fu's F s (Clement, Posada, & Crandall, 2000) using Arlequin ver. 3.5.1.2 under the infinite site model with 10,000 coalescent simulations. Second, population demographic changes were also examined by estimating the Harpending's raggedness index (HR) based on mismatch distribution for each population with significance assessed using Arlequin by 1,000 permutations. A significant HR value (p < .05) is taken as evidence for rejecting the sudden population expansion model (Schneider & Excoffier, 1999). Third, the demographic history of whether D. citri experienced a range expansion was investigated by the spatial expansion model in Arlequin (Schneider & Excoffier, 1999). Mismatch distributions of pairwise differences among sequences were calculated with parametric bootstrapping (1,000 replicates). According to simulations, demographically stable or admixed populations must present a multimodal distribution, whereas populations that have undergone a recent expansion generally show a unimodal distribution (Rogers & Harpending, 1992).

| Genetic variability and differentiation
We succeeded in amplifying a 760-bp COI region for the specimens from China and found only one polymorphic site, which was nonsynonymous mutation site. Two unique haplotypes were identified T A B L E 1 Region, country of collections, number of individuals, haplotypes (number of individuals), and nucleotide and haplotype diversity of Diaphorina citri in each sampled region

| Global relationships
A haplotype network analysis revealed the genetic structuring among global D. citri populations ( Figure 2 The probabilities of each individual belonging to a lineage were plotted by locality (Figure 3a), indicating geographic structuring.
The STRUCTURE analysis (Figure 3b) confirmed the presence of two lineages genetically differentiated.

| Demographic history
The  (Fu, 1997;Tajima, 1989). Furthermore, haplotype network analyses using the genetic data for D. citri further suggested rapid population expansion in that two universal haplotypes located at the central positions of two connected star-like clusters of rare haplotypes for both lineages (Figure 2). The whole sample and both lineages showed unimodal patterns of mismatch distribution curves and the strong biases toward low divergence values as distinguished by 0-and 1-nucleotide changes ( Figure 4). Furthermore, none of the Harpending's raggedness (HR) index of mismatch distribution was significant ( Figure 4). These findings illustrated a relatively recent and rapid population expansion from a small number of colonizers (Rogers & Harpending, 1992) or a range expansion with high levels of migration between neighboring demes (Excoffier, 2004;Ray, Currat, & Excoffier, 2003).

Source of variation df Sum of squares
Variance components

Percentage of variation
The whole sample

| Genetic diversity and differentiation
We present the most inclusive analysis of the mtDNA phylogeographic structure of D. citri to date, collating previously published COI data (Boykin et al., 2012;de León et al., 2011;Guidolin, Fresia, & Cônsoli, 2014;Lashkari et al., 2014) (Avise, 2000). However, our findings seem more likely to reflect colonization histories and are consistent with genetic diversity of introduced populations being low due to the small size of propagules and bottleneck effects (Nei et al., 1975). South America, we hypothesize, was colonized recently by D. citri, yet it displays high mtDNA diversity.
This could be the result of multiple invasions and/or a relative longer expansion history, consistent with prior studies (Guidolin et al., 2014;Silva et al., 1968). The genetic diversity of D. citri has been studied in Iran and Pakistan (Lashkari et al., 2014) well as between lineages. We also found strong genetic differentiation between the two lineages even though they do co-occur in some localities. These results indicate limited or no mtDNA gene flow between the lineages, and limited gene flow within lineages among geographically isolated populations. On this basis, we hypothesize that this species has limited natural (nonhuman mediated) dispersal capacity (at least females, Kobori, Nakata, Ohto, & Takasu, 2011), and therefore, local populations may be unique and amenable to sitespecific containment strategies to reduce the citrus greening disease spread. This observed geographic structure of D. citri is consistent with a process of range expansion followed by isolation of populations, as demonstrated in Excoffier, Foll, and Petit (2009) studies.
The observed significant positive selection (Ka/Ks ≫ 1) here is interesting as invasion success may occur after severe population bottleneck events or even with singly mated female invasion depending on its heterozygosity level (Zayed, Constantin, & Parker, 2007). In recent years, massive pesticides were used to control D. citri in the citrus orchards (Lopes, Frare, Yamamoto, Ayres, & Barbosa, 2007;Yamamoto et al., 2009). Therefore, the pest has been undergoing the intense selective pressure. Additionally, rapidly changing environments in invasive regions, such as climate change, also produce strong selective pressure. These directional selection can quickly lead to the loss of haplotypes and the expansion of rare haplotypes with different fitness attributes, affecting management strategies to be adopted for controlling the pest D. citri and/or the citrus greening disease.

| Dispersal history
Our mtDNA phylogeographic analysis suggests two haplotype lineages, each characterized by a common "ancestral" haplotype and several derivatives of it ( Figure 2). Although several potential processes might yield such a pattern (as outlined above), we hypothesize based on these results that the two lineages have independent colonization histories from two separate ancestral haplotypes. We hypothesize three dispersal routes across the globe, in part showing parallel dispersal patterns for the two lineages in Southeast Asia, Africa, and South America. The genealogies based on network can provide insight into ancestral haplotypes, as these haplotypes are expected to be common, widely distributed and highly connected (Chen et al., 2010). In the study by Boykin et al. (2012) probably India (Hall, 2008). Our findings are consistent with the latter; a southern and/or southeastern Asia origin. The Asian citrus psyllid is known to occur in China likely as an invasion within the last 80 years (Hoffmann, 1936). It colonized South America in the 1940s (Costa Lima, 1942). During the 1990s, D. citri invaded North America (French, Kahlke, & De Graca, 2001;Halbert & Núñez, 2004;Pluke, Qureshi, & Stansly, 2008;Tsai, Wang, & Liu, 2000).

| Future directions
The continent-global scale phylogeographic history of the psyllid D. citri has to date exclusively been studied with mtDNA represented by the COI locus. These studies have allowed understanding of the matrilineage structure of this species and generated hypotheses regarding dispersal and colonization of this important pest. However, while mtDNA data have proven highly useful in phylogeographic studies, there are several shortcomings of this approach (e.g., Avise, 2000).
Not the least of these is the potential mismatch between the dispersal history of males and females and likely mismatch between a species tree and any single gene tree. With this paper, we summarize available COI data and add new data re-analyses of global patterns, and discuss prior and novel hypotheses on dispersal and colonization history of D. citri. Robustly testing these hypotheses, and gaining a deeper understanding of population structure of this species, will require data drawn at random from the whole genome, utilizing methods such as RADseq or targeted sequence capture. Such data will tease apart male and female dispersal, testing the "female-based" dispersal hypotheses