Faecal samples (Table 1) were kept in absolute ethanol (99%) and stored at −20°C to preserve the DNA in the faces (Aifat et al. 2021). DNA was extracted from 0.5–1.0 g of faeces using the InnuPREP Stool DNA Kit (Analytik Jena, Germany) according to manufacturer’s protocol (Md-Zain et al. 2019). An approximately 600-bp region of the D-loop gene was amplified by PCR using a species-specific primer pair (GIBDLF3 5′-CTTCACCCTCAGCACCCA AAGC-3′ and GIBDLR4 5’-GGGTGATAGGCCTGTGATC-3’), as described by Whittaker et al. (2007). Amplifications were carried out using MyTaq ™ Red Mix from Bioline in a total reaction volume of 25 µl consisting of 12.5 µl of MyTaq™ Red Mix, 3.0 µl of DNA template, 1.0 µl of each 20 µM primer and 7.5 µl ddH₂O. The thermal cycling programme consisted of a predenaturation at 95°C for 3 min, followed by 30 cycles of denaturation at 95°C for 15 s, annealing at 50°C for 30 s and elongation at 72°C for 10 s and final elongation at 72°C for 10 min. The PCR products were electrophoresed using 1.5% agarose gel in 1× TAE buffer.

Amplified DNA was purified using the InnuPREP Double Pure Kit (Analytik Jena). The purified PCR product was sequenced by Apical Scientific Sdn. Bhd., Malaysia and the resulting sequences were edited using BioEdit software v.7.2.6.1 (Alzohairy 2011). The sequence chromatogram was visually proofread,and, then, the validity of the sequence was checked against other GenBank sequences using BLAST to evaluate quality and degree of match of our DNA sequences to the expected species and gene locus.

The phylogenetic trees were constructed using a number of methods for a robust comparison: distance-based method (Neighbour-joining), character-based method (Maximum Parsimony), likelihood (Maximum Likelihood), Bayesian Inference and molecular clock estimation. Approximately 348 bp was used for phylogenetic and haplotype analyses. Of approximately 255 variable sites from sequence analyses, 217 were parsimony informative. The Neighbour-joining (NJ) and Maximum Parsimony (MP) trees were constructed using Mega X software (Kumar et al. 2018) and Maximum Likelihood (ML) was constructed using PAUP 4.0b10 (Kimura 1980). The tree was constructed with 1000 bootstrap replications to obtain the bootstrap confidence level, which estimates the robustness of the tree topologies. The Bayesian tree was constructed using an appropriate substitution model for D-loop sequences, which was selected by hLRT in ModelTest 3.7. The best model was HKY + G, with a gamma shape parameter of 0.7293 and base frequencies of 0.3495 for adenosine, 0.3136 for cytosine, 0.1103 for guanine and 0.2346 for thymine. A metropolis-coupled Markov Chain Monte Carlo model was run with 10 million generations and the tree was sampled every 100 generations.

The molecular clock was constructed using BEAUti 1.8.4 and BEAST 1.8.4 (Drummond et al. 2012). We defined Siamang datasets in BEAST 1.8.4, an outgroup dataset (Papio papio, Macaca fascicularis and Pongo pygmaeus) was used as a calibration point and the other dataset included members of the family Hylobatidae (Nomascus, Hoolock, Symphalangus and Hylobates). The molecular divergence phylogenetic tree was constructed using the uncorrelated log-normal relaxed clock model (Drummond et al. 2006) and a Yule model before estimating the substitution rate for all nodes in the tree with uniform priors on the mean (0, 100) and standard deviation (0, 10). Published dates with a mean of 33.65 Ma and standard deviation of 3.45 Ma ([CI]: 30–37.1 Ma) were utilised for M. fascicularis and P. papio (Israfil et al. 2011). The tree was run for 10 million generations and sampled every 100 generations, with 10% of the initial trees discarded as burn-in. Tracer version 1.5 was used to assess the estimated sample size (ESS) from the log files produced by BEAST and the ESS of all parameters exceeded 200. The maximum-clade-credibility tree topologies were calculated from the posterior distribution and TreeAnnotator version 1.8.4 was employed to summarise the trees and view them in FigTree v.1.4.3. Finally, DnaSP v.4.0 was used to determine the number of haplotypes, haplotype diversity and segregation sites in the population (Rozas et al. 2003) and the relationships between haplotypes were illustrated in the form of minimum-spanning network (MSN) through Network 4.6.1.2.

The NJ phylogeny tree (Fig. 2) was generated using the Kimura-2 parameter with 1000 bootstrap replications. All Siamang samples remained in the monophyletic clade with a 100% bootstrap value. Siamangs formed two clades. The first clade represented the monophyletic group of Siamangs from Peninsular Malaysia supported by a 97% bootstrap value. Amongst this clade, ECSS644 and NASSC671 formed a monophyletic subclade with 86% bootstrap support that branched out from other Peninsular Siamangs. The remaining Peninsular Siamangs formed another monophyletic subclade supported by a high bootstrap value (97%). Siamangs from Fraser’s Hill (ARSSC529) had the closest relationship with Siamangs from Genting Highlands (ECSSC643), although other individuals from Genting Highlands (ECSSC643, ECSSC644, ECSSC674, NASSC670 and NASSC671) were geographically sympatric.

Figure 2.  

Neighbour-joining (NJ) phylogeny tree. Bootstrap values are indicated above the branch.

The second clade depicted the separation of the Siamang monophyletic subclade formed by DQ862110 and DQ862112.1 from other Sumatran Siamangs distinguished by a 67% bootstrap value. The remaining individuals (DQ862114.1, DQ862115.1, DQ862116.1 and DQ862117.1) formed a monophyletic subclade supported by a high bootstrap value of 98%. The tree topology also depicted two major clades: Hylobates + Hoolock (supported with an 89% bootstrap value) and Symphalangus + Nomascus (supported by a 44% bootstrap value). The Nomascus clade was formed by N. siki, N. leucogenys and N. gabriellae with a 99% bootstrap value. The Hoolock clade supported by a 94% bootstrap value was represented by H. hoolock. Within clade Hylobates, further subclades were defined. The first subclade characterised the klossii–agilis relationships supported by a 73% bootstrap value. The second subclade represented three Hylobates species, which were H. muelleri, H. pileatus and H. lar, supported by a high bootstrap value of 97%. H. muelleri was the first to diverge from two other species, followed by splitting of H. pileatus from H. lar population.

The MP tree (Fig. 3) agreed with all Siamang groups in the same clade supported by a 99% bootstrap value. The differentiation between Peninsular and Sumatran Siamangs shown by MP analysis was parallel to the NJ tree, except for the subclade formed by Siamang members. The two clades structured within Symphalangus represented the Peninsular and Sumatran clades, which were supported by 97% and 92% bootstrap values, respectively. In the Peninsular Siamang clade, ECSSC643 was the earliest to branch out from the clade comprising ECSS644, ARSSC529, NASSC670, NASSC671 and NASSC674 (supported by a 61% bootstrap value). In the Sumatran clade, the branch topology within the members was equal to that of the NJ tree, but was supported with a higher bootstrap value. An 80% bootstrap value had grouped DQ862112.1 and DQ862101.1 together and formed a sister subclade to the remaining Sumatran Siamangs that the subclade supported by 97%. The topology pattern of hylobatids positively reflected the NJ tree where Nomascus formed a monophyletic clade with Symphalangus (53% bootstrap value) and the Hoolock clade together with Hylobates (88% bootstrap value). Each genus further branched into its own monophyletic clade: Nomascus, Symphalangus, Hoolock and Hylobates supported by 99%, 99%, 86% and 99% bootstrap values, respectively. For the Hylobates clade, H. klossii and H. agilis remained as the basal species.

Figure 3.  

Maximum parsimony (MP) phylogeny tree. Bootstrap values are indicated above the branch.

The ML tree (Fig. 4) supported the divergence between the two populations of Symphalangus. Within the Symphalangus major clade, two clades were formed, which were the Peninsular Malaysia Siamang clade (98% bootstrap value) and Sumatran Siamang clade (83% bootstrap value). The Sumatran clade had a typical branching pattern to the NJ and MP trees. However, in the Peninsular clade, four different branches were formed. The first branch, represented by NASSC670 and NASSC674, grouped together with a 58% bootstrap value. The second branch cladding of NASSC671 and ECSSC644 with a 69% bootstrap value, while ARSSC529 and ECSSC643 branched individually. The ML tree illustrated almost a typical branching pattern as NJ and MP analysis for hylobatids unless for the Nomascus position. A 98% bootstrap value supported the clustering of all Hylobatidae members. In contrast to NJ and MP, the Nomascus clade (52% bootstrap value) diverged from the Symphalangus + Hoolock + Hylobates clade (95% bootstrap value). Later, Hylobates and Hoolock were grouped in the monophyletic clade with a bootstrap value (95%) separated from the Symphalangus clade. This clade later diverged into two subclades represented by Hylobates and Hoolock, both of which were supported by 100% and 87% bootstrap values, respectively.

Figure 4.  

Maximum Likelihood (ML) phylogeny tree. Bootstrap values are indicated above the branch.

The BI tree (Fig. 5) agreed with the topological branch with ML analysis. A high posterior probability (PP) value (1.00) supported the grouping of all Hylobatidae genera in the same clade. However, in the Hylobatidae clade, Nomascus was recognised as a basal genus with 0.84 PP, which diverged the earliest and formed a sister clade to Symphalangus, Hoolock and Hylobates. Three remaining genera were grouped in the same clade with a 1.00 PP value. Then, it was followed by divergence of Symphalangus forming its own clade supported by 1.00 PP. Hoolock and Hylobates still branched together with 1.00 PP and later separated and formed a different clade with 0.95 and 1.00 PP values, respectively. Similar to other phylogenetic trees, BI showed absolute distinction of two clades of Symphalangus with high PP support. The first Siamang clade represented the Peninsular Malaysia Siamang population supported by 1.00 PP and the second clade comprised all Sumatran Siamangs with a 0.98 PP value. Within the Peninsular Malaysian clade, ECSSC644 and NASSC671 had the closest genetic relationships as a result of the NJ tree. For the Sumatran clade, the first individual to diverge was DQ862101.1, followed by DQ862112.1.

Figure 5.  

Summary of Bayesian inference and molecular clock phylogenetic tree. Values above the branch represent the posterior probability for BI. Values on the nodes represent the divergence times of the clade.

The estimated evolutionary divergence dates were illustrated and supported by Bayesian analysis, as shown in Fig. 5. The analysis indicated that S. syndactylus diverged approximately 11.1 million years ago (MYA), splitting the species from the Nomascus clade. Later on, S. syndactylus was separated into two clades. The molecular clock separation between the Peninsular Malaysia and Sumatran populations of S. syndactylus occurred approximately 3.12 MYA. Splitting event occurred much earlier within the Peninsular Siamang population (1.92 MYA) than in the Sumatran Siamang population (1.85 MYA). The tree genus group of the family Hylobatidae (Symphalangus, Hoolock, Hylobates and Nomascus) diverged from P. pygmaeus, a split that was estimated to occur 32.17 MYA. Nomascus radiation occurred approximately 3.86 MYA, splitting from Symphalangus, Hoolock and Hylobates at 18.39 MYA. Hoolock and Hylobates diverged from Symphalangus 11.1 MYA. Then, this was followed by speciation of Hylobates in several events at 4.29 MYA. Amongst Hylobates, the clade formed by H. lar, H. pileatus and H. muelleri was the earliest group to branch out from the other Hylobates clades, which diverged approximately 3.65 MYA. Another lineage branched out and separated the clade (comprising two species, namely, H. agilis and H. klossii) approximately 3.21 MYA.

MSN was used to visualise the associations between haplotypes, based on the haplotype analysis (Fig. 6). The network analysis revealed that the Peninsular Malaysia Siamang has a different haplotype from the Sumatran Siamang. Symphalangus contained six unique haplotypes: Hap_1 to Hap_3 characterised the haplotype of the Peninsular Siamang population and Hap_4 to Hap_6 represented the Sumatran Siamang population. The analysis depicted Hap_1 (Peninsular Malaysia Siamang) as the connecting point to Sumatran Siamangs with at least eight mutational steps as the fewest steps to indicate the relationship between these two populations. The MSN also revealed three haplotypes of Nomascus (Hap_26, Hap_27 and Hap_28), three haplotypes represented Hoolock (Hap_29, Hap_30 and Hap_31) and 18 haplotypes of Hylobates (Hap_12 to Hap_25) varied according to the species. The result depicted the fewest mutational steps (9) illustrated between Symphalangus and Nomascus haplotypes compared to other genera. Thus, we drew the closest relationship between Symphalangus and Nomascus amongst hylobatids. Hoolock and Hylobates haplotypes were adjacent to each other with at least 21 mutational steps.

Figure 6.  

The minimum spanning network (MSN) illustrating the relationships between the Symphalangus population and other hylobatids. Each circle represents a haplotype.

Phylogenetic trees indicate the genetic relationships amongst the members of Hylobatidae. Two different branch topologies of the hylobatid clade were obtained to explain the genetic relationships amongst these genera. In the NJ and MP trees, Nomascus and Symphalangus branched together, forming a sister clade to the Hylobates–Hoolock clade. For ML and BI trees, alternatively exclude Nomascus from being in a clade together with Symphalangus, which formed a monophyletic clade with Hoolock and Hylobates. Moderate support in bootstrap value was recorded for the topology of the Nomascus–Symphalangus clade, but there was a high support for the Hylobates–Hoolock clade formation. However, the PP value showed positive support for each topography clade formed amongst the hylobatids. The Symphalangus position is prone to form a clade together with Hylobates + Hoolock in bootstrap analysis, specifically NJ and MP trees. This can be explained by the fact that the earliest primatologists classified Symphalangus as a subgenus of Hylobates (Kloss 1929). Regardless, Kloss (1929) discovered that the dwarf gibbon is the absolute intermediate species between ordinary gibbons and siamangs, excluding Symphalangus as an out-group and grouping all gibbons into one genus. This is parallel to this study as Symphalangus later diverged from the Hylobates clade. In addition, the results showed high support of bootstrap and PP values for the branching of Hylobates and Hoolock in the monophyletic clade (without Symphalangus and Nomascus).

Hoolock gibbons were previously classified as members of the genus Bunopithecus (Prouty et al. 1983) and Hylobates and Hoolock once formed a monophyletic clade, as proposed by Groves (2001) and Brandon-Jones et al. (2004). Mootnick and Groves (2005) reclarified the Hoolock gibbons as having different karyotypes from Hylobates, renaming the genus as Hoolock. Thus, this study supports the distinction of Hoolock from Hylobates as mentioned by a previous researcher. Groves (2001) acknowledged Nomascus as a subgenus of Hylobates, which contradicted our study, which suggested Nomascus as the basal genus of hylobatids. However, because of contrasting support between bootstrap and PP values in this study, the relationship of these four genera to each other remains obscure as no previous research revealed typical branching patterns (Carbone et al. 2014). The relationships amongst the four genera were previously unresolved, despite the efforts of experts, such as Chan et al. (2010) and Matsudaira and Ishida (2010), who successfully constructed the evolutionary tree that supported Nomascus as the pioneer genus (without adding Hoolock), but the structure collapsed when all genera were analysed (Roos 2016).

However, the current classification of Nomascus as a basal taxon of Hylobatidae has been accepted, as illustrated in the Maximum Likelihood tree by Israfil et al. (2011), which is consistent with our study with stronger support as more samples of hylobatids were included. The groupings of different Hylobates species (H. lar, H. muelleri, H. agilis, H. pileatus and H. klossii) are congruous for all four phylogenetic trees. According to Roos (2016), the branching pattern amongst Hylobates species may be controversial, but there is a common agreement on H. pileatus forming the basal lineage of this group. However, our study indicates contradictory results, where H. klossii and H. agilis (in a monoclade) diverged the earliest, followed by the divergence of H. muelleri, H. lar and H. pileatus. The relationships between the species comprising clades H. klossii–H. agilis and H. lar–H. muelleri remain unresolved as both bootstrap and PP values indicated low support, which agrees with the results published by Carbone et al. (2014). This ambiguity is common in situations where sympatric species may migrate, leading to hybridisation (Roos 2016).

In terms of divergence times, the Hylobatidae mitochondrial D-loop sequences estimate divergences in this study is slightly near to those proposed by Israfil et al. (2011) and these occurred approximately 19.7–24.1 MYA. However, the divergence of Syndactylus in this study did not match the divergence date of Hylobates syndactylus from Nomascus, estimated to be approximately 5.6–7.2 MYA as proposed by Israfil et al. (2011). This mismatch in divergence age might be due to the high number of samples and the type of locus used. This study only covered a part of the mitochondrial region, which is the D-loop, compared with Israfil et al. (2011), who used the whole mitochondrial and sex-linked gene. However, the D-loop region is still considered more relevant to the phylogenetic study of rapidly changing gene that radiate over a short period of time (Abdul-Latiff et al. 2019). In addition, different topological branches elucidate the different time ranges, as Israfil et al. (2011) proposed Hoolock as the basal genus instead of Nomascus.

The results connotated Siamangs from Peninsular Malaysia as the basal taxon of radiation because the population divergence date was much earlier than that of Sumatran Siamangs. This is comparable to studies on Trachypithecus cristatus (Roos et al. 2008) and Hylobates (Whittaker et al. 2007, Israfil et al. 2011) that the mainland population is considered as ancestral stock, while the southern population is derived species. Thus, it can be presumed that the speciation of Siamangs occurred first in Peninsular Malaysia, after which Siamangs migrated and diversified in Sumatra. As argued by Chatterjee (2006) and Chatterjee et al. (2009), Hylobates radiated from mainland Southeast Asia southwards and eastwards (via Peninsular Malaysia and Sumatra) ~ 8.3 MYA. During the Pleistocene, the early diversification of hylobatids was followed by a period of rapid speciation, notably within the concolour and lar gibbon genera (Carbone et al. 2014). Symphalangus was claimed to have diverged from the hylobatid lineage through southward expansion (Chatterjee et al. 2009) as it can be traced in the west of the Mekong River (Thinh et al. 2010).

Past scenarios might have caused the small population of different species to thrive in small refugees during the Pliocene (~ 3.17–0.31 MYA), which consequently drove the species apart and formed different populations (Wang 1994). At this age, Presbytis and Macaca had colonised Sumatra from the Malay Peninsula throughout the Late Pliocene and Early Pleistocene, where cold climate occurred at approximately 1.8 Ma (Harrison et al. 2006). Later, sea level falls formed a route that connected the Javan and Peninsular Malaysia mainland first, without connecting to Sumatra and Borneo during the Pleistocene (~ 1.25 MYA) (Hall et al. 1998). This was followed by drastic falling of sea levels (Gorog et al. 2004) causing the mainland to expand and being connected to Sumatra (Roos et al. 2008). The connection between these two regions via a land corridor has allowed the radiation of Hylobates (Chatterjee 2006) and T. cristatus (Roos et al. 2008) from the mainland to the southward region. Thus, we deduced the baseline of the same land bridge that caused the migration of Siamangs from Peninsular Malaysia to Sumatra.

Compared to this study, splitting of Siamangs occurred 3.12 MYA, the age of which is in line with the late Pliocene era when Siamang from Peninsular might move along with Presbytis and Macaca migrating from Malay Peninsular to Sumatra where Peninsular Malaysia and Sumatra were already connected to each other. It is possible Siamang that colonised Peninsular Malaysia migrate after a vicariance event such as the rise and fall of sea level (Abdul-Latiff et al. 2014) resulting in the existence of Siamang in Sumatra. Thus, we proposed the Peninsular Malaysia Siamang (S. s. continentis) positioned in the mainland as the ancestor and the Sumatran Siamang (S. s. syndactylus) that diverged later as derived species.

On the basis of the previous result, the taxonomic position of the two proposed subspecies is accepted because the D-loop analysis illustrated significant subspecies distinction. Raffles (1821) was the first to describe Symphalangus, based on the sample found in Sumatra. At that point, he concluded that the Siamang was a single species and no further classification was made. None of the fossil records was reviewed to clarify the ages of the Siamang population in Peninsular Malaysia and Sumatra to confirm which population existed the earliest. However, Thomas (1908) compared Sumatran skulls with British Museum skull specimens to explore the proclivity for minor variations in Peninsula specimens. This result indicates that the skull size of the Peninsula Siamang is smaller than that of the Sumatran Siamang. Thomas (1908) referred to Siamangs from Peninsular Malaysia as S. s. continentis and Sumatran Siamangs as S. s. syndactylus.

The subspecies distinction was embraced by another morphological description, recorded by Kloss (1929). In terms of colouration, the face of S. s. syndactylus is loosely covered with white fur and faintly tinted with brown, whereas S. s. continentis has less pale hair on the face and fur that is frequently less thick than that of Sumatran Siamangs, which invade colder climates (Kloss 1929). Both populations are physically comparable in terms of size, although S. s. continentis seems larger than S. s. syndactylus (Kloss 1929). The morphological characteristics in terms of skull size (Thomas 1908) and body size (Kloss 1929) for each subspecies seem to be contradicting. It cannot be confirmed that the skull size of S. s. syndactylus is larger than that of S. s. continentis as described by Thomas (1908) because no evidence has been provided on which sex of Siamangs for each subspecies has been used for comparison. This is because different skull sizes have been documented for different sexes as indicated by island rule (Welch 2009). Hooijer (1960) argued that the size and proportions of the skulls of S. s. syndactylus are remarkably varied in males and females. The average male skull is larger than the female skull with a thinner brain case, more projected orbiting and deeper postorbital constriction. Symphalangus s. continentis also has sexual variations of the skull that are comparable to those of S. s. syndactylus (Hooijer 1960).

Hylobatids have been postulated to be a mechanism for decrease in size or dwarfing (Pilbeam 1996). Body size changes in response to a colonising continental species may be adaptable to an island depending on the island mass and its existing animals and vegetation (Welch 2009, Burness et al. 2001). Hylobates have evolved with specific features, such as a reduced body size and well-developed fore-limb suspensory locomotion (Reichard et al. 2016). However, Siamangs are enigmatic as their size is one-half to two times larger than that of other gibbons. More recently, other morphological features, such as nose and toe structures (Fig. 7), were reported by Mootnick (2006). Symphalangus s. syndactylus has a shallow and long tapering point pattern over its nostrils and webbing between the fourth and fifth toes, whereas S. s. continentis has a wider nose and second and third toes extending up to the distal interphalangeal joint (Mootnick 2006).

Figure 7.  

Symphalangus s. continentis (a and b) and S. s. syndactylus (c and d) nose structure and toe structure (Mootnick 2006).

These features add evidence for classifying the Peninsular Malaysian Siamangs as S. s. continentis and Sumatran Siamangs as S. s. syndactylus. Other researchers have continued to classify Siamangs according to morphological appearance as the molecular approach requires genetic data that were previously unavailable, thus hindering molecular studies at the species level (Roos 2016). However, this study successfully reclarified the existence of two subspecies of Symphalangus.

The primary threat to Siamangs is forest fragmentation, which limits their ranging area and the availability of food sources (Marshall 2009). As with other primates, the Siamang habitats in Malaysia are being disrupted either directly or indirectly by human activities. For instance, in a survey of 14 forest fragments that include the Virgin Jungle Reserves and adjacent disturbed forests in Peninsular Malaysia, S. syndactylus was present in two out of six large fragments and absent from all small- and medium-sized forest fragments (Laidlaw 2000). Habitat fragmentation leads to an increased risk of predation and hunting of this species, resulting in an increased probability of local extinction. To protect this species, both in situ and ex situ conservation are needed. Molecular studies on Siamangs to reveal taxonomical classification have shown that Peninsular and Sumatran Siamangs have different gene pools. It is important to acknowledge the genetic relationship up to the subspecies level as this may affect the decision for translocation as similar gene pools should be grouped together to encourage survival. Our results confirmed subspecies classifications within S. syndactylus and this information will be useful in identifying the most suitable region for Siamang translocation for the purpose of conservation.