Mitonuclear discordance is a frequently encountered pattern in phylogeographic studies and occurs when mitochondrial and nuclear DNA display conflicting signals. Discordance among these genetic markers can be caused by several factors including confounded taxonomies, gene flow, and incomplete lineage sorting. In this study, we present a strong case of mitonuclear discordance in a species complex of toads (Bufonidae: Incilius coccifer complex) found in the Chortís Block of Central America. To determine the cause of mitonuclear discordance in this complex, we used spatially explicit genetic data to test species limits and relationships, characterize demographic history, and quantify gene flow. We found extensive mitonuclear discordance among the three recognized species within this group, especially in populations within the Chortís Highlands of Honduras. Our data reveal nuclear introgression within the Chortís Highlands populations that was most likely driven by cyclical range expansions due to climatic fluctuations. Though we determined introgression occurred within the nuclear genome, our data suggest that it is not the key factor in driving mitonuclear discordance in the entire species complex. Rather, due to a lack of discernible geographic pattern between mitochondrial and nuclear DNA, as well as a relatively recent divergence time of this complex, we concluded that mitonuclear discordance has been caused by incomplete lineage sorting. Our study provides a framework to test sources of mitonuclear discordance and highlights the importance of using multiple marker types to test species boundaries in cryptic species.

We collected ddRADseq data for 84 individuals following the protocol described in Peterson et al. (2012) and following parameters specified in Streicher et al. (2014). Our final library was analyzed on one Illumina HiSeq2500 lane (150 bp single end reads) at the Genomic Sequencing and Analysis Facility (GSAF) at The University of Texas (https://www.wikis.utexas.edu/display/GSAF). The workflow for data processing, filtering, and formatting was automated using scripts available from Portik et al. 2017 (https://github.com/dportik/Stacks_pipeline). In brief, the raw Illumina reads were demultiplexed using stacks v1.35 (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013), the restriction site overhangs were removed using the fastx_trimmer module of the fastx-toolkit (www.hannonlab.cshl.edu/fastx_toolkit), and the sequencing quality was examined on a per sample basis using fastqc v0.10.1 (www.bioinformatics.babraham.ac.uk/projects/fastqc). Loci were created, catalogued, and identified using ustacks, cstacks, and sstacks, respectively. populations was then used to generate alleles for loci present in 70% of all individuals, which resulted in 2,211 loci. Custom filtering removed invariant loci (n=150), non-biallelic loci (n=2), and loci containing at least one individual with more than two alleles (n=854). For loci containing multiple SNP sites (average number of SNPs per locus was 2.13 (± 1.14)), we randomly chose a single SNP to be used for subsequent analyses. Any samples missing data for more than 60% of loci were removed. After completing the above filtering steps, our final SNP dataset consisted of 64 samples and 1,207 loci.