Phylogenetic evidence that both ancient vicariance and dispersal have contributed to the biogeographic patterns of anchialine cave shrimps

Cave shrimps from the genera Typhlatya, Stygiocaris and Typhlopatsa (Atyidae) are restricted to specialised coastal subterranean habitats or nearby freshwaters and have a highly disconnected distribution (Eastern Pacific, Caribbean, Atlantic, Mediterranean, Madagascar, Australia). The combination of a wide distribution and a limited dispersal potential suggests a large-scale process has generated this geographic pattern. Tectonic plates that fragment ancestral ranges (vicariance) has often been assumed to cause this process, with the biota as passive passengers on continental blocks. The ancestors of these cave shrimps are believed to have inhabited the ancient Tethys Sea, with three particular geological events hypothesised to have led to their isolation and divergence; (1) the opening of the Atlantic Ocean, (2) the breakup of Gondwana, and (3) the closure of the Tethys Seaway. We test the relative contribution of vicariance and dispersal in the evolutionary history of this group using mitochondrial genomes to reconstruct phylogenetic and biogeographic scenarios with fossil-based calibrations. Given that the Australia/Madagascar shrimp divergence postdates the Gondwanan breakup, our results suggest both vicariance (the Atlantic opening) and dispersal. The Tethys closure appears not to have been influential, however we hypothesise that changing marine currents had an important early influence on their biogeography.

problematic, as evidenced by the aforementioned studies 21,22 , which highlight some methodological issues, in particular alternative fossil calibrations, which can yield younger estimates and thus imply dispersal.
Although adult atyids have never been reported from open marine habitats, it is well established that some members of the family require saline water for larval development 13,23 , which suggests a way in which seemingly unlikely oceanic dispersal could be possible. Anchialine cave water provides a salinity gradient interphase between fresh ground water and marine water and may offer suitable, vertically segregated, environmental conditions for both adults and juveniles 24 . Some anchialine species, including other atyids, occur on oceanic islands thousands of kilometers apart, sometimes with the same mitochondrial haplotypes 15,[25][26][27] . The divergent phylogeographic patterns exhibited among anchialine shrimp have been linked to differences in the duration of their respective planktonic larval (dispersive) phases 26,28 . However, it is the planktonic characteristics of ancestral . Nodes with maximum nodal support (i.e. 1.00/100) are marked with asterisks. Nodes defining the TST phylogeny are numbered for discussion in the main text (nodes 41 to 53). Fossil calibration points are indicated with codes F1 to F4. Inset multi-plot charts show the age density distributions for two relevant nodes (the root of the tree and the origin of the TST clade) resulting from analysing the data set using each single-fossil separately or all-fossils-at-a-time as calibration points. Coloured circles on the left of the TST complex species names refer to their respective geographical distributions according to the map on Fig. 1. populations from millions of years ago, not contemporary ones, that may have generated some of these patterns, and so all we can do is infer and hypothesise. Even so, the fact of a wide current distribution suggests that trans-oceanic dispersal could also have played an important role in the generation of the extreme disjunct pattern exhibited by the TST shrimp-complex.
Previous phylogenetic studies on the TST shrimp-complex have relied on a limited number of species and on relatively short DNA sequences 15,[29][30][31][32] . These studies were inconclusive with regard to the causes and time-frame of the extreme disjunct distribution of the complex. The most comprehensive of these studies 15 was based on mitochondrial and nuclear DNA sequences from representatives of distinct areas, but did not include the Malagasy Typhlopatsa 33 or the newly discovered Zanzibari Typhlatya, and thus could not consider Gondwanan vicariance. The inferred molecular phylogeny was calibrated by constraining three nodes based on palaeogeographic considerations, which provided node age estimates that were inconsistent with a simple vicariant scenario by plate tectonics 15 .
Here we have taken a whole mitochondrial genome approach to resolve phylogenetic and biogeographic relationships by sequencing fourteen species from three genera from the TST complex that represent all of the speciation rates (cool colours = slow, warm = fast) along each branch of the TST complex. Each unique colour section of a branch represents the mean of the marginal posterior density of speciation rates on a localised segment of the tree. Note the temporal deceleration in speciation rates toward tips (warm colours at the root, cool colours at the tips). (b) Evolutionary net diversification rates through time (i.e., speciation minus extinction rates). Colour density shading denotes the 10% through 90% Bayesian credible regions on the distribution of rates at any time.
Scientific RepoRts | 7: 2852 | DOI:10.1038/s41598-017-03107-y geographic regions where it is known to occur. In particular, we test the ancient vicariant origin hypothesis, which if true would be supported by the following predictions: (1) divergence times of sister-species on the two shores of the Atlantic Ocean should be compatible with the opening of this oceanic basin (ca. 110 Ma 34,35 ); (2) divergence times between sister lineages on opposite sides of the Indian Ocean should be consistent with the establishment of deep-water conditions associated with the separation of the Australian and the Indo-Madagascar landmasses of Gondwana (that started at 132 Ma 36 ); and (3) the divergence of Indian Ocean and Mediterranean taxa should  trace back to the occlusion of the Tethys sea (20-16 Ma 37 ). We assess these predictions using fossil-calibrated phylogenetic trees, together with lineage diversification analyses, to appraise the potential roles played by both vicariance and long distance dispersal in the biogeographic history of this group.

Results
Molecular phylogenetics. Sixteen new mitogenomes were produced (Table 1). All sequenced mitogenomes matched the ancestral putative pancrustacean gene order and included 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs) and two ribosomal genes. The final concatenated alignment of the 13 mitochondrial PCGs for the 30 studied taxa (outgroups included) comprised 10,953 bp (Supplementary Table 1). Xia's 38 test showed low levels of substitution saturation even under the assumption of an extremely asymmetric tree (ISS = 0.405, ISSc = 0.818, P two-tailed <0.00001). Similar results were found when testing for saturation at each codon position separately, except for third coding positions where the occurrence of significant levels of saturation was detected under the hypothesis of an extremely asymmetric tree (ISS = 0.741, ISSc = 0.554, P twotailed <0.00001). However, the balance of the phylogenetic tree obtained under the best partitioning scheme (see below) was TCI = 1414 from a set of possible values ranging from 338 to 3654, thus suggesting that our phylogeny is far from being asymmetrical. The best partitioning scheme for the nucleotide-based analyses consisted of subdividing the dataset by codon position and coding DNA strand (Supplementary Tables 2 and 3). Other partitioning strategies analysed resulted in an identical topology (Supplementary Fig. 1a-p) except for the single partition scheme (Supplementary Table 2), which differed in the phylogenetic placement of Typhlatya garciai and T. pearsei, although this partition scheme was the worst ranked based on the Bayesian Information Criterion (BIC). Codon-based substitution models, exclusion of the third coding positions or protein data resulted in the same topology or topologies showing minor differences ( Supplementary Fig. 1h-m and 1p). Phylogenies based on single genes were both less resolved and supported than the tree based on the 13 PCGs (see Supplementary  Fig. 1q for the tree obtained with the widely used cytochrome oxidase subunit 1 marker as an example).
The inferred tree was fully congruent with most recent phylo-mitogenomic proposals on decapod systematics 39 and relationships amongst atyid taxa 32 (Fig. 2). The TST complex was consistently recovered as a monophyletic lineage with maximum nodal support ( Fig. 2 -node 41), rendering the genus Typhlatya paraphyletic as it currently stands since both Stygiocaris and Typhlopatsa nest within it. Within the TST clade, two well-supported monophyletic sub-lineages occur (Fig. 2). One (node 42) includes the Australian genus Stygiocaris and the Malagasy Typhlopatsa, the Mediterranean species T. arfeae and T. miravetensis, and the Caribbean species T. monae. The analysis recovers Stygiocaris and Typhlopatsa as sister-groups with maximum support, and the same holds for the two Mediterranean Typhlatya species. The phylogenetic relationships within the other sub-lineage, which includes Typhlatya species from the western Atlantic (Mexico, Cuba, Bermuda) and eastern Africa (Zanzibar) (node 49), were also identical in all the analyses performed except for the unsupported relative position of the Cuban T. garciai in the protein analysis based on the CAT profile mixture model 40 . The systematic position of T. iliffei (Bermuda) and Typhlatya sp. (Zanzibar) could not be established with confidence within this sub-lineage.
Divergence time and diversification. The best fitting clock model resulted from setting an independent relaxed log-normal clock model for each partition plus a Yule diversification pattern (Supplementary Table 4). Node ages and their respective 95% highest posterior densities (HPD95%) are shown in Supplementary  Table 5 for the analyses with the best Bayes factors using fossil calibrations. The four fossil calibration points (Supplementary Text 1) showed age congruence with each other along the phylogeny as deduced from the cross validation test ( Supplementary Fig. 2, and Fig. 2 for the age density distributions in the root of the tree and for the origin of the TST clade), but were largely incompatible with biogeographic calibrations 15 . The latter rendered much younger node ages and unexpectedly high mitochondrial substitution rates, and their combination with fossil data in the same analysis resulted in a high variation in the substitution rates across branches (i.e. high ucld. stdev values) (Supplementary . The selected chronogram that used only fossil constraints placed the origin of the TST complex in the early Cretaceous (133.6 Ma, HPD95% = 148.9-118.8; Fig. 2 -node 41). Within the TST complex, two main ancestral lineages were both inferred to date to relatively soon afterwards, with the western Atlantic/east African lineage (node 47) first diverging 119.  Table 5).
Diversification took place early in the TST complex as inferred from the significantly negative gamma value, thus rejecting the null hypothesis of a constant diversification rate (γ = -2.71, P one-tailed = 0.003, P two-tailed = 0.007). The same result was obtained when correcting for incomplete taxon sampling (critical γ value = -1.9; P = 0.008) and using the lineages-through-time (LTT) plot ( Supplementary Fig. 6). Accordingly, both LASER 41 and DDD 42 analyses independently resulted in the selection of a diversity-dependent diversification process as the best fitting model for our dataset (Supplementary Tables 6 and 7). Moreover, BAMM 43 analyses revealed that the diversification rate has progressively decreased over time without causing any significant shift in the evolutionary history of the TST complex (posterior probability = 0.92; Fig. 3a). These results are consistent with mean speciation and extinction rates defining a low net diversification rate for this lineage of cave shrimp (0.0088 events/Ma/lineage) (Fig. 3b).
Biogeographical history of the TST complex. LAGRANGE analyses based on three matrices of migration probabilities resulted in similar log-likelihood values: -lnL = 19.24, 19.69 and 21.49 for three, two and a single time intervals, respectively. Hence, we selected the three and two time slices models following the standard cut-off value of two log-likelihood units difference 44 . Both models resulted in the same hypothesis regarding historical biogeography (Fig. 4). Thus, the most recent common ancestor (mrca) of the TST complex (node 41) was inferred to live along the western shores of Tethys, in the area corresponding to the Proto-Caribbean and the Proto-Atlantic, and there it diverged into two primary lineages (Supplementary Table 8). One of them (node 47) remained in the same ancestral area except for the lineage leading to the current Zanzibari species (node 48). Some members of the second main lineage (node 42) remained in the Proto-Caribbean (T. monae), but others became established in what is now the W Mediterranean (node 45) and Indian Ocean (node 43). In both cases the very disjunct distributions within the two lineages (nodes 45 and 48) occurred very early in the diversification of the genus (mid-Cretaceous), before the regional groups of species formed. S-DIVA analyses implied the same biogeographical hypotheses as LAGRANGE (Supplementary Fig. 7) and allowed for hypotheses of nodes in the phylogeny where dispersal and/or vicariance events likely took place. Dispersal was inferred to have occurred in the mrca of the TST clade (node 41) and also in the more weakly supported node 47, while vicariance episodes were inferred for nodes 42, 43, 45, and 48.

Discussion
Our mitogenomic analyses show the TST complex is divided into two main ancestral lineages. A strongly supported clade (node 42) comprises the taxa from NW Australia and Madagascar, the W Mediterranean and one of the species from the Caribbean (T. monae). The other lineage (node 47) contains species from the likely ancestral Colours defining the distribution areas refer to the map in Fig. 1. Coloured rectangles on the left of each node show the most probable inherited ranges for each of the two daughter branches (see Supplementary Table 8 for a complete report on the relative probabilities for all the splits). Grey pie charts on the right of each node reflect the relative probability of occurrence of that split. Histogram on the left of the tree depicts the four most probable biogeographic episodes for the origin of the TST clade. Nodes are numbered according to area of the western Atlantic (Yucatan Peninsula and Cuba, plus Bermuda), as well as an undescribed new species from east Africa (Zanzibar).
The first plate tectonic vicariant scenario of particular interest in TST's evolutionary history relates to sister-species on opposite shores of the Atlantic, and whether this aligns with the establishment of deep-water conditions in the opening Atlantic at 110 Ma 34,35 . The two main lineages present in our phylogeny (nodes 42 and 47) were inferred to date back to the early Cretaceous (148.9-118.8 Ma), a period when the Tethyan Seaway was at its largest. Our estimated time frame for the split between the species separated by the Atlantic Ocean in both ancestral major lineages is compatible with the age of the first establishment of deep-water conditions between the east and west shores of this oceanic basin 34,35 . Thereby, the separation of the lineages currently present on Bermuda and Zanzibar (node 48 in Fig. 2) was estimated to occur at 121-91 Ma, and the split between T. monae (Caribbean) and the western Mediterranean species (T. arfeae, T. miravetensis) (node 45 in Fig. 2) also falls within the same age range, suggesting that the two splits may have been triggered by the disconnection between the east and west Atlantic shores, and thus relate to vicariance via tectonic movements.
The second vicariant scenario on which our phylogenetic data may be informative is that associated with the separation of the Australian and the Indo-Madagascar landmasses of Gondwana. Our phylogenetic analysis recovers a strongly supported sister relationship between the NW Australian Stygiocaris and the monotypic, SW Malagasy genus Typhlopatsa. Both taxa are strictly subterranean and are known to occur only in anchialine karstic caves separated by a 7,000 km span of deep-ocean. Interestingly, the two genera were considered to be closely related by their original descriptor 45 . Moreover, in a previous study using nuclear and mitochondrial markers from numerous surface and cave-dwelling atyids from the Indo-Pacific (but not Typhlopatsa), Typhlatya was recovered as sister to Stygiocaris 31 .
The chronogram with the best score suggests that the Australian-Malagasy clade (node 43) is Late Cretaceous in age (85.1 Ma; 97.9-72.7 Ma). At that time, Antarctica/Australia was already well separated from the Indo-Madagascar landmass as part of the sequential break-up of the supercontinent Gondwana. The estimated age for the split between the Australian and the Malagasy shrimp genera would invalidate a vicariant origin by fragmentation of Gondwana as this would postdate the disconnection of the shallow continental shelf waters by 30-35 Ma on average 36,46 , suggesting instead an origin associated with ancient trans-oceanic dispersal. Recently an ancient vicariant origin by plate tectonics was proposed by Chakrabarty et al. 20 to explain the sister relationship found between the obligate cave-dwelling gobiid fishes Milyeringa (NW Australia) and Typhleotris (Madagascar), which mirror exactly the distribution displayed by the Australian-Malagasy shrimp clade. However reanalyses of this relationship suggests that it also postdates any continental connection by tens of millions of years 22 . Interestingly our S-DIVA analysis inferred vicariance for the same node, but this algorithm has been shown to favour vicariance over early dispersal 47 . In any event, a recent long-distance trans-oceanic dispersal is an unlikely explanation for the Australoid-Malagasy shrimp or fish relationship given the deep genetic divergence between the sister lineages.
The third geologically derived vicariant event that we have considered is the final occlusion of Tethys by the collision of the Arabic and the Anatolian Plates at ca. 20-16 Ma 48 . This scenario is different from the previous two, in which the developing barriers which may have sundered populations were represented by expanding deep waters. In this case, the barrier is represented by land as the eastern Mediterranean was blocked from the Arabian Sea. Our age estimates for the divergences between the species from the Indian Ocean (Zanzibar, Madagascar, W Australia) and the rest of taxa (western Atlantic, Mediterranean) exceeds the date established elsewhere for the occlusion of the connection between the Mediterranean and the Indian Ocean (20-16 Ma 37 ) by at least 70 million years (nodes 42,48). This strongly implies these divergences are unrelated to the closure of the Tethyan Seaway.
However, our deep estimates for divergences at either side of the Tethyan Seaway broadly overlap (node 42: Caribbean/Mediterranean vs Madagascar/Australia 134-104 Ma; node 48: Bermuda vs Zanzibar 121-91 Ma). This suggests that these divergences could have been caused by the same historical biogeographic event, but one much more ancient than the closure of the Tethyan Seaway. One possibility is the early Cretaceous disruption of the circum-equatorial west-flowing current system established through the Tethys since its opening in the Late Jurassic times 49 . This disruption was mainly due to the establishment of new connections of the Tethys to other oceanic basins, plus the broadening of those formerly established. We suggest that this might be the basis of the differentiation that took place within the two major lineages of the TST shrimp complex.
The calibration of molecular clocks and their application to infer the biogeographic history of evolutionary lineages are often contentious and complex 22 . Molecular clock calibrations are usually enforced by assigning well-dated and confidently identified fossils to particular nodes of phylogenies 50 . Node age estimates in our TST phylogeny were explored using four decapod fossils as calibration priors and compared with the results obtained using three geological events hypothesised to have had biogeographic influence on contemporaneous TST populations. Fossil ages were largely congruent in cross validation analyses, which has not always been the case for atyid time-trees 32 . However they were largely incompatible with biogeographic-based geological calibrations ( Supplementary Figs 3-5). It has been argued that calibrations solely based on geological events can hinder biogeographic hypotheses as they can lend to a ´circular reasoning by presupposing the very speciation mode they are trying to test´5 1 , while geological information is advised to be assessed independently of a paleontological one 52 . In addition, calibrations near the root of phylogenetic trees are important, although deep calibrations can lead to spurious results when inferring the age of relatively shallow nodes 53 . Therefore, we have favoured the estimates based on fossil calibrations, at least for deeper nodes not close to the tips of the tree.
Our results on species diversification show that the ancestors of the current species in the TST complex diversified early in the phylogeny. This was followed by a gradual slowdown in diversification in a diversity-dependent manner. The latter might be related to a life-style transition from the open sea to the highly specialised anchialine and subterranean fresh water environments. This is a process presumably associated with the progressive reduction in dispersal abilities 15 . Overall, the troglomorphic adaptations to the subterranean habitat may have ultimately driven these species to become virtually landlocked in coastal aquifers, separated by broad oceanic basins 15 .
Scientific RepoRts | 7: 2852 | DOI:10.1038/s41598-017-03107-y Accordingly, for a formal reconstruction of the biogeographic history of the TST complex we made the assumption that the ancestral shrimp lineages populated -either as adults, larvae or both -a shallow tropical sea that progressively broke-up. We conservatively modelled migration probabilities among shallow waters according to historical plate tectonic information 46 . The pattern obtained suggests the existence of an ancestral shrimp lineage distributed along the Proto-Caribbean and Proto-Atlantic (node 41). This ancestor diverged into the two above-mentioned main lineages, with one remaining largely tied to this ancestral area, except for the branch leading to the species from Zanzibar. The other main lineage resulted in a broad distribution that encompassed not only the Proto-Caribbean, but also the current W Mediterranean and Indian Ocean basins, with implications of both dispersal and vicariance.
In conclusion, we examined the influence of vicariance and dispersal on the evolutionary history of a group of anchialine shrimp through phylo-mitogenomics, biogeographic/diversification analyses and calibration of molecular clocks. Our results reflect the complexity of inferring very ancient, widespread events, since they are distant in time, and are unlikely to be the result of a single, simple process. We found evidence for both major processes in the biogeographic patterns evident in this group. Vicariance via tectonic forces is our current best explanation for the divergence of sister-species separated by the Atlantic, whereas ancient long distance dispersal seems more likely to explain sister-species separated by the Indian Ocean. Further, while we did not find evidence of vicariance due the influence of the relatively recent closure of the Tethyan Seaway, we may have uncovered a more ancient vicariant event. This may be the disruption or shift of marine currents in ancient oceans. The combination of these geological and biological processes may explain the fascinatingly complex patterns of isolation on distant regions exhibited by this group of cave shrimp.

Methods
Taxon sampling. We sequenced the mitogenomes of 14 of the 20 species in the TST complex, plus the mitogenomes of the related Halocaridines fowleri and two specimens of T. galapagensis from different islands of the Galapagos. We also included one existing TST mitogenome (T. miravetensis 54 ), and retrieved the mitogenomes of twelve decapod genera and an euphausiid from GenBank to be used as more divergent outgroups ( Table 1). Selection of outgroup taxa was based mainly upon the existence of fossils identified with high confidence and to assure a reasonably broad age coverage.
DNA purification and mitogenome sequencing. Total DNA was purified and mitogenomes were amplified in large amplicons through long-range PCR as described in Pons et al. 55 . Amplicons were sequenced using a Junior NGS 454 platform and assembled into contigs in CodonCode Aligner (CodonCode Corporation, Dedham, MA, USA). Gene annotation was performed with DOGMA 56 and MITOS 57 WebServers.
Phylogenetic analyses. The aligned dataset comprised the nucleotide sequences of the 13 protein-coding mitochondrial genes (PCGs). Nucleotide substitution saturation was assessed using the Xia's test 38 . The balance of the phylogenetic tree obtained under the best partitioning scheme was measured using the total cophenetic index 58 . PartitionFinder 59 was used to explore the partitioning schemes and evolutionary models for the standard DNA-based and the protein-based analyses. Tree topologies were analysed using maximum likelihood (ML) in IQTREE multicore v. 1.3.12 60 specifying 1000 bootstrap replicates, and compared using BIC. We also explored the effect of excluding third coding positions from the aligned DNA matrix and the implementation of codon-based substitution models 61 , and analysing them at the protein level. In addition, the Bayesian CAT profile mixture model 40 was also explored. Phylogenetic trees were rooted with the euphausiid included in the dataset as an outgroup.
Estimation of divergence times. Molecular dating analyses were performed on the nucleotide sequence dataset after removing the divergent euphausiid outgroup. Tree topology, model parameter values and node ages were co-estimated and optimised in BEAST v. 2.3.2 62 . Bayesian analyses were run for 300 million generations. Mean values and confidence intervals of parameters and ages were estimated using TreeAnnotator in the BEAST package. Node ages were estimated under different clock models such as strict clock, uncorrelated log-normal clock and random local clocks (RLC). Analyses implementing RLC were run for 500-1000 million generations with the first 250 million discarded as burnin. For RLC, we implemented narrower priors as specified in Supplementary Table 9. The ages of four fossils were implemented as constraints to calibrate the relaxed molecular clock. Each calibration point was specified as a log-normal distribution where the limits of its 95% confidence interval was defined by the maximum and the minimum fossil age (Supplementary Text 1 for detailed information); namely: the oldest known dendrobranchiate shrimp (the penaeid Aciculopoda mapesi; Late Devonian/Famennian 374.5-359.2 Ma), the oldest known decapod Palaeopalaemon newberryi (Famennian 382.7-358.9 Ma); the oldest representative of the Stenepodidea (Devonostenopus pennsylvaniensis; Late Devonian Carboniferous-Mississippian 358.9-323.2 Ma) and finally the presumed oldest Gebiidea (Upogebia obscura; Early Triassic 252. . In order to test congruence among fossil calibration points, we performed a cross validation analysis, i.e. comparing four additional BEAST runs, each calibrated with a single fossil. Alternatively, we ran the same analyses but, instead of the above four fossil age constraint priors, used three younger biogeographical events assumed to have affected the diversification of particular TST shrimp lineages as described in Botello et al. 15 (Supplementary Text 1).
Different combinations of partition schemes and clock models were compared based on Bayes Factors calculated with marginal likelihood values estimated from path-sampling analyses in BEAST (40 steps with 10 million generations each). The final state files from regular BEAST analyses were used as the starting state of each step to speed up calculations. RLC analyses failed to converge even when performing steps of 50 million generations and stronger priors. Hence, we also compared models and partitioning schemes using Akaike's information criterion (AIC) through MCMC as implemented in Tracer (AICM) 63 .
Scientific RepoRts | 7: 2852 | DOI:10.1038/s41598-017-03107-y Lineage diversification analyses. Analyses were performed on the selected BEAST dated tree. Deviations from a constant rate pure birth process were tested in the TST complex by calculating the gamma statistic 64 using LASER v. 2.3 41 . The MCCR test 64 was implemented to account for incomplete taxon sampling performing 500 simulations. LASER was also used to compare the fit of various models of diversification with constant and variable diversification rates using AIC. A LTT plot was generated using the R package APE 3.0 65 . BAMM v. 2.0 43 and BAMMtools 66 were used to estimate branch-specific rates of speciation and extinction and to determine the occurrence of distinct evolutionary processes in the phylogeny, with a global sampling fraction of 0.74. We also investigated diversification patterns using the likelihood-based method implemented in the R package DDD 42 , with parameters estimated in BAMM.
Analysis of historical biogeography. The biogeographic history of the TST complex was explored under the dispersal-extinction-cladogenesis model (DEC) in LAGRANGE 67 , and the statistical dispersal-vicariance analysis (S-DIVA) 68 implemented in RASP 3.2 47 . Outgroup taxa were excluded in the preferred BEAST chronogram. We constructed a biogeographic model consisting of coding terminals as occurring in any of the following five areas: Australia (A), Caribbean (including Bermuda) (C), Europe (E), Madagascar (M) and Zanzibar (Z). Migration probabilities were defined to model the changes in connectivity among marine shallow waters as described in Leprieur et al. 46 , with no restriction on the number of regions in which a lineage may be present. Probability values ranged from 0.1 for well-isolated areas to 1.0 for contiguous shallow waters based on  Table 10). The LAGRANGE online configuration tool 67 was used to design the migration probability matrices in different time lapses. In order to find the best-fitting time slice model, we calculated the global maximum likelihood in models with three (134-110, 110-20, and 20 Ma to the present), two (134-110 and 110 Ma to the present), and only one time interval (134 Ma to the present). A cut-off value of two log-likelihood units was applied to determine the existence of significant differences between models 44 . The same five areas scheme were used in S-DIVA analyses, accounting for phylogenetic uncertainty by performing the analyses on 1000 random post-burnin BEAST trees.