Genetic variation and phylogeographic structure of Spodoptera exigua in western China based on mitochondrial DNA and microsatellite markers.

The beet armyworm, Spodoptera exigua, is a significant agricultural pest of numerous crops and has caused serious economic losses in China. To effectively control this pest, we analyzed its genetic variation, population genetic structure and demographic history. We used mitochondrial DNA (mtDNA) fragments of the cytochrome oxidase subunit I (COI) and eight nuclear microsatellite loci to investigate genetic diversity and population genetic structure of S. exigua populations at 14 sampling sites in western China. Both mtDNA and microsatellite data indicated low levels of genetic diversity among all populations. A moderate genetic differentiation among some S. exigua populations was detected. Neighbor-joining dendrograms, STRUCTURE, and principal coordinate analysis (PCoA) revealed two genetically distinct groups: the KEL group and the remaining population group. Isolation by distance (IBD) results showed a weak significant correlation between geographic distance and genetic differentiation. Haplotype networks, neutrality testing, and mismatch distribution analysis indicated that the beet armyworm experienced a recent rapid expansion without a recent genetic bottleneck in western China. Thus, the results of this population genetic study can help with the development of strategies for managing this highly migratory pest.


Introduction
The beet armyworm, Spodoptera exigua (Lepidoptera: Noctuidae), is a significant polyphagous pest on vegetables, maize, cotton, soybeans, and ornamental plants [1,2]. Generally, S. exigua larvae feed on host leaves, decreasing crop yields and leading to the death of plants. Originating from South Asia, this species is already widely distributed in the tropical and temperate regions of Europe, Africa, North America, and Asia [3]. In China, S. exigua was first recorded in Beijing in the 1890s. It has a wide distribution in the main crop-producing areas of China and has caused severe economic losses in recent years. For example, S. exigua has spread to some provinces in North and East China and has infested a total area of over 2.7 million hectares [4]. This pest especially damages welsh onions in Northern China and has infested more than 8000 hectares in Tianjin, reducing the annual welsh onion production by 30% [5,6]. Amplification products were purified and then sequenced by Sangon Biotech (Shanghai) Co., Ltd. Chromatograms, including sense and antisense, were edited and assembled using DNASTAR 5.0 (DNASTAR, Madison, WI, USA) to obtain a single consensus sequence. Microsatellite alleles were analyzed using GeneMapper 4.0 software (Applied Biosystems).  Table. The figure is modified from the graphic of Song et al. [23]. n in the pie chart represents the sample size of each population. https://doi.org/10.1371/journal.pone.0233133.g001

PLOS ONE
Phylogeographic structure of Spodoptera exigua

Microsatellite amplification and genotyping
In this study, individuals were genotyped using eight microsatellite loci, and technical details were provided by Kim et al. [25]. These microsatellite loci were assigned unique fluorophores for the fluorescent tagging of DNA in a PCR.

Data analyses
Mitochondrial data. The mtDNA COI sequences were aligned manually using the CLUS-TAL X 1.81 program [26] and were subsequently checked manually. We used MEGA 5.2 to analyze the sequence variation information and genetic distance between DNA sequences based on Kimura's two-parameter method [27]. For each population, haplotype diversity (π) and nucleotide diversity (h) were estimated using DnaSP 4.0 software [28]. The median-joining network of the mtDNA COI haplotypes was conducted in NETWORK 2.0 [29]. We calculated population genetic differentiation (F ST ) statistics between populations using Arlequin 3.0 [30]. Analysis of molecular variance (AMOVA) based on F ST and haplotype frequencies were completed using Arlequin 3.0. Isolation by distance (IBD) was determined by testing the correlation between geographical distance and pairwise F ST /1-F ST using the Mantel test with Arlequin 3.0 [30]. Demographic history changes were analyzed for S. exigua using two neutrality tests, Tajima's D [31] and Fu's F S [32], across the 14 geographic populations. Finally, mismatch distributions were calculated between the observed and expected mismatch distributions by Arlequin 3.0 [29].
Microsatellite data. The software Micro-Checker 2.2.3 [33] was used to test for errors and null alleles. Genotypic linkage disequilibrium (LD) and Hardy-Weinberg equilibrium (HWE) were computed using GenePop 3.4 [34]. Corrections for multiple tests were performed via Bonferroni corrections. Genetic diversity indexes, such as the mean number of alleles (Na), the effective number of alleles (Ne), Shannon's information index (I), observed heterozygosity (H O ), expected heterozygosity (He), and unbiased expected heterozygosity (uHe) were determined using GenAlEx 6.41 [35]. Allelic richness (A R ), the fixation index (F ST ), and the inbreeding coefficient (F IS ) were calculated with FSTAT 2.9.3.2 [36]. Polymorphism information content (PIC) was calculated in CERVUS 2.0 [37]. The levels of genetic differentiation between pairs of populations were also estimated in Arlequin 3.0 [30].
We used POPTREE 2 to construct an unrooted phylogenetic tree by the neighbor-joining method with 1,000 bootstrap replicates [38]. The software STRUCTURE 2.3 [39] was used to identify clusters of genetically similar populations based on the Bayesian approach. Markov chain Monte Carlo (MCMC) values were set for a burn-in period of 30,000 and a run length of 105 iterations under an admixture model with correlated allele frequencies within populations. To identify the optimal number of groups (K), measures of change in the likelihood function (ΔK) and F ST (ΔF ST ) were calculated using the package CORRSIEVE [40], and the plateau criterion was applied [41]. The software CLUMPP 1.1 [42] was used to conduct model averaging of individual ancestry coefficients across 10 independent runs. Next, clusters were visualized using DISTRUCT 1.1 [43]. PCoA, based on the covariance of the genetic distance matrix, was implemented in GenAlex 6.41 [35]. The significance of the hierarchical partitioning of the genetic structure among the geographic groups was tested using an AMOVA, which was conducted in Arlequin 3.0 [30]. We used the program Bottleneck 1.2.02 to investigate whether any of the populations experienced genetic bottlenecks [44]. We performed the test using the infinite allele (IAM), strict stepwise mutation model (SMM), and two-phase (TPM) models of microsatellite evolution. We assessed the significance of the tests using Wilcoxon's test, which is the most appropriate test in cases where fewer than 20 microsatellite loci are used [45]. Asymmetric pairwise measures of a recent gene flow were estimated using a genetic approach, assignment tests, and detection first-generation migrants for each population implemented in GENECLASS 2.0 software. This software was also used to identify first-generation migrant individuals of S. exigua using a partial Bayesian method [46].

Genetic diversity
In this study, 312 individuals were collected from 14 populations in western China and genotyped with eight microsatellite loci. The null allele frequency per locus was mostly less than 0.1 (S2 Table). The Hardy-Weinberg expectations (HWE) test showed that 57 of the 112 locuspopulation combinations had significantly deviated, and no significant linkage disequilibrium was found. The average F ST , which was not used, and the excluding null alleles (ENAs) per the locus correction method, was 0.184 and 0.171, respectively, and these F ST estimations did not differ significantly (P = 0.462). Therefore, the presence of null alleles did not influence F ST estimation (S3 Table). Modern genetic variation among eight microsatellite loci of S. exigua in western China was found (S4 Table). Overall, the middle level of genetic diversity of all of the populations in the sampled regions was obtained, and the mean observed heterozygosity was equal to the mean expected heterozygosity (S5 Table). The estimates of microsatellite genetic variation differed among these populations. For example, the observed number of alleles (Na) across microsatellite loci ranged from 4.250 in Xilengrad (NMXM) to 8.875 in Kunming (KM). The effective number of alleles (Ne) ranged from 2.908 in Haikou (DLH) to 4.166 in Guiyang (GY). Shannon's information index (I) across loci ranged from 1.130 to 1.545. The inbreeding coefficients (F IS ) of the eight microsatellite loci in all 14 populations of S. exigua were nearly zero (F IS = 0.131). Kunming (KM) population had the highest genetic diversity, whereas the DLH population had the lowest genetic diversity.
A 576 bp fragment of the mtDNA COI gene sequence was sequenced for 291 S. exigua individuals (deposited in GenBank under accession nos. KY069907-KY069921). Of the 576 total characters, 566 were conserved sites and 10 were variable sites. We observed 10 haplotypes in 291 individuals. The average number of haplotypes in each population was 2.286 ± 0.914 and ranged from 1 to 4, and the Korla (KEL) population had the highest number of haplotypes. Hap 2 was the most common haplotype, present in all populations, and was shared by 247 individuals. Hap 4 was the second most frequent haplotype, being shared by 20 individuals from two populations. Genetic diversity indices haplotype diversity (h) and nucleotide diversity (π) are shown in S6 Table. On average, h and π were 0.273 and 0.00069, respectively. The highest genetic diversity (h = 0.489, π = 0.00085) was detected in Taiyuan (GSTY). In contrast, the genetic diversity value of Xilengrad (NMXM) and Hanoi (HN) was zero because no variable sites were found (S6 Table).

Median-joining network and genetic differentiation
These 10 mtDNA COI haplotypes were from 14 locations, which covered more than 4,000 km. These results indicated the random distribution of mtDNA COI haplotypes in S. exigua, which may reflect a high level of gene flow in this species. Furthermore, the median-joining network of mtDNA COI haplotypes indicated no apparent geographical clustering of haplotypes ( Fig  1). A star-shaped network suggested that many haplotypes differed from Hap 2 by only one or two mutations, which suggested that S. exigua experienced a population expansion event.
A global F ST value of 0.217 was obtained based on the microsatellite data from all populations. The pairwise F ST values between populations ranged from −0.053 to 0.484, and the exact tests detected that 68 of the 91 pairwise comparisons differed significantly (Table 1). Overall, a moderate level of genetic differentiation among populations was determined. We found significant differentiation between most populations, such as between Korla (KEL), Delingha (DLH), Kunming (KM), Chengdu (SC), Multan (BM), Faisalabad (BF), and Hanoi (HN). We also found significant differentiation between Yinchuan (YINC), Guiyang (GY), and other populations. Based on the mitochondrial DNA data, the pairwise F ST values computed from ranged from -0.20 to 0.94, with an average F ST of 0.109 (Table 1). We identified a significant genetic differentiation in 24 of the 105 comparisons (P < 0.05). Although a significant difference existed between most populations including Korla (KEL) and Lanzhou (GSTY), we found no significant differentiation among the other remaining populations.
Population genetic structure POPTREE analysis based on microsatellite data. Based on microsatellite data, the unrooted neighbor-joining tree, including the 14 S. exigua populations defined two major clades (Fig 2). One clade corresponded to the Korla (KEL) population, which originated from Xinjiang, and the second clade was composed of the remaining 13 populations in western China. Due to the absence of a strong phylogeographic structure in the latter clade, frequent gene flow in the vast western region of China is suggested.
Bayesian clustering. Based on microsatellite data, we used a clustering algorithm in STRUCTURE 2.3.3 to infer the relationships between the 14 S. exigua populations in western China (Fig 3). The mean LnP(D) values increased slowly starting from K = 2, likely representing the most appropriate number of major clusters. The maximum ΔK was reached at Principal coordinate analysis (PCoA). The population-based PCoA was performed based on a Nei's genetic distance matrix using allele frequencies of the eight microsatellite markers in the 14 S. exigua populations (Fig 4). The first and second axes account for 33.46% and 55.85% of the total variance, respectively. The above neighbor-joining tree and Bayesian clustering, based on these data collected from 312 individuals, indicated two distinct groups, which is supported by the efficacy of the PCoA approach.
Analysis of molecular variance (AMOVA). The global AMOVA indicated that greater genetic variance was partitioned within individuals (microsatellite data: 78.31%, mtDNA COI: 62.46%), and the remaining genetic variation (microsatellite data: 21.69%, mtDNA COI: 10.85%) could be explained by the variation among populations ( Table 2). The hierarchical analysis demonstrated that only 28.28% and 6.87% of the variation was explained by the variation among two groups based on microsatellite data and mtDNA COI data, respectively, whereas the most variation (microsatellite data: 62.46%, mtDNA COI: 88.77%) was explained by variation within populations. Consistent with the results of the NJ tree, PCoA and STRUC-TURE, the AMOVA analyses were also supported one genetic group.
Isolation By Distance (IBD). A Mantel test for IBD was conducted to estimate the correlation between the genetic distance and geographic distance of S. exigua. If the dispersal of S. exigua is limited by distance, genetic and geographical distances should be positively correlated, producing a pattern of isolation by distance. Applying the Mantel test, the weak significant positive correlation between genetic and geographic distance was observed among any of the S. exigua populations in western China (microsatellites: R 2 = 0.134, P = 0.010, Fig 5A;

PLOS ONE
Phylogeographic structure of Spodoptera exigua mtDNA: R 2 = 0.050, P = 0.030, Fig 5B). In contrast to the mtDNA data, the microsatellite data were more in line with the IBD pattern, perhaps as a result of the limited mtDNA gene sequence variable information and the relatively slower evolution rate.
Demographic history. Based on the mtDNA COI data, neutrality tests were conducted using Tajima's D [31] and Fu's F S [32] statistics. Significantly negative values were obtained (Tajima's D = -1.678, P < 0.05, Fu's F S = -7.987) (S6 Table). The widespread single haplotype (Hap 2) and the star-like haplotype network map provide strong evidence of expansion (Fig  1). The unimodal mismatch distribution of all populations also indicated a sudden demographic expansion in these populations (S1 Fig). The bottleneck effect analysis based on the microsatellite data showed no significant heterozygote excess in 13 of the 14 populations (except for the DLH population) under three models (IAM, TPM, and SMM), which is consistent with a normal L-shaped distribution of allelic  Table. https://doi.org/10.1371/journal.pone.0233133.g003 frequency. Therefore, these results indicated that S. exigua had expanded demographically without experiencing a severe population bottleneck in western China (S7 Table).
Individual assignment and migrant detection. In this study, we detected a high individual assignment success, 191 of 312 individuals (61.0%) were definitively classified at a P value of 0.05 based on the partial Bayesian method. Our results showed 130 migrant individuals from all of the 312 individuals, and the majority of putative migrant genotypes (110/135) were restricted to three populations of southwestern China, including Kunming (KM), Guiyang (GY) and Chengdu (SC). It also provides a highly directional from southwest China to other western regions populations. In addition, we also observed 11 individuals as potentially firstgeneration migrants (F0), of which 5 individuals migrated from three southwestern populations to almost all the other populations of S. exigua (S8 Table).  Table. https://doi.org/10.1371/journal.pone.0233133.g004

Genetic variation and genetic differentiation
In the present study, S. exigua was found to have low haplotype diversity and low nucleotide diversity based on the mtDNA COI gene (Hd = 0.273 ± 0.034, Pi = 0.00069 ± 0.00010), and a moderate expected heterozygosity (He) of 0.605 based on the microsatellite loci, demonstrating a low level of genetic diversity in S. exigua in western China. This result is consistent with previous studies [17,20,21]. Similarly, several studies suggested a low level of genetic diversity d.f., degree of freedom; � P < 0.05, ��� P < 0.001: significance level. Two groups including eastern region group and western region group based on STRUCTURE analysis (See Fig 3).
https://doi.org/10.1371/journal.pone.0233133.t002 in other migratory Lepidoptera, such as monarch butterflies (Danaus plexippus) [47]. In general, species that are capable of dispersal show little genetic differentiation between populations. S. exigua is a significant agriculture pest, but little information has been published about its dispersal ability [48]. Overall, the low genetic differentiation among most populations implies that the relatively active dispersal capacity might result in increased gene flow. However, distinct ecological climates in different regions might also contribute to the limitation of gene flow among regions. In the present study, we found a significant genetic differentiation between Korla (KEL) and all remaining populations based on the microsatellite and mitochondrial DNA data. This was caused by the long-term independent diversification of these two lineages, along with other factors such as geographic barriers and temperature limitations, and these factors may play important roles in maintaining the present phylogeographic patterns. This result is consistent with findings regarding the diamondback moth, Plutella xylostella, where sampling locations throughout temperate and subtropical Chinese regions were based on the climate division in China [49].

Population genetic structure
Knowledge of the population genetic structure can provide insight into the evolutionary and ecological processes of species. In the present study, we found no distinct distribution pattern of this pest based on mtDNA haplotypes phylogenetic tree. In line with that, the mtDNA haplotypes network analysis showed that the dominant haplotype (Hap2) was evenly distributed in the western regions of China, indicated that gene flow occurred among the populations of this species in recent history and blurred the original relationship between different populations. These genetic features are evidence of migration, as demonstrated in other migratory species [47,[50][51][52]. POPTREE, PCoA, STRUCTURE, and AMOVA analyses based on microsatellite data divided all individuals into two clusters, the first cluster including the Korla (KEL) population, and the second clade composed of the remaining 13 populations in western China. To the best of our knowledge, geographic barriers and climatic factors might contribute to the isolation of the studied population and they might play an important role in preventing gene flow among groups [53,54]. For example, mountains shaped the population structure of the terrestrial species Locusta migratoria [55] and Paeonia rockii [56]. We speculate that the differing climates and large mountain systems in western China might have prevented or slowed nuclear gene flow among the western populations of S. exigua, resulting in lineage differentiation. For example, Mt. Tianshan might have acted as a substantial geographical barrier, producing significant population differentiation between KEL and the remaining populations. Similarly, genetic structuring on a local scale has also been reported in the migratory moth Helicoverpa armigera in Australia, where it is thought that drought conditions impacted its migratory capacity and resulting genetic structure [57]. We found weak significant positive correlation between genetic and geographic distance among all of the S. exigua populations in western China, and this result is consistent with those reported in previous studies [17,21,22]. IBD affects are pronounced in moderately mobile species or are weak in both low-mobility and high-mobility species since extensive gene flow homogenizes populations across the large scale of geographic regions [58].

Pest management implications
S. exigua is a polyphagous species that has been described as feeding on more than 300 plants, which implies its considerable adaptive potential. This pest originates from South Asia. It is currently widely distributed in most of the main crop-producing areas in China. We collected five individuals of S. exigua using sex pheromone traps for the first time in Yichun, Heilongjiang, in 2012. Therefore, this species has not accumulated additional genetic variation due to its short history of expansion in Northern China. In the present study, we found a high level of genetic diversity in the Kunming (KM) population. This high level of diversity can probably be explained by S. exigua not experiencing serious founder effects, genetic bottlenecks, or strong selection (e.g., from insecticides) in this region. Resistance to insecticides in insects is an example of evolutionary adaptation to environmental changes [59]. Since the past several decades, cultural and chemical controls have been used to prevent the spread and damage of S. exigua. S. exigua has developed resistance to some insecticides. For example, this species has a long history of exposure to carbamate pesticides and has developed a high level of resistance to methomyl in California, USA, since 1989 [9]. Therefore, quantifying the level of resistance to insecticides and investigating the distribution of resistance genes relative to the genetic structure and the level of gene flow among Chinese S. exigua populations are essential. Understanding the dispersal ability, genetic structure, and population demography of this pest is critical for determining both theoretical aspects of its evolution and the effective implementation of pest forecasting systems. In the future, we will explore population genetic differentiation and population genetic structure based on the genomic level of S. exigua, and deeply reveal the evolution relationship, reconstruct population history, it will help us to gain insight into the adaptation to climates and ecology factors at the genomic level. These works will be carried out before a reasonable control strategy is developed for managing S. exigua.

Conclusions
We used mitochondrial DNA (mtDNA) fragments of the cytochrome oxidase subunit I (COI) and eight nuclear microsatellite loci to investigate genetic diversity and population genetic structure of S. exigua populations at 14 sampling sites in Western China. Low levels of genetic diversity and a moderate genetic differentiation among S. exigua populations were detected. Population genetic structure analysis revealed two genetically distinct groups: the KEL group and the remaining population group. Isolation by distance (IBD) results showed no significant correlation between geographic distance and genetic differentiation. Population demographic history analysis indicated that the beet armyworm experienced a recent rapid expansion without a recent genetic bottleneck in western China.