Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

The Phylogeography and Population Demography of the Yunnan Caecilian (Ichthyophis bannanicus): Massive Rivers as Barriers to Gene Flow

  • Hui Wang ,

    Contributed equally to this work with: Hui Wang, Xia Luo

    Affiliation Anhui Key Laboratory of Eco-engineering and Bio-technique, School of Life Sciences, Anhui University, Hefei, 230601, Anhui, China

  • Xia Luo ,

    Contributed equally to this work with: Hui Wang, Xia Luo

    Affiliation Anhui Key Laboratory of Eco-engineering and Bio-technique, School of Life Sciences, Anhui University, Hefei, 230601, Anhui, China

  • Shaoquan Meng,

    Affiliation College of Life Science & Technology, Yulin Normal University, Yulin, 537000, Guangxi, China

  • Yongjian Bei,

    Affiliation College of Life Science & Technology, Yulin Normal University, Yulin, 537000, Guangxi, China

  • Tao Song,

    Affiliation Anhui Key Laboratory of Eco-engineering and Bio-technique, School of Life Sciences, Anhui University, Hefei, 230601, Anhui, China

  • Tao Meng,

    Affiliation Guangxi Forestry Inventory and Planning Institute, Nanning, 530011, Guangxi, China

  • Guifen Li ,

    ylsylgf@163.com (GL); zhangbw@ahu.edu.cn (BZ)

    Affiliation College of Life Science & Technology, Yulin Normal University, Yulin, 537000, Guangxi, China

  • Baowei Zhang

    ylsylgf@163.com (GL); zhangbw@ahu.edu.cn (BZ)

    Affiliation Anhui Key Laboratory of Eco-engineering and Bio-technique, School of Life Sciences, Anhui University, Hefei, 230601, Anhui, China

Abstract

Ichthyophis bannanicus is the only caecilian species in China. In this study, the phylogeography and population demography of I. bannanicus were explored, based on the mitochondrial DNA genes (cyt b and ND2) and 15 polymorphic microsatellite loci. Altogether 158 individuals were collected from five populations in Yunnan province, Guangxi province, Guangdong province, and Northern Vietnam. Phylogeographical and population structure analysis identified either two groups (Xishuangbanna, Northern Vietnam-Yulin-Yangchun-Deqing) or three groups (Xishuangbanna, Northern Vietnam-Yulin-Yangchun, and Deqing), indicating that the Red River and Pearl River systems may have acted as gene-flow barriers for I. bannanicus. Historical population expansion that happened 15–17 Ka ago was detected for mtDNA data and was possibly triggered by warmer weather after the Last Glacial Maximum. However, the Bayesian simulations of population history based on microsatellite data pinpointed population decline in all populations since 19,123 to 1,029 years ago, demonstrating a significant influence of anthropogenic habitat alteration on I. bannanicus.

Introduction

Caecilians (Gymnophiona or Apoda) are the least-known group of the three orders of living amphibian on account of their subterranean habits [1]. They are readily distinguished from frogs and salamanders by their worm-like, annulated bodies with tentacles and reduced eyes [2]. This group is found across much of the wet tropics and some subtropical regions apart from Madagascar, Australasia, and Southeast Asia east of Wallace’s line [2]. Caecilian taxonomy has been problematic even at the family level due to the lack of morphological and molecular data [3,4]. Currently, approximately 190 caecilian species in 10 families are recognized in Amphibiaweb [5] and the number is rapidly increasing due to rigorous regional surveys. In the face of global amphibian decline, more and more attention is being paid to amphibian species. However, compared with frogs and salamanders, little is known about caecilians. Among amphibian population genetic studies published between 2001 and 2010, about 68% of the studies were about frogs and 32% about salamanders, and caecilians were not represented at all [6].

The Yunnan caecilian (Ichthyophis bannanicus) is the only species of Gymnophiona in China. It was first found in the Yunnan Province in 1976 and initially identified as Ichthyophis glutinosus [7]. In 1984, based on more samples found in the Guangxi Province, it was classified as a new species and called I. bannanicus [8]. In 1985, it was also found in the Guangdong Province [9]. Since then, a few sporadic new occurrences have been reported in China, all restricted to the above three provinces. Recently, I. bannanicus was also found in Vietnam [10], Thailand and Laos [3]. Given the few records, the distribution of I. bannanicus is still unclear, and even less is known about its population size. The first survey of national terrestrial wild animals organized by the State Forestry Administration of China in 1995 reported that the population size of I. bannanicus was approximately 10,000 individuals (http://www.china.com.cn/chinese/PI-c/603534.htm). A pessimistic view was presented by the China Red Data Book of Endangered Animals, indicating that the number of I. bannanicus was only 400 and that it should be classified as an endangered species in China [11]. However, follow-up field surveys suggested that the population size was far bigger than the numbers proposed earlier. For example, an investigation on the population size of I. bannanicus in Beiliu City (Guangxi Province, China) showed approximately 10,800 individuals living in a single village in spite of the obvious population decline in this area [12]. In addition, Chou Wenhao reported to the Global Amphibian Assessment (GAA) workshop that about 200 individuals were collected in a single night in the Xishuangbanna National Reserve in 2000 [1]. Now, I. bannanicus is classified as near threatened (NT) in the China Species Red List [13]. All these data demonstrated that at least in some areas I. bannanicus still maintained a relatively large population size. However, little was known about the population genetics studies of the species needed to inform conservation planning.

It is well known that the genetic structure and population history of existing wildlife were profoundly affected by landscape features and/or historical climate fluctuations. In the range of I. bannanicus, well-developed river systems are one of the dominant landscape features, including the Yangtze, Pearl, Red, Mekong, Salween and Irrawaddy rivers and their criss-crossed tributaries. With typical amphibian features, the life history of I. bannanicus consists of an aquatic larval stage of about two years and a subterranean post-metamorphic stage in which the animals rarely enter the water [14,15]. Therefore, the distribution and phylogeographic patterns of I. bannanicus are thought to be influenced by river systems. Earlier studies indicated that I. bannanicus is strongly influenced by low temperatures, and when the temperature drops below 18°C, hibernation and even death may occur [15]. It is generally acknowledged that glaciers were never present in Southern China [16,17]. However, the fluctuating temperatures and water levels accompanying glaciations and deglaciations may have greatly affected the migration and dispersal of I. bannanicus. Nowadays, human activities are generally thought to be one of the major causes of habitat fragmentation and population declines for many endangered species [1820]. Human impact may occur directly on species by hunting for instance, or indirectly by anthropogenic habitat changes. Human capture of caecilians was rarely reported due to their elusive behavior and lack of economic value. Therefore, anthropogenic habitat changes were regarded as the major cause of population decline or even local extinction for some caecilians [21,22]. On the other hand, certain agricultural practices such as crop plantations and fish farming are actually beneficial for caecilians, by providing suitable habitat (wet and loose soil), rather than having a negative impact. [2326]. Although a trend of population decline was observed in field surveys [12], the abundance of I. bannanicus in synanthropic habitats (such as farmland) was higher than in primary habitats. However, whilst we know much more about how caecilians can survive or thrive in disturbed habitats, little is known about demographic changes due to primary habitat loss and adaptation to new habitat [1]. Therefore, it is necessary to estimate the impact of anthropogenic habitat changes on I. bannanicus.

In this study, we used multiple molecular markers including two mitochondrial DNA encoded genes, cytochrome b (cyt b) and NADH dehydrogenase 2 (ND2), and fifteen nuclear DNA microsatellite loci aiming: I. to evaluate population genetic diversity of I. bannanicus; II. to investigate the effects of landscape features and historical climate fluctuations on genetic structure and population evolutionary history; III. to assess the impacts of anthropogenic habitat changes on population demography.

Materials and Methods

Sample collection and DNA extraction

One hundred and fifty-eight tissue samples of I. bannanicus were collected in southern China and northern Vietnam from 1998 to 2010 (Fig 1 and S1 Table). Samples were from seven localities in the five currently known distribution regions (defined here as populations): Xishuangbanna (BN, N = 34) from the Yunnan Province; Yulin (YL, N = 66) and Yangchun (YC, N = 19) from the Guangxi Province; Deqing (DQ, N = 25) from the Guangdong Province; Northern Vietnam (VN, N = 14) (S1 Table). Samples were preserved in 95% ethanol or stored at -20°C in the laboratory. Genomic DNA was isolated using the proteinase K digestion followed by SDS-phenol/chloroform method [27].

thumbnail
Fig 1. The study area and sampling localities of I. bannanicus.

1, Xishuangbanna, Yunnan province (BN); 2, Tam Dao, Vietnam (VN); 3–5, Yulin, Guangxi province (YL); 6, Yangchun, Guangxi province (YC); 7, Deqing, Guangdong province (DQ); S1, Tam Dao, Vietnam; S2, Laos; S3, Ubon Ratchathani Province, Thailand; S4, Phu Wan, Thailand; S5, Chiang Mai Province, Thailand; S6, Mae Hong Son, Thailand; S7, Longling, Yunnan, China.

https://doi.org/10.1371/journal.pone.0125770.g001

Ethical statement

In the present study, collection of samples was performed within a long-term investigation project on I. bannanicus and all samples were from individuals that died naturally and were found during field work. This investigation project and the sample collection were approved by the Forestry Administration of Guangdong, Guangxi and Yunnan provinces. The specimen collection in Vietnam was conducted under the ethics license from the Vietnam National Museum of Nature, and was accepted by the Animal Research Ethics Committee of Anhui University and Yulin Normal University. Our experimental procedures complied with the current laws on animal welfare and research in China, and were specifically approved by the Animal Research Ethics Committee of Anhui University and Yulin Normal University.

MtDNA sequencing and data analysis

Mitochondrial DNA fragments encoding cyt b and ND2 genes were amplified for 158 individuals. Amplification of the cyt b fragment was performed using the primers L14014 (5’-GGA TCT AAC CAA GAC TAA TGG TC-3’, forward) and H15198 (5’-GCT AAT GCT TTC TGA TAA GCT AC-3’, reverse) [28]. Amplification of the ND2 fragment was performed using the primers ND2-F (5’-CAC CTT CTA AGT AAG GTC AGC-3’) and ND2-R (5’-TCC TTC TTT GGG CTT TGA AGG-3’). All PCRs were performed with the same conditions in 50 μL: 20 to 80ng of genomic DNA, 25 μL 2×EasyTaq PCR SuperMix polymerase (TransGen Biotech, containing 1.25U Ex Taq, 0.4mM dNTP, 4mM Mg2+) and 0.4 μM of primers. Reactions were performed with the following profile: 5 min denaturing step at 95°C followed by 35 cycles of 30s at 95°C, 40s at 55°C and 40s at 72°C, and a final 10 min extension period at 72°C. PCR products were purified and subsequently sequenced on the ABI 3730 semi-automated Sequencer (PE Applied Biosystems). The cyt b DNA fragment was sequenced with the primers L14014 and L14500 (5’-GAG CCA CAG TAA TTA CAA ATC TC-3’), and ND2 was sequenced with the primers ND2-F and ND2-R. Nucleotide sequences were checked by eye and assembled using the program SeqMan in DNAStar [29], and open-reading frame for the two genes were verified. Cyt b and ND2 genes were combined together for analysis and aligned in MEGA 5.0 [30]. Nucleotide diversity (π, Nei 1987), haplotype diversity (h, Nei 1987) and haplotypes were identified using DnaSP v5.10.1 [31]. Genetic differentiation (Fst) between populations was calculated as implemented in Arlequin 3.5.1.3 [32] based on haplotype frequency differences with 10,000 permutations. The best-fitting nucleotide substitution model for sequence data was identified using the software MrModeltest 2.3 [33] and the GTR+G substitution model was selected. We estimated phylogeographic relationships using Bayesian Inference on the combined cyt b and ND2 haplotype sequences, implemented in BEAST 1.7.4. [34]. The homologous sequence of Ichthyophis glutinosus (AY456251.1) was used as outgroup. Two independent runs of Markov Chain Monte Carlo (MCMC) analyses for 30 million generations were conducted with sampling every 1,000 generations, and 10% of the initial samples were discarded as “burn-in”. The convergence of chains was confirmed until average standard deviation of split frequency was below 0.01, and Bayesian Posterior Probabilities (BPP) were used as an indicator of node credibility (≥ 95% was considered significant) [35]. A median-joining network was constructed to depict relationships among haplotypes using Network 4.6.1.2 [36]. In addition, an independent phylogenetic analysis was repeated, combining the haplotypes found in the present study with seven published cyt b sequences of I. bannanicus (S1–S7, S1 from Northern Vietnam, S2 from Laos, S3–S6 from Thailand and S7 from Yunnan, China, Fig 1 and S1 Table) [3,37].

Historical demographic inferences were carried out using the Bayesian Skyline Plot (BSP) as implemented in the software BEAST with GTR+G substitution model. Two independent runs of Markov Chain Monte Carlo (MCMC) analyses for 30 million generations were conducted with sampling every 1,000 generations and 10% of the initial samples were discarded as “burn-in”. Additionally, to scale the time axis on BSPs, we used data estimates for the most recent common ancestor (MRCA) of all haplotypes in a lineage. Lacking fossil evidence, we assumed a range of substitution rate of 0.50–1.00% per Myr for cyt b or ND2 based on evolutionary rates commonly proposed for amphibians [3843]. Plots for each analysis were visualized using Tracer 1.4 [44]. Tajima’s D [45] and Fu’s Fs [46] were calculated using Arlequin 3.5.1.3.

Microsatellite genotyping and data analysis

Fifteen polymorphic tetranucleotide microsatellite loci for I. bannanicus were screened in the present study [47]. PCR reaction was performed in 20 μL system containing 20–50 ng of template DNA, 10 μL 2×EasyTaq PCR SuperMix polymerase (TransGen Biotech, containing 1.25U Ex Taq, 0.4mM dNTP, 4mM Mg2+) and 0.2 μM of each primer. Thermal cycling was performed with the following profile: 5 min denaturing step at 95°C; followed by 35 cycles of 30s at 95°C, 1 min at optimal annealing temperature [47] and 1 min at 72°C; and a final extension at 72°C for 10 min. All PCR products were visualized on an ABI 3730 semi-automated sequencer (PE Applied Biosystems) with GS500 marker and analyzed using GENEMARKER 1.85 (version 1.3, SoftGenetics LLC).

Micro-Checker [48] was used to detect the presence of null alleles and genotyping errors in microsatellite genotyping. Genetic variation was assessed by summary statistics, including mean number of alleles per locus (MNA), observed heterozygosities (HO), expected heterozygosities (HE) and inbreeding coefficients (FIS), which were performed for each population using GENETIX 4.02 [49]. Deviation from Hardy-Weinberg equilibrium (HWE) across all loci for each population was assessed using an exact probability test implemented in GENEPOP 4.2.1 [50]. Significance values for multiple comparisons were adjusted using the Bonferroni correction [51].

Genetic differentiation between populations was quantified using FST values which were calculated in GENETIX. In contrast to the preceding analyses (based on region-defined populations), a Bayesian clustering approach was used to estimate the number of potential clusters present in the microsatellite data and to assign individuals to inferred clusters using the STRUCTURE software [52]. The posterior probability for different values of K (the number of population) tested was set from 1 to 8, and ten independent runs were performed at the length of 1,000,000 MCMC with a 100,000 burn-in period using no prior information and assuming admixture and correlated allele frequencies. The true K was selected using the maximal value of the log likelihood [Ln Pr(X/K)] of the posterior probability of the date for a given K [52]. Further, the statistic ΔK, the second-order rate of change in the log probability of the data between successive values of K was also derived [53]. The results were graphically displayed by the software DISTRUCT [54].

Demographic history based on microsatellites was assessed using the following methods. First, the Wilcoxon’s sign rank test [55] was used to test for heterozygosity excess under the stepwise mutation model (SMM) and two-phase mutation model (TPM), with 95% single step mutations and 5% multi-step mutations and a variance of 12 following the recommendations of Piry et al. [56]. Second, a mode-shift test [57] was carried out to detect any distortion of the expected L-shaped distribution of allele frequency. Both Wilcoxon’s sign rank test and mode-shift test were performed in BOTTLENECK 1.2.02 [56]. Finally, we inferred past demographic change of populations from microsatellite data using MSVAR 1.3 [58] with a coalescent simulation-based Bayesian likelihood method. This approach employs coalescent MCMC simulations to estimate the posterior probability distribution of population parameters, based on the observed distribution of microsatellite alleles and their repeat numbers. The main parameters from MSVAR are: current effective population size (N0), historical or ancestral effective population size (N1) and the time since the beginning of population change (T). Five independent simulations were run on subsamples of five populations. In every simulation, we ran each chain with a thinning interval of 10,000 steps, leading to a total number of Monte Carlo searches of 1×109 steps with the first 10% of total iterations discarded as burn-in. The remaining data were used to obtain the lower (5%), the median (50%), and the upper (95%) quantiles of the posterior distributions. Different means for the mean N1 were used to represent three demographic histories: a stable population (N0 = N1), an expanding population (N0>N1) and a decreasing population (N0<N1). The generation time of I. bannanicus was set to 5 years [59,60]. We estimated the marginal posterior distributions of the model parameters using the LOCFIT package [61] implemented in R v2.11.1 (R Development Core Team, 2010).

Results

Genetic diversity

A total of 2,181 base pairs of mtDNA (1038 for ND2 and 1,143 for cyt b) were obtained from 158 individuals after alignment. A total of 85 variable nucleotide sites were found with no insertions or deletions and 23 different haplotypes (H01-H23) were identified (S2 Table, GenBank accession numbers KM513544-KM513566). All haplotypes were restricted to a single population except for H08, which was the most common haplotype, and shared by four populations (VN/YL/YC/DQ) (S2 Table). Summary statistics, including haplotype richness, haplotype diversity, and nucleotide diversity within each population are given in Table 1. On the whole, a moderate genetic diversity was detected in I. bannanicus with an average haplotype diversity of h = 0.703 and a nucleotide diversity of π = 0.0077 (Table 1). The value of both haplotype diversity (h = 0.681) and nucleotide diversity (π = 0.0021) were remarkable higher for BN than the other populations (Table 1).

thumbnail
Table 1. Genetic variation and neutral test based on mitochondrial control region and fifteen microsatellite loci in five populations of I. bannanicus.

https://doi.org/10.1371/journal.pone.0125770.t001

For microsatellites, 158 individuals were successfully genotyped at fifteen loci. Microsatellite genotyping raw data are given in S3 Table. Micro-Checker indicated no evidence of null alleles nor genotyping errors such as large allele dropout or stuttering. The MNA for the five populations ranged from 3.64 to 8.29 with a mean value of 12.43. Observed heterozygosity (HO ± SD) ranged from 0.462 ± 0.17 to 0.666 ± 0.218, and expected heterozygosity (HE ± SD) ranged from 0.479 ± 0.255 to 0.715 ± 0.202. In accordance with mitochondrial data, the BN population had the highest MNA, HO and HE (Table 1). There was no deviation in HWE apart from five loci in the VN population prior to Bonferroni correction. There was no linkage disequilibrium at any locus in any population.

Phylogeography and divergence dating

The Bayesian tree revealed that there were two distinct haplotype groups with high bootstrap support: Clade A included seven haplotypes (H01–H07) all from the BN population, and Clade B included all the other sixteen haplotypes (H08–H23) (Fig 2). Within Clade A, H04 and H06 were divided from the other haplotypes and formed a distinct subclade. Clades A and B shared their most recent common ancestor (MRCA) about 1.39 million years before present (Ma BP) (95% HPD, 0.84–1.7 Ma). The phylogenetic tree of the above 23 haplotypes mixed with seven published cyt b sequences constructed with 714 homologous base pairs also showed two major lineages as above (S1 Fig). As for the seven new sequences, sample S1 from northern Vietnam joined Clade B and S2–S7 from Thailand and Laos merged into Clade A (S1 Fig).

thumbnail
Fig 2. The Bayesian tree of 23 haplotypes from mitochondrial data of I. bannanicus with I. glutinosus as outgroup.

Value above each branch refers to the Bayesian posterior probability (BPP).

https://doi.org/10.1371/journal.pone.0125770.g002

The median-joining network showed a clear pattern with obvious phylogeographic structure in accordance with the Bayesian tree (Fig 3). Two haplogroups were distinguished from the network analysis: Clade A was composed of H01–H07, and Clade B contained H08–H23. Within Clade A, H04 and H06 were far from the other five closely linked haplotypes. In Clade B, haplotypes clustered together closely around dominant haplotype H08 and presented a typical star-like distribution. H22 was another dominant haplotype in group B, containing most individuals (23) from the DQ population (25 individuals overall).

thumbnail
Fig 3. Median-joining network based on maximum parsimony among mtDNA haplotypes of I. bannanicus.

Nodes contain the haplotype name and are proportional to haplotype frequencies. Length of branches is proportional to the number of changes from one haplotype to the following, with a number next to the branch representing more than one mutation step. Black nodes indicate undetected intermediate haplotype states.

https://doi.org/10.1371/journal.pone.0125770.g003

Genetic structure

Genetic differentiation analysis produced similar results for mtDNA and microsatellite data: significant FST values (range 0.4548–0.5690 for mtDNA data and 0.3097–0.3358 for microsatellite data) were detected between BN and the other populations, and between DQ and all the other populations as well (range 0.5569–0.8964 for mtDNA data and 0.0878–0.3332 for microsatellite data). The FST values between VN, YL and YC ranged 0.0078–0.0554 for mtDNA (non-significant) and 0.0282–0.0811 for microsatellites, indicating weak to moderate genetic differentiation (Table 2). Based on STRUCTURE analysis, although no obvious maximum log likelihood of posterior probability was found [LnP (X/K) = −5325.2], the ΔK statistics output showed a clear maximum at K = 2 (Fig 4a and 4b). When K = 2, two clear clusters were identified: one cluster (Yellow) was composed of individuals from BN, and the other cluster (Red) of individuals from all the other populations (Fig 4c). The ΔK value based on the STRUCTURE output showed a large value for K = 3, too (Fig 4c). Correspondingly, a new independent cluster (Blue), composed of all individuals from the DQ population appeared, which may indicate possible substructuring within clade B (Fig 4c).

thumbnail
Table 2. Pairwise FST estimates based on mtDNA sequence (below diagonal) and microsatellite data (above diagonal), respectively.

https://doi.org/10.1371/journal.pone.0125770.t002

thumbnail
Fig 4. Bayesian STRUCTURE clustering results based on microsatellite genotypes among five I. bannanicus populations.

(a) The linear relationship between LnP(D) and K, (b) ΔK values as a function of K based on 5 runs, (c) STRUCTURE output at K = 2 and 3.

https://doi.org/10.1371/journal.pone.0125770.g004

Demography

Demographic analysis was carried out in five populations and in the two clades recovered by phylogenetic analysis (A and B). Significant negative Fu’s Fs test and Tajima’s D test values were detected in YL, YC and DQ populations, supporting population expansion (Table 1). For BN population (equivalent to Clade A), Tajima’s D test provided a significant negative value, but Fu’s Fs test showed conflicting results. Bayesian skyline plots (BSP) indicated that I. bannanicus maintained a relatively constant population size for quite a long time before a gradual decrease around 0.2 Ma BP, and then a strong growth occurred since approximate 16 thousand years before present (Ka BP) for Clade A (Fig 5). In contrast, Clade B maintained a constant population size before a strong growth approximately 14 Ka BP (Fig 5).

thumbnail
Fig 5. Bayesian skyline plot of effective population size.

Time is expressed in million years.

https://doi.org/10.1371/journal.pone.0125770.g005

Based on microsatellite data, BOTTLENECK indicated no significant signs of a recent bottleneck in all populations under both TPM and SMM. In addition, the mode-shift test showed a normal L-shaped distribution of allele frequencies in all populations. However, Bayesian coalescent simulation provided compelling evidence for population decline (Table 3, Fig 6). Summary statistics for the posterior distributions of N0, N1 and T were generated using MSVAR 1.3 with five independent replicates showing concordant results in all populations (Table 3, Fig 6). The BN population had the highest N0 (20366) and the YL population had the highest N1 (65574) while the VN population had the lowest N0 and N1 (170 and 8615). Overall, decline severity was low (N0/N1 ≥0.1) for BN, meanwhile moderate or severe population declines (N0/N1 ≥0.01) were detected in the VN, YL, YC and DQ populations. The medians of T for the five populations were 1029 (BN), 2140 (VN), 19123 (YL), 8182 (YC) and 9432 (DQ) respectively (Table 3, Fig 6).

thumbnail
Table 3. Medians of current (N0) and ancestral effective population sizes (N1) and time since population decline (T) for five populations based on a generation time of 5 years, using Bayesian MSVAR simulations.

https://doi.org/10.1371/journal.pone.0125770.t003

thumbnail
Fig 6. Estimated posterior distributions of current (N0, blue curve) and ancestral (N1, red curve) effective population sizes and time since population change (T, black curve) using MSVAR on a logarithmic scale based on a generation time of 5 years from five populations.

(a) BN, (b) VN, (c) YL, (d) YC and (e) DQ.

https://doi.org/10.1371/journal.pone.0125770.g006

Discussion

Genetic diversity

In the present study, genetic diversity of I. bannanicus was analyzed based on mtDNA and microsatellite data. Average, mtDNA haplotype diversity (h = 0.703) was higher than in some extremely endangered amphibians, such as Ranodon sibiricus (h = 0.367) [62], and lower or close to some non-endangered species, such as Rana nigromaculata (h = 0.4046–0.9454) [63], and Feirana quadranus (h = 0.6599–0.9498) [64]. Based on microsatellite data, the average MNA and HE were 12.43 and 0.674, respectively, suggesting high genetic diversity in I. bannanicus. Among the five studied populations, BN showed the highest genetic diversity based on either mtDNA or microsatellite data. The high ancestral effective population size N1 and low population decline severity of BN (Table 3, Fig 6) may explain its higher genetic diversity, and may be due to its larger ancestral population and longer population stablity.

Gene flow barriers

Up to now, little is known about gene flow barriers or dispersal pattern for caecilian species. In the present study, deep genetic divergences and remarkable differentiations were detected among the five populations (Table 2, Fig 4). These results may point to possible gene flow barriers for I. bannanicus. The geographic distance between DQ population and YL/YC population (191 km) is much less than the distance between YL/YC and VN population (530 km). However, a remarkable differentiation was detected between DQ and YL/YC/VN but not within YL/YC/VN. Considering that there is neither obvious elevation difference nor other geographic barriers between DQ and its adjacent populations except the Pearl River, we inferred that the Pearl River may constitute an effective geographic barrier between DQ and the other populations. Generally speaking, the effectiveness of rivers as barriers to gene flow is obvious in small mammals [65,66]. But for amphibians, the barrier effect of rivers is uncertain. For example, Li et al. found major rivers to be barriers for the alpine stream frog Scutiger boulengeri [67], and the Chagres River was shown to be a barrier for the túngara frog Physalaemus pustulosus [68]. In contrast, other studies have shown that mountain ridges rather than rivers constituted genetic barriers for amphibians [69,70]. Well-developed river systems are one of the dominant landscape features throughout the range of I. bannanicus, including Yangtze River, Pearl River, Red river, Mekong River, Salween River, Irrawaddy River and their crisscrossed tributaries. It seems that the Red River may act as an effective gene flow barrier between BN and other populations in this study. Because of insufficient data, we do not know whether other massive rivers such as the Mekong River or Salween River are also acting as gene flow barriers.

The conclusion that massive rivers may act as gene flow barriers was in conflict with our original hypothesis. We originally thought that the optimal migration path for caecilian species was through drainage systems because of weak movement ability in subterranean life stage [28]. It would appear that crossing a large river is very difficult for I. bannanicus.

However, the genetic barrier effect of rivers is more unstable than that of mountain ridges due to changing hydrological factors like water level, width, flow rate, bed height and channel position which directly affect the effectiveness of the barrier. For example, the barrier effect on the plateau wood frog Rana kukunoris was much stronger at the downstream of the Yalong River, where the water flow was higher and the river valley steeper and drier [71]. During historical periods, these hydrological factors, especially water levels, fluctuated with climate fluctuation and so were obviously affected by the Quaternary glaciations in the Pleistocene. In the glacial periods, temperatures lowered and massive sheets of ice locked water, resulting in low water levels in river systems, weakening the barrier effect. In China during glacial periods, the mean annual temperatures dropped by about 10–12°C, the sea level fell by about 140 m and the paleo-coastline pushed about 600 km away from the present coastline [16,17], and the water level of the Yangtze River fell by about 20–45 m during the Last Glacial Maximum [72]. During interglacial periods, the temperature rose and more water entered the rivers, resulting in wider and larger flowing rivers that strengthened their barrier effect. In the present study, the divergences between Clades A and B were estimated to have occurred about 1.39 Ma BP, corresponding to the Donau—Günz Interglacial that occurred 1.7–1.3 Ma BP [73] which supported the proposed hypothesis viewing the Red River as the main geographic barrier between the two I. bannanicus clades.

Demographic history

A strong historical population expansion signal was detected in both Clade A and B approximately 16 and 14 Ka BP respectively (Fig 5), which was just after the Last Glacial Maximum (LGM, 18–22 Ka BP). This post-LGM expansion event is not coincident with those of many Chinese species whose expansion events happened earlier than the LGM (e.g. [7479]), but consistent with species in North America and Europe (e.g. [8083]). We thought that this discordance was not about the geography or climate differences between I. bannanicus and other Chinese species. It may because that the I. bannanicus is more sensitive to low temperature than the other species [15]. Hence, we considered that the warmer weather after LGM might have triggered this rapid population expansion.

Although BOTTLENECK and the model-shift test did not support population declines, the MSVAR results proved strong population declines in all populations (Table 3, Fig 6). This inconsistent results can be ascribed to poorer ability of detecting neither too weak nor too recent population declines in BOTTLENECK than in MSVAR [84]. The times since population decline (T) were from 1,029 to 19, 123 years ago (Table 3, Fig 6), throughout the whole Neolithic Era when farming began. In southern China, rice cultivation appeared 10–12 Ka BP according to the fossil record [85,86]. Overall, the recent population declines of I. bannanicus could be attributed to human agricultural activities in the Neolithic Era. Farming in the Neolithic Era was mainly slash-and-burn cultivation which still survives in some minorities in southern China. Compared to modern farming, slash-and-burn had more severe impacts on the natural environment, and is still considered as an important threat in many areas nowadays [87]. Slash-and-burn practice and the resulting deforestation would have had various negative impacts on caecilians: first, individuals were swept away by the felling and fire because of their weak migration ability; second, because of poor soil fertility in slash-and-burned lands, the soil dwelling animals that the caecilians feed on would also decline [25,88]; third, because of the low soil fertility resulting from slash-and-burn agriculture, lands were abandoned after only a few years and new lands exploited, so spreading the negative impact. This land use pattern caused persistent and high-intensity habitat destruction for caecilians, thus explaining the population decline evidenced by MSVAR. For the BN population, the decline started later and was less severe than for the other populations (Table 3, Fig 6). Accordingly, in the early Holocene the whole of southern China was covered by evergreen broad-leaf forest [89]. However, at present, primeval forests remain only in Xishuangbanna (BN). We therefore suggest that these differences in the time and severity of decline are attributable to the less severe deforestation impact in Xishuangbanna than in other areas. Our data therefore strongly suggest that anthropogenic habitat alteration and deforestation led to I. bannanicus population declines with decline severity directly correlated to the degree of primary habitat destruction.

Fortunately, there is a view that current agricultural practices are not necessarily harmful to caecilians and may even be beneficial to them in some cases [2326]. Throughout the range of I. bannanicus, the staple crop is rice and mechanization is relatively low. The moist and loose soil in paddy fields may benefit I. bannanicus. In addition, some researchers reported that caecilians can tolerate agrochemicals [1,25]. Therefore, the current success of I. bannanicus in the agricultural environment may be just a remarkable recovery from historical population declines.

Supporting Information

S1 Fig. The Bayesian tree of 23 haplotypes from mitochondrial data combined with seven published cyt b sequences of I. bannanicus with I. glutinosus as outgroup.

Value above each branch refers to the Bayesian posterior probability (BPP).

https://doi.org/10.1371/journal.pone.0125770.s001

(EPS)

S1 Table. Samples used for analysis in this study together with the information on localities, coordinates, GenBank accession numbers for published sequences, populations and sample size.

https://doi.org/10.1371/journal.pone.0125770.s002

(DOC)

S2 Table. Twenty-three haplotypes identified in this study with information on numbers of individual (N) and sample numbers.

https://doi.org/10.1371/journal.pone.0125770.s003

(DOC)

S3 Table. Microsatellite genotyping raw data of fifteen loci. Abbreviations are the same as those in Fig 1.

https://doi.org/10.1371/journal.pone.0125770.s004

(XLS)

Acknowledgments

The authors would like to thank Prof. Shiping Gong, Prof. Jinping Chen (South China Institute of Endangered Animals, Guangzhou, China) and Nguyen Thien Tao (Vietnam National Museum of Nature, Hanoi, Vietnam) for their help in collecting samples. We would also like to thank Dr. John Bailey (University of Leicester) for his help on manuscript writing. We would also like to thank Tao Pan, Chaochao Hu and Ling Ding for their advice.

Author Contributions

Conceived and designed the experiments: BZ GL. Performed the experiments: HW XL. Analyzed the data: TS TM HW XL. Contributed reagents/materials/analysis tools: YB SM. Wrote the paper: HW XL.

References

  1. 1. Gower DJ, Wilkinson M (2005) Conservation biology of caecilian amphibians. Conserv Biol 19: 45–55.
  2. 2. Taylor EH (1968) The caecilians of the world: a taxonomic review. Kansas: University of Kansas Press.
  3. 3. Nishikawa K, Matsui M, Yong HS, Ahmad N, Yambun P, Belabut DM, et al. (2012) Molecular phylogeny and biogeography of caecilians from Southeast Asia (Amphibia, Gymnophiona, Ichthyophiidae), with special reference to high cryptic species diversity in Sundaland. Mol Phylogenet Evol 63: 714–723. pmid:22387289
  4. 4. Nussbaum RA, Wilkinson M (1989) On the classification and phylogeny of caecilian (Amphibia: Gymnophiona) a critical review. Herpetol Monogr 3: 1–42.
  5. 5. AmphibiaWeb: Information on amphibian biology and conservation. [web application]. (2014) California: AmphibiaWeb: Berkeley.
  6. 6. Emel SL, Storfer A (2012) A decade of amphibian population genetic studies: synthesis and recommendations. Conserv Genet 13: 1685–1689.
  7. 7. Yang D, Su C (1976) Ichthyophis glutinosus found from Xishuanbanna of Yunnan. Acta Zool Sin China 22: 175.
  8. 8. Yang D (1984) Finding and description of a new species Ichthyophis bannanicus. J Herpetologica China 3: 73–74.
  9. 9. Pan J, Liu C, Qian X (1985) Ichthyophis bannanicus found from Guandong. J South China Normal University (Natural Science Edition) China 1: 111–113.
  10. 10. Nguyen VS, Ho TC, Nguyen QT (2009) Herpetofauna of Vietnam. Frankfurt am Main: Edition Chimaira.
  11. 11. Zhao E (1998) China Red data Book of Endangered Animals (Amphibians and Reptiles). Beijing: Science Press.
  12. 12. Meng S, Bei Y, Li Y, Long S, Xie W, Li G (2006) An investigation on the quantity and distribution of Ichthyophis bannanicus in the town of Liuma in Beiliu city. J Yulin Teachers college (Natural Science) China 127: 114–117.
  13. 13. Wang S, Xie Y (2006) China Species Red List. Vol. 1. High Education Press.
  14. 14. Meng S, Li G (2006) Research Advances of Ichthyophis bannanicus. Sichuan J Zool China 25: 423–425.
  15. 15. Li G, Bei Y, Tan Y, Meng S, Xie W, Li J (2010) Observations of the reared caecilian, Ichthyophis bannanica. Sichuan J Zool China 29: 81–84.
  16. 16. Li J, Shu Q, Zhou S, Zhao Z, Zhang J (2004) Review and prospects of Quaternary Glaciation research in China. J Glaciol Geocryol China 26: 235–243.
  17. 17. Shi Y, Zhao J, Wang J (2011) New understanding of Quaternary Glaciations in China. Shanghai: Shanghai Popular Science Press.
  18. 18. Zhang B, Li M, Zhang Z, Goossens B, Zhu L, Zhang S, et al. (2007) Genetic viability and population history of the giant panda, putting an end to the "evolutionary dead end"? Mol Biol Evol 24: 1801–1810. pmid:17513881
  19. 19. Agudo R, Rico C, Vilà C, Hiraldo F, Donázar JA (2010) The role of humans in the diversification of a threatened island raptor. BMC Evol Biol 10: 384. pmid:21144015
  20. 20. Burney DA, Flannery TF (2005) Fifty millennia of catastrophic extinctions after human contact. Trends Ecol Evol 20: 395–401. pmid:16701402
  21. 21. Daniels RJR (2002) Some aspects of the threatened amphibians of India; Sharma BD, Kumari T, editors. New Delhi: Anmol Publications.
  22. 22. Kupfer A, Müeller H (2004) On the taxonomy of ichthyophiid caecilians from southern Thailand: a reevaluation of the holotype of Ichthyophis supachaii Taylor 1960 (Amphibia: Gymnophiona: Ichthyophiidae). Amphibia-Reptilia 25: 87–97.
  23. 23. Oommen OV, Measey GJ, Gower DJ, Wilkinson M (2000) Distribution and abundance of the caecilian Gegeneophis ramaswamii (Amphibia: Gymnophiona) in southern Kerala. Curr Sci India 79: 1386–1389.
  24. 24. Hofer D (2000) A short note about the status and abundance of caecilian populations. Froglog 42: 1.
  25. 25. Measey GJ (2004) Are caecilians rare? An East African perspective. J East African Nat Hist 93: 1–21.
  26. 26. Measey GJ, Gower DJ, Oommen OV, Wilkinson M (2003) Quantitative surveying of endogeic limbless vertebrates-a case study of Gegeneophis ramaswamii (Amphibia: Gymnophiona: Caeciliidae) in southern India. Appl Soil Ecol 23: 43–53.
  27. 27. Sambrook JEF, Fritsch EF, Maniatis T (1989) Molecular Cloning: A Laboratory Manual, 2nd edn. Cold Spring Harbor, New York: Cold Spring Harbor Laboratory Press.
  28. 28. Li G, Luo X, Bei Y, Meng S, Zhang B (2010) Genetic diversity and population structure of four populaitons of Ichthyophis bannanicus (Gymnophiona, Ichthyophiidae), based on mtDNA cytochrome b gene sequence varation. Acta Zootax Sin China 35: 494–503.
  29. 29. Burland TG (1999) DNASTAR's Laser gene sequence analysis software. Bioinformatics Methods and Protocols: Humana Press. pp 71–91.
  30. 30. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S (2011) MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol 28: 2731–2739. pmid:21546353
  31. 31. Librado P, Rozas J (2009) DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451–1452. pmid:19346325
  32. 32. Excoffier L, Lischer HE (2010) Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour 10: 564–567. pmid:21565059
  33. 33. Nylander JAA (2004) MrModeltest v2. Program distributed by the author. Evolutionary Biology Centre, Uppsala University.
  34. 34. Drummond AJ, Suchard MA, Xie D, Rambaut A (2012) Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol 29: 1969–1973. pmid:22367748
  35. 35. Leache AD, Reeder TW (2002) Molecular systematics of the eastern fence lizard (Sceloporus undulatus): a comparison of parsimony, likelihood, and Bayesian approaches. Syst Biol 51: 44–68. pmid:11943092
  36. 36. Bandelt HJ, Forster P, Röhl A (1999) Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol 16: 37–48. pmid:10331250
  37. 37. Gower DJ, Kupfer A, Oommen OV, Himstedt W, Nussbaum RA, Loader SP, et al. (2002) A molecular phylogeny of ichthyophiid caecilians (Amphibia: Gymnophiona: Ichthyophiidae): out of India or out of South East Asia? P Roy Soc B-Biol Sci 269: 1563–1569. pmid:12184826
  38. 38. Lin L, Qu Y, Li H, Zhou K, Ji X (2012) Genetic structure and demographic history should inform conservation: Chinese cobras currently treated as homogenous show population divergence. PloS one 7(4): e36334. pmid:22558439
  39. 39. Macey JR, Schulte JA II, Ananjeva NB, Larson A, Rastegar-Pouyani N, Shammakov SM, et al. (1998) Phylogenetic relationships among Agamid Lizards of the Laudakia caucasia species group: testing hypotheses of biogeographic fragmentation and an area cladogram for the Iranian Plateau. Mol Phylogenet Evol 10(1): 118–131. pmid:9751922
  40. 40. Macey JR, Schulte JA II, Larson A, Fang Z, Wang Y, Tuniyev BS, et al. (1998) Phylogenetic relationships of toads in the Bufo bufo species group from the eastern escarpment of the Tibetan Plateau: a case of vicariance and dispersal. Mol Phylogenet Evol 9(1): 80–87. pmid:9479697
  41. 41. Tan AM, Wake DB (1995) MtDNA phylogeography of the California newt, Taricha torosa (Caudata, Salamandridae). Mol Phylogenet Evol 4(4): 383–394. pmid:8747295
  42. 42. Weisrock DW, Macey JR, Ugurtas IH, Larson A, Papenfuss TJ (2001) Molecular phylogenetics and historical biogeography among salamandrids of the “True” salamander clade: rapid branching of numerous highly divergent lineages in Mertensiella luschani associated with the rise of Anatolia. Mol Phylogenet Evol 18(3): 434–448. pmid:11277635
  43. 43. Yan F, Zhou W, Zhao H, Yuan Z, Wang Y, Jiang K, et al. (2013) Geological events play a larger role than Pleistocene climatic fluctuations in driving the genetic structure of Quasipaa boulengeri (Anura: Dicroglossidae). Mol Ecol 22(4): 1120–1133. pmid:23216961
  44. 44. Rambaut A, Drummond A (2007) Tracer: a program for analysing results from Bayesian MCMC programs such as BEAST & MrBayes.
  45. 45. Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123: 585–595. pmid:2513255
  46. 46. Fu Y (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147: 915–925. pmid:9335623
  47. 47. Li G, Luo X, Bei Y, Wang H, Meng S, Zhang B (2010) Isolation and characterization of polymorphic tetranucleotide microsatellites in the Ichthyophis bannanicus (Gymnophiona: Ichthyophiidae). Conserv Genet Resour 2:129–131.
  48. 48. Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes 4: 535–538.
  49. 49. Belkhir K, Borsa P, Chikhi L (2001) 1996–2001 GENETIX 4.02, logiciel sous Windows TM pour la génétique des populations.
  50. 50. Rousset F (2008) GENEPOP'007: a complete re-implementation of the genepop software for Windows and Linux. Mol Ecol Resour 8: 103–106. pmid:21585727
  51. 51. Rice WR (1989) Analysing tables of statistical tests. Evolution 43: 223–225. pmid:2561774
  52. 52. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155: 945–959. pmid:10835412
  53. 53. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol 14: 2611–2620. pmid:15969739
  54. 54. Rosenberg NA (2004) DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes 4: 137–138.
  55. 55. Cornuet JM, Luikart G (1996) Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics 144: 2001–2014. pmid:8978083
  56. 56. Piry S, Luikart G, Cornuet JM (1999) BOTTLENECK: a computer program for detecting recent reductions in the effective size using allele frequency data. J Hered 90: 502–503.
  57. 57. Luikart G, Cornuet JM (1998) Empirical evaluation of a test for identifying recently bottlenecked populations from allele frequency data. Conserv Biol 12: 228–237.
  58. 58. Storz JF, Beaumont MA (2002) Testing for genetic evidence of population expansion and contraction: an empirical analysis of microsatellite DNA variation using a hierarchical Bayesian model. Evolution 56: 154–166. pmid:11913661
  59. 59. Kupfer A, Nabhitabhata J, Himstedt W (2004) Reproductive ecology of female caecilian amphibians (genus Ichthyophis): a baseline study. Biol J Linn Soc 83: 207–217.
  60. 60. Wen Y (1985) Observations on feeding habits of captive Ichthyophis glutinosus. Acta Herpetol Sin China 4: 147–148.
  61. 61. Loader C (2007) Locfit: Local regression, likelihood and density estimation. R package version: 1–5.
  62. 62. Chen S, Zhang Y, Wang X, Sun J, Xue Y, Peng Z, et al. (2012) Extremely low genetic diversity indicating the endangered status of Ranodon sibiricus (Amphibia: Caudata) and implications for phylogeography. PLoS One 7: e33378. pmid:22428037
  63. 63. Yang Y, Zhang D, Li Y, Ji Y (2004) Mitochondrial DNA diversity and preliminary biogeographic inference of the evolutionary history of the black—spotted pond frog Rana nigromaculata populations in China. Acta Zool Sin China 50: 193–201.
  64. 64. Wang B, Jiang J, Xie F, Li C (2012) Postglacial colonization of the Qinling Mountains: phylogeography of the Swelled Vent frog (Feirana quadranus). PloS one 7: e41579. pmid:22848532
  65. 65. Eriksson J, Hohmann G, Boesch C, Vigilant L (2004) Rivers influence the population genetic structure of bonobos (Pan paniscus). Mol Ecol 13: 3425–3435. pmid:15488001
  66. 66. Carmichael LE, Nagy JA, Larter NC, Strobeck C (2001) Prey specialization may influence patterns of gene flow in wolves of the Canadian Northwest. Mol Ecol 10: 2787–2798. pmid:11903892
  67. 67. Li R, Chen W, Tu L, Fu J (2009) Rivers as barriers for high elevation amphibians: a phylogeographic analysis of the alpine stream frog of the Hengduan Mountains. J Zool 277: 309–316.
  68. 68. Lampert KP, Rand AS, Mueller UG, Ryan MJ (2003) Fine-scale genetic pattern and evidence for sex-biased dispersal in the túngara frog, Physalaemus pustulosus. Mol Ecol 12: 3325–3334. pmid:14629349
  69. 69. Funk WC, Blouin MS, Corn PS, Maxell BA, Pilliod DS, Amish S, et al. (2005) Population structure of Columbia spotted frogs (Rana luteiventris) is strongly affected by the landscape. Mol Ecol 14: 483–496. pmid:15660939
  70. 70. Lougheed SC, Gascon C, Jones DA, Bogart JP, Boag PT (1999) Ridges and rivers: a test of competing hypotheses of Amazonian diversification using a dart-poison frog (Epipedobates femoralis). P Roy Soc B-Biol Sci 266: 1829–1835. pmid:10535104
  71. 71. Zhao S, Dai Q, Fu J (2009) Do rivers function as genetic barriers for the plateau wood frog at high elevations? J Zool 279: 270–276.
  72. 72. Yang D (1986) The paleoenvironment of the mid-lower regions of Yangtze River in the full-glacial period of late pleistocene. Acta Geogr Sin China 41: 302–310.
  73. 73. Kettle AJ, Heinrich D, Barrett JH, Benecke N, Locker A (2008) Past distributions of the European freshwater eel from archaeological and palaeontological evidence. Quaternary Sci Rev 27: 1309–1334.
  74. 74. Li Z, Yu G, Rao D, Yang J (2012) Phylogeography and demographic history of Babina pleuraden (Anura, Ranidae) in southwestern China. PloS one 7(3): e34013. pmid:22448286
  75. 75. Qu Y, Lei F, Zhang R, Lu X (2010) Comparative phylogeography of five avian species: implications for Pleistocene evolutionary history in the Qinghai-Tibetan plateau. Mol Ecol 19(2): 338–351. pmid:20002586
  76. 76. Shi W, Wang H, Zhu L, Zhu Q, Han D, Chang Q, et al. (2010) The genetic divergence and gene flow pattern of two muntjac deer (Muntiacus reevesi) populations, Wannan and Dabie Mountains, from the effect of Yangtze River and the late Pleistocene glacial oscillations. Acta Theriol Sin China 30(4): 390–399.
  77. 77. Song G, Qu Y, Yin Z, Li S, Liu N, Lei F (2009) Phylogeography of the Alcippe morrisonia (Aves: Timaliidae): long population history beyond late Pleistocene glaciations. BMC Evol Biol 9(1): 143.
  78. 78. Yu G., Zhang M, Rao D, Yang J (2013) Effect of pleistocene climatic oscillations on the phylogeography and demography of red knobby newt (Tylototriton shanjing) from Southwestern China. PloS one 8(2): e56066. pmid:23424644
  79. 79. Zhang D, Chen M, Murphy RW, Che J, Pang J, Hu J, et al. (2010) Genealogy and palaeodrainage basins in Yunnan Province: phylogeography of the Yunnan spiny frog, Nanorana yunnanensis (Dicroglossidae). Mol Ecol 19(16): 3406–3420. pmid:20666999
  80. 80. Avise JC (1998) Pleistocene phylogeographic effects on avian populations and the speciation process. P Roy Soc B-Biol Sci 265(1395): 457–463. pmid:9569664
  81. 81. Hewitt GM (1999) Post-glacial re-colonization of European biota. Biol J Linn Soc 68(1–2): 87–112.
  82. 82. Hewitt GM (2000) The genetic legacy of the Quaternary ice ages. Nature 405: 907–913. pmid:10879524
  83. 83. Hewitt GM (2004) Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc Lond B 359(1442): 183–195. pmid:15101575
  84. 84. Girod C, Vitalis R, Leblois R, Fréville H (2011) Inferring population decline and expansion from microsatellite data: a simulation-based evaluation of the Msvar method. Genetics 188: 165–179. pmid:21385729
  85. 85. Zhang W, Yuan J (1998) A preliminary study on the ancient rice excavated from Yuchanyan, Daoxian, Hunan Province. Acta Agron Sin China 24: 416–420.
  86. 86. Xiang A (2005) Archaeological research of prehistoric rice cultivation in Guangdong Province. Agr Archaeology China 1: 149–155.
  87. 87. Cornell JD (2011) Slash and burn. Available: http://www.eoearth.org/view/article/156045. Accessed 27 October 2011.
  88. 88. Deng X, Zou S, Fu X, Yao T, Sheng C, Bai Z (2003) The impacts of land use practices on the communities of soil fauna in the Xishuangbanna rainforest, Yunnan, China. Acta Ecol Sin China 23: 130–138.
  89. 89. Zheng Z, Deng Y, Zhang H, Yu B, Chen Z (2004) Holecene environmental changes in the tropical and subtropical areas of the south China and the relation to human activities. Quaternary Sciences China 24: 387–393.