Vanishing refuge? Testing the forest refuge hypothesis in coastal East Africa using genome‐wide sequence data for seven amphibians

High‐throughput sequencing data have greatly improved our ability to understand the processes that contribute to current biodiversity patterns. The “vanishing refuge” diversification model is speculated for the coastal forests of eastern Africa, whereby some taxa have persisted and diversified between forest refugia, while others have switched to becoming generalists also present in non‐forest habitats. Complex arrangements of geographical barriers (hydrology and topography) and ecological gradients between forest and non‐forest habitats may have further influenced the region's biodiversity, but elucidation of general diversification processes has been limited by lack of suitable data. Here, we explicitly test alternative diversification modes in the coastal forests using genome‐wide single nucleotide polymorphisms, mtDNA, spatial and environmental data for three forest (Arthroleptis xenodactyloides, Leptopelis flavomaculatus and Afrixalus sylvaticus) and four generalist (Afrixalus fornasini, A. delicatus, Leptopelis concolor and Leptopelis argenteus) amphibians. Multiple analyses provide insight about divergence times, spatial population structure, dispersal barriers, environmental stability and demographic history. We reveal highly congruent intra‐specific diversity and population structure across taxa, with most divergences occurring during the late Pliocene and Pleistocene. Although stability models support the existence of some forest refugia, dispersal barriers and demographic models point towards idiosyncratic diversification modes across taxa. We identify a consistent role for riverine barriers in the diversification of generalist taxa, but mechanisms of diversification are more complex for forest taxa and potentially include topographical barriers, forest refugia and ecological gradients. Our work demonstrates the complexity of diversification processes in this region, which vary between forest and generalist taxa, but also for ecologically similar species with shared population boundaries.


| INTRODUCTION
Biodiversity is unequally distributed across the earth, with the greatest concentration occurring in tropical regions (Gaston, 2000).
Understanding the processes that have generated and maintained this pattern has been a major question in biology for the best part of a century (Brown, 2013). Landscape changes are hypothesized to have driven vicariant evolution by fragmenting species distributions that were formerly continuous, generating congruent spatial and temporal patterns of genetic differentiation across co-distributed taxa (Rosenzweig, 1995;Smith et al., 2014). Vicariance may occur due to a number of different processes, including the formation of dispersal barriers by rivers (the riverine barrier hypothesis; Gascon et al., 2000;Haffer, 1997;Moritz, Patton, Schneider, & Smith, 2000;Plana, 2004;Voelker et al., 2013;Wallace, 1852), by mountains (Fjeldså & Lovett, 1997) or by areas of unsuitable habitat (Kirschel et al., 2011;Schneider, Smith, Larison, & Moritz, 1999). In some regions, landscape changes over time may have been so severe that particular areas would have functioned as refugia, while diversity in surrounding areas was entirely lost (e.g., the forest refuge hypothesis: Endler, 1982;Haffer, 1969Haffer, , 1997Mayr & O'Hara, 1986;Moreau, 1954;Moritz et al., 2000;Plana, 2004). In some cases, diversification may have been triggered by climate change, with taxa making an ecological switch from forest to non-forest habitats in order to persist (the vanishing refuge hypothesis, Vanzolini & Williams, 1981). Finally, ecotones can facilitate disruptive selection on phenotypes despite the presence of gene flow, which may lead to dispersal with incomplete genetic barriers (e.g., range expansion), ultimately driving parapatric divergence (Moritz et al., 2000;Smith, Wayne, Girman, & Bruford, 1997). To understand how and why biodiversity accumulates in tropical regions, the underlying diversification processes need to be tested (e.g., Charles et al., 2018;Ntie et al., 2017;Portik et al., 2017). However, this has remained difficult for most tropical biodiversity hotspots due to a lack of thorough taxonomic, geographic and molecular (genome-wide) sampling.
The Coastal Forests of Eastern Africa (CFEA) are one of Africa's foremost centres of diversity and a designated global biodiversity hotspot (Burgess, D'Amico Hales, Underwood, Dinerstein, & Ecoregion, 2004). Along with the adjacent Eastern Arc Mountains, the CFEA form an important area of endemism highly threatened by anthropogenic impacts (Barratt et al., 2014;Burgess & Clarke, 2000;Burgess, Clarke, & Rodgers, 1998;Burgess, Mwasumbi, Hawthorne, Dickinson, & Doggett, 1992). Distributional data from plants, vertebrates and invertebrates demonstrate a high number of narrow-ranged endemics and a pronounced north-south biogeographic divide caused by the rain shadow of Madagascar (Burgess & Clarke, 2000;Burgess et al., 1998). Diverse communities of taxa, many endemic, are found in the remaining forest fragments of Kenya (Tana River, Kwale, Arabuko-Sokoke), Tanzania (East Usambara, Pemba island, Uluguru, Udzungwa, Pugu hills and Lindi), Mozambique (Bazaruto archipelago), Malawi (Mulanje massif) and Chimanimani and Haroni-Rusitu in Zimbabwe (Burgess et al., 1998;Figure 1a). Historical environmental change across tropical East Africa has been frequent since the Miocene, and the current CFEA are considered to be the remnants of a once continuous forest that has expanded and contracted for the past 40 million years (Axelrod & Raven, 1978;Demenocal, 1995;Maslin et al., 2014;Mumbi, Marchant, Hooghiemstra, & Wooller, 2008). Combined knowledge of endemism patterns and environmental change have led to the assumption that current CFEA biodiversity mainly originated from the isolation and persistence of ancient lineages in forest refugia, with local extinctions and in some cases adaptation to non-forest habitats across the rest of the region (Azeria, Sanmartín, Ås, Carlson, & Burgess, 2007;Barratt et al., 2017a;Burgess et al., 1998). The coastal forests have thus been described as a "vanishing refuge" (Burgess et al., 1998), although to date forest refugial processes have not been thoroughly tested against alternative modes of diversification (Damasceno, Strangas, Carnaval, Rodrigues, & Moritz, 2014;Kirschel et al., 2011;Schneider et al., 1999;Zhen et al., 2017).
To test competing hypotheses of ecological association over long time periods (i.e., millions of years), it is helpful to use taxa that retain ancestral variation at small spatial scales. Amphibians are ideal candidates due to their poor dispersal abilities and highly deme-structured species in East Africa (Bittencourt-Silva et al., 2017;Blackburn & Measey, 2009;Lawson, 2013;Zimkus et al., 2017). Here, we use nuclear DNA (nucDNA) single nucleotide polymorphism (SNP) data acquired from restriction site-associated DNA sequencing (RAD-seq), mitochondrial DNA (mtDNA) and georeferenced occurrences for seven widespread amphibians found in forest and non-forest habitats across the CFEA. Although these seven taxa represent only a small proportion of the 55 known amphibians in this region (Barratt et al., 2017a), they cover a variety of life histories and ecological associations (from generalists to forest specialists), which should reflect the evolutionary processes that have occurred in this region. Our sampling represents almost the complete range of each of the study taxa, allowing us to test the forest refuge hypothesis against alternative modes of diversification, including allopatric divergence across landscape barriers (rivers and topographic) and parapatric divergence across ecological gradients. Using our analytical framework, we make an important first step towards understanding fine-scale diversification patterns and processes in this highly threatened biodiversity hotspot which is a suitable model for similar studies in other tropical regions. For each focal taxon, we assess (a) How many distinct populations are there, and how are they related? (b) Do population boundaries coincide with geographic features and dispersal barriers, and are these boundaries shared across taxa? (c) Which demographic mechanisms have played a role in population diversity and divergence?
For forest taxa which may be older in terms of their evolutionary history, we predict concordant population structure across taxa (Burgess et al., 1998), with distributions that have remained locally stable over multiple time periods (i.e., refugia), and demographic signals of allopatric divergence and size expansion/secondary contact which would represent forest contraction and expansions. Generalist taxa present in non-forest habitats such as Miombo woodland and savannah, on the other hand, should demonstrate signals of migration between populations and incomplete dispersal barriers, indicating the potential role of ecological gradients in line with the vanishing refuge hypothesis. Given the large number of major rivers, mountains and raised plateaus intersecting our CFEA sampling, we additionally investigate the role of landscape barriers, which would F I G U R E 1 Study region in East Africa, encompassing the Coastal Forests of Eastern Africa and surrounding areas, with sampling locations marked (forest and generalist taxa). (a) Major hydrological and topographical features, including forest refugia demarcated by red polygons based on data from Burgess et al. (1998). (b) Terrestrial ecoregions indicating the location of known biogeographic realms, associated with climate gradients, and sampling localities used in this study (black dots) [Colour figure can be viewed at wileyonlinelibrary.com] result in clear population structure without subsequent migration, size changes or secondary contact. We make a distinction between forest and non-forest (Miombo woodland and savannah) to aid interpretation of results with regard to the long-term (pre-Pleistocene) environmental stability of these areas and to the ecological preferences of the study taxa.
We firstly investigate population structure within taxa to define putative geographic boundaries of populations. We then use explicit demographic model selection to reach a conclusion on whether refugial models of diversification are applicable for forest and generalist taxa. We support demographic results by modelling evolutionary relationships, geographic distributions, effective migration and diversity, and long-term environmental stability based on ecological niche models (ENMS). This approach allows us to assess the congruence of patterns across taxa and evaluate their consistency with putative forest refugia, geographic features (river and topographic barriers) and known biogeographic regions caused by ecological (rainfall) gradients.

| DNA sequencing and data filtering
Genomic DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen) following manufacturer's instructions. A partial fragment of the mitochondrial 16S gene was amplified via polymerase chain reaction to verify species identity using the NCBI BLAST tool against our own barcoding database of amphibians across the region (full details can be found in Barratt et al., 2017a; GenBank Accession nos included in Supporting Information Table S1). DNA was quantified prior to RAD-seq library preparation using a Qubit fluorometer (Invitrogen). We followed the Etter, Bassham, Hohenlohe, Johnson, and Cresko (2011) laboratory protocol to prepare RAD-seq libraries using the SbfI restriction enzyme (Supporting Information Appendix S1).
We sequenced individuals across five RAD-seq libraries (45-51 samples each), with a unique barcode adapter per individual in each library to demultiplex sequences bioinformatically. The final eluted products were sequenced in a single run on an Illumina Hi-seq 2500 (100-bp single-end reads) at the D-BSSE sequencing facility in Basel, Switzerland (Supporting Information Table S2).
We used STACKS 1.41 (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013) to process RAD-seq data and produce single nucleotide polymorphism (SNP) data sets. We used the process_radtags module to filter out low-quality reads, demultiplexing individuals into their own fastq file. The standard workflow of ustacks, cstacks and sstacks modules was used to align reads into stacks, to build a catalogue of consensus loci by merging alleles across individuals and to match individuals to the catalogue, respectively (data sets summarized in Supporting Information

| Population structure
We filtered STACKS "haplotype" files to remove loci that were invariant between samples, loci with at least one individual with more than two alleles (i.e., potentially paralogous loci), and loci that were not bi-allelic. We investigated population structure per taxon using discriminant analysis of principal components (DAPC) in the ADEGENET R package (Jombart & Ahmed, 2011), after converting STACKS output files into Fstat format using PGDSPIDER 2.1.0.3 (Lischer & Excoffier, 2012). Unlike model-based clustering methods, the DAPC method is free of assumptions regarding Hardy-Weinberg equilibrium (Jeffries et al., 2015;Jombart & Ahmed, 2011) and less sensitive to minor allele frequency thresholds (Linck & Battey, 2017). We defined values of k between 1 (i.e., a single panmictic population) and 8, using Bayesian information criterion (BIC, Schwarz, 1978) scores across tested k values to infer the number of populations. Due to the taxonomic uncertainty of some samples of A. delicatus/A. sylvaticus and L. argenteus/L. concolor, we conducted an initial analysis of the full data set to confirm species memberships, followed by analyses of each taxon separately. To complement our DAPC analyses, we also ran ADMIXTURE (Alexander, Novembre, & Lange, 2009) with formatted bed files converted using PLINK 1.07 (Purcell et al., 2007). As with DAPC analyses, k ranged between 1 and 8, and we used the 10-fold cross-validation procedure to estimate the number of population clusters. Population membership of each individual for chosen k values was verified by inspecting the clustering analysis plots in DAPC and ancestry coefficients in ADMIXTURE barplots (Supporting Information Tables S4 and S5). DAPC clustering for multiple values of k is shown in Supporting Information Figure S1.

| Demographic model selection
To evaluate the likelihoods of alternative demographic models within each taxon based on RAD-seq data, we used the diffusion approxi- and represent divergence due to ecological gradients. Although these broad categorizations are simplistic, they enable comparisons across taxa to be more readily made, although it should be noted that other processes which are not explicitly captured by our model parameters may also contribute to the demographic patterns observed (e.g., range expansion, local adaptation, recent anthropogenic impacts).
We ran a total of 15 alternative 2D models that differed in parameters for migration rates, periods of isolation and population size changes (visually represented in Supporting Information Figure  | 4293 optimizations for each model before performing the final model selection. We did not transform obtained parameters into absolute migration rates and divergence times because our primary aim was to perform model selection, and these parameter values should ideally be estimated using an accurate mutation rate which we lack for our study taxa (Gutenkunst et al., 2009). We therefore compare the relative time intervals of population divergences obtained from δaδi with those obtained from a Bayesian coalescent-based approach using bi-allelic SNPs as described in the next section. To determine the best fitting models, the AIC and log likelihoods were inspected, and ΔAIC scores were used to calculate relative likelihoods and Akaike weights (ω i ). The model with the highest Akaike weight was selected as the most likely for each divergence event (Burnham & Anderson, 2002). To explore the possible influence of recent anthropogenic impacts, we ensured that a variation in the top-ranked model per taxon was also tested that included a size change step.
Two-dimensional models already had these models within the original fifteen tested, and we built an additional three 3D models to test for each of A. fornasini, L. flavomaculatus and A. xenodactyloides. To verify that our models were reasonable explanations of the JSFS, we performed goodness of fit tests. For each taxon, we fit the topranked model using our optimized parameters, scaled the resulting model spectrum by the inferred theta value and used the model spectrum to generate 100 Poisson-sampled frequency spectra. We then optimized each simulated frequency spectrum to obtain a distribution of log-likelihood scores and Pearson's chi-squared test statistic and subsequently determined whether our empirical values were contained within these distributions. A more detailed description of demographic model selection and goodness of fit tests is shown in Supporting Information Appendix S1. To support our demographic model selection, we also estimated evolutionary relationships, effective migration and diversity, and the long-term stability of taxa using ENMS, as described below.

| Evolutionary relationships
We explored evolutionary relationships with mtDNA and RAD-seq data separately. Sequences of 16S mtDNA for all samples included in RAD-seq libraries were edited in GENEIOUS 6 (Kearse et al., 2012) and aligned with the RAxML tree estimator using a GTRCAT model in SATé-II (Liu et al., 2012) before analyses in BEAST 2.4.8 (Bouckaert et al., 2014). We used BMODELTEST (Bouckaert & Drummond, 2017)  We reconstructed phylogenomic relationships using SNAPP 1.3 (Bryant, Bouckaert, Felsenstein, Rosenberg, & Roychoudhury, 2012) implemented in BEAST 2.4.8 (Bouckaert et al., 2014). SNAPP is a package that infers species or population trees from unlinked markers such as bi-allelic SNPs, implementing a coalescent model with an algorithm to integrate all possible gene trees rather than explicitly sampling them. To reduce computational requirements and run times, we selected 2-6 representative individuals per population (based on population structure results) with at least 50% complete data matrices. Backwards (u) and forwards (v) mutation rates (expected mutations per site per generation) were estimated from the data in SNAPP, with the birth rate (λ) of the Yule prior (indicating the rate at which populations diverge from one another) based on the number of samples used (Supporting Information Table S6). The run for each data set used a chain length of 1,000,000 generations, sampling every 1,000 trees. We inspected final log files and created maximum clade credibility trees (median node heights) by combining two independent runs per taxon in TREEANNOTATOR 2.4.6 after discarding 10% as burn-in. To verify that the selected individuals did not severely bias SNAPP results, we repeated each analysis using a different random selection of individuals per population (non-overlapping where possible, Supporting Information Table S6).
Given that an accurate mutation rate for amphibians is unavailable, we refrained from inferring the absolute timing of divergence events based on bi-allelic SNP data. Furthermore, as the SNAPP model is coalescent-based, it can account for incomplete lineage sorting, but the presence of high gene flow can cause underestimates of node ages (Leaché, Harris, Rannala, & Yang, 2014a). Additionally, the assembly and filtering of our RAD-seq data sets may also adversely affect dating estimates due to high numbers of retained loci being under selection or linked, which would potentially reduce calculated genetic diversity (Huang & Knowles, 2016). Despite the uncertainties of absolute dating, our SNAPP and δaδi analyses enabled a relative comparison of the divergence time intervals between populations of each species. We investigated divergence estimates for each of these analyses alongside mtDNA estimates and conducted Pearson's correlation tests between the relative divergence time intervals to aid our discussion of forest and generalist taxa.
This program identifies geographic areas where genetic similarity is greater than expected under isolation by distance using spatial and SNP data, without the need for environmental and topographic information. The effective migration and effective diversity estimates are interpolated across geographic space and provide a visual representation of observed genetic dissimilarities, including regions with higher or lower gene flow (i.e., barriers) and effective (i.e., genetic) diversity than expected. We converted filtered STACKS "haplotype" files into PLINK format (.bed files, Purcell et al., 2007) and used the BED2DIFFS program to calculate dissimilarity matrices for each taxon. Based on the size of the habitat and the number of demes required to fill that habitat, we defined a deme size (the density of populations) of 700 in EEMS, but also explored results using smaller deme sizes of 250 and 500 using the SNP version of EEMS (runeems_snps).
We used a MCMC length of 25,000,000 with a burn-in of 1,000,000, each for three replicates. We combined results using the EEMS R plotting (rEEMSplot) package and plotted surfaces of effective migration (m) and effective diversity rates (q).  (Brown, 2014). To create stability models, we calculated the mean of the negative log suitability per epoch, using the exponent of this value to create continuous stability models per taxon ranging between 0 (absent) and 1 (present, Rosauer, Catullo, Vanderwal, & Moussalli, 2015). As the Last Glacial

| Ecological niche and stability models
Maximum data were only available at 2.5 arc-seconds (approximately 5 km 2 ), continuous stability models are resampled to this resolution.
The palaeoclimatic data we used for the ENM projection only cover the period until the Last Interglacial (120 kybp). We therefore use the stability models to reflect recent major climatic events, and as a proxy for deeper time in the absence of accurate palaeoclimate data further back to the Miocene.  Figure S1. Population structure analyses across all taxa revealed that thirteen of the twenty populations were restricted to six allopatric areas which are non-overlapping. The remaining seven populations exhibited wider ranges associated with either the northern or southern parts of the CFEA region (Figure 3).

| Demographic model selection
The consistency of log likelihoods and parameters across replicates of each model increased after the first, second and third optimization rounds and indicated convergence for the best-ranked models across replicates (Supporting Information Table S9, Figure S4). Null models (no divergences) were consistently lowest ranked for 2D models, supporting the genetic distinctiveness of populations. For forest taxa, demographic model selection was inconsistent across taxa (Table 2; Figure 4a), with the best-ranked model selected as historical gene flow for A. xenodactyloides, and allopatric divergence models selected for L. flavomaculatus (historical isolation followed by secondary contact) and A. sylvaticus (divergence and isolation). For all generalist taxa, allopatric models were consistently found as the best-ranked (   Table S10, and full goodness of fit results are shown in Supporting Information Table S11.

| Evolutionary relationships
Mitochondrial DNA-inferred relationships ( Figure 5a) were mostly congruent with those inferred using bi-allelic SNPs from RAD-seq data, despite some discordance in population membership and the evolutionary relationships between populations, which we suspect is

| Effective migration and diversity surfaces
The EEMS analyses based on smaller deme sizes (250, 500, not shown) merged numerous localities that we consider important to F I G U R E 7 Stability models representing areas of persistent suitable habitat across from the Last Interglacial period (120 kya) until the present. (a) Forest taxa. (b) Generalist taxa. Black and white colours represent topography (white = higher elevation, black = lower elevation), and yellow-red colours represent highest stability (ranging between 0 and 1). Rivers are also shown, along with refugia identified by Burgess et al. (1998) shown as red polygons, matching the labelled refugia in Figure 1a  Mulanje massif and surrounding areas).

| Ecological niche and stability models
The best ENMS per taxon had AUC scores between 0.89 and 0.99 (Supporting Information Table S8), with Bio14 (precipitation of driest month) consistently having the highest contribution to the models for six taxa along with input from Bio4 (temperature seasonality) and Bio2 (mean diurnal range). Together, these three variables accounted

| DISCUSSION
At present, our understanding of diversification in the CFEA is far from complete. We used analyses of population structure, demographic model selection, evolutionary relationships and stability models to evaluate diversification patterns and processes with genome-wide data for the first time within this highly threatened biodiversity hotspot. Utilizing such data for multiple co-distributed taxa provides insight into the spatial distribution of biodiversity and the underlying diversification processes (e.g., Bell et al., 2017), which smaller spatial and genetic data sets have previously been unable to address in detail (see Davey & Blaxter, 2010;Lexer et al., 2013).
Using a multi-faceted analytical strategy, our findings suggest that, although biodiversity patterns across the CFEA appear consistent with a forest refugia-driven model of diversification, this hypothesis alone is insufficient to fully explain the region's biological diversity even for forest taxa. We show that multi-taxon approaches can help to develop a more comprehensive understanding of the biotic history of the region.
Most divergence events observed in our focal taxa occurred during the Pliocene-Pleistocene based on mtDNA analysis. Differences BARRATT ET AL.
| 4303 in topologies with SNP data and difficulties in inferring absolute dates without an accurate mutation rate for amphibians make direct comparisons of node ages across these data sets difficult. Despite this, we note that almost all divergences are temporally consistent between SNAPP and δaδi, with time interval parameters (in genetic units) being highly correlated. Although gene flow can have adverse effects on divergence date estimates using Bayesian coalescent models (Leaché, Fujita, Minin, & Bouckaert, 2014b), we found evidence for a lack of gene flow in six of our seven population comparisons, indicating age estimates are not likely to be underestimated from this cause. However, we do note that divergence times may be overestimated when dealing with closely related lineages due to inappropriate models of nucleotide evolution using genome-wide data (Lischer, Excoffier, & Heckel, 2014

| Understanding tropical diversification using multiple taxa
Taxa found across the CFEA, as in other heterogeneous tropical biodiversity hotspots, are a rich mixture of old and young taxa, each with unique ecological characteristics. Inferring diversification processes in tropical biodiversity hotspots is therefore challenging, as even congruent biodiversity patterns are likely to be generated by highly complex processes that vary both temporally and spatially.
With this study, we managed to capture some of these complex processes for forest and generalist taxa and were able to evaluate potential diversification processes at work across taxa using several high-throughput data sets. In doing so, we were able to explicitly test existing hypotheses against alternative diversification modes for the first time in this region and make a broad evaluation of the forest refuge model. Forest refugial processes appear to be only partially responsible for the current diversity in the CFEA, and a range of other diversification mechanisms support the idiosyncrasy of these processes across taxa. Although we found some clear differences between forest specialists and generalist taxa, counter to our predictions the forest specialists were less consistent in diversification mechanisms than generalists.
While the conceptual framework we employ is by no means the only available option to address such hypotheses, approaches such as ours can help to reveal the nuances of diversification which lie behind the correct but rather coarsely defined forest refuge hypothesis. Building upon this work, new statistical tools and techniques for handling high-throughput data, such as minimizing allelic dropout while maximizing amounts of informative data and building more accurate demographic models will be key to leveraging optimal information from RAD-seq data sets in particular. Moreover, newly published model-based methods for genome-wide SNP data Xue & Hickerson, 2017) as well as those already existing for mtDNA (Huang, Takebayashi, Qi, & Hickerson, 2011;Oaks, 2014) will enable further work for assessing concordant patterns of co-diversification amongst co-distributed taxa in tropical hotspots.
These methods allow empirical investigation of the synchronicity of divergences within and between species and populations, detecting if co-distributed taxa exhibit shared responses to major historical events. We strongly encourage future work in tropical regions to take advantage of these emerging tools as well as those we utilized here, and as burgeoning quantities of genome-wide spatial and molecular data become available, we advocate that research collaborations incorporating multiple taxonomic groups will further improve our understanding of the evolutionary processes occurring in biodiversity hotspots.

DATA ACCESSIBILI TY
All raw, unprocessed sequences are deposited in NCBI GenBank (Accession nos MG871749-MG871982 for 16S mtDNA) and the Sequence Read Archive (Accession no. SRP150605 for RAD-seq).
We provide a large package on DRYAD (https://doi.org/10.5061/ dryad.315km76) that includes our final RAD-seq filtered "haplotypes" files, input files for analyses, including mtDNA (also submitted to GenBank). We include data and results for our analyses in this package (ADMIXTURE, DAPC, BEAST, SNAPP, EEMS, δaδi, STACKS and ENMS). All newly created 3D demographic models, along with scripts for performing goodness of fit tests, are freely available in an updated version of the δaδi analysis pipeline found at https://github.com/ dportik/dadi_pipeline.

AUTHOR CONTRI BUTIONS
The