Phylogeography and Demographic History of Chinese Black-Spotted Frog Populations (Pelophylax nigromaculata): Evidence for Independent Refugia Expansion and Secondary Contact
© Zhang et al; licensee BioMed Central Ltd. 2008
Received: 06 August 2007
Accepted: 24 January 2008
Published: 24 January 2008
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.
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.
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.
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 . 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 . 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. . 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 . 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 .
Collection data for P. nigromaculata specimens included in this study.
County, Province (Abbr.)
GenBank Accession No.
Guangde, Anhui (GD)
AY803813–AY803817, AY803831, DQ006243–DQ006245
Yuexi, Anhui (YX)
Ganyu, Jiangsu (GU)
Cixi, Zhejiang (CX)
Xiangshan, Zhejiang (XS)
Laiyang, Shandong (LA)
Linyi, Shandong (LI)
Zaozhuang, Shandong (ZZ)
Luquan, Hebei (LQ)
Huangchuan, Henan (HC)
3, 41–48, 77
AY803823, AY803848–AY803852, DQ006246–DQ006249
Shashi, Hubei (SS)
AY803823, AY803851, AY803853–AY803855, DQ006258–DQ006261
Nanxian, Hunan (NX)
53, 62–68, 77
AY803844, AY803856–AY803858, DQ006253–DQ006256
Gao'an, Jiangxi (GA)
Yiyang, Jiangxi (YY)
Xiajiang, Jiangxi (XJ)
Baoshan, Yunnan (BS)
Guilin, Guangxi (GL)
Dushan, Guizhou (DS)
AY803873–AY803874, AY803877, DQ006237
Guiyang, Guizhou (GY)
Renhuai, Guizhou (RH)
Dazhu, Sichuan (DZ)
Guanghan, Sichuan (GH)
Wugang, Hunan (WG)
Qiqiha'er, Heilongjiang (QQ)
14, 19, 69
Dongliao, Jilin (DL)
Dongliao, Jilin (DL)
Dongliao, Jilin (DL)
Dongliao, Jilin (DL)
Dongliao, Jilin (DL)
Liuhe, Jilin (LH)
Tonghua, Jilin (TH)
Liaoyang, Liaoning (LY)
AY803889, AY803894–AY803895, DQ006252
Cytochrome b variation
Distribution of Cyt b haplotypes shared in different local populations of P. nigromaculata.
Cyt b phylogenetics
Cyt b sequence distances and the standard deviations among major lineages under the HKY model.
1.22% ± 0.002
4.06% ± 0.005
4.09% ± 0.005
3.84% ± 0.005
7.72% ± 0.008
7.07% ± 0.008
7.73% ± 0.008
7.08% ± 0.008
7.69% ± 0.008
7.69% ± 0.008
7.71% ± 0.008
6.86% ± 0.008
1.10% ± 0.003
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.
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
Clades with significant associations at the 0.05 level in nested contingency analysis.
Permutational χ2 statistic
Total cladogram A
Total cladogram B
Inferences from the nested clade distance analysis in Chinese range of P. nigromaculata.
Contiguous range expansion
Contiguous range expansion
Restricted gene flow with isolation by distance
Fragmentation or isolation by distance
Past fragmentation or long-distance colonization
Range expansion, long-distance colonization or past fragmentation
Fragmentation or isolation by distance
Long-distance colonization or past fragmentation
Range expansion, long-distance colonization or past fragmentation
Restricted gene flow with isolation by distance
Total cladogram B (3)
Contiguous range expansion
Restricted gene flow with isolation by distance
Demographic parameters of lineages for which range expansion were detected in mismatch distributions.
Tajima' D (P-value)
Fu's Fs (P-value)
0.975 ± 0.006
0.028 ± 0.014
0.972 ± 0.008
0.009 ± 0.005
0.631 ± 0.080
0.001 ± 0.001
0.799 ± 0.043
0.003 ± 0.002
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 , 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) .
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, A o = 1.80, A e = 1.45, H e = 0.27, and I = 0.40, much higher than the mean values of all populations, especially the PPB (data not shown).
Evolutionary history and independent refugia
Amphibians are known to exhibit a higher degree of population subdivision than any other major animal taxon . 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) and other vertebrate groups (12.5% in Microtus oeconomus). Zhang et al. had also detected a deep haplotype split within P. nigromaculata on mitochondrial 5' control region . 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 . 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 . 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 , environmental diversity of tropics and subtropics in lower elevations are still maintained . This region has retained many relic floras of Tertiary and was also a refugium for temperate and subtropical floras in cold periods of Quaternary . 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 , which lies between China and North Korea. But the Korean Peninsula, characterized by subtropical mountain climate , 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
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 , 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 , 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 . 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  when northeastern China was cold and the Changbai Cordillera was covered with ice caps . 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.
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.
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] . 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 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.
We used PAUP*4.0b10  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 . 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  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 , and the D test statistic proposed by Tajima  implemented in DnaSP 3.0  for testing the hypothesis that all mutations are selectively neutral .
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  cladogram networks were constructed using TCS version 1.13 . The connection for haplotypes at 95% confidence level were nested into higher level clades followed the algorithm given by Templeton et al. . We used GEODIS version 2.0  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 (D c and D n ), the interior-tip distances (I-T c and I-T n ) were calculated. The observed D c and D n 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  was used to infer possible causation.
Mismatch distributions were created using ARLEQUIN version 2.000  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  to confirm evidence of demographic expansion. We estimated the time since expansion with the equation τ = 2ut, 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 . 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 (θ = 4N e μ), migration rate (M = 2N e m), time of divergence (T = t/2N e ) . 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 .
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 (H e ), the observed number of alleles per locus (A o ), the effective number of alleles per locus (A e ), Shannon's information index (I), total gene diversity, the level of gene flow using POPGENE version 1.31 .
Furthermore, to illustrate the relationship among populations, genetic divergence (F st ) 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 . 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.
Inter simple sequence repeat
Polymerase chain reaction-restriction fragment length polymorphism
Unweighted Pair Group Method with Arithmetic mean
Nested clade analysis
Isolation by distance
Last glacial maximum.
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.
- Hewitt GM: Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc. 1996, 58: 247-276.View ArticleGoogle Scholar
- Hewitt G: The genetic legacy of the Quaternary ice ages. Nature. 2000, 405 (6789): 907-913. 10.1038/35016000.View ArticlePubMedGoogle Scholar
- Avise JC, Walker D, Johns GC: Speciation durations and Pleistocene effects on vertebrate phylogeography. Proc Biol Sci. 1998, 265 (1407): 1707-1712. 10.1098/rspb.1998.0492.PubMed CentralView ArticlePubMedGoogle Scholar
- Avise JC, Walker D: Pleistocene phylogeographic effects on avian populations and the speciation process. Proc Bio Sci. 1998, 265 (1395): 457-463. 10.1098/rspb.1998.0317. 1998/05/07View ArticleGoogle Scholar
- Zink RM, Klicka J: The tempo of avian diversification: a comment on Johnson and Cicero. Evolution. 2006, 60 (2): 411-412. 2006/04/14View ArticlePubMedGoogle Scholar
- Cicero C, Johnson NK: The tempo of avian diversification: Reply. Evolution. 2006, 60 (2): 413-414.View ArticleGoogle Scholar
- Lovette IJ: Glacial cycles and the tempo of avian speciation. Trends Ecol Evol. 2005, 20 (2): 57-59. 10.1016/j.tree.2004.11.011. 2006/05/17View ArticlePubMedGoogle Scholar
- Zhang RZ: Zoogeography of China. 1999, Beijing , Science PressGoogle Scholar
- Zhang RZ: Geological events and mammalian distribution in China. Acta Zoologica Sinica. 2002, 48 (2): 141-153.Google Scholar
- Liu DS, Li ZG: Geography of Asia. 1996, Beijing , Commercial PressGoogle Scholar
- Fei L, Ye CY, Huang YZ: Key to Chinese Amphibia. 1990, Chongqing , Chongqing Branch, Science and Technology Literature Publishing HouseGoogle Scholar
- 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.Google Scholar
- 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.Google Scholar
- 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.Google Scholar
- 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. 10.1016/j.ympev.2006.11.032.View ArticlePubMedGoogle Scholar
- 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.Google Scholar
- 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.PubMedGoogle Scholar
- Templeton AR: Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history. Mol Ecol. 1998, 7 (4): 381-397. 10.1046/j.1365-294x.1998.00308.x.View ArticlePubMedGoogle Scholar
- Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992, 9 (3): 552-569.PubMedGoogle Scholar
- Avise JC: Phylogeography: The History and Formation of Species. 2000, Cambridge , Harvard University PressGoogle Scholar
- 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. 10.1266/ggs.75.79.View ArticlePubMedGoogle Scholar
- Gillespie A, Molnar P: Asynchronous maximum advances of mountain and continental glaciers. Rev Geophys. 1995, 33 (3): 311-364. 10.1029/95RG00995.View ArticleGoogle Scholar
- 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. 10.1111/j.1365-294X.2004.02269.x.View ArticlePubMedGoogle Scholar
- Barber PH: Phylogeography of the canyon treefrog, Hyla arenicolor (Cope) based on mitochondrial DNA sequence data. Mol Ecol. 1999, 8 (4): 547-562. 10.1046/j.1365-294x.1999.00593.x.View ArticlePubMedGoogle Scholar
- Galbreath KE, Cook JA: Genetic consequences of Pleistocene glaciations for the tundra vole (Microtus oeconomus) in Beringia. Mol Ecol. 2004, 13 (1): 135-148. 10.1046/j.1365-294X.2004.02026.x.View ArticlePubMedGoogle Scholar
- 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. 10.1016/j.ympev.2004.03.006.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- Zhou YL: Geography of the Vegetation in Northeast China. 1997, Beijing , Science PressGoogle Scholar
- 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.Google Scholar
- Jin CZ, Xu QQ, Zheng JJ: On the dispersal events of Mammuthus during the late Late Pleistocene. Vertebrate Palastatica. 1998, 36 (1): 47-53.Google Scholar
- Zhang RZ: Relict distribution of land vertebrates and Quaternary glaciation in China. Acta Zoologica Sinica. 2004, 50 (5): 841-851.Google Scholar
- 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.Google Scholar
- 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.Google Scholar
- Edwards RL, Cutler KB, Cheng H, Gallup CD: Geochemical Evidence for Quaternary Sea-level Changes. Treatise on Geochemistry. 2003, 6: 343-364.View ArticleGoogle Scholar
- 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. 10.1038/382241a0.View ArticleGoogle Scholar
- 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. 10.1046/j.1365-2699.1999.00346.x.View ArticleGoogle Scholar
- Liu JQ, Ni YY, Chu GQ: Main Palaeoclimatic events in the Quaternary. Quatern Sci. 2001, 21 (3): 239-248.Google Scholar
- 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.Google Scholar
- Remington CL: Suture-zones of hybrid interaction between recently joined biotas. Evol Biol. 1968, 2: 321-428.Google Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- 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. 10.1111/j.1365-294X.2006.03194.x.View ArticlePubMedGoogle Scholar
- 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. 10.1266/ggs.76.311.View ArticlePubMedGoogle Scholar
- 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. 10.1093/nar/25.24.4876.PubMed CentralView ArticlePubMedGoogle Scholar
- Swofford DL: PAUP* — Phylogenetic Analysis Using Parsimony (* and Other Methods. Beta version 4.0b10). Sunderland , Sinauer Associates.
- Posada D, Crandall KA: MODELTEST: testing the model of DNA substitution. Bioinformatics. 1998, 14 (9): 817-818. 10.1093/bioinformatics/14.9.817.View ArticlePubMedGoogle Scholar
- Felsenstein J: Confidence limits on phylogenetics: an approach using the bootstrap. Evolution. 1985, 39: 783-791. 10.2307/2408678.View ArticleGoogle Scholar
- Fu YX, Li WH: Statistical tests of neutrality of mutations. Genetics. 1993, 133 (3): 693-709.PubMed CentralPubMedGoogle Scholar
- Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123 (3): 585-595.PubMed CentralPubMedGoogle Scholar
- Rozas J, Rozas R: DnaSP version 3: an integrated program for molecular population genetics and molecular evolution analysis. Bioinformatics. 1999, 15 (2): 174-175. 10.1093/bioinformatics/15.2.174.View ArticlePubMedGoogle Scholar
- Kimura M: The neutral theory of molecular evolution. 1983, Cambridge, Massachusetts , Cambridge University PressView ArticleGoogle Scholar
- 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 CentralPubMedGoogle Scholar
- Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000, 9 (10): 1657-1659. 10.1046/j.1365-294x.2000.01020.x.View ArticlePubMedGoogle Scholar
- 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. 10.1046/j.1365-294x.2000.00887.x.View ArticlePubMedGoogle Scholar
- Templeton AR: Statistical phylogeography: methods of evaluating and minimizing inference errors. Mol Ecol. 2004, 13 (4): 789-809. 10.1046/j.1365-294X.2003.02041.x.View ArticlePubMedGoogle Scholar
- Schnerider S, Roessli D, Excoffier L: ARLEQUIN ver. 2.000 A Software for Population Genetic Data Analysis. 2000, Geneva , Genetics and Biometry Laboratory, University of GenevaGoogle Scholar
- Tajima F: The effect of change in population size on DNA polymorphism. Genetics. 1989, 123 (3): 597-601.PubMed CentralPubMedGoogle Scholar
- Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997, 147 (2): 915-925.PubMed CentralPubMedGoogle Scholar
- Nielsen R, Wakeley J: Distinguishing migration from isolation: A Markov chain Monte Carlo approach. Genetics. 2001, 158 (2): 885-896.PubMed CentralPubMedGoogle Scholar
- Yeh FC, Yang RC, Boyle T: POPGENE, the user friendly shareware for population genetic analysis. 1999, Edmonton, Canada , Molecular Biology and Biotechnology Center, University of AlbertaGoogle Scholar
- Rohlf FJ: NTSYS-pc, Numerical Taxonomy and Multivariate Analysis System, version 2.1. 2000, Setauket, New York , Exeter SoftwareGoogle Scholar