Bayesian divergence-time estimation with genome-wide SNP data of sea catfishes (Ariidae) supports Miocene closure of the Panamanian Isthmus

The closure of the Isthmus of Panama has long been considered to be one of the best defined biogeographic calibration points for molecular divergence-time estimation. However, geological and biological evidence has recently cast doubt on the presumed timing of the initial isthmus closure around 3 Ma but has instead suggested the existence of temporary land bridges as early as the Middle or Late Miocene. The biological evidence supporting these earlier land bridges was based either on only few molecular markers or on concatenation of genome-wide sequence data, an approach that is known to result in potentially misleading branch lengths and divergence times, which could compromise the reliability of this evidence. To allow divergence-time estimation with genomic data using the more appropriate multi-species coalescent model, we here develop a new method combining the SNP-based Bayesian species-tree inference of the software SNAPP with a molecular clock model that can be calibrated with fossil or biogeographic constraints. We validate our approach with simulations and use our method to reanalyze genomic data of Neotropical army ants (Dorylinae) that previously supported divergence times of Central and South American populations before the isthmus closure around 3 Ma. Our reanalysis with the multi-species coalescent model shifts all of these divergence times to ages younger than 3 Ma, suggesting that the older estimates supporting the earlier existence of temporary land bridges were artifacts resulting at least partially from the use of concatenation. We then apply our method to a new RAD-sequencing data set of Neotropical sea catfishes (Ariidae) and calibrate their species tree with extensive information from the fossil record. We identify a series of divergences between groups of Caribbean and Pacific sea catfishes around 10 Ma, indicating that processes related to the emergence of the isthmus led to vicariant speciation already in the Late Miocene, millions of years before the final isthmus closure.

with a molecular clock model that can be calibrated with fossil or biogeographic 148 constraints. SNAPP is well suited for analyses of genome-wide data as it infers the species 149 tree directly from single-nucleotide polymorphisms (SNPs), through integration over all absolute divergence times when properly calibrated with fossil or biogeographic evidence. 163 We evaluate the accuracy and precision of our approach using an extensive set of 164 simulations, and we compare it to divergence-time estimation based on concatenation. We  Table 1. Based on the results of experiments 1-5, 188 we developed recommendations for divergence-time estimation with SNP data, and we then 189 applied this approach to infer timelines of evolution for Neotropical army ants and sea  Table S1). All simulated data sets were based 196 on the same set of 100 species trees generated with the pure-birth Yule process (Yule 1925) 197  for all branches. These population sizes were set to N = 25 000 diploid individuals for most 208 analyses, but we also used the larger population sizes N = 100 000 and N = 400 000 in the 209 simulations conducted for experiment 2 (Table 1). For each simulated gene tree, between 2, 210 4, or 8 terminal lineages were sampled per species, corresponding to 1, 2, or 4 diploid individuals per species (Table 1)  resulting data sets of close to 10 000 unlinked SNPs were further subsampled randomly to 228 generate sets of 300, 1 000, and 3 000 bi-allelic SNPs for each species tree (see Table 1).

229
For the analyses in experiments 1-4, each of the 100 data sets of 300, 1 000, and 230 3 000 SNPs was translated into the format required for SNAPP, where heterozygous sites 231 are coded with "1" and homozyguous sites are coded as "0" and "2". Per site, the codes "0" 232 and "2" were randomly assigned to one of the two alleles to ensure that the frequencies of 233 these codes were nearly identical in each data set. For experiment 4 in which we tested for 234 the effect of ascertainment bias in SNAPP analyses, the data sets of 1 000 SNPs were also 235 modified by adding invariant sites. To each set of 1 000 SNPs, between 12 184 and 32 740 236 invariant sites (alternating "0" and "2") were added so that the proportion of SNPs in 237 these data sets matched the mean proportion of variable sites in the alignments initially 238 generated for the respective species tree. Finally, for analyses using concatenation in 239 experiment 5, we added the same numbers of invariant sites to the data sets of 1 000 SNPs; 240 however, in this case we used the untranslated versions of these data sets with the original 241 nucleotide code, and also used nucleotide code for the added invariant sites (randomly 242 selecting "A" , "C", "G", or "T" at each site). This assumption was met in our simulated data sets but may often be violated by  Finally, as we were interested in SNAPP's ability to infer divergence times rather 284 than the species-tree topology (which has been demonstrated previously; Bryant et al.

285
2012), we fixed the species-tree topology to the true topology. We provide a script written 286 in Ruby, "snapp_prep.rb", to generate XML input files for SNAPP corresponding to the 287 settings described above (with or without a fixed species tree). Note that these settings, 288 including the use of scale-invariant prior distributions, were deliberately not tailored 289 towards our simulated data sets, but were instead intended to be generally applicable for 290 divergence-time estimation with any SNP data set. As a result, the XML files produced by 291 our script should be suitable for analysis without requiring further adjustments from the 292 user. Our script is freely available at https://github.com/mmatschiner/snapp_prep.

293
Details on operators used in our analyses are provided in Supplementary Text S1.

294
As SNAPP is specifically designed for the analysis of bi-allelic SNPs, its algorithm for or invariant sites were added to data sets of 1 000 SNPs (see Table 1). This option did 301 not apply to the analyses of concatenated data in experiment 5 as these were not  analyses of 100 data sets of 300, 1 000, and 3 000 SNPs with node-age constraints on either 317 the root or a younger node, is shown in Figure 1 and summarized in containing the true node age was always slightly higher in analyses with root-node 324 constraints even though the width of these HPD intervals was generally smaller. Notes: Accuracy was measured as the percentage of 95% HPD intervals containing the true node age. Precision was measured as the mean width of 95% HPD intervals for node-age estimates. Both measures are presented separately for young (true node age < 10 myr) and old (true node age > 10 myr) nodes. Ex.   Results are based on 100 species trees and 300 to 3 000 SNPs generated per species tree. a) Node ages estimated with an age constraint on the root. b) Node ages estimated with an age constraint on a node that is approximately a third as old as the root. Mean age estimates of constrained and unconstrained nodes are marked with red and gray circles, respectively, and vertical bars indicate 95% HPD intervals.   Table 3). While both parameters were 355 underestimated roughly by a factor of three when ascertainment bias was corrected for, 356 leaving this bias unaccounted led to parameter overestimation by more than an order of 357 magnitude. Importantly, however, when ascertainment bias was accounted for, the Results are based on data sets of 1 000 SNPs generated for each of 100 species trees, analyzed with and without SNAPP's ascertainment-bias correction or after adding invariant sites to the data sets. Gray circles indicate mean estimates and 95% HPD intervals are marked with vertical bars. The visualization of node-age estimates in a) is equivalent to the illustration in Fig. 1, except that only unconstrained nodes are shown. Note that logarithmic scales are used for estimates of the clock rate (a) and Θ (b).

325
resulting estimates of the population size N (calculated as N = Θ/4µ with µ being the mutation rate per generation, i.e., the estimated clock rate divided by the number of 360 generations per myr) accurately recovered the true population size used for simulations 361 (N = 25 000 in all simulations conducted for experiment 4; see Table 1), as 95% of the 95% 362 HPD intervals included the true parameter value (Fig. 2c, Table 3). In contrast, the 363 population size was underestimated when ascertainment bias was not corrected for: Mean 364 estimates were on average 17.4% lower than the true population size and 35% of the 95% 365 HPD intervals did not include the true parameter value (Fig. 2c, Table 3).

366
Our results of experiment 4 also showed that when invariant sites were excluded,  True node age (myr) Results are based on analyses of 100 data sets of 1 000 SNPs, simulated with population sizes N = 25 000, N = 100 000, and N = 400 000. Gray and red dots indicate node-age estimates obtained with the MSC implemented in SNAPP and with BEAST analyses of concatenated data, respectively. Node-age error is measured as the ratio of the estimated node age over the true node age. Solid lines represent mean node-age errors in bins of 0.2 myr. Only nodes with true ages up to 10 myr are shown to highlight differences between the two methods. Note that a logarithmic scale is used for node-age error.
recovered reliably in these analyses and were included in 97.2% of the 95% HPD intervals 378 (Fig. 2d).  Notes: Mean node-age error is presented separately for nodes with young (true node age < 10 myr) and old (true node age > 10 myr) nodes. Note that very young nodes (true node age < 0.5 myr) are excluded from this comparison.
concatenation, the mean age estimate for a node with a true age around 3 Ma (±0.2 myr)  Fresh fin tissues were preserved in 96% ethanol for subsequent DNA extraction.
following the manufacturer's instructions. RNase treatment after digestion (but before precipitation) was performed in order to improve the purity of the samples. DNA   All fossils used for phylogenetic analyses are summarized in Supplementary Table S6.   almost exclusively from the width of the calibration density (Fig. 1). In addition to data-set 600 size, the placement of the node-age calibration also had an effect on the precision of 601 divergence-time estimates, which was improved when the root node was calibrated instead 602 of a younger node. This suggests that future studies employing divergence-time estimation with SNAPP should make use of constraints on the root node if these are available from 604 the fossil record, from biogeographic scenarios, or from previously published time-calibrated 605 phylogenies (as in our analyses of empirical SNP data of Neotropical army ants and sea 606 catfishes). While we did not test the performance of multiple calibration points with 607 simulated data, the use of additional calibration points can be expected to further improve 608 the precision of divergence-time estimates; therefore these should be used if available.

609
It should be noted that even though all our analyses of both simulated and 610 empirical data sets were calibrated through node-age constraints, this so-called "node implemented. This means that particularly in clades that may be expected to have 644 different mutation rates in different lineages, the precision of divergence-time estimates 645 may be exaggerated, which should be considered in the interpretation of such results.

646
Our experiment 4 revealed that when SNP data sets are used without the addition 647 of invariant sites, SNAPP's estimates for the clock rate and Θ did not match those used in 648 simulations (Fig. 2a,b, Table 3). While this mismatch might appear as a weakness of our 649 approach, we do not consider it unexpected that these estimates change when