Cryptic Genetic Diversity Is Paramount in Small-Bodied Amphibians of the Genus Euparkerella (Anura: Craugastoridae) Endemic to the Brazilian Atlantic Forest

Abstract

Morphological similarity associated to restricted distributions and low dispersal abilities make the direct developing “Terrarana” frogs of the genus Euparkerella a good model for examining diversification processes. We here infer phylogenetic relationships within the genus Euparkerella, using DNA sequence data from one mitochondrial and four nuclear genes coupled with traditional Bayesian phylogenetic reconstruction approaches and more recent coalescent methods of species tree inference. We also used Bayesian clustering analysis and a recent Bayesian coalescent-based approach specifically to infer species delimitation. The analysis of 39 individuals from the four known Euparkerella species uncovered high levels of genetic diversity, especially within the two previously morphologically-defined E. cochranae and E. brasiliensis. Within these species, the gene trees at five independent loci and trees from combined data (concatenated dataset and the species tree) uncovered six deeply diverged and geographically coherent evolutionary units, which may have diverged between the Miocene and the Pleistocene. These six units were also uncovered in the Bayesian clustering analysis, and supported by the Bayesian coalescent-based species delimitation (BPP), and Genealogical Sorting Index (GSI), providing thus strong evidence for underestimation of the current levels of diversity within Euparkerella. The cryptic diversity now uncovered opens new opportunities to examine the origins and maintenance of microendemism in the context of spatial heterogeneity and/or human induced fragmentation of the highly threatened Brazilian Atlantic forest hotspot.

About the Authors:

Luciana A. Fusinatto

* E-mail: [email protected]

Affiliation: Departamento de Ecologia, Instituto de Biologia Roberto Alcantara Gomes, Universidade do Estado do Rio de Janeiro, Rio de Janeiro, Rio de Janeiro, Brasil

João Alexandrino

Affiliations Departamento de Ciências Biológicas, UNIFESP - Universidade Federal de São Paulo, Campus Diadema, Diadema, São Paulo, Brasil, Departamento de Zoologia, Instituto de Biociências, Universidade Estadual Paulista, Rio Claro, São Paulo, Brasil

Célio F. B. Haddad

Affiliation: Departamento de Zoologia, Instituto de Biociências, Universidade Estadual Paulista, Rio Claro, São Paulo, Brasil

Tuliana O. Brunes

Affiliations CIBIO/UP, Centro de Investigação em Biodiversidade e Recursos Genéticos da Universidade do Porto, Vairão, Portugal, Departamento de Biologia, Faculdade de Ciências da Universidade do Porto, Porto, Portugal

Carlos F. D. Rocha

Affiliation: Departamento de Ecologia, Instituto de Biologia Roberto Alcantara Gomes, Universidade do Estado do Rio de Janeiro, Rio de Janeiro, Rio de Janeiro, Brasil

Fernando Sequeira

Affiliation: CIBIO/UP, Centro de Investigação em Biodiversidade e Recursos Genéticos da Universidade do Porto, Vairão, Portugal

Introduction

The Neotropics harbor one of the highest levels of biodiversity on earth [1], but only a fraction of this biodiversity has been described [2]. Amphibians are among the vertebrates in which species discovery has increased considerably in the last decades [3–7]. The search for new species in poorly known regions, the use of molecular data and of integrative approaches (combining natural history with genetic, morphological, and ecological data) all have contributed to increasing knowledge on amphibian biodiversity [8–12]. Species discovery has often been accompanied by challenging systematics and problematic taxonomies for several groups of amphibians. Morphological conservatism often associated to highly homoplastic traits is commonly mentioned as the cause for discordance between phylogenetic relationships inferred from molecular and morphological data [13–15]. The case of the Neotropical direct developing anurans of the taxon Terrarana is paradigmatic of this taxonomic complexity [5]. The taxonomy of this clade, composed by almost a thousand species distributed in four families (Ceuthomantidae, Brachycephalidae, Eleutherodactylidae and Craugastoridae) [16], has been controversial due to scarcity of diagnostic traits and high levels of trait plasticity and convergence [5,16–22]. Recent phylogenetic reconstructions with multimarker DNA data examined the general taxonomy of Terrarana [5,16,19], but the evolutionary history and phylogenetic relationships of several clades within this taxon remain largely unknown.

The genus Euparkerella (Craugastoridae) is one of the least known of the Terrarana [16,19]. It currently consists of four morphologically-based species - E. brasiliensis (Parker, 1926) (Figure 1), E. cochranae Izecksohn, 1988, E. robusta Izecksohn, 1988, and E. tridactyla Izecksohn, 1988. These species, distinguishable only by cryptic morphological traits [23], are restricted to small areas of Atlantic forest in southeastern Brazil [23,24]. One of the most impressive characteristic of Euparkerella species is their small body size (maximum SVL around 22 mm) with extremely reduced digits compared to those of their larger relatives [19], features that are commonly associated with the process of miniaturization [26]. The reduction of adult body size is a recurrent tendency in the evolution of amphibians, being strongly associated with terrestrial life history [25,26]. Small body size and associated traits (e.g., digit reduction) were recently suggested to be innovative morphological features resulting from microhabitat adaptation, that would reduce the vagility and the physiological tolerance of species, ultimately concurring to geographic range restriction and fragmentation, and higher population structure and diversification rates [26,27]. So, morphological similarity and restricted distributions with predicted low dispersal abilities make the genus Euparkerella a good model for examining diversification processes.

[Figure omitted. See PDF.]

Figure 1. Euparkerella brasiliensis. From Cachoeiras de Macacu, Rio de Janeiro State. Locality 5d. (Photo: L. A. Fusinatto) .

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

In the last decades several new genera and species of amphibians were recognized and described for the Neotropics as a result of increased sampling effort and taxonomic reevaluations using molecular data [5,28–30], but knowledge remains very poor on within genus species genetic diversification and evolutionary history. We here contribute to reduce this knowledge gap by examining patterns and levels of genetic diversity and inferring the phylogenetic relationships within the genus Euparkerella. To this purpose, we use DNA sequence data from both mitochondrial and nuclear genes coupled with traditional Bayesian phylogenetic reconstruction methods and more recent coalescent approaches of species tree inference. In addition, we investigate species-level biological diversity applying several methods [31] to delimit putative novel evolutionary units within the most extensively sampled Euparkerella species (E. brasiliensis, E. cochranae) using Bayesian cluster assignment analyses, the Genealogical Sorting Index, and in particular a recent Bayesian coalescent species delimitation approach.

Materials and Methods

Ethics Statement and Sampling

We analyzed 39 individuals from all the four known Euparkerella species, from 15 localities (28 collected in this study and 11 from scientific collections; Fig. 2; Table S1 in file S1), and one sample from the outgroup species Barycholos ternetzi [19]. All field collections did not involve endangered or protected species, and were obtained under appropriate permits (Instituto Chico Mendes de Conservação da Biodiversidade - ICMBio, permit number 18887; and Instituto Estadual do Ambiente - INEA, permit number 040/2010). Techniques used to capture, tissue sampling and euthanasia sought to minimize animal suffering and were in accordance with recommendations of the Herpetological Animal Care and Use Committee (HACC) of the American Society of Ichthyologists and Herpetologists (available at: http://www.asih.org/publications), as well as were accepted by ICMBio and INEA. When collected, individuals were euthanized using an anesthetic application over the skin (5 % Lidocaine), whereas other samples were obtained by toe clipping and followed by specimen release in the field. Prior to toe clipping, a suitable level of anesthetic was applied in the local of incision (2% Lidocaine).

[Figure omitted. See PDF.]

Figure 2. Samples localities of Euparkerella.

Numbers corresponds to populations of: Euparkerella tridactyla (1); E. robusta (2); E. cochranae (3,4); E. brasiliensis (5–7) and Euparkerella sp. (8).

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

Following the morphological criteria proposed by Izecksohn [23], we analyzed two specimens of E. tridactyla, one of E. robusta, 10 of E. cochranae, and 22 of E. brasiliensis. Since the four individuals from Guapimirim, RJ (Locality 8, see Figure 2 and Table S1 in File S1), were not fully diagnosable as E. brasiliensis or E. cochranae, we designated them as Euparkerella sp.

DNA amplification, sequencing and polymorphism

Whole genomic DNA was extracted from tissue samples using the QUIAGEN DNeasy Tissue kit (QUIAGEN). We amplified one mitochondrial gene fragment, cytochrome c oxidase subunit 1 (CO-I), and the nuclear β- fibrinogen intron 7 gene (β-fibint7), a segment with exon 2 and intron 2 of the cellular myelocytomatosis gene (C-myc2), the recombination activating protein 1 gene (RAG-1), and a segment of the exon 1 of tyrosinase gene (TYR). Protocols for DNA extraction, PCR (Polymerase Chain Reaction) gene amplification and sequencing were described in Supporting Information (Protocol S2), and primers used in PCR reactions are listed in Table S2 in File S1 and described in the literature [32–36]. Sequences were edited with the program ChromasPro v. 1.5 (Technelysium Pyt Ltd.). Sequence data were deposited in GenBank and are available under accession numbers KF625052 - KF625166 (Table S1 in File S1), except for β-fibint7 for which amplified sequences were shorter than 200 bp (Table 1), the limit size accepted for this database. Only for this gene, raw sequences are available from the Dryad Digital Repository: http://doi.org/10.5061/dryad.690jk. Alignments were made and manually edited in BioEdit Sequence Alignment Editor [37]. For nuclear genes, the polymorphic nucleotide positions of heterozygous individuals were coded with the IUPAC ambiguity codes. For each heterozygous genotype sequence, we inferred the most likely phased haplotypes with the Bayesian algorithm implemented in PHASE v2.1.1[38]. All known haplotypes were incorporated for haplotype inference. We ran PHASE three times with different random seeds and checked if haplotype estimation was consistent across runs. Each run was conducted using default values. Standard genetic diversity indices were estimated using the program DNAsp v. 5.10.1[39]. Genetic divergence estimated from mtDNA sequences within groups and between pairs of populations were measured by the p-uncorrected distance as implemented in the software MEGA v. 5.1 [40].

Fragment Size (bp) N S h %Hd (SD) %π (SD)

CO-I 543 34 194 14 92 (2.3) 11 (0.9)

β-fibint7 157 - 159 70 24 30 97 (0.6) 3.5 (0.1)

C-myc2 458 - 473 70 45 29 94 (1.6) 1.7 (0.1)

RAG-1 824 42 38 14 93 (1.5) 0.8 (0.1)

TYR 507 42 46 25 97 (1.1) 1.9 (0.2)

Table 1. DNA Polymorphism for five loci used in this study.

N – number of sequences (phased in nuclear genes); S – number of polymorphic positions; h – number of haplotypes; %Hd – haplotype diversity percentage and respective standard deviation (SD); %π nucleotide diversity percentage and respective standard deviation (SD).

CSV

Download CSV

Phylogenetic analysis and divergence time estimates

Phylogenetic analyses were based on independent gene tree inference for the five gene fragments, on traditional multilocus tree inference using sequences from all genes concatenated into a single data matrix, and on the coalescent-based approach implemented in *BEAST, where genes trees are simultaneously estimated to find the best multilocus species tree. The most appropriate model of nucleotide evolution and the best-fitting partitioning scheme were selected using PartitionFinder v. 1.0 [41] under the Akaike information criterion [42] (for adopted models, see Table S3 in File S1). For each partition, we set the number of distinct relative substitution rates and the rate of variation among sites (gamma distribution and/or proportion of invariable sites) that best adjusted to the nucleotide evolution model. The prior on the substitution rate was set as variable to accommodate rate differences across partitions.

The best partition strategy and the most appropriate model of nucleotide evolution (Table S3 in File S1) were used for Bayesian inference with the program MrBayes v.3.1.2 [43], for both gene trees and concatenated multilocus tree estimates. Four MCMC (Markov chain Monte Carlo) were run simultaneously in each analysis, each of one with four chains sampled every 1000 generations, until split frequencies reached 0.002. We removed the first 25% of trees as burn-in and tested convergence of chains plotting log-likelihood values versus number of generations in program Tracer v.1.5 [44]. The species tree was generated in *BEAST [45] implemented in BEAST v.1.7.4 [46]. A priori assignment of individual alleles to a “species”, was based on information from well-supported and geographically structured clades in the multilocus tree (concatenated dataset), the Bayesian population structure analysis, and genealogical indices of population divergence (see the results section).

To estimate divergence times to the most recent common ancestral (tMRCA) we used the mean nucleotide substitution rates of RAG-1 and TYR obtained by Heinicke et al. [18] for Terrarana (0.0015 mutations/site /million years). Although divergence times inferred on the basis of nucleotide substitution rates are not free of drawbacks, and are susceptible to errors due to substitution pattern variability across organisms and genes [47], the lack of independent calibration points (e. g., fossil records, biogeographic events) for the studied group prevents the application of an alternative approach. Although divergence times should be interpreted with caution, we consider the use of mutation rates still valuable to infer a timeframe for hypotheses of historical biogeography and patterns of diversification within Euparkerella. Models of nucleotide evolution were simplified with respect to codon positions to avoid overparameterization. Divergence times were estimated using a lognormal relaxed clock. We adopted the Yule speciation process and followed default settings for other parameters. Due to computational limitations, data were analyzed by four independent runs of 50 million generations (total of 200 million) sampled each 5000 generations, and 10% removed as burn-in. We assessed chain convergence with the program Tracer v. 1.5 [44], assembled different runs outputs with LogCombiner v. 1.7.4 [46] (Drummond et al. 2012), and generated the final tree with TreeAnnotator v. 1.7.4 [46].

Delimitation of genetic evolutionary units

In addition to phylogenetic analysis we used several methods to assess taxonomic diversity within the morphology-based species from Rio de Janeiro (E. brasiliensis, E. cochranae, and Euparkerella sp.): 1) cluster assignment analyses based on allele frequencies; 2) Bayesian coalescent-based species delimitation, and; 3) Genealogical Sorting Index.

Genetic structure analysis were investigated with no a-priori assignment of individuals to clusters using a Bayesian model-based algorithm implemented in the program STRUCTURE v. 2.3.3 [48] and examined independently for nuclear markers only and for mtDNA combined with nuclear markers. Although the use of a Bayesian clustering analysis in our system is not free of drawbacks, especially regarding the low sample size for each of the three entities herein analyzed, which may decrease the power of reliably inferring the optimal number of clusters [49] we still consider the use of such approach as useful for non-phylogenetic discovery tool of species-level diversity and for comparative purposes across different type of methods [50]. Sequence data were converted to Structure input file format using the program xmfa2struct (available at: http://www.xavierdidelot.xtreemhost.com/clonalframe.htm). The STRUCTURE analysis was performed under the admixture ancestry model, with five independent runs for each K ranging from 1 to 10, using 5x10⁵ MCMC repetitions and discarding the first 5x10⁴ iterations as burn-in. To choose the K value that best fitted our data we used the on-line program Structure Harvester v.0.6.93 [51] to monitor the estimated log posterior probability of the data (ln Pr (X/K)) [48], and estimate the second-order rate of change of the likelihood function (ΔK) [49]. We used the program CLUMPP v.1.1.2 [52]. to assemble results of the five independent runs generated by Structure in an average admixture matrix Q, which was represented graphically with the program DISTRUCT v. 1.1 [53].

.Delimitation of putative evolutionary units as inferred by phylogenetic and Bayesian clustering analyses were tested using a multilocus, coalescent-based species delimitation method implemented in the program “Bayesian Phylogenetics and Phylogeography” (BP & P; [54]). This method, that assumes that shared polymorphism is attributable to incomplete sorting of ancestral polymorphism, combines species phylogenies and gene genealogies via ancestral coalescent processes, providing statistical support (posterior probabilities) for different species delimitations. Following most users of BP&P [50,55,56] we used the species tree estimated with *BEAST as guide tree. Since the guide tree may play a critical role in the outcome of the species delimitation model, phylogenetic uncertainties on its topology may constitute an important drawback for accurate species delimitations [50,55]. Given the low nodal support (pp < 0.80) for the relationships among some of the evolutionary lineages recovered by our species tree estimates (E. cochranae 3 and E. brasiliensis 5, 6 and 7, see results), we followed the recommendation of [55], comparing the results from the use of different guide trees representing all possible phylogenetic resolutions of those evolutionary lineages. So, for the BP&P analysis we considered three different guide trees: (E. robusta 2, (E. tridactyla 1, ((Euparkerella sp. 8 + E. cochranae 4), ((E. brasiliensis 7 + E. cochranae 3), (E. brasiliensis 6 + E. brasiliensis 5))))); (E. robusta 2, (E. tridactyla 1, ((Euparkerella sp. 8 + E. cochranae 4), (E. cochranae 3, (E. brasiliensis 7, (E. brasiliensis 6 + E. brasiliensis 5)))))); and (E. robusta 2, (E. tridactyla 1, ((Euparkerella sp. 8 + E. cochranae 4), (E. brasiliensis 7, (E. cochranae 3, (E. brasiliensis 6 + E. brasiliensis 5)))))). We followed the same methodology used by [55] for adjusts of gamma prior (G) in population size (θ) and age of the root in the species tree (τ₀) parameters. Accordingly, we tested three scenarios to examine the influence of priors in the posterior probabilities: i) assuming relatively large ancestral population sizes and deep divergences (θ = 1, 10, τ₀ = 1, 10); ii) assuming relatively small ancestral population sizes and shallow divergences among species (θ = 2, 2000, τ₀ = 2, 2000); and iii) assuming large ancestral population sizes and relatively shallow divergences among species (θ = 2, 2000, τ₀ = 2, 2000). We ran two rjMCMC analyses for each trial using 500000 generations, 10000 as burn-in, with sample interval of five.

We also used the Genealogical Sorting Index method (gsi) [57] to test the level of genealogical divergence in our nuclear gene trees for each of the evolutionary units inferred by our multilocus phylogenetic analysis and Bayesian clustering analyses within E. cochranae, E. brasiliensis and Euparkerella sp.. This index allows us to take a measure of the degree of exclusive ancestry for single locus gene trees in pre-defined groups [57,58]. To estimate the gsi, we first generated rooted trees topologies for each individual nuclear locus. For each, we generated 1000 trees using the Maximum Likelihood (ML) approach implemented in the program RAxMLGUI v.1.2 [59], a graphic interface of RAxML [60], adopting the GTR gamma as model of nucleotide evolution. Subsequently, we generated a Majority Rule consensus tree keeping branches with posterior probabilities higher than 50%. Those 1000 individual gsi measurements were used to estimate an ensemble gsi statistic (gsi_T) for each locus. All values and their significance were obtained at the web server of Genealogical Sorting Index (available in: http://www.genealogicalsorting.org/).

All alignments and main input files used in this work are available from the Dryad Digital Repository: http://doi.org/10.5061/dryad.690jk

Results

Sequence variation

Standard genetic diversity summary statistics across loci are summarized in Table 1. For the mtDNA gene (CO-I), we obtained an alignment of 543 base pair sequences, corresponding to 14 unique haplotypes and 194 segregating sites (Table 1). Neither stop nor nonsense codons were detected. Within the genus Euparkerella, the pairwise mtDNA sequence divergence ranged from 26.6%, between E. tridactyla and E. robusta to 6.1 %, between Euparkerella. sp. 8 and E. cochranae 4; whereas within groups p-distance ranged from 0 % to 2.1 % (Table 2). Within E. brasiliensis, the sequence divergence ranged from 6.3% to 11.7%, while within E. cochranae this divergence was 13.6% (Table 2).

N (1) (2) (3) (4) (5) (6) (7) (8)

(1) E. tridactyla 2 0.0 125 285 365 320 420 440 450

(2) E. robusta 1 20.63 -- 160 240 190 300 315 225

(±1.73)

(3) E. cochranae 2 18.78 18.05 0.0 105 60 155 170 90

(±1.66) (±1.56)

(4) E. cochranae 7 19.71 17.22 13.63 2.1 50 60 75 15

(±1.61) (±1.54) (±1.41)

(5) E. brasiliensis 9 17.82 16.94 9.30 12.48 0.6 105 120 35

(±1.60) (±1.51) (±1.13) (±1.30)

(6) E. brasiliensis 6 19.21 16.08 9.64 12.89 6.26 0.5 15 75

(±1.71) (±1.49) (±1.19) (±1.35) (±0.95)

(7) E. brasiliensis 5 20.44 18.78 11.23 14.55 11.23 11.66 0.0 90

(±1.78) (±1.66) (±1.31) (±1.39) (±1.33) (±1.37)

(8) Euparkerella sp. 4 20.44 17.68 13.63 6.08 12.38 13.32 14.55 0.0

(±1.71) (±1.60) (±1.46) (±0.91) (±1.35) (±1.42) (±1.46)

Table 2. Genetic divergence (%) within and among populations, and geographical distance among populations/species of Euparkerella.

Genetic divergence estimated for the mitochondrial gene fragment CO-I. Uncorrected p-distance (%) and standard error among and within populations/species below and in the the diagonal (underlined), respectively Geographic distance in kilometers above the diagonal. N, number of specimens analyzed.

CSV

Download CSV

For nuclear loci, we obtained fragments of 157–159 bp for β-fibint7, 458–473 bp for C-myc2, 824 bp for RAG-1, and 507 for TYR (Table 1). The number of haplotypes ranged from 14 for RAG-1, to 30 for β-fibint7 (Table 1).

Gene trees

The Bayesian inference of the mtDNA CO-I gene tree resulted in an overall highly resolved and strongly supported topology for all major relationships (Figure 3). Contrasting with the previously morphologically-defined species E. robusta and E. tridactyla, both E. brasiliensis and E. cochranae were not recovered as monophyletic groups. Euparkerella cochranae 3 (from locality 3) and E. brasiliensis (from localities 5, 6 and 7) form a well-supported clade, and E. cochranae 4 (from locality 4) together with individuals belonging to Euparkerella sp. (locality 8) form another clade with high posterior probability. The former clade appears to be highly substructured with the occurrence of two sub-clades with high posterior probability, corresponding to individuals of E. brasiliensis from localities 5 and 6, but their relationships are basically unresolved (Figure 3). Although resolution varied among nuclear loci, gene trees showed an overall lower resolution when compared to the mtDNA tree, and some phylogenetic non-concordance (Figure 4). The placement of E. cochranae 3 was non-concordant across markers, by clustering with E. cochranae 4 and E. robusta on the β-fibint7 tree, but with its position unresolved for other nuclear markers. Another non-concordant pattern between mtDNA and nuclear gene trees was the monophyly of Euparkerella sp., which was recovered by all nuclear, but not mtDNA, gene trees (Figure 4).

[Figure omitted. See PDF.]

Figure 3. Mitochondrial (CO-I) gene tree of Euparkerella.

Bayesian phylogenetic inference of mitochondrial haplotypes of Euparkerella. Colors refer to populations and combinations of numbers-letters indicate localities (Figure 2 and Table S1 in File S1) corresponding to haplotypes. Posterior probabilities are indicated to the left of nodes. Asterisks represents posteriors equal or higher than 0.95.

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

[Figure omitted. See PDF.]

Figure 4. Nuclear genes trees of Euparkerella.

Bayesian phylogenetic inferences of nuclear haplotypes of four gene fragments (β-fibint7, C-myc, RAG-1 and TYR) of Euparkerella. Colors refer to populations and combinations of numbers-letters indicate localities (Figure 2 and Table S1 in File S1) corresponding to haplotypes. Posterior probabilities are indicated left to nodes. Asterisks represents posteriors equal or higher than 0.95.

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

Multigene trees and divergence time estimates

Both the multilocus tree resulting from concatenated data analysis and the species tree showed very similar topologies, but with better resolution and higher level of node support for the concatenated multilocus tree (Figure 5). Both E. brasiliensis and E. cochranae grouped in two different clades, a pattern also seen on the mtDNA tree. However, only the mtDNA and the concatenated multilocus tree showed high posterior probability for the clade formed by all E. brasiliensis individuals (localities 5, 6 and 7) and E. cochranae 3. Estimates of tMRCA indicated that the diversification of Euparkerella started in Miocene (11.2 and 5.5 million years ago). Subsequent diversifications in the genus appear to have occurred in the Pliocene and at the Plio-Pleistocene boundary (Figure 5).

[Figure omitted. See PDF.]

Figure 5. Multilocus trees of Euparkerella.

Bayesian phylogenetic inference of concatenated sequence data (left) and *BEAST Species Tree (right) of Euparkerella. Colors refer to populations indicated on the map (see Table S1 in File S1 for details). Posterior probabilities are indicated near nodes. Asterisks represent posteriors equal or higher than 0.95. Grey bars in Species Tree indicate 95% interval of tMRCA estimated in million years.

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

Delimitation of evolutionary units

The level of genetic structure recovered by the Bayesian clustering methods for both nuclear and combined nuclear and mtDNA datasets of E. brasiliensis, E. cochranae, and Euparkerella sp. was K = 6, as suggested by both the highest log posterior probability of the data (ln Pr (X/K) and the ΔK score (Figure 6). Each of the six sampled localities of Rio de Janeiro was classified into distinct genetic clusters with high individual membership coefficient (0.95 to 1.0). Moreover, there is an overall concordance between these six distinct clusters with terminal clades inferred by phylogenetic analyses.

[Figure omitted. See PDF.]

Figure 6. Structure results for Euparkerella populations.

Assignment proportion of individuals to populations of Euparkerella from Rio de Janeiro State. Structure results for number of groups K=6. Left bar includes both mitochondrial and nuclear genes, right bar includes only nuclear genes.

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

The Bayesian species delimitation using BPP analyses are shown in Figure 7. BPP analyses produced speciation probabilities of 1.0 on all nodes across all three combinations of prior distributions for θ and τ₀ and alternative guide trees (Figure 7). Estimates of genealogical sorting index (gsi) were highly significant (all p-values < 0.02) for each nuclear locus across all six evolutionary units as defined by phylogenetic inferences and Bayesian clustering analysis, and across loci within those units, ranging from 0.55 to 1.0, and from 0.89 to 1.0, respectively (Table 3).

[Figure omitted. See PDF.]

Figure 7. BP&P species delimitation.

Bayesian species delimitation assuming three alternative phylogenetic arrangement based in *BEAST species tree (see Figure 5). Numbers below the nodes represent speciation probabilities values. Colors and numbers refers to populations indicated on the map (Figure 2, see Table S1 in File S1 for details).

https://doi.org/10.1371/journal.pone.0079504.g007

E. cochranae E. cochranae E. brasiliensis E. brasiliensis E. brasiliensis Euparkerella sp.

(3) (4) (5) (6) (7) (8)

gsi P gsi P gsi P gsi P gsi P gsi P

β-fibint7 0,5803 0,0001 0,8327 0,0001 1 0,0001 0,7647 0,0001 1 0,0001 1 0,0001

C-myc2 1 0,0001 1 0,0001 0,8582 0,0001 1 0,0001 1 0,0001 1 0,0001

RAG-1 1 0,0152 0,5486 0,0001 1 0,0001 1 0,0001 1 0,0001 1 0,0002

TYR 1 0,017 1 0,0001 1 0,0002 1 0,0001 1 0,0001 1 0,0001

total (gsi_T) 0,8951 0,0001 0,8453 0,0001 0,9646 0,0001 0,9412 0,0001 1 0,0001 1 0,0001

Table 3. Genealogical sorting index for Euparkerella populations/species.

CSV

Download CSV

Discussion

The first study examining the genetic diversity within the genus Euparkerella uncovered high levels of genetic diversity within the two morphological species E. cochranae and E. brasiliensis, which is in line with previous studies of Terrarana [5,11,12,61]. The gene trees at five independent loci and trees from combined data (concatenated dataset and the species tree) allowed us to uncover at least six deeply diverged and geographically coherent evolutionary units: two within E. cochranae, three within E. brasiliensis, and another corresponding to Euparkerella sp., which may have diversified between the Miocene and the Pleistocene. The delimitation of these six units was supported by all methods herein used, including Bayesian clustering analysis, Bayesian coalescent-based species delimitation, and the Genealogical sorting index, providing thus strong evidence for a current underestimation of diversity in Euparkerella.

Phylogenetic analyses

Topologies of nuclear DNA gene trees were somewhat discordant, revealing less power than mitochondrial DNA (CO-I) in resolving phylogenetic relationships. Discordance between nuclear and mitochondrial markers is expected in analyses of closely related species [62,63]. Our results notably showed however the concordance between mtDNA and most nDNA loci in recovering of several clades within E. cochranae (two clades) and E. brasiliensis (three clades). The mtDNA CO-I tree additionally showed high support for four main groups that correspond to E. robusta and E. tridactyla, which are basal relative to other two well supported sister clades: one composed by E. cochranae 4 and Euparkerella sp. 8, and the other including E. brasiliensis 5, E. brasiliensis 6, E. brasiliensis 7, and E. cochranae 3.

Phylogenetic relationships within Euparkerella are better supported in the tree inferred using the multilocus concatenated dataset when compared either to the mtDNA tree or the multilocus species tree, as observed in several other case studies [30,64–67]. This is likely explained by differences in the sensitivity of distinct methods to gene tree incongruences. While species trees are more sensitive to the resolution of independent gene trees, reconstructions from concatenated datasets enable that lesser resolution given by a particular gene may be compensated by the higher resolution of another gene [45,67,68]. The higher resolution of the CO-I data would compensate the lesser resolution of the nuclear loci in the tree inferred from concatenated data, and incongruences between the several loci would decrease the resolution of the species tree. This appears to be exactly the case of clade E. cochranae 3 - E. brasiliensis, which shows distinct degrees of resolution depending on the tree reconstruction method (well resolved for the CO-I and concatenated datasets, and a polytomy for the species tree). Thus, we suggest that the conflict between the tree inferred from concatenated data and the species tree is explained by gene tree incongruences affecting the latter and not necessarily by a wrong tree inference from the concatenated dataset [69]. The degree of population and gene sampling undertaken in this initial study of genetic diversification within Euparkerella limits our ability to discuss this matter any further. Future work sampling more individuals and localities of all species of Euparkerella, and covering more of the genome diversity perhaps will result in a more complete understanding of phylogenetic relationships in this genus.

Species delimitation: diversification within the genus Euparkerella

The combined nuclear and cytoplasmic DNA data uncovered several genetic units in Euparkerella, especially within the previously recognized species E. cochranae and E. brasiliensis. There is general concordance between gene trees and trees from combined data in recovering two clades within E. cochranae (E. cochranae 3 and E. cochranae 4), three clades within E. brasiliensis (E. brasiliensis 5, E. brasiliensis 6, and E. brasiliensis 7), and Euparkerella sp., all corresponding to specific localities or regions. Population structure inferred from Bayesian analysis of multilocus nuclear data, and from combined cytoplasmic and nuclear data, suggests six genetically differentiated units that correspond to the six clades recovered from phylogenetic analyses. Each geographic region corresponds to a distinct genetic group with probabilities between 0.95 and 1.0. This overall congruence between results from both phylogenetic and Bayesian clustering methods give additional confidence for the delimitation of six distinct evolutionary units within the three previously defined Euparkerella species. So, despite the limitations of our data (reduced sample size and meaningful species-level genetic diversity), this result is in line with the observed in other studies that used Bayesian clustering analysis to identify and/or delimit evolutionary units, even within widespread species complexes [56,58,70]

The Bayesian species delimitation method (BPP) consistently supported the same six distinct evolutionary units as inferred by both phylogenetic and Bayesian clustering analysis, even accounting for variation on population sizes and divergence times. The BPP analysis may have the drawback of relying on an accurate guide tree [55] which in our case (*BEAST species tree) had some ambiguous relationships estimates. Previous studies have indeed shown that BPP can produce misleading results due to incorrectly prior-specified trees [50,55,71]. However, Euparkerella species delimitation using BPP consistently recovered posterior probabilities of 1.0 for every speciation event when we used different guide trees representing all possible branching patterns in regions of the species tree with low support. Moreover, using an independent method (GSI), these six genetic evolutionary units are characterized by high levels of genealogical exclusivity both for individual nuclear gene trees (gsi) and in combined multilocus analysis (gsi_T) (Table 3), which is consistent with a history of independent divergence for these units [57,58,72].

Whether our genetic units conform or not to a particular species concept is less relevant than the processes of old evolutionary diversification uncovered for E. brasiliensis and E. cochranae. The observed pairwise mtDNA genetic divergences among the three units (clades) of E. brasiliensis (6,3% - 11,2%) and between the two units (clades) of E. cochranae (13,6%) are generally similar or larger than genetic divergences observed among other amphibian species [4,25,73,74]. Divergence times inferred from the Euparkerella species tree suggest that diversification of those genetic units occurred between the Miocene and the Plio-Pleistocene boundary (1.4 - 5.1 mybp), which is also concordant to times of divergence inferred for other species of amphibians [27,34], including species of Terrarana [12].

Geographical patterns of genetic diversity

High levels of genetic substructure were observed within Euparkerella, with all know genetic units occurring in single localities or geographically restricted locality groups, without any records of sympatry. The six genetic units uncovered within Euparkerella from Rio de Janeiro State correspond to populations that are between 15 and 100 km away from the closest congeneric populations. Diversification in relatively small geographic scales seems relatively common in small-bodied, and particularly in miniaturized vertebrates [25,27,30,75,76]. Organismal body size has several functional, ecological, and evolutionary implications, in particular reducing physiological tolerance (e.g., higher sensitivity to dehydration) and dispersal in the case of extremely small-bodied species [27,77]. The lower levels of gene flow between populations induced by those morphophysiological traits may lead to range fragmentation at relatively small geographic scales and posterior genetic diversification, even without apparent barriers to gene flow [27]. We hypothesize the same general processes in small-bodied or miniaturized species determined the geographical structure and divergence observed in E. brasiliensis and E. cochranae, perhaps combined with habitat heterogeneity. Within the range of Euparkerella in Rio de Janeiro, the landscape is formed by mountains that promote altitudinal climatic variation (from sea level to 1000 meters) accompanied by considerable physiognomic diversity (Restinga - a coastal vegetation with marine influence; Lowland rainforest, 5–50m a.s.l..; Submontane rainforest, 50–500m a.s.l.; and Montane rainforest, 500–1500m a.s.l.) [78]. The three genetic units recovered within E. brasiliensis occur in three distinct physiognomies: E. brasiliensis 7 occurs in the Restinga, whereas E. brasiliensis 5 and 6 occur in Submontane and Montane rainforest. Our results are indicative of some association between small scale habitat heterogeneity and the diversification of Euparkerella. This study cannot rule out factors that are commonly mentioned as promoters of geographic isolation, such as geographic barriers, climate change, and other sources of habitat disturbance [79–82], to have been involved in the diversification of Euparkerella. Alternative hypotheses explaining diversification in Euparkerella will only be properly examined when studies combining higher geographic sampling with analyses of physiological, behavioral, ecological, and environmental variables are undertaken in the future.

Taxonomical and conservation implications

The present study was the first to examine genetic diversity within the genus Euparkerella. Genetic diversity was very high and six distinct genetic units have likely independently diversified since the late Pliocene. Our results imply in some cases that morphological traits previously used to distinguish species are not truly diagnostic. For example, E. cochranae 3 is morphologically indistinguishable from E. cochranae 4 (the two are currently recognized as the same species), whereas both are not sister clades in our phylogenetic analyses. Discordance between molecular and morphological data in diagnosing small-bodied or miniaturized species is relatively common and may be associated to the reduction of visible characters [75]. One way to overcome this limitation will be to undertake detailed morphological, anatomical, and osteological analyses using high resolution techniques, such as geometrical morphometry, radiometric imaging and high-contrast 3D imaging throughout X-ray micro and nanotomography [25,30,83]. An integrative analysis of morphology, physiology, behavior, and ecology will likely be instrumental in the future to properly undertake a much needed taxonomic revision of the genus Euparkerella.

Taken together, the high levels of genetic diversity uncovered within Euparkerella at small geographic scales suggest that the spatial scale of sampling may be essential to find cryptic diversity in reduced size amphibians. Given the intrinsic morphoecological characteristics of Euparkerella combined with habitat heterogeneity it is likely that genetic units in this genus may increase with higher sampling, especially across the ranges of the species E. tridactyla and E. robusta. We predict that such levels of cryptic diversity, whether corresponding to incipient or fully-fledged speciation, may be excellent models to examine the origins and maintenance of microendemism in the context of spatial heterogeneity and/or human induced fragmentation of the highly threatened Brazilian Atlantic forest hotspot.

Supporting Information

[Figure omitted. See PDF.]

Supporting Information Tables S1, S2, S3.

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

(DOCX)

Protocol S1.

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

(DOCX)

Acknowledgments

We are grateful to numerous colleagues that contributed to this work, especially to Bruno Lamy and Andreza Siqueira for fieldwork assistance and to Mariana lyra and Tiago S. Neves for laboratory and/ or data analysis assistance. Instituto Chico Mendes de Conservação da Biodiversidade (ICMBio) and Instituto Estadual do Ambiente (INEA) for collecting permits (18887 and 040/2010, respectively). To curators of zoological collections José P. Pombal Jr. (MNRJ) and Leonora Pires Costa (UFES), and colleague Felipe B. S. Telles for tissue samples provided.

Author Contributions

Conceived and designed the experiments: LAF JA CFBH CFDR FS. Performed the experiments: LAF TOB. Analyzed the data: LAF TOB FS. Contributed reagents/materials/analysis tools: LAF JA CFBH CFDR FS. Wrote the manuscript: LAF JA CFBH TOB CFDR FS.

Citation: Fusinatto LA, Alexandrino J, Haddad CFB, Brunes TO, Rocha CFD, Sequeira F (2013) Cryptic Genetic Diversity Is Paramount in Small-Bodied Amphibians of the Genus Euparkerella (Anura: Craugastoridae) Endemic to the Brazilian Atlantic Forest. PLoS ONE 8(11): e79504. https://doi.org/10.1371/journal.pone.0079504

References

1. Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J (2000) Biodiversity hotspots for conservation priorities. Nature 403: 853–858. doi:10.1038/35002501. PubMed: 10706275.

2. Wilson EO (2003) The future of life. New York: Vintage Books.

3. Faivovich J, Haddad CFB, Garcia PCA, Frost DR, Campbell JA et al. (2005) Systematic review of the frog family Hylidae, with special reference to Hylinae: phylogenetic analysis and taxonomic revision. Bull Am Museum Nat Hist 294: 1–240. doi:10.1206/0003-0090(2005)294[0001:SROTFF]2.0.CO;2.

4. Fouquet A, Gilles A, Vences M, Marty C, Blanc M et al. (2007) Underestimation of species richness in neotropical frogs revealed by mtDNA analyses. PLOS ONE 2: e1109. doi:10.1371/journal.pone.0001109. PubMed: 17971872.

5. Hedges SB, Duellman WE, Heinicke MP (2008) World N, Molecular T. New World Directdev Frogs (Anura: Terrarana): molecular phylogeny, classification, biogeography, and conservation. Zootaxa 1737: 1–182.

6. De La Riva I (2007) Bolivian frogs of the genus Phrynops, with the description of twelve new species (Anura: Brachycephalidae). Herpetol Monogr 21: 241–277. doi:10.1655/07-011.1.

7. Vieites DR, Wollenberg KC, Andreone F, Köhler J, Glaw F et al. (2009) Vast underestimation of Madagascar’s biodiversity evidenced by an integrative amphibian inventory. Proc Natl Acad Sci U S A 106: 8267–8272. doi:10.1073/pnas.0810821106. PubMed: 19416818.

8. Hanken J, Wake D (1994) Five new species of minute salamanders, genus Thorius (Caudadata: Plethodontidae), from Northern Oaxaca, Mexico. Copeia: 573–590.

9. Campbell JA, Savage JM (2000) Taxonomic reconsideration of the Middle American Frogs of the Eleutherodactylus rugulosus group (Anura: Leptodactylidae): a reconnaissance of subtle nuances among frogs. Herpetol Monogr 14: 186–292. doi:10.2307/1467048.

10. Köhler J, Vieites DR, Bonett RM, García FH, Glaw F et al. (2005) New amphibians and global conservation: a boost in species discoveries in a highly endangered vertebrate group. BioScience 55: 693–696. doi:10.1641/0006-3568(2005)055[0693:NAAGCA]2.0.CO;2.

11. Heinicke MP, Duellman WE, Trueb L, Means DB, Macculloch RD et al. (2009) A new frog family (Anura: Terrarana) from South America and an expanded direct-developing clade revealed by molecular phylogeny. Zootaxa 2211: 1–35.

12. Rodríguez A, Vences M, Nevado B, Machordom A, Verheyen E (2010) Biogeographic origin and radiation of Cuban Eleutherodactylus frogs of the auriculatus species group, inferred from mitochondrial and nuclear gene sequences. Mol Phylogenet Evol 54: 179–186. doi:10.1016/j.ympev.2009.08.023. PubMed: 19706330.

13. Cherry LM, Case SM, Wilson AC (1977) Relative rates of morphological evolution in primates, carnivores and frogs. Am Zool 17: 910–910.

14. Stuart BL (2008) The phylogenetic problem of Huia (Amphibia: Ranidae). Mol Phylogenet Evol 46: 49–60. doi:10.1016/j.ympev.2007.09.016. PubMed: 18042407.

15. Wiens JJ (2008) Systematics and Herpetology in the Age of Genomics. BioScience 58: 297–307. doi:10.1641/B580405.

16. Pyron RA, Wiens JJ (2011) A large-scale phylogeny of Amphibia including over 2800 species, and a revised classification of extant frogs, salamanders, and caecilians. Mol Phylogenet Evol 61: 543–583. doi:10.1016/j.ympev.2011.06.012. PubMed: 21723399.

17. Frost DR, Grant T, Faivovich J, Bain RH, Haas A et al. (2006) The amphibian tree of life. Bull Am Museum Nat Hist: 1–370. doi:10.1206/0003-0090(2006)297[0001:TATOL]2.0.CO;2.

18. Heinicke MP, Duellman WE, Hedges SB (2007) Major Caribbean and Central American frog faunas originated by ancient oceanic dispersal. Proc Natl Acad Sci U S A 104: 10092–10097. doi:10.1073/pnas.0611051104. PubMed: 17548823.

19. Canedo C, Haddad CFB (2012) Phylogenetic relationships within anuran clade Terrarana, with emphasis on the placement of Brazilian Atlantic rainforest frogs genus Ischnocnema (Anura: Brachycephalidae). Mol Phylogenet Evol 65: 610–620. doi:10.1016/j.ympev.2012.07.016. PubMed: 22842090.

20. Lynch JD (1976) The species groups of the South American frogs of the genus Eleutherodactylus (Leptodactylidae). Kansas: Occas Pap Museum Nat Hist Univ. pp. 1–24.

21. Savage JM (1987) Systematics and distribution of the Mexican and Central American rainfrogs of the Eleutherodactylus gollmeri group (Amphibia: Leptodactylidae). Fieldiana Zool New Ser 33: 1–57.

22. Hedges SB (1989) Evolution and biogeography of West Indian frogs of the genus Eleutherodactylus: slow-evolving loci and the major groups. In: CA Woods. Biogeography of the West Indies: past, present and future. Gainesville: Sandhill Crane Press. pp. 305–370.

23. Izecksohn E (1988) Algumas considerações sobre o gênero Euparkerella, com descrição de três novas espécies (Amphibia, Anura, Leptodactylidae). Rev Bras Biol 48: 59–74.

24. Frost DR (2013) Amphibian Species of the World: an online reference. version 5.6 (9 January 2013); Am Museum Nat Hist. Available: http://research.amnh.org/vz/herpetology/amphibia/.

25. Rittmeyer EN, Allison A, Gründler MC, Thompson DK, Austin CC (2012) Ecological guild evolution and the discovery of the world’s smallest vertebrate. PLOS ONE 7: e29797. doi:10.1371/journal.pone.0029797. PubMed: 22253785.

26. Zimkus BM, Lawson L, Loader SP, Hanken J (2012) Terrestrialization, miniaturization and rates of diversification in African puddle frogs (Anura: Phrynobatrachidae). PLOS ONE 7: e35118. doi:10.1371/journal.pone.0035118. PubMed: 22509392.

27. Wollenberg KC, Vieites DR, Glaw F, Vences M (2011) Speciation in little: the role of range and body size in the diversification of Malagasy mantellid frogs. BMC Evol Biol 11: 1–15. doi:10.1186/1471-2148-11-217. PubMed: 21194491.

28. Hanken J, Wake DB (2001) A seventh species of minute salamander (Thorius: Plethodontidae) from the Sierra de Juárez, Oaxaca, México. Herpetologica 57: 515–523.

29. Harvey MB, Almendáriz ANA, JB M, R DB (2013) A new species of Noblella (Anura: Craugastoridae) from the Amazonian Slopes of the Ecuadorian Andes with comments on Noblella lochites (Lynch). Zootaxa 3635: 1–14. doi:10.11646/zootaxa.3635.1.1.

30. Clemente-Carvalho RBG, Klaczko J, Perez SI, Alves ACR, Haddad CFB et al. (2011) Molecular phylogenetic relationships and phenotypic diversity in miniaturized toadlets, genus Brachycephalus (Amphibia: Anura: Brachycephalidae). Mol Phylogenet Evol 61: 79–89. doi:10.1016/j.ympev.2011.05.017. PubMed: 21693192.

31. Carstens BC, Pelletier Ta, Reid NM, Satler JD (2013) How to fail at species delimitation. Mol Ecol, 22: 4369–83. doi:10.1111/mec.12413. PubMed: 23855767.

32. Sequeira F, Ferrand N, Harris DJ (2006) Assessing the phylogenetic signal of the nuclear β-Fibrinogen intron 7 in salamandrids (Amphibia: Salamandridae). Amphibia-Reptilia 27: 409–418. doi:10.1163/156853806778190114.

33. Crawford AJ (2003) Relative rates of nucleotide substitution in frogs. J Mol Evol 57: 636–641. doi:10.1007/s00239-003-2513-7. PubMed: 14745532.

34. Brunes TO, Sequeira F, Haddad CFB, Alexandrino J (2010) Gene and species trees of a Neotropical group of treefrogs: genetic diversification in the Brazilian Atlantic Forest and the origin of a polyploid species. Mol Phylogenet Evol 57: 1120–1133. doi:10.1016/j.ympev.2010.08.026. PubMed: 20813192.

35. Chiari Y, Vences M, Vieites DR, Rabemananjara F, Bora P et al. (2004) New evidence for parallel evolution of colour patterns in Malagasy poison frogs (Mantella). Mol Ecol 13: 3763–3774. doi:10.1111/j.1365-294X.2004.02367.x. PubMed: 15548289.

36. Bossuyt F, Milinkovitch MC (2000) Convergent adaptive radiations in Madagascan and Asian ranid frogs reveal covariation between larval and adult traits. Proc Natl Acad Sci U S A 97: 6585–6590. doi:10.1073/pnas.97.12.6585. PubMed: 10841558.

37. Hall TA (1999) BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser: 95–98. PubMed: 10780396.

38. Stephens M, Donnelly P (2003) A comparison of bayesian methods for haplotype reconstruction from population genotype data. Am J Hum Genet 73: 1162–1169. doi:10.1086/379378. PubMed: 14574645.

39. Librado P, Rozas J (2009) DnaSP v5 : a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451–1452. doi:10.1093/bioinformatics/btp187. PubMed: 19346325.

40. Tamura K, Peterson D, Peterson N, Stecher G, Nei M et al. (2011) MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol 28: 2731–2739. doi:10.1093/molbev/msr121. PubMed: 21546353.

41. Lanfear R, Calcott B, Ho SYW, Guindon S (2012) PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol 29: 1695–1701. doi:10.1093/molbev/mss020. PubMed: 22319168.

42. Akaike H (1973) Information theory as an extension of the maximum likelihood principle. In: BN PetrovF. Csaki. Second International Symposium of Information Theory. . Budapest: Akademiai Kiado. pp. 267–281.

43. Ronquist F, Teslenko M, Van der Mark P, Ayres DL, Darling A et al. (2012) MrBayes 3.2: efficient bayesian phylogenetic inference and model choice across a large model space. Syst Biol 61: 539–542. doi:10.1093/sysbio/sys029. PubMed: 22357727.

44. Rambaut A, Drummond AJ (2009) Tracer v1.5. Available http://tree.bio.ed.ac.uk/software/tracer/.

45. Heled J, Drummond AJ (2010) Bayesian inference of species trees from multilocus data. Mol Biol Evol 27: 570–580. doi:10.1093/molbev/msp274. PubMed: 19906793.

46. Drummond AJ, Suchard MA, Xie D, Rambaut A (2012) Bayesian phylogenetics with BEAUti and the BEAST. p. 1.7. Mol Biol Evol 29: 1969–1973 doi:10.1093/molbev/mss075.

47. Pulquério MJF, Nichols RA (2006) Dates from the molecular clock: how wrong can we be ? Trends Ecol Evol 22: 180–184. doi:10.1016/j.tree.2006.11.013. PubMed: 17157408.

48. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155: 945–959. PubMed: 10835412.

49. 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. doi:10.1111/j.1365-294X.2005.02553.x. PubMed: 15969739.

50. Carstens BC, Satler JD (2013) The carnivorous plant described as Sarracenia alata contains two cryptic species. Biol J Linn Soc 109: 737–746. doi:10.1111/bij.12093.

51. Earl DA, VonHoldt BM (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour 4: 359–361. doi:10.1007/s12686-011-9548-7.

52. Jakobsson M, Rosenberg NA (2007) CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23: 1801–1806. doi:10.1093/bioinformatics/btm233. PubMed: 17485429.

53. Rosenberg NA (2004) DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes 4: 137–138. doi:10.1046/j.1471-8286.2003.00566.x.

54. Yang Z, Rannala B (2010) Bayesian species delimitation using multilocus sequence data. Proc Natl Acad Sci U S A 107: 9264–9269. doi:10.1073/pnas.0913022107. PubMed: 20439743.

55. Leaché AD, Fujita MK (2010) Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus). Proc R Soc Lond B Biol Sci 277: 3071–3077. doi:10.1098/rspb.2010.0662. PubMed: 20519219.

56. Welton LJ, Siler CD, Oaks JR; Diesmos a C, Brown RM (2013) Multilocus phylogeny and Bayesian estimates of species boundaries reveal hidden evolutionary relationships and cryptic diversity in Southeast Asian monitor lizards. Mol Ecol, 22: 3495–510. doi:10.1111/mec.12324. PubMed: 23710945. PubMed: 23710945.

57. Cummings MP, Neel MC, Shaw KL (2008) A genealogical approach to quantifying lineage divergence. Evolution (N Y) 62: 2411–2422. doi:10.1111/j.1558-5646.2008.00442.x.

58. Weisrock DW, Rasoloarison RM, Fiorentino I, Ralison JM, Goodman SM et al. (2010) Delimiting species without nuclear monophyly in Madagascar’s Mouse Lemurs. PLOS ONE 5: e9883. doi:10.1371/journal.pone.0009883. PubMed: 20360988.

59. Silvestro D, Michalak I (2012) raxmlGUI: a graphical front-end for RAxML. Org Divers Evol 12: 335–337. doi:10.1007/s13127-011-0056-0.

60. Stamatakis A (2006) RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics 22: 2688–2690. doi:10.1093/bioinformatics/btl446. PubMed: 16928733.

61. Kieswetter CM, Schneider CJ (2013) Phylogeography in the northern Andes: complex history and cryptic diversity in a cloud forest frog, Pristimantis w-nigrum (Craugastoridae). Mol Phylogenet Evol. doi:10.1016/j.ympev.2013.08.007.

62. Edwards SV, Beerli P (2000) Perspective: gene divergence, population divergence, and the variance in coalescence time in phylogeographic studies. Evolution (N Y) 54: 1839–1854.

63. Degnan JH, Rosenberg NA (2009) Gene tree discordance, phylogenetic inference and the multispecies coalescent. Trends Ecol Evol 24: 332–340. doi:10.1016/j.tree.2009.01.009. PubMed: 19307040.

64. Brumfield RT, Liu L, Lum DE, Edwards SV (2008) Comparison of species tree methods for reconstructing the phylogeny of Bearded Manakins (Aves: Pipridae, Manacus) from multilocus sequence data. Syst Biol 57: 719–731. doi:10.1080/10635150802422290. PubMed: 18853359.

65. Wiens JJ, Kuczynski CA, Stephens PR (2010) Discordant mitochondrial and nuclear gene phylogenies in emydid turtles: implications for speciation and conservation. Biol J Linn Soc 99: 445–461. doi:10.1111/j.1095-8312.2009.01342.x.

66. Leaché AD (2010) Species trees for spiny lizards (Genus Sceloporus): Identifying points of concordance and conflict between nuclear and mitochondrial data. Mol Phylogenet Evol 54: 162–171. doi:10.1016/j.ympev.2009.09.006. PubMed: 19751840.

67. Townsend TM, Mulcahy DG, Noonan BP, Sites JW Jr., Kuczynski CA et al. (2011) Phylogeny of iguanian lizards inferred from 29 nuclear loci, and a comparison of concatenated and species-tree approaches for an ancient, rapid radiation. Mol Phylogenet Evol 61: 363–380. doi:10.1016/j.ympev.2011.07.008. PubMed: 21787873.

68. Pasachnik SA, Echternacht AC, Fitzpatrick BM (2010) Gene trees, species and species trees in the Ctenosaura palearis clade. Conserv Genet 11: 1767–1781. doi:10.1007/s10592-010-0070-3.

69. Belfiore NM, Liu L, Moritz C (2008) Multilocus Phylogenetics of a Rapid Radiation in the Genus Thomomys (Rodentia: Geomyidae). Syst Biol 57: 294–310. doi:10.1080/10635150802044011. PubMed: 18432550.

70. Thomé MTC, Zamudio KR, Haddad CFB, Alexandrino J (2012) Delimiting genetic units in Neotropical toads under incomplete lineage sorting and hybridization. BMC Evol Biol 12: 1–13. doi:10.1186/1471-2148-12-242. PubMed: 22214193.

71. Miralles A, Vences M (2013) New metrics for comparison of taxonomies reveal striking discrepancies among species delimitation methods in Madascincus lizards. PLOS ONE 8: e68242. doi:10.1371/journal.pone.0068242. PubMed: 23874561.

72. Sakalidis ML, Hardy GESJ, Burgess TI (2011) Use of the Genealogical Sorting Index (GSI) to delineate species boundaries in the Neofusicoccum parvum-Neofusicoccum ribis species complex. Mol Phylogenet Evol 60: 333–344. doi:10.1016/j.ympev.2011.04.026. PubMed: 21575732.

73. Smith MA, Poyarkov Na, Hebert PDN (2007) CO1 DNA barcoding amphibians: take the chance, meet the challenge. Mol Ecol Notes 8: 235–246. doi:10.1111/j.1471-8286.2007.01964.x.

74. Alonso R, Crawford AJ, Bermingham E (2012) Molecular phylogeny of an endemic radiation of Cuban toads (Bufonidae: Peltophryne) based on mitochondrial and nuclear genes. J Biogeogr 39: 434–451. doi:10.1111/j.1365-2699.2011.02594.x.

75. Glaw F, Köhler J, Townsend TM, Vences M (2012) Rivaling the world’s smallest reptiles: discovery of miniaturized and microendemic new species of leaf chameleons (Brookesia) from northern Madagascar. PLOS ONE 7: e31314. doi:10.1371/journal.pone.0031314. PubMed: 22348069.

76. Pabijan M, Wollenberg KC, Vences M (2012) Small body size increases the regional differentiation of populations of tropical mantellid frogs (Anura: Mantellidae). J Evol Biol 25: 2310–2324. doi:10.1111/j.1420-9101.2012.02613.x. PubMed: 22998688.

77. Hanken J, Wake DB (1993) Miniaturization of body size: organismal consequences and evolutionary significance. Annu Rev Ecol Syst 24: 501–519. doi:10.1146/annurev.es.24.110193.002441.

78. Veloso HP, Rangel Filho ALR, Lima JCA (1991) Classificaçao da vegetação brasileira, adaptada a um sistema universal. Rio de Janeiro: IBGE, Departamento de Recursos Naturais e Estudos Ambientais.

79. Hewitt G (2000) The genetic legacy of the Quaternary ice ages. Nature 405: 907–913. doi:10.1038/35016000. PubMed: 10879524.

80. Moritz C, Patton JL, Schneider CJ, Smith TB (2000) Diversification of rainforest faunas: an integrated molecular approach. Annu Rev Ecol Evol Syst 31: 533–563. Available: http://arjournals.annualreviews.org/doi/abs/10.1146/annurev.ecolsys.31.1.533. doi:10.1146/annurev.ecolsys.31.1.533.

81. Ghalambor CK, Huey RB, Martin PR, Tewksbury JJ, Wang G (2006) Are mountain passes higher in the tropics? Janzen’s hypothesis revisited. Integr Comp Biol 46: 5–17. doi:10.1093/icb/icj003. PubMed: 21672718.

82. Carnaval AC, Moritz C (2008) Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic forest. J Biogeogr 35: 1187–1201. doi:10.1111/j.1365-2699.2007.01870.x.

83. Metscher BD (2009) MicroCT for comparative morphology: simple staining methods allow high-constrast 3D imaging of diverse non-mineralized animal tissues. BMC Physiol 9: 11. doi:10.1186/1472-6793-9-11. PubMed: 19545439Retrieved onpublished at whilst December year 1111 from. . doi:10.1186/1472-6793-9-11.