Stop and Go – Waves of Tarsier Dispersal Mirror the Genesis of Sulawesi Island

The Indonesian island of Sulawesi harbors a highly endemic and diverse fauna sparking fascination since long before Wallace’s contemplation of biogeographical patterns in the region. Allopatric diversification driven by geological or climatic processes has been identified as the main mechanism shaping present faunal distribution on the island. There is both consensus and conflict among range patterns of terrestrial species pointing to the different effects of vicariant events on once co-distributed taxa. Tarsiers, small nocturnal primates with possible evidence of an Eocene fossil record on the Asian mainland, are at present exclusively found in insular Southeast Asia. Sulawesi is hotspot of tarsier diversity, whereby island colonization and subsequent radiation of this old endemic primate lineage remained largely enigmatic. To resolve the phylogeographic history of Sulawesi tarsiers we analyzed an island-wide sample for a set of five approved autosomal phylogenetic markers (ABCA1, ADORA3, AXIN1, RAG1, and TTR) and the paternally inherited SRY gene. We constructed ML and Bayesian phylogenetic trees and estimated divergence times between tarsier populations. We found that their arrival at the Proto-Sulawesi archipelago coincided with initial Miocene tectonic uplift and hypothesize that tarsiers dispersed over the region in distinct waves. Intra-island diversification was spurred by land emergence and a rapid succession of glacial cycles during the Plio-Pleistocene. Some tarsier range boundaries concur with spatial limits in other taxa backing the notion of centers of faunal endemism on Sulawesi. This congruence, however, has partially been superimposed by taxon-specific dispersal patterns.


Introduction
Situated at the triple junction of the Australian, Eurasian, and Pacific plates the Indonesian island of Sulawesi is part of the tectonically most active region in the world. Geologically, Sulawesi has been formed by converging micro-terranes that successively emerged due to subduction zone processes since the Miocene and fused into one another not before the early Pliocene period around 5 MYA [1]. Today, the island covers more than half of the Wallacean terrestrial biosphere (Fig 1A), and therefore it reflects the hotspot of intermingling Asian and Australian biota [2]. Sulawesi hosts a rich and highly endemic fauna with overlapping geographic ranges between disparate taxa [3][4][5][6][7][8]. It is not surprising that this pattern of regional endemism has been associated with the complex geological evolution of Sulawesi, where Neogene and Quaternary tectonic and climatic changes are assumed to be the key factors driving radiations/divergence [4,7]. Apart from the central mountains several fault systems and basins are prominent relief features on Sulawesi [1,2] (Fig 1B). At least some of these physiographic landmarks determine geographic borders of currently known endemism areas [4,5,7,9] (Fig 1C). It is thus conceivable that each region flanked by these suture lines and depressions temporarily formed its own biogeographic entity, likely isolated from each other through water influx during interglacial periods. However, phylogeographic patterns also depend on time and location of colonization and subsequent range expansion [4]. Likewise, lineage-specific habitat preferences, dispersal abilities and sex-biased dispersal behavior can affect population structuring [4,5,[10][11][12][13], thus leading to phylogenetic signals not matching major geological events. Therefore it remains challenging to subdivide Sulawesi into general biotic units, both for conservation purposes [4] and for the generation of testable biogeographic hypotheses [14].
We here elucidate colonization and diversification of a group of widely distributed arboreal island endemics, the Sulawesi tarsier, whose phylogeography is thought to reflect some very general phenomena related to the island´s biogeography. Modern tarsiers represent the oldest lineage of extant haplorhine primates [15]. Their geographical range is restricted to insular Southeast Asia, where they fall into the three spatially and evolutionary distinct clades of Western (Tarsius bancanus), Philippine (Tarsius syrichta), and Sulawesi tarsiers (Fig 1A). The latter group inhabiting Sulawesi and smaller surrounding islands is the most speciose tarsier clade with Tarsius dentatus, T. fuscus, T. lariang, T. pumilus, and T. wallacei being endemic to mainland Sulawesi [9,16]. Further, several other yet unclassified Sulawesi tarsier taxa have been recognized [9]. Those being part of this study are referred to as Tarsius sp. and/or termed according to their respective sample site. Sulawesi tarsiers show a high degree of regional endemism reported on the level of both genetic makeup and vocalization [7,9]. However, neither phenotypic nor molecular data (autosomal STR genotypes, Y chromosomal and mitochondrial haplotypes) sampled so far allowed a thorough and comprehensive inference of a pan-Sulawesi tarsier phylogeography, either because studies were focused on Central-Sulawesian populations [7,17], or because the applied molecular marker was not adequate for resolving phylogenetic relationships among closely related taxa [9]. Therefore, we analyzed the sex-determining region of the Y-chromosome (SRY) and a set of five autosomal and taxonomic widely applicable phylogenetic markers (ABCA1, ADORA3, AXIN1, RAG1, and TTR) [18] in a tarsier sampling with a broad geographic coverage reflecting both geological and vocal features (Fig 1B-1D and Table 1). On this basis we inferred phylogenetic relationships of distinct tarsier populations and, furthermore, tested the plausibility of tectonic and climatic forces driving speciation on Sulawesi. To this end we correlate paleo-environmental and paleo-geographic data with divergence times among Sulawesi tarsiers. In addition we asked whether and to what extent range fragmentation or lineage-specific life history traits shaped current phylogeographic patterns of these group-living nocturnal primates comparing the tarsier´s spatial distribution to that of other island endemics.

Sample collection and preparation
We collected tissue samples of 65 tarsiers at nine sites in seven study areas located throughout Sulawesi between 2009 and 2010 ( Fig 1D, blacks dots). Tarsiers were localized by tracing their  [1]. Arrows point to topographically significant regions. C) Species ranges of distinct lineages. White continuous lines: Macaque and toad hybrid zones, white dotted lines: Toad ranges deviating from the nearest macaque hybrid zone [4]; Continuous black lines: Distribution of tarsier acoustic forms [9,17], thick black lines and circle indicate the discontinuous range of T. wallacei [17], in black dotted areas, tarsier species boundaries are yet to be determined; White dashed line: population differentiation in a Sulawesian bat species, Thoopterus nigrescens [5] To avoid injuring small blood vessels in the translucent pinna of the ear, we checked their course with the help of a flashlight. We then obtained small ear biopsies (triangles with 2 mm side lengths) from the thin tip of the pinna using scissors and tweezers (cleaned with water, ethanol, and wipes) and applied iodine solution to the tarsier ear. Animals did not required sedation prior to the biopsy. All tarsiers were subsequently released alive at their respective capture site. Permits to capture and sample wild tarsiers (no. SK. 198/IV-SET/2009 and SK. 50/IV-SET/2010), and to enter protected areas (no. SI. 86/Set-3/2009, SI 17/Set-3/2010, and SI. 025/BTNBABUL-1/PK/2010) were issued by the Indonesian CITES management authority, the Directorate General of Forest Protection and Nature Conservation (PHKA). We extracted DNA from tissue samples using a DNeasy Blood and Tissue Kit (Qiagen). All specimens were subjected to whole genome amplification (WGA) using the GenomiPhi DNA Amplification Kit (GE Healthcare) to increase the amount of limited genomic material and to meet national and CITES export and import requirements in terms of CITES Appendix-II species. CITES permits covered export (no. 0517/IV/SATS-LN/2010) and import (E-0179/10) of WGA samples from Indonesia to Germany. Additionally, our study comprised WGA-amplified DNA probes of 95 central Sulawesi tarsiers from seven populations (Fig 1D) sampled during previous studies [7,17]. DNA samples of the Philippine and of the Western tarsier were provided by J. Brosius and J. Schmitz (University of Muenster, Germany), and by Y. Rumpler (Les Hôpitaux Universitaires de Strasbourg, France), respectively.

Sequence data collection
In our analyses we included five nuclear autosomal loci (exonic: ADORA3, AXIN1, RAG1; intronic: ABCA1, TTR) of the "Phylogenomic toolkit" [18], in total yielding around 3400 bp of sequence information. Primers and conditions were partly modified (S1 Table) and applied to a pruned sample set comprising 28 specimens from Sulawesi, and one specimen from each the Philippine and the Western tarsier (Table 1, GenBank accession numbers KP642169-KP642408 and KP642434-KP642493) for PCR amplification and sequencing (for more details see S1 Text). Sulawesian populations/study areas were represented by two individuals each constituting a terminal taxon of a preliminary Cytochrome b maximum likelihood tree estimated from sequence data of 151 specimens ( We further PCR amplified and sequenced (for more details see S1 Text) a fragment of the Y-chromosomal SRY gene for 24 males (630 bp, Genbank accession numbers KPKP642409-KP642433, for more details see Table 1 and S1 Text) adopting PCR primers and conditions from Merker et al. [7]. Data sets from former studies and comprising 35 males [7,17] (Table 1, Genbank accession numbers FJ614510-13, FJ614517-21, FJ614523-34, FJ614562-68, HM115985-91) were incorporated for further analyses.

Sequence data analyses
Models of DNA sequence evolution were selected based on Akaike´s information criterion corrected for small sample sizes (AICc) using Treefinder v. March 2011 [24]. We included substitution models for each locus (S2 Table) into subsequent phylogenetic analyses as outlined below.
Multilocus species trees were estimated from combined nuclear markers of the "Phylogenomic Toolkit" [18] employing Ã Beast v. 1.6.2 [25]. We consulted the ucld.stdev and the coefficient of variation parameter of preliminary simulations under an uncorrelated relaxed lognormal clock to reveal deviations from a global clock. Both parameters indicated low rate variation among lineages Thus, we used a strict molecular clock model to analyze multilocus sequence data. We ran ten independent chains of 1x10E7 length sampling every 10,000 generations and leaving the other priors at their default values. For MCMC diagnostics we examined ESS values of the Beast output in Tracer v. 1.5 [26] and Awty [27]. In all analyses, ESS values were generally above 1,000, chains converged and mixed well. Awty´s plots for comparisons of split frequencies across all paired MCMC chains further confirmed convergence. We combined two of the ten analyses to create a final multilocus species tree.
We reconstructed phylogenies for Y-chromosomal data conducting maximum likelihood (ML) and Bayesian analyses. ML trees were estimated using Garli v. 2.0 [28] with two search replicates and stepwise-addition starting trees. Bootstrap support was assessed by 100 replicates of which a majority rule consensus tree was inferred with the program Consense of the Phylip package v. 3.69 [29,30]. Bayesian phylogenies were generated by MrBayes v. 3.2 [31]. We ran four independent analyses of 5x10E6 with four Markov chains, sampling every 1,000 generations. We verified convergence of the four runs by examining the average standard deviation of split frequencies (ASDSF), the stationarity of log likelihood values (LnL) of the cold chain, the chain swap acceptance rate, and the potential scale reduction factor (PSRF). ASDSF decreased to 0.01 within the first 800,000 generations, LnL reached stationarity distribution in the first thousand generations, the swap acceptance rate between chains was between 20-60%, and PSRF values were close or equal to 1. On this basis, we concluded that all four runs converged, discarded the burn-in (25%) and built a majority rule tree from the resulting trees.
Moreover, we reconstructed haplotype networks of each autosomal marker and of the SRY gene based on non-identical sequences and using the TCS method [32] as implemented in the PopART package v. 1.7 [33].

Estimating divergence times
We also applied Ã Beast v. 1.6.2 on the five combined nuclear gene markers for dating divergence within Tarsius. Divergence time estimates recently published by Perelman et al. [34] were used to calibrate seven nodes outside the genus Tarsius. For this we complemented our multiple sequence alignments with anthropoid and strepsirhine species (sequence resources are listed in S3 Table) and created eight taxon sets: Homininae, Hominidae, Hominoidea, Catarrhini, Anthropoidea, Tarsius, Haplorhini (monophyletic), and Strepsirrhini. Models of nucleotide evolution (S2 Table) and substitution rates were determined as described above. Accordingly, divergence times (Table 2 and S4 Table) were estimated under an uncorrelated lognormal relaxed clock using a gamma distribution as prior on branch-specific substitution rates (shape 0.001, scale 1000). We applied normal distribution priors to the seven calibration points and added information for the mean time to the most recent common ancestor (tmrca) in million years ago (MYA) and a standard deviation, both derived from node ages and 95% credibility intervals estimated by Perelman et al. [34]. For node calibration we set divergence times as follows: 1) Homininae 6.68/0. 64 69. Three independent Beast analyses were run for 5x10E7 generations saving every 10,000th tree. Trace files were analyzed in Tracer v. 1.5 [26] and Awty [27] to confirm convergence of the Markov chain.

Nuclear DNA Sequence Data
Multilocus Bayesian species tree inference revealed two evolutionary lineages (pp = 1.0) on Sulawesi (Fig 2A) hereafter referred to as lineage 1 and lineage 2 (Figs 2A and 1E). The   lariang define species-and lineage-specific clusters pointing to a recent common ancestor. In ADORA3-and TTR-networks (S3 and S6 Figs) species-clusters only fairly correspond to those inferred from the multilocus species tree analysis. The comparatively short lengths of ADOR-A3-sequences could be responsible for the low information value regarding phylogenetic relationships. Masking gap positions (see PopART documentation) led to the collapsing of several alleles into one node and, therefore, may also to the loss of relevant phylogenetic information in the TTR-network. Overall, parsimony networks of autosomal alleles show that, despite patterns of incomplete lineage sorting, speciation processes in Sulawesi tarsiers are well advanced.
The pan-Sulawesi sample comprised ten SRY-haplotypes found in 59 males. Six haplotypes were unique to their sample location (BAN, DUA, KEN, LAB, LUW, OGA), four haplotypes were species-specific and/or shared between adjacent populations (Figs 1D and 2B). We reconstructed ML and Bayesian gene trees based on one Philippine and ten Sulawesi SRY haplotypes (Fig 2B). Phylogenetic analyses revealed that Sulawesi tarsiers split into two major paternal lineages (ML bootstrap value = 100, pp = 1.0). As inferred by Bayesian species tree inference KEN was grouped together with T. wallacei, DUA, LAB, and OGA. T. fuscus, T. lariang, T. dentatus, KOR, and LUW constituted the second lineage. Within lineage 2 T. dentatus clustered with KOR and LUW, with one haplotype being common to all six study sites and another haplotype being unique to LUW (Fig 2B). The haplotype network (Fig 3) broadly confirms these results by separating northern populations of lineage 1 from lineage 2-species and by the slightly lower sequence divergence between the KEN-haplotype and haplotype-clusters of northern populations. Within lineage 2 all three species (T. dentatus, T. fuscus, and T. lariang) form contiguous clusters.

Divergence Times
Tarsiers have a scarce fossil record [35,36] notoriously impeding the choice of calibration points for a robust estimation of divergence times on the tarsier lineage. We therefore used mean node age estimates of a recently published primate phylogeny based on comprehensive genomic and fossil data [34]. Furthermore, we applied a relaxed clock method to increase the accuracy of divergence time estimation [37].

Tarsier Phylogeography and Island Formation
Over the past two decades paleo-geographic maps of the Wallacea region have been continually revised and adapted to the current state of knowledge [1,[38][39][40][41][42]. Reconstructing areas of subaerial land from the absence of marine deposits is still a geologists' method of choice. The call for cross-fertilization of geological and biological inferences, however, became louder over the last years [40]. Thus, tracing the dispersal history of Sulawesi's old endemics might add new tesserae to reconstruct the geological and environmental setting of past epochs. Arboreal tarsiers and their extinct relatives occupied tropical rainforests in Southeast Asia since the Eocene epoch [35]. They are strong habitat specialists, for which reason the deduction of their dispersal routes may be informative about forested areas at times of dispersal. In order to identify drivers of faunal diversification, we interpret estimated divergence times among extant tarsiers in relation to past geological and paleoclimatic events. We further linked our divergence time estimates to proposed land and sea distributions at the paleo-Sulawesi archipelago [40][41][42][43] and to present-day topographical maps indicating interglacial shallow seas and presumably flooded low lying land areas at times of historical sea level increase. We used this information to locate potential land bridges and contemporary reproductive barriers in the geological past and thus inferred the most-parsimonious albeit hypothetical dispersal pattern resulting in modern species distribution (Fig 4).
Our genetic analyses do not contradict a monophyletic origin of Sulawesi tarsiers and point to their separation from other crown tarsiers (also including the last common ancestor of Philippine and Western tarsiers) between late Oligocene and early Miocene (Fig 1E). This period was marked by a notable drop in sea level [19] and might therefore have facilitated the colonization of Wallacea via dispersal [43]. Moreover, the estimated median node age at 22.3 MYA relates well to the mid-Miocene collision of the Sula-Spur with the Sundaland margin [1,40,43] (Fig 1E) and the resulting emergence of land in the paleo-Sulawesi archipelago. Ophiolite emplacements are indicative for areas of emergent land forming southeast Sulawesi from early Miocene onwards [40,43], thus offering a conceivable point of arrival for the progenitor of Sulawesi tarsiers. However, paleo-geographic maps also indicate the existence of a small area of subaerial land in the region of today´s southwest peninsula of Sulawesi since the early Miocene [1,38,40] allowing for colonizing the island from this direction.
Our broad (although not all-encompassing) sampling of extant taxa suggests that, assuming no other yet undetected lineage split off before (see next section), diversification in crown Sulawesi tarsiers had begun at earliest 13 Ma after initial colonization of Sulawesi, most likely reflecting the emergence of scattered land masses in the region. At the Plio-Pleistocene border crown Sulawesi tarsiers split up into two lineages (Fig 1E). It is remarkable that this divergence event roughly corresponds to a glacial maximum at 2.5 MYA [19]. Assuming that tarsiers first populated Sulawesi´s south (either south-east or south-west), episodic exposure of shelves during the ice-age at 2.5 MYA may have enabled tarsiers to cross the Gulf of Bone [1,13,38,44]. This in turn could have induced the formation of the two lineages, one expanding its range from southeast to north Sulawesi (lineage 1), and the other dispersing in a southwest-northeast direction (lineage 2, Fig 4).
The relatively deep divergence between southeastern and northern populations of lineage 1 might be the result of strike-slip faulting and the related formation of the Malili-Lake-System in southeast Sulawesi beginning in the Pliocene [45]. Assessing the phylogeographic structure in northern Sulawesi, our data did not clearly resolve evolutionary relationships between populations ranging from Ogatemuku (OGA) to the most northeastern tip of the island. Interglacial oceanic inundations, especially during the Pleistocene, are thought to have provoked allopatric speciation in several taxa inhabiting the northern peninsula [2][3][4]. However, sea-level induced range fragmentations may have been insufficently long-lasting for speciation to be completed in northern tarsier populations. Alternatively, reduced distribution ranges, either naturally or anthropogenically effected, or decreased population densities can promote interspecific hybridization [46,47]. Relaxed assortative mating in order to avoid inbreeding among recently evolved tarsier populations inhabiting the narrow northern peninsula could therefore have led to the collapse of newly acquired reproductive barriers. A third option would be that extant populations of LAB and DUA are descendants from OGA and only recently colonized the central and most northeastern part of the northern peninsula.
Based on nuclear sequence data lineage 2-populations fall into two distinct clades, with Tarsius lariang being sister to T. fuscus, T. dentatus and the most eastern populations of KOR and LUW (Fig 1E). Based on our divergence time estimates we hypothesize two dispersal waves out of southwest Sulawesi, which was likely a separate island until mid-Pleistocene [48]. In this period (Proto-) Tarsius lariang could have reached west Sulawesi in a first wave (Fig 4). In a second wave, individuals may have colonized central parts of the island following an ecological gradient from dry to moist climatic conditions [49], which might have facilitated differentiation between T. fuscus and T. dentatus. Considering present cross-taxon congruence in endemism around Lake Tempe ( [4], this study), recurrent flooding of the Tempe Depression during the Pleistocene [2] likely played a central role in driving and maintaining species divergence. In the course of fusing land masses, T. dentatus crossed the Palu-Koro fault [7] and superseded lineage 1-species from lowland rainforests in central and east Sulawesi (Figs 1C-1D and 4). Molecular evidence presented here confirms previous acoustic studies assigning tarsiers of the eastern peninsula to T. dentatus [9,17]. In case of the still hybridizing tarsier species in central Sulawesi, the paleo-environmental setting around the Palu Valley seems to have maintained rather than initiated reproductive isolation [7] between the two recently diverged species T. dentatus and T. lariang.
Our reconstruction of a crown tarsiid phylogeny again strengthens a sister group relationship between Western and Philippine tarsiers suggesting their divergence in late Miocene [9,50]. The median node age at 9.8 MYA perfectly matches the period of lowest tertiary sea level [1,19,51] possibly facilitating tarsier dispersal from Sundaland to Mindanao (Philippines) [52].
Aware that the use of secondary calibration points bears the risk of erroneous divergence time estimates [53], we are confident that relatively narrow confidence intervals, the obvious correspondence between tarsier divergence times and major geological events and climate fluctuations, as well as good node support for several sub-clades build a reliable framework for our inferences about tarsier diversification. Moreover, median node ages of crown tarsiers, the Western-/Philippine tarsier split and the initial speciation event in Sulawesi tarsiers are consistent with estimates by [9,50] further strengthening the reliability of internal split times obtained in this study.

Pathways to Sulawesi and Radiation to Mountain and Offshore Regions
Tarsiers likely crossed Wallace´s Line (see Fig 1A) taking advantage of episodically exposed dispersal corridors connecting the Sundaic region and the paleo-Sulawesi archipelago. Passages across the Makassar Strait and along the Sunda Arc cover the most conceivable routes (see Fig  1A). The Strait of Makassar forms a deep sea barrier between Borneo and Sulawesi since the Eocene. Intermittent land connections especially in its central regions are considered unlikely [41]. Further south temporarily emergent volcanic chains and carbonate platforms extending from Java over the East Java Sea to South Sulawesi could have served as stepping stones during Paleo-Neogene periods [1,48,54]. But, even though offering dispersal opportunities over geological times, the tectonic and paleoenvironmental instability of this region seems to have significantly decreased dispersal success of tarsiers. Our data suggest a single entry of a founder population into Sulawesi, likely going along with the separation of crown Sulawesi tarsiers from other crown tarsiers during the Miocene (see previous section). The long time gap between island colonization and diversification of Sulawesi tarsiers, however, allows alternative scenarios. In the first scenario, crown Sulawesi tarsiers made several attempts to colonize Sulawesi from the Proto-Java islands. One single lineage succeeded at the Plio-Pleistocene border and subsequently underwent further diversification while other lineages became extinct. In the second scenario, Tarsiers´arrival on Sulawesi indeed falls into the Miocene epoch, but initial cladogenesis within Sulawesi tarsiers predates the corresponding bifurcation in our inferred species tree as a result of incomplete taxon sampling [55]. Although our sampling comprises a wide geographic coverage of Sulawesi lowland tarsiers, it lacks data on the montane pygmy tarsier [56] and on lowland species inhabiting some offshore islands near Sulawesi [16,[56][57][58][59] (Fig 4). We therefore cannot exclude insufficient taxon coverage as a possible source for the time gap between tarsier arrival and initial (observed) speciation on Sulawesi. However, it remains unclear at this point whether offshore and montane species are phylogenetically embedded within or evolved independently from the two major Sulawesi tarsier lineages described here.

Units, Shifts, and Conservation
The Wallacea region is a hotspot of biodiversity [60]. However, as complex as its geological history is the spatial distribution of species on Sulawesi, the largest Wallacean island. Regional range overlaps of independent evolutionary lineages point to a shared biogeographic history, while at the same time highlighting discrepancies to other local species assemblages [8,12,[61][62][63][64]. Nevertheless, there is a wide consensus for dispersal as the prevailing mode of island colonization-contrary to the hypothesis of ancient vicariance via micro-continental drift-and Plio-/Pleistocene diversification of most of Sulawesi´s extant terrestrial and limnic fauna [4,7,13,14,43,63,65] (see Results of this study). Cross-taxon congruence of distantly related terrestrial species like Sulawesi toads, macaques, and tarsiers ( Fig 1C) [4] is strong evidence that at least parts of the island temporarily constituted isolated biotic entities in the more recent geological past. Current geographical and genetic structuring of Sulawesi tarsiers suggests that post-speciation range shifts may have moved species boundaries thus blurring positional signals of past reproductive isolation. Therefore, ancient (physical) barriers to gene flow do not necessarily qualify for defining general areas of conservation, although they provide useful guidance in identifying regions of genetic endemism [4]. Finally, further phylogeographic studies are needed to fully understand the history and progress of land formation on Sulawesi and its consequences for species diversity. Future work should in particular take advantage of powerful molecular tools rather than analyzing single (usually mitochondrial) genetic loci with higher susceptibility to incomplete lineage sorting, gender bias, and accelerated evolutionary rates that can lead to overestimates of divergence times [66][67][68][69][70][71]     the Primate Research Center at Bogor Agricultural University) for assistance with authorities and administration. Our thanks also go to numerous field assistants, who helped collecting samples.