People recognize that the impact of Pleistocene glaciation cycles on floral and faunal distributions is a major force shaping population diverging patterns in many organisms [1,2]. For phytophagous insects feeding on a limited range of host trees, distribution data for hosts define likely distributions for the insects [3]. Therefore, phytophagous insects that are widespread in southern tropical and subtropical zones are ideal models for the study of how past glaciations had affected population differentiation and speciation. For the agricultural insects, it is important to understand the geographical origin of cryptic species and invasive species, to reveal the population distribution of the pest and design the effective control measures and prevent the further diffusion [4,5].

Locust is one of the most important agricultural pest insects. Phylogeographic and population genetic studies of the migratory locust Locusta migratoria of Asia have been reported [6-8]. Zhang and Kang [7] chosen microsatellite loci to look into the genetic diversity among geographical populations of the migratory locust. Zhang et al. [8] chosen eight nuclear microsatellite loci as markers to study relationships of substructured populations in Chinese L. migratoria, suggesting the Locust populations in China should be divided into three distinct population groups : Hainan population, Tibetan population and the North China population. Ma et al. [9] explored the worldwide genetic structure and phylogeography of the locust populations based on the sequence information of complete mitochondrial genomes of 65 individuals and three mitochondrial genes of 263 individuals from 53 sampling sites; although the migratory locust can migrate over long distances, their results revealed high genetic differentiation among geographic populations. They believed that historical climatic fluctuations played a primary role. Additionally, Lovejoy and the rest have used mitochondrial DNA to explain the biogeography of the desert locust Schistocerca gregaria [10].

The yellow-spined bamboo locust, Ceracris kiangsu, is endemic to China. C. kiangsu distributes in a broad range across southern China, from Yunnan to Jiangsu [11,12]. This bamboo locust is one of the most economically important pest insects affecting Chinese bamboo forest, with outbreaks commonly causing near total mortality of bamboo within the susceptible age classes [11-14]. They mainly feed in large groups on the leaves of bamboo plants, often causing new culms to die and decrease in producing new shoots. Also, when lack of food for this species damages maize, rice corn, and so on [13,15]. However, there are no available phylogeographical data on C. kiangsu. Such data urgently need to unfold the evolutionary history and contemporary population genetic structure of this locust, and to guide developing proper management and control tactics for the bamboo locust [16-19].

In this study, we chosen mitochondrial DNA as a molecular marker to build the phylogeographical structure and phylogenetic relationships of C. kiangsu among the geographically separated populations. The reason of selecting mitochondrial ND2 gene is the third variable gene, after ATPase 8 and ND6 genes in the mitochondria [20]. And recent studies demonstrate that ND2 gene is a powerful molecular marker for phylogenetic analysis and the population structure research [21-24].

Here, we used ND2 sequence data to examine: 1) the population genetic structure of C. kiangsu among its distribution regions; 2) whether phylogeographical structure exists in the bamboo locust populations in China; 3) the important factor to affect historical population demography of this locust on native scale.

Specimen collection and DNA extraction

The specimens used in this study were collected by some collaborators during July 2008 to August 2009, from 13 locations in China, shown in Figure 1 and Table 1. The sampling scheme has almost covered the species’ distribution ranges. To avoid bias in our inferences about the phylogeography [25], each sample site was analyzed with a sufficient number of individuals (1220), except Guizhou sample site with only two individuals collected. The specimens were individually sealed in 50 mL modified centrifuge tubes, and preserved in 100% ethanol until DNA isolation.

Figure 1: Map showing the approximate locations of the sampling sites. Sampling localities and code are listed in Table 1.

Total genomic DNA was extracted from one hind leg using a standard proteinase K/phenol extraction protocol [26], and stored at -30 until ready for use.

Polymerase chain reactions were performed with primers ND2F118 (5′-CTCTCATTTATTCCC ATATTAG-3) and ND2R908 (5′- GGTTTAAGTGTCATTGATAGTG-3), designed to amplify a fragment of about 790 bp of the ND2 gene. Amplifications were performed on a DNA Engine® Peltiter Thermol Cycler (BIO-RAD) in 50 μL reaction volume composed of: 29 μL of sterilized distilled water, 5 μL of LA PCR Buffer II (Takara), 5 μL of 25 mM MgCl₂, 6 μL of dNTPs mix (2.5 mM each), l.5 μL of each primer (10 μM), 1.5 μL of DNA template and 0.5 μL (2.5 U) of TaKaRa rTaq polymerase (TaKaRa Bio Inc, Dalian, China). The PCR program was: initial denaturation at 94 for 4 min followed by 35 cycles of 30 s at 94, 30 s annealing at 50, 1 min at 72 and a subsequent 10 min final extension at 72.

PCR products were tested by electrophoresis on an agarose gel. After a single and bright target band was observed in the agarose gel, the PCR product was purified by using the V-gene PCR Clean-up purification Kit. If more than one band present, the appropriately sized PCR product was cut from the gel and extracted using a Biospin Gel Extraction Kit (Shanghai Sangon Biotech Co., Ltd.). Some specific products were ligated to the pGEM-T Easy Vector (Promega, Madison, WI, USA). Big Dye Terminator Cycle Sequencing Kits (Applied Biosystems, USA) were used for circularly sequencing the all PCR products, by using the same primers as for PCR amplification.

The ND2 sequences were checked, assembled and aligned by using the software BioEdit version 7.0.9.0 [27] and SeqMan (DNASTAR Inc., Madison, Wisconsin, USA). The assembled sequences were blasted at the National Center for Biotechnology Information (http:// www.ncbi.nlm.nih.gov/guide/) to corroborate if they were target fragments.

All sequences were aligned using MEGA 6.0 [28] with the Clustal W option and default parameters. Some basic sequence statistics were also computed using MEGA. Haplotype identification, haplotype diversity (Hd), nucleotide diversity (π) [29] and the number of polymorphic sites (S) were estimated in DnaSP version 5.00.04 [30].

We performed analyses of molecular variance (AMOVA) in ARLEQUIN 3.11 (Excoffier L. Zoological Institute, University of Berne, Switzerland) with 1,000 permutations to examine partitioning of genetic diversity within and among populations, to test for hierarchical population genetic structure in C. kiangsu. We calculated FST using haplotype frequencies only to evaluate genetic partitioning among populations.

Two tests of selective neutrality, Tajima’s [31] and Fu’s [32], and mismatch distribution analyses were performed using ARLEQUIN 3.11 to detect recent demographic signatures of population expansions [33]. To compare observed distributions with those expected under the expansion model, we calculated the sum of square deviation (SSD) and the Harpending’s raggedness index (HRI) [34].

To represent phylogeographical structure among the haplotypes, a median-joining network was generated for all haplotypes by using software Network version 4.5.1.6 (Fluxus Technology Ltd.), with the median-joining method [35].Phygenetic analyses of all haplotypes were reconstructed by using neighbour-joining (NJ) in MEGA 6.0 [28] and maximum-likelihood (ML) in PAUP* 4.0 b 10 [36]. The NJ tree was constructed with both bootstrap and interior branch support tests. ML analysis was carried out by a heuristic search of 10 random addition analyses with treebisection- reconnection (TBR) branch swapping. The HKY substitution model [37] were obtained by modeltest version 3.7 [38] based on the akaike information criterion (AIC). The confidence levels of the nodes in the NJ and ML trees were estimated using 1000 bootstrap pseudo replicates. Two individuals from the known closely related species Ceracris nigricornis [39] were chosen as outgroups.

Genetic diversity

We obtained 736 bp fragments of the ND2 gene from 211 individuals after doing alignments. We detected 43 polymorphic sites and discovered that 17 were parsimony informative in the 736 bp fragment, while 26 sites were singleton variable sites. The nucleotide frequencies of this gene were 0.371, 0.143, 0.393 and 0.093 for T, C, A and G, respectively. From all the 211 sequences analyzed, 53 different haplotypes were detected, of which 40 were unique. The most common haplotype was widespread and represented 48.3% of the sequenced individuals. Haplotype and nucleotide diversity were 0.7508% and 0.198%, respectively, for the whole samples. The genetic diversity of all samples examined varied from 0.100% (CH population) to 0.297% (HR population) for nucleotide diversity and from 0.386 (CH population) to 1.000 (HG population) for haplotype diversity. Nucleotide and haplotype diversities at each sample site, as well as the number of haplotype in each sample, were listed in Table 1.

Table 1: Sampling sites, latitude/longitude.

The results showed the lack of genetic structure among the C. kiangsu populations (FST=0.02). When we regarded each sample site as one population, the AMOVA based on haplotype frequencies revealed that variation among populations accounted for 2.23% of the total variations, while the variation within populations was 97.77% of the total variations. We divided samples of one province into one group, the analysis results also showed that almost all the differences in ND2 sequences were between individuals within a population (var=97.62%, FST=0.02), while none of the variation could be attributed to differences between populations within a group (var=-0.83%, FSC=-0.01) or between groups (var=3.21%, FCT=0.03). All population pairwise FST for ND2 were small or zero with significant P values except the 14 pairs (CH and QZ, CH and AH, CH and HR, CH and FJ,CH and TY, CH and GX, CH and JX, CH and HG, CH and ZJ, GN and HR, RA and HR, HR and SC, SC and FJ, and FJ and GX) (marked with * in Table 2).

Table 2: Pairwise ND2 FST (below diagonal) and corresponding P values (above diagonal).

Mismatch distributions for the all-individuals of ND2 consisted of a distinct unimodal curve. The other indices were Sum of Squared deviation (SSD), 0.001 (P=0.709), and Harpending's Raggedness index (HRI), 0.038 (P=0.740). Tajima’s D analysis revealed negative values (D= -2.349), and the Fu’s test of neutrality based on 1000 simulating samplings was significantly negative (FST= -28.037). All these results suggest that a recent colonization and population expansion in C. kiangsu (Figure 2) has happened. When the mismatch analyses were separately implemented for the 13 locations, expansions were found in the 10 locations (FJ, GN, RA, HR, JX, QZ, SC, CH, TY and ZJ) (the figure not presented here).

Figure 2: Distribution of the number of pairwise differences between haplotypes of ND2 of C. kiangsu . HRI-Harpending’s raggedness index.

The median-joining network displayed a star-like pattern, where the one most common haplotype is located on the star's center and the derivatives are connected to it by short branches. The most common haplotype was found in all the sampling sites, and the second common haplotype was found only in the ten sampling sites. Other haplotypes were derived directly or indirectly from these shared common haplotypes, so the haplotypes could not be clearly separated into the 13 geographical populations (Figure 3).

Figure 3: Network of ND2 haplotypes for C. kiangsu constructed in NETWORK. Each open circle represents one haplotype with size of circle proportional to abundance of each haplotype. Small open circles with a label mv represent intermediate haplotypes. Haplotypes from the same location had the same color(s).

The phylogenetic reconstruction based on ND2 data (Figure 4) also supported this result. The phylogenetic trees resulting from ML and NJ analyses exhibited some shallow genealogy (Figure 4). The phylogenetic analysis results highly supported monophyly of C. kiangsu (bootstrap value=100). Within this lineage, haplotype 3 (Hap 3) and haplotype 36 (Hap 36) from Chongqing City are sister to the remaining haplotypes, but monophyly for the remainders of the lineages was supported with lower bootstrap values. Specimens with identical geographic origin did not necessarily group together. No distinct phylogenetic groups were detected. These results were corresponding to the median-joining network of haplotypes for C. kiangsu.

Figure 4: Phylogenetic tree of ND2 gene of C. kiangsu . Values above branches represent ML bootstrap values; values below branches represent NJ bootstrap values. Bootstrap values of 50% or greater are shown.

Polymorphism and population genetic structure

The average A-T content of the ND2 gene was 0.711, which is consistent in general to be highly A+T biased in insects [40-43]. The overall nucleotide and haplotype diversities of ND2 gene in C. kiangsu were 0.198% and 0.7508, respectively, in which the haplotype diversity was similar to many other insects, but the nucleotide diversity was much lower than those reported for the many other insects [44-46]. In the present study, high haplotype diversity and low nucleotide diversity reveal that a rapid demographic expansion from a small effective population size occurred in the bamboo locust [47]. High haplotype diversity may also show multiple refuges and secondary contact of populations from different refuges [48]. The FJ, HR and TY populations had the higher nucleotide and haplotype diversities, so we assumed that these locations may be the refuges in the quaternary glaciations [48]. It is worth to pay attention to the highly genetic diversity of HG population, and we speculate that a few samples may be the most important cause.

AMOVA analyses revealed the majority of the total variation of ND2 was found in individuals within a population, and that variation among populations accounted for 2.23% of the total variation. With other analysis results (mismatch distributions and neutrality tests) for C. kiangsu, these results are the signature of no significant population genetic structure. When the populations were treated separately, AMOVA analyses of all location pairwise FST of ND2 showed small or zero with significant P values except the 14 pairs (CH and QZ, CH and AH, CH and HR,CH and FJ,CH and TY, CH and GX, CH and JX, CH and HG, CH and ZJ, GN and HR, RA and HR, HR and SC, SC and FJ, and FJ and GX). AMOVA results also showed a lack of genetic structure among populations .The low-level of overall and pairwise population differentiation values suggested that these populations have not been geographically and genetically isolated to some extent. The reason may be the total samples experienced rapid dynamic and the secondary contact of haplotypes following the expansion [48].

Although the sampling locations in this study separate at least 200 km apart, no obvious phylogeographical pattern was revealed for the bamboo locust in the all-sample-ing sites. Both phylogeny and median joint network analyses did not divide the bamboo locust haplotypes into distinct groups coinciding with any of the other populations in the all-sampling locations. Shared haplotypes were found from all sampling locations, with a lack of interpopulation genetic variation in this locust, suggesting there are no effective barriers to gene flow among the populations. The median joint networks also supported the hypothesis of a recent expansion in size from a smaller number of founders, since the common ancestral haplotype located on the centre of the star-like network, and other haplotypes derived directly or indirectly from these shared common haplotypes. Coalescence theory [49] predicts the probability that a haplotype is the oldest is equal to its frequency in the sample, and the expected rank of haplotypes by age is equal to their rank by frequency,. In the present study, the most common haplotype is identified as base by the derived haplotypes from median joint networks. The phylogenetic trees exhibited some shallow genealogy, suggesting coalescence to the most common recent ancestor. This suggestion can be interpreted because of a low long-term effective population size, which most likely result from changes in population size through time, or extinction and decolonization of subpopulations [50].

The most common haplotype of ND2 sequences were geographically widespread. Other rare haplotypes were similar to the common haplotypes. Occurring of low frequency haplotypes in the widely separated locations implies high gene flow [51]. Usually, highly migratory species are expected to have slight phylogeographical structure among their distributional ranges [52], because strong gene flow has the effect of homogenizing genetic variation over geographical populations, counteracting random drift, selection and mutation [52-55]. Although the bamboo locust is widespread in the southern China, our results also revealed that populations separated over 1,000 km in East China did not show genetic differentiation and there is no isolation by distance within this region. According to classical population genetic theory, these results indicate there are strong gene flows [7]. These results reflect that Wuyi and Nanling mountains coupled with some other mountains in southern China were not effective barriers limiting gene exchange between neighbouring populations on both sides of these mountain ranges. As a forestry and agricultural pest, the utility of pesticide will play a significant affection on the phylogeographical structure of C. kiangsu.

Mismatch distribution analyses found the observed distributions of mutation differences among haplotypes for the whole ND2 sequences as a whole fitted the expected distribution under a model of sudden population expansion (Figure 2). Tajima’s D analysis and Fu’s test of total samples of molecular marker revealed negative values (Table 1). Results showed that selective neutrality was rejected for the DNA marker. These negative values are not necessarily signatures of nonneutral molecular evolution since many authors have noted that Tajima’s test is also sensitive to other, such as recent population bottlenecks or population growth, which will drive the value of Tajima’s D towards more negative values [57,58]. Indeed, significatively negative values for Tajima’s D statistic and Fu’s test in the present study most likely resulted from recent population growth. All these results suggest there is a recent colonization and population expansion in C. kiangsu. We speculate the great fluctuations in climatic and ecological or environmental conditions during the Quaternary glaciations played an important role in the population demography.

Moderate to relative high haplotype diversities and low nucleotide diversity, the wide geographic distribution of common haplotypes and not strait geographic distribution of rare haplotypes, and the absence of phylogeographical structure based on the ND2 data, suggest that high gene flow is the most important factor infecting the genetic structure of the C. kiangsu populations. A recent and sudden population expansion of this locust is also supported by the medianjoining network. Populations of this species are most likely originated in FJ, HR and TY populations, because the populations in these regions harbor greater genetic diversities and a higher number of haplotypes. Shallow phylogenetic trees are evidences of high movement rate between neighbouring habitats, colonization and longdistance migration. Therefore, we think that changes on the climatic and ecological or environmental conditions during the Quaternary glaciations play an important role in the genetic diversity and population demography for C. kiangsu. This finding is important to expand the perspective of control measures for this bamboo locust, because it suggests that controls might be effective only if did simultaneously for the whole distribution range. In the future, we prone to sampling more specimens and selecting more rapidly evolving markers such as mitochondrial DNA control region or microsatellites, which may be useful to confirm these hypotheses.

We would like to thank the collectors for providing specimens: Yueming Cheng, Xiu-jin Chen, Yue-jin Liu, Li-xia Dai, Li-jun Dong, Zhaohui Tang, Ya-cheng Wang, Dai-lin Wang, Shu-lin Guo, Jian-wen Liu, Jian-ming Luo, Ru-ming Liu, Dan-yang Shi, Gan-jian Qi, Heng-qing Xie and Wei-ning Liu, and also thank Wen-Juan Xuan,Wan Chen, Zhou Fan and Jin-Lliang Zhao for their assistances on experiments. This work was supported by grants from the National Natural Sciences Foundation of China (No. 30970339 and No. 31572246) to Guo-Fang Jiang, Nanjing Normal University Innovative Team Project (No. 0319PM0902), A Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions, and the Guangxi Natural Science Foundation (2012GXNSFBA053056).