The white-backed planthopper, Sogatella furcifera (Horváth) (Hemiptera, Delphacidae), is a long-range migrant and is a major rice insect pest in Asia (Kisimoto & Sogawa, Reference Kisimoto, Sogawa, Drake and Gatehouse1995; Sogawa, Reference Sogawa, Heong, Cheng and Esclada2014). As Chinese hybrid rice varieties have spread throughout China and Vietnam, outbreaks of S. furcifera have occurred subsequently in these regions. Furthermore, S. furcifera has become problematic in Japan (Sogawa, Reference Sogawa, Heong, Cheng and Esclada2014). In temperate Asian countries, including Korea, Japan, and northern China, S. furcifera cannot overwinter and is thought to migrate from northern Vietnam or southern China in June and July (Matsumoto et al., Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013). S. furcifera overwinters in northern Vietnam, and migrates to southern China from April to early May (Otuka et al., Reference Otuka, Matsumura, Watanabe and Ding2008; Hu et al., Reference Hu, Lu, Tuan, Liu, Xie, McInerney and Zhai2017), and is also able to overwinter in southern China (Hu et al., Reference Hu, Xie and Wang1988, Reference Hu, Liu, Fu, Huang, Wang, Liu, Lü and Ye2015). Therefore, populations of S. furcifera in southern China are likely to be mixed. It is unknown whether northern Vietnam populations are homogeneous or mixed. The sources of S. furcifera populations in temperate Asian countries, in particular Korea and Japan, remain unclear.
To elucidate the source regions of S. furcifera in these regions, the population structure of S. furcifera has been examined with various molecular markers, which include cytochrome c oxidase subunit I (COI), inter-simple sequence repeats (ISSR), and microsatellites. Mun et al. (Reference Mun, Song, Heong and Roderick1999) found that COI sequence variation of S. furcifera was low in Korea, China, Malaysia, Philippines, and Vietnam. Matsumoto et al. (Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013) used cox1–trnL2-cox2 data to examine S. furcifera populations from ten Asian regions (Japan, China, Taiwan, northern and southern Vietnam, Laos, Thailand, northern and southern Philippines, and Papua New Guinea) and found little genetic differentiation among these Asian populations. In contrast, Liu et al. (Reference Liu, Gui and Li2010) found a high level of genetic variation among populations of S. furcifera in China, Laos, Myanmar, and Vietnam by ISSR analysis. The S. furcifera populations were grouped into two clusters; one from northern Vietnam and Laos to southeastern and central Yunnan, and the other from Myanmar to southwestern Yunnan. Sun et al. (Reference Sun, Jiang, Wang and Hong2014) examined genetic variation of five S. furcifera populations 1314 km apart in southern China using 19 microsatellite markers. The level of population differentiation was very low, forming one genetic cluster.
In general, population genetic studies on S. furcifera in Asian countries (Mun et al., Reference Mun, Song, Heong and Roderick1999; Matsumoto et al., Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013; Sun et al., Reference Sun, Jiang, Wang and Hong2014; Qiu et al., Reference Qiu, Jiao, Hu, Liu, Huang and Li2016) show high levels of genetic connectivity and strong gene flow in S. furcifera, which makes it difficult to identify the exact sources of S. furcifera in temperate Asian countries. Northern Vietnam or southern China are considered the primary sources for the Far East Asian populations. While there is a possibility that all the Asian populations of S. furcifera may be one panmictic population, it is still very important to more specifically elucidate migration patterns, gene flow, and genetic connectivity among various geographic populations of S. furcifera in Asia because characterization of biotypes, such as insecticide resistance and susceptibility to resistant rice varieties, is crucial for planthopper management (Matsumoto et al., Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013).
In the present study, we sampled S. furcifera from 43 locations in seven countries in Asia (Korea, Bangladesh, China, Laos, Nepal, Thailand, and Vietnam). Also, using the microsatellites described by Nam et al. (Reference Nam, Coates, Kim, Park and Lee2015), we characterized the genetic differentiation and gene flow among these populations, and identified possible migration routes of S. furcifera to Korea.
Materials and methods
S. furcifera was collected in 2012, 2013, and 2014 from various countries in Asia (Supplementary Table S1). These comprised one site each in Laos, Korea, Nepal, Thailand, and Vietnam, four sites in Bangladesh in 2012; one site each in China, Nepal, and Thailand, two sites in Bangladesh, and 15 sites in Korea in 2013; and four sites in China and ten sites in Korea in 2014. Samples from Bangladesh, Vietnam, Thailand, Nepal, and Laos were collected as part of the Asian Food and Agriculture Cooperation Initiative and other samples were collected by us. All the samples were collected using sweeping net and put in vials with 95% ethanol and stored at −20°C until DNA extraction.
DNA was extracted from individuals without abdomen using the AccuPrep DNA Extraction Kit (Bioneer, Korea). For genotyping, we used 12 microsatellite loci developed for S. furcifera by Nam et al. (Reference Nam, Coates, Kim, Park and Lee2015). Polymerase chain reaction (PCR) was performed in two separate multiplex groups. Multiplex group 1 included WBPH_T5, WBPH_T7, WBPH_T9, WBPH_T11, WBPH_T13, and WBPH_T16 markers. Multiplex group 2 included WBPH_T3, WBPH_T4, WBPH_T7, WBPH_T8, WBPH_T15, and WBPH_T18 markers. We used Takara Taq ™ (TaKaRa, Japan) in a total volume of 10 µl, with 4.9 µl distilled water, 1.0 µl 10 × PCR buffer, 1.0 µl 10 mM dNTP mixture, 0.5 µl of each primer (final concentration, 10 pmol µl−1), 0.1 µl of Taq polymerase, and 2.0 µl template DNA. PCR was conducted using an initial denaturation of 4 min at 94°C, followed by 35 cycles of 94°C for 30 s, annealing at 61°C for 30 s, 72°C for 40 s, and a final extension at 72°C for 15 min.
Genetic variation and genetic structure
The microsatellite genotype data of S. furcifera were analyzed by Micro-Checker (Van Oosterhout et al., Reference Van Oosterhout, Hutchinson, Wills and Shipley2004) for the presence of null alleles, short allele dominance, and the appearance of non-functional gene regions. The Microsatellite Toolkit (Park, Reference Park2001) was used to estimate the mean number of alleles (N A) per locus, and the observed (H O) and expected (H E) heterozygosities. The inbreeding coefficient at each locus, F IS, was calculated using FSTAT v. 2.9.3 (Goudet, Reference Goudet2001). Deviations from the Hardy–Weinberg equilibrium (HWE) were tested by Genepop v. 4.2.1 (Raymond & Rousset, Reference Raymond and Rousset1995).
The population structure of S. furcifera was estimated by the Bayesian clustering procedure using STRUCTURE 2.3.3 (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000). As different runs can produce different likelihood values, ten runs were conducted in order to quantify the amount of variation of the most likely number for each K. The range of possible clusters (K) tested was set from one to ten, with ten iterations. The lengths of the Markov Chain Monte Carlo iterations and burn-in were set to 100,000 and 200,000, respectively. For the K-value, STRUCTURE estimates the maximal value of the log likelihood (ln Pr(X/K)) of the posterior probability of each K-value (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000). We calculated the mean posterior probability for each K value and the second-order rate of change in the log probability of the data between successive values of △K, the true number of populations, was calculated by using the program structure harvester (Earl & vonHoldt, Reference Earl and vonHoldt2012).
We conducted an analysis of molecular variance (AMOVA) for the hierarchical partitioning of genetic variation among and within populations. Using the Genalex program (Peakall & Smouse, Reference Peakall and Smouse2006), we conducted a principal coordinate analysis (PCoA) and then plotted the scatter diagram based on factor scores along the two PCoA revealing most variation. The PCoA analysis was conducted separately for each year.
Gene flow measures
The historical rates of gene flow between populations (N em) were calculated according to the relationship N em = (1 − F ST)/4F ST (Wright, Reference Wright1931). N em infers the effective number of migrants per generation, N e is the effective population size, and m is the migrant rate. Pairwise estimates of genetic differentiation (F ST) between populations were calculated using FSTAT v. 2.9.3 (Goudet, Reference Goudet2001). Moreover, FreeNA (Chapuis & Estoup, Reference Chapuis and Estoup2007) was used to estimate F ST, considering null alleles, and the results of F ST adjusted for null alleles and F ST assuming no null alleles were compared.
We plotted F ST/(1 − F ST) values against geographic distances between all pairs of sample locations to examine isolation by distance (IBD) using the Mantel test and the Genalex program. By using AMOVA for populations and individuals, the hierarchical partitioning of genetic variation was analyzed, which provided an estimate of the proportion of genetic variation within and between populations.
To detect genetic signatures of dispersal and immigration of S. furcifera (Rannala & Mountain, Reference Rannala and Mountain1997), we conducted population assignment/exclusion tests using the Geneclass2 program (Piry et al., Reference Piry, Alapetite, Cornuet, Paetkau, Baudouin and Esop2004). The direct assignment test (without P-value) calculates the proportion of individuals correctly allocated to the most probable population of origin, even if the true population of origin is not clear among the reference populations. The exclusion test calculates the likelihood of a genotype occurring in the population, which is compared with the distribution of likelihoods of simulated genotypes for each reference population. If the likelihood (α) of an individual is below a predetermined threshold (e.g., α = 0.01), the population is excluded as the possible origin of the individual (Cornuet et al., Reference Cornuet, Piry, Luikart, Estoup and Solignac1999). The exclusion method does not presume that the true population of origin has been sampled because this method does not focus on the specific task of immigrant identification and each population is treated independently (Cornuet et al., Reference Cornuet, Piry, Luikart, Estoup and Solignac1999). In the exclusion test, the frequency possibilities of multilocus genotypes in each reference population were determined using Monte Carlo simulations of 10,000 individuals for the population (Paetkau et al., Reference Paetkau, Slade, Burden and Estoup2004). We followed a Bayesian statistical approach (Rannala & Mountain, Reference Rannala and Mountain1997) using the Monte Carlo resampling method.
To infer the migration pathway of individuals between populations, we applied the detection of first-generation migrants in Geneclass2 (Piry et al., Reference Piry, Alapetite, Cornuet, Paetkau, Baudouin and Esop2004). We computed the statistical criterion (the ratio L home/L max) for likelihood estimation of migrant detection (L) (Paetkau et al., Reference Paetkau, Slade, Burden and Estoup2004). The assignment criterion values of the simulated independent individuals were computed using a simulation of 10,000 independent individuals at probability thresholds of α = 0.05 and α = 0.01.
A total of 118 alleles were identified across the 12 microsatellite loci for 1644 S. furcifera individuals from the 43 locations in Asia. Allelic richness varied from 4.855 to 8.482, and H O and H E ranged from 0.275 to 0.615 and from 0.506 to 0.790, respectively (table 1). All 43 populations had positive F IS values across loci.
Number of alleles, expected heterozygosity (H E), observed heterozygosity (H O), inbreeding coefficient (F IS) and probability (P-value) of being in HWE.
1 Hardy–Weinberg exact test (Raymond & Rousset, Reference Raymond and Rousset1995) with Bonferroni correction (P = 0.000123).
The genetic differentiation between each pair of populations (corrected pairwise F ST) and the effective number of migrants exchanged per generation (N em) are shown in Supplementary Tables S2–S4. Uncorrected estimates of pairwise F ST values ranged from 0.0528 for the Bangladesh (B3) and Laos (LA) populations (ENA corrected F ST = 0.0555; B3 and LA populations) to 0.2477 for the Vietnam (V) and Thailand (TH) populations (ENA corrected F ST = 0.2513; V and TH populations) (Supplementary Table S2). Uncorrected estimates of pairwise F ST values ranged from 0.0400 for the Korea (CN) and Korea (MY) populations (ENA corrected F ST = 0.0434; CN and MY populations) to 0.2558 for the Korea (WD) and Korea (CG) populations (ENA corrected F ST = 0.2489; WD and CG populations) (Supplementary Table S3). Uncorrected estimates of pairwise F ST values ranged from 0.0088 for the Korea (SA3) and Korea (CW) populations (ENA corrected F ST = 0.0114; SA3 and CW populations) to 0.3492 for the Korea (BA) and Korea (MY) populations (ENA corrected F ST = 0.3222; WD and CG populations) (Supplementary Table S4). Both estimates of F ST (F ST adjusted for null alleles and F ST assuming no null allele results) were similar. Global estimates of F ST across all loci and all populations were high and significant (uncorrected F ST = 0.1762, 95% confidence interval (CI) = 0.1387–0.2137; ENA corrected F ST = 0.1610, 95% CI = 0.1285–0.1944). In 2012, the N em calculated from the uncorrected F ST ranged from 0.759 (Vietnam (V) and Thailand (TH)) to 4.4848 (Laos (LA) and Bangladesh (B3)), and in 2013, N em ranged from 0.7273 (Korea (CG) and Korea (WD)) to 6.000 (Korea (MY) and Korea (CN)). In 2014, N em ranged from 0.8017 (China (CH4) and Korea (MY)) to 28.2530 (Korea (CW) to Korea (SA3)). The range of N em implied an intermediate to high level of gene flow.
AMOVA among S. furcifera populations from 43 locations in Asia manifested that most genetic variation was partitioned among populations and individuals within populations (table 2). Approximately 76% of the total genetic variation was accounted for by the individuals within populations, and 24% of the total genetic variation was among populations. Significant correlation was found between genetic and geographical distance among the populations using the Mantel test of IBD (r 2 = 0.4585, P = 0.01), indicating the role of geographic isolation in the genetic structure of S. furcifera (fig. 1).
Results of PCoA analysis are shown in fig. 2. Total variance explained by the first two factors in 2012, 2013, and 2014 were 48% (25.8% for axis 1 and 22.2% for axis 2), 37% (21.5% for axis 1 and 15.8% for axis 2), and 64% (40.9% for axis 1 and 23.1% for axis 2), respectively. In 2012, the population from Thailand (TH) was conspicuously divergent from the others, and in 2014, the population from Milyangshi (MY) in Korea was divergent from the others. In 2013, dispersion of all the population was observed.
Bayesian clustering indicated three clusters (K = 3); structure output revealed a maximum value of ΔK of 79.54 (Supplementary Fig. S1). The results of a Bayesian cluster analysis of multilocus microsatellite genotypes in 2012, 2013, and 2014 are displayed as pie graphs in figs 3–5, respectively.
Assignment/exclusion test and investigation of first-generation immigrants
The percentage of S. furcifera individuals sampled from each population assigned to and excluded from each reference population, and the mean assignment log-likelihood for each potential donor population in each year, are presented in Supplementary Tables S5–S7. Populations from most locations included members whose possible origins in other populations could be excluded with ≥99% certainty. In Supplementary Tables S5–S7, individuals from all populations could not be assigned to their own population, since the percentage was too low. Among the rest of the population, individuals from the BA (Korea) population in 2013 showed the highest result of percentage, which could be assigned to their own population with 26.3%. In the exclusion test, B1 (2012) population could be excluded with 100% and MY (2014) population could be excluded with 94.7–100% certainty (0.01 threshold) as a putative origin of all populations. The mean assignment likelihood for individuals in 2012 showed that the highest value of the Nepal (Ne) population came from the Laos (La) population (mean assignment log-likelihood = −18.02), while the highest assignment likelihood of the Laos population came from the Bangladesh (B3) (−21.48). Also, the highest assignment likelihood of individuals of the B3 came from the La (−19.75) (Supplementary Table S5). In the mean assignment likelihood for individuals in 2013, the Korea population of JC&CW, WD&JD, BA&CG, GR&KP, and Bangladesh (B5) and Nepal (NE) population showed connection between sites and they were most likely from the same source (Supplementary Table S6). In 2014, the mean assignment likelihood of individuals of the SA3 (Korea) population (except from itself) came from the CH1 (China) population (−14.26), the population of SA3 and CH1 was most likely from the same source. Moreover, there were connections between TA (Korea) population and CH3 (China), NYJ (Korea) population and CH2 (China), in which China was expected to be a possible source of Korea population (Supplementary Table S7).
The number of first-generation migrant individuals for each year is summarized in Supplementary Tables S8–S10. We considered the L home/L max ratio as for the detection of first-generation immigrant between sample locations. In 2012, a total of 63 and 30 individuals were observed as possible first-generation immigrants at thresholds of α = 0.05 and α = 0.01, respectively. In 2013, 184 and 82 individuals were identified, respectively, and 81 and 39 individuals were verified, respectively, in 2014 as potential first-generation immigrants. The probable dispersal pathway among populations is displayed in figs 6–8, and was based on the first-generation immigrants at a threshold α = 0.01. We eliminated the geographical barrier and reverse movement from Korea. We speculate that there may be several migration pathways of S. furcifera to Korea. One putative route is from Nepal (NE) to Northern Vietnam (V) and from Northern Vietnam directly to Southwestern Korea (SA) (e.g., in 2012). In fig. 6, various migration pathways were revealed among Vietnam, Nepal, and Bangladesh, and these relations may sustain the result of direct dispersal pathway from Nepal to Southwestern Korea. The second putative route is from Southern China to Korea (e.g., in 2013 and 2014). It would be also possible that S. furcifera would migrate from Northern Vietnam to Korea via Southern China. Moreover, there were frequent migrations of S. furcifera among Bangladesh, Laos, and Nepal. It would be plausible that high dispersal rate may show in Southeast and Southwest Asia. In figs 7 and 8, various dispersal pathways are indicated from China to Korea, with a definite genetic connection among them. We verified the movement in Korea individually for 2013 and 2014, high rate of dispersion occurred and movement to adjacent locations appeared among the sites in Korea (Supplementary Fig. S2).
In this study, we investigated the genetic diversity and connectivity of S. furcifera populations in various locations in Asia. Low level of expected heterozygosity and high level of gene flow (N em) were observed in these regions. S. furcifera has been speculated to migrate from the Indo-China Peninsula, in particular northern Vietnam to China (Sun et al., Reference Sun, Jiang, Wang and Hong2014), and its migration to Korea and Japan is more likely from southern and/or south-eastern China (Mun et al., Reference Mun, Song, Heong and Roderick1999; Otuka et al., Reference Otuka, Matsumura, Watanabe and Ding2008). This long-range migration may cause the low genetic differentiation of S. furcifera (Sun et al., Reference Sun, Jiang, Wang and Hong2014). Several studies have addressed the low genetic differentiation among migratory insect populations over long distances (Malausa et al., Reference Malausa, Dalecky, Ponsard, Audiot, Streiff, Chaval and Bourguet2007; Krumm et al., Reference Krumm, Hunt, Skoda, Hein, Lee, Clark and Foster2008). Migration route and origins of S. furcifera populations in Korea have been based mainly on circumstantial evidence (see Mun et al., Reference Mun, Song, Heong and Roderick1999), but continuous efforts have been made applying molecular markers to verify the migration route and origins of S. furcifera populations in East Asia since Mun et al. (Reference Mun, Song, Heong and Roderick1999). Matsumoto et al. (Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013) utilized mitochondrial COI sequence and reported no genetic structure among ten regions, supporting the hypothesis that northern Vietnam area is considered as the primary source region of S. furcifera to East Asia, which includes Japan, Korea, and Northern China (Sogawa, Reference Sogawa1995; Otuka et al., Reference Otuka, Matsumura, Watanabe and Ding2008). However, mitochondrial sequences did not demonstrate much differentiation of local populations and it was difficult to distinguish regional populations of S. furcifera based on the difference in mitochondrial sequences. In another study, Sun et al. (Reference Sun, Jiang, Wang and Hong2014) applied 21 microsatellites to S. furcifera samples collected from five locations in China. A high degree of gene flow and genetic connectivity was evident with a low value of F ST, in which five geographic populations, distanced by 414.1–1314.2 km, in China represented a single panmictic population. In contrast, Liu et al. (Reference Liu, Gui and Li2010) showed a higher level of genetic differentiation (F ST = 0.72) among countries in Asia (11 geographic regions in Yunnan Provinces of China, Vietnam, Laos, and Myanmar) using ISSR markers. Most of the S. furcifera in Yunnan migrated from the adjacent southern and southwestern countries based on PCoA analysis.
In the genetic structure analysis based on the structure harvester, we revealed three genetic clusters in S. furcifera populations in Asia. The relationship among populations in Asia consisted of genetic connectivity and genetic divergence (figs 3–5). In the genetic cluster, there was genetic connectivity between Vietnam and Korea in 2012 and between China and Korea in 2013 and 2014. It is possible that S. furcifera migrates from Northern Vietnam to Korea via Southern China. In the mean assignment log-likelihood for each possible donor population in 2013, the highest value was between SA3 (Korea) and CH1 (China), and were most likely from the same source. Moreover, the genetic connection showed between TA (Korea) and CH3 (China), NYJ (Korea), and CH2 (China) in 2014, in which China might be the possible source of Korea population. Assignment testing and first-generation immigrants (L home/L max) revealed the possible dispersal pathway of S. furcifera in Asia. The analysis indicated that the source of S. furcifera in Korea may be China and Vietnam. However, the exact origin of the population could not be verified since there were too few samples from various countries and there was too much gene flow among the populations. We could assume that S. furcifera migrates from northern Vietnam to Korea via southern China. Furthermore, the winter distribution of S. furcifera may be considered a key factor affecting the spatial and temporal change of population density (Hu et al., Reference Hu, Liu, Fu, Huang, Wang, Liu, Lü and Ye2015). Several studies indicated that S. furcifera overwinters in Yunnan (Liu et al., Reference Liu, Yang, Lin and Kong1991; Tao & Sogawa, Reference Tao and Sogawa2000) and the temperature in Yunnan has been increasing since 1987 due to global climate change, with expansion of northern tropical and southern subtropical areas in China (Cheng et al., Reference Cheng, Song, Huang and Xu2014). The geographical range of S. furcifera overwintering population may increase in the future and migratory populations of S. furcifera may expand in Korea.
In conclusion, genetic diversity and connectivity among the S. furcifera populations in Asia based on 12 microsatellite loci indicated that, although the number of samples taken across the Asian countries in the same year was limited, there were recurrent movements of S. furcifera among adjacent locations. The genetic structure of S. furcifera populations in Korea was diverse, indicating that Korean populations have originated from diverse sources. Our results reveal that southern China and northern Vietnam are most likely the main source of S. furcifera populations in Korea. However, some other regions, such as Nepal and Bangladesh, might be potential sources via interconnection with Vietnam populations. Further studies are needed to examine the possibility of direct flight of S. furcifera from Nepal and Bangladesh to Korea. Also, future studies should examine more geographic populations in Asia collected in the same year with powerful analytical methods, such as single nucleotide polymorphism markers. These markers have higher genotyping efficiency, data quality, and level of genetic diversity compared with other molecular markers (Morin et al., Reference Morin, Luikart and Wayne2004).
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485318000755.
This research was supported by grants from Rural Development Administration in Korea (PJ00922907) and was partially supported by the Brain Korea Plus.