Finding complexity in complexes: assessing the causes of mitonuclear discordance in a problematic species complex of Mesoamerican toads

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.

Methods
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.

线粒体-核基因组冲突（mitonuclear discordance）是系统地理学研究中屡见不鲜的格局，当线粒体DNA与核DNA呈现冲突信号时便会出现该现象。这类遗传标记间的冲突可由多种因素导致，包括分类学混淆、基因流以及不完全谱系分选（incomplete lineage sorting）。本研究针对中美洲乔蒂斯地块（Chortís Block）分布的蟾蜍科（Bufonidae）Incilius coccifer复合群，报道了一例典型的线粒体-核基因组冲突案例。为解析该复合群内线粒体-核基因组冲突的成因，我们采用空间显性遗传数据开展物种界定与系统发育关系检验、种群历史动态特征解析以及基因流量化分析。研究发现，该类群内3个已认定物种间存在广泛的线粒体-核基因组冲突，尤以洪都拉斯乔蒂斯高地的种群最为显著。数据显示，乔蒂斯高地种群存在核基因渐渗现象，其最可能的驱动因素为气候波动引发的周期性分布范围扩张。尽管我们证实核基因组内发生了基因渐渗，但数据表明其并非驱动整个物种复合群出现线粒体-核基因组冲突的关键因素。鉴于线粒体DNA与核DNA间未呈现可辨识的地理分布格局，且该复合群的分化时间相对较晚，我们最终判定线粒体-核基因组冲突由不完全谱系分选所致。本研究为检验线粒体-核基因组冲突的来源提供了分析框架，并强调了利用多标记类型验证隐存物种种群边界的重要性。

方法
我们参照Peterson等（2012）描述的实验流程与Streicher等（2014）指定的参数，为84个个体获取了双酶切位点相关DNA测序（ddRADseq）数据。最终构建的文库于德克萨斯大学基因组测序与分析中心（GSAF）的Illumina HiSeq2500测序通道（150 bp单端读段）中完成测序，相关平台链接为https://www.wikis.utexas.edu/display/GSAF。数据处理、过滤与格式化流程采用Portik等（2017）公开的自动化脚本，脚本仓库链接为https://github.com/dportik/Stacks_pipeline。简言之，原始Illumina读段通过Stacks v1.35软件进行去多重测序（Catchen等，2013），利用fastx-toolkit（官网为www.hannonlab.cshl.edu/fastx_toolkit）的fastx_trimmer模块去除限制性酶切位点粘性末端，并通过FastQC v0.10.1软件对每个样本的测序质量进行评估，FastQC官网为www.bioinformatics.babraham.ac.uk/projects/fastqc。分别通过ustacks、cstacks和sstacks工具完成基因座的构建、编目与鉴定。随后使用populations程序为在70%所有个体中均存在的基因座生成等位基因，最终得到2211个基因座。通过自定义过滤步骤移除了不变基因座（n=150）、非双等位基因座（n=2）以及至少包含1个个体存在多于2个等位基因的基因座（n=854）。对于包含多个单核苷酸多态性（SNP）位点的基因座（每个基因座平均SNP数量为2.13 ± 1.14），我们随机选取1个SNP用于后续分析。剔除了在超过60%基因座上存在数据缺失的样本。完成上述过滤步骤后，最终的SNP数据集包含64个样本与1207个基因座。