Hierarchies of evolutionary radiation in the world’s most species rich vertebrate group, the Neotropical Pristimantis leaf litter frogs

Abstract

The Neotropical leaf litter frog genus Pristimantis is very species-rich, with 526 species described to date, but the full extent of its diversity is much higher and remains unknown. This study explores the phylogenetic processes and resulting evolutionary patterns of diversification in Pristimantis. Given the well-recognised failure of morphology- and community-based species groups to describe diversity within the genus, we apply a new test for the presence and phylogenetic distribution of higher evolutionary units. We developed a phylogeny based on 260 individuals encompassing 149 Pristimantis presumed species, sampled at mitochondrial and nuclear genes (3718 base pair alignment), combining new and available sequence data. Our phylogeny broadly agrees with previous studies, both in topology and age estimates, with the origin of Pristimantis at 28.97 (95% HDP =21.59 – 37.33) million years ago (MYA). New taxa that we add to the genus, which had not previously been included in Pristimantis phylogenies, suggest considerable diversity remains to be described. We assessed patterns of lineage origin and recovered 14 most likely (95% CI: 13–19) phylogenetic clusters or higher evolutionary significant units (hESUs) within Pristimantis. Diversification rates decrease towards the present following a density-dependent pattern for Pristimantis overall and for most hESU clusters, reflecting historical evolutionary radiation. The timing of diversification suggests that geological events in the Miocene, such as Andes orogenesis and Pebas system formation and drainage, may have had a direct or indirect impact on the evolution of Pristimantis and thus contributed to the origins of evolutionary independent phylogenetic clusters.

Key words:

• anurans
• diversification rate
• evolutionary radiation
• higher evolutionary significant units (hESUs)
• Neotropics
• species richness

Introduction

The tropical forests of South America boast a disproportionately high percentage of the world’s biodiversity (Myers, Mittermeier, Mittermeier, da Fonseca, & Kent, 2000; Wilson, 1988), with mountainous regions being particularly speciose (Cracraft, 1985; Gentry, 1992; Hernandez & Young, 2007). Investigation into the high levels of diversity and endemism in Neotropical amphibians, facilitated largely by traits like low vagility and tendency toward philopatry, has the potential to greatly increase our understanding of the diversification of all South American biota (Duellman, 1999; García-R et al., 2012; Loyola, Kubota, da Fonseca, Lewinsohn, 2009; Zeisset & Beebee, 2008). To date, phylogenetic and biogeographical studies have supported the hypothesis that Andean orogeny (4.5–23 MYA), which changed the course of the Amazon River and formed ridges, as well as the formation of the Pebas wetlands system (8–23 MYA), caused rapid diversification across many taxonomic groups, including amphibians (Antonelli et al., 2010; Castroviejo-Fisher, Guayasamin, Gonzalez-Voyer, & Vilà, 2014; Ceccarelli et al., 2016; De-Silva, Elias, Willmott, Mallet, & Day, 2016; Hoorn et al., 2010; Lougheed, Gascon, Jones, Bogart, & Boag, 1999; Muñoz-Ortiz, Velásquez-Álvarez, Guarnizo, & Crawford, 2015; Noonan & Wray, 2006; Sanín et al., 2016; Santos et al., 2009). Furthermore, disproportionately high species richness is found in the montane regions, for example in dart-poison frogs (Santos et al., 2009) and glass frogs (Castroviejo-Fisher et al., 2014; Hutter, Guayasamin, & Wiens, 2013), as well as birds (Kattan & Franco, 2004).

The widespread and understudied Neotropical amphibian genus Pristimantis (Anura, Terraranae, Craugastoridae) is the most speciose of terrestrial vertebrates at approximately 526 currently described species (Frost, 2018). There has been an increase in the description of new Pristimantis species in recent years due to molecular techniques revealing cryptic diversity (e.g. Elmer, Dávila, & Lougheed, 2007b; Elmer and Cannatella, 2008; Hutter & Guayasamin, 2015; Ortega-Andrade et al., 2015) and increased sampling in previously uncharted areas (Barrio-Amorós, Mesa, Brewer-Carías, & McDiarmid, 2010; Kaiser et al., 2015). The taxonomy of this group is complex and problematic due to the prevalence of high genetic variation within putative species (Crawford, 2003) combined with a lack of diagnostic morphological characters (Duellman & Lehr, 2009; Padial, Grant, & Frost, 2014). Several Pristimantis species are considered to be widespread over large areas, which conflicts with both current diversification hypotheses and the ecological traits amphibians exhibit, which have tended to facilitate speciation at more restricted geographic scales (Crawford 2003; Fouquet et al., 2007; Funk, Caminer, & Ron, 2012). Therefore, researchers have suggested that these widely distributed species more likely represent cryptic species complexes that have not been identified due to a lack of sampling (Elmer, Dávila, & Lougheed, 2007a; Funk et al., 2012). Recognising true species richness and cryptic species has major implications for fundamental understandings of diversification and biogeography, and critical utility for conservation management of biodiversity and of landscapes (Fouquet et al., 2007; Funk et al., 2012).

The basis for the exceptional species richness within Pristimantis and other terraranans has not been identified rigorously. Several hypotheses have been proposed, including specific phenotypic characteristics (Gonzalez-Voyer, Padial, Castroviejo-Fisher, & De la Riva, 2011), and the less restrictive direct development mode of reproduction allowing populations to disperse more widely, occupy a range of habitats (Hedges, Duellman, & Heinicke, 2008), and establish a rapid accumulation of population genetic isolation across the landscape (Elmer et al., 2007b; Prohl, Ron, & Ryan, 2010). Phylogenetic studies have placed the geographic origin of Pristimantis within South America (Heinicke, Duellman, & Hedges, 2007; Pinto-Sánchez et al., 2012), with multiple independent dispersal events into Central America (Hedges et al., 2008; Pinto-Sánchez et al., 2012). More recently, Mendoza, Ospina, Cárdenas-Henao, & Garcia-R (2015) recognised a mid-elevational band (1000–3000 masl) in the North-western Andes (Colombia and Ecuador) as the focal point for the origin and radiation of Pristimantis, giving further support to the hypothesis that the mountainous north-western part of the continent acts as a species pump. However, the relationship of evolutionary rates within and across these Pristimantis clades has not been assessed.

The genus Pristimantis was divided into species series and species groups based on extensive morphological work by Lynch & Duellman (1997). Later, many of these groups were not found to be monophyletic in subsequent molecular analyses (Hedges et al., 2008; Mendoza et al., 2015; Padial et al., 2014; Pinto-Sánchez et al., 2012). However, researchers were reluctant to revise the taxonomy of this genus largely due to a lack of comprehensive taxon sampling (Hedges et al., 2008). Padial et al. (2014) reclassified Pristimantis into 11 species groups based on a parsimony-based phylogeny; yet given the extremely high species richness, this was based on less than a quarter of current species within the genus, suggesting there are still many unknown assignments. Without the cohesion of designated species groups, Pristimantis currently lacks any empirically derived units of diversity to focus and contextualise taxonomy and systematics. Such delimitations are key to further research on speciation, biogeography, and processes promoting diversification.

Our study explores the phylogenetic processes associated with diversification of the world’s most species-rich vertebrate genus, Pristimantis. We accomplish this by combining new multi-locus sequence data and new species with published datasets to achieve three goals. First, our study produced a dense, time calibrated phylogeny of the genus Pristimantis. Second, given the well-recognised failure of species groups to delineate natural genetic groups within the genus, we test for the presence and phylogenetic distribution of key evolutionary units within Pristimantis. This approach de-emphasises both the role of described species, which we know are far from complete, and the unevenness in geographical sampling, which may be biasing our current understanding of species richness across groups of Pristimantis. Finally, we highlight some taxonomic issues within our Pristimantis phylogeny with regards to geographic and species sampling, mis-identification, and likely instances of cryptic diversity, noting some taxonomic and geographic areas for further investigation.

Materials and methods

Species and genes analysed

The dataset comprises 260 individuals encompassing 149 recognised Pristimantis species, which represents around 30% of the currently described species (526 species, Frost, 2018). We combined 107 new and unpublished sequences (28 12S, 44 16S, 28 cytb and seven Rag1) for 50 individuals from Ecuador (sourced from collections at Museo de Zoología of Pontificia Universidad Católica del Ecuador (QCAZ); molecular methods follow) and 646 published sequences (170 for 12S, 210 for 16S, 70 for COI, 25 for cytb, 84 for Rag1 and 87 Tyr) accessed from GenBank (Table S1, see supplemental material online, which is available from the article’s Taylor & Francis Online page at https://doi.org10.1080/14772000.2018.1503202). These 210 individuals were originally sourced from the following countries Ecuador (n = 62), Colombia (n = 45), Peru (n = 34), Panama (n = 23), Costa Rica (n = 7), Guyana (n = 6), Trinidad and Tobago (n = 3), Honduras (n = 1) and Venezuela (n = 1), with 28 individuals having unknown location. We rooted the phylogeny using molecular data from four outgroups from the closely related genera Craugastor and Diasporus (following Pinto-Sánchez et al., 2012).

Molecular techniques for new data

Genomic DNA was extracted from ethanol-preserved tissue using either standard phenol-chloroform methods or a Qiagen DNEasy kit, following the standard protocol. The following primers were used for PCR: 12S (12Sa-L, 12Sb-H; Goebel Donnelly & Atz, 1999), 16S (16Sar-L, 16Sbr-H; Goebel et al., 1999), cytb (MVZ15L, MVZ16H; Goebel et al., 1999) and Rag1 (R182_Rag, R270_Rag; Heinicke et al., 2007).

Cocktails for each mitochondrial PCR reaction contained: 0.3 μM of the relevant forward and reverse primer, 5μL PCR enhancing buffer (2.5 mM MgCl₂, 10 mM Tris pH 8.4, 50 mM KCl, 0.02 mg BSA, 0.01% gelatin), 0.3 mM dNTP, 0.625 units Taq DNA polymerase (Fermentas Life Sciences, UK) and approximately 1 to 3 ng DNA (in 1:10 dilutions). An Applied Biosystems 2700 or 9700 thermocycler was used for PCR. Cycling parameters for 16S were: 94° C for 2 min, 35 cycles of denaturing 94° C for 30 sec, annealing 60° C for 20 sec and extending 72° C for 20 sec, a final extension at 72° C for 5 min followed by extended cooling at 10 or 4° C. The parameters for 12S were the same as 16S, except with an annealing temperature of 50° C. The cytochrome b parameters were: initial denaturing at 92° C for 3 min, 38 cycles denaturing 92° C for 1 min, annealing 51° C for 1 min and extension 72° C for 1 min, a final extension of 72° C for 5 min followed by extended cooling at 4 or 10° C. PCR product was purified using Pall AcroPrep 96 Filter Plates.

Rag1 PCR cocktail contained: 2µl of 10 x PCR Rxn Buffer (Invitrogen), 1µl of MgCl₂ (50mM, Invitrogen), 2µl of dNTPs (10mM), 0.2µl of both forward and reverse primer (10 pmole/µl), 0.2µl of Taq polymerase (Invitrogen) and 1 to 3 ng of DNA. The PCR protocol was as follows: an initial denaturating at 94 °C for 5 min, 35 cycles of denaturating at 94 °C for 30 sec, annealing at 60 °C for 30 sec and extension at 72 °C for 60 sec, a final extension at 72 °C for 7 min, followed by an extended cooling at 10 °C. The PCR products were purified using an ExoSap kit (Thermo Fisher) following the standard protocol. All mtDNA and nuclear reaction sets included a negative control.

All genes were sequenced using Applied Biosystems Big-Dye Ver 3.1 chemistry on an Applied Biosystems 3730 automated capillary DNA sequencer. Samples were sequenced in the forward and reverse directions, except for a subset of six cytb, one 12S, three 16S and five Rag1 samples that were only sequenced one direction.

Phylogenetic analyses

The sequence chromatograms were cleaned and compiled into consensus (when sequenced bi-directionally) using Geneious v6.1.4 (Kearse et al., 2012, http://www.geneious.com). The sequences for all six genes were aligned separately using the ‘ClustalW’ alignment option within Geneious with the default settings, and then concatenated. Before concatenation, 12S and 16S sequences were separately analysed with Gblocks v0.91 on the server (Castresana, 2000, http://molevol.cmima.csic.es/castresana/Gblocks_server.html) to remove unalignable regions using the following parameters: the minimum length for a block set at 5, allowed gap positions in half positions and the rest of settings kept as default (Massana et al., 2004). The final matrix had a total of 58% missing data.

Models of evolution were evaluated to select the best nucleotide substitution models to be used in further analysis. PartitionFinder v1.1.1 (Lanfear, Calcott, Ho, & Guindon, 2012) was used to identify the best substitution model and partition scheme on the concatenated dataset for each gene and codon position separately using the following parameters: branch lengths = linked, models = beast (models specific to the BEAST v1.7.5 program), model selection = Bayesian Information Criterion (BIC), and schemes = greedy. A six-partition scheme was found to be the best overall strategy: 16S, 12S and cytb_1 (GTR + I + G); cytb_2 (TrN + G); cytb_3 (TrN + G); COI_1, COI_2, COI_3 (GTR + I + G); Rag1_1, Rag1_2, Tyr_1, Tyr_2 (GTR + I + G) and Rag1_3 and Tyr_3 (HKY + G) and used for the phylogenetic analysis.

A time-calibrated phylogeny was inferred using Bayesian methods in BEAST v1.8 (Drummond, Suchard, Xie, & Rambaut, 2012) using secondary calibrations with the three age estimates used in Mendoza et al. (2015), sourced originally from Heinicke et al. (2007) and Pinto-Sánchez et al. (2012). A normal prior distribution was used for all calibrations. BEAUti v.1.8 (Drummond et al. 2012) was used to prepare the xml file used by BEAST. A mean of 24.45 and standard deviation (SD) of 5.0 (2.5% quantile = 14.65, 75% quantile = 34.25) Ma was selected for the crown age of Pristimantis (following Mendoza et al. 2015). For the crown age of the monophyletic P. pardalis group (Wang, Crawford, & Bermingham, 2008) we used a mean of 8.6 and SD of 1.8 (2.5% quantile = 5.07, 75% quantile = 12.12) Ma, and for the crown age of P. taeniatus we chose a mean of 8.3 and SD of 1.76 (2.5% quantile = 4.85, 75% quantile = 11.75) Ma (following Mendoza et al. 2015). We used an uncorrelated lognormal relaxed clock model and the Yule species tree prior, following the methods of Pinto-Sánchez et al. (2012). All other priors were left as default for BEAUti v.1.8 (ucld.mean = fixed value of 1, ucld.stdev = exponential with bounds 0 to ∞). Gaps were treated as missing data. The BEAST analysis was run for a chain length of 30 million with samples taken every 3,000 generations. From the resulting 10,000 trees, the first 10% was discarded as burn-in.

Higher evolutionary significant units (hESUs)

To test for the presence of evolutionary significant units above the species level (hESUs) across the tree, we employed the Generalised Mixed Yule Coalescent (GMYC) model (Fontaneto et al., 2007; Fujisawa & Barraclough, 2013; Pons et al., 2006), using the R package (R Core Team, 2017) ‘SPLITS’ 1.0 (Ezard, Fujisawa, & Barraclough, 2009). GMYC can determine evolutionary clusters above the species level using a time-calibrated tree, with tips representing species and not individuals (per Humphreys & Barraclough, 2014). The GMYC model estimates the time at which a change in the rate of phylogenetic branching cannot be explained by a single coalescent process. We employed the multiple threshold method (Monaghan et al., 2009) for increase in branching rate implemented in ‘SPLITS’, which has been suggested to be more accurate and robust (Humphreys & Barraclough, 2014; Postaire, Magalon, Bourmaud & Bruggemann, 2016). The tree used in this analysis was our maximum clade credibility (MCC) tree, pruned to include only one individual per species, excluding “aff”, “cf” and “sp” specimens as well as any uncertain identifications (total = 139 taxa included). Even though it was developed for single locus trees, the GMYC model has been shown to work effectively for combined datasets (e.g. Barej, Penner, Schmitz & Rödel, 2015; Jörger, Norenburg, Wilson & Schrödl, 2012; Low et al., 2016). This analysis tests for the presence of hESUs (Humphreys & Barraclough, 2014; Humphreys et al., 2016), which could represent phylogenetic species groups, as opposed to any phenetic species group (e.g. as proposed by Lynch & Duellman, 1997).

Rates of diversification

We applied multiple assessments of diversification rate, for Pristimantis overall and for the identified hESUs. First, we tested specific models of diversification through time for all of Pristimantis and for each hESU using LASER 2.3 (Rabosky, 2006), selecting the best fitting model using Akaike’s Information Criterion. Only hESUs with more than four taxa were used in the analysis (n = 10), as it requires more than three branching events. Second, lineage through time plots for the pruned MCC tree inferred in BEAST were generated in APE 3.09 (Paradis, Claude & Strimmer, 2004), using 1000 trees sampled from the posterior distribution. Third, to test whether diversification rates remained constant over time we applied the gamma statistic (Pybus & Harvey, 2000) with the Monte Carlo Constant Rate (MCCR) test (Pybus & Harvey, 2000) in LASER, which also accommodates incomplete species sampling (i.e., there are 526 Pristimantis species described, while our pruned tree has 139 species) by simulating full phylogenies, which are then randomly pruned to mimic incomplete sampling. Fourth, we used BAMM (Rabosky, 2014) to estimate speciation rates for the whole Pristimantis and each hESU using the MCC tree. BAMM is able to detect shifts in evolutionary rates within a phylogeny, using reversible jump Markov chain Monte Carlo to explore various models of lineage diversification. In addition, BAMM allows users to incorporate knowledge of incomplete taxon sampling in the likelihood calculation, so for the globalSamplingFraction parameter we used a value of 0.30, since our MCC tree includes 139 of the 526 species described for Pristimantis. We used a value of 1 for the prior on the expected number of shifts. Priors for BAMM were generated using the function setBAMMpriors in the R package BAMMtools v 2.1.6 (Rabosky et al., 2014) by providing the MCC tree (priors: lambdaInitPrior = 1.36594224825722; lambdaShiftPrior = 0.0397449192506994; muInitPrior = 1.36594224825722). We ran 5 million generations, sampling parameters every 1000 generations. The output of BAMM was then analysed in R using BAMMtools to calculate mean speciation and extinction rate, bayes factor (BF), and marginal odds-ratios. We discarded the first 30% of generations as burn-in and assessed convergence using CODA v 0.16–1 (Plummer, Best, Cowles & Vines, 2006). All parameters had effective sample sizes >400. We identified the most likely diversification scenarios for Pristimantis using Bayes Factor (BF) and marginal posterior-to-prior odd ratios, calculated speciation rates, and reconstructed rate-through-time curves for the genus and each hESU.

Results

Phylogenetic analyses

The total alignment of six genes was 3718 base pairs (bp) in length, made up of 715 bp for 12S, 503 bp for 16S, 666 bp for cytb, 673 bp for COI, 630 bp for Rag1 and 531 bp for Tyr. Posterior probabilities of ≥0.95 were obtained for 63% of the nodes (163 of 262) in our Bayesian consensus tree, though some nodes are weakly supported (Fig. 1). The origin of Pristimantis was estimated to be 28.97 MYA (95% HDP = 21.59 – 37.4 MYA) (Fig. 1; Fig. S1, see supplemental material online).

Fig. 1 Bayesian consensus time-tree for Pristimantis based on a 6-gene analysis of 264 frogs. The posterior probabilities are shown for each node. Time in millions of years before present and geological eras are noted. New specimens with this study are marked with ^. Full specimen information in Table S1 (see supplemental material online). 95% HPD for dating at nodes are shown in Fig. S1 (see supplemental material online). Calibration nodes noted with black dot (see Methods).

Our phylogeny did not recover as monophyletic any of the 15 species groups or series included in our dataset, with the exception of the peruvianus group (Table S1, see supplemental material online). The largest species group, unistrigatus, was represented by 97 individuals in our dataset and had individuals that appeared in most clades across our tree (Fig. 1).

Our phylogenetic analysis included 37 species that were represented by at least two exemplar individuals. After revisiting the voucher notes and assigning some vouchers with ‘aff’, 22 of these (60%) were found to be monophyletic, including five species for which we added new individuals (P. altamnis, P. lanthanites, P. librarius, P. luscombei, and P. quaquaversus). Eight species that were part of our dataset of new sequences were not recovered as monophyletic (P. altamazonicus, P. conspicillatus, P. croceoinguinis, P. diadematus, P. kichwarum, P. malkini, P. martiae and P. omevividis), along with a further seven species already on GenBank (P. bogotensis, P. eriphus, P. pardalis, P. peruvianus, P. unistrigatus, P. viejas and P. zophus). These likely represent misidentified specimens and cryptic diversity within putative species rather than incomplete lineage sorting (Table S2, see supplemental material online).

Higher evolutionary significant units

The GMYC model test for phylogenetic clusters showed evidence for the presence of 14 clusters (95% confidence interval of 13–19 clusters), or separately evolving units, within the Pristimantis tree (p = 0.004) (Fig. 2, Table 1). This confidence interval is taken from the 95% confidence set of alternative models that identified the phylogenetic clusters, with 14 being the most likely value. Phylogenetic branching rates (a transition from clusters to coalescence; Humphreys and Barraclough 2014) increased at ca. 24.5, 20.5, 16.8 and 15.7 MYA based on the tree and the multiple-threshold GMYC model (Fig. S2, see supplemental material online). Two taxa, P. rozei and P. caprifer placed in a sister relationship, were not included in any cluster. Post hoc assessment suggests that one cluster, hESU A, corresponds to a previously described species group, the peruvianus species complex (as per Hedges et al., 2008; Table S1, see supplemental material online).

Fig. 2 Higher Evolutionary Significant Units (hESU) within Pristimantis. a) The 14 hESUs (A-N) recovered by the GMYC model with the multiple thresholds method (95% CI 13–19 clusters). Two species in black, P. caprifer and P. rozei, were not assigned (n/a). Node numbers indicate support based on the confidence set of models.

Table 1. Summary of higher evolutionary units (hESUs) recovered from time-calibrated Bayesian phylogeny of the pruned Pristimantis tree (139 in-group species; two non-assigned taxa not shown). Number of taxa per hESU, and time of origin (MYA) drawn from the time-calibrated tree.

hESU no.	No. of taxa	Time of origin MYA (95% CI)
A	6	11.75 (7.13–17.2)
B	4	24.51 (16.71–33.32)
C	6	20.53 (13–27.44)
D	11	18.54 (12.5–25.55)
E	3	11.9 (7.11–18.23)
F	15	16.8 (11.63–22.61)
G	4	17.87 (11.42–24.97)
H	22	15.29 (9.85–21.35)
I	17	16.65 (11.3–22.44)
J	7	12.96 (9.93–21.11)
K	7	15.34 (7.89–18.67)
L	3	10.54 (6.59–15.13)
M	11	14.12 (9.81–18.7)
N	21	15.67 (11.25–20.92)

Rates of diversification

Rate-variable models of diversification, which allow diversification rate to change overtime, were preferred over rate-constant models for Pristimantis overall and for all hESUs except one (hESU C pure birth: r = 0.053, AIC 23.90) (Table S3, see supplemental material online) by the LASER analysis. Of those rate-variable models, the density-dependent logistic (DDL) model was the best fit for Pristimantis overall and for the majority of hESUs (A, D, F, I, M, N; hESUs with fewer than five species were not included). These results concur with inference from the gamma statistic and the MCCR test, which suggested the diversification rate in Pristimantis was decreasing (γ = −5.67, p = 0.0008). Lineage through time plots for Pristimantis also suggest that species accumulation declined towards the present, with a plateau from ca. 11 MYA relative to a pure-birth model of diversification (Fig. S3, see supplemental material online).

BAMM analyses additionally suggested a density-dependent speciation process for the genus Pristimantis, with a speciation rate (mean λ = 0.18 species/My) that started high and declined towards the present (Fig. S4, Fig. S5, see supplemental material online). BAMM detected positive but low (mean μ = 0.03 species/My) extinction rate, but interpretation of extinction rates from trees without fossil data should be taken cautiously. This analysis does not support a shift in diversification rates for Pristimantis (Bayes Factors of 1 to 3 shifts each <1). For the marginal odds ratio, we retained rate shifts with a ratio greater than or equal to 5 (per default in BAMMtools); this gave two scenarios, one with no shifts being by far most likely (Fig. S4, f = 0.9, see supplemental material online) and a much less supported scenario with a single shift in speciation rate in hESU H (Fig. S4, f = 0.065, see supplemental material online). This is the same hESU for which the LASER analysis recovered a different rate-variable model from the other hESUs, i.e. a yule3rate as the best fitting diversification model (Table S3, see supplemental material online). A yule-n-rate model allows for shifts in speciation rates over time; this may explain congruence with BAMM identifying a shift in speciation within Pristimantis only in that clade. The rate-through-time plot for each hESU from the BAMM analysis shows declining speciation rates for Pristimantis and each hESU (A-N) (Fig. S5, see supplemental material online). Speciation rates and declines were similar between hESUs and Pristimantis, indicating support for a zero-shift scenario within the genus. We interpret the results of BAMM cautiously given possible misestimates, especially with incomplete lineage sampling (Meyer & Wiens, 2018), but note that the patterns of diversification rate are overall concordant across the multiple approaches we explored here.

Discussion

Phylogeny of Pristimantis

The time-calibrated phylogeny places the origin of Pristimantis at 28.97 MYA (95% HDP = 21.59 – 37.33 MYA), in the Oligocene epoch. Most of the phylogenetic relationships we recovered here are in general agreement with previous studies. However, because a number of basal nodes have low support (Fig. 1), as in previous studies (Hedges et al., 2008; Mendoza et al., 2015; Pinto-Sánchez et al., 2012), the ordering of larger clades differs amongst published trees due to these splits being unresolved. The phylogeny of Mendoza et al. (2015) described five main clades (A, B, C, D1, D2) with their A, B, C clades generally resolved in our phylogeny, but with some taxa being unstable (e.g. P. caprifer, P. rozei). Their D1 and D2 clades appear the same in both topologies. Here we do not assess the concordance of relationship of each taxon across our full tree and the tree of Mendoza et al. (2015) as such comparisons overemphasize final topologies, which as we have noted in some cases have weak support at deeper nodes. Rather, we focus on the statistically supported groupings at shallower scales inferred by hESUs and the taxonomic implications of the species groupings associated with our new taxa.

Cryptic species and undescribed diversity

In many instances for which multiple exemplars of a species were included in the phylogeny, the “species” was not resolved as monophyletic and/or showed substantial phylogenetic distinctiveness; for example, P. altamazonicus-diadematus, P. croceoinguinis, P. malkini, P. variabilis, P. taeniatus (sensu lato), P. ockendeni, P. librarius, P. martiae (sensu lato). P. zophus, and P. bogotensis-frater. We hypothesise that cryptic diversity or misidentification of specimens may be the explanation in most of these cases. We highlight several incidences of these two factors occurring within our phylogeny in Table S2 (see supplemental material online), along with suggested revisions to aid in inference of the true relationship and to draw attention to these vouchers. These groups require further investigation of additional specimens and molecular information to advance alpha taxonomy.

Misidentification can easily occur as comprehensive resources may not always be available and taxa often have conserved morphology. In Pristimantis, many species are very morphologically similar and there are challenges for identification without intrusive investigations (Hedges et al., 2008) or genetic analysis. Pristimantis has been understudied in terms of cryptic diversity compared with more charismatic tropical amphibians such as dendrobatids and treefrogs (Funk et al., 2012; Lötters, Schmitz, Reichle, Rödder, & Quennet, 2009; Vaz-Silva & Maciel, 2011). Species delimitation is particularly important in highly diverse yet threatened ecosystems such as Neotropical forests, and the current number of described species is undoubtedly a gross underestimation of the true diversity within Pristimantis.

Evolutionary patterns of diversification: Rate and distribution

We report four periods of increased phylogenetic branching since the origin of the Pristimantis. The first increase at ca. 24 MYA may correspond with the start of the Andes formation in the northern regions. The extreme geological activity that characterised this period drove speciation in some South American biota (Hoorn et al., 2010), and it has been previously suggested to explain this first burst of diversification in Pristimantis (Heinicke et al., 2007; Mendoza et al., 2015; Pinto-Sánchez et al., 2012). The other increases in the early to mid-Miocene, ca. 20.5, 16.8 and 15.7 MYA, may correspond to the extensive landscape changes in the era of the Pebas system (Hoorn et al., 2010). This vast wetland across Amazonia is thought to have driven the radiations of several taxa, including molluscs and ostracods (Wesselingh and Salo, 2006), plants (Antonelli et al., 2010), crocodilians and turtles (Hoorn et al., 2010), salamanders (Elmer, Bonett, Wake & Lougheed, 2013), and dendrobatid frogs (Santos et al., 2009), and may have directly or indirectly promoted diversification for several of the Pristimantis phylogenetic hESU clusters we have identified here.

Multiple lines of evidence in our analyses suggest that diversification rates within the genus Pristimantis were faster at the time of origin and slowed down towards the present. The LTT plots approximate the expected “decrease in diversification rates with early radiations” pattern of an upward sloping curve ending with a plateau (e.g. Couvreur, Forest & Baker, 2011). A decrease in diversification rates, in particular density-dependent ones, is a common pattern of evolutionary radiations (De-Silva et al., 2016; Etienne & Haegeman, 2012; Rabosky & Lovette, 2008b). This may represent density-dependent speciation in evolutionary radiations due to complex interaction between species traits and geographical range size, which in turn could facilitate spatial drivers of speciation on continents (Rabosky & Lovette, 2008a). The gamma statistic and the MCCR test we applied in the speciation rates analysis accommodate incomplete taxonomic sampling (Pybus & Harvey, 2000) but we acknowledge that these rate estimates may be affected by unequal and incomplete diversity throughout our tree (Cusimano & Renner, 2010). Our results here, in agreement with previous studies (Mendoza et al. 2015; Pinto-Sánchez et al. 2012), suggest that Pristimantis started diversifying in an era that is associated with Andean uplift and the formation of the Pebas system, which may have provided greater habitat complexity favouring isolation (Hoorn et al. 2010). As Pristimantis species are direct developers, this group may have been able to exploit this new opportunity to colonise and diversify in new areas. Hutter, Lambert & Wiens (2017) looked at the relationship between elevation and species richness in the Andes, using Hyloidae as model system, and found strong association between rates of diversification and elevational change in Pristimantis, suggesting a close link between Andean uplift and Pristimantis diversity.

Evolutionary radiation processes have been suggested for several taxa as a result of Andean uplift (e.g. Ceccarelli et al., 2016; Sanín et al., 2016). We find support for an evolutionary radiation scenario for Pristimantis, with faster speciation rates at the beginning of the radiation and slowdowns towards the present. However, continental evolutionary radiations do not necessarily show these diversification patterns. For example, Schweizer, Hertwig & Seehausen (2014) studied the radiation of Neotropical parrots and did not identify either a rapid diversification rate at the origin of the group nor a decrease in diversification rate towards the present, suggesting a lack of density-dependent speciation across the continent. Furthermore, the diversification rates of phenotypic evolution in neotropical ovenbirds (family: Furnariidae) suggested a different evolutionary history compared to the speciation rates, suggesting that morphological differentiation and speciation processes may not be tightly linked (Derryberry et al., 2011). Significant correlations between some phenotypic traits, such as skin morphology, and species richness were found for new world direct developing frogs (Gonzalez-Voyer et al., 2011). Further studies on differentiation in morphological traits within and between Pristimantis hESUs would be valuable for investigating whether the patterns in speciation and phylogenetic diversification processes correspond to morphological evolution in this group.

Species groups and diversity in Pristimantis

Species are traditionally considered as the fundamental units of evolution, with higher taxonomic levels (e.g. genera and families) functioning as more arbitrary categories in classifying the diversity of life (Barraclough & Humphreys, 2015). However, the complex interplay between microevolutionary and macroevolutionary diversity suggests that fine-scale processes within and among closely related species contribute to higher-level patterns (Humphreys & Barraclough 2014). For example, geographical isolation and ecological divergence promoting species diversity can also lead to the formation of discrete phylogenetic clusters such as hESUs. These hESUs are groups of taxa that share a phylogenetic history, phylogenetic rates and patterns of diversification. Our analysis recovered between 13 and 19 phylogenetic clusters within our broad phylogeny of Pristimantis, with the most evidence for the 14 hESUs we explored in detail (Fig. 2, Table 1). We hypothesise that significant phylogenetic clusters that have undergone shared evolutionary histories and rates underlie the high extant diversity of Pristimantis. Simulations and empirical analyses suggest that incomplete sampling can either increase, decrease or not affect the inferred number of hESUs, dependant on the phylogenetic history and diversity of the specific clades (Humphreys and Barraclough, 2014); well-supported hESUs are generally less sensitive to these issues (Humphreys and Barraclough, 2014). Our sampling of Pristimantis is far from complete (approx. 30% of species), which may ultimately revise the number or complement of hESUs. Nonetheless, given the information at hand, we consider that these clusters represent evolutionarily informative phylogenetic species groupings, which would be a more accurate biological representation of taxon similarities than the phenetic species groupings established in the pre-molecular era.

The hESUs we report here for the first time do not relate to any previously described species group, apart from hESU A, which corresponds to the peruvianus species complex (as per Hedges et al., 2008). The issue of non-monophyly among current species groups within Pristimantis has been raised by many authors (e.g. Hedges et al., 2008; Mendoza et al., 2015; Padial et al., 2014; Pinto-Sánchez et al., 2012), with Heinicke et al. (2015) defining this problem as “unwieldy”. After reclassifying the species groups based on further molecular work, Padial et al. (2014) remarked that there is a lack of morphological synapomorphies for these groups, which could be a cause for concern given the increasing number of Pristimantis being described and great uncertainty around their placement. Thus, a critical evaluation of morphological characteristics within robustly statistically supported molecular groupings, such as our hESUs proposed here, is an important next step to clarify the many unresolved species groups, cryptic species, and evolutionary relationships in Pristimantis. Given the extremely high level of diversity in Pristimantis, combining phylogenetic history (as done here), biogeography/ecology, and alpha taxonomy will be critical to advancing our understanding of systematics for this group and evolutionary processes more generally. Our findings suggest that efforts at revising groups of taxa in Pristimantis would benefit from a comprehensive assessment of specimens and species within each hESU, including their genetic, ecological, morphological and geographic patterns, rather than based on any other grouping. Our study and many others show that approaches based on deep exploration of single localities, non-monophyletic species groups, or alpha taxonomy based on morphology (e.g. Heinicke et al. 2015), will provide important but scattered information across the genus Pristimantis. We propose that focusing on hESUs to work toward establishing new phylogenetic species groups would be an efficient way to inform both evolutionary process and taxonomic outcome simultaneously.

The number of species within the hESUs varied from three to 22. These are likely under sampled and do not necessarily reflect the full species complement in each group. There are representative samples from eleven countries in this dataset, but most samples come from just two countries (Colombia and Ecuador), with almost half of all specimens (n = 122) coming from Ecuador, including all the new DNA sequence data. Records of Pristimantis on the Global Biodiversity Information Facility (GBIF) show the same pattern, with 80% of all records from Ecuador and Colombia (accessed June 2017). These two countries have the highest species richness of the genus Pristimantis (Colombia = 220, Ecuador = 211 species; Frost, 2018); however, other South American countries also have high species richness but many fewer records. For example, Peru has 141 known species (Frost, 2018) but only 7% of the records on GBIF; in addition, Brazil has both low numbers of records on GBIF (4%) and low species richness (n = 30). More sampling is needed from across the full distribution of Pristimantis to capture the true diversity of this group.

Conclusions and further implications

We inferred a multi-gene phylogeny of the species-rich leaf litter frog genus Pristimantis, adding 50 new exemplars (representing eight known species, nine known but problematic [or unresolved] taxa, and two potentially new species) from Ecuador for a total for 260 individuals. Rates of diversification were faster early in Pristimantis, which had an origin ca. 28 MYA, and slowed approaching present times, which is a pattern usually associated with evolutionary radiations. Our study suggests high levels of cryptic species diversity and a lack of consistency in identification across a number of wide-ranging species. In an effort to better resolve taxonomic and evolutionary groupings, we identified 14 hESU clusters that have experienced independent and shared evolutionary histories, which could be interpreted as new monophyletic and evolutionarily informed species groups. To better understand the evolutionary history of this group, molecular studies should be conducted on more species, so that phylogenetic relationships, hESUs delimitation, and diversification rates are more clearly resolved. Further phylogenetic and population genetics research on ecological contextualised and co-distributed lineages of Neotropical amphibians is key to resolving the primary historical and contemporary drivers of diversification.

Supplemental data

Supplemental data for this article can be accessed at: http://dx.doi.org/10.1080/14772000.2018.1503202.

Acknowledgements

We thank L. Coloma for facilitating sampling, identifying specimens, and support in the QCAZ museum; G. Vigle, R.O. de Sá, M. Bustamante, Z. Aguilar, S. Ron, S. Padilla, JC Santos, C. Proaño, and N. Acosta for specimens contributed to QCAZ; and F. Ayala-Varela, T. Sugahara, P. Menéndez, S. Padilla, and I. Tapia for assistance to KRE in the field. This project was funded by: NSERC PGS, Queen’s Field Travel grant, NERC NE/N003942/1, and University of Glasgow support to KRE, NSERC Discovery grant to SCL, and Amphibian Tree of Life with DCC. We thank A. Crawford and an anonymous reviewer for comments that improved the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.