Speciation in a biodiversity hotspot: Phylogenetic relationships, species delimitation, and divergence times of Patagonian ground frogs from the Eupsophus roseus group (Alsodidae).

The alsodid ground frogs of the Eupsophus genus are divided into two groups, the roseus (2n = 30) and vertebralis (2n = 28), which are distributed throughout the temperate Nothofagus forests of South America. Currently, the roseus group is composed by four species, while the vertebralis group consists of two. Phylogenetic relationships and species delimitation within each group are controversial. In fact, previous analyses considered that the roseus group was composed of between four to nine species. In this work, we evaluated phylogenetic relationships, diversification times, and species delimitation within the roseus group using a multi-locus dataset. For this purpose, mitochondrial (D-loop, Cyt b, and COI) and nuclear (POMC and CRYBA1) partial sequences from 164 individuals were amplified, representing all species. Maximum Likelihood (ML) and Bayesian approaches were used to reconstruct phylogenetic relationships. Species tree was estimated using BEAST and singular value decomposition scores for species quartets (SVDquartets). Species limits were evaluated with six coalescent approaches. Diversification times were estimated using mitochondrial and nuclear rates with LogNormal relaxed clock in BEAST. Nine well-supported monophyletic lineages were recovered in Bayesian, ML, and SVDquartets, including eight named species and a lineage composed by specimens from the Villarrica population (Bootstrap:>70, PP:> 0.99). Single-locus species delimitation analyses overestimated the species number in E. migueli, E. calcaratus, and E. roseus lineages, while multi-locus analyses recovered as species the nine lineages observed in phylogenetic analyses (Ctax = 0.69). It is hypothesized that Eupsophus diversification occurred during Mid-Pleistocene (0.42-0.14 Mya), with most species having originated after the Last Southern Patagonian Glaciation (0.18 Mya). Our results revitalize the hypothesis that the E. roseus group is composed of eight species and support the Villarrica lineage as a new putative species.

Introduction

From an operational point of view, the notion of biodiversity encompasses several different levels of biological organization, from the genetic make up of the species to ecosystems and landscapes, in which the species is the most significant unit. Species are used for comparisons in almost all biological fields including ecology, evolution, and conservation [1-3], and there is no doubt the central unit for systematics is also the species [4]. Furthermore, biodiversity hotspots are selected on the basis of the species they possess, conservation schemes are assessed on how many species are preserved, and conservation legislation and politics are focused on species preservation [5,6].

Despite the importance of the species concepts debate [7,8], and since the species as taxonomic hierarchy is also considered a fundamental topic in biology [9], it is broadly accepted that species are best conceptualized as dynamic entities connected by "grey zones" where their delimitation will remain inherently ambiguous [4,10]. Under this perspective, species delimitation, i.e. the act of identifying biological diversity at species-level [11], is particularly challenging in actively radiating groups composed of recently diverged lineages. The difficulty lies in the fact that recently separated species are less likely to possess all or even many of the diagnosable characters such as phenetic distinctiveness, intrinsic reproductive incompatibility, ecological uniqueness, or reciprocal monophyly, that constitute operational criteria for their delimitation [4,12]. This becomes more complex when hybridization and introgression among related species are considered common and major contributors to speciation and diversification [13]. Genealogical discordance obtained with different markers is a result of these processes, but also of incomplete lineage sorting, selection, or demographic disparities [14]. Thus, hypotheses of the boundaries of recently diverged species may remain unclear due to incomplete lineage sorting, introgression, complex of cryptic species that cannot be distinguished by morphology alone, sampling deficiencies, or different taxonomic practices [2,4].

Ever since genetic data became easier and less expensive to gather, the field of species delimitation has experienced an explosion in the number and variety of methodological approaches [3,11,15-17]. These new approaches proceed by evaluating models of lineage composition under a phylogenetic framework that implements a coalescent model to delimit the species [11,18]. In this regard, these approaches estimate the phylogeny while allowing for the action of population-level processes, such as genetic drift in combination with migration, expansion, population divergence, or combinations of these processes [19-21]. Thus, the species delimitation models can involve population size parameters (i.e. θs for the extant species and common ancestors), parameters for the divergence times (τ), and coalescent models specifying the distribution of gene trees at different loci [22-26].

Some methodological approaches to species delimitation use single-locus sequence information itself as the primary information source for establishing group membership and defining species boundaries [27-29]. Other methods are designed to analyze multi-locus data sets and require a priori assignment of individuals to species categories [21,30,31]. The performance of species delimitation methods are quantified by the number of different species recognized in each case, and the congruence with data at hand such as life history, geographical distribution, morphology, and behavior [15,32]. Although there is difficulty to integrate genetic and non-genetic data to increase the efficacy of species detection [33], there are available methods to measure the congruence and resolving power among species delimitation approaches [34].

The history of the Patagonian landscape offers exceptional opportunities to investigate diversification and promote conservation strategies by studying the past, present, and future of evolutionary processes using amphibians as a model of study. In this region, the amphibian fauna of Chile is not particularly diverse (60 species; [35]) but includes 10 endemic genera, some of them having one, a few species (e. g. Calyptocephalella, Chaltenobatrachus, Hylorina, Insuetophrynus, Rhinoderma), or as many as 18 (Alsodes). Among these amphibians are the frogs of the genus Eupsophus Fitzinger 1843. This taxon currently includes six species distributed almost throughout the temperate Nothofagus forest of South America [35]. Nevertheless, Eupsophus have puzzled frog systematics for decades [36-39], and a clear consensus has not yet been reached regarding the number of species that make up this genus [40-42]. In fact, the genus Eupsophus was classically divided into two groups with following species [36,43]: 1) roseus group, composed of E. altor, E. roseus, E. calcaratus, E. contulmoensis, E. insularis, E. septentrionalis, E. migueli and E. nahuelbutensis, all of them with 30 chromosomes, with individuals of 34–42 mm body size (snout-vent distance) [44]; and 2) the vertebralis group, composed of E. vertebralis and E. emiliopugini, both species with 28 chromosomes and individuals with a body size of 50–59 mm (snout-vent distance) [44]. Nevertheless, recently molecular analyses within the roseus group synonymized E. altor with E. migueli as well as E. contulmoensis, E. septentrionalis, and E. nahuelbutensis with E. roseus [37]. Therefore, the roseus group is currently composed by four species: E. migueli, E. insularis, E. roseus and E. calcaratus [35].

In this study, we present phylogenetic and species delimitation of the roseus group, using 164 new samples from all species covering most of their distribution range. We used three mitochondrial and two nuclear markers, three of them are different to those used by Blotto et al. [36] and Correa et al. [37] [Control Region (D-loop), Propiomelanocortin (POMC), and β Crystallin A1 (CRYBA1)]. These molecular datasets were used to carry out phylogenetic reconstructions and an extensive number of single- and multi-locus species delimitation methods. Species trees and diversification times were estimated to support phylogenetic and species boundaries inferences. New samples, different markers, and multiple bioinformatic techniques allowed us to test, in an independent way, phylogenetic and species delimitation hypothesis of the roseus group.

Materials and methods

Ethics statement

This study was carried out under supervision and approval of the Bioethics and Biosecurity Committee of the Universidad Austral de Chile (UACh, Resolutions No. 236/2015 and 61/15), and the Servicio Agrícola y Ganadero (SAG, Resolution No. 9244/2015). After capture, animals were kept in the dark in fabric bags for a maximum of two hours. Euthanasia was carried out in the field by an overdose of benzocaine 50 mg/mL in a humid chamber. The Corporación Nacional Forestal, Ministerio de Agricultura, Gobierno de Chile allows to collect buccal swabs samples of Eupsophus species from wild protected areas (CONAF, Permit No. 11/2016.-CPP/ MDM/jcr/ 29.02.2016).

Sample collection

A total of 164 samples of Eupsophus from 45 localities in Chile were analysed (Fig 1, S1 Table). Each sampling site was geo-referenced with a GPS Garmin GPSmap 76CSx. Two E. emiliopugini individuals, three E. vertebralis, and one Alsodes norae were used as outgroups (S1 Table, gray cells). Although mostly samples were obtained from buccal swabs according to Broquet et al. [45], some animals were euthanized. Liver tissue was extracted, conserved in 100% ethanol, and stored at -20°C. The specimens were deposited in herpetological collection of the Institute of Marine and Limnological Sciences, Universidad Austral de Chile (ICMLH). The voucher and isolate numbers were included in the sequences information.

Fig 1

Map depicting 45 localities of Eupsophus samples from Chile (listed in S1 Table).

E. roseus: localities 1–16 (red), E insularis: locality 17 (purple), E. migueli: localities 18–20 (blue), E. calcaratus: localities 21–43 (yellow). Localities of outgroup were: E. emiliopugini: 44 and 45 (white), E. vertebralis: 12, 19, 22, Alsodes norae: 19.

DNA extraction, amplification, and sequence alignment

Whole genomic DNA was extracted using Chelex following Walsh et al. [46]. We amplified via the polymerase chain reaction (PCR) three mitochondrial regions: a segment of D-loop [47], Cytochrome b (Cyt b; [48, 49]), and Cytochrome oxidase subunit I (COI; [50]), and two nuclear regions: POMC [51], and CRYBA1 [52]. We mixed reaction cocktails for PCR using 100 ng DNA, 10 μmol of each oligonucleotide primer, 2X of Platinum Taq DNA Polymerase master mix (Invitrogen, Cat. No. 10966), and nuclease-free water to final volume of 25 μL. We verified successful PCR qualitatively by viewing bands of appropriate size following electrophoresis on 1.0% agarose gels. PCR products were sequenced in Macrogen Inc. (Seoul, Korea). Electropherograms were visualized and aligned with Geneious v.9.1.3 (GeneMatters Corp.) using the iterative method of global pairwise alignment (Muscle and ClustalW) implemented in the same software [53,54]. An inspection of aligned sequences by eye as well as manual corrections was also carried out. We expanded our dataset with sequences of E. calcaratus reported in Nuñez et al. [55], and mitochondrial sequences reported by Suárez-Villota et al. [49]. All newly generated sequences from Eupsophus and Alsodes were submitted to GenBank (MK180849-MK181499).

Phylogenetic analyses

Phylogenetic trees were constructed with concatenated dataset using Maximum Likelihood (ML) and Bayesian inference (BI). Evolutionary models and partitioning strategies were evaluated with Partitionfinder v2.1.1 [56] and the best partition was identified using the Bayesian information criterion [57]. ML trees were inferred using GARLI v2.0 [58] with branch support estimated by nonparametric bootstrap (1000 replicates) [59]. Bayesian analyses were performed using MrBayes v3.2 [60]. Each Markov chain Monte Carlo (MCMC) was started from a random tree and run for 5.0x107 generations with every 1000th generation sampled from the chain. MCMC stationarity was checked as suggested in Nylander et al. [61]. All sample points prior to reaching the plateau phase were discarded as “burn-in”, and the remaining trees were combined to find the a posteriori probability of phylogeny. The analyses were repeated four times to confirm that they all converged on the same results [62].

Species trees were reconstructed using the Singular Value Decomposition Scores for Species Quartets (SVDquartets) [63] and species tree reconstruction in BEAST v2.4.8 (*BEAST) [30,64].

The SVDquartets method infers relationships among quartets of taxa under a coalescent model and then estimates the species tree using a quartet assembly method [63,65]. We evaluated all the possible quartets from the concatenated data set using SVDquartets module implemented in PAUP* v4.0a [66]. Quartet’s Fiduccia and Mattheyses algorithm [67] and multispecies coalescent options were used to infer species tree from the quartets. We used nonparametric bootstrap with 100 replicates to assess the variability in the estimated tree [59].

For *BEAST, multi-species coalescent module implemented in BEAST [30,64] and concatenated dataset were used. We set the partition scheme and models found by Partitionfinder. Mutation rates, clock models, and tree priors were the same as detailed in divergence time estimates section (see below). MCMC were run three times for 5.0x107 generations each, logging tree parameters every 50,000 generations. Posterior distribution was summarized with Densitree v2.01 [64]. Chain mixing, convergence, and a posteriori probability were estimated in the same way as the Bayesian analyses described above.

Species delimitation analyses

Two single-locus analyses, Bayesian General Mixed Yule Coalescent model (bGMYC; [27,68]) and multi-rate Poisson Tree Processes (mPTP; [69]) were performed on mitochondrial dataset. The GMYC model distinguishes between intraspecific (coalescent process) and interspecific (Yule process) branching events on a phylogenetic tree [29]. We used the last 100 trees sampled from the posterior distribution of a Bayesian analysis for mitochondrial sequences (detailed in next section). Bayesian GMYC analyses were assessed using the R package bGMYC, where each tree was ran for 50,000 generations, discarding the first 40,000 generations as burn-in and using thinning intervals of 100 generations (as recommended by Reid and Carstens [70]). The threshold parameter priors (t1 and t2) were set at 2 and 170, and the starting parameter value was set at 25.

mPTP is a phylogeny-aware approach that delimits species assuming a constant speciation rate with different intraspecific coalescent rates [69]. For this analysis, a tree obtained with mitochondrial dataset in MrBayes was used as input on the web server (http://mptp.h-its.org//tree).

Four multi-locus coalescent-based methods were applied to species delimitation: Tree Estimation using Maximum Likelihood, (STEM; [18,21]), Bayesian Species Delimitation (BPP; [26,71]), Multi-locus Species Delimitation using a Trinomial Distribution Model (Tr2; [72]), and Bayes Factor Delimitation (BFD; [73]). As required by these methods, a set of analyses assigning individuals to a series of species categories was performed (delimitation scenarios).

STEM analysis followed Harrington and Near [31]. ML scores for each species tree were generated with STEM v2.0 [21] and evaluated using the information-theoretic approach outlined by Carstens and Dewey [18].

BPP analysis was applied using Bayesian Phylogenetics and Phylogeography software (BPP v.2.2; [26,71]). We used A10 mode, which delimits species using a user-specified guide tree (species delimitation = 1, species tree = 0). The species tree obtained with *BEAST was used as guide tree. Population size parameters (θs) and divergence time at the root of the species tree (τ0) were estimated using A00 mode [71], while the other divergence time parameters were considered as the Dirichlet prior ([24]: equation 2). Each analysis was run four times to confirm consistency among runs. Following a conservative approach, only the speciation events supported by probabilities larger or equal to 0.99 were considered for species delimitation.

The Tr2 analysis followed Fujisawa et al. [72]. Gene trees were obtained in GARLI and the polytomies were resolved using internode branch lengths of 1.0x10-8 in Mesquite v2.75 [74].

For the BFD analysis, we reconstructed a species tree for each delimitation scenario using BEAST, as it was detailed in phylogenetic analyses section (see above). After the standard MCMC chain has finished, a marginal likelihood estimation (MLE) was performed for each species tree, using both path sampling and stepping-stone via an additional run of ten million generations of 100 path-steps (1,000 million generations). Subsequently, the Bayes factor between delimitation scenarios was calculated using MLEs [73] and evaluated using the framework of Kass and Raftery [75].

The taxonomic index of congruence (Ctax) between pairs of species delimitation methods was estimated following the Miralles and Vences’ protocol [34]. In order to access most congruent species delimitation approaches, the mean Ctax value for each method was also estimated.

Divergence time estimates

Divergence times were estimated with concatenated mitochondrial and nuclear dataset using the Bayesian method (BEAST v2.4.8; [64]). We used Neobatrachian mutation rates of 0.291037% and 0.374114% per million years for COI and POMC, respectively [76]. The mutation rates from the other markers were estimated using as prior nuclear or mitochondrial rates for all genes as reported by Irrisarri et al. [76] (0.379173% and 0.075283%, respectively). Partitionfinder provided nucleotide substitution models. LogNormal relaxed clock model and birth-death process as tree prior were used. Bayes factor analysis [77] indicated that this setting received decisive support compared with other models and tree priors available in BEAST. Markov chains in BEAST were initialized from the tree obtained from species tree analyses to calculate posterior parameter distributions, including the tree topology and divergence times. We ran this analysis for 5x107 generations, and sampling every 1000th generation. The first 10% of samples were discarded as “burn-in”, and we estimated convergence to the stationary distribution and acceptable mixing using Tracer v1.6 [78]. An additional BEAST analysis was carried out only with mitochondrial dataset using the same setting to obtain the last 100 trees. These trees were used as input in bGMYC (see section above).

Results

Phylogenetic patterns in E. roseus group

We aligned the five DNA markers for a total of 2576 sites, 858 were variable and 700 were phylogenetically informative. Three of these markers corresponded to mitochondrial dataset with a total of 1799 nucleotide sites, 750 variable, and 629 phylogenetically informative (see information for each marker in S2 Table). Evolutionary models and partitioning strategy obtained in Partitionfinder are also indicated in supplementary data (S2 Table).

Phylogenetic analyses using concatenated mitochondrial and nuclear sequences recovered three main well-supported clades corresponding to Clade A (including E. insularis and E. migueli), Clade B (E. roseus), and Clade C (E. calcaratus) (Fig 2). Although ML and Bayesian analyses recovered to B and C were sister clades, phylogenetic relationships among these clades received low support (Fig 2). Within these clades it is possible to recognize nine highly supported monophyletic lineages (Fig 2; Bootstrap >70, PP>0.99, lineages 1–9). The phylogenetic relationships among Eupsophus species using only mitochondrial dataset recovered the same pattern described for the concatenated dataset (S1A Fig), while nuclear dataset analyses showed a basal polytomy where only lineages 1 and 8 were resolved (Bootstrap >80, PP>0.99; S1B Fig, blue arrows). Other clades exhibited only high posterior probability support in nuclear analyses. For example the clade composed by four individuals from lineage 4, and 21 individuals from lineage 9 (PP>0.99; S1B Fig, red arrow). Nevertheless, low bootstrap support and low variation detected for nuclear markers (S2 Table) prevent us to discuss such a mito-nuclear discordance.

Fig 2

Phylogenetic relationships among Eupsophus species.

This maximum likelihood (ML) tree was reconstructed using concatenated nuclear and mitochondrial data set. Topologies obtained by ML and Bayesian inference were similar. Numbers above branches represent bootstrap scores and Bayesian posterior probabilities. Isolate numbers consist of the species abbreviation (E. roseus: ER, E. migueli: EM, E. insularis: EI, and E. calcaratus: EC), locality abbreviation listed in S1 Table, and field number. Major clades (A, B, and C) and lineages (1–9) of Eupsophus are indicated.

The most congruent result among single- and multi-locus analyses recognized nine monophyletic lineages as different species (Fig 3; mean Ctax = 0.69, see all Ctax values in S3 Table). These nine lineages were the same ones recovered in the phylogenetic analyses and were also supported in the consensus tree from the SVDquartets analysis (Fig 3; Bootstrap >70). Taking into consideration both the geographical distribution (Fig 1) and phylogenetic analyses of Blotto et al. [36], these lineages corresponded to the formerly eight Eupsophus species of the roseus group: E. altor, E. migueli, E. insularis, E. contulmoensis, E. nahuelbutensis, E. septentrionalis, E. roseus, E calcaratus, plus a lineage composed by specimens from the locality of Villarrica, hereafter referred to as Eupsophus sp. (Fig 3).

Fig 3

SVDquartets and species delimitation analyses.

Majority-rule consensus tree from the SVDquartets analysis. Nodal support values are bootstrap proportions. Bars on the right of the tree indicate the species limits as proposed by bGMYC, mPTP, STEM, BPP, Tr2 and BFD analyses. All analyses were carried out with mitochondrial and nuclear loci, except bGMYC and mPTP which used only mitochondrial data set. Limits of formerly Eupsophus species and putative species from Villarrica (Eupsophus sp.) are indicated with different colors on the branches of the tree and with square bracket on the right of the bars. These limits correspond to the most congruent species delimitation scenario (see S3 Table).

Bayesian GMYC analyses detected more than one species in these nine lineages except in E. insularis and E. contulmoensis (Fig 3). Multi rate PTP detected six species corresponding to E. altor, E. migueli, E. insularis, E. contulmoensis, E. nahuelbutensis, and E. septentrionalis lineages, and more than one species in E. roseus and E calcaratus lineages (Fig 3). The nine-species scenario (Fig 4A, gray cell) was the highest supported in BPP and Tr2 analyses (Fig 4B, black arrows, scenario 12). For the STEM analysis the eight-species scenario, where Eupsophus sp. and E. roseus represent a single species, was the highest supported (Fig 4A, scenario 11). Nevertheless, among the other species delimitation scenarios, the STEM analysis greatly favored a nine-species delimitation scenario (Fig 4B, S4 Table). Highest MLEs in BFD analysis were obtained for eight-species scenario, where E. altor and E. migueli corresponded to one species (Fig 4, scenario 10). In this case, Bayes factor comparisons were greater than two, which allowed us to choose that better scenario (S5 Table). Nevertheless, comparisons with some scenarios including that of nine-species were around four, which indicate non-strong or decisive support to the best model (S5 Table). Other possible scenarios, including that proposed by Correa et al. [37] (scenario 3), were lowly supported for all multi-locus analyses (Fig 4).

Fig 4

Multi-locus species delimitation analyses.

A) species delimitation scenarios. Specimens were assigned to the delimited species indicated in Fig 3. Abbreviations within parenthesis indicate the grouping tested in each scenario. E. roseus: ER, E. migueli: EM, E. insularis: EI, and E. calcaratus: EC, E. altor: EA, E. contulmoensis: ECO, Eupsophus sp.: EV, E. nahuelbutensis: EN, E. septentrionalis: ES. Some abbreviated localities from S1 Table were added to species abbreviation to indicate a specific locality grouping. The most congruent scenario is indicated in gray. B) probability, marginal likelihood (MLE), or score values generated for each scenario using different species delimitation approaches. Black arrow indicates the credible species hypotheses. For Tr2 lowest score indicates the better-delimited scenario. For STEM and BFD were plotted model probabilities and MLE values using stepping-stone sampling, respectively (see S4 and S5 Tables).

Species tree and divergence time estimates among Eupsophus species

Species tree reconstructions in *BEAST and SVDquartets, using the nine lineages (= species), recovered similar phylogenetic relationships to the Bayesian and ML analyses (Fig 5). Under this scenario, E. calcaratus diverged early in Eupsophus radiation for both the species tree and the divergence time tree. Overlaying posterior sets of trees generated in BEAST and plotted by DensiTree supported this topology (Fig 5). Thus, we decided to used consensus species tree as a prior to estimate the divergence times among Eupsophus species (Fig 5, in blue).

Fig 5

Species tree and divergence times of Eupsophus.

This cladogram illustrates the posterior distribution of the species trees inferred with BEAST based on the most congruent species delimitation scenario (Figs 3 and 4, S2 Table). High color density is indicative of areas in the species trees with high topology agreement. Different colors represent different topologies. Consensus species tree are colored in blue. Nodal values are Bayesian posterior probability (BEAST) and bootstrap proportions (SVDquartets). Mean divergence dates in million years and 95% credible intervals are indicated (below the support values).

The age of crown-group Eupsophus and the origin of E. calcaratus are estimated at 0.396 (0.351–0.442) Myr. Eupsophus insularis diverged at 0.268 (0.230–0.308) Myr, while E. altor and E. migueli at 0.096 (0.077–0.116) Myr (Fig 5). The split between E. roseus and Eupsophus sp. /E. contulmoensis, E. nahuelbutensis, and E. septentrionalis was around 0.134 (0.114–0.154) Myr. The divergence between E. roseus and Eupsophus sp. was estimated at 0.088 (0.072–0.106) Myr. Eupsophus septentrionalis diverged at 0.111 (0.193–0.131) Myr, followed of E. contulmoensis and E. nahuelbutensis at 0.054 (0.041–0.067) Myr (Fig 5).

Discussion

Species delimitation in the Eupsophus roseus group

The most congruent species delimitation results detected nine species in the E. roseus group, and eight of them (namely E. altor, E. calcaratus, E. contulmoensis, E. insularis, E. migueli, E. nahuelbutensis, E. roseus, and E. septentrionalis) were concordant with taxonomic proposals of the last decades [36,38,49,79-83]. Although, gaps in morphological, geographic, cytogenetic, bioacoustic, and behavioral information prevent us to carry out a protocol for integrative taxonomy, our molecular approach is concordantly with integrative studies available for some species as E. altor [81]. Moreover, we provided molecular evidences for separation of nine evolutionary lineages, a key step to carry out future work protocols under an integrative taxonomical approach [84].

The highest level of congruence was obtained with BPP and Tr2 methods (mean Ctax = 0.69; nine species), followed by STEM, and BFD (mean Ctax = 0.63; eight species; Figs 3 and 4, S3 Table). Although, Eupsophus sp. and E. roseus clades were recovered as a single species by STEM, they were recovered as different species by BPP, Tr2, mPTP, and BFD analyses, similar to the case of E. migueli and E. altor which were recovered as a single species by BFD but as two different species in the other analyses. Therefore, the greatest congruence indicated that Clade B is composed by five different species (Eupsophus sp., E. roseus, E. nahuelbutensis, E. contulmoensis and E. septentrionalis), while Clade A is composed by three (E. altor, E. migueli, and E. insularis) as it was suggested in previous works [36,81]. The differences among the results of these species delimitation methods could be derived from their different sensibility to the ratio of population size to divergence time, as reported between BPP and bPTP [17]. Hence the importance of carrying out several species delimitation methods to examine whether the proposed groups are consistently recovered with different algorithms [17,11]. This was evident when we compared results from multi-locus analyses with bGMYC result (mean Ctax = 0.27), which overestimated the number of species in all lineages except in E. insularis and E. contulmoensis (Fig 3). It is known that bGMYC has shortcomings when datasets consist of few putative species [85] and cannot be used as sufficient evidence for evaluating the specific status without additional data or analyses [86]. Moreover, this method tends to overestimate the number of species when the ancestral polymorphism is low [87]. Therefore, rather than using this method as a species delimitation approach, we used it to obtain alternative scenarios to be tested with multi-locus analyses (e.g. scenario 13, Fig 4).

Our delimitation results did not agree with a recent hypothesis [37], which would be related to the use of different molecular markers and species delimitation analyses. Three of our markers were found to be highly variables (Cyt b, COI, D-loop), while two were conservative (POMC and CRYBA1; see S2 Table). Thus, we use at least three strong markers (sequences with many polymorphic sites), a key aspect to carry out coalescent analyses when less than ten markers are used [88]. On the other hand, we used several multi-locus coalescent methods to delimitate species (BPP, STEM, R2, and BFD), while Correa et al. [37] based its inferences in single-locus analyses (bGMYC, mPTP, and Automatic Barcode Gap Discovery, ABGD). In this sense, the two groups of synonymized species were recovered as two species in analyses of mPTP (using mitochondrial data set) and ABGD (using mitochondrial + nuclear data set) performed by these authors [37]. The ABGD method is based on genetic distances computed from a single-locus (COI) and requires a priori specification of an intraspecific distance threshold [89]. The robustness and accuracy of coalescent approaches over distance methods are well known, partly because the latter do not appeal to an explicit species concept [17,90]. Therefore, we decided not to include ABGD in our main species delimitation analyses. Nevertheless, we conducted ABGD analyses using our COI and concatenated datasets, obtaining different results (see S1 File). In this regard, the use of two potential barcode gaps allowed us to detected nine and five groups with COI, while seven and four groups were obtained with concatenated dataset. Consequently, ABGD results can be influenced by the application of a method designed for single-locus (DNA barcoding) to concatenated dataset, as well as by the a priori election of distance threshold. Moreover, ABGD analysis underestimated species diversity among species with low divergence [89,91]. Thus, ABGD tool is recommended as a first grouping hypothesis but not as robust and definitive species delimitation proof [89].

Phylogenetic relationships and divergence time in the Eupsophus roseus group

Monophyly of E. roseus group and its nine delimited species was strongly supported, concordant with previous analyses (Fig 2; [36,49]). Although the early divergence of E. calcaratus was not strongly supported in Bayesian, ML, and SVDquartet approaches; our analyses resolved all other interspecific relationships among delimited species (Figs 2 and 3). In fact, the plot of overlying posterior sets of species trees (Fig 5) showed few alternative interspecific relationships. One example of this, is the early divergence of E. septentrionalis within Clade B, which was also recovered by Blotto et al. [36] and Suárez-Villota et al. [49] (Fig 5, in red).

Phylogenetic and species delimitation analyses recognized to Eupsophus sp. as a distinct species (Figs 3 and 4). In fact, SVDquartet analysis detected this clade with greater support than other well-defined species such as E. insularis (Fig 3; bootstrap: 95%), and high probabilities were detected in single- and multi-locus species delimitation analyses (Figs 3 and 4). These results are concordant with previous works suggesting a species-level for this lineage [55]. Although Correa et al. [37] also detected a close phylogenetic relationship between Villarrica and E. roseus specimens, they considered the three specimens from this locality within the E. roseus diversity. We sampled 17 specimens from this locality and they were monophyletic with high support (Fig 2; Bootstrap: 94.4, PP: 0.99). Additionally, we did not detect syntopy instances in Villarrica, which could result in the recovery of specimens from other localities within the Villarrica clade (i.e. interpopulational paraphyly). This paraphyletic pattern is common for localities within the E. roseus lineage, an additional support to consider the possibility that Villarrica specimens do not belong to E. roseus species. For example, specimens from Fundo Santa María (FS) are recovered with specimens from other localities [e.g. Mafil (MA), Llancahue (LA)], in several highly supported clades within the E. roseus lineage (Fig 2).

We used the divergence rate in agreement with estimates for several other Neobatrachian species [76]. We fully recognize that this approach is far from ideal with several potential sources of error [92], but a beginning exploration of evolutionary histories of these endemic Patagonian species will in our view benefit from provisional estimates. Under this assumption, most of the delimited species from the E. roseus group diverged from 0.134 to 0.054 Mya during the Valdivian interglacial [93], except E. calcaratus and E. insularis whose origin is older (before of the Last Southern Patagonian glaciation, 0.18 Mya). The oldest deposits of Mocha Island are dated from the Eocene and Miocene [94] whereas extensive terraces from Pliocene and Pleistocene characterize more recent settings [95]. Although the origin of E. insularis in the Mocha Island remain unknown, these large terraces might have been a suitable habitat for both its settlement and differentiation from the continental Eupsophus species. Anyway, it is possible that all species lived during Valdivia interglacial and were subsequently affected by the Last Glacial Maximum (LGM, 0.020–0.014 Mya; [96,97]). Valdivia interglacial was characterized by the presence of North Patagonian forests and Valdivian rainforests [98], which are habitats associated to Eupsophus species [44,81]. These suitable Late Pleistocene habitats for Eupsophus species were probably contracted during periods of glacial advance, whereas distributional range shifted during glacial retreats and warming. Therefore, it is possible to hypothesize a wide distribution of Eupsophus species during the interglacial, followed by restricted distribution in refugia during the LGM. Consequently, current restricted distribution of some Eupsophus species (e. g. E. migueli, E. altor, E. contulmoensis, E. nahuelbutensis, Eupsophus sp., E. septentrionalis) could be related with Pleistocene cycling events. In fact, geographical isolation effect of Quaternary cycling events over other vertebrate species has been hypothesized [99-101].

Finally, the lineage represented by the Villarrica specimens (Eupsophus sp.) diverged from E. roseus at ~ 0.088 Mya (Fig 5). Under this temporal scenario it is possible that this lineage lived during the interglacial and was subsequently affected by LGM. A central east colonization of an ancestral E. roseus population could have given rise to Eupsophus sp. during warmer interglacial conditions. In this sense, this putative species probably represents a remnant lineage left behind in central-west Chilean refugia present during LGM. In short, isolation during LGM, the monophyly, and coalescent species delimitation suggest taxonomic differentiation of the Villarrica specimens.

Using new molecular datasets and coalescent analyses, our approach revitalizes in an independent way the hypothesis that the E. roseus group is composed of eight species. Moreover, we suggest the taxonomic differentiation for the Villarrica specimens. Finally, we suggest filling bioacoustic, morphological, behavioral, and karyotypic data gaps for a deep Eupsophus revision.

Sampling locations of Eupsophus species.

Coordinates, sample size (N), corresponding species according to Frost [35] and map number from Fig 1 are indicated. Species used as outgroup are also listed (gray cells).

(DOC)

Click here for additional data file.

Sites characterization, partitioning schemes, and nucleotide substitution models for sequences used in this study.

Conservative (C), variable (V), informative (I) and total sites for each marker are indicated. Partitioning schemes and nucleotide substitution models were determined using Partitionfinder, version 2.1.1 [56].

(DOC)

Click here for additional data file.

Taxonomic index of congruence (Ctax) calculated for each pair of approaches.

Mean of all the Ctax values obtained involving a given approach (Mean Ctax) and total number of species supported by each approach (sp.) is indicated. Species delimitation approaches: Bayesian General Mixed Yule Coalescent model (bGMYC), multi-rate Poisson Tree Processes (mPTP), Tree Estimation using Maximum likelihood, (STEM), Bayesian Species Delimitation (BPP), Multi-locus Species Delimitation using a Trinomial Distribution Model (Tr2), and Bayes factor delimitation (BFD).

(DOC)

Click here for additional data file.

Likelihood scores and Akaike’s information criterion (AIC) results for STEM analysis (see Carstens and Dewey [18]).

Species delimitation scenarios for Eupsophus species are indicated in Fig 4A. Species number (sp.), Log-likelihood of the species tree (−lnL), number of parameters (k), AIC, AIC difference (Δi), the relative likelihood of model given the data (L), and the model probabilities (w) are indicated. Note the proximity between–lnL from scenario 11 and 12.

(DOC)

Click here for additional data file.

Bayes factor delimitation results.

Marginal likelihood (MLE) and Bayes factor estimates for species delimitation scenarios indicated in Fig 4A. Species number (sp.) as well as values using path (PS) and stepping-stone (SS) sampling are indicated.

(DOC)

Click here for additional data file.

Maximum likelihood trees using A) mitochondrial ( Bayesian analyses recovered similar topologies in both cases. Support values (bootstrap and posterior probabilities) are shown above branches. Blue arrows indicate lineages 1 and 8 shown in Fig 2, red arrow indicates a clade supported only in Bayesian analyses.

(TIF)

Click here for additional data file.

Automatic Barcode Gap Discovery (ABGD) results using A) COI and B) concatenated (D-loop, Cytb, COI, POMC and CRYBA1) datasets.

(DOCX)

Click here for additional data file.