Introduction
The white-backed planthopper (WBPH), Sogatella furcifera (Horváth) (Hemiptera: Delphacidae) is a serious rice pest that is distributed in South, Southeast and East Asia, the South Pacific and Australia (Mun et al., Reference Mun, Song, Heong and Roderick1999). WBPH not only causes direct sucking damage, but also is the main insect vector of southern rice black-streaked dwarf virus (SRBSDV), which is a notorious member of Fijivirus (He et al., Reference He, Liu, He, Wang, Chen, Guo, Correll and Song2013). WBPH is unable to overwinter in high latitudes such as in North China, Korea and Japan, and can only undergo northwards or north-eastwards long-distance migration, following north-easterly flowing air current each year from tropical and subtropical regions (National coordinated research group for white-backed planthoppers, 1981).
Much effort has been devoted to measuring the dispersal of WBPH, including high-altitude aerial netting, recapture of color-labeled individuals, dissection of female ovaries, radar monitoring and atmospheric current analysis, etc. (National coordinated research group for white-backed planthoppers, 1981). These methods are logistically difficult to set up, and thus are often limited in space and time. WBPH populations migrating from Indo-China Peninsula are considered as the primary sources of WBPH in Northern China (Shen et al., Reference Shen, Shang and Liu2003), while WBPH populations migrating from Northern Vietnam or Southern China are considered as the primary sources of the WBPH in Korea and Japan (Matsumoto et al., Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013). However, the exact primary resource and migration routes of the WBPH are still controversial. In addition, the genetic background of different WBPH geographic populations has been unclear. Whether all the Asian populations represent one large panmictic population is under debate.
Several studies have used molecular markers (microsatellites or mitochondrial genes) to investigate the dispersal of pest insects (Loxdale & Lushai, Reference Loxdale, Lushai, Woiwod, Reynolds and Thomas2001; Llewellyn et al., Reference Llewellyn, Loxdale, Harrington, Brookes, Clark and Sunnucks2003; Endersby et al., Reference Endersby, McKechnie, Ridland and Weeks2006; Wei et al., Reference Wei, Shi, Gong, Jin, Chen and Meng2013). Molecular markers can be used in an indirect approach, such as for comparing allele frequency differences between populations or for constructing of gene trees. Another way is to use them in individual-based assignment tests, allowing individuals to be assigned to their population or origin. The individual-based assignment tests can estimate dispersal parameters directly (Broquet & Petit, Reference Broquet and Petit2009). Both methods require that the molecular markers be highly polymorphic. To date, mitochondrial COI gene is the main molecular marker used in WBPH (Mun et al., Reference Mun, Song, Heong and Roderick1999, Matsumoto et al., Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013). Using partial COI gene (823 bp), Mun et al. (Reference Mun, Song, Heong and Roderick1999) have studied the genetic variation of WBPH among five Asian countries (Korea, Philippines, China, Malaysia and Vietnam). Their result showed that the genetic differentiations among the five countries were in low level, and implied that the entire Asian populations may represent one large panmictic population. Similar to Mun et al.(Reference Mun, Song, Heong and Roderick1999), Matsumoto et al. (Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013) studied the population structure of WBPH sampled from seven Asian countries and regions (Japan, China, Taiwan, Vietnam, Laos, Thailand and Philippines) using a longer mitochondrial DNA sequence [cytochrome oxidase subunit I (cox1) – transfer RNA for Leu (trnL2) – cox2, 1, 927 bp]. However, significant genetic differentiations were detected, especially between Taiwan and southern Vietnam (pairwise F ST=0.236). This result seems not support the idea of the entire Asian populations representing one large panmictic population. Besides the COI gene, inter-simple sequence repeat (ISSR) markers have also been used for population genetic inference of WBPH (Liu et al., Reference Liu, Gui and Li2010). However, ISSR markers are increasingly questioned by researchers concerning their reproducibility, dominance and homology.
A disadvantage of using mtDNA marker is that only female movements can be inferred, and this may result in biased estimates of dispersal rates (Burg et al., Reference Burg, Lomax, Almond, Brooke and Amos2003). In addition, COI gene has been demonstrated less robust for distinguishing WBPH populations (Matsumoto et al., Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013). In recent years, microsatellites, also called simple sequence repeats (SSRs), have been widely used in studies of population genetics and dispersal because of their high levels of polymorphism, codominant Mendelian inheritance, and ease of detection by polymerase chain reaction (PCR). Although microsatellite DNA has been reported in the WBPH transcriptome (Xu et al., Reference Xu, Zhou, Zhou, Wu and Zhou2012), useable microsatellite markers have not been characterized from the WBPH transcriptome. Presently, only four microsatellite markers have been identified in WBPH and these are transferred from Laodelphax striatellus (Liu et al., Reference Liu, Zhang and Hou2014). And the four microsatellites have not been used in population structure study of WBPH. A potential problem of transferring markers is that there may be a decrease in allelic diversity and polymorphism when loci are used in non-source species (Ge et al., Reference Ge, Sun, Cui and Hong2013). The lack of microsatellite markers developed specifically for WBPH has caused a bottleneck in studies of population genetics in WBPH (Liu et al., Reference Liu, Gui and Li2010). In this study, we developed a new set of microsatellite markers for WBPH, and developed efficient multiplex PCRs for detecting them. We also demonstrate the utility of the new markers by conducting a preliminary population structure study on five geographic populations. The high polymorphism and ease of use of these loci suggest they will be powerful tools for population genetic and dispersal studies in the future.
Materials and methods
WBPH sampling and DNA extraction
One hundred and fifty-eight adult females of WBPH, representing five populations, were sampled from five provinces of China during the summer of 2010 (fig. S1). The locations, rice cropping system strategies, sampling dates and number of individuals sampled are given in table 1. These populations were sampled by randomly collecting adults from 20 rice plants in a 5×5 m square. Genomic DNA was extracted from individual adult females using the Wizard®SV Genomic DNA Purification System (Promega) according to the manufacturer's instructions.
Populations, locations, rice cropping system strategies and dates for the analyzed samples of WBPH in China and basic genetic diversity indices calculated from 19 microsatellite loci.
N, sample sizes, N A, number of alleles per locus; A R, mean allelic richness over 19 microsatellites (for standardized samples of 24 individuals); H O, mean observed heterozygosity over 19 microsatellites; H E, mean expected heterozygosity over 19 microsatellites.
Microsatellites development
Transcriptome sequences were downloaded from the Short Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) with accession numbers of SRP009194 (http://www.ncbi.nlm.nih.gov/sra) (Xu et al., Reference Xu, Zhou, Zhou, Wu and Zhou2012). These reads were assembled using TRINITY software with default parameters (Grabherr et al., Reference Grabherr, Haas, Yassour, Levin, Thompson, Amit, Adiconis, Fan, Raychowdhury and Zeng2011). Short transcripts (<350 bp) were discarded from the resulting assembly. After assembling, a total of 21,147 unigenes were obtained with a total length of 14.1 Mb. The obtained unigenes were analyzed using SciRoKo 3.4 software (Kofler et al., Reference Kofler, Schloetterer and Lelley2007) to identify reads containing simple tandem repeats. We used the Mismatched; Fixed Penalty mode with mismatch penalty value five to identify microsatellites. Only dinucleotide to pentanucleotide motifs with more than four repeats were selected for further study. A total of 1014 microsatellites-containing sequences were identified. We designed the primers using PRIMER 3 software (Rozen & Skaletsky, Reference Rozen, Skaletsky, Misener and Krawetz2000) with the following parameters, a 40–60% GC content, 100–300 bp final product size, primer T m s of 52–58 °C (optimum 55 °C) and a 3 °C maximum difference in melting temperature between paired primers. Primers were successfully designed for 663 of these microsatellite-containing sequences. We randomly selected 220 of the 663 primer pairs for PCR validation. The selected ESTs were also compared to the NCBI nr protein database using the BLASTX program (http://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastx&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome). A cut-off value of 1e-7 was chosen to assign a potential homologue for each EST sequence.
To examine the polymorphism by fluorescent genotyping method, we first amplified single locus by PCR. The DNAs of six individuals were used as the PCR templates. The forward primer of each pair was labeled with a fluorescent dye FAM or HEX (table 2). PCR amplifications were conducted with the Type-it Microsatellite PCR Kit (QiagenTM) in a 10 μl reaction volume containing 0.2 μM of each primer, 10–100 ng template DNA, 1×Type-it Multiplex PCR Master Mix (Qiagen). The thermal profile used an initial denaturation step of 95 °C for 3 min followed by 35 cycles of denaturing at 94 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 30 s. A 15 min final extension at 72 °C was added at the end of cycle. Negative controls (ddH2O) were included in PCR reactions. PCR products were analyzed using ABI 3130 sequencer (Applied Biosystems) according to the manufacturer's instructions. Allele sizes were determined using GENEMAPPER version 4.0 (Applied Biosystems), using LIZ-500 as size standard. The loci which showed polymorphism and clear peaks were selected for the subsequent multiplex PCR development.
Characteristics of 21 microsatellite loci isolated from WBPH.
Primer, primer concentration for both forward and reverse primers; N A, number of alleles per locus; H O, observed heterozygosity; H E, expected heterozygosity; F IS, inbreeding coefficient, the bold indicated significant departure from HWE after FDR correction for multiple tests (P<0.05); r, null allele frequencies calculated with the Oosterhout algorithm (van Oosterhout et al., Reference van Oosterhout, Hutchinson, Wills and Shipley 2004 ); annealing temperature of each pair of primers is at 55 °C.
Optimal combinations of the polymorphic microsatellite loci for multiplex PCR reactions were determined with MULTIPLEX MANAGER V1.0 (Holleley & Geerts, Reference Holleley and Geerts2009). The amplifications protocol for multiplex PCR was consistent with the protocol for single locus PCR, except for the number of primer pairs contained in the PCR reaction. The DNA of six individuals used in the initial polymorphism examination was also used as PCR template. To check the pairwise compatibility of the primers, the multiplex PCR profiles were compared to those produced in single locus PCR. Amplification profiles were further optimized by adjusting the relative concentration of each SSR primer (initial with equimolar concentrations: forward primer: 0.2 μM; reverse primer: 0.2 μM) to obtain a homogeneous signal across all loci. After the multiplex PCR development, these new microsatellites were subsequently characterized using 20 WBPH DNA samples randomly picked from Wenshan population.
Null allele frequencies were determined with Micro-Checker version 2.2.3 with the Oosterhout algorithm (van Oosterhout et al., Reference van Oosterhout, Hutchinson, Wills and Shipley2004). The number of alleles per locus (N A), observed heterozygosity (H O) and expected heterozygosity (H E) were assessed using GenAlEx 6.5 (Peakall & Smouse, Reference Peakall and Smouse2012). Hardy–Weinberg equilibrium (HWE) for each locus and genotypic linkage disequilibrium between pairs of microsatellites (210 pairs) were calculated with GENEPOP V 3.4 (Raymond & Rousset, Reference Raymond and Rousset1995). False discovery rate (FDR) correction was applied for 210 tests for linkage disequilibrium and 21 tests for HWE using R 3.0.2 software (R Development Core Team, 2005). This method will preclude the false positive in multiple tests. When the loci deviated from HWE, tests were performed using GENEPOP V3.4 to find out whether deviations were the result of a deficit or an excess of heterozygotes.
Population structure analyses
A total of 158 individuals sampled from five provinces of China were genotyped at the 19 new microsatellites which are free of null alleles. The population genetic diversity indices such as number of alleles (N A), observed heterozygosity (H O), expected heterozygosity (H E), and the fixation index at each locus in each population were assessed using GenAlEx 6.5 software. We used FSTAT 2.9.3.2 software (Goudet, Reference Goudet1995) to calculate the allelic richness (A R) based on a minimum population size of 24 individuals. To check loci potentially under selection, we used the program LOSITAN (http://popgen.eu/soft/lositan/) to generate 50,000 simulated loci distributions with two options ‘neutral mean F ST’ and ‘force mean F ST’.
Pairwise F ST values for each population comparison were calculated with FSTAT 2.9.3.2 software (Goudet, Reference Goudet1995) to determine the genetic differentiation of the five WBPH populations. We used STRUCTURE V2.3.4 software (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000) to infer patterns of genetic structure. The simulation was run using the model combination of admixture ancestry model with sampling locations as prior (LOCPRIOR), and the correlated allele frequency model. The Markov chain Monte Carlo simulation was run 10 times for each value of K (1–5) for 106 iterations after a burn-in period of 105. The most probable number of K in the data was determined by comparing the log likelihood (lnPr(X/K). In addition, discriminant analysis of principle components (DAPC) was also used to explore the population structure using the ADEGENET software (Jombart, Reference Jombart2008; Jombart et al., Reference Jombart, Devillard and Balloux2010). Unlike the STRUCTURE software, ADEGENET software uses the non-model-based multivariate approach (do not rely on HWE, nor do they suppose the absence of linkage disequilibrium). To identify the optimal number of clusters, the k-means clustering algorithm was employed. We calculated the Bayesian information criterion (BIC) for K=1–40, where K=number of clusters. The optimal clustering solution should correspond to the lowest BIC. This method is more suitable for populations that violated HWE and linkage equilibrium assumptions.
Results
Microsatellite characteristics
Of the 220 primer pairs, 95 primer pairs amplified the expected products. Of the 95 successful primer pairs, 21 produced clean band patterns and revealed polymorphism among the tested WBPH genotypes. For these 21 selected ESTs, homology searches of the NCBI nr database with BLASTX found three ESTs with significant hits to insect genes (table 2). The 21 microsatellites were successfully amplified in four multiplex PCR sets and binned into three runs of an automated DNA sequencer (fig. 1). The primer sequences and the concentration for each locus were given in table 2. The 21 microsatellites were subsequently characterized by 20 individuals from Wenshan population.
Examples of multiplex PCR sets and multiloading runs. All of the 21 loci present clear peaks and generally homogeneous signal.
No linkage disequilibrium between loci was detected (P>0.05 for all pairs of loci after FDR corrections). The number of alleles per locus ranged from 3 to 9, with an average of 4.8 (table 2). The observed and expected heterozygosity ranged from 0.200 (SF9) to 0.900 (SF17 and SF18) and from 0.184 (SF13) to 0.799 (SF5), respectively. The inbreeding coefficient (F IS) ranged from −0.227 (SF17) to 0.541 (SF6), with an average of 0.018. After FDR correction for multiple tests, three (SF2, SF6 and SF18) of the 21 loci significantly deviated from HWE. Further analyses revealed that SF2 and SF6 were caused by homozygote excess, and SF18 was caused by heterozygote excess. Null alleles were identified by MICRO-CHECKER 2.2.3 software for SF2 (0.213) and SF6 (0.255) (table 2).
Population structure
Table 1 shows the genetic diversity indices averaged over all loci for each WBPH population sample. The number of alleles and allelic richness for each population ranged from 4.5 (Nanning) to 5.3 (Wenshan) and 4.2 (Guiyang) to 4.8 (Longyan), respectively. The observed and expected heterozygosity ranged from 0.436 (Sanya) to 0.494 (Wenshan) and from 0.454 (Sanya) to 0.482 (Wenshan), respectively (table 1). After FDR correction for multiple tests, only one (SF18 in Wenshan) of the 95 locus-population combinations displayed a significant departure from HWE (table S1), which was caused by heterozygote excess. Details of single-locus indices of genetic diversity can be found in supplementary table S1. For the neutral test, all of the 19 loci were under the neutral expectations (fig. 2).
Comparison of F ST and H E in 19 loci to identify outliers and potential candidates for selection using LOSITAN. All of the 19 loci are under the neutral expectations.
Pairwise estimates of F ST calculated between pairs of populations were very low, ranging from 0 (Wenshan vs. Nanning, Wenshan vs. Guiyang and Guiyang vs. Nanning) to 0.0063 (Sanya vs. Longyan), with an average F ST of 0.002. Permutation tests showed that none of these pairwise F ST values were significantly different from zero (table 3). The STRUCTURE analysis result showed that the log-likelihood of the multilocus genotypic data was maximal and with a low variance for K=1 (fig. 3). The DAPC result showed that the 158 individuals fall into three clusters (fig. 4A). The scatterplot of the first two components of the discriminant analysis (DA) showed that the three clusters were set apart from each other clearly (fig. 4B). The membership probability of each individual was presented in fig. 4C. However, the cluster 3 in Wenshan had relatively higher proportions than the remaining two clusters. And the cluster 2 in Nanning and Sanya had slightly higher proportions than cluster 1 and cluster 3 (fig. 4C).
Inference of the number of genetic clusters (K) from STRUCTURE simulations for five Chinese WBPH populations. L(K), mean log likelihood over 10 runs (error bars=standard deviations). L(K) values indicate that the most probable number of clusters is one.
DAPC analytic results. (A) Scatter plots of BIC values calculated by adegenet software. BIC values indicate that the most probable number of clusters is three. (B) Plot of DAPC for three assigned genetic clusters, each indicated by different colors. Dots represent different individuals. (C) Cluster membership probabilities of each individuals based on the discriminant functions of DAPC. Each individual is represented by a vertical bar. The distributions of the three clusters in the five populations reflect strong gene flow between populations.
Pairwise F ST (lower-left matrix) and population differentiation P-values (upper-right matrix) between all populations.
Pairwise F ST values between all populations were calculated from 19 microsatellites according to Weir & Cockerham's (Reference Weir and Cockerham1984). Two hundred random permutations were used to test for significant genetic differentiation between populations.
Discussion
Microsatellite characteristics
The 21 microsatellites that we developed in this study had high degree of polymorphism. Most of the loci were in HWE, which suggested that the WBPHs were under sexual random mating. Only three loci (SF2, SF6 and SF18) significantly deviated from HWE. For SF2 and SF6, Micro-Checker suggested the presence of null alleles as the most probable cause of the deviations. For SF18, the deviation of HWE was caused by heterozygote excess. Usually, heterozygote excess is not as common as the heterozygote deficiency, and therefore has not been fully theoretically explored. Overdominant selection favoring heterozygotes, associative overdominance, and negative assortative mating are generally used to explain heterozygote excess in natural populations (Stevens et al., Reference Stevens, Salomon and Sun2007). However, heterozygote excess at SF18 was not detected in other populations when we performed population structure analyses. Based on this fact, we speculated that the heterozygote excess in Wenshan population may be resulted from the biased sampling. A potential problem with the EST-microsatellites is that they are located in coding or UTR regions, so that they are probably linked to sites under selection. Therefore, estimates of population structure and genetic diversity based on these markers may not reflect the neutral expectation. However, neutrality test demonstrated that all of the loci were under neutral expectations. These microsatellites should be powerful tools for inferring the dispersal rate from the distribution of genetic variation among populations of the WBPH. In addition, we succeeded in defining multiplex SSR sets and multi-loading combinations, which makes it possible to fingerprint each genotype at 21 loci using only four PCR runs and three electrophoresis tracks or injections. This will reduce the cost and labor for analyzing a large number of samples.
Population structure
The level of expected heterozygosity in WBPH was similar to that in the brown planthopper Nilaparvata lugens (Stål) (Sun et al., Reference Sun, Zhang, Ge and Hong2011), and was lower than that in the small brown planthopper L. striatellus (Fallén) (Liu et al., Reference Liu, Zhang and Hou2014). This may be due to that both the WBPH and N. lugens have weaker overwinter ability than the L. striatellus, most of the WBPH and N. lugens in China are immigrants from the Indo-China Peninsula. Founder effect may result in the relatively lower genetic diversity in WBPH and N. lugens than in L. striatellus. Although the rice farming systems are not exactly the same among the five sampling locations (table 1), the five WBPH populations had similar level of genetic diversity. This could be as a result of the strong gene flow, which allowed homogenizing genetic diversity among regions.
The five WBPH populations showed a high degree of genetic connectivity. This population structure pattern was consistent with those of other migratory insects (Endersby et al., Reference Endersby, McKechnie, Ridland and Weeks2006; Franklin et al., Reference Franklin, Ritland and Myers2011; Lyons et al., Reference Lyons, Pierce, Barribeau, Sternberg, Mongue and De Roode2012), and contrasted with insects that are widely distributed, but have limited dispersal ability, such as the striped riceborer Chilo suppressalis (Walker) (Meng et al., Reference Meng, Shi and Chen2008). The low level of population differentiation indicated that the high dispersal ability of this pest allowed genetic mixing between populations from remote geographical origins. This population genetic structure pattern also suggested that the five geographic populations constituted one panmictic population. However, whether all the Asian populations represent one large panmictic population is still unclear. Results of the previous studies may have been restricted by small sample size or less robust markers, and cannot give a good resolution on this issue. Further study is needed to analyze more geographic populations of Asia with the powerful microsatellite markers. The DAPC analytic result suggested that the 158 individuals fell into three clusters. In accordance with the low F ST values and the cluster result conducted by STRUCTURE software, the distributions of the three clusters in the five populations also reflect strong gene flow between populations. The reason of existing three clusters was unclear here. Whether the three clusters are associated with different ecotypes or come from three geographic resources, these speculations are needed to be carefully verified in the future studies.
Compared with the previous studies, the pairwise F ST reported here was much lower than in Matsumoto's study using mitochondrial DNA sequence. Two reasons may contribute to this inconsistent result between the two studies. First, the mtDNA is haploid and maternally inherited, and therefore it has a fourfold smaller effective population size. So, greater F ST-values are expected for mtDNA than nuclear DNA at drift-migration equilibrium because of the smaller effective population size (Birky et al., Reference Birky, Maruyama and Fuerst1983). Another reason is the different scale of geographic populations employed in the two studies. The scale of geographic populations was larger in Matsumoto's study (Matsumoto et al., Reference Matsumoto, Matsumura, Sanada-Morimura, Hirai, Sato and Noda2013) than in our study. In contrast to Matsumoto's study, Liu et al. (Reference Liu, Gui and Li2010) revealed much higher level of population differentiation (F ST=0.72) among geographic regions (11 geographic regions in Yunnan Provinces of China and three Southeast Asian countries, including Vietnam, Laos and Myanmar) using ISSR markers. This unexpected high F ST value may be caused by the dominant character of ISSR markers. Dominant molecular markers tend to overestimate genetic differentiation among populations (Aagaard et al., Reference Aagaard, Krutovskii and Strauss1998).
In conclusion, we developed 19 high-quality microsatellite markers for the WBPH. The multiplex PCR sets for the 19 microsatellites will facilitate large-scale population genetic studies of WBPH with less cost and labor intensity. Preliminary population genetic structure analyses of WBPH individuals showed that the entire population comprised three multi-locus genotypes spreading across the five provinces sampled. These microsatellites also provide a new approach to inferring the dispersal rates of this pest.
The supplementary material for this article can be found at http://www.journals.cambridge.org/BER
Acknowledgements
We thank Jia-Fei Ju and Xiang-Fei Zhang of the Department of Entomology, Nanjing Agricultural University (NJAU) for help with the collection of WBPH. We are grateful to Dr. Gao Hu of NJAU for making the collection map. This work was supported by a grant-in-aid from the National Key Basic Research Program (973 Program, No. 2009CB119200) from the Ministry of Science and Technology of China, a grant-in-aid (No. 200903051) from the Science and Technology Research Program of the National Agricultural Public Welfare Fund from the Ministry of Agriculture of China and a grant-in-aid from the Specialized Research Fund for the Doctoral Program of Higher Education (SRFDP) from the Ministry of Education of China (Priority Development Area, 20110097130005).
