Incorporating Topological and Age Uncertainty into Event-Based Biogeography of Sand Spiders Supports Paleo-Islands in Galapagos and Ancient Connections among Neotropical Dry Forests

Event-based biogeographic methods, such as dispersal-extinction-cladogenesis, have become increasingly popular for attempting to reconstruct the biogeographic history of organisms. Such methods employ distributional data of sampled species and a dated phylogenetic tree to estimate ancestral distribution ranges. Because the input tree is often a single consensus tree, uncertainty in topology and age estimates are rarely accounted for, even when they may affect the outcome of biogeographic estimates. Even when such uncertainties are taken into account for estimates of ancestral ranges, they are usually ignored when researchers compare competing biogeographic hypotheses. We explore the effect of incorporating this uncertainty in a biogeographic analysis of the 21 species of sand spiders (Sicariidae: Sicarius) from Neotropical xeric biomes, based on a total-evidence phylogeny including a complete sampling of the genus. Using a custom R script, we account for uncertainty in ages and topology by estimating ancestral ranges over a sample of trees from the posterior distribution of a Bayesian analysis, and for uncertainty in biogeographic estimates by using stochastic maps. This approach allows for counting biogeographic events such as dispersal among areas, counting lineages through time per area, and testing biogeographic hypotheses, while not overestimating the confidence in a single topology. Including uncertainty in ages indicates that Sicarius dispersed to the Galapagos Islands when the archipelago was formed by paleo-islands that are now submerged; model comparison strongly favors a scenario where dispersal took place before the current islands emerged. We also investigated past connections among currently disjunct Neotropical dry forests; failing to account for topological uncertainty underestimates possible connections among the Caatinga and Andean dry forests in favor of connections among Caatinga and Caribbean + Mesoamerican dry forests. Additionally, we find that biogeographic models including a founder-event speciation parameter (“+J”) are more prone to suffer from the overconfidence effects of estimating ancestral ranges using a single topology. This effect is alleviated by incorporating topological and age uncertainty while estimating stochastic maps, increasing the similarity in the inference of biogeographic events between models with or without a founder-event speciation parameter. We argue that incorporating phylogenetic uncertainty in biogeographic hypothesis-testing is valuable and should be a commonplace approach in the presence of rogue taxa or wide confidence intervals in age estimates, and especially when using models including founder-event speciation.


Introduction
The reconstruction of the biogeographic history of organisms is one of the main aims of systematic biology. Based on a known phylogeny, researchers may attempt to glimpse into the ancestral area where a particular clade originated [1], infer the number and direction of dispersal events [2], or estimate the number of vicariant events along the evolutionary history of the group [3]. More often than not, the geographic distribution of a group is known only from its extant species, occasionally accompanied by a few fossil forms. This represents but a fraction of the diversity of a group throughout its history, and we cannot directly assess the distribution range of unknown extinct species. Thus, biogeographers interested in such questions must resort to methods that attempt to estimate past distribution ranges from the data available from known species.
Initially, such estimates of ancestral geographical ranges relied on algorithms such as Fitch or Camin-Sokal parsimony optimization [1,4]. This approach treats distribution ranges as discrete characters and optimizes them as such along the phylogenies, but has some drawbacks. For instance, only tips of the phylogeny may be 'polymorphic' and occur in more than one area simultaneously, and important biogeographic processes such as vicariance are not modeled at all. This changed with the advent of event-based biogeography methods, pioneered by dispersal-vicariance analysis (DIVA) [3]. Such methods attempt to explain the current distribution of organisms by modeling biogeographic events such as dispersal (range expansion into a new area), range contractions (extinction of a species in a particular area) and allopatry (allopatric speciation leading to each descendant inheriting only part of the range of the ancestor species). DIVA brought a fundamental advance with respect to previous methods: changes in a lineage's geographic distribution may happen as anagenetic events (dispersal or extinction along branches) or cladogenetic events (allopatry at nodes). Furthermore, the states in such estimates are geographic ranges, which may consist of one or more areas; thus, both tips and ancestors may be 'polymorphic'. By using a parsimony framework, DIVA assigns costs to dispersal and extinction events, and treats allopatry as the null expectation for explaining biogeographic history. A similar rationale was used later to elaborate a likelihood-based approach to ancestral range estimation named dispersal-extinction-cladogenesis analysis (DEC) [5,6]. DIVA's parsimony-based optimization is agnostic to branch lengths. On the other hand, DEC models anagenetic events as a time-continuous function that takes branch lengths into account; thus, longer branches are more likely to contain anagenetic events such as dispersal or extinction. More importantly, DEC can accommodate dated trees, and thus prior information on geological history can be utilized for explicit testing of biogeographic hypothesis [7]. Modifications of DIVA and DEC have been implemented and further elaborated in software packages for biogeography, such as RASP [8] and BioGeoBEARS [9,10]. The latter has become particularly popular due to its flexible implementation of biogeographic models allowing for different types of cladogenetic events, such as allopatry, subset sympatry, and founderevent speciation [10,11]; in addition, all models are implemented in a likelihood framework, allowing direct model comparison. It also allows incorporating prior information, such as dispersal probabilities among areas, to put biogeographic hypotheses to test in an explicit manner (e.g., [12]).
A unifying feature of all methods discussed above is that they rely on the knowledge of the geographic distribution of each of the taxa, and of the phylogenetic tree describing their interrelationships. As in any comparative method, estimates of ancestral ranges are only as reliable as the primary data underlying it. The knowledge of the geographic distribution of tips may be affected by the Linnean and Wallacean shortfalls [13], i.e., gaps in the data about existing species and their distributions. However, it is arguable that estimates of ancestral ranges are usually carried out by systematists that specialize in a particular taxon, who strive to include all or most known species and distribution records in their sampling, so as to mitigate any negative effects of such shortfalls. On the other hand, the imperfect knowledge of the underlying phylogenetic tree can be potentially more problematic. In recent years, the accumulation of genomic-scale studies has shown that the phylogenetic relationships of some clades are elusive even with massive quantities of data due to e.g., incomplete lineage sorting and/or very short internodes (e.g., [14,15]), and some portions of such trees are shrouded in topological uncertainty. There is uncertainty not only regarding the topology, but also in the estimates of divergence times. This uncertainty in age estimates is also resistant to the accumulation of genetic markers (see [16]). Rates of molecular evolution rarely conform to a strict molecular clock, and branch lengths estimated from molecular sequences are a product of substitution rate and time, such that each of these two parameters are individually non-identifiable [17]. In addition, estimation of dated trees builds upon several assumptions, some of which may not be met by empirical datasets [18], and require fossil or other external calibrations, whose selection and justification is not free of difficulties [19,20]. In this scenario, analyses that depend on phylogenetic trees might benefit from accounting for any uncertainties regarding their topologies and/or ages [21]. In this paper, we refer to "phylogenetic uncertainty" as any uncertainty regarding the topology and/or node ages.
Efforts have been made to incorporate such phylogenetic uncertainty into ancestral range estimation. Nylander et al. [22] used multiple trees from the posterior distribution of a Bayesian analysis to estimate ancestral ranges, and thus infer the biogeographic history of a clade of birds whose phylogeny had proven difficult to solve; they called this approach Bayes-DIVA. Similar pipelines have been incorporated into RASP (S-DIVA; [23]) and BioGeoBEARS (using the run_bears_optim_on_multiple_trees function). These functions summarize estimates of ancestral ranges of several trees on a target tree, such as a majority-rule consensus (as output by e.g., MrBayes) or a maximum clade credibility tree (MCC tree; as output by e.g., BEAST) and are routinely employed by researchers interested in incorporating phylogenetic uncertainty (e.g., [24,25]). However, the results are summarized in a single tree, and the underlying variability in the estimates is difficult to grasp. A more elaborate approach is the joint estimation of the phylogeny and ancestral ranges (e.g., [26,27]), although it can be computationally intensive in complex systems. Furthermore, the individual underlying trees usually are not taken into account when comparing competing biogeographic hypotheses regarding particular events. We argue that it is important to understand the effect of analyzing multiple trees during hypothesis-testing. This is especially relevant in time-stratified analyses, where different time periods can have different biogeographic possibilities (e.g., different availability of areas or dispersal probabilities), as confidence intervals of age estimates may cross boundaries of time slices.
In addition to phylogenetic uncertainty, there is the uncertainty associated with stochastic time-continuous models such as DEC. Because of this, estimates for ancestral nodes might include several possible states, each with its own probability. This hampers counting biogeographic events or estimating their ages. Dupin et al. [2] solved this by introducing biogeographic stochastic mapping (BSM), a method to count biogeographic events while accounting for uncertainty in ancestral range estimation. Several replicates are run, and in each one the states at each node are resolved by taking into account the probabilities of each state. This allows, for instance, counting the number of dispersal events among areas, or estimating the time when such transitions took place. This approach is elegant, but it is usually employed on estimates based on a single tree, and thus incorporates uncertainty on ancestral range estimates conditional on a single topology. If there is topological uncertainty leading to different ancestral range estimates, this approach will underestimate the uncertainty in the biogeographic reconstruction ( Figure 1). Furthermore, using a single tree as input ignores the confidence intervals in age estimates, which are usually large. We argue that it would be productive to run stochastic maps over a sample of trees to incorporate both phylogenetic and stochastic uncertainty simultaneously. This approach has been successfully implemented to estimate dispersal rates among areas through time in some studies (e.g., [27][28][29]). Biogeographic stochastic maps are insufficient to account for uncertainty in ancestral range estimates. Each column represents a single tree from the posterior distribution of the Bayesian analysis of our dataset, with the most likely estimated states in the top row and different stochastic maps in the bottom rows. While estimates in each stochastic map are different, the most remarkable differences are found among trees with different topologies. Note the tip marked with a grey star, which is a rogue taxon (Sicarius andinus). As it shifts its position across the different trees of the stationary phase of the chain, it leads to substantially different ancestral range estimates in each tree.
To illustrate our argument, we explore the effect of uncertainty in the biogeographic history of Neotropical sand spiders (Sicarius). These spiders represent an ideal system for this test because they are moderately diverse (21 species) and all species have been included in a dated total-evidence phylogeny [30]. The genus has a disjunct distribution, and each species is restricted to one or two arid areas surrounded by mesic habitats; phylogeographic and phylogenetic patterns suggest they are very poor dispersers [31][32][33]. Thus, their areas of distribution are clearly delimited and could be interpreted as "islands" of dry biomes inserted in a matrix of unsuitable humid habitats. Because of their distribution and moderate diversity, most of the biogeographic transitions should be straightforward to interpret, thus allowing us to easily measure the effects of phylogenetic uncertainty in biogeographic estimates.
Specifically, we focus on two particularly pressing questions. The first is the timing of arrival of Sicarius in the Galapagos archipelago. The Galapagos are currently inhabited by a single species of sand spider, Sicarius utriformis (Butler), which is sister to S. peruensis (Keyserling) from the Peruvian coastal deserts [30,33]. Although the oldest emerged islands are~3.5 million years (Myr) old [34], geological evidence suggests that the archipelago existed for at least 14.5 Myr, when it was composed of paleo-islands that are now submerged [35,36]. The age of divergence between S. utriformis and S. peruensis has a 95% confidence interval between 1.2 and 22.2 Myr (median 9.7; [30]), and thus it is possible that this pair of species split before the current islands were formed, but during the time the archipelago was formed by paleo-islands. This makes this system ideal to test the effect of the uncertainty of age estimates in biogeographic inference.
The second question is the estimation of the ancestral distribution range of a clade endemic to the Brazilian Caatinga. This area is one of the largest and most diverse tropical dry forests in the world [37], and six Sicarius species inhabit the region [33]. These six species form a well-supported monophyletic group [30], suggesting a single dispersal into this area. In this case, from which area did the ancestor of this clade come? The American tropical dry forests and other xeric biomes such as deserts and scrublands currently have a disjunct distribution, but the similarity in their biota suggests they have been connected in the past (see [37]). Based on plant distributions, some argued that the Caatinga could have been connected to dry forests in Bolivia and Argentina, or to dry forests in the Caribbean coast of northern South America [38,39]; these connections would have taken place by expansion of dry forests over areas that now are covered by mesic biomes. Connections between the Caatinga and the Caribbean dry forests imply a northern route passing through present-day Amazon (a rainforest); alternatively, the Caatinga could have been connected with southern formations, such as the Monte or the Chiquitano dry forests, passing through present-day Cerrado (a savannah). We thus aim to identify the most likely route for the occupation of the Caatinga by sand spiders. However, this is hampered by the fact that one species, Sicarius andinus Magalhaes et al. from the Peruvian Andes, is a rogue taxon in the phylogeny and does not have a well-resolved phylogenetic position [30]. Different positions of this species may yield different biogeographic estimates for the Caatinga clade, and thus we must take this into account.
In this paper, we test the effect of taking phylogenetic uncertainty into account to address the two biogeographic questions mentioned above. We combine processing of several trees of the stationary phase of the Markov chains of a Bayesian analysis with biogeographic stochastic maps for each tree, so that both phylogenetic and biogeographic uncertainties are considered simultaneously. Specifically, we test (1) whether data from sand spiders support dispersal to Galapagos in the last 3.5 Myr (age of oldest emerged island), in the last 14.5 Myr (age of oldest known submerged paleo-island) or an unconstrained model (representing the possibility of older, yet-undetected paleo-islands), and (2) whether the Caatinga was connected to northern (Caribbean, Mesoamerican or Andean dry forests) or southern (Chiquitano dry forests or Monte) biomes, as well as the age of such connections. We anticipate that analyzing a sample of trees, instead of a single target tree, provides invaluable insights for testing competing biogeographic hypotheses. Finally, we provide scripts for R [40] to replicate the analyses described below, in the hope they will be useful for further studies.

Summaries of Biogeographic Inferences in Face of Uncertainty
We prepared an R script to (1) sample n trees randomly from BEAST output files and prune them to the taxa of interest, (2) estimate ancestral ranges using BioGeoBEARS for each of these trees, (3) run biogeographic stochastic maps for each of these estimates, and finally (4) parse the results and summarize them. The summaries include: (1) maximum likelihood parameter estimates, log-likelihoods and AICc scores for ancestral range estimates of each tree, (2) tables containing the transitions between geographic ranges for each stochastic map of each of the n sampled trees, along with the age of such transitions; (3) a graphic of lineages through time per area averaged over all trees and stochastic maps; (4) tables with all possible geographic ranges in both rows and columns, and average counts of how many times a transition between a particular pair of ranges took place; such tables are broken down by each type of transition (dispersal, extinction, allopatry, etc.) and summarized in a table containing all types of transitions; and (5) a summary of the most common biogeographic transitions (as the mean number of transitions per BSM replicate). The script (Online Supplementary File S1) uses functions from the packages phytools [41], ape [42], sjmisc [43] and BioGeoBEARS [9] and has been written as to be easily adapted to most datasets. The necessary input files are the same as those used in BioGeoBEARS, except that users may provide multiple trees instead of a single target tree.

Model Selection
BioGeoBEARS implements three different biogeographic models that differ mainly in the events that may take place during cladogenesis when the ancestor has a widespread range (i.e., its range consists of two or more areas; see [10] for a summary). DEC has an identical implementation to the model described in [6] and allows narrow vicariance (one descendant inherits exactly a single area, the other inherits the rest of the range) or subset sympatry (one descendant inherits exactly a single area, while the other inherits the whole distribution range of the ancestor). DIVA-like allows the same cladogenetic events as the original implementation of DIVA [3]: narrow allopatry (as in DEC) and wide allopatry (each descendant may inherit two or more areas from the ancestor). BAYAREA-like does not allow changes in distribution ranges during cladogenesis, and each descendant inherits exactly the same range as the ancestor. Each of these three models can be modified by the addition of a "jump-dispersal" free parameter ("+J") that allows for founder-event speciation during cladogenesis, i.e., one of the descendants occupies a single area that is not part of the range of the ancestor, while the other descendant inherits the same range as the ancestor [11]. Thus, models including founder-event speciation are unique in that they allow dispersal to take place during cladogenetic events. Please note that founder-event speciation is also called jump dispersal, which leads some researchers to interpret it as long-distance dispersal; it should be noted that this notion is incorrect, as this parameter merely models dispersal taking place instantaneously during cladogenetic events (hence the "jump"). It is affected by dispersal matrices in the same way as anagenetic dispersal is, and thus the distances among areas are not more important for jump dispersal than they are to other parameters.
Biogeographic history may be estimated under each of these six models and the fit of the data to each of them can be compared by using the Akaike information criterion (AIC; [44]) and Akaike weights (AICw; [45]). This procedure has been used to guide the selection of models that best explain the data and for testing biogeographic hypotheses [10,11]. We here estimate the fit of the data to six different models (DEC, DIVAlike, BAYAREA-like, DEC + J, DIVA-like + J, BAYAREA-like + J) for initial exploration of the behavior of the estimates. The script used for model selection can be found as Online Supplementary File S5.
We wanted to investigate whether the choice of a particular biogeographic model interacts in any way with the decision to run analyses over a sample of posterior trees. Thus, we ran the aforementioned script under two models with (DIVA-like + J and BAYAREA-like + J) and without (DIVA-like and DEC) founder-event speciation. The choice of these four models was made to include those that are a better fit to the data (DIVA-like and the +J variant), a model that is frequently used in empirical studies (DEC) and a model excluding the possibility of allopatry, and thus more similar to simple parsimony reconstructions (BAYAREA-like + J). Each of these four models was used in both unconstrained and timestratified analysis (see below), and using either the MCC tree or 100 posterior trees as source. In this latter case, to make these analyses directly comparable, the same sample of 100 trees was used for all runs. The combination of four biogeographic models, three time-stratification scenarios and two sources of trees resulted in 24 different runs (see Online Supplementary Figure S10). For each individual tree, 100 biogeographic stochastic maps were estimated to count biogeographic events and estimate the number of lineages through time by area (see below).

Phylogeny and Distribution Data
To reconstruct the biogeographic history of Sicarius, we used a recently published phylogeny estimated using morphology and DNA sequences and dated using a combination of fossil calibrations and substitution rates for the histone H3, subunit A gene [30]. To understand the effect of incorporating uncertainty into biogeographic inference, analyses were run both on the maximum clade credibility (MCC) tree, and on a sample of 100 trees randomly drawn from the posterior distribution of the analysis, after removing the first 10% samples as burn-in (see above). Hexophthalma is the African sister group of Sicarius and was included in the analyses; the remaining terminals (Loxosceles and non-sicariid outgroups) were pruned from the trees. Species geographic ranges and trees can be found as Online Supplementary Files S2-S4.
The distribution of each species has been fully mapped in recent taxonomic publications on the genus including ca. 1800 adult specimens from natural history collections and recent field expeditions [33,46]. We classified species distribution in ten areas: southern Africa deserts and xeric scrublands (F), Argentinean Monte (O), Atacama Desert and neighboring Chilean xeric scrublands (T), Sechura desert in the Peruvian coast (S), Andean dry forests (D), Chiquitano dry forests in Bolivia (C), Mesoamerican dry forests (M), dry forests in the Caribbean coast of Colombia (B), Caatinga dry forest in Brazil (I), and the Galapagos Islands (G) (Figure 2). Most of these areas are clearly delimited by geographic barriers such as oceans or mountain ranges, or correspond to well-recognized ecoregions or phytogeographic units [37,47].

BioGeoBEARS Parameters and Time Stratification
Sand spiders are poor dispersers [30] and most species are restricted to a single area, with only two species occurring in two areas [33]. For this reason, and to speed up calculations, we restricted the maximum range size to include three areas. Likelihood calculations were carried out with optimx [48]. We compared ancestral ranges estimates under three scenarios: (1) unconstrained, allowing dispersal to the islands in any time, and (2 & 3) two different time-stratified scenarios, each with two time slices, where the Galapagos Islands were only available for occupation in the more recent slice. The boundary between the two slices was set to (2) 15 Myr, as geological evidence points to submerged islands that are at least 14.5 Myr old [36], or (3) 3.5 Myr, representing the approximate age of the oldest emerged island [34]. Files for implementing the time-stratified model are available as Online Supplementary Files S6-S8.

Model Selection and Estimates of Ancestral Ranges
Regarding runs on MCC trees, overall DEC and DIVA-like resulted in similar ancestral range estimates among them, as did all models including a +J parameter. Parameter estimates and fit of data under different models are summarized in Online Supplementary  Table S1. AIC values indicate that DIVA-like (log-likelihood: −51.29, AIC: 107. 16) is the favored model among those not including a founder-speciation free parameter, while DIVA-like + J (log-likelihood: −44.24, AIC: 95.69) is favored among models including this parameter. Because it has been demonstrated that models including founder-event speciation are prone to over-fitting [5]; but see [49], we here show results of the ancestral range estimates for the MCC tree under DIVA-like ( Figure 2); results under DIVA-like + J can be seen in Online Supplementary Figure S11.
We briefly investigated the effect of taking topological uncertainty into account during model selection by comparing the AICc values of the estimates of each of the 100 trees for the models DIVA-like, DIVA-like + J and DEC. A histogram of such values displays some overlap among values of the different models (Online Supplementary Figure S12). However, for each individual tree the relationship of the preferred model is maintained (i.e., DIVA-like + J is preferred to DIVA-like, which is preferred to DEC), and is identical to the order found by running the analysis using the MCC tree. Thus, taking uncertainty into account has not affected the results of model selection.

Lineages through Time by Area
We estimated the number of lineages occupying each individual area through time. For this, the tree has been divided in 1 Myr slices, and we counted the number of species in each area in each slice; widespread species are counted once in each area of their distribution range. The number of lineages in each area increases with time ( Figure 3) as a combination of within-area speciation and new dispersals into the area. When using a single tree, changes in diversity associated with cladogenetic events (e.g., within-area speciation or founder-event speciation) appear as abrupt increases in the plot (Figure 3a,b). This effect is stronger in the model including founder-event speciation (Figure 3a), but disappears when topological and age uncertainty is taken into account (Figure 3c,d). It should be noted that now extinct lineages do not contribute to the plots, so the graphs only portray minimal (not actual) regional species richness over time (see Discussion).
The results (Figure 3) indicate that the ancestor of Hexophthalma + Sicarius lived in a range composed of (1) southern African deserts and scrublands, and (2) either the Atacama desert, Sechura desert, or Argentinean Monte. These latter three have similar probabilities of being part of the ancestral range of sand spiders due to uncertainty in topology, divergence times, and biogeographic estimates. There has been a slow and steady increase in diversity in these temperate desert areas for the last~100-80 Myr. On the other hand, lineages only started occupying tropical dry forests later, around 50 Myr. The Caatinga has become the most species-rich region relatively rapidly due to within-area speciation, mainly during the Miocene.  ) show abrupt changes in the number of lineages that are related to cladogenetic events, whose age does not vary because there is no uncertainty associated to the tree. This effect is more pronounced in the model with founder-event speciation (a) because it relies more on cladogenetic events for estimating ancestral ranges. Using several trees (each with 100 stochastic maps) to account for topological and age uncertainty removes this effect and smooths the curves (c,d), reducing differences between models with and without founder-event speciation.

Comparing Time-Stratified vs. Unconstrained Models in the Face of Age Uncertainty
We compared the fit of the data to unconstrained models (occupation of Galapagos possible at any time) to time-stratified models (occupation of Galapagos possible only in the last 15 Myr, or in the last 3.5 Myr). First, we investigated the inferred age of dispersal to the islands in the unconstrained analysis. In the analyses using a MCC tree as input, dispersal to the Galapagos has been inferred to occur more recently than 15 Myr in 72% (DIVA-like; mean 13.8 Myr) or 97% (DIVA-like + J; mean 8.8 Myr) of the stochastic maps (Figure 4a,b). Using a MCC tree results in sharper, overconfident distributions of the inferred ages of dispersal to the islands. This is especially notable in the case of the model with founderevent speciation, where 88% of the replicates inferred that dispersal into the islands is a cladogenetic jump-dispersal with age equal to that of the S. utriformis-S.peruensis node in the MCC tree (Figure 4a). In analyses using 100 posterior trees as input, dispersal to the Galapagos has been inferred to occur more recently than 15 Myr in 57.8% (DIVA-like; mean 16 Myr) or 78.2% (DIVA-like + J; mean 12 Myr) of the replicates (Figure 4c,d). Thus, when using 100 posterior trees, inferred ages of the dispersal to Galapagos are older in average. In addition, the distribution of inferred ages is flatter, as it takes into account the uncertainty in node ages; this reduces the difference between models with and without founder-event speciation. Figure 4. Histograms with the inferred age of dispersal of Sicarius to the Galapagos Islands. Analyses based on a maximum clade credibility tree (MCC tree) (a,b) produce sharper estimates that are closely tied to the age of split of Sicarius utriformis (from Galapagos) and its sister species. This effect is especially pronounced in the model including founder-event speciation (a), where 88% of the estimates have exactly the same age of that split. Accounting for uncertainty in topology and age estimates (c,d) reveals that the age of dispersal is much more uncertain. The dashed vertical line indicates the boundary between time slices in the time-stratified analysis and corresponds to the geological evidence of the oldest submerged paleo-islands of the Galapagos archipelago.
We then compared the fit of the data to unconstrained model vs. time-stratified model allowing dispersal only in the last 15 Myr. Using the MCC tree, the data fit better to a stratified model under DIVA-like, DIVA-like + J, and BAYAREA-like + J, while it fits better to an unconstrained model under DEC (stars in Figure 5). In all cases, however, the support for the preferred model is very weak, as the ratio between AICc weights of the preferred model ranges only between 1.29 to 2.58, and thus we cannot decisively reject any of the two models in favor of the other. When comparing models over 100 posterior trees, we observed an interesting pattern. Again, in most trees the data fit slightly better a time-stratified scenario under DIVA-like, DIVA-like + J, and BAYAREA-like + J (in 58, 68, and 58 of the 100 trees, respectively), and an unconstrained scenario under DEC (in 91 of the trees), but relative supports in these cases are once again ambiguous (1.67-2.91). However, in the fewer cases of trees where the data fit better an unconstrained model, the median relative support is higher and shows decisive support for this model (5.00-151.74) ( Figure 5). We suspected that this pattern might be due to the age of split between S. utriformis-S. peruensis. The confidence interval of this age spans a wide range (1.2-22.2 Myr) and actually crosses the boundary between the two time slices of the time-stratified model (15 Myr). To investigate this, we plotted the AICc weight of the unconstrained model against the age of the split ( Figure 5). The plot indicates that when the age of the split is younger than 15 Myr, there is weak support for the time-stratified model (DIVA-like and DIVA-like + J), to the unconstrained model (DEC), or the support to either of them is ambiguous (BAYAREA-like + J). On the other hand, the data fit better to the unconstrained model in all the 17 trees where the S. utriformis-S. peruensis split is older than the boundary of the time slice (15 Myr). While this is true for the four biogeographic models employed, the effect is much stronger in models including a founder-event speciation parameter, which tend to favor the unconstrained model more strongly: the relative support for the unconstrained model is higher in DIVA-like + J (776.2 ± 369) and BAYAREA-like + J (85.7 ± 110.86) than in DIVA-like (37.81 ± 37.97) and DEC (8.84 ± 4.73). Thus, in at least some of the trees from the posterior distribution, there is strong support for an unconstrained colonization of the Galapagos taking place before 15 Myr.
Finally, we compared the fit of the data to a model allowing dispersal to the Galapagos in the last 15 Myr (age of the oldest recorded paleo-islands) to a model allowing dispersal only in the last 3.5 Myr (age of the oldest emerged island). Under all biogeographic models, the model allowing dispersal in the last 15 Myr is strongly favored (Figure 6). In the few trees where this split is younger than 3.5 Myr, the 15 Myr model is still strongly favored by DIVA-like and DEC, while support to either scenario is ambiguous in DIVA-like + J and BAYAREA-like + J. Thus, our data indicate that dispersal of sand spiders to the Galapagos took place before the appearance of the oldest emerged island.

Ancestral Range Estimates of Particular Nodes in the Face of Topological Uncertainty
We estimated the most likely state in the root node (corresponding to the split between African Hexophthalma and Neotropical Sicarius). For simplicity, we only report the results under DIVA-like, which are similar to those of other biogeographic models explored here. Analyses using the MCC tree yield estimates for an ancestral range most likely including southern African scrublands and one of the temperate American deserts (Sechura, Atacama or Monte): most likely ranges are FOT (30.7%), FOS (30.3%), FTS (24.6%), FO (4.8%), FT (3.8%) and FS (3.8%), summing to a total of 98.3%. The most likely ranges across the 100 trees are FOT (27.8%), FOS (24.4%), FTS (22.4%), FO (5.5%), FT (4.2%) and FS (3.7%), total 88.1%. These latter ranges have their likelihoods slightly diminished because of an increase in the likelihood of ranges including Andean dry forests, namely FOD (2.3%), FSD (2%), and FTD (1.5%). Nonetheless, these changes are rather small, and, even in the face of phylogenetic uncertainty, we can be fairly confident that the ancestor of Sicarius + Hexophthalma lived in a range including deserts and xeric scrublands of southern Africa and southern or western South America.
We used stochastic maps to estimate the ancestral range that sourced species to the Brazilian Caatinga. In more than 99.9% of the maps, occupation of the Caatinga is the result of a single dispersal event, either jump-dispersal (DIVA-like + J) or anagenetic dispersal to a wide range including the Caatinga immediately followed by vicariance (DIVA-like). Using 100 stochastic maps resolved from the ancestral range estimates on the MCC tree, the most likely candidates for the area that originated this single dispersal are Mesoamerican dry forests, Caribbean dry forests, Argentinean Monte, Atacama desert or Andean dry forests (Table 1). When accounting for topological uncertainty using 100 stochastic maps for each of the 100 posterior trees, the results are similar but there is a substantial increase in the likelihood of Andean dry forests to be the source area, and this area actually becomes the most likely source (Table 1); using only the MCC tree underestimates this possibility. This increase in the likelihood is related to the presence of the rogue taxon Sicarius andinus, which inhabits Andean dry forests and is resolved as the sister taxon to the Caatinga clade in several trees of the posterior distribution.  We inferred the age of the dispersal event to the Caatinga and discovered a similar pattern to what we observed regarding the Galapagos. When using a MCC tree, the estimates are sharp and overconfident (Figure 7a,b), especially in the analysis including founder-event speciation, which shows four peaks closely tied to the ages of the cladogenetic events immediately leading to the Caatinga clade. Using 100 posterior trees to estimate the age of the dispersal event incorporates the uncertainty in the node ages, and reduces the differences between DIVA-like and DIVA-like + J (Figure 7c,d).

Summary of Biogeographic Events in the Face of Uncertainty
We were able to count the most frequently inferred biogeographic events (dispersals, extinctions, allopatry and within-area speciation) by taking into account uncertainty in topology and ancestral range estimates. The results are summarized in Figure 2 and Table 2. A substantial fraction of Sicarius diversity has been generated by within-area speciation, mainly in the Caatinga and the Atacama and, to a lesser extent, in the Monte scrubland, Sechura desert and Andean dry forests. Three of the dispersal events are robust to both types of uncertainty: one from the Sechura to Galapagos, one from the Atacama to Sechura, and one from the Andes to the Chiquitano dry forest. Other dispersal events are sensitive to uncertainty in topology and/or in ancestral range estimates, but generally involve geographically close areas, such as Mesoamerican and Caribbean dry forests, Atacama and Monte, and Atacama and Andes. The Caatinga has exchanged lineages with a single other area whose identity is uncertain, but candidates are Mesoamerican, Caribbean or Andean dry forests and, less likely, the Atacama desert and the Monte scrubland. Figure 7. Histograms with the inferred age of dispersal of Sicarius to the Caatinga. Analyses based on a maximum clade credibility tree (MCC tree) (a,b) produce sharper estimates closely tied to the ages of cladogenetic events in this tree. This is especially notable in the model including founder-event speciation (a), which displays four peaks, each associated with the age of the four successive nodes in the tree that could have originated the founder-event dispersal to the Caatinga. Accounting for uncertainty in topology and age estimates (c,d) reveals that the age of dispersal is much more uncertain.

Inference of Biogeographic Events in the Face of Uncertainty
In recent times, we have learned that tree topology may reach stability with more sequence data, but some clades remain elusive even in massive phylogenomic datasets [14,15]. In addition, estimates of ages of divergence are based on limited data and many assumptions (e.g., the placement of fossil calibrations) and will always be uncertain; as Bromham [18] neatly pointed out, "paleontological evidence and molecular dates paint history with a broad brush, not fine pen work". In this scenario, it seems unwise to put too much faith in a single tree, which represents only one of many similarly probable topologies and age estimates. Several studies have successfully incorporated such uncertainty in biogeographic inference [27,29]. Our results corroborate this: when combining biogeographic stochastic maps with a sample of trees, there are valuable insights to be gained.
We illustrate this concept by studying particular biogeographic events. Biogeographic maps can be used to estimate both transitions among areas and ages of biogeographic events [2]. However, they are tied to the particular tree that is used as a base, which may bias the biogeographic inferences. When we estimated biogeographic transitions between the Caatinga and other dry Neotropical areas, analyses using only the MCC tree underestimated the possibility of dispersal coming from the Andes when compared to the analyses using 100 posterior trees (Table 1). Perhaps more importantly, the inferred ages of biogeographic events are estimated with overconfidence when stochastic maps are run on a single tree (Figures 3a,b, 4a,b and 7a,b). This is because stochastic maps can only infer the ages of events along the branches and nodes of that particular tree. This overconfidence is potentially problematic, because many biogeographic studies are aimed at linking phylogenetic and geological history; the correlation (or lack thereof) of node ages (or biogeographic events) with geological events is often used to reach biological conclusions (e.g., [50]). If such uncertainty is not taken into account, researchers may achieve inaccurate results. We show that running stochastic maps on a single tree does not fully capture all the possible biogeographic histories of a clade, and thus strongly encourage researchers to use this method in combination with a sample of trees to take advantage of its full potential.
Additionally, it seems that using several trees reduces the differences between different biogeographic models (e.g., DIVA-like and its +J variant). We find that lineages through time by area and inferred ages of dispersal are more similar among models when phylogenetic uncertainty is taken into account (Figures 3c,d, 4c,d and 7c,d). Thus, at least in some cases, phylogenetic uncertainty might be more influential than uncertainty in the choice of a particular biogeographic model.
On a lighter note, our results suggest that phylogenetic uncertainty is not crucial for model selection. We show that the outcome of model selection (DEC, DIVA-like, BAYAREAlike, and their +J variants) is the same regardless if the comparison is performed using the MCC tree or 100 posterior trees (Online Supplementary Figure S12). Thus, it seems that model selection can be performed using the MCC tree, and biogeographic stochastic maps can be run on a sample of posterior trees using only the preferred model.

Cladogenetic Events, Founder-Event Speciation, and Uncertainty
We find that overconfidence in age estimates is stronger for cladogenetic events (relative to anagenetic events) when a single tree is used. Models including founder-event speciation disproportionately favor cladogenetic events [5], and thus are more prone to overconfidence, especially when inferring ages of biogeographic events. Such models have been introduced by Matzke [10,11] as a way to model the possibility of dispersal to a new area being followed by speciation over the course of very few generations-almost instantaneously in an evolutionary timescale. In practice, this parameter allows dispersal to happen simultaneously with cladogenetic events, i.e., at tree nodes. This is opposed to anagenetic dispersal, which happens along the branches. Stochastic maps resolve the age of an anagenetic dispersal at any point of such a branch, but founder-event dispersal is always tied to the age of a particular node. Because of this, ages of events inferred from models including founder-event speciation are estimated with more overconfidence and closely tied to the ages of nodes when a single tree is used to run the stochastic maps (Figures 3a, 4a and 7a). This overconfidence is partly smoothed by the variability in ages that can be incorporated by using a sample of posterior trees (Figures 3c,d, 4c,d and 7c,d). Thus, we recommend incorporating phylogenetic uncertainty when running stochastic maps especially when using models including founder event-speciation. It is especially important to have this in mind since model selection often favors such models [5,11,49]. It should also be noted, however, that even models favoring anagenetic dispersal also suffer from overconfidence when using a single tree.

Lineages through Time by Area
We here explore the concept of "lineages through space and time" plots. Former implementations of such plots have been independently conceived by Spriggs et al. [51] and Ceccarelli et al. [52], who were interested in tracking changes in diversity of two different areas through time. These earlier approaches equated the most likely area estimate at nodes as the "true" ancestral range. A later development by Skeels [53] takes advantage of biogeographic stochastic maps to take the uncertainty in ancestral range estimations into account, but this script was designed to work on the results based on a single tree. We here use a similar approach that also takes phylogenetic uncertainty into account. By dividing the tree into user-defined time slices, it is possible to count the number of lineages occupying each area in each time period. Diversity increases through time as a combination of in-situ speciation and new dispersals into the area. Additionally, diversity can also decrease if the model infers extinctions in a particular area (Supplementary Figure S13), or if the taxon sampling includes fossil tips (ILFM, unpublished data).
These plots may be used as rough approximations to detect areas which might be under special processes. For instance, our data on sand spiders indicate that the Caatinga rapidly accumulated diversity after the initial dispersal to this area (Figure 3), which is consistent with observations that this is one of the most diverse dry forests in the Neotropical region [37]. Nevertheless, it should be noted that these plots present limitations. It is known that DEC and its variants vastly and consistently underestimate extinctions [6]; this has been observed in our data, as extinction was estimated to be close to 0 under some models (Online Supplementary Table S9). Known but unsampled species are effectively "extinct" for the purposes of the model and will not be considered, which may affect the results of the plot, particularly if the sampling is uneven across areas. This is clear when we take a look at the diversification of lineages from Africa: only four in-situ speciation events are inferred ( Figure 3, Table 2), even though the genus Hexophthalma currently includes 8 species [54]. It is clear that diversification in this area is underestimated by our sampling, which only includes three species. Thus, we recommend that such plots are only considered when the sampling of the group is complete, or at least when the unsampled species are not biased to a particular area. Second, even if the sampling is complete, these plots will not account for truly extinct species. Researchers interested in comparing diversification rates among areas, while accounting for extinct lineages, should refer to methods specifically designed for modelling this (e.g., BAMM; [55]). On a lighter note, because the approach we take can accommodate phylogenetic uncertainty, one can include taxa with unknown phylogenetic placement by incorporating them to a sample of trees based only on taxonomic information, as undertaken by [29], thus possibly alleviating the effect of missing taxa.

Sand Spiders Dispersed to Galapagos Paleo-Islands
The Galapagos Islands are volcanic in origin, and while the oldest emerged island is 3-4 million years old [34], there is evidence of submerged paleo-islands to the east [35] and north [36] of the current archipelago. The particular biota of the islands has affinities with those of South America, the Greater Antilles and other Pacific islands [56,57]. It is clear that the only sand spider of the islands, S. utriformis, is of South American origins, as its closest relative lives on the coast of Peru (Figure 2). The estimated age of split of these two species (1.2-22 Myr, 95% highest posterior density interval, median 9.7) exceeds the age of the oldest emerged islands. Accordingly, a biogeographic model allowing for dispersal to the islands in the last 15 Myr is strongly favored in a relation to a model allowing dispersal only in the last 3.5 Myr (Figure 6). In most trees the S. utriformis-S. peruensis split is older than 3.5 Myr; thus the second scenario requires that S. utriformis had originated in coastal Peru through in-situ speciation before 3.5 Myr, dispersed to the islands after their formation, and then went extinct in the continent (Online Supplementary Figure S13). In contrast, the first scenario only requires one dispersal event from the ancestor of the pair of species to the islands, followed by a (costless) allopatric event. Thus, our data fit better a model in which Sicarius reached the Galapagos when the archipelago was formed by currently submerged paleo-islands.
Interestingly, the oldest paleo-island reported by Christie et al. [35] is only~530 km distant from the continent, while the current islands lie at~940 km away from the continent. Thus, the paleo-islands were closer to the continent, which could help explaining how a group with poor dispersal capabilities reached a volcanic archipelago. This scenario is similar to the metapopulation vicariance model that has been proposed to explain the presence of ancient, poorly dispersing groups in recent volcanic archipelagos, particularly in the Galapagos [57,58]. It could also explain other instances of clades whose estimated ages are older than the islands they occupy, such as sheet-web weavers in the Juan Fernández Islands [59].
Is it possible that there were even older paleo-islands in the Galapagos? Christie et al. [35] suggested that it is very possible that the history of the archipelago could be as old as that of the volcanic hotspot, spanning 80-90 million years, and thus other uncharted paleo-islands might exist. Further evidence for this hypothesis has been recently reviewed by Heads & Grehan [57]. To account for this possibility, we compared a model where dispersal to the Galapagos was only possible in last 15 Myr to an unconstrained model where dispersal could happen at any moment-thus, considering the possibility of even older paleo-islands. The data never provide definite support to the 15 Myr model, and in those trees where the split between S. utriformis and S. peruensis is older than 15 Myr, it strongly supports the unconstrained model ( Figure 5). In addition, the inferred age of dispersal to Galapagos in the unconstrained model has some probability of being older than 15 Myr (Figure 4). Thus, our data are unable to reject dispersal to the Galapagos before the age of the oldest recorded paleo-islands. This is in line with the opinion that the archipelago must be very old and uncharted paleo-islands exist [35].
Interestingly, we found a correlation between the support of alternative models and the ages of split ( Figure 5). There is a positive correlation between age of split and AIC weight of the unconstrained model, stronger in the models without founder-event speciation (DEC, r = 0.86; DIVA-like, r = 0.83; DIVA-like + J, r = 0.63, BAYAREA-like + J, r = 0.50). This is probably because these models favor anagenetic events, and thus the amount of time passed is correlated with opportunity for dispersal. This means that even if the split is younger than 15 Myr, but very close to the limit between time slices, the window of opportunity for a dispersal is very narrow, and thus even in these cases the unconstrained model is favored.

Ancient Connections among Neotropical Dry Forests
Our data indicate that the Caatinga clade reached this area as a result of a single dispersal. Such dispersal most likely originated from northern dry forests in the Caribbean and Mesoamerica, or from Andean dry forests. Interestingly, this second possibility only appears as a likely candidate when topological uncertainty is taken into account (Table 1). Such connections seem to support the hypothesis by Pennington et al. [38] that some areas of present-day Amazonia have been replaced by xeric vegetation in the past. They compiled detailed evidence that compellingly indicates dry conditions in the region during the Quaternary. Phylogenetic patterns of other organisms, such as birds, support the hypothesis of recent connections among Neotropical dry forests (e.g., [60]). Our data, on the other hand, indicate that dispersal of sand spiders to the Caatinga happened as early as in the Oligocene (Figure 7), with no recent dispersals to other areas. Other groups of organisms display similarly ancient histories in this biome (e.g., geckos; [61]). Thus, it may be plausible that currently disjunct Neotropical dry areas might have gone through many periods of connections over the last 30 million years, and that groups with different dispersing capabilities could have responded idiosyncratically to such connections.

Script Availability
The script for replicating the analyses is available as Supplementary Materials S1, or from GitHub (https://github.com/ivanlfm/BGB_BSM_multiple_trees, accessed on 20 June 2021). The code is thoroughly commented and documented, with explanations for each of the options. It has been successfully tested with two additional datasets, one of them including fossil tips, and thus we expect it to be easily adaptable to other researchers' needs. These datasets varied between~30 and~100 taxa,~10 areas (with maximum range size set to~3) and runs included 2-3 free parameters. In each case, analyses run over 100 posterior trees were usually completed between 8 and 12 h on a standard personal computer (Intel ® i5-5200U 2.20 GHz with 4 GB of RAM). Thus, we expect that the analyses outlined above are not too demanding computationally. As it relies on BioGeoBEARS, the script can be used with dated trees produced by any phylogenetic software.

Concluding Remarks
(1) Biogeographic hypothesis-testing, inferences of transitions among areas and age of biogeographic events using biogeographic stochastic maps benefit from running the analysis over a sample of trees, instead of a single, target tree, since they incorporate uncertainty in topology and age estimates that might be relevant to the questions at hand. We provide a broadly customizable R script to run such analyses.
(2) Including phylogenetic uncertainty is especially important in models including a founder-event speciation parameter; ages of biogeographic events estimated under these models are tightly tied to cladogenetic events, such that using a single tree results in overconfident estimates that disregard the uncertainty in age estimates present in trees from the posterior distribution.
(3) Our data strongly suggest Sicarius most likely dispersed to the Galapagos Islands before the formation of the oldest emerged island. This is congruent with geological findings that indicate that seamounts along the volcanic hotspot are former paleo-islands of the archipelago that are now submerged.
(4) Our data indicate Sicarius dispersed into the Caatinga around 30 Mya, suggesting an ancient colonization of this area. The route of dispersal is unclear due to topological uncertainty, but most likely consisted of a northern route connecting the Caatinga to the Caribbean and Mesoamerican dry forests, or of a southern route connecting the Caatinga to the Andean dry forests.