Email updates

Keep up to date with the latest news and content from BMC Evolutionary Biology and BioMed Central.

Open Access Research article

Phylogeography and Demographic History of Chinese Black-Spotted Frog Populations (Pelophylax nigromaculata): Evidence for Independent Refugia Expansion and Secondary Contact

Hua Zhang, Jie Yan, Guoqiang Zhang and Kaiya Zhou*

Author Affiliations

Jiangsu Key Laboratory for Biodiversity and Biotechnology, College of Life Sciences, Nanjing Normal University, Nanjing 210046, China

For all author emails, please log on.

BMC Evolutionary Biology 2008, 8:21  doi:10.1186/1471-2148-8-21


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2148/8/21


Received:6 August 2007
Accepted:24 January 2008
Published:24 January 2008

© 2008 Zhang et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Pleistocene glaciations had considerable impact on phylogeographic patterns within and among closely related species of many vertebrates. Compared to Europe and North America, research on the phylogeography of vertebrates in East Asia, particularly in China, remains limited. The black-spotted frog (Pelophylax nigromaculata) is a widespread species in East Asia. The wide distribution of this species in China makes it an ideal model for the study of palaeoclimatic effects on vertebrates in East Asia. Our previous studies of P. nigromaculata revealed significant subdivisions between the northeast China populations and populations in other regions of the mainland. In the present study, we aim to see whether the deepest splits among lineages and perhaps subsequent genealogical divisions are temporally consistent with a Pleistocene origin and whether clade geographic distributions, with insight into expansion patterns, are similarly spatially consistent with this model.

Results

Using 1143 nucleotides of the mitochondrial cytochrome b gene from 262 individuals sampled from 28 localities, two main clades (clade A and clade B) differing by c. 7.72% sequence divergence were defined from parsimony analyses. The corresponding timing of lineage divergence, 0.92 Mya, indicates a most likely Pleistocene split. The A clade is further subdivided into two sub-clades, A1 and A2 with 1.22% sequence divergence. Nested clade phylogeographical and population demographic analyses suggested that the current distribution of this frog species was the result of range expansion from two independent refugia during the last interglacial period. We discovered a population within which haplotype lineages A and B of P. nigromaculata coexist in the Dongliao area of China by nucleotide sequences, PCR-RFLP and ISSR (inter simple sequence repeat) patterns. The ISSR result in particular supported divergence between the mitochondrial clades A and B and implied introgressive gene flow between the two divergent lineages.

Conclusion

Nested clade phylogeographical and population demographic analyses indicate that the current distribution of P. nigromaculata is the result of range expansion from two independent refugia during the last interglacial period in late Pleistocene. One refugium was in east China and the lower elevations of south-western plateau. The distribution of the other mitochondrial clade is consistent with the presence of a refugium in the Korean Peninsula. The gene flow as detected by ISSR markers suggests a range expansion of the two refugia and a secondary contact between the two highly divergent lineages in the Dongliao (DL) area of northeast China.

Background

During Pleistocene glaciations, Northwest Europe, Siberia and the northern most regions in North America were covered with ice sheets [1,2]. Pleistocene glaciations had extensive impact on phylogeographic patterns within and among closely related species of many vertebrates [3]. In Europe and North America, various species have been studied to determine the climatic and geological effects on phylogeographic patterns and population structures. In Europe and America, species dispersed to southern locations to survive in refugia and then expanded northward again during interglacial periods when climate recovered to or exceeded current mean global temperatures. Although there are discussion and debates regarding the degree of Pleistocene effects on speciation in North American birds [4-6], Pleistocene conditions appear to have played an active role both in initiating intraspecific and species-level divergences of the most closely related North American taxa [7]. However, climate recovery in East Asia during interglacial periods did not seem to resemble that of Europe or North America primarily due to the continuing uplift of the Tibet Plateau during late Tertiary. In China, the glacial advance was not as extensive as in Europe and North America because of the monsoons in East Asia. Biotic zones of Asia were also located at higher northern latitude than on other continents [8-10]. Despite the relatively mild climate, species' distribution were affected by climatic fluctuations during Pleistocene. However, research on the phylogeography of vertebrates in East Asia, especially in China, remains limited both in terms of the number of model species and geological timescales.

The black-spotted frog (Pelophylax nigromaculata) is a widespread species found in East Asia covering both Palaearctic and Oriental realms. The wide distribution of this species in Northeast, North, East, Central and Southwest China renders it an ideal model for investigating the palaeoclimatic effects on vertebrates in East Asia. P. nigromaculata is variable in size and color pattern over its wide geographic range; but shows no distinct morphological differentiation among populations and no subspecies have been designated. Traditionally, the black-spotted frog was treated as a member of genus Rana, but was recently placed into Pelophylax, a genus reintroduced by Fei et al. [11]. Recent studies on the phylogeny of ranid frogs based on mitochondrial and nuclear DNA sequences consistently supported Pelophylax as a well resolved monophyletic unit [12-15].

Studies of molecular phylogeography of P. nigromaculata have been independently reported by two groups. Yang et al. proposed that during the last interglacial period, P. nigromaculata experienced a rapid population expansion and that its distribution range fluctuated latitudinally in response to climatic oscillations [16]. But details of the evolutionary history and population structure were not inferred because of the small number of sampling localities (eight) and the relatively few Cyt b sequence data (704 bp). In a previous study, we examined the genetic structure of this frog in mainland China using 5'control region sequences and detected significant subdivision between Jilin-Liaoning (Northeast China) populations and other mainland populations [17].

In this report, we tested the hypothesis that the Quaternary glaciations were the causes of this significant genetic structure. We analyzed the complete mitochondrial Cyt b sequences and generated ISSR fingerprint data in order to reconstruct the phylogeographic patterns in black-spotted frog populations sampled from 28 localities across the Chinese range (Fig. 1, Table 1). These data allowed us to assess the evolutionary history (demography, age and origins) of various Chinese populations and to address the hypothesis that the observed phylogeographic divisions are consistent with past range fragmentation by Pleistocene glaciations.

thumbnailFigure 1. Sampling sites and mitochondrial haplotype distributions of P. nigromaculata populations. The numbers indicate local populations (see Table 1 for details). KR and JP stand for haplotypes from Korea and Japan, respectively. Lower case letters "a" to "d" indicate distinct but geographically proximate locales within the Dongliao area.

Table 1. Collection data for P. nigromaculata specimens included in this study.

Results

Cytochrome b variation

Of the entire 1143 bp of Cyt b sequences, 197 nucleotides are variable and 130 are parsimony-informative. Nine, 29, and 159 variable sites are in the first, second, and third codon positions, respectively. No indels or premature stop codons were observed. This result, together with a strong bias against guanine (mean G = 13.3%, A = 25.8%, T = 27.4% and C = 33.4%), implies that the target fragment is mitochondrial Cyt b rather than its nuclear homologue. We identified 118 haplotypes among the 216 P. nigromaculata complete Cyt b sequences [GenBank: AY803813~AY803895 and DQ006233~DQ006267]. Eighty-seven of these were unique haplotypes, 19 were shared within local populations and 12 were shared among local populations (Table 2). The sequence distances (based on the HKY model) vary from 0.09 to 8.65% (mean 2.2%) among ingroup taxa, distances among different haplotype clades vary from 1.10% to 7.73% (Table 3). The overall haplotype diversity (h) and nucleotide diversity (π) were 0.975 and 0.028, respectively.

Table 2. Distribution of Cyt b haplotypes shared in different local populations of P. nigromaculata.

Cyt b phylogenetics

Maximum parsimony (MP) analysis from 122 ingroup haplotypes resulted in 6884 most parsimonious trees, each 566 steps in length (CI = 0.7621, RI = 0.9027). Two major monophyletic haplotype groups were revealed in the Chinese population of P. nigromaculata: designated here as clade A and clade B (as shown in Fig. 2a). Clade A included all haplotypes present in most sampling regions of this frog species: from northeast China (locality 24, 47°20'N) to southwest China (locality 16, 25°08'N); clade B was found limited to eastern part of northeast China (locality 25, 26, 27 and 28). Clade A was further subdivided into two sub-clades: A1, with a broad geographical distribution and weak phylogenetic structure; A2, with higher bootstrap support of 85% and distributed to mountainous regions of southwest China (locality 17–23, see Fig. 1). When included the Korean and Japanese haplotypes, this phylogenetic topology showed that the single Korean haplotype and clade B clustered with 100% bootstrap support, and that the three Japanese haplotypes formed a sister relationship with clade A and showed 100% bootstrap value (Fig. 2).

thumbnailFigure 2. Phylogenetic trees based on Cyt b haplotypes. (a) Fifty percent majority-rule consensus of 6884 parsimonious trees based on 122 Cyt b unique haplotypes of P. nigromaculata (CI = 0.7621, RI = 0.9027). Numbers above branches are bootstrap values (1000 replicates). Letters or numbers in parenthesis indicate the population abbreviation and the quantity of individuals sharing the haplotype within each local population. The asterisk followed by an integer indicates the number of individuals sharing the haplotype among different populations (see Table 2 for details). (b) The concordant topology of maximum likelihood, maximum parsimony and neighbor-joining trees resulting from the reduced 49 exemplar ingroup taxa. Numbers above branches are ML bootstrap values (100 replicates); numbers below branches are MP/NJ bootstrap support values (1000 and 5000 replicates).

Maximum parsimony (MP), maximum likelihood (ML) and neighbor-joining (NJ) analysis of our reduced dataset resulted in concordant topologies (topology and bootstrap support value of major clades were shown in Fig. 2b). ML analysis resulted in a tree with a likelihood score of -ln L = 3396.27; MP analysis resulted in 3372 most parsimonious trees, each 528 steps in length (CI = 0.7670, RI = 0.8388). The most obvious feature of the three topologies was consistent with those revealed by MP tree based on 122 haplotypes (Fig. 2a). The sequence distances among these haplotype clades also yielded obvious insight of such phylogenetic relationships (Table 3). The division pattern of haplotypes into several clades (A1, A2, B and Japan) was concordant in three analyses, but support for sub-clade A2 was lower in ML (65%) than in MP (80%) and NJ (88%) analyses. Lower support in ML compared to the MP and NJ analyses was likely the result of a more complex model of sequence evolution used for ML reconstruction.

Table 3. Cyt b sequence distances and the standard deviations among major lineages under the HKY model.

For all the sequence data, Fu and Li's D* and F* were negative and not significant, and Tajima's D = -1.22, which was not significantly different from 0 (P > 0.10). Patterns of variation within our mitochondrial data are thus consistent with the neutrality hypothesis.

PCR-RFLP analysis

Most haplotypes are unique to local populations. Haplotypes shared across local populations (12 instances, Table 2) occur entirely in geographically adjacent areas except for one case (haplotype 3), which is shared between southwest China (BS) and central China (HC). To our surprise, locality 25 (northeast China) is composed of haplotypes grouping of the two highly divergent lineages (A and B). This observation suggests that the Dongliao area (locality 25, DL) may represent a secondary contact zone of the two deeply split lineages. Accordingly, we widened our sampling scales in this region and collected samples from four contiguous sites of 25-a, 25-b, 25-c, and 25-d with the maximum pairwise distances about 28 km (Fig. 1). We used 56 individuals in the PCR-RFLP analysis on Cyt b. When the PCR products digested with Taq I and Xba I, the expected fragment sizes were 1090 and 151 bp for haplotype A, and were 704 and 537 bp for haplotype B. The PCR product-digestion patterns allow us to readily distinguish which haplotype group each product (sample) belongs to. All 12 samples from site 25-a are identified as haplotype A, while 6 and 4 in 25-b are respectively A and B. In 25-c, 7 and 8 belongs to A and B respectively and total 19 samples from 25-d are all B haplotype. Such a distribution pattern strongly supports that secondary contact and introgression between lineages A and B had occurred in this region.

Nested-clade phylogeographic analysis

Haplotypes separated by up to eight mutational steps were connected following the hierarchical nesting procedure. The two haplotype clades were grouped into two independent networks (A and B). One hundred and six haplotypes and 5 nesting levels were included in cladogram A (Fig. 3), 12 haplotypes and 2 nesting levels in cladogram B (Fig. 4). For both cladograms, fifteen nested groups showed significant geographical associations in nested contingency test (Table 4). Further, the nested distance analysis showed that 34 nested groups in both haplotype cladograms have geographic distributions that differ significantly from random expectations (results not shown). Contiguous range expansion, isolation by distance, long-distance colonization and past fragmentation were inferred at different levels (Table 5), which facilitated the identification of additional sampling areas.

thumbnailFigure 3. Nested-clade design on Cyt b haplotype parsimony network of P. nigromaculata clade A. Each line represents a single mutational change, irrespective of length. Black circles represent unsampled or extinct haplotypes. Numbers in italic identify clade-groups on each nesting level.

thumbnailFigure 4. Nested-clade design on Cyt b haplotype parsimony network of P. nigromaculata clade B. Each line represents a single mutational change, irrespective of length. Black circles represent unsampled or extinct haplotypes. Numbers in italic identify clade-groups on each nesting level.

Table 4. Clades with significant associations at the 0.05 level in nested contingency analysis.

For both cladograms, range expansion was inferred at higher levels (Table 5). Tip (young) clades are more widespread geographically than interior (ancestral) haplotypes, which is characteristic of range expansion [18]. For total cladogram A, there was no conclusive inference due to the lack of tip-interior differentiation. Restricted gene flow with isolation by distance is a repeated pattern inferred in clade-groups 2-30, 3-2, 4-1 and 4-2. These four groups correspond to lineage A1 and distribute most sampling localities. For group 3-1 and 3-12 which include haplotypes from coastal population (locality 4), south-western population (locality 16), and most central or eastern populations, we do not have sufficient evidence to discriminate between long-distance movement and the combined effects of gradual movement during a past range expansion and fragmentation. Long-distance colonization (or past fragmentation) is interpreted for group 3-13, which is located in the south-western plateau and corresponds to lineage A2. Inadequate geographical sampling is identified within group 3-11 which corresponds to the regions between northeast and north China (Fig. 1 and 3). Allopatric fragmentation is predicted within group 4-5 and 4-6. The former group is located across southwest China and eastern most regions, while the latter is composed of haplotypes from north-eastern populations and the east China populations.

Table 5. Inferences from the nested clade distance analysis in Chinese range of P. nigromaculata.

Historical demography

Mismatch distributions for the entire group exhibited two distinct modes (Fig. 5). We assume that the left peak represents a stationary expansion corresponding to lineage A, but the "older" peak on the right is unlikely due to demographic expansion. Rather, may result from the deep split between the lineages, perhaps caused by bottlenecking and lineage sorting in allopatry (initial reduction of variation). For the three independent lineages, the pairwise differences all show unimodal patterns, which suggests that demographic expansions occurred within each. Mismatch distributions for lineage B and sub-lineage A1 did not differ significantly (P > 0.05) from expectations under the sudden-expansion model [19]. Whereas, sub-lineage A2 was weakly significant (P = 0.043), but Tajima's D and Fu's Fs are both negative (Table 6). Our estimated time since expansion for lineage A1, A2 and B were different. Expansions of sub-lineage A1 took place contemporarily about 0.08 Mya, whereas expansion was about 0.012 Mya for sub-lineage A2 and 0.028 Mya for lineage B. The sequence variation (Table 6) of each lineage differs in our estimated h and π. The presence of high h and high π (0.975 and 2.8%) in the total populations seems to reflect the fact that different mtDNA lineages were analyzed together. Similarly, the same genetic diversity pattern was observed for sub-lineage A1 (h = 0.972,π = 0.9%), which means this lineage is a large population with admixed local populations. Relatively low h and extremely low π (0.631 and 0.13%) were found in sub-lineage A2, suggesting that a prolonged bottleneck was present in this population. In lineage B, π is low (0.3%) but h is otherwise (0.799), which can be explained by rapid population growth from ancestral populations [20].

thumbnailFigure 5. Mismatch distributions for each major haplotype clade of P. nigromaculata. The abscissa shows the number of pairwise differences between compared haplotypes. The ordinate shows the proportion for each value. Black shows observed frequency distribution while the white bars show the distribution expected under the sudden-expansion model.

Table 6. Demographic parameters of lineages for which range expansion were detected in mismatch distributions.

Based on the results from the MDIV analysis, the posterior distribution for T revealed a sharp peak at 1.16. With substitution rate of 4.4% per Myr [21], our estimate of the timing of lineage divergence, 0.92 Mya (95% credibility internal 1.5 to 0.26 Mya), supports a Pleistocene split. Further within the lineage A, divergence time between the two sub-lineages was estimated to be 0.18 Mya (95% credibility internal 0.10 to 0.32 Mya), roughly corresponding to the Riss glaciation (210,000–135,000 years ago) [22].

ISSR polymorphism and population structure

The 10 ISSR primers generated a total of 127 bands, of which 103 were polymorphic thus showing 81.10% polymorphism and an average genetic diversity of 0.27. The number of bands generated by individual primers varied from 9 to 21, with the average being 12.7. The genetic variation indices at the species level were PPB = 81.10, Ao = 1.80, Ae = 1.45, He = 0.27, and I = 0.40, much higher than the mean values of all populations, especially the PPB (data not shown).

The UPGMA (Unweighted Pair Group Method with Arithmetic mean) dendrogram showed similar topology to that of the mitochondrial phylogeny, showing divergence between the mitochondrial clade A and clade B, and also between sub-clade A1 and sub-clade A2 (data not shown). As expected, we detected potential gene flow between lineage A and lineage B in DL area. Fig. 6 was the PCR banding pattern with primer UBC826. As shown (right side), band L. A can be readily recognized as a lineage A-special band, while band L. B was unique in lineage B. In the 3 individuals (DL1, 2, 6) of the mitochondrial linage B, L. B was detected, while in 5 individuals (DL3, 4, 5, 7 and 8) of lineage A, both bands (L. A and L. B) were amplified (Fig. 6 left part). Interestingly, in individual DL9 which belongs to mitochondrial lineage A (haplotype 14), only band L. B. was detected. This phenomenon was not found in all other local populations, and provided nuclear evidence to support the secondary contact and introgression between lineage A and B at Dongliao area.

thumbnailFigure 6. Amplified patterns generated by using ISSR primer UBC826. Letters and numbers above the figure are codes for individuals from different localities (letters are abbreviations same to Table 1, numbers are corresponding to individuals). "M" represents the molecular marker of 100-bp ladder, with fragment length of each band indicated. L. A and L. B are bands unique to lineage A and lineage B, respectively.

Discussion

Evolutionary history and independent refugia

Amphibians are known to exhibit a higher degree of population subdivision than any other major animal taxon [23]. Our analyses revealed that there are two reciprocally monophyletic lineages within Chinese range of P. nigromaculata: the mainland lineage (A) and the northeast lineage (B). The deep intraspecific split (7.72%) in Cyt b is not of the same magnitude as to some anuran amphibians (maximum of 11.35% in Mantella bernhardi; 13.7% in Hyla arenicolor [24]) and other vertebrate groups (12.5% in Microtus oeconomus [25]). Zhang et al. had also detected a deep haplotype split within P. nigromaculata on mitochondrial 5' control region [17]. The interpretation was that the relatively long-term hibernation of the north-eastern population limited its ability of gene flow with other mainland populations. In other words, genetic drift under natural selection in northeast China resulted in a deep divergence of this frog species. In such a case, selective sweeps would be expected to affect the evolutionary pattern of Cyt b, but this conjecture was not supported by the neutral tests of our dataset which indicated that its patterns of variation were consistent with the hypothesis of neutrality. According to the estimated divergence time of Cyt b haplotypes, the two lineages of P. nigromaculata diverged approximately 0.92 Mya (credibility internal 1.5 to 0.26 Mya), similar to that of other anurans isolation events [26,27]. This divergence time was roughly congruent with the Gunz glaciation (1.2 – 0.9 Mya), which may be the cause of allopatric isolation and the lineage split.

Geographically, lineage A (A1 and A2) covers almost the entire sampling region of this study. Lineage B is restricted to the eastern part of northeast China (Fig. 1). The patterns of diversification suggest complex histories involving both allopatric isolation among refugial areas and prominent patterns of dispersal. We proposed that there were two independent refugia during Pleistocene glaciations. One refugium was the eastern monsoon region and the lower elevations of the south-western plateau. Unlike Europe and North America, China was not covered by continental ice sheets during the Quaternary [28,29], even though the climate was cold and dry. The permafrost had expanded southward by about 10 degrees latitude to the present day Great Wall line and the ice-age "mammoth fauna" had roamed in north China and had even reached the Yangtze River Estuary [30]. Conversely, the climate in eastern China was relatively moist and warm; therefore, eastern China was likely an alternative for this water-dependent frog to live through the cold periods. Zhang also suggested that eastern China was a refugium for temperate and tropical-subtropical faunas during cold stages of the Quaternary [31]. The fact that most haplotypes from eastern local populations are present in interior clades of the parsimony network A again substantiates our east China refugium contention. Although the south-western plateau had been identified with 3–5 glaciations that are dated to late and middle Pleistocene [32], environmental diversity of tropics and subtropics in lower elevations are still maintained [9]. This region has retained many relic floras of Tertiary and was also a refugium for temperate and subtropical floras in cold periods of Quaternary [33]. It is possible for P. nigromaculata to have survived the lower elevations in cold periods and expand afterwards.

We suggest that the other refugium was situated in the present day Korean Peninsula. This Peninsula is situated between 33°–43° N and 124°–132° E, and is adjacent to northeast China. Recent studies on the Quaternary indicated that there were mountain glaciations developed in the Changbai Cordillera [32], which lies between China and North Korea. But the Korean Peninsula, characterized by subtropical mountain climate [10], was not deeply affected by the global climate changes starting in the Pleistocene. The sustained stable environment was presumably favorable to P. nigromaculata. The Korean haplotype clustered with clade B and differed by only 0.009 on Cyt b sequence from B further strengthens the notion of a Korean refugium.

Expansion of lineage A

Using NCA on the Cyt b data, we obtained further insights into the phylogeographic pattern in lineage A of P. nigromaculata. First, contiguous range expansion was predicted for cladogram A in the highest nesting level (5-1, 5-2), which represents an older event in the evolutionary history of P. nigromaculata. This conclusion helps explain the species' current distribution in the north-eastern areas and in the high elevations of south-western plateau where former glaciations occurred. The tip haplotypes are interpreted as being younger than the interior ones [18]. We found that the majority of younger haplotypes distributed in the most peripheral populations of the species' range (north-eastern, south-western and some costal regions). This pattern, together with the unimodal mismatch distributions, supports the range expansion inference of lineage A. Most ancient haplotypes distributed in eastern China also suggest that this area was the origin of the last interglacial range expansion. From the high genetic variation of this lineage and the fact that most common haplotypes are shared among eastern and south-western populations, we can conclude that the eastern and south-western regions were the ice-age refugium for lineage A. From this refugium, P. nigromaculata expanded in all directions to form the current distribution pattern (Fig. 7). Following range expansion, restricted gene flow isolation by distance (IBD) at lower nesting levels (2-30, 3-2, 4-1 and 4-2) took place repeatedly, which seems to represent more recent dynamics of this species. These results likely indicate that gene flow decreased with geographic distance, despite the fact that there were rivers distributed across the frogs' range. It also reflects the reduced vagility of amphibian species. The close relationship between Japanese and clade A haplotypes might be an indication of its origin. Throughout the Quaternary, sea level has risen and fallen as continental ice sheets have waned and waxed [34], and it was about 70–125 m lower than present during the last glacial period when the Eurasia continent was in connection with Japanese archipelago by a land bridge [34-36]. Additionally, the deepest Korean strait makes it controversial of whether a land-bridge connection existed in that place during Late Quaternary [36]. Via the land bridge, we suggest, P. nigromaculata colonized Japan from Chinese eastern coast (Fig. 7). However, further samples from the Japanese population are required in future studies.

thumbnailFigure 7. Map of major dispersal patterns of Chinese population of P. nigromaculata identified in this study. Points show the sampling locations in Fig. 1. Solid arrows indicate dispersal pattern of lineage A. Dotted arrows indicate patterns of dispersal of lineage B. Arrows with star "*" indicate that the pattern is not confirmed. Gray arrows from Chinese eastern coast to Japanese archipelago indicate the dispersal pattern via past land bridge. The symbol "?" indicates that the origin of the population has not been resolved. Gray areas represent past glacial refugia.

For group 2-19 (north-eastern populations QQ and DL were included) in lineage A1, we could not distinguish past fragmentation\long-distance colonization from contiguous range expansion due to the inadequate sampling design. Again, in group 3-11, the two north-eastern populations and two populations near the Shandong Peninsula (locality 3, 8) were included, for which the patterns of fragmentation and isolation by distance were not identified due to sampling gaps. Nonetheless, it may imply that, when recolonizing northeast China during the last interglacial period, these peninsula-adjacent populations were founders. The recolonizing route might have been via the Shandong region to northeast China (Fig. 7). This is also consistent with the topology of MP tree (Fig. 2), in which the north-eastern haplotypes are clustered with the peninsula-adjacent haplotypes. The route through which lineage A1 recolonized Baoshan (Location 16) is somewhat puzzlings. The Baoshan population is an isolated population hundreds of kilometres apart from the nearest P. nigromaculata population in Sichuan. It is recognized not only by the specimens we obtained in 2002, but also by 63 specimens collected in 1998 by the researchers from Chengdu Institute of Biology, Chinese Academy of Sciences (Fei and Jiang, personal communication). However, due to the relatively high diversity level (both in mitochondrial haplotypes and ISSR profiles), and the relatively few individuals studied, we were unable to unequivocally identify the origin of the Baoshan population.

Group 3-13, corresponding to lineage A2 is composed of all south-western plateau populations characterized by long-distance colonization, and possibly couple with subsequent fragmentation. The unimodal mismatch distribution in this lineage suggests recent rapid colonization. The estimated time since expansion was about 0.012 Mya, after the last glacial maximum (LGM, 0.018 Mya). According to the palaeoclimatic records of this region [9], ice sheet on the mountain had shifted vertically in response to the oscillations of ice ages, natural zones have also moved vertically since the Pleistocene, but retained environmental diversity of tropics and subtropics in lower elevations. P. nigromaculata likely survived in these areas in ice-ages and rapidly expanded subsequently. Based on the dispersal theory of Hewitt [2], once a place is occupied by a leading edge expansion population through dispersal, it will be much more difficult for those behind them to follow. The imprint of this is reduced genomic variability. The strikingly low levels of genetic variation, especially the observation that multiple individuals from distant geographic regions share one single haplotype (No. 20 in Table 2: 5 from DS, 3 from WG, 10 from GH, 5 from RH and 4 from DZ) strongly supports rapid expansion.

Recolonization to northeast China of lineage B

Contiguous range expansion is also inferred for total cladogram B, which restricted to northeast China. Mismatch distribution also lend strong support to the conclusion of population expansion of this lineage. The last global glaciation corresponds respectively to the Würm, Wisconsin and Weichsel glaciation in Alps, North America and western Europe [37]. In China the same period is called Dali glaciation, which began at about 0.05 Mya. In Europe and North America, the most extensive glacial period was during the last glacial maximum (LGM) beginning 0.020 ~ 0.018 Mya, whereas in China the most extensive glacial period was not LGM, but the early stage (0.054 ~ 0.044 Mya) of the Dali glaciation [38] when northeastern China was cold and the Changbai Cordillera was covered with ice caps [32]. The Korean Peninsula was assumed to be a refugium not only because the relatively mild climate of the Peninsula, but also because of the phylogenetic basal position of the Korean haplotype. Our estimated time of expansion was about 0.028 Mya for the P. nigromaculata lineage B and it was likely after the early extensive glacial period in Dali glaciation, that this lineage colonized northeast China via the Changbai Cordillera, but not after the last glacial maximum when many lineages in Europe and North America underwent expansions. The current distribution of lineage B suggests that after recolonizing northeast China, this lineage remained in situ rather than expanding southward, at least in the present interglacial period.

Secondary contact and introgressive gene flow

Secondary contact of previously isolated populations after the last glacial maximum in Europe and North-America are well established [2,39-42]. In the present study, we uncovered a population where haplotype lineage A and B of P. nigromaculata coexisted in Dongliao area (locality 25). The DL population, located in northeast China, is the only one composed of both lineages, a conclusion confirmed by nucleotide sequences, PCR-RFLP, and ISSR patterns. Hence, the DL region (especially 25-b and 25-c with a distance only 2.2 km) was likely a rather narrow secondary contact zone during the current interglacial period. The full range of the secondary contact zone remains to be defined.

In general, hybrid or contact zones between populations or species of varying degrees of differentiation provide a potentially important source of information about how taxa diverge and populations interact. We could not conclude that in the area of sympatry, genetic isolation exists between the two lineages based solely on high mtDNA sequence divergence of 7.7%. Indeed our nuclear ISSR data suggest some inter-lineage gene flow in the DL population. Two lines of evidence support this contention. First, in 5 individuals of mitochondrial lineage A, both bands of L. A. and L. B. were amplified. Second, one lineage A individual only possessed L. B band. The former may arise from the hybridization between males of lineage B and females of lineage A, and the later might be caused by backcross of female hybrids of the former with males of lineage B. In other words, these two divergent lineages may have experienced a recent hybridization resulting in introgressive gene flow.

Conclusion

Our phylogeographic study suggests that, two mitochondrial lineages diverged about 0.9 Mya beginning in independent refugia during the Pleistocene glaciation. Two main clades, A and B, differing by c. 7.72% sequence divergence were detected. The A clade is further subdivided into two sub-clade, A1 and A2 differing by 1.22%. The two refugia seemed to be located in east China and the lower elevations of south-western plateau, and in Korean Peninsula, respectively. Contiguous range expansion from the two refugia during last interglacial period resulted in the current distribution pattern of this frog species. Lineage A covers almost the entire sampling region of this study. Lineage B is restricted to the eastern part of northeast China. Secondary contact was detected in DL area where both lineages coexisted, as confirmed by nucleotide sequences, PCR-RFLP and ISSR analyses. ISSR results provided further evidence for introgressive gene flow between the two lineages in this area. Although climate events of the Pleistocene in East Asia did not seem to resemble that of Europe or North America, they do have had marked effect on the historical distribution and intraspecific divergence of amphibians as in Europe and North America.

Methods

Sample collection and Molecular techniques

Sampling was designed to cover most of the range of P. nigromaculata across the temperate and subtropical zones in China. We obtained 262 individuals from 28 localities during the year 2002 and 2003. Two hundred and six samples were used in cytochrome b (Cyt b) sequence analyses, the other 56, collected from 4 sites (25-a, 25-b, 25-c and 25-d) in Dongliao area (locality 25, DL, Fig. 1) were used in PCR-RFLP analyses to identify secondary contact between haplotype group A and B. The distance between sites 25-a and 25-b is 11.9 km, and that of 25-b and 25-c, 25-c and 25-d is 2.2 and 28.1 km, respectively. For inter simple sequence repeat (ISSR) analyses, six or nine (just in DL) individuals were sampled from the 28 localities mentioned above. Rana catesbeiana [GenBank: AF205089] and Rana rugosa [GenBank: AF205093] were used as outgroup taxa. Four conspecific haplotypes from Korea (KR) and Japan (JP-1, 2, 3) [GenBank: AF205087, AY315755~AY315757] were also included in the phylogenetic analyses to help infer relationships with Chinese haplotypes. The sampling code and population localities are given in Table 1, and the map of collecting localities is shown in Fig. 1.

Tissue was taken from leg muscle and stored in 95% ethanol. Extraction of DNA followed standard phenol/chloroform extraction methods. The complete 1143 bp of mitochondrial Cyt b gene was amplified with primers Glu-f: GAC TCT AAC CTG GAC CAA TAG-3' and contr5'r: TAA ATT TAT GCT CTA TAC A-3', which were designed based on the published mtDNA complete sequence of Rana nigromaculata [GenBank: AB043889] [43]. The polymerase chain reaction (PCR) was carried out in 30 μL volumes containing 3.0 μL 10 × Buffer, 2.0 mM MgCl2, 0.2 mM each dNTP, 0.25 μM each of primers, 5–10 ng of template DNA and 1.0 U of Taq DNA polymerase (Promega). The PCR cycling parameters were 95°C for 4 min, 35 cycles of 95°C for 40 s, 50°C for 45 s, 72°C for 1 min, and 1 cycle of 72°C for 7 min. Templates for sequencing were purified and sequenced in both directions by United Gene Holdings (Shanghai). Nucleotide sequences were initially aligned using CLUSTAL X Version 1.81[44] and corrected manually. Complete sequences were assembled using Seqman II (DNASTAR). The program DNAClub which could predict cutting sites in a DNA sequence was used to search for restriction endonucleases that would yield unique, specific restriction digestion profiles for each haplotype group. For the 56 unsequenced samples, we used two restriction enzymes of Taq I and Xba I to digest the Cyt b products simultaneously, by which we could distinguish the two sharply diverged haplotype groups rapidly and unambiguously.

We also investigated diversification and potential gene flow using inter simple sequence repeat PCR (ISSR). Ten 3'-anchored primers (UBC807, 808, 809, 811, 817, 818, 825, 826, 827 and 835) from University of British Columbia (ISSR set #9) were used for the study. The reaction mixture (25 μL) contained 2.5 μL of 10 × reaction buffer, 25 ng of DNA template, 0.5 μM of a single primer, 2.0 mM MgCl2, 0.2 mM of each dNTP and 1.0 U of Taq DNA Polymerase (Promega). The PCR products were separated on 2% agarose gels buffered with 0.5 × TBE, detected by staining with ethidum bromide, and photographed under ultraviolet light. Molecular weights were estimated using a 100-bp DNA ladder.

Phylogenetic analysis

We used PAUP*4.0b10 [45] for phylogenetic analyses using maximum parsimony (MP), maximum likelihood (ML) and neighbor-joining (NJ) methods. Most local populations showed low level of genetic differentiation and were composed of haplotypes with only 1 to 3 nucleotide changes. For ML phylogenetic reconstruction, these similar haplotypes would considerably increase the computation time. Therefore, the dataset comprised of all haplotypes was only used for MP reconstruction. To save the computing time, we cut the sequences with less than four-step substitutions for ML analysis to a 49 ingroup exemplar dataset including 45 from China, 1 from Korea and 3 from Japan. The best-fit likelihood model of sequence evolution was HKY + I + G, which was chosen on the basis of hierarchical likelihood-ratio tests (hLRTs) as implemented in MODELTEST 3.06 [46]. The parameters of the model were calculated using PAUP*. The base frequencies were 0.2582, 0.3340, 0.1334 and 0.2743 for A, C, G and T, respectively, transition/transversion ratio ti/tv = 4.2823, proportion of invariable sites Pinvar = 0.5841, gamma distribution shape parameter α = 1.9961. We also conducted MP and NJ analyses for our reduced exemplar taxa. ML heuristic search was conducted using the following options: addition sequence = 'as-is', with 1 tree held at each step, TBR swapping algorithm, Collapse and Multrees options in effect, steepest descent option not in effect. Nodal support of the ML tree was estimated using bootstrap method [47] with 100 pseudoreplicates; Maximum parsimony (MP) analyses were conducted using heuristic search with 100 random sequence additions and tree-bisection-reconnection (TBR) branch swapping. Robustness of the MP trees was assessed by 1000 bootstrap replicates; NJ reconstruction was conducted using distances calculated under the HKY model, support for the topology was evaluated by bootstrapping 5000 replicates. We also used the chosen model to calculate the diversity indices and the sequence distances.

We calculated the statistical tests D* and F* proposed by Fu and Li [48], and the D test statistic proposed by Tajima [49] implemented in DnaSP 3.0 [50] for testing the hypothesis that all mutations are selectively neutral [51].

Nested-clade phylogeographic analysis

We applied nested clade analysis (NCA [18,52]) to gain further insight into the demographic history of P. nigromaculata. The statistical parsimony [52] cladogram networks were constructed using TCS version 1.13 [53]. The connection for haplotypes at 95% confidence level were nested into higher level clades followed the algorithm given by Templeton et al. [52]. We used GEODIS version 2.0 [54] to test the null hypothesis of no association between haplotypes and their geographical locations by randomly permuting the observations within a nesting clade across geographical locations. After 1000 random permutations, the clade and nested clade distances (Dc and Dn), the interior-tip distances (I-Tc and I-Tn) were calculated. The observed Dcand Dn were then contrasted to the null distribution to infer the statistically significantly large (L) and small (S) distances at 5% deviation level. For the clades with significant geographical associations, the recent published inference key [55] was used to infer possible causation.

Demographic history

Mismatch distributions were created using ARLEQUIN version 2.000 [56] to detect historical demographic expansions for the major lineages (with or without range expansion in NCA analysis). We also calculated Tajima's D [49,57] and performed Fu's Fs test [58] to confirm evidence of demographic expansion. We estimated the time since expansion with the equation τ = 2ut [19], where τ is the mode of the mismatch distribution, an index of time since expansion expressed in units of mutational time, u is the mutation rate for the whole sequence, and t is the time since expansion. Finally, the comparisons of haplotype diversity (h) and nucleotide diversity (π) within different lineages were used to infer historical demography and examine the biological inference from our nested clade analysis.

Divergence time between lineages

We analyzed the divergence time between the two lineages with a coalescent-based approach using the program MDIV [59]. We ran the program under the 'finite model', which does not assume that only one mutation occurs at each site but takes into account the possibility of multiple hits, differences in the nucleotide frequencies and the presence of transition/transversion bias. MDIV estimated values for theta (θ = 4Neμ), migration rate (M = 2Nem), time of divergence (T = t/2Ne) [59]. We then ran each simulation 5 × 106 times with a 10% burn- in period as suggested. Likelihood values for θ, M and T were calculated and the value with the highest posterior probability accepted as the best estimate. Values for T were calculated using an evolution rate estimated at 4.4% per Myr [21].

ISSR analysis

Unequivocally scorable and consistently reproducible amplified ISSR bands were scored as present (1) and absent (0), each of which was treated as an independent character regardless of its intensity. Fragments of the same molecular weight were considered as the same locus.

The binary matrix was used under the Hardy-Weinberg equilibrium to calculate the percentage of polymorphism (PPB), Nei's gene diversity (He), the observed number of alleles per locus (Ao), the effective number of alleles per locus (Ae), Shannon's information index (I), total gene diversity, the level of gene flow using POPGENE version 1.31 [60].

Furthermore, to illustrate the relationship among populations, genetic divergence (Fst) were calculated using software ARLEQUIN version 2.000. We generated UPGMA (unweighted pair group method with arithmetical averages) trees using the SAHN- clustering and TREE programs from the NTSYS-pc, version 2.1 package [61]. In this tree, individuals were enforced to their populations except for those from the DL population (site 25), which were subdivided into DLA and DLB according to their ISSR patterns.

Abbreviations

ISSR: Inter simple sequence repeat; PCR-RFLP: Polymerase chain reaction-restriction fragment length polymorphism; MP: Maximum parsimony; ML: Maximum likelihood; NJ: Neighbor-joining; UPGMA: Unweighted Pair Group Method with Arithmetic mean; NCA: Nested clade analysis; IBD: Isolation by distance; LGM: Last glacial maximum.

Authors' contributions

HZ isolated genomic DNA, and conducted the amplification and following RFLP analysis of the complete Cyt b gene. She collected, analyzed and summarized the data, and drafted the manuscript; JY participated in analysis and interpretation of data, and helped draft the manuscript; GQZ performed ISSR analysis; KYZ conceived the study and participated in its design and data interpretation, and preparing the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We would like to thank all colleagues who kindly helped us during samples collecting. We also thank Xiongfei Zhang, Keran Bi, Xiaojun Yang and Xiuling Lü for kind assistance in the laboratory and helpful suggestions on data analyses; Liang Fei and Jianping Jiang for information about specimens of P. nigromaculata collected from Baoshan, and for identification of some P. nigromaculata individuals. We are grateful to two anonymous reviewers for valuable comments and language corrections that greatly improved the manuscript. Support for this research was provided by the National Natural Science Foundation of China (No. 30470249) to KYZ.

References

  1. Hewitt GM: Some genetic consequences of ice ages, and their role in divergence and speciation.

    Biol J Linn Soc 1996, 58:247-276. OpenURL

  2. Hewitt G: The genetic legacy of the Quaternary ice ages.

    Nature 2000, 405(6789):907-913. PubMed Abstract | Publisher Full Text OpenURL

  3. Avise JC, Walker D, Johns GC: Speciation durations and Pleistocene effects on vertebrate phylogeography.

    Proc Biol Sci 1998, 265(1407):1707-1712. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Avise JC, Walker D: Pleistocene phylogeographic effects on avian populations and the speciation process.

    Proc Bio Sci 1998/05/07 edition. 1998, 265(1395):457-463. Publisher Full Text OpenURL

  5. Zink RM, Klicka J: The tempo of avian diversification: a comment on Johnson and Cicero.

    Evolution 2006/04/14 edition. 2006, 60(2):411-412. PubMed Abstract OpenURL

  6. Cicero C, Johnson NK: The tempo of avian diversification: Reply.

    Evolution 2006, 60(2):413-414. OpenURL

  7. Lovette IJ: Glacial cycles and the tempo of avian speciation.

    Trends Ecol Evol 2006/05/17 edition. 2005, 20(2):57-59. PubMed Abstract | Publisher Full Text OpenURL

  8. Zhang RZ: Zoogeography of China. Beijing , Science Press; 1999.

  9. Zhang RZ: Geological events and mammalian distribution in China.

    Acta Zoologica Sinica 2002, 48(2):141-153. OpenURL

  10. Liu DS, Li ZG: Geography of Asia. Beijing , Commercial Press; 1996.

  11. Fei L, Ye CY, Huang YZ: Key to Chinese Amphibia. Chongqing , Chongqing Branch, Science and Technology Literature Publishing House; 1990.

  12. Jiang JP, Zhou KY: Evolutionary relationships among Chinese ranid frogs inferred from mitochondrial DNA sequences of 12S rRNA gene.

    Acta Zoologica Sinica 2001, 47(1):38-44. OpenURL

  13. Jiang JP, Zhou KY: Phylogenetic relationships among Chinese ranids: inferred from sequence data set of 12S and 16S rDNA.

    Herpe J 2005, 15(1):1-8. OpenURL

  14. Frost DR, Grant T, Faivovich J, Bain RH, Haas A, Haddad CFB, De Sá RO, Channing A, Wilkinson M, Donnellan SC, Raxworthy CJ, Campbell JA, Blotto BL, Moler P, Drewes RC, Nussbaum RA, Lynch JD, Green DM, Wheeler WC: The amphibian tree of life.

    Bulletin of the AMNH 2006, 297:1-370. OpenURL

  15. Che J, Pang J, Zhao H, Wu GF, Zhao EM, Zhang YP: Phylogeny of Raninae (Anura: Ranidae) inferred from mitochondrial and nuclear sequences.

    Mol Phylogenet Evol 2007, 43(1):1-13. PubMed Abstract | Publisher Full Text OpenURL

  16. Yang YH, Zhang DX, Li YM, J. JY: Mitochondrial DNA diversity and preliminary biogeographic inference of the evolutionary history of the black-spotted pond frog Rana nigromaculata populations in China.

    Acta Zoologica Sinica 2004, 50(2):193-201. OpenURL

  17. Zhang XF, Zhou KY, Chang Q: Population genetic structure of Pelophylax nigromaculata in Chinese mainland, based on mtDNA control region sequences.

    Acta Genetica Sinica 2004, 31(11):1232-1240. PubMed Abstract OpenURL

  18. Templeton AR: Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history.

    Mol Ecol 1998, 7(4):381-397. PubMed Abstract | Publisher Full Text OpenURL

  19. Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences.

    Mol Biol Evol 1992, 9(3):552-569. PubMed Abstract | Publisher Full Text OpenURL

  20. Avise JC: Phylogeography: The History and Formation of Species. Cambridge , Harvard University Press; 2000.

  21. Sumida M, Kaneda H, Kato Y, Kanamori Y, Yonekawa H, Nishioka M: Sequence variation and structural conservation in the D-loop region and flanking genes of mitochondrial DNA from Japanese pond frogs.

    Genes Genet Syst 2000, 75(2):79-92. PubMed Abstract | Publisher Full Text OpenURL

  22. Gillespie A, Molnar P: Asynchronous maximum advances of mountain and continental glaciers.

    Rev Geophys 1995, 33(3):311-364. Publisher Full Text OpenURL

  23. Palo JU, Schmeller DS, Laurila A, Primmer CR, Kuzmin SL, Merila J: High degree of population subdivision in a widespread amphibian.

    Mol Ecol 2004, 13(9):2631-2644. PubMed Abstract | Publisher Full Text OpenURL

  24. Barber PH: Phylogeography of the canyon treefrog, Hyla arenicolor (Cope) based on mitochondrial DNA sequence data.

    Mol Ecol 1999, 8(4):547-562. PubMed Abstract | Publisher Full Text OpenURL

  25. Galbreath KE, Cook JA: Genetic consequences of Pleistocene glaciations for the tundra vole (Microtus oeconomus) in Beringia.

    Mol Ecol 2004, 13(1):135-148. PubMed Abstract | Publisher Full Text OpenURL

  26. Austin JD, Lougheed SC, Boag PT: Discordant temporal and geographic patterns in maternal lineages of eastern north American frogs, Rana catesbeiana (Ranidae) and Pseudacris crucifer (Hylidae).

    Mol Phylogenet Evol 2004, 32(3):799-816. PubMed Abstract | Publisher Full Text OpenURL

  27. Nielson M, Lohman K, Sullivan J: Phylogeography of the tailed frog (Ascaphus truei): Implications for the biogeography of the Pacific Northwest.

    Evolution 2001, 55(1):147-160. PubMed Abstract OpenURL

  28. Zhou YL: Geography of the Vegetation in Northeast China. Beijing , Science Press; 1997.

  29. Kahlke HD: On the complex of the Stegodon-Ailuropoda-Fauna of Southern China and the chronological position of Gigantopithecus blacki v. Koenigswald.

    Vertebrate Palastatica 1961, 5(2):83-108. OpenURL

  30. Jin CZ, Xu QQ, Zheng JJ: On the dispersal events of Mammuthus during the late Late Pleistocene.

    Vertebrate Palastatica 1998, 36(1):47-53. OpenURL

  31. Zhang RZ: Relict distribution of land vertebrates and Quaternary glaciation in China.

    Acta Zoologica Sinica 2004, 50(5):841-851. OpenURL

  32. Li JJ, Shu Q, Zhou SZ, Zhao ZJ, Zhang JM: Review and prospects of Quaternary glaciation research in China.

    J Glaciol Geocryol 2004, 26(3):235-243. OpenURL

  33. Zeng J, Liu YC: The research of the geographic distribution of the rare and endangered plants in Sichuan and their characteristics of flora.

    Journal of Chongqing Teachers College 1995, 12(4):39-47. OpenURL

  34. Edwards RL, Cutler KB, Cheng H, Gallup CD: Geochemical Evidence for Quaternary Sea-level Changes.

    Treatise on Geochemistry 2003, 6:343-364. OpenURL

  35. Bard E, Hamelin B, Arnold M, Montaggioni L, Cabioch G, Faure G, Rougerie F: Deglacial sea-level record from Tahiti corals and the timing of global meltwater discharge.

    Nature 1996, 382(6588):241-244. Publisher Full Text OpenURL

  36. Millien-Parra V, Jaeger JJ: Island biogeography of the Japanese terrestrial mammal assemblages: an example of a relict fauna.

    J Biogeography 1999, 26(5):959-972. Publisher Full Text OpenURL

  37. Liu JQ, Ni YY, Chu GQ: Main Palaeoclimatic events in the Quaternary.

    Quatern Sci 2001, 21(3):239-248. OpenURL

  38. Cui ZJ, Zhang W: Discussion about the glacier extent and advance/retreat asynchrony during the last glaciation.

    J Glaciol Geocryol 2003, 25(5):510-516. OpenURL

  39. Remington CL: Suture-zones of hybrid interaction between recently joined biotas.

    Evol Biol 1968, 2:321-428. OpenURL

  40. Burbrink FT, Lawson R, Slowinski JB: Mitochondrial DNA phylogeography of the polytypic North American rat snake (Elaphe obsoleta): A critique of the subspecies concept.

    Evolution 2000, 54(6):2107-2118. PubMed Abstract OpenURL

  41. Zamudio KR, Savage WK: Historical isolation, range expansion, and secondary contact of two highly divergent mitochondrial lineages in spotted salamanders (Ambystoma maculatum).

    Evolution 2003, 57(7):1631-1652. PubMed Abstract OpenURL

  42. Leaché AD, Cole CJ: Hybridization between multiple fence lizard lineages in an ecotone: locally discordant variation in mitochondrial DNA, chromosomes, and morphology.

    Mol Ecol 2007, 16(5):1035-1054. PubMed Abstract | Publisher Full Text OpenURL

  43. Sumida M, Kanamori Y, Kaneda H, Kato Y, Nishioka M, Hasegawa M, Yonekawa H: Complete nucleotide sequence and gene rearrangement of the mitochondrial genome of the Japanese pond frog Rana nigromaculata.

    Genes Genet Syst 2001, 76(5):311-325. PubMed Abstract | Publisher Full Text OpenURL

  44. Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG: The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools.

    Nucleic Acids Res 1997, 25(24):4876-4882. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  45. Swofford DL: PAUP* — Phylogenetic Analysis Using Parsimony (* and Other Methods. Beta version 4.0b10). Sunderland , Sinauer Associates.;

  46. Posada D, Crandall KA: MODELTEST: testing the model of DNA substitution.

    Bioinformatics 1998, 14(9):817-818. PubMed Abstract | Publisher Full Text OpenURL

  47. Felsenstein J: Confidence limits on phylogenetics: an approach using the bootstrap.

    Evolution 1985, 39:783-791. Publisher Full Text OpenURL

  48. Fu YX, Li WH: Statistical tests of neutrality of mutations.

    Genetics 1993, 133(3):693-709. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  49. Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism.

    Genetics 1989, 123(3):585-595. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  50. Rozas J, Rozas R: DnaSP version 3: an integrated program for molecular population genetics and molecular evolution analysis.

    Bioinformatics 1999, 15(2):174-175. PubMed Abstract | Publisher Full Text OpenURL

  51. Kimura M: The neutral theory of molecular evolution. Cambridge, Massachusetts , Cambridge University Press; 1983.

  52. Templeton AR, Crandall KA, Sing CF: A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation.

    Genetics 1992, 132(2):619-633. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  53. Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies.

    Mol Ecol 2000, 9(10):1657-1659. PubMed Abstract | Publisher Full Text OpenURL

  54. Posada D, Crandall KA, Templeton AR: GeoDis: a program for the cladistic nested analysis of the geographical distribution of genetic haplotypes.

    Mol Ecol 2000, 9(4):487-488. PubMed Abstract | Publisher Full Text OpenURL

  55. Templeton AR: Statistical phylogeography: methods of evaluating and minimizing inference errors.

    Mol Ecol 2004, 13(4):789-809. PubMed Abstract | Publisher Full Text OpenURL

  56. Schnerider S, Roessli D, Excoffier L: ARLEQUIN ver. 2.000 A Software for Population Genetic Data Analysis. Geneva , Genetics and Biometry Laboratory, University of Geneva; 2000.

  57. Tajima F: The effect of change in population size on DNA polymorphism.

    Genetics 1989, 123(3):597-601. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  58. Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection.

    Genetics 1997, 147(2):915-925. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  59. Nielsen R, Wakeley J: Distinguishing migration from isolation: A Markov chain Monte Carlo approach.

    Genetics 2001, 158(2):885-896. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  60. Yeh FC, Yang RC, Boyle T: POPGENE, the user friendly shareware for population genetic analysis. Edmonton, Canada , Molecular Biology and Biotechnology Center, University of Alberta; 1999.

  61. Rohlf FJ: NTSYS-pc, Numerical Taxonomy and Multivariate Analysis System, version 2.1. Setauket, New York , Exeter Software; 2000.