Delimitation despite discordance: Evaluating the species limits of a confounding species complex in the face of mitonuclear discordance

Abstract The delimitation of species is an essential pursuit of biology, and proper taxonomies are crucial for the assessment and conservation management of organismal diversity. However, delimiting species can be hindered by a number of factors including highly conserved morphologies (e.g., cryptic species), differences in criteria of species concepts, lineages being in the early stages of the speciation or divergence process, and discordance between gene topologies (e.g., mitonuclear discordance). Here we use a taxonomically confounded species complex of toads in Central America that exhibits extensive mitonuclear discordance to test delimitation hypotheses. Our investigation integrates mitochondrial sequences, nuclear SNPs, morphology, and macroecological data to determine which taxonomy best explains the divergence and evolutionary relationships among these toads. We found that a three species taxonomy following the distributions of the nuclear SNP haplotypes offers the best explanation of the species in this complex based off of the integrated data types. Due to the taxonomic instability of this group, we also discuss conservation concerns in the face of improper taxonomic delimitation. Our study provides an empirical and integrative hypothesis testing framework to assess species delimitation hypotheses in the face of cryptic morphology and mitonuclear discordance and highlights the importance that a stable taxonomy has over conservation‐related actions.

DNA (mtDNA) and/or a limited number of nuclear DNA (nuDNA) markers, though in recent years large genomic datasets that include hundreds to thousands of loci have become much more commonplace (Fennessy et al., 2016;Hofmann et al., 2019;Thielsch et al., 2017;Victor, 2015). This has provided greater taxonomic resolution for many groups of organisms and has given us a true picture of the evolutionary relationships of species (Funk & Omland, 2003;Spinks et al., 2014). While adding more molecular data to taxonomic studies can be useful, it can also sometimes confound the taxonomic and evolutionary relationships of species even further. This is often demonstrated when you have incongruent evolutionary histories, such as between mitochondrial and nuclear loci ("mitonuclear discordance"; Toews & Brelsford, 2012). Mitonuclear discordance can arise from a number of factors including introgression, incomplete lineage sorting (ILS), and sex-biased dispersal (Firneno et al., 2020;Funk & Omland, 2003;Ivanov et al., 2018;Phuong et al., 2016;Toews & Brelsford, 2012). While mtDNA markers have proven useful through molecular barcoding for the detection of cryptic species, there have been a number of cases where mitochondrial-based methodologies conflict with traditional taxonomy and genomic inferences (Toews & Brelsford, 2012). The addition of larger genomic datasets, along with a growing prevalence of discordance between genetic datasets, has caused for the need to incorporate more lines of evidence into species delimitation in an integrative taxonomic framework, as well as the reevaluation of cryptic species and/or confounding species complexes (Padial et al., 2010).
Here, we address a controversial issue of species delimitation in a confounding complex of Central American true toads (Anura: Bufonidae: Incilius coccifer complex). Since the taxonomic revision of the species in this complex, which was based on morphology and mtDNA by Mendelson et al. (2005), the nominal species have been the subject of continuous debate and reevaluation as a complex (Firneno & Townsend, 2019;McCranie, 2015;McCranie & Castañeda, 2007;Mendelson et al., 2011). The I. coccifer complex is currently comprised of three distinct, yet closely related mitochondrial lineages that exhibit a highly conserved morphology: I. coccifer in the Pacific lowlands from southern Mexico to northern Costa Rica, I. ibarrai in the highlands from western Guatemala to western Honduras, and I. porteri in the highlands of central Honduras (Firneno & Townsend, 2019;Mendelson et al., 2005Mendelson et al., , 2011. The recent addition of genomic SNP data found extensive mitonuclear discordance between the two genetic datasets, which was determined to be the cause of incomplete lineage sorting (ILS) of mitochondrial haplotypes (Firneno et al., 2020). The conflicting genetic evidence has raised a complicated example of species delimitation where a stable and accurate taxonomy is of particular concern, as these organisms have been constantly reassessed as being data deficient, endangered, or least concern by the IUCN Red List of Threatened Species and exist in some of the most disturbed and threatened habitats in Central America (Mendelson et al., 2004;Whitfield et al., 2016).
The mitonuclear discordance found in the I. coccifer complex involves differences in tree topology and the distribution of populations/species (Firneno et al., 2020;Firneno & Townsend, 2019;Mendelson et al., 2005Mendelson et al., , 2011, therefore making it unclear exactly how many species there are in the complex and what their geographic range distribution is. Here we integrate and reevaluate previously collected molecular data (mtDNA-Firneno & Townsend, 2019;Mendelson et al., 2005; genome-wide SNPs- Firneno et al., 2020), previously and newly collected morphology (Mendelson et al., 2005), and macroecological modeling in a robust statistical and comparative integrative taxonomic framework. We use three hypotheses to test and reevaluate the species limits of the I. coccifer complex in the face of mitonuclear discordance: (a) there are three species in the complex that follow the distributions proposed in Mendelson et al. (2005) according to mtDNA and morphological data (currently recognized taxonomy); (b) there are three species in the complex that follow the distributions proposed in Firneno et al. (2020) according to SNP data and biogeographic barriers of the region; or (c) there are some other number of species in the complex that will be revealed by integrating molecular, morphological, and macroecological data that have different range distributions than the currently proposed ones. We test these hypotheses and define species in this paper using the general lineage concept (GLC), due to its pluralistic approach and the ability to use numerous criteria to define a species under its definition (Mayden, 1997;Queiroz, 1998Queiroz, , 1999Sangster, 2014). Our study highlights the importance of using multiple genetic markers and multiple lines of evidence in an integrative hypothesis testing framework to resolve complex species delimitation scenarios where different lines of evidence may conflict, as well as the importance of a stable taxonomy to guide proper conservation assessment measures.
populations was then used to generate alleles for loci present in 80% of all individuals, which resulted in 3,225 loci. Custom filtering removed "blank" loci (n = 87), invariant loci (n = 1,120), nonbiallelic loci (n = 0), and loci containing at least one individual with more than two alleles (n = 0). For loci containing multiple SNP sites, we randomly chose a single SNP to be used for subsequent analyses. Samples missing data for more than 50% of loci were removed. After completing the above filtering steps, our final SNP dataset consisted of 64 samples and 2,018 loci.
We also created a complimentary mtDNA dataset of the cytochrome oxidase I (COI) gene using available sequences on GenBank (Firneno & Townsend, 2019;Firneno et al., 2020). Individuals in this dataset either had complementary RADseq data associated with them or were from roughly the same sampling sites (Table S1; Figure 1a,b). We generated a multiple alignment of 610 bp for 59 samples (including two outgroups; Table S1) in mega7  with the muscle algorithm (Edgar, 2004) using the default parameters.

| Population structure
We determined the number of discrete populations present across the sampled range of the I. coccifer complex with our RADseq dataset using a combination of Bayesian and likelihood clustering analyses and multivariate methods. We used structure v2.3.4 (Falush et al., 2003;Pritchard et al., 2000) to examine the number of population clusters and potential admixture between populations in our dataset using MCMC. Hierarchical analyses were performed for 10 runs per K, up to a maximum of six populations, and used the admixture model with a burn-in of 10,000 steps followed by 100,000 steps. We summarized our results using structureharvester (Earl & vonHoldt, 2012) and evaluated the number of populations based on inspection of likelihood plots and following Evanno et al. (2005). To complement our structure analysis, we used a maximum likelihood approach with admixture (Jombart, 2008). We performed ten replicate analyses to evaluate up to seven populations. To assess the best K value, we performed 10-fold cross-validation and determined the K value with the lowest cross-validation error. We also evaluated the number of discrete populations using a discriminant analysis of principal components (DAPC) with adegenet v2.0.0 (Jombart, 2008;Jombart & Ahmed, 2011). A maximum of 10 clusters were investigated using the k-means algorithm. The preferred number of clusters was evaluated using BIC scores. We explored a range of three to five clusters to describe using DAPC. To minimize overfitting, an initial DAPC was used to find the a-score for each set of clusters and this value was used to select the number of principal components to retain in a subsequent reanalysis (Jombart, 2008;Jombart & Ahmed, 2011). Group membership probabilities were then examined for each cluster. To independently assess the validity of population differentiation and assignment, we used the fineradstructure software package (Malinsky et al., 2018) to construct a coancestry matrix from our RADseq data. We used a 100,000 burn-in followed by 100,000 MCMC steps sampling every 1,000 steps, and the tree was constructed with 10,000 hill-climbing iterations. The results were visualized using the FINERADSTRUCTUREPLOT.R and FINESTRUCTURELIBRARY.R scripts (included in the fineradstructure package file).

| Phylogenetic analyses
We estimated the phylogenetic relationships independently for our mtDNA and SNP datasets. We used a maximum likelihood approach carried out in RAxML v8.0 (Stamatakis, 2014) on both genetic datasets with 1,000 bootstrap replicates under a GTR substitution model. We used a Bayesian approach for our mtDNA dataset by employing beast v2.5.1  using an HKY model of nucleotide substitution, relaxed clock model, coalescent constant population, and a random starting tree with a Markov Chain Monte Carlo (MCMC) run for 2 × 10 7 generations, sampling every 1,000 generations producing a total of 10,000 trees. We assessed the run using Tracer v1.6 (Rambaut & Drummond, 2009) to examine convergence. A burn-in of 10% was discarded, and a maximum clade credibility tree with median heights was created from the remaining 9,000 trees. We then estimated the species tree for our SNP dataset using snapp v1.3.0 (Bryant et al., 2012) implemented in beast2.
To reduce run times, we subsampled each population to include 3-6 representatives, for a total of 20 individuals. We estimated the mutation rates (u and v) from the data (1.06 and 0.94, respectively) within beauti. We assigned a gamma distribution to our birth rate ( ) of the Yule prior, with and alpha of 2.0 and a beta of 2.0. Our snapp prior was assigned an alpha of 11.75, a beta of 109.73, and a lambda of 0.01. We performed two independent runs with a chain length of 1,000,000 generations, sampling every 1,000 generations. Runs were assessed using tracer v1.6 (Rambaut & Drummond, 2009) to examine convergence, and tree topologies and node heights were visualized using densitree (Bouckaert, 2010).

| Species delimitation
To estimate the best-fit number of species units based on our two genetic datasets, we used a suite of different species delimitation analyses. We applied three single locus delimitation analyses to our mtDNA dataset including (a) the General Mixed Yule Coalescent (GMYC) model, which establishes thresholds between the branching patterns of ultrametric gene trees from inter-and intraspecific branches in order to define speciation events (Fujisawa & Barraclough, 2013); (b) the Poisson Tree Process (PTP; Kapli et al., 2016), which, like GMYC, aims to identify the transition between inter-and intraspecific processes, but specifically requires that the phylogenetic gene trees branch lengths are proportional to the number of substitutions, rather than to time as for GMYC; (c) and the Automated Barcode Gap Discovery (ABGD) method (Puillandre et al., 2012), which partitions sequences into groups based on comparisons of pairwise distances and compared to tree-based delimitation methods (GMYC and PTP) that are suspected to over split, ABGD offers a more conservative approach to estimate the number of species given comprehensive sampling (Puillandre et al., 2012).
For the GMYC analysis, we applied a single model using the R package splits to our ultrametric tree generated in BEAST. Next, we F I G U R E 1 Distribution of haplotypes in the I. coccifer complex (I. coccifer in red, I. ibarrai in yellow, and I. porteri in blue) and sampling for our (a) mitochondrial and (b) SNP datasets. The Chortís/Guatemala Highland boundary is indicated by the dashed line. (c) Coancestry matrix from fineradstructure based on our SNP dataset (coefficients of coancestry are color-coded from low (yellow) to high (black), and the dendrogram depicts a clustering of individual sampled based on the pairwise matrix of coancestry coefficients) indicating main populations (bolded squares) with substructuring identified within the sampling east of the Chortís/Guatemala Highlands boundary (dashed squares). (d) Discriminant analysis of principal components of SNP dataset. (e) Maximum likelihood cladograms of mtDNA (left) and SNP (right) datasets exhibiting extent of mitonuclear discordance. Nodes with high support (bootstrap ≥90) are indicated by black dots, and asterisks (*) indicate samples with haplotypes that do not follow the exact pattern of mitonuclear discordance shown by dashed lines. (f) Summary of species delimitation analyses, with collapsed cladogram inferred from BEAST analysis of mtDNA (left) and SNAPP species tree of SNP dataset (right). Nodes with high support (posterior probability ≥0.9) are indicated by black dots. Blocks represent species units that have been estimated for each analysis. Gray areas connecting blocks indicate that the support for the split between those taxa is low to moderate (gdi between 0.1 and 0.7) implemented the bPTP analysis using our maximum likelihood tree on the bPTP server (https://speci es.h-its.org; Zhang et al., 2013). For this analysis, we ran 100,000 MCMC generations, with a thinning of 100 and burn-in of 0.1. Finally, we applied the AGBD analysis to our data using the web interface (https://bioin fo.mnhn.fr/abi/publi c/abgd). For our analysis of each locus dataset, we used the default maximum intraspecific distance values (P max = 0.1), a relative barcode gap width of X = 1.5, and a distance correction of Jukes-Cantor (JC69). Default settings were used for all of the remaining parameters.
We then applied four species delimitation methods to our e.g., implementing BFD*). We performed our analyses testing four models, which are presented in Table S2. We compared and ranked models to select the best-supported species hypothesis. We calculated the Bayes Factor (BF) by subtracting the value of the log MLE for the model representing the current taxonomic classification from each alternative model and multiplying the difference by two (BF = 2(model 1 − model 2)). We ranked all of the models and selected the model with the highest BF (Table 1).
We used the same parameters as our snapp analysis (see above).
We performed 48 path sampling steps, with 100,000 MCMC generation and a preburn-in of 1,000 generations. The BF was then calculated for each alternative model, the models were ranked, and a best model was chosen (Table S2).
Next, we used Bayes species delimitation as implemented in BPP v4.2 (Flouri et al., 2018(Flouri et al., , 2020. We implemented the A10 (species delimitation using a fixed guide tree) analysis under the MSC accommodating introgression (MSCi) model using relationships derived from phylogenetic analyses as a guide tree. The mapping of individuals to their identifiers was done in the imap file supplied to BPP.
Following this, we used the gdi to determine whether putative species boundaries correspond to species-level divergences between populations of the I. coccifer complex [25][26]. We calculated gdi values in BPP v4.2 (Flouri et al., 2018(Flouri et al., , 2020) under the MSCi model. The MSCi approach provides a measure of uncertainty in population estimates, and gdi can be calculated using the posterior distributions for θ and τ (Leaché et al., 2019). The posterior probability distributions for θ and τ were estimated in the A00 analysis (Jackson et al., 2017) using a fixed species tree containing our three lineages (I. coccifer, I. ibarrai, and I. porteri, based on SNP haplogroups). We used the same diffuse inverse-gamma priors for θ and τ that we used in our BPP species delimitation analysis (see above).
To assess convergence, we compared the posterior distributions from three independent runs (same parameters as BPP above). We then calculated gdi for each species by combining all samples from the posterior distributions using the equation: gdi = 1 − e −2τ/θ (Chan & Grismer, 2019;Leaché et al., 2018). Lineages were considered distinct species when gdi > 0.7, the same species when gdi < 0.2, and/ or ambiguous species status if 0.2 ≥ gdi ≤ 0.7 (Jackson et al., 2017;Pinho & Hey, 2010).
Finally, we used delineate (Sukumaran et al., 2021) to determine the number of species that are within the I. coccifer complex under the PSM. We used the species tree from our SNAPP analysis as a guide tree, and within our control file, we constrained the I. coccifer lineage (based on previous descriptions of morphological and ecological distinctiveness) but left the I. ibarrai and I. porteri (highland) lineages as unconstrained to determine if they are separate species from the I. coccifer or each other.

| Morphological data acquisition and processing
To determine which genetic species assignment (mtDNA or SNP) and

| Species distribution modeling
To determine the extent of macroecological differentiation between genetically discordant populations, we created species distribution models (SDMs) for all three species based on their distributions according to their respective genetic datasets. We compiled all available museum locality information from vertnet (Constable et al., 2010) and our own sampling for toads of the I. coccifer complex. We created two datasets corresponding to our genetic hypotheses (mtDNA or SNP) and assigned species into their respective genetic datasets either based on genetic confirmation or being within the distribution of a genetic haplotype. This resulted in 410 records that were used to create SDMs (Table S3) and BIO19 = precipitation of coldest quarter) and six envirem variables (annual potential evapotranspiration (PET), thermicity index, climatic moisture index, aridity index, terrain roughness index, and topographic wetness), bringing our final dataset to a total of seventeen variables. We constrained our environmental variable layers to the range spanned by all of our occurrence records.
To build SDMs, we used maxent (Phillips & Dudík, 2008)  We inferred the amount of niche overlap between species for the mtDNA distribution and the SNP distribution SDMs using by calculating Schoener's D and I by pairwise comparison using the nicheOverlap functions in the R package dismo (Hijmans et al., 2011).
These metrics give an output value from 0 to 1, where a value of 0 indicates no overlap between niches/low niche similarity and a value of 1 indicates that the niches completely overlap/are identical.

| Genetic and species delimitation analyses
In agreement with previous studies (Firneno et al., 2020), we see extensive mitonuclear discordance among our two genetic datasets of individuals within the Incilius coccifer complex. We estimated four populations from our SNP dataset, with a single lowland population, a highland population that exist west of the Chortís/Guatemala Highlands boundary, and two populations with admixture east of the Chortís/Guatemala Highlands boundary (Figures 1d and S1). Our fineradstructure analysis identified three main populations (as defined above) with some substructuring among individuals in the population east of the Chortís/Guatemala Highlands boundary, which may be due to nuclear introgression or isolation-by-distance due to the mosaic nature of the organisms' highland distributions (Figure 1c).

| Morphological analyses
The morphological variation across the species complex best supports the two species hypothesis, with an overall correct classification rate of 92.4% for both mtDNA and SNP assignments. The three species and four species hypotheses did have high classification rates but were 10.1%-12% lower than the DFAs using a two species classification (Table 1). We found that I. coccifer is a morphologically distinct taxon, which had a correct classification rate of 94.9% across all six of the DFAs and showed a high degree of separation in the PCAs (Table 1; Figure 2a). The decreased rate of correct classification in the three and four species DFAs was associated with delineating the Highland taxon into multiple species (Table 1). These populations, which are currently recognized as I. porteri and I. ibarrai, show a high degree of morphological similarity. Specifically, it is the Honduran populations that are very similar.
The VP analysis found that the nuDNA and mtDNA assignments both significantly explained the dataset (Figure 2b). The nuDNA and mtDNA assignments explained 36% and 32% of the morphological variation, respectively. The VP analysis revealed that 27% of the explained variation was shared by both assignments (Figure 2b).

F I G U R E 2
Analyses of the association between species delimitation and morphometry. (a) Visualization of the first two PCA loadings of the morphological variation. Convex polygons were produced using the mtDNA and SNP species assignments for across our molecular species delimitation hypotheses to show the degree of morphological variation and overlap between the species. (b) Venn diagram of the results of the VP analysis that depicts the adjusted R 2 values and significance levels of the mtDNA-and SNP-based delimitations for discordant highland populations (I. ibarrai and I. porteri). Values outside of the shaded areas represent marginal effects (R 2 adj; e.g., the amount of variation explained when testing each delimitation separately). Values in the intersection of the shaded area represent variation explained in common by both delimitations. Values outside of the intersection of the shaded area represent the conditional effects (variation uniquely explained by each delimitation). *p <.001

| Ecological niche modeling and niche overlap
The top candidate SDMs that constituted our final models had high AUC scores ranging from 0.81 to 0.95. Based on the known geographic distributions of the species under either phylogenetic hypothesis, very little geographic overestimation occurred in the models. The only model that we see some overestimation into more lowland regions is the I. ibarrai SNP model ( Figure S4), which we attribute to the low number of samples (<20) that have been used for that model. It should be noted that there is evident estimation of the I. ibarrai and I. porteri models of both datasets ( Figure S4)

| D ISCUSS I ON
Delimiting species can be hindered by a number of factors including highly conserved morphologies (e.g., cryptic species), differences in criteria of species concepts, lineages being in the early stages of the speciation or divergence process, and discordance between gene tree topologies (e.g., mitonuclear discordance). Here we provided an example of a taxonomically confounding species complex of toads from Central America that exhibits extensive mitonuclear discordance and presented an integrative hypothesis testing framework to assess differentiation of integrated data using the different phylogenetic hypotheses caused by mitonuclear discordance. Below we summarize our findings and evidence for the taxonomic relationships of the species within the I. coccifer complex, along with a discussion of the conservation implications for this species complex.

| Species delimitation
While strides have been made in species delimitation using genetic data (Fujita et al., 2012;Leaché et al., 2018;Luo et al., 2018), taxonomists are still required to make subjective judgments based on multiple separate lines of evidence to decide on the boundary between populations and species (Sites & Marshall, 2004). Following single lines of evidence, such as a single genetic locus, can lead to discordant results that do not reflect an accurate taxonomy (Jackson et al., 2017). We chose to implement single locus species delimitation methods on our mtDNA dataset, and they suggest that there are between 3 and 5 species within the complex. Our tree-based single locus methods (GMYC and PTP) delineate a species boundary between I. ibarrai populations in Guatemala and those in Honduras, which have been previously suggested as two potentially distinct species (Firneno & Townsend, 2019;Mendelson et al., 2011). This Single locus species delimitation methods rely on the topologies inferred from a single locus and assume that it represents the current species relationship. However, incorrect topologies due to introgression and/or ILS can mislead these methods (Dupuis et al., 2012;Knowles & Carstens, 2007). Because ILS was found to be the cause of mitonuclear discordance in the I. coccifer complex and that nuclear introgression was present within populations of I. porteri within the Chortís Highlands (Firneno et al., 2020), we also used species delimitation methods that incorporate multilocus data in both a MSC and PSM framework. These methods have been shown to take into account ILS and are robust to low levels of introgression (Zhang et al., 2011). Our multilocus delimitation methods widely varied estimating 1-4 species units. The methods that apply the MSC (BFD* and BPP) tended to delimit species along the lines of our population inferences, whereas our other methods (gdi and delineate) were more conservative in their estimates, indicating only one or two species units present. An often-cited issue with species delimitation methods that use the MSC framework (BFD* and BPP) is that they tend to delimit species boundaries based on population structure, therefore overestimating the number of species within a focal group (Sukumaran & Knowles, 2017). We see this in our data and is the reason why we used comparative methods that implement other frameworks other than the MSC. It has been remarked that the use of gdi can have a number of weaknesses including that the criterion depends on the population divergence time relative to the population size (small populations and recent divergence times may skew results), that ambiguity may arise when two populations being compared are of very different sizes, and that the metric has a wide range of indecision and reflects a more subjective nature of species delimitation that these methods are trying to avoid (Leaché et al., 2019). Likewise, it has been argued that the PSM implemented by delineate is unrealistic because it models instantaneous speciation in a single generation, which does not seem to reflect how natural species convert from being an incipient species to a true species (Leaché et al., 2019;Sukumaran & Knowles, 2017). These may be reasons why our heuristic and PSM-based analyses yielded such conservative results.
The addition of more genetic markers can often aid in the clarity of delimiting species. However, there are cases, as exemplified here, where adding more genetic data can reveal more complexity within and between species than originally assumed or hinder the ability of taxonomists to make objective decisions concerning the delimitation of species (Hinojosa et al., 2019;Ivanov et al., 2018;Papakostas et al., 2016;Pedraza-Marrón et al., 2019). Due to the discordance between our genetic datasets and the wide array of species units estimated by our species delimitation analyses, it is difficult to discern the structure and number of species in this complex based on genetic data alone. However, we were able to use these competing phylogenetic and delimitation hypotheses to further integrate and evaluate other lines of evidence (e.g., morphology and macroecology) to determine which hypothesis best describes the variation seen in these other data types and most accurately reflects the proper taxonomic relationships and distributions of the species in this complex.

| Taxonomic implications
Under the GLC, the only necessary property of a species is its existence as an independently evolving lineage (de Queiroz, 2005(de Queiroz, , 2007.
Following our molecular phylogenetic and species delimitation methodologies under the GLC, we cannot clearly and objectively make an inference of the number of species based on our genetic data.
Also, due to mitonuclear discordance it is still not clear which dataset (mtDNA or SNP) gives us a more concrete indication of where the species range limits lie in this complex. Therefore, we integrated morphological and macroecological analyses to further test hypotheses of species number and distribution. These combinations of data in an integrative taxonomic framework have long been shown to clarify problematic species complexes such as the I. coccifer complex (Padial et al., 2010). It should be noted that, though macroecological modeling is not often used for species delimitation, when combined with genetic and/or morphological information, species delimitation hypotheses can improve in objectivity and strength (Hidalgo-Galiana et al., 2014;Leaché et al., 2009;Razkin et al., 2016;Rissler & Apodaca, 2007).
Though Firneno et al. (2020) suggest that there may be more than three species in this complex, a hypothesis that was based only on genetic data, we do not find evidence within our genetic, morphological, and macroecological data for this claim. Within our morphological data, we see little separation between highland populations within Honduras ( Figure 2a) and a low amount of support for discerning these populations (Table 1) when they are treated as separate species; however, we do see good support for a separate highland population in Guatemala. Instead, we see high support for only two species (a lowland entity and a highland entity) and moderate support for three species (a lowland entity and two highland entities) with relatively high overlap between the two highland entities. Our macroecological analyses clearly show separation of suitable habitat between at least two species (Figure 3), but with higher niche overlap when splitting highland entities into two species/populations. It is also worth noting that, while we could not test a four species hypothesis within our macroecological analyses (the locality datasets upon spatial rarefication become too small to generate reliable models), we would expect high levels of niche overlap as we see in our current three species models (Figure 3). Along with the species limits, we were also concerned with determining which distribution (mtDNA or SNP) reflected the taxonomy as well. We found support for a more parsimonious distributional  Rovito et al., 2015). It may also be worth noting that I. porteri's distribution may extend into the highlands of Nicaragua (still comprising the Chortís Highlands), due to the lack of any biogeographic separation. However, no confirmed samples have ever been reported from this region.
Ultimately, we were able to use integrated lines of evidence in a hypothesis testing framework to evaluate the species limits and distributions of the three species found within this highly confounding species complex. A major factor that may also be confounding this species complex is that it is a relatively young complex (divergence was estimated at ~790,000 years ago; Firneno et al. (2020). Young lineages may still be undergoing the process of speciation/divergence, which can prove challenging for the types of delimitation methods that we bring to bear (Fujita et al., 2012;Padial et al., 2010). In the future, more taxon-specific lines of evidence (e.g. bioacoustics, osteology, etc.) may be useful in clarifying and solidifying this complex's taxonomy further, as well as using a process-based delimitation approach that integrates processes of speciation into delimitation practices (Smith & Carstens, 2020). Finally, due to the long-standing taxonomy and distributions that were originally proposed by Mendelson et al. (2005), there have been existing conservation recommendations that may now need to be reevaluated based upon these new findings.

| Conservation implications
A stable (Smith & Carstens, 2020) taxonomy is necessary to guide conservation assessment measures and practices for species.
Central American amphibians, such as those within the I. coccifer complex, represent some of the world's most vulnerable organisms due to rapid habitat loss and zoonotic disease (e.g., chytridiomycosis) in this region (Gallant et al., 2007;James et al., 2015;Lips, 1998Lips, , 2016Whitfield et al., 2016) While we sought to evaluate the species limits of the Incilius coccifer complex, we also wanted to provide an integrative hypothesis testing framework to aid in the delimitation of species complexes in the face of mitonuclear discordance. Our results highlight the importance of using comprehensive and integrated data for addressing complex species delimitation problems. Our findings further emphasize the necessity of an accurate taxonomy to guide conservation measures, especially among organisms that inhabit ecosystems that are highly imperiled. With this integrative taxonomic approach, we hope to fuel a novel perspective on pursuing species delimitation and phylogeographic work using complex and cryptic species and work within this region of Central America where amphibians (and other organisms) are facing grave conservation concerns due to anthropogenic changes to their ecosystems.

ACK N OWLED G M ENTS
We thank S. Lainez (Section of Protected Areas and Wildlife, Instituto Rivera, for their time and help with some of the analyses that were utilized in this study, and the revising and editing process for this manuscript.

CO N FLI C T O F I NTE R E S T
We declare that we have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
We have included a large data package on DRYAD (https://doi. org/10.5061/dryad.j9kd5 1cb5) that includes our final ddRADseq filtered "haplotypes" files, COI alignment files, resulting input files for a number of analysis programs (structure, admixture, DAPC, fineradstructure, BFD*, BPP, gdi, delineate, beast, and snapp), an excel spreadsheet of raw morphological measurements split into the two distributional datasets used, and a.csv file of both locality datasets used for SDMs.