Extensive gene flow of white-backed planthopper in the Greater Mekong Subregion as revealed by microsatellite markers

The white-backed planthopper (WBPH), Sogatella furcifera (Horváth), is a destructive pest of rice in the Greater Mekong Subregion (GMS) countries including Cambodia, Laos, Myanmar, Thailand, Vietnam, and China’s Yunnan Province. Our previous study not only confirmed the immigration sources of the WBPH in China’s Yunnan Province were from Myanmar, Vietnam, and Laos, but also indicated that Cambodia was likely an additional migration source. To further clarify the migration sources and patterns of the WBPH in the GMS, we investigated the genetic structure of 42 WBPH populations using microsatellite loci markers. The analysis of genetic diversity, heterozygosity deficit, and heterozygosity excess based on the nuclear markers suggest that there is extensive gene flow between the 42 sampled populations from the GMS. The genetic structure confirmed the immigration sources of WBPH as revealed by mitochondrial markers and trajectory analyses methods in previous studies. These findings will aid in the sustainable regional management of this insect pest in the GMS.

also indicated that Cambodia was a likely additional migration source. As nuclear genetic markers, microsatellite loci have been widely used in elucidating the genetic structure of insect populations, because they are inherited codominantly and have a broad distribution and high abundance throughout the genome [15][16][17][18] . Determining the genetic structure of the WBPH in the GMS based on nuclear markers will provide further insights into the gene flow and migration patterns of the WBPH in this region.
In the present study, we investigated the genetic structure of 42 WBPH populations using nuclear (microsatellite loci) markers, to reveal the gene flow and migration patterns of the WBPH in the GMS. These results will benefit future sustainable management programs of this insect pest in the GMS. Analyses of genetic structure within populations. The estimator of the fixation index, Fis, was significantly different in 20 of the 42 populations, demonstrating that the presence of sub-structure within the populations was common (Table 1). In testing for deviation from mutation-drift equilibrium in BOTTLENECK, we detected a significant heterozygosity deficit (Wilcoxon test P < 0.05) in only three populations (CY, JP, and V3). The significant heterozygosity deficit in the three populations may result from demographic expansion 18,19 because there were no significant departures from Hardy-Weinberg equilibrium (Fis) in these populations ( Table 1), suggesting that the significant deviation from mutation-drift equilibrium was not due to sub-structure (the Wahlund effect) within these localities. In testing the deviation from mutation-drift equilibrium in BOTTLENECK software, we did not detect a significant heterozygosity excess in any population under the TPM or SMM models, although under the IAM model, a significant heterozygosity excess (Wilcoxon test P < 0.05) was detected in twelve of the populations ( Table 2), indicating that these twelve populations might have experienced a genetic bottleneck.
Analyses of genetic structure among populations. When considering each pairwise Fst, 512 of 861 (59.5%) Fst values were associated with a significant exact test (Table 3). Analyses using STRUCTURE software identified two genetic clusters overall (K = 2) ( Fig. 1): one cluster consisted mainly of individuals from the two populations in Thailand (T1 and T2), and a few individuals from three of the populations from Cambodia (C1, C2, and C3); another cluster consisted mainly of all individuals from the other 37 populations, a few of individuals from the two populations in Thailand (T1 and T2), and most individuals from the three populations in Cambodia (C1, C2, and C3). When K = 3, the individuals from the two populations from Thailand (T1 and T2), and a few of individuals from three of the Cambodian populations (C1, C2, and C3) could also be differentiated from the other individuals.

Sum of squares
Variance components  Table 4. Analysis of molecular variance (AMOVA) for structures of Sogatella furcifera collections. The Mantel test results produced an r value of 0.0807 for microsatellite alleles (P = 0.8810) (Fig. 2), indicating that no correlations were found between genetic distance and geographical distance among the populations of the WBPH in the GMS countries, indicating that extensive gene flow exists among these WBPH populations.

Percentage of variation
Significant genetic structure of the WBPH was observed at two hierarchical levels (among populations and within populations) ( Table 4). Most of the variation was at the within populations level (91.53%). Although the variation among populations (8.47%, P < 0.01) was small, it was significant. These results demonstrated that the variations of genetic differentiation in the WBPH are mainly from inter-populations.
Gene flow based on microsatellite data. Based on microsatellite data, the average values of the numbers of migrants in the different countries were similar except in Thailand (Table 5). In China's Yunnan Province, the average number of migrants was the highest in southern Yunnan, while the lowest numbers were found in the western region; this is similar to previously published results based on mitochondrial COI data 9

Discussion
Evidence for extensive gene flow in the WBPH within the GMS. This study showed that the level of the genetic diversity in most populations originating from different countries was similar, suggesting that extensive gene flow occurs between the WBPH populations within the GMS. The heterozygosity deficit is used to test for population expansion, whereas the heterozygosity excess test is used to provide evidence of a genetic bottleneck. In this study, only three populations from China's Yunnan Province had a significant heterozygosity deficit. The results in testing the heterozygosity excess were completely inconsistent under the TPM, SMM, and IAM models. Di Rienzo et al. showed that most microsatellites fit TPM better than SMM or IAM 20 . Based on the TPM model, there is no significant heterozygosity excess in any of the tested populations, which suggests that no severe bottleneck effects exist in the GMS populations. Our study indicate that bottleneck effects have not played an essential role during the genetic differentiation of the WBPH. This may be due to the bottleneck effects on heterozygosity being transient and observable for only a few generations 21 .

Migration sources of the WBPH within the GMS. The extensive gene flow between WBPH populations
within the GMS is consistent with our previous study 9 showing that there are multiple immigration sources of the WBPH in China's Yunnan Province including Myanmar, Vietnam, Laos, and Cambodia. Although the populations from central Thailand (T1 and T2 populations) were shown to have an extensive gene flow in this study, both had limited gene flow with neighboring countries, limiting the probability of immigration of these populations into China's Yunnan Province. These results support previous results based on the trajectory analyses methods 12,22,23 that the populations in central Thailand would be incapable of immigrating to Yunnan due to the lack of sufficient wind currents, incorrect wind direction and the excessive distance involved, although gene flow with the populations in Vietnam, Laos, and Cambodia would be possible 23 .
The Cambodian populations (C1-C3 populations, especially C2) do have a limited gene flow with those from central Thailand but have extensive gene flows with other populations from other countries in the GMS. Our previous study showed that the specific mtCOI haplotype from Cambodia is only found outside the country in China's Yunnan Province (BS and GM populations). The extensive nuclear gene flow also substantiates the probable occurrence. However, the trajectory analyses methods demonstrated that the emigrant population from Cambodia would not be able to migrate to Yunnan Province 23 . The immigration of the WBPH individuals from Cambodia into China's Yunnan may be indirect from Vietnam or Laos, which have extensive gene flow with those in Cambodia. Whether the immigration of WBPH from Cambodia into China's Yunnan is direct or/and indirect remains to be determined.

Future research on the source population of S. furcifera within the GMS. The immigration sources
and patterns of the WBPH as demonstrated by mitochondrial and nuclear markers are helpful in devising future, sustainable, regional management programs for this important pest in the GMS. However, the genetic basis for the migration should be further explored, including the influence of such factors as wing polyphenism 24,25 . For Table 5. Numbers of effective migrants per generation (N e m) of Sogatella furciferain the Greater Mekong Subregion. : mutation-scaled population size, which is effective population size × mutation rate per site per generation; M: mutation-scaled immigration rate, which is the immigration rate divided by the mutation rate. example, Xu et al. showed that two insulin receptors in the migratory brown planthopper (Nilaparvata lugens) (Stål) (Hemiptera: Delphacidae) play an important role in controlling long versus short wing development 25 , providing the first evidence of a molecular basis for the regulation of wing polyphenism in insects. We had elucidated the genetic diversity of the WBPH in the GMS countries based on mtCOI and SSR markers. In future studies, it will be necessary to analyze the reliability and significance of these molecular markers relative to their consistency with WBPH biological data. With the exception of the molecular markers, the genome difference and transcriptome analyses also should be considered in a followup study. Although the 42 populations that were collected from the GMS countries help to explain the genetic diversity within somewhat limited areas, additional populations from different regions and from different seasonal occurrence in the GMS should be collected and analyzed in future studies. Additional attention needs to also be paid to more widely distributed populations, such as those from Malaysia, Indonesia, the Philippines, Bangladesh, Pakistan, India and other known occurrences of S. furcifera within Asia and outside of Asia to further explore and clarify the source population of S. furcifera in the GMS.

Conclusions
Based on the nuclear (microsatellite) markers, the analysis of the genetic diversity, heterozygosity deficit, and heterozygosity excess suggested that there is extensive gene flow between the WBPH populations in the GMS. The genetic structure confirmed the immigration sources of the WBPH as revealed by mitochondrial markers. There is a certain gene flow between the populations in Thailand and Cambodia. It should be further explored whether the immigration of WBPH from Cambodia into China's Yunnan Province is direct or/and indirect. These results will be helpful to the sustainable regional management of this insect pest in the GMS.

Materials and Methods
Field sampling and DNA extraction. Adult Table 7) were screened from 40 pairs of newly designed primers based on the WBPH microsatellite sequences in GenBank (until November 11, 2014) and were then used to amplify the loci using WBPH DNA as the template. The primers and the annealing temperature are described in Table 2. The PCR reactions were performed in 20 μL of a solution containing 2 μL 10 × buffer, 1.5 mM MgCl 2 , 0.2 μM dNTPs, 1 unit Taq DNA polymerase, 2 μL template DNA, and 0.2 μM of each primer. PCR amplification was carried out as follows: initial denaturation at 94 °C for 4 min, followed by 35 cycles of 30 s at 94 °C, 90 s at the primer-specific annealing temperature (Table 1) and 60 s at 72 °C, and a final elongation step at 72 °C for 30 min. The products were run on an ABI 3730xl DNA analyzer (Sangon, ShangHai, China) and the allele size was determined by comparing the mobility of the PCR products to the GeneScan ™ 400HD size standard using GeneMapper software version 3.2 (Applied Biosystems, ShangHai, China).
Based on the microsatellite alleles, the average number of alleles per locus (Na), the effective number of alleles (Ne), the observed heterozygosity (Ho), the expected heterozygosity (He), and Nei's expected heterozygosity (Nei) of each of the 42 WBPH populations were calculated using POPGENE v.1.31 26 . The estimator of the fixation index, Fis, was performed to detect deviation from neutrality using GENEPOP v.4.2 27 . Wilcoxon test P value for heterozygosity deficit compared to expectations at mutation-drift equilibrium (Pwil) was calculated using ARLEQUIN v.3.5 software 28 .
Analyses of genetic structure within populations based on microsatellite data. Deviation of the mutation-drift equilibrium in each population was tested using the BOTTLENECK software 18 . The heterozygosity deficit was evaluated using the Wilcoxon test under the two-phase mutation model (TPM) recommended for microsatellite data 20  Analyses of genetic structure among populations based on microsatellite data. The traditional population differentiation approach, Weir and Cockerham's estimator of the fixation index Fst 29 , was calculated using GENEPOP v.3.4 software 27 . The correlation between genetic differentiation and geographic distance was examined by Mantel test using IBDWS v.3.15 software 28 . The distribution of genetic variation was investigated by the analysis of molecular variance (AMOVA) using ARLEQUIN v.3.5 software 28 , and by calculating allelic diversity, heterozygosity, and pairwise values of Fst among the 42 populations. The genetic clustering of samples were examined using STRUCTURE v.2.3.2 software 30 , using the Bayesian clustering approach with a burn-in period of 50,000 iterations and one million Markov chain Monte Carlo (MCMC) repetitions under the admixture ancestry model. Twenty independent runs were performed for each testing K value, ranging from K = 1 to 42, and ΔK was used to calculate the optimal number of genetic clusters (K) 31 .  Table 7. Sequence of microsatellite primers designed in this study.
populations in the GMS, the effective numbers of migrants per generation N e m was calculated using microsatellited data respectively. N e m is M ( = N e μ, where μ is the mutation rate per site per generation; M = m/μ, where m is the migration rate) calculated using Bayesian search strategies in MIGRATE v. 3.2.16 32 .