DOI: 10.5061/dryad.q573n5tfw
title: Finding complexity in complexes: assessing the causes of mitonuclear discordance in a problematic species complex of Mesoamerican toads

description: 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.
description: 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.