Pongo tapanuliensis, Alexander Nater & Maja P. Mattle-Greminger & nton Nurcahyo & Matthew G. Nowak & Marc de Manuel & Tariq Desai & Colin Groves & Marc Pybus & Tugce Bilgin Sonay & Christian Roos & Adriano R. Lameira & Serge A. Wich & James Askew & Marina Davila-Ross & Gabriella Fredriksson & Guillem de Valles & Ferran Casals & avier Prado-Martinez & Benoit Goossens & Ernst J. Verschoor & Kristin S. Warren & Ian Singleton, 2017

Alexander Nater, Maja P. Mattle-Greminger, Anton Nurcahyo, Matthew G. Nowak, Marc de Manuel, Tariq Desai, Colin Groves, Marc Pybus, Tugce Bilgin Sonay, Christian Roos, Adriano R. Lameira, Serge A. Wich, James Askew, Marina Davila-Ross, Gabriella Fredriksson, Guillem de Valles, Ferran Casals, avier Prado-Martinez, Benoit Goossens, Ernst J. Verschoor, Kristin S. Warren, Ian Singleton, David A. Marques, Joko Pamungkas, Dyah Perwitasari-Farajallah, Puji Rianti, Augustine Tuuga, Ivo G. Gut, Marta Gut, Pablo Orozco-terWengel, Carel P. van Schaik, Jaume Bertranpetit, Maria Anisimova, Aylwyn Scally, Tomas Marques-Bonet, Erik Meijaard & Michael Krüzen, 2017, Morphometric, Behavioral, and Genomic Evidence for a New Orangutan Species, Current Biology 27, pp. 1-40 : 2-7

Pongo tapanuliensis sp. nov. Nurcahyo, Meijaard, Nowak, Fredriksson & Groves.

Tapanuli Orangutan.

Etymology

The species name refers to three North Sumatran districts (North, Central, and South Tapanuli) to which P. tapanuliensis is endemic.

Holotype

The complete skeleton of an adult male orangutan that died from wounds sustained by local villagers in November 2013 near Sugi Tonga, Marancar, Tapanuli (Batang Toru) Forest Complex (1 35’54.1’’N, 99 16’36.5’’E), South Tapanuli District, North Sumatra, Indonesia. Skull and postcranium are lodged in the Museum Zoologicum Bogoriense, Indonesia, under accession number MZB39182 . High-resolution 3D reconstructions of the skull and mandible are available from MorphoBank, http://morphobank.org/permalink/?P2591.

Paratypes

Adult individuals of P. tapanuliensis (P2591 - M435788 – P2591 - M435790 ) photographed by Tim Laman in the Batang Toru Forest Complex (1 41’9.1’’N, 98 59’38.1’’E), North Tapanuli District, North Sumatra, Indonesia. Paratypes are available from MorphoBank at http://morphobank.org/permalink/?P2591.

Differential Diagnosis

We compared the holotype to a comprehensive comparative dataset of 33 adult male orangutans from ten institutions housing osteological specimens. Summary statistics for all measurements are listed in Tables S1–S 3 . Pongo tapanuliensis differs from all extant orangutans in the breadth of the upper canine (21.5 versus <20.9 mm), the shallow face depth (6.0 versus>8.4 mm), the narrower interpterygoid distance (at posterior end of pterygoids 33.8 versus>43.9 mm; at anterior end of pterygoids, 33.7 versus>43.0 mm), the shorter tympanic tube (23.9 versus>28.4 mm, mostly>30 mm), the shorter temporomandibular joint (22.5 versus>24.7 mm), the narrower maxillary incisor row (28.3 versus>30.1 mm), the narrower distance across the palate at the first molars (62.7 versus>65.7 mm), the shorter horizontal length of the mandibular symphysis (49.3 versus>53.7 mm), the smaller inferior transverse torus (horizontal length from anterior surface of symphysis 31.8 compared to>36.0 mm), and the width of the ascending ramus of the mandible (55.9 versus>56.3 mm).

Pongo tapanuliensis differs specifically from P.abelii by its deep suborbital fossa, triangular pyriform aperture, and angled facial profile; the longer nuchal surface (70.5 versus <64.7 mm); the wider rostrum, posterior to the canines (59.9 versus <59 mm); the narrower orbits (33.8 versus <34.6 mm); the shorter (29.2 versus>30.0 mm) and narrower (23.2 versus>23.3 mm) foramen magnum; the narrower bicondylar breadth (120.0 versus>127.2 mm); the narrower mandibular incisor row (24.4 versus>28.3 mm); and the greater mesio-distal length of the upper canine (19.4 versus <17.6 mm). The male long call has a higher maximum frequency range of the roar pulse type (>800 versus <747 Hz) with a higher ‘‘shape’’ (>952 versus <934 Hz/s).

Pongo tapanuliensis differs from P. pygmaeus by possessing a nearly straight zygomaxillary suture and lower orbit (orbit height 33.4 versus>35.3 mm); the male long call has a longer duration (>111 versus <90 s) with a greater number of pulses (>52 versus <45 pulses), and is delivered at a greater rate (>0.82 versus <0.79 pulses per 20 s).

Pongo tapanuliensis differs specifically from Pongo ‘‘pygmaeus’’ palaeosumatrensis in the smaller size of the first upper molar (mesio-distal length 13.7 versus>14.0 mm, buccolingual breadth 11.4 versus>12.1 mm, crown area 155.2 versus>175.5 mm2; Figure S1 ).

Description

Craniometrically, the type skull of P. tapanuliensis ( Figure 1B ) is significantly smaller than any skull of comparable developmental stage of other orangutans; it falls outside of the interquartile ranges of P. abelii and P. pygmaeus for 24 of 39 cranio-mandibular measurements ( Table S1 ). A principal-component analysis (PCA) of 26 cranio-mandibular measurements commonly used in primate taxonomic classification [ 5, 6] shows consistent differences between P. tapanuliensis and the two currently recognized species ( Figures 1C and S 2 ).

The external morphology of P. tapanuliensis is more similar to that of P. abelii in its linear body build and more cinnamon pelage than that of P. pygmaeus. The hair texture of P. tapanuliensis is frizzier, contrasting in particular with the long, loose body hair of P. abelii. Pongo tapanuliensis has a prominent moustache and flat flanges covered in downy hair in dominant males, whereas flanges of older males resemble more those of Bornean males. Females of P.tapanuliensis have beards,unlike those of P. pygmaeus.

Distribution

Pongo tapanuliensis occurs only in a small number of forest fragments in the districts of Central, North, and South Tapanuli, Indonesia ( Figure 1A ). The total distribution covers approximately 1,000 km2, with an estimated population size of fewer than 800 individuals [ 7]. The current distribution of P. tapanuliensis is almost completely restricted to medium elevation hill and submontane forest (300–1300 m above sea level) [ 7 – 9]. Although densities are highest in primary forest, it does occur at lower densities in mixed agroforest at the edge of primary forest areas [ 10, 11]. Until relatively recently, P. tapanuliensis was more widespread to the south and west of the current distribution, although evidence for this is largely anecdotal [ 12, 13].

DISCUSSION

Other hominoid species and subspecies were previously described using standard univariate and multivariate techniques to quantify morphological character differences. The elevation of bonobos (P. paniscus) from a subspecies to a species dates back to Coolidge [ 14] and was based on summary statistics of primarily morphological data from a single female specimen of P. paniscus, five available P. paniscus skulls, and comparative data of what is now P. troglodytes. Groves and colleagues [ 5] and Shea et al. [ 15] supported Coolidge’s proposal using larger sample sizes and discriminant function analyses. Shea et al. [ 15] remarked that the species designation for P. paniscus, which was largely based on morphological comparisons, was ultimately strengthened by genetic, ecological, and behavioral data, as we attempted here for Pongo tapanuliensis . For the genus Gorilla, Stumpf et al. [ 16] and Groves [ 17] used craniomandibular data from 747 individuals from 19 geographic regions, confirming a classification of the genus into two species (G. gorilla and G. beringei), as proposed earlier by Groves [ 1]. Other recent primate species descriptions primarily relied on an inconsistent mix of data on pelage color, ecology, morphology, and/or vocalizations [ 18 – 23], with only a few also incorporating genetic analyses [ 24, 25].

Here, we used an integrative approach by corroborating the morphological analysis and behavioral and ecological data with whole-genome data of 37 orangutans with known provenance, covering the entire range of extant orangutans including areas never sampled before ( Figure 2A ; Table S4 ). We applied a model-based approach to statistically evaluate competing demographic models, identify independent evolutionary lineages and infer levels of gene flow and the timing of genetic isolation between lineages. This enabled us to directly compare complex and realistic models of speciation. We refrained from directly comparing genetic differentiation among the three species in the genus Pongo with that of other hominoids, as we deem such comparisons problematic in order to evaluate whether P. tapanuliensis constitutes a new species. This is because estimates of genetic differentiation reflect a combination of divergence time, demographic history, and gene flow and are also influenced by the employed genetic marker system [ 26, 27].

A PCA ( Figure 2B ) of genomic diversity highlighted the divergence between individuals from Borneo and Sumatra (PC1) but also separated P. tapanuliensis from P. abelii (PC2). The same clustering pattern was also found in a model-based analysis of population structure ( Figure 2C ) and is consistent with an earlier genetic study analyzing a larger number of non-invasively collected samples using microsatellite markers [ 28]. However, although such clustering approaches are powerful in detecting extant population structure, population history and speciation cannot be inferred, as these methods are not suited to distinguish between old divergences with gene flow and cases of recent divergence with isolation [ 29, 30]. To address this problem and further investigate the timing of population splits and gene flow, we therefore employed different complementary modeling and phylogenetic approaches.

We applied an approximate Bayesian computation (ABC) approach, which allows inference and comparison of arbitrarily complex demographic modes based on the comparison of the observed genomic data to extensive population genetic simulations [ 31]. Our analyses revealed three deep evolutionary lineages in extant orangutans ( Figures 3A and 3B ). Colonization scenarios in which the earliest split within Pongo occurred between the lineages leading to P. abelii and P. tapanuliensis were much better supported than scenarios in which the earliest split was between Bornean and Sumatran species (model 1 versus model 2, combined posterior probability: 99.91%; Figure 3A ). Of the two best scenarios, a model postulating colonization of both northern Sumatra and Borneo from an ancestral population most likely situated south of Lake Toba on Sumatra had the highest support (model 1a versus model 1b, posterior probability: 97.56%; Figure 3A). Our results supported a scenario in which orangutans from mainland Asia first entered Sundaland south of what is now Lake Toba on Sumatra, the most likely entry point based on paleogeographic reconstructions [ 32]. This ancestral population, of which P. tapanuliensis is a direct descendant, then served as a source for the subsequent different colonization events of what is now Borneo, Java, and northern Sumatra.

We estimated the split time between populations north and south of Lake Toba at 3.4 Ma ( Figure 3B ; Table S5 ). Under our best-fitting model, we found evidence for post-split gene flow across Lake Toba (0.3–0.9 migrants per generation; Table S5 ), which is consistent with highly significant signatures of gene flow between P. abelii and P. tapanuliensis using D statistics (CK, BT, WA, Homo sapiens: D = 0.2819, p <0.00001; WK, BT, LK, Homo sapiens: D = 0.2967, p <0.00001). Such gene flow resulted in higher autosomal affinity of P. tapanuliensis to P. abelii compared to P. pygmaeus in the PCA ( Figure 2B ), explaining the smaller amount of variance captured by PC2 (separating P. tapanuliensis from all other populations) compared to PC1 (separating P. pygmaeus from the Sumatran populations). The parameter estimates from a Bayesian full-likelihood analysis implemented in the software G-PhoCS were in good agreement with those obtained by the ABC analysis, although the split time between populations north and south of Lake Toba was more recent (2.27 Ma; 95% highest posterior density [HPD]: 2.21– 2.35; Table S5 ). The G-PhoCS analysis revealed highly asymmetric gene flow between populations north and south of the Toba caldera, with much lower levels of gene flow into the Batang Toru population from the north than vice versa ( Table S5 ).

The existence of two deep evolutionary lineages among extant Sumatran orangutans was corroborated by phylogenetic analyses based on whole mitochondrial genomes ( Figure 4A ), in which the deepest split occurred between populations north of Lake Toba and all other orangutans at 3.97 Ma (95% HPD: 2.35–5.57). Sumatran orangutans formed a paraphyletic group, with P. tapanuliensis being more closely related to the Bornean lineage from which it diverged 2.41 Ma (1.26–3.42 Ma). In contrast, Bornean populations formed a monophyletic group with a very recent mitochondrial coalescence at 160 ka (94–227 ka).

Due to strong female philopatry [ 33], gene flow in orangutans is almost exclusively male mediated [ 34]. Consistent with these pronounced differences in dispersal behavior, phylogenetic analysis of extensive Y chromosome sequencing data revealed a comparatively recent coalescence of Y chromosomes of all extant orangutans 430 kya ( Figure 4B ). The single available Y-haplotype from P. tapanuliensis was nested within the other Sumatran sequences, pointing at the occurrence of male-mediated gene flow across the Toba divide. Thus, in combination with our modeling results, the sex-specific data highlighted the impact of extraordinarily strong male-biased dispersal in the speciation process of orangutans.

Our analyses revealed significant divergence between P. tapanuliensis and P. abelii ( Figures 3B and 4A ) and low levels of male-mediated gene flow ( Figures 3B and 4B ), which, however, completely ceased 10–20 kya ( Figure 3C ). Populations north and south of Lake Toba on Sumatra had been in genetic contact for most of the time since their split, but there was a marked reduction in gene flow after 100 ka ( Figure 3C ), consistent with habitat destruction caused by the Toba supereruption 73 kya [ 35]. However, P. tapanuliensis and P. abelii have been on independent evolutionary trajectories at least since the late Pleistocene/early Holocene, as gene flow between these

populations has ceased completely 10–20 kya ( Figure 3C ) and

is now impossible because of habitat loss in areas between

the species’ ranges [ 7].

Nowadays, most biologists would probably adopt an operational species definition such as ‘‘a species is a population (or group of populations) with fixed heritable differences from other such populations (or groups of populations)’’ [ 36]. With totally allopatric populations, a ‘‘reproductive isolation’’ criterion, such as is still espoused by adherents of the biological species concept, is not possible [ 37, 38]. Notwithstanding a long-running debate about the role of gene flow during speciation and genetic interpretations of the species concept [ 39, 40], genomic studies have found evidence for many instances of recent or ongoing gene flow between taxa that are recognized as distinct and well-established species. This includes examples within each of the other three hominid genera. A recent genomic study using comparable methods to ours revealed extensive gene flow between Gorilla gorilla and G. beringei until 20–30 ka [ 41]. Similar, albeit older and less extensive, admixture occurred between Pan troglodytes and P. paniscus [ 42] and was also reported for Homo sapiens and H. neanderthalensis [ 43]. Pongo tapanuliensis and P. abelii appear to be further examples, showing diagnostic phenotypic and other distinctions that had persisted in the past despite gene flow between them.

Due to the challenges involved in collecting suitable specimens for morphological and genomic analyses from critically endangered great apes, our description of P. tapanuliensis had to rely on a single skeleton and two individual genomes for our main lines of evidence. When further data become available, a more detailed picture of the morphological and genomic diversity within this species and of the differences to other Pongo species might emerge, which may require further taxonomic revision. However, is not uncommon to describe species based on a single specimen (e.g., [ 44 – 46]), and, importantly, there were consistent differences among orangutan populations from multiple independent lines of evidence, warranting the designation of a new species with the limited data at hand.

With a census size of fewer than 800 individuals [ 7], P. tapanuliensis is the least numerous of all great ape species [ 47]. Its range is located around 100 km from the closest population of P. abelii to the north ( Figure 2A ). A combination of small population size and geographic isolation is of particularly high conservation concern, as it may lead to inbreeding depression [ 48] and threaten population persistence [ 49]. Highlighting this, we discovered extensive runs of homozygosity in the genomes of both P. tapanuliensis individuals ( Figure S3 ), pointing at the occurrence of recent inbreeding.

To ensure long-term survival of P. tapanuliensis , conservation measures need to be implemented swiftly. Due to the rugged terrain, external threats have been primarily limited to road construction, illegal clearing of forests, hunting, killings during crop conflict, and trade in orangutans [ 7, 11]. A hydroelectric development has been proposed recently in the area of highest orangutan density, which could impact up to 8% of P. tapanuliensis ’s habitat. This project might lead to further genetic impoverishment and inbreeding, as it would jeopardize chances of maintaining habitat corridors between the western and eastern range ( Figure 1A ), as well as smaller nature reserves, all of which maintain small populations of P. tapanuliensis .

Figures

Figure 1. Morphological Evidence Supporting a New Orangutan Species (A) Current distribution of Pongo tapanuliensis on Sumatra. The holotype locality is marked with a red star. The area shown in the map is indicated in Figure 2 A. (B) Holotype skull and mandible of P. tapanuliensis from a recently deceased individual from Batang Toru. See also Figure S 1 and Tables S 1 and S 2. (C) Violin plots of the first seven principal components of 26 cranio-mandibular morphological variables of eight north Sumatran P. abelii and 19 Bornean P. pygmaeus individuals of similar developmental state as the P. tapanuliensis holotype skull (black horizontal lines). See also Figure S 2.

Figure 2. Distribution, Genomic Diversity, and Population Structure of the Genus Pongo (A) Sampling areas across the current distribution of orangutans. The contour indicates the extent of the exposed Sunda Shelf during the Last Glacial Maximum. The black rectangle delimits the area shown in Figure 1 A. n indicates the number of sequenced individuals. See also Table S 4. (B) PCA of genomic diversity in Pongo. Axis labels show the percentages of the total variance explained by the first two principal components. Colored bars in the insert represent the distribution of nucleotide diversity in genome-wide 1 - Mb windows across sampling areas. (C) Bayesian clustering analysis of population structure using the program ADMIXTURE. Each vertical bar depicts an individual, with colors representing the inferred ancestry proportions with different assumed numbers of genetic clusters (K, horizontal sections).

Figure 3. Demographic History and Gene Flow in Pongo (A) Model selection by ABC of plausible colonization histories of orangutans on Sundaland. The ABC analyses are based on the comparison of 3,000 non-coding 2 - kb loci randomly distributed across the genome with corresponding data simulated under the different demographic models. The numbers in the black boxes indicate the model’s posterior probability. NT, Sumatran populations north of Lake Toba; ST, the Sumatran population of Batang Toru south of Lake Toba; BO, Bornean populations. (B) ABC parameter estimates based on the full demographic model with colonization pattern inferred in (A). Numbers in gray rectangles represent point estimates of effective population size (Ne). Arrows indicate gene flow among populations, and numbers above the arrows represent point estimates of numbers of migrants per generation. See also Table S 5. (C) Relative cross-coalescent rate (RCCR) analysis for between-species pairs of phased high-coverage genomes. A RCCR close to 1 indicates extensive gene flow between species, and a ratio close to 0 indicates genetic isolation between species pairs. The x axis shows time scaled in years, assuming a generation time of 25 years and an autosomal mutation rate of 1.5 3 10 8 per site per generation. See also Figure S

Figure 4. Sex-Specific Evolutionary History of Orangutans Bayesian phylogenetic trees for (A) mitochondrial genomes and (B) Y chromosomes. The mitochondrial tree is rooted with a human and a central chimpanzee sequence and the Y chromosome tree with a human sequence (not shown). ** Posterior probability = 1.00. (C) Genotype-sharing matrix for mitogenomes (above the diagonal) and Y chromosomes (below the diagonal) for all analyzed male orangutans. A value of 1 indicates that two males have identical genotypes at all polymorphic sites; a value of 0 means that they have different genotypes at all variable positions