- Research article
- Open Access
Phylogeography and Demographic History of Chinese Black-Spotted Frog Populations (Pelophylax nigromaculata): Evidence for Independent Refugia Expansion and Secondary Contact
BMC Evolutionary Biology volume 8, Article number: 21 (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 .
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.
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.
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).
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.
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
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.
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 . 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.
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 . 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 .
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).
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.
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
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 . 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 , 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 . 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.
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.
Hewitt GM: Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc. 1996, 58: 247-276.
Hewitt G: The genetic legacy of the Quaternary ice ages. Nature. 2000, 405 (6789): 907-913. 10.1038/35016000.
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.
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/07
Zink RM, Klicka J: The tempo of avian diversification: a comment on Johnson and Cicero. Evolution. 2006, 60 (2): 411-412. 2006/04/14
Cicero C, Johnson NK: The tempo of avian diversification: Reply. Evolution. 2006, 60 (2): 413-414.
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/17
Zhang RZ: Zoogeography of China. 1999, Beijing , Science Press
Zhang RZ: Geological events and mammalian distribution in China. Acta Zoologica Sinica. 2002, 48 (2): 141-153.
Liu DS, Li ZG: Geography of Asia. 1996, Beijing , Commercial Press
Fei L, Ye CY, Huang YZ: Key to Chinese Amphibia. 1990, Chongqing , Chongqing Branch, Science and Technology Literature Publishing House
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.
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.
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.
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.
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.
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.
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.
Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992, 9 (3): 552-569.
Avise JC: Phylogeography: The History and Formation of Species. 2000, Cambridge , Harvard University Press
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.
Gillespie A, Molnar P: Asynchronous maximum advances of mountain and continental glaciers. Rev Geophys. 1995, 33 (3): 311-364. 10.1029/95RG00995.
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.
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.
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.
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.
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.
Zhou YL: Geography of the Vegetation in Northeast China. 1997, Beijing , Science Press
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.
Jin CZ, Xu QQ, Zheng JJ: On the dispersal events of Mammuthus during the late Late Pleistocene. Vertebrate Palastatica. 1998, 36 (1): 47-53.
Zhang RZ: Relict distribution of land vertebrates and Quaternary glaciation in China. Acta Zoologica Sinica. 2004, 50 (5): 841-851.
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.
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.
Edwards RL, Cutler KB, Cheng H, Gallup CD: Geochemical Evidence for Quaternary Sea-level Changes. Treatise on Geochemistry. 2003, 6: 343-364.
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.
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.
Liu JQ, Ni YY, Chu GQ: Main Palaeoclimatic events in the Quaternary. Quatern Sci. 2001, 21 (3): 239-248.
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.
Remington CL: Suture-zones of hybrid interaction between recently joined biotas. Evol Biol. 1968, 2: 321-428.
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.
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.
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.
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.
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.
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.
Felsenstein J: Confidence limits on phylogenetics: an approach using the bootstrap. Evolution. 1985, 39: 783-791. 10.2307/2408678.
Fu YX, Li WH: Statistical tests of neutrality of mutations. Genetics. 1993, 133 (3): 693-709.
Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123 (3): 585-595.
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.
Kimura M: The neutral theory of molecular evolution. 1983, Cambridge, Massachusetts , Cambridge University Press
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.
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.
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.
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.
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 Geneva
Tajima F: The effect of change in population size on DNA polymorphism. Genetics. 1989, 123 (3): 597-601.
Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997, 147 (2): 915-925.
Nielsen R, Wakeley J: Distinguishing migration from isolation: A Markov chain Monte Carlo approach. Genetics. 2001, 158 (2): 885-896.
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 Alberta
Rohlf FJ: NTSYS-pc, Numerical Taxonomy and Multivariate Analysis System, version 2.1. 2000, Setauket, New York , Exeter Software
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.
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.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
About this article
Cite this article
Zhang, H., Yan, J., Zhang, G. et al. Phylogeography and Demographic History of Chinese Black-Spotted Frog Populations (Pelophylax nigromaculata): Evidence for Independent Refugia Expansion and Secondary Contact. BMC Evol Biol 8, 21 (2008). https://doi.org/10.1186/1471-2148-8-21
- Inter Simple Sequence Repeat
- Last Glacial Maximum
- Mismatch Distribution
- Interglacial Period
- Secondary Contact