Among direct-developing rain frogs of the genus Craugastor is a clade of 19 described species (bocourti series) that occur in Mexico and northern Central America. Many of these 19 species have been described based on subtle morphological differences and have never been examined using molecular data. Here, we used a multilocus dataset (one mitochondrial (mtDNA) and four nuclear (nDNA) gene fragments, totalling 3,048 concatenated base pairs) to investigate species limits and phylogenetic relationships among 60 northern rain frogs referable to 12 species, with a focus on species from Guatemala. We inferred phylogenies using maximum likelihood and Bayesian analyses on separate mtDNA and nDNA datasets. Concatenated and coalescent species-tree analyses support the monophyly of multiple species, with interspecific relationships mostly unresolved. These mtDNA and nDNA trees were often incongruent with morphology-based taxonomy. For example, two genetically shallow clades contained individuals referable to at least five described species, whereas a single described species contained deep divergences indicative of multiple cryptic species. These findings indicate that morphology-based taxonomy has both overestimated and underestimated actual species diversity (depending on the species), an interpretation supported by two molecular species-delimitation procedures. Based on these findings, we synonymise C. glaucus (Lynch, 1967) and C. stuarti (Lynch, 1967) with C. xucanebi (Stuart, 1941). We also synonymise C. nefrens (Smith, 2005) and C. cyanochthebius McCranie & Smith, 2006 with C. campbelli (Smith, 2005). The molecular data also support multiple undescribed species, notably within C. decoratus (Taylor, 1942). Overall, we show how morphology-based species delimitation can both underestimate and overestimate species richness in morphologically conservative groups.
Introduction
Using DNA sequence variation to infer species boundaries has greatly increased our understanding of amphibian evolution. For many anurans (frogs and toads), these data have revealed vast underestimates of species richness (Stuart et al. 2006, Fouquet et al. 2007, Vieites et al. 2009, Funk et al. 2012, Rojas et al. 2018, Scherz et al. 2022). This underestimation has been particularly well-documented in a large radiation of New World direct-developing frogs called Brachycephaloidea (Padial et al. 2014) or Terraranae (Heinicke et al. 2018). In this clade, molecular data have revealed high levels of morphologically cryptic species diversity (Miyamoto 1983, Crawford 2003, Crawford et al. 2007, Padial & De La Riva 2009, Kieswetter & Schneider 2013, Fusinatto et al. 2013, Pie et al. 2018). Interestingly, the discovery of ‘cryptic’ taxa via molecular methods even extends to the familial level in this group, with many species traditionally placed in the genus Eleutherodactylus now distributed across multiple families (Heinicke et al. 2009). While the increasing affordability and ease of collecting DNA sequence data has shed much light on terraranan diversity, extensive molecular studies have yet to be conducted on many groups.
Northern rain frogs (Craugastor bocourti Species Series; sensu Hedges et al. 2008, hereafter ‘bocourti series’) occur as far north as southern Tamaulipas in Mexico, as far west as central Guerrero in Mexico, as far south as Guatemala, and as far east as Honduras (Fig. 1). Most members of the bocourti series were previously included in the Eleutherodactylus alfredi group (Lynch 1966, 1967a, Campbell et al. 1989, Smith 2005). These frogs now belong to the family Craugastoridae (sensu Barrientos et al. 2021) and the subgenus Hylactophryne (sensu Hedges et al. 2008, Padial et al. 2014, Heinicke et al. 2018). Morphologically, species in the bocourti series are easily distinguished from congeners by the presence of large, expanded pads on fingers III and IV (Campbell et al. 1989; Fig. 2). The bocourti series currently contains 19 described species: C. alfredi, C. batrachylus, C. bocourti, C. campbelli, C. cyanocthebius, C. decoratus, C. galacticorhinus, C. glaucus, C. guerreroensis, C. megalotympanum, C. nefrens, C. polymniae, C. silvicola, C. spatulatus, C. stuarti, C. taylori, C. uno, C. xucanebi and C. yucatanensis (Frost 2023, Hedges et al. 2008). Conducting research on the bocourti series is challenging because most species are rarely encountered in the field (Streicher et al. 2011, Palacios-Aguilar 2017, Carbajal-Márquez et al. 2019). For example, 74% of these 19 species were described based on three or fewer specimens (Boulenger 1898, Taylor 1940, 1942, Stuart 1941, Shannon & Werler 1955, Lynch 1966, 1967a, b, Savage 1985 ‘1984’, Campbell et al. 1989, Canseco-Márquez & Smith 2004). Furthermore, obtaining tissue samples for molecular analysis can be difficult because some species are considered vulnerable by the IUCN (e.g. Craugastor uno, Angulo et al. 2020; Craugastor alfredi, Luría-Manzano & Ramirírez-Bautista 2017). Despite these limitations, previous molecular studies have included up to ten described species from the bocourti series and supported its monophyly (Crawford & Smith 2005, Heinicke et al. 2007, Hedges et al. 2008, Pyron & Wiens 2011, Padial et al. 2014, Portik et al. 2023).
In addition to small sample sizes, another serious issue in the systematics of the bocourti series is that many species are distinguished by only subtle morphological differences. These differences include traits like palmar and plantar tubercle depth (Smith 2005), colour pattern (Canseco-Márquez & Smith 2004), and the relative widths of finger and toe pads (Shannon & Werler 1955, Lynch 1967a, McCranie & Smith 2006). Because of this overall morphological uniformity, assigning individuals to species, particularly juveniles and poorly preserved natural history specimens, is often a nontrivial task. This issue is exacerbated in locations where poorly known taxa are putatively sympatric (e.g. Central Mexico and southern Guatemala; Fig. 1).
In this study, we build on previous work by generating a multilocus dataset (one mitochondrial and four nuclear gene fragments) for 12 described species of the bocourti series with a focus on species from Guatemala and Honduras (Fig. 1). We used these data to produce gene and species trees and estimate species limits to compare to morphology-inferred species boundaries.
Material and Methods
Assigning specimens to species
As noted above, species in the bocourti series are challenging to identify because morphological differences among them are subtle, and many poorly known taxa occur in sympatry in some regions (e.g. Central Mexico and southern Guatemala; Fig. 1). Because we sampled several specimens from these challenging regions, we used comparisons of voucher specimens with type specimens to assign specimens to species. In cases where we could not confidently match vouchered material to types, we describe specimens using open nomenclature (sensu Bengston 1988). Specifically, we used the designator affinis (aff.) when a specimen had a morphology that was highly similar to (but not completely congruent with) a described species.
Fig. 1.
Geographic distribution of samples used in the molecular analysis (circles) and type localities (stars) for the Craugastor bocourti series. Filled stars correspond to species with topotypic sampling.

We were able to confidently match specimens used for our molecular analyses to the following species based on type comparisons (type specimens used for comparisons in brackets; see acknowledgements for definitions of museum abbreviations): C. alfredi (BMNH 1947.2.15.54-55), C. bocourti (MNHNP 6413), C. campbelli (UTA A-33452), C. glaucus (TCWC 21463), C. nefrens (UTA A-45279), C. stuarti (UMMZ 126738), C. uno (UTA A-7984), and C. xucanebi (UMMZ 89914). We applied open nomenclature to several specimens originating from the Sierra Madre Oriental and Sierra Madre del Sur of Mexico (C. aff. decoratus, C. aff. polymniae, C. aff. silvicola, and C. aff. spatulatus) and eastern Honduras (C. aff. nefrens) because they did not convincingly match type specimens. These assignments resulted in 12 species names being associated with our molecular sampling, either through confident identification based on type specimens (eight species) or substantial morphological similarity to the type material (four species). Two of these species, C. glaucus and C. nefrens, have not previously had molecular data generated for them.
Taxonomic, genetic, and geographic sampling
We collected specimens from the field in Mexico, Guatemala, and Honduras and generated molecular data for 49 individuals assigned to eight morphology-based species and for 11 individuals with uncertain taxonomic affinities (open nomenclature). The sampling included specimens of C. alfredi (n = 4), C. bocourti (n = 3), C. campbelli (n = 3), C. aff. nefrens (n = 4), C. aff. decoratus (n = 7), C. glaucus (n = 1), C. nefrens (n = 1), C. aff. polymniae (n = 2), C. aff. silvicola (n = 1), C. aff. spatulatus (n = 2), C. stuarti (n = 3), C. uno (n = 4), and C. xucanebi (n = 25). Specific collection localities for all specimens are listed in Table 1. We obtained topotypic (or nearly topotypic) material for seven of the eight confidently identified species (C. alfredi, C. bocourti, C. campbelli, C. glaucus, C. nefrens, C. stuarti, and C. uno; Fig. 1).
We selected molecular markers used extensively in terraranan phylogenetics (e.g. Heinicke et al. 2007, Padial et al. 2014). For all 60 individuals, we sequenced a 460 base pair (bp) segment of the mitochondrial 12S ribosomal subunit gene (12S). We obtained a multilocus nuclear DNA (nDNA) dataset for a subset of these individuals (n = 50), including representatives of all eight confidently identified species. Our nDNA dataset consisted of fragments of four genes: 351 bp of rhodopsin (exon 1; Rho), 653 bp of the recombination-activating protein 1 (RAG-1), 560 bp of tyrosinase precursor (Tyr), and 981 bp of the cellular-myelocitomatosis gene (c-myc). The c-myc gene fragment consisted of separate 403 bp and 578 bp sections that correspond to exons two and three, respectively. The outgroup taxa sampled included: 1) the sister group of the bocourti series (the C. augusti series; Hedges et al. 2008, Streicher et al. 2014), including one individual of C. augusti and one of C. tarahumaraensis; and 2) the species C. daryi, a representative of the Craugastor subgenus Campbellius (Hedges et al. 2008). The subgenus Campbellius was targeted as a distant outgroup because it has been supported as the earliest branching lineage of Craugastor in multiple studies (Crawford & Smith 2005, Hedges et al. 2008, Pyron & Wiens 2011, Padial et al. 2014, Portik et al. 2023). A complete list of voucher specimens, locality information, and GenBank numbers are provided in Table 1.
Fig. 2.
Examples of the Craugastor bocourti series. Clockwise from top left, palmar surface of the holotype of Eleutherodactylus hidalgoensis (FMNH 100094) demonstrating the expanded pads on fingers III and IV that is characteristic of the series; C. uno from Guerrero, Mexico; C. aff. decoratus from San Luis Potosí, Mexico; C. aff. decoratus from Hidalgo, Mexico; C. aff. nefrens from Cortés, Honduras; C. campbelli from Izabal, Guatemala.

Laboratory methods and DNA sequencing
We obtained genomic DNA from muscle and liver tissues using Qiagen® kits. Before DNA extraction, tissue samples were stored in either 70-100% ethanol or an SDS-based lysis buffer. We used Promega reagents (e.g. polymerases, dNTPs, MgCl2) to amplify DNA fragments. PCR amplification was performed using Green Taq Master Mix (Promega) and previously designed oligonucleotide primers (Table 2). To amplify 12S we used a standard thermal cycling profile with an annealing temperature of 50 °C. We amplified nuclear loci (c-myc, RAG-1, Rho, Tyr) using the touchdown thermal cycling protocol described in Streicher et al. (2009). Cycle sequencing was performed using Big-Dye terminator chemistry, and samples were sequenced by either SeqWright (Houston, Texas, USA) or the University of Texas at Arlington (UTA) genomics core facility (Arlington, Texas, USA). We used the PCR primers listed in Table 2 as sequencing primers and sequenced gene fragments in both directions for all loci except the c-myc exons, which were only sequenced in one direction.
We used the programs Sequencher v 4.1 (Gene Codes Corp; Ann Arbor, Michigan, USA) or Geneious v 7.1.3 (Biomatters; http://www.geneious.com) to clean the resulting DNA sequences. For nDNA loci, we assumed that a site was heterozygous (coded as a degenerate base) if equal chromatogram peaks were present for both bases (Hare & Palumbi 1999). We did not phase nuclear loci because putatively heterozygous sites were rare (< 10 observed) and thus unlikely to have a large influence on topological reconstructions. Instead, we coded these sites as degenerate bases for analysis and GenBank submission.
Table 1.
Specimen information and GenBank accession numbers for samples of Craugastor used in this study. MX = Mexico, GT = Guatemala, and HD = Honduras.

Continued

Continued

Continued

Continued

Table 2.
Primers used to amplify and sequence Craugastor mitochondrial (mtDNA) and nuclear (nDNA) gene fragments.

Sequence alignments were performed using MUSCLE (Edgar 2004) software executed in MEGA 5.1 (Tamura et al. 2011). We identified reading frames for protein-coding genes (RAG-1, Tyr, Rho, c-myc) and performed multilocus concatenation using Geneious v 7.1.3. All alignments used in this study were deposited in the Natural History Museum Data Portal ( http://data.nhm.ac.uk).
Separate phylogenetic analyses of mitochondrial and nuclear DNA
We separately analysed mtDNA and concatenated nDNA data using maximum likelihood (ML) and Bayesian phylogenetic methods. To select an appropriate model of nucleotide evolution, we used MEGA 5.1 for mtDNA (12S), treating the entire gene fragment as a single partition. For nDNA, we used PartitionFinder v1.1 (Lanfear et al. 2012) with codon partitions specified for all loci (including for exons 1 and 2 of c-myc separately). The nucleotide models selected for each partition are listed in Table 3.
We performed ML analyses using MEGA 5 using a heuristic tree search criterion with nearest-neighbour interchange. We assessed ML branch support using bootstrapping analysis (1,000 pseudoreplicates). We performed Bayesian Markov Chain Monte Carlo (MCMC) analyses using MrBayes 3.2.1 (Ronquist & Huelsenbeck 2003). For Bayesian MCMC analyses we unlinked all parameters between the identified partitions and used a variable rate prior. We used default settings for Bayesian MCMC concatenated analyses except that we increased the sampling frequencyto1,000generationsandrantheanalysesfor 10 million generations. We used the online software ‘Are we there yet?’ (AWTY; Wilgenbusch et al. 2004) to select an appropriate number of sampled trees to discard as burn-in before summarising the paired Bayesian MCMC runs and calculating posterior probabilities, using the majority rule consensus of the post-burnin trees. We visualised trees using the software FigTree v 1.4 (available at http://tree.bio.ed.ac.uk/software/figtree/).
Combined Bayesian species-tree inference
We performed a Bayesian species-tree analysis using all five markers (1 mtDNA + 4 nDNA) with the software *BEAST v 1.6.0 (Heled & Drummond 2010) and BEAST 1.7.5 (Drummond & Rambaut 2007). As in our MrBayes 3.2.1 concatenated analysis of nDNA, we treated the two c-myc exons as separate data partitions. While recombination within c-myc is unlikely, we chose to treat the exons as separate partitions because model selection identified different codon model schemes for each exon (Table 3). We used the seven matched mtDNA and nDNA clusters of individuals as our a priori ‘species’ for the *BEAST analysis (Table 4). Clusters were identified using the following criteria: 1) they contained the same individuals, 2) they were monophyletic in both nDNA and mtDNA analyses, and 3) their monophyly received high Bayesian statistical support. We also included two individuals of C. augusti and one individual of C. tarahumaraensis as outgroups.
Table 3.
Models of nucleotide evolution selected for alignments of Craugastor mitochondrial (mtDNA) and nuclear (nDNA) gene fragments. Partition column indicates regions of nDNA that were included in the same partition because of Model overlap.

We used the same model-partitioning scheme and models employed for the separate analyses of mtDNA and nDNA when using MrBayes (Table 3). We used a Yule tree prior and unlinked loci and substitution models across the dataset. We used an uncorrelated relaxed-clock model to account for rate variation among branches. Species-trees were selected and scored based on a 10 million generation MCMC analysis sampled every 1,000 generations. We used Tracer v 1.5 (Rambaut & Drummond 2007) to ensure that all parameters had reached effective sample sizes over 200. We used a burn-in of 10% (1 million generations) before summarising species trees to determine the 95% Highest Posterior Density (HPD) of posterior distributions. To visualise the 95% HPD tree set, we used the program DensiTree v 2.1 (Bouckaert 2010).
Table 4.
Matched mitochondrial (mtDNA) and nuclear (nDNA) DNA clusters of Craugastor used in species tree analysis.

Species delimitation procedures
To infer species boundaries in each dataset, we used two approaches: ML Poisson Tree Process (PTP) species delimitation (Zhang et al. 2013) and a generalised mixed Yule-coalescent (GMYC) model (Fujisawa & Barraclough 2013). We ran PTP and GMYC separately on the mtDNA ML tree and the combined nDNA ML tree (both rooted with C. daryi) using the online species delimitation server ( https://species.h-its.org/). This server implements a Bayesian version of PTP where Bayesian support values are added to the best species-partitioning schemes to quantify the relative support for each inferred species delimitation. We used the following parameter settings to generate these support values: a search across 100,000 generations, a thinning value of 100, and burn-in of 10%. We ran a single-threshold model of GMYC. To time-calibrate each phylogeny for the GMYC analysis, we used treePL (Smith & O'Meara 2012) and two calibrations from Portik et al. (2023); 1) the divergence between the subgenera Campbellius and Hylactophryne (estimated at 39.1 million years ago, mya) and 2) the divergence between the C. augusti series and C. bocourti series (estimated at 20.3 mya). For the treePL analyses, we used the thorough setting with a smoothing value of 0.1, as determined by cross-validation analyses run on both trees.
Fig. 3.
Mitochondrial phylogram for the Craugastor bocourti series resulting from Maximum Likelihood analysis. Nodal support values from ML bootstrapping appear below branches (grey) and Bayesian inference appears above branches (black). NS indicates branches with no statistical support.

Fig. 4.
Nuclear phylogram for the Craugastor bocourti series resulting from concatenated Maximum Likelihood analysis of four nuclear fragments. Nodal support values from ML bootstrapping appear below branches (grey) and Bayesian inference appear above branches (black). NS indicates branches with no statistical support.

Results
Mitochondrial DNA phylogeny
Using the Bayesian Inference Criterion (BIC) in MEGA 5.1, the GTR+Γ model was selected as the most appropriate model of nucleotide evolution for the mtDNA dataset (Table 3). The AWTY examination of topological congruence (using both symmetrical difference and agreement scores) suggested that paired MrBayes 3.2.1 searches in the mtDNA analysis either did not converge (symmetrical-difference scores) or converged after 5 million generations (agreement scores). Thus, to be conservative we discarded 50% of the retained trees as burn-in (leaving 5,000 trees in the final sample). As in previous studies (Crawford & Smith 2005, Hedges et al. 2008), our Bayesian analysis strongly supported the monophyly of the bocourti series and the subgenus Hylactophryne (Fig. 3). Although deep and shallow nodes in the Bayesian consensus tree were well supported (posterior probability > 0.95), nodes at intermediate depths were often weakly supported. Two exceptions to this pattern were the strongly supported nodes uniting 1) C. uno with C. aff. spatulatus and 2) C. aff. decoratus with C. aff. polymniae. Hereafter, we refer to the clade of C. aff. decoratus + C. aff. polymniae as the C. decoratus complex (Fig. 3). Interestingly, we did not recover a monophyletic C. decoratus complex in our ML analysis. However, the ML analysis did show moderate support for the clade of C. uno + C. aff. spatulatus (bootstrap support (BS) = 78, Fig. 3).
We found that two genetically shallow clades contained all individuals of five putative species: 1) a clade containing C. campbelli, C. aff. nefrens, C. nefrens, and some sampled populations previously identified as C. xucanebi (from the Guatemalan department of Izabal) and 2) a clade containing C. glaucus, C. stuarti, and C. xucanebi (from central and western Guatemala). Hereafter, we refer to the first clade as the C. campbelli complex and the second clade as the C. xucanebi complex. We selected names for these complexes using the oldest available name in the clade, which also included samples from the type locality. Although individuals of C. xucanebi are present in the C. campbelli complex clade, those individuals are from Izabal, Guatemala, far away from the type locality of C. xucanebi in Alta Verapaz, Guatemala (Fig. 1). For the C. xucanebi complex, several samples of C. xucanebi are from Alta Verapaz, Guatemala.
Fig. 5.
Coalescent species tree estimated using *BEAST. Posterior probabilities appear below branches. Dark lines correspond to the consensus tree, and light blue trees in the background are the posterior distribution of species trees, all visualised using DensiTree. Only posterior probabilities > 0.90 are indicated on the consensus tree.

Importantly, these results render the described species C. xucanebi as non-monophyletic. The C. campbelli complex and C. xucanebi complex were found to be sister taxa (although with low support values, BS < 70). We also observed deeper levels of divergence within the C. decoratus complex than expected, with C. aff. polyminae rendering C. aff. decoratus non-monophyletic. The mitochondrial tree supported the monophyly of the sampled individuals of the species C. uno, C. aff. spatulatus, C. alfredi, and C. bocourti.
Nuclear DNA phylogeny
The best-fitting PartitionFinder scheme (lnL = –6749.92; BIC score = 14464) for the concatenated nuclear genes had five partitions with four models (Table 3). In total our concatenated nDNA alignment contained 2,545 bp. The AWTY examination of topological congruence (using both symmetrical difference and agreement scores) suggested that the paired Bayesian nDNA runs converged between 3 and 4 million generations. Therefore, we discarded 50% of the trees (first 5 million generations) as burn-in. As in the mtDNA analysis, our nDNA Bayesian analysis strongly supported the monophyly of the bocourti series and of the subgenus Hylactophryne (Fig. 4). However, similar to the mtDNA dataset, ML bootstrap analysis did not show high support for these groupings. ML bootstrap analysis did show strong support for two clades: 1) the sister relationship between C. uno and C. aff. spatulatus, and 2) the monophyly of the C. decoratus complex. Two additional patterns in the nDNA tree were similar to the mtDNA tree: 1) individuals of the C. campbelli and C. xucanebi complexes were placed within two shallow clades that each contained multiple species and 2) nodes at intermediate depths (i.e. nodes between species) were weakly supported. While there were fewer individuals of C. aff. decoratus included in the nDNA phylogeny, deep levels of divergence among the three individuals were evident, as in the mtDNA topology.
Fig. 6.
Results of Poisson Tree Process (PTP) and generalised mixed Yule coalescent (GMYC) species delimitation indicating that the mitochondrial (mtDNA) dataset supported 18-22 species. Boxes adjacent to tip labels reference the delimitation scheme described in Table 5.

Fig. 7.
Results of Poisson Tree Process (PTP) and generalised mixed Yule coalescent (GMYC) species delimitation indicating that the concatenated nuclear dataset (nDNA) supported 12-14 species. Boxes adjacent to tip labels reference the delimitation scheme described in Table 5.

Table 5.
Results of Poisson Tree Process (PTP) species delimitation among different species of the Craugastor bocourti series and three outgroup species. Support for different species delimitations is in Bayesian support (BS) as determined by PTP. Sample sizes (n) are the number of individuals included in each putative species. GMYC column described if Generalised mixed Yule coalescent delimitation detected the same entity as PTP. Results in bold indicate those species delimited consistently in all delimitation analyses of mitochondrial (mtDNA) and nuclear (nDNA). See text for definitions of C. campbelli complex and C. xucanebi complex.

Species-tree analysis
The 95% HPD set of species trees (Fig. 5) was congruent with the separate mtDNA and concatenated nDNA analyses (Figs. 3 and 4) in that the only relationships supported by most of the trees (i.e. posterior probability > 0.90) were the clade uniting C. uno and C. aff. spatulatus and the monophyly of the bocourti series relative to the outgroups (C. augusti and C. tarahumaraensis).
Species-delimitation analysis
The PTP delimitation procedure inferred 18 species in the mtDNA dataset and 14 species in the nDNA dataset (Table 5). The GMYC delimitation procedure inferred 22 species in the mtDNA dataset (confidence interval, CI; 16-24) and 12 species (CI; 10-14) in the nDNA dataset (Figs. 6 and 7). Comparisons between the delimitation procedures revealed many similarities. First, all delimitation procedures supported recognising the C. campbelli complex as a single species (Figs. 6 and 7, Table 5). The C. xucanebi complex was recognised as a single species in all but the GMYC mtDNA delimitation (Fig. 6; see explanation below). In PTP analyses of these two complexes, Bayesian support for a single species within each complex was low in the mtDNA delimitation but higher in the nDNA delimitation (PTP; C. xucanebi complex, mtDNA = 0.09, nDNA = 0.42; C. campbelli complex, mtDNA = 0.16, nDNA = 0.78). Second, all delimitation procedures supported the recognition of each of the following taxa as single species: C. bocourti (PTP; mtDNA = 0.80, nDNA = 0.68), C. uno (PTP; mtDNA = 0.47, nDNA = 0.39), and C. aff. spatulatus (PTP; mtDNA = 0.47, nDNA = 0.49). Finally, all delimitations procedures recognised the three outgroup taxa as species: C. daryi (PTP; mtDNA = 1.00, nDNA = 1.00), C. augusti (PTP; mtDNA = 0.96, nDNA = 0.97), C. tarahumaraensis (PTP; mtDNA = 0.96, nDNA = 0.97).
There were also differences between PTP and GMYC species delimitation results (Table 5). The PTP and GMYC mtDNA delimitations were similar, except for the C. xucanebi complex, where GMYC inferred five species. A similar result was observed for the PTP and GMYC nDNA delimitations, where GMYC inferred a single species in C. alfredi, whereas PTP inferred three different species within the C. alfredi clade. A related and difficult-to-interpret result is that two of the three individuals of C. alfredi in the PTP delimitation with nDNA had zero (or close to zero) differences in branch length. We would not expect two individuals with such shallow divergence to be different species (Fig. 7).
Discussion
Analyses of mtDNA and nDNA revealed four noteworthy patterns: 1) deep levels of divergence within the C. decoratus complex, 2) low levels of genetic divergence among named species within the C. campbelli complex and C. xucanebi complex, 3) evidence for the distinctiveness of three species (C. alfredi, C. bocourti, and C. uno), and 4) an inability to strongly resolve interspecific relationships near the base of the bocourti series. Importantly, the first two patterns, coupled with the results of the species delimitation analyses, indicate that morphology-based taxonomy has both over-delimited and under-delimited species in this group, depending on the taxon. Below we discuss our results and their implications for the taxonomy and future study of northern rain frogs.
Craugastor decoratus complex: under-delimited
Taylor (1942) described C. decoratus (as Eleutherodactylus decoratus) from near ‘Banderia’ (Banderilla) in the Mexican state of Veracruz, and Eleutherodactylus hidalgoensis from Tianquistengo, Hidalgo. Lynch (1967c) described latitudinal variation in C. decoratus, synonymised E. hidalgoensis with C. decoratus, and designated two subspecies, C. d. decoratus (from Hidalgo and Veracruz) and C. d. purpurus (from Gomez Farias, Tamaulipas, México). These were originally described as E. d. decoratus and E. d. purpurus, respectively. Campbell et al. (1989) described C. polymniae (as Eleutherodactylus polymniae) from the Sierra Juárez near Vista Hermosa in the Mexican state of Oaxaca. Bayesian MCMC (MrBayes 3.2.1) analysis of mtDNA recovered three well-supported clades within our C. decoratus complex (D1-D3; Fig. 3). While clade D2 contains individuals that we assigned to C. aff. polymniae, clades D1 and D3 (which are not sister clades in our mtDNA analysis) contain specimens that are morphologically similar to the C. decoratus material examined by Lynch (1967c), indicating that there are likely multiple species within the current concept of C. decoratus. Indeed, the species delimitation results for this complex support this interpretation (Table 5). Notably, clade D3 contains an individual from approximately 6 km from the type locality of E. hidalgoensis, suggesting that this name may be available for the inferred species.
Our nDNA dataset lacked many individuals in the mtDNA dataset, so future nDNA sequencing will provide valuable information for delimiting candidate species. Although we lacked genetic samples to include it here, we suspect that C. batrachylus is also a member of the C. decoratus complex given the geographic proximity of the type locality in southern Tamaulipas to the type locality of C. decoratus purpurus and populations of C. decoratus in Hidalgo and San Luis Potosí (Fig. 1). This possible relationship is also suggested by the overall morphological similarity between C. batrachylus and C. decoratus (see Taylor 1940).
Craugastor campbelli complex: over-delimited
Smith (2005) described C. campbelli and C. nefrens (as Eleutherodactylus campbelli and E. nefrens) from Guatemala in the montane wet forests of the Montañas del Mico and Sierra de Caral, respectively. He also examined specimens of C. xucanebi from the Sierra de Santa Cruz in the western Izabal department and concluded that they differed from C. campbelli and C. nefrens. Our results suggest that all northern rain frogs from southern Izabal (south of Lake Izabal and the Polochic and Dulce rivers), adjacent eastern Honduras, and some populations in the Sierra de Santa Cruz, form a clade with shallow within-clade divergence, not including C. xucanebi from central and western Guatemala (also inhabiting the Sierra de Santa Cruz) and the adjacent state of Chiapas in Mexico (Figs. 3 and 4). The occurrence of geographic variation in the morphology of young northern rain frog lineages is not unprecedented (Streicher et al. 2011) and is a likely contributor for the incongruence between morphological and molecular variation in the C. campbelli complex.
The C. campbelli complex occurs in the Guatemalan department of Izabal across the Motagua-Polochic fault system, an important biogeographic boundary for many terrestrial vertebrates (Castoe et al. 2009, Daza et al. 2010). However, much like a study on Bolitoglossa salamanders (Rovito et al. 2012), we found that the C. campbelli complex ranges across the faults, indicating either that 1) these features have not limited dispersal or 2) recent range expansions have obscured the historical importance of this barrier.
Based on our phylogenetic results (Figs. 3 and 4) and species delimitation results (Figs. 6 and 7), we recommend that C. nefrens be considered a junior synonym of C. campbelli (the first species described in Smith (2005)). We also recommend that some populations of C. xucanebi inhabiting the Sierra de Santa Cruz and Puerto Barrios regions of Izabal be referred to as C. campbelli, those with an expanded toe V tip. We also consider the sampled individuals of C. aff. nefrens (frogs from eastern Honduras) to be C. campbelli. McCranie & Smith (2006) reported that C. cyanocthebius from western Honduras was closely allied to C. campbelli and C. nefrens. Although we lack molecular data from the type locality of C. cyanocthebius (Fig. 1), the synonymisation of C. nefrens with C. campbelli renders C. cyanocthebius morphologically undiagnosable as there is now overlap in the total number of vomerine teeth (0–8 in C. cyanochthebius versus 0–14 in C. campbelli), the shape of the proximal subarticular tubercle of Toe V (rounded in C. cyanochthebius versus pointed to rounded in C. campbelli), and the ratio between proximal and distal subarticular tubercles of Toe V (see McCranie & Smith 2006). As such, we recommend that C. cyanocthebius also be considered a junior synonym of C. campbelli.
Craugastor xucanebi complex: over-delimited
Stuart (1941) described C. xucanebi(as Eleutherodactylus xucanebi) from the department of Alta Verapaz in Guatemala. Lynch (1967a) described C. glaucus (as E. glaucus) and C. stuarti (as E. stuarti) from the state of Chiapas in Mexico and the department of Huehuetenango in Guatemala, respectively. Lynch (1967a) compared these two species with other Craugastor species, but he did not compare C. glaucus and C. stuarti to one another. Furthermore, both species described by Lynch (1967a) exclusively used specimens collected by L.C. Stuart and the descriptions relied heavily on colour pattern. Based on our phylogenetic results (Figs. 3 and 4) and three out of four species delimitation results (Figs. 6 and 7), we recommend that C. glaucus and C. stuarti be considered junior synonyms of C. xucanebi. The mtDNA GMYC species delimitation analysis inferred five species within the C. xucanebi complex, including a species containing a single individual of C. glaucus and a species containing a clade of mostly C. stuarti (Fig. 6). Nevertheless, we consider this delimitation scheme to be ‘over-split’ given 1) the shallow genetic divergence levels within the C. xucanebi clade and 2) the results of the three other molecular species delimitation analyses which all inferred a single species for this clade.
We suspect that C. taylori, although not included in our molecular analysis, may be closely related to (or a member of) the C. xucanebi complex, given its geographic distribution (Fig. 1) and given overlap in the diagnostic morphological characters that are supposed to separate C. taylori from other Craugastor. Specifically, Lynch (1966) reported that C. taylori was separated from other species by having large tympana (3/4 of the eye diameter), a smooth dorsum, pallid venter, tarsal fold, and lack of vocal slits. However, other species in the C. xucanebi complex often share these morphological characteristics. For example, C. glaucus lacks vocal slits in males, has a tympanum that is 3/4 of the eye, and has smooth skin (Lynch 1967a).
Evidence of rapid diversification and comparisons to previous studies
In systematics, a pattern of long branches among very short internodes is often interpreted as evidence of an ancient rapid diversification (e.g. Rothfels et al. 2012). In this study, we observed this pattern among clearly differentiated species of the bocourti series using mtDNA (Fig. 3) and concatenated nDNA (Fig. 4). Viewing the species-trees simultaneously via DensiTree demonstrates this pattern particularly well (Fig. 5). Thus, the species-tree inferences amongst species in the bocourti series are effectively a cloud-like polytomy. Given our phylogenetic results, one interpretation of this pattern is that northern rain frogs underwent a rapid radiation (Schluter 2000) early in their evolution.
Another study that included most species in the bocourti series was recently published by Portik et al. (2023). Their anuran-wide tree included ten species from the bocourti series. While the tree of Portik et al. (2023) strongly supports the monophyly of the bocourti series, it also has low branch support for relationships amongst species, as in our results here. Thus, at least two studies suggest that early, rapid radiation may have occurred in the bocourti series, leading to low branch support for most interspecific relationships. A possible way to improve support in the future would be to obtain phylogenomic data for these species.
There are other consistencies between our study and the tree of Portik et al. (2023), including the sister relationship of C. alfredi and C. silvicola (which supports that our C. aff. silvicola is indeed C. silvicola), and a shallow, strongly supported branch uniting C. xucanebi and C. stuarti. The most notable difference between our trees and the tree of Portik et al. (2023) tree is their placement of C. spatulatus and C. polymniae as sister taxa (with strong support), whereas in our mtDNA analysis we found these species to be distantly related. However, further work is needed to determine if what we call C. aff. polymniae and C. aff. spatulatus here is equivalent to the C. polymniae and C. spatulatus of Portik et al. (2023).
In a molecular study of many terraranan species, Padial et al. (2014) did not support the reciprocal monophyly of the bocourti and augusti series. Yet, these two groups were supported by many other molecular studies (Crawford & Smith 2005, Hedges et al. 2008, Portik et al. 2023). We find the recommendation of Padial et al. (2014) to ‘lump’ the bocourti and augusti series to be untenable, given the large amount of morphological and molecular data that support the reciprocal monophyly of these two series.
Misdiagnosed diversity?
We found that previous taxonomy based on morphological data overestimated species richness in some northern rain frog lineages. This finding contrasts with many previous molecular analyses of anuran species, in which morphology-based species contained one or more cryptic species (e.g. Padial & de la Riva 2009, Kieswetter & Schneider 2013). Specifically, we found that individuals from five species recognised based on morphology are contained within two shallow molecular clades (C. campbelli complex and C. xucanebi complex). Based on these low levels of genomic divergence and based on two explicit molecular species delimitation analyses (Table 5, Figs. 6 and 7), we conclude that the populations referable to C. campbelli, C. glaucus, C. nefrens, C. stuarti, and C. xucanebi, should be considered only two species (C. campbelli and C. xucanebi) and not five. Except for C. aff. nefrens, which we identified as C. campbelli, and the implied undiagnosability of C. cyanochthebius, which we also identify as C. campbelli, the taxonomy of other samples to which we applied open nomenclature remains uncertain and will need to be examined in future research.
Biologists familiar with terraranan systematics may not be surprised by our report of overestimated biodiversity because many species possess intraspecific polymorphisms in morphology (see Savage 1987, Lynch 1993). As such, delimiting species in terraranans may involve greater taxonomic uncertainty than in many other vertebrate groups. Indeed, Stuart (1941) referenced the uncertainty felt by terraranan systematists when he wrote, “I have overcome my hesitancy to further multiply Eleutherodactylid names and herein describe them…” in his description of C. xucanebi. While molecular data have led to ‘multiplied names’ in many terraranan groups (e.g. Crawford et al. 2010) and have revealed several undescribed species in the bocourti series, our study is also an important reminder that molecular data can also play an important role in countering overestimation of species diversity, especially in groups that are morphologically conservative, polymorphic, and rarely encountered.
Acknowledgements
We are indebted to J.A. Campbell for assistance with many aspects of this study related to fieldwork and financial support. We greatly appreciate the use of samples and/or photographs provided by M.E. Acevedo, R. García Anleu, S.M. Rovito, T.J. Devitt, and J.R. Mendelson III. We thank the following individuals (and their respective institutions) for allowing access to type specimens in their care: A. Resetar, T.F. Lian, and R. Inger (FMNH; Field Museum of Natural History), G. Schneider (UMMZ; University of Michigan Museum of Zoology), T. Hibbits (TCWC; Texas Cooperative Wildlife Collection), A. Wynn, R. Heyer, and R. Wilson (USNM; Smithsonian Institution Museum of Natural History), B. Clarke, D. Gower, and M. Wilkinson (BMNH; British Museum, Natural History, currently Natural History Museum, London), A. Ohler (MNHNP; Museum National d'Histoire Naturelle), G. Pandelis (UTA; University of Texas at Arlington), C. Phillips and M. Dreslik (INMH; Illinois Natural History Survey), and R. Brown and A. Campbell (KU; University of Kansas). Specimens used in this study were previously collected under the ethical and legal oversight of UNAM, UTA and Operation Wallacea. Funding for this research was provided in part by NSF grants to J.A. Campbell (Mexico, DEB-0613802 and 0102383; Guatemala, DEB-9705277), UTA start-up funds, a Bioclon grant, and NSF grant (DEB-0416160) awarded to E.N. Smiths and support from Operation Wallacea to M. Jocqué. We thank B. Hedges and A. Crawford for providing important information on Craugastor samples. For helpful comments on the manuscript, we thank A. Crawford and an anonymous reviewer. We thank P. Kok for the invitation to submit this work to the Special Issue and L. Polačiková for assistance with editing the accepted manuscript.
This is an open access article under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits use, distribution and reproduction in any medium provided the original work is properly cited.
Author Contributions
J. Streicher and E. Smith designed the study and examined museum specimens; J. Streicher generated and analysed the molecular data; J. Streicher, J. Wiens and E. Smith contributed to the analytical design; O. García-Vázquez, M. Jocqué, J. Streicher and E. Smith participated in fieldwork; J. Streicher, J. Wiens and E. Smith wrote the paper with input from O. García-Vázquez and M. Jocqué. All authors approved of the final version of the manuscript.
Data Availability Statement
DNA alignments and phylogenetic trees are available on the NHM Data Portal at https://doi.org/10.5519/7mmve3n2.