Conservation Genetics of the Philippine Tarsier: Cryptic Genetic Variation Restructures Conservation Priorities for an Island Archipelago Primate

Establishment of conservation priorities for primates is a particular concern in the island archipelagos of Southeast Asia, where rates of habitat destruction are among the highest in the world. Conservation programs require knowledge of taxonomic diversity to ensure success. The Philippine tarsier is a flagship species that promotes environmental awareness and a thriving ecotourism economy in the Philippines. However, assessment of its conservation status has been impeded by taxonomic uncertainty, a paucity of field studies, and a lack of vouchered specimens and genetic samples available for study in biodiversity repositories. Consequently, conservation priorities are unclear. In this study we use mitochondrial and nuclear DNA to empirically infer geographic partitioning of genetic variation and to identify evolutionarily distinct lineages for conservation action. The distribution of Philippine tarsier genetic diversity is neither congruent with expectations based on biogeographical patterns documented in other Philippine vertebrates, nor does it agree with the most recent Philippine tarsier taxonomic arrangement. We identify three principal evolutionary lineages that do not correspond to the currently recognized subspecies, highlight the discovery of a novel cryptic and range-restricted subcenter of genetic variation in an unanticipated part of the archipelago, and identify additional geographically structured genetic variation that should be the focus of future studies and conservation action. Conservation of this flagship species necessitates establishment of protected areas and targeted conservation programs within the range of each genetically distinct variant of the Philippine tarsier.


Introduction
Biodiversity-rich tropical forests are being degraded worldwide, and the pace of forest destruction is exceptionally rapid in insular Southeast Asia [1].With only 4-8% of its original forest remaining [2], the Philippines has been designated as both a global conservation biodiversity hotspot [3] and a Megadiverse nation [4]-a distinction shared only with Madagascar.Within this archipelago, the Philippine tarsier, a small endemic nocturnal primate, has been enlisted as a flagship species for an emerging societal conservation movement and an expanding ecotourism industry [1,5].
Traditionally, taxonomy and conservation have been inextricably linked and most conservation strategies have targeted formally named taxonomic units: species or subspecies [6].Although most conservation efforts have targeted these taxonomic entities, conserving finer-grained genetic diversity across a species' range [7,8] is essential to preserving metapopulation dynamics, preventing inbreeding depression, avoiding population collapses, and ultimately ensuring against extinction [9][10][11].Unique evolutionary lineages or genetically defined ''Evolutionarily Significant Units'' [12] are appropriate targets of conservation programs aimed at preserving genetic diversity among and within species; targeting empirically defined distinct evolutionary lineages has the added benefit of potentially removing the subjectivity sometimes associated with traditional taxonomy [13].
Despite being the focus of a disproportionate number of intensely focused studies, multiple lines of evidence suggest that numerous primate taxa await discovery and formal taxonomic description.For example, between the 1975 and 1999 nocturnal primate species diversity grew worldwide by 2.85-fold increase [14][15][16][17][18], and has since climbed by an additional 1.69-fold increase [19].Unique among nocturnal primates, tarsiers are found only in insular Southeast Asia (See Appendix S1).Ten species are recognized currently, with several new taxa recently proposed through bioacoustic analysis and molecular data [15].The Philippine tarsier (Fig. 1) has been the focus of recent attempts to understand morphological variation [16], clarify taxonomy [17], and establish the conservation status of populations [18,19]; these efforts have met with limited success and left unanswered the questions of appropriate targets for conservation action.
In the Philippines, many morphologically indistinguishable yet genetically differentiated ''cryptic'' species [20][21][22] have recently been discovered, suggesting that the codistributed Philippine tarsier might also harbor hidden evolutionary lineages (unrecognized species, genetic variants, putative taxa, etc.), which may be unprotected and not incorporated into current conservation planning.
We assessed the conservation genetics of Philippine tarsier to determine (1) how many distinct evolutionary lineages can be identified, (2) whether genetic structure conforms to biogeographical predictions and/or (3) expectations derived from current taxonomy (three subspecies).Due to limited natural history data or consensus from other sources of information (morphology, bioacoustics, ecology), we argue that genetic data should be used to distinguish evolutionarily distinct tarsier population groups as objective, empirically defined conservation priorities.

Mitochondrial data collection
Genetic material from T. syrichta was nondestructively sampled (ear and tail-tip biopsies) from throughout as much of the species' range as feasible (targeting ranges of all described subspecies; Appendix S1) on the large islands of Mindanao, Samar, Leyte, Bohol, and Dinagat, Philippines (Figs. 1, 2; Appendix S2).Deceased animals were salvaged from illegal roadside animal dealers, local Provincial or City Environmental Natural Resources (PENRO, CENRO) officers or university administrators (confiscations of poached animals by students, animal traders, and bush meat hunters), and from indigenous hunters in forested areas.Salvaged animals were prepared as museum specimens and deposited at the Leyte State University (ViSCA) Natural History Museum, the National Museum of the Philippines, or the University of Kansas Biodiversity Institute (Appendix S2).We sampled 77 individuals of T. syrichta from 17 localities for genetic material from most of the large islands within the Mindanao Pleistocene Aggregate Island Complex (PAIC) [2], including representatives of all recognized subspecies.We sequenced the mitochondrial 12S ribosomal RNA (12S), Cytochrome B (CytB), and NADH Dehydrogenase Subunit 2 (ND2) gene fragments and nine nuclear microsatellites.
Sequencing proved highly problematic for our many degraded tissue samples and approximately half our samples would not amplify for any targeted gene regions.Of the half that could be amplified for .2target regions, initial sequence alignments were produced in Muscle [24] and manual adjustments were made in MacClade 4.08 [25].We estimated the phylogeny for each mitochondrial gene fragment independently using likelihood and concatenated mitochondrial sequences following no observation of statistically significant incongruence between datasets.Exploratory analyses of the combined dataset of 36 individuals inferred relationships that did not strongly conflict with topologies inferred from individual gene fragments.

Microsatellite data collection
We sampled a total of 66 individuals from 17 localities of Tarsius syrichta across the southern Philippines.Each locality was represented by 1-27 individuals (mean = 3.88).DNA was extracted using QIAGEN DNeasy extraction kits and we used nine microsatellite makers previously developed for Tarsius syrichta: T5, T6, T22, T34, T35, T43, T69 [26] and Tarsius spp.: T42, T50 [27].PCR reaction volumes and conditions for the PCR amplifications followed the protocol of Schuelke [28].Conditions of the PCR amplification were as follows: 94uC for 30 s, 56uC for 45 s, 72uC for 45 s, followed by 8 cycles at 94uC for 30 s, 53uC for 45 s, 72uC for 45 s, and a final extension at 72uC for 10 min.Primers T42, T43 and T45 were fluorescently labeled with 6-FAM by Integrated DNA Technologies (Coralville, Iowa) for PCR.Primer annealing temperatures for these three loci were as follows: T42 at 53uC [27], T43 at 54uC and T45 at 57uC [26].Subsequently, 2.1 ml of the PCR product was added to 9 ml of formimide and .15ml was electrophoresed with the LIZ-500 size standard and analyzed on a 3730 Genetic Analyzer (Applied Biosystems).Fragments were sized with GeneMapper version 4.1 (Applied Biosystems).We replicated data collection for all samples a minimum of two times per locus in order to monitor for allelic dropout and false alleles.In cases where we detected allelic dropout or inconsistent genotypes (n = 2), the sample was replicated once more in order to minimize genotyping errors.In a few instances of ambiguous genotypes, despite these efforts, samples were simply excluded from subsequent analyses (n = 8).Microsatellite data were archived at Dryad (doi:10.5061/dryad.r7468).

Phylogenetic analyses and lineage delimitation: mitochondrial data
Partitioned Bayesian analyses of the combined dataset were conducted in MrBayes v3.1.2[29].As the mtDNA genes sampled in this study included incomplete sections for some targeted genes, datasets were partitioned by gene only.The Akaike Information Criterion (AIC) implemented in jModeltest [30] was used to select the model of nucleotide substitution for each subset (Table 1).We ran four independent Metropolis-coupled MCMC analyses, each with four chains and the default heating scheme (temp = 0.2).All analyses were run for 20 million generations, sampling every 1000 generations.To assess stationarity, all sampled log-likelihoods and parameter values from the cold Markov chain were plotted against generation time and compared among independent runs using Tracer v1.4 [31].These runs demonstrated patterns consistent with stationarity after 4 million generations, hence the first 20% of samples were discarded as burn-in.
Partitioned maximum likelihood (ML) analyses were conducted in RAxMLHPC v7.0 [32] for the combined dataset, using the same partitioning strategy as implemented in the Bayesian analyses.The more complex model (GTR + I + C) was used for all subsets, and 100 replicate ML inferences were performed for each run.Each inference was initiated with a random starting tree, employed the rapid hill-climbing algorithm and support was assessed with 1000 bootstraps [33].
To generate an ultrametric phylogeny for our Generalized Mixed Yule-Coalescent (GMYC) analysis, we ran a strict clock divergence time analysis of the mtDNA data, partitioning the data by gene, and fixing the mean substitution rate to 1.0 in BEAST v1.7.5 [34].We first removed all individuals with identical sequences, however, since zero-length branches can cause the GMYC analysis to over-partition the dataset.BEAST analyses were run for 100 million generations, sampling every 10,000.Convergence was assessed with Tracer as described above.We then used the ultrametric phylogeny with the ''single threshold'' GMYC model [35] to determine the number of evolutionary lineages in the mtDNA dataset.

Population structure and gene flow: nuclear data
As currently recognized, Tarsius syrichta spans the Mindanao PAIC in the southern Philippines.As an initial starting point, our a priori hypothesis of population structure, given our understanding that all these islands are part of the same faunal region [2], were that genetic diversity might be partitioned among individual islands.In order to test these expectations, we used the Bayesian clustering method of the program Structure v2.3.3 [27,28,[36][37][38][39] with allelic data from nine nuclear loci, to estimate the number of populations among T. syrichta samples, as well as infer the probabilities of individuals belonging to each of the estimated populations.We considered individuals with q values between 0.10 and 0.90 to be admixed [39][40][41].Due to the absence of a priori knowledge of inter-locus relationships, analyses were run under the correlated allele frequency model and the admixture ancestry model [36,40].All other parameters were left with default settings.Taking a conservative approach to evaluating population boundaries for this study, we considered the number of possible populations, represented by K, to range from one (a single geographically non-structured species) to seven (each island sampled in this study, and allowing for additional population genetic complexity that may be present on the large island of Mindanao [2]).We conducted 10 independent runs in Structure for each of the seven values of K (1-7).We ran each analysis for 10 million iterations, with a burn-in of 500,000.The rate of change of K (DK) was evaluated in Structure Harvester (Web v0.6.93)[42] following the method of Evanno et al. [43] to determine the number of preferred genetic clusters.The results of structure analyses were visualized using the program Distruct v1.1 (Fig. 2) [44].
To assess whether there is evidence of recent migration among populations, we estimated the level of historical gene flow within a Bayesian framework in the program BAYESASS v1.3 [45] for the four evolutionary lineages supported in phylogenetic analyses (Fig. 1).Four independent runs were conducted, each with 20,000,000 iterations, a sampling frequency of 2,000, and a burn-in of 2,000.Trace files were assessed using the program Tracer v1.5 [31] for evidence of convergence prior to result summary.
To visualize population genetic structure, we generated phylogenetic networks for the mitochondrial dataset by employing the Neighbor-Net algorithm [46] in the program SplitsTree v4.10 [47] (Fig. 2).To assess the support for inferred splits in the network, a bootstrap analysis was conducted with 1000 pseudoreplicates.

Genetic variation and gene networks
Of the 77 samples collected, 66 could be characterized for microsatellite variation (Tables 2, 3) but, due to sample degradation, only 36 could be sequenced for mitochondrial DNA.As anticipated, the results of network analyses corroborate the major results observed from phylogenetic analyses (Figs. 1, 2); these same mitochondrial clades were identified as distinct in the GMYC lineage delineation procedure.

Phylogenetic and lineage delimitation analyses
GMYC analysis of the mitochondrial phylogeny identified three putative evolutionary lineages; each of these three lineages was supported in Bayesian and ML phylogeographic analyses and phylogenetic networks of mitochondrial sequence data, and were provisionally adopted as the following distinct evolutionary groups (Fig. 1): (1) Bohol, Samar, and Leyte island, (2) Dinagat and Northeast Mindanao (Caraga Region) island, and (3) Mindanao island, composed of (a) eastern and (b) western subclades; substantial genetic divergence was found among groups (2.1-4.7% in mtDNA; Table 4), with markedly less genetic variation within groups (0.0-1.1%) and minimal gene flow between most groups (Table 5).

Assignment of individuals to population clusters
The ad hoc DK test of the results of Structure analyses of microsatellite loci preferred a K of two (LnP[K] = 22341.58;Dinagat-Caraga vs. all remaining samples), suggesting an initial estimate of only two distinct populations as defined by our nuclear data alone.All individuals were assigned to one of these two genetic clusters with q values above 0.90.With the exception of populations from Dinagat Island and the northeastern peninsula of Mindanao Island (Surigao del Norte Province; Fig. 2), all sampled populations were inferred to be part of one cohesive genetic group.Evaluating a three-population model in structure reveals a similar pattern of unambiguous individual assignments to populations; however, the model of K = 3 which was close to maximally preferred included additional support for the distinctiveness of populations on Bohol (Fig. 2).All Structure analyses provided support for the cryptic Dinagat-Caraga population (Fig. 2).

New Philippine tarsier evolutionary lineages: novel conservation priorities
The unambiguous support across all analyses for a genetically distinct Dinagat-Caraga tarsier lineage, mirrored in all analyses of all loci, is novel, and identifies a range-restricted ''cryptic'' evolutionary entity with unique conservation concerns.Although moderately forested, Dinagat and northeast Mindanao are impoverished economically, lack low-elevation protected areas (Fig. 2), and have become the focus of particularly intensive mining operations-all of which threaten the remaining suitable habitat of this newly documented evolutionary lineage.The Dinagat-Caraga tarsier should therefore be regarded as tantamount to the conservation importance of celebrated Philippine flagship species (e.g., Philippine eagle, tamaraw, golden-spotted monitor lizard, etc.: [2]).Nearby Siargao Island is protected and may harbor the same genetic variant identified here from Dinagat; intensive studies of these populations are urgently needed.
We anticipate future studies will further refine the known distribution of the Dinagat-Caraga tarsier and extend its range to nearby Siargao Island.Additionally it is not clear if the Surigao del Norte (NE Mainland Mindanao, Caraga Region) samples identified here occur naturally on Mindanao or are the result of Dinagat and/or Siargao island animals confiscated from smugglers by regional wildlife enforcement and subsequently released on Mindanao (PSO, personal communication with Surigao del Norte Provincial Environmental Natural Resources Office [PENRO] staff).

Geographical distribution of genetic variation
Our empirical findings are more complex than is reflected in current taxonomy [18,48], stand in contrast to past biogeographic paradigms [2,49], uncover a novel conservation target (Dinagat-Caraga lineage; Figs. 1, 2), and identify a conservation research priority of urgent concern (potential differentiation between western and eastern Mindanao).The moderate level of sequence divergence (2.4%; Table 4) detected between eastern and western Mindanao (both supported as distinct lineages in the GMYC analysis, but with marginally significant migration rates detected in analysis of microsatellite data; Table 4) may indicate that separate tarsier conservation programs for these populations are warranted if future studies confirm their distinctiveness.An east-west sampling transect across Mindanao will be necessary to investigate the genetic relationships of these populations which, given the gap in our sampling, may simply be the extremes of natural geographically based genetic structure.
Our results from analyses of mitochondrial data do not differentiate the populations on Bohol from those on Samar-Leyte (currently recognized as separate subspecies: [48]).Considering the suboptimal K = 3 results, Structure analyses of microsatellites identified Bohol potentially as distinct.However, given that these analyses potentially are sensitive to unequal sampling, and that we possessed many more samples from Bohol than from Samar-Leyte, we suspect this result may be an artifact; future studies with better sampling from Samar and Leyte islands, and additional microsatellite loci will be necessary to investigate further genetic relations within the Bohol-Samar-Leyte clade.At present, in our efforts to adhere to objective, empirically defined evolutionary lineages, we stop short of defining Bohol as a distinct conservation target, and

Combining mitochondrial and nuclear gene signals to establish targets for conservation action
Under current implementation of international conservation status assessment (IUCN, 2014), prioritization of populations for conservation action follows the recognition of named taxonomic units, i.e., species and subspecies.In the face of differing taxonomic arrangements, biogeographic expectations, and differing geographic patterns of nuclear versus mitochondrial genetic variation, how should the Philippine tarsier be prioritized for implementation of applied conservation measures?Currently, one tarsier sanctuary has been established on Bohol Island and another is under consideration for construction on Leyte.Would these two, and only these two efforts adequately conserve genetic components of Philippine tarsier diversity?We argue that they would not.Would establishment of conservation programs based on current tarsier taxonomy (one effort on Bohol, another on Samar-Leyte, and a third on Mindanao) properly conserve the genetic variation we have documented?Similarly, we argue that such an approach would fail to conserve the genetic variation elucidated here.
Brandon-Jones et al. [17] characterized the three subspecies (one from Bohol, another from Leyte-Samar, and a third from Mindanao) as ''dubious'' taxa.They wrote: '' Groves (2001) recognized no subspecies for Tarsius syrichta.Hill (1955) recognized T. s. syrichta/carbonarius/fraterculus as a poorly defined subspecies, perhaps synonymous, with T. s. syrichta.Museum specimen variation seemed insignificant to Niemitz (1984), but an inadequate basis for judgment, according to Musser, and Dagosto (1987).''Our study sheds some light on the uncertainty of Brandon-Jones et al. [17] and suggests that a detailed study employing other lines of evidence is still needed to revise Philippine tarsier taxonomy.One possible outcome of further study could be that our evolutionary groups are not sufficiently distinct based on new lines of evidence (morphology, bioacoustics, whole genome data) to warrant taxonomic separation, and that Philippine tarsiers should therefore be classified as a single taxon, T. syrichta [16,50,51].Alternatively, further investigation might find two or more of our conservation priority groups warrant taxonomic separation, in which case, the epithets applied by Hill [52] will be available (as species or subspecies).Three names are available for a subset of the evolutionary lineages we identify: (1) Samar-Leyte, (2) Bohol, and (3) eastern Mindanao (but not including western Mindanao [Zamboanga], or the novel Dinagat-Caraga lineage).For the immediate future, we argue that applied conservation efforts based upon the combined intersection of our variable genetic results are superior to those based upon existing taxonomy.
A number of genetic and social system phenomena could conceivably account for the differences in lineage-specific support we have inferred between mitochondrial and nuclear loci.Possibilities include persistence of genetic polymorphism, lineage sorting, nuclear gene flow between distinct mitochondrial lineages, sex-biased dispersal, mitochondrial gene sweeps, or an undocumented Philippine tarsier mating system.Given the present sampling, we are unable to distinguish between these possibilities, which may provide intriguing questions for future research.However, from a practical and applied conservation perspective, we argue that a single, clearly optimal, and objective solution results from the combination of patterns observed here in mitochondrial and nuclear gene loci.
Given our variable inferences, we advocate a combined, cumulative approach towards identification and recognition of tarsier evolutionary lineages for conservation planning.Thus, we argue that in order to maximally preserve genetic variation, it is the combined results of our variable sources of information that likely will have the best chance of maximally preserving genetic variation in the Philippine tarsier.As such, we advocate a multitiered conservation program involving conservation programs and protected areas in the ranges of each distinctive population lineages identified here, minimally including (1) Bohol-Samar-Leyte, (2) Dinagat-Caraga, and (3) Mindanao.Until a comprehensive taxonomic study using multiple lines of evidence (molecular, morphological and/or bioacoustic data) can be undertaken, we find that the existing evidence provides an inadequate basis for distinguishing between taxonomic alternatives.Instead, we emphasize that the archipelago's tarsier populations are partitioned into at least the three genetic variants, which we have empirically defined here for practical conservation purposes.We caution primatologists from taking taxonomic action until multiple lines of evidence converge on a meaningful solution [20,22], ideally with information from throughout the genome [53] and more robust sampling along multiple transects from throughout the range of this species (i.e., eastern versus western Mindanao).Nevertheless, our results necessitate a refined conservation strategy in order to effectively preserve the geographical distribution of genetic variation in this flagship species.Such an approach will greatly enhance the prospects for continued survival of this endemic primate and, combined with many other recent discoveries in the country, will contribute to the recognition of the Table 5. Identities by geographical region and numbers (in parentheses) of Philippine tarsier (Tarsius syrichta) evolutionary lineages, putative taxa, or conservation targets inferred (new data and analyses) or predicted (biogeography, taxonomy) from each of the five available sources of information.archipelago as a globally significant biodiversity conservation priority [2,5,49].

Table 1 .
Models of evolution selected by AIC and applied for partitioned, phylogeographic analyses 1 .