Evolutionary pattern of Macaca fascicularis in Southeast Asia inferred using Y-chromosomal gene
BMC Ecology and Evolution volume 21, Article number: 26 (2021)
We analyzed a combined segment (2032-bp) of the sex-determining region and the testis-specific protein of the Y-chromosome (Y-DNA) gene to clarify the gene flow and phylogenetic relationships of the long-tailed macaques (Macaca fascicularis) in Southeast Asia. Phylogenetic relationships were constructed using the maximum likelihood, Bayesian inference, and the median-joining network from a total of 164 adult male M. fascicularis from 62 localities in Malaysia, including sequences from the other regions from previous studies.
Based on Y-DNA, we confirm the presence of two lineages of M. fascicularis: the Indochinese and Sundaic lineages. The Indochinese lineage is represented by M. fascicularis located northwards of the Surat Thani-Krabi depression region and is introgressed by the Macaca mulatta Y-DNA. The Sundaic lineage is free from such hybridization event, thus defined as the original carrier of the M. fascicularis Y-DNA. We further revealed that the Sundaic lineage differentiated into two forms: the insular and the continental forms. The insular form, which represents the ancestral form of M. fascicularis, consists of two haplotypes: a single homogenous haplotype occupying the island of Borneo, Philippines, and southern Sumatra; and the Javan haplotype. The more diverse continental form consists of 17 haplotypes in which a dominant haplotype was shared by individuals from southern Thai Peninsular (south of Surat Thani-Krabi depression), Peninsular Malaysia, and Sumatra. Uniquely, Sumatra contains both the continental and insular Y-DNA which can be explained by a secondary contact hypothesis.
Overall, the findings in this study are important: (1) to help authority particularly in Malaysia on the population management activities including translocation and culling of conflict M. fascicularis, (2) to identify the unknown origin of captive M. fascicularis used in biomedical research, and; (3) the separation between the continental and insular forms warrants for the treatment as separate management units.
The most common and abundant Macaca species in Malaysia, Macaca fascicularis (Raffles, 1821) (also known as long-tailed macaques, crab-eating macaques, and cynomolgus macaques), are widely distributed in nature and occupy vast areas of mainland Southeast Asia (SEA) (Bangladesh, Thailand, Cambodia, Vietnam, Laos, Myanmar, Peninsular Malaysia and Singapore), the Greater and Lesser Sunda Islands (Indonesia, Brunei, and the Malaysian Borneo), and the Philippines, extending across 30° of latitude and 35° of longitude . This species is highly adaptable to disturbed environments and secondary forests, and commonly populate low elevation habitat types favouring seashores and mangrove forests, riverbanks, and swamp forests [2, 3]. In Peninsular Malaysia, they are widespread and populate areas in sympatry with human settlements [4, 5]. At the same time, in the Malaysian Borneo (Sabah and Sarawak), they are distributed throughout the lowlands populating the coastal regions .
Macaca fascicularis like all the other species within the genus Macaca demonstrates a high level of male-mediated dispersal while females are sedentary in nature [7, 8]. This implies that the transfer of genes between groups and populations occurs almost exclusively through the dispersal of males . Moreover, the dominant male exhibits a polygynous mating system, which results in a homogenous gene pool of offspring. Because of this characteristics, particular M. fascicularis groups or populations would consist of homogenized Y-chromosomal DNA (Y-DNA) .
Y-DNA is uniparentally inherited and as a haploid genetic system, is passed down paternally from male adults to their male progenies. Due to the lack of recombination, Y-DNA has one-fourth of the effective population size (Ne) of that of nuclear DNA [9, 10]. Additionally, the polygynous mating behaviour of dominant males further lowers the Ne of Y-DNA loci [9, 11,12,13,14]. Nevertheless, the paternal inheritance and the highly dispersal nature of male macaques mean that the information from Y-DNA is important to infer the male-mediated gene flow and to distinguish between contemporary gene flow and ancestral polymorphism .
The sex-determining region (SRY) and the testis-specific protein (TSPY) are two unlinked gene segments within the Y-DNA which have been used to infer male-mediated gene flow in primates, particularly in Cercopithecinae [9, 15,16,17,18,19,20,21,22,23]. Specifically, Tosi, Morales and Melnick  and Bunlungsup et al.  used both the SRY and TSPY segments coupled with mitochondrial (mtDNA) to investigate the contact zone between M. fascicularis and M. mulatta in the Indochinese region while Tosi, Morales, and Melnick [9, 16] used a small number of representatives from various region of SEA to investigate the phylogeography of M. fascicularis.
The widespread distribution of M. fascicularis makes it a prime candidate for the study of phylogeography and evolutionary studies. During the Quaternary Periods (Pleistocene and Holocene), the currently separated landmasses in SEA were intermittently joined and formed Sundaland that permits the dispersal of land mammals, including primates [24,25,26]. The climatic oscillations, especially during the Pleistocene, including major volcanic eruptions of Mount Toba, would have induced vegetation changes causing the expansions or contractions of habitat that influenced the dispersal of animals [27,28,29]. Moreover, major river systems, lakes, and mountain ranges may also have acted as a barrier to gene flow [26, 28].
In Peninsular Malaysia, the rising number of cases of human-Macaca conflict, particularly involving M. fascicularis has prompted the Department of Wildlife and National Parks (DWNP) to implement several population control measures including translocations and culling [5, 30]. Conflict M. fascicularis mainly male individuals, have been translocated from conflict urban areas to other locations with a less inhabited human population . Therefore, in this study, by investigating the SRY and TSPY segments, information on whether Y-DNA of M. fascicularis consists of continuously related or segregated haplotypes would provide information on whether the translocated individuals posed any risks of contamination to the existing population. Furthermore, taking into account the historical zoogeography and phylogeography of SEA, investigation into the male-mediated gene flow would also provide insight into the hypothesis on the existence of historical or current barrier as well as evidence of recent gene flow among the M. fascicularis populations in SEA, particularly in Malaysia.
Haplotype mapping, sequence characteristics, and genetic diversity
A total of 654- and 1370-bp of sequence lengths were obtained from the sex-determining region (SRY) and the testis-specific protein (TSPY) segments, respectively, from all 164 male M. fascicularis samples. Since both segments are closely linked on the Y-chromosome, the datasets from 287 sequences were combined to produce an alignment of 2032-bp (including indels). All the sequences in this study were registered with the National Center for Biotechnology Information (NCBI) and were given accession numbers, as shown in Additional file 1: Table S1.
Considering all the M. fascicularis sequences (N = 232), 20 haplotypes were observed. Out of the 20 haplotypes, the Malaysian M. fascicularis were represented by 18 haplotypes in which 16 were unique to Peninsular Malaysia (Additional file 1: Table S1; Fig. 1a). Peninsular Malaysia constitutes 17 haplotypes where Haplotype 1 is the most dominant and can be found in all the states and major islands (except on Tioman Island). Interestingly, Haplotype 1 is also shared with all nine M. fascicularis sequences from the southern Thai Peninsular (south of Surat Thani-Krabi depression), four sequences from Sumatra (Fig. 1b), and the introduced population to Mauritius (N = 10). The Indochinese M. fascicularis sequences from Vietnam, Cambodia, Thailand (north of Surat Thani-Krabi depression) formed a single haplotype, Haplotype 20 (Additional file 2: Table S2; Fig. 1b) and were uniquely shared with M. mulatta from Burma and South China. The Malaysian Borneo sequences (Sarawak, Labuan Island, and Sabah) formed a single haplotype (Haplotype 18) shared with sequences from Kalimantan and the Philippines including two sequences from southern Sumatra while Java was represented by Haplotype 19.
Based on the haplotype mapping, the M. fascicularis sequences were assigned into three geographical groups (except for sequences from Mauritius) for further sequence characterization: (1) Indochinese (Haplotype 20); (2) Continental Sundaic consisting of sequences from the Thai peninsula, Peninsular Malaysia, and Sumatra (Haplotype 1–17), and; (3) Insular Sundaic (Haplotype 18–19). Four indels were observed in M. fascicularis, distinguishing between the Indochinese and the Sundaic (continental and insular) haplotypes where two events of; (1) three-bp indels (ACA) at nucleotide position (np) 596–598 and (2) one-bp thymine (T) indel at np 974 occurred (Table 1). Furthermore, within the Sundaic lineage, a single transversion mutation (adenine to thymine) at np 159 was observed on the SRY segment which distinguished the insular Sundaic from the continental Sundaic haplotypes.
Table 2 summarized the sequence characterization, genetic diversity indices, genetic structure, and differentiation among the geographical groups of M. fascicularis. The continental Sundaic is the most variable group displaying the highest number of variable sites (VS) and parsimony-informative sites (PIS). Estimation of the Y-DNA nucleotide diversity (π) among the Sundaic lineage revealed that the continental Sundaic haplotypes had a significantly higher π (4.4 × 10–4) compared to the insular Sundaic (0.4 × 10–4). The genetic diversity for the Indochinese lineage was incalculable since only a single haplotype was observed. The pairwise fixation index (FST) calculated among the groups revealed that the Indochinese M. fascicularis is highly divergent from both the Sundaic forms (FST of 0.90 and 0.99 to the continental and insular forms, respectively) while within the Sundaic lineage, the insular and continental forms are moderately differentiated from each other (FST of 0.69). The genetic distances among the haplotypes range from 0.05 to 0.35% (data not shown). Among the geographical groups, the genetic distance between continental Sundaic and insular Sundaic is at 0.15%, while the genetic distances between Indochinese with the continental Sundaic and insular Sundaic is at 0.27% and 0.17%, respectively.
Phylogenetic relationships and estimation of divergences
The combined SRY and TSPY phylogenetic relationships constructed by using the maximum likelihood (ML) (− lnL = 3809.0846) and Bayesian inference (BI) methods produced similar topologies and thus were represented by the BI tree (Fig. 2). In general, the tree grouped the Macaca into their respective species groupings consistent with the four major lineages of Macaca: (1) sylvanus-silenus, (2) sinica, (3) arctoides, and (4) fascicularis.
Within the fascicularis species group, the tree formed two distinct clades: (1) a clade consisting of M. cyclopis, M. fuscata, M. mulatta, and the Indochinese M. fascicularis, and (2) M. fascicularis haplotypes of Sundaic origin, with high support value (99% and 100% for ML and BI, respectively). Furthermore, the Sundaic M. fascicularis lineage further bifurcated into the continental and insular forms but with low support values (below 50%). Similarly, the median-joining network (MJN) (Fig. 3) produced similar groupings, as observed in the phylogenetic tree. Moreover, MJN showed that the continental Sundaic haplotypes were all derived from the insular form, particularly from Haplotype 18. Furthermore, all the other haplotypes within the continental Sundaic haplotypes were derived from the dominant Haplotype 1.
The divergence date estimations calculated by calibrating the BI using the proposed date for the last common ancestor (LCA) of Macaca at 5.5 million years ago (mya) are shown in Fig. 2. The fascicularis species group (cyclopis, fascicularis, fuscata, and mulatta) shared the LCA at around 2.35 mya. Next, the bifurcation between the continental and insular Sundaic forms was estimated to have occurred at ~ 2.00 mya. After that, the continental form shared the LCA at around 1.65 mya while the insular form shared the LCA at around 0.71 mya.
By analyzing 164 male samples from 62 localities throughout Malaysia (including five major islands), this work represents the most comprehensive effort to infer male-mediated gene flow of M. fascicularis in Malaysia based on the TSPY and SRY segments of the Y-DNA gene. Combining this dataset with the sequences of M. fascicularis from the other regional populations, a complete and broad analysis of the Y-DNA gene flow were able to be estimated, covering the entire range of the species in SEA. In general, the combined segments of the SRY and TSPY in this study were able to provide several paternal insights into the gene flow, phylogenetic relationships at the population level as well as divergence time estimation of M. fascicularis.
In this study, the phylogenetic relationship of Macaca based on the Y-DNA supported the four major species groupings congruent with the morphological classification by Delson : sylvanus, silenus, sinica, and fascicularis species groups. Although supported with low bootstrap values (except for fascicularis species group which is supported by high bootstrap values of 99 and 100% for BI and ML, respectively), all Macaca species were clustered into their respective species groups. Similarly, previous nuclear DNA studies also supported the four major groupings [9, 16, 32,33,34].
Within the fascicularis species group, two major clades were observed: (1) a species complex group comprising of M. cyclopis, M. fuscata, M. mulatta, and the Indochinese M. fascicularis and (2) M. fascicularis haplotypes of Sundaic origin. Both clades shared an LCA at ~ 2.35 mya, which is in proximity to the previous estimation of 2.30 mya by Tosi, Morales, and Melnick . This estimation coincides with the multiple major climate cooling events during Late Pliocene and early Pleistocene occurring between 2.8 and 2.4 mya [35, 36] that influenced the paleoenvironments of SEA [24, 25]. The cold and dry climate during the glacials gave rise to savannah-like vegetation or semi-open woodland [25, 27, 29, 37], causing restriction of dispersal for primates, thus could lead to the separation between both clades.
Hybridization of M. mulatta Y-DNA into M. fascicularis in the Indochinese region
The haplotype sharing between M. mulatta (N = 12) and M. fascicularis (N = 34) represents the male-mediated gene flow of M. mulatta into M. fascicularis populations from the Indochinese region. Several authors have previously reported on the possible occurrence of the hybridization events between the Indochinese M. fascicularis and M. mulatta [9, 15, 16, 23, 38,39,40]. The extent of the contamination of M. mulatta Y-DNA into the Indochinese M. fascicularis population is not only restricted to the contact zones between both species but possibly more widespread across the Indochinese region [23, 38, 39].
In this study, we confirmed the absence of the male-mediated gene flow of M. mulatta into M. fascicularis from the Sundaic region, which is restricted at the Surat Thani-Krabi depression region, as shown in Fig. 1b. This is similar to the separation observed between the southern and northern pig-tailed macaques . It is possible that the Surat Thani-Krabi depression played a significant role as a barrier that limits the downwards dispersal of the contaminated Indochinese M. fascicularis into the Sundaic lineage. Thus, based on the current finding, we defined the Sundaic M. fascicularis as the original carrier of Y-DNA in M. fascicularis.
The region of Isthmus of Kra and Kangar-Patani marks the biogeographical barrier that separated the Indochinese and Sundaic zoogeographical regions . These regions display the shift from mixed deciduous forest to the north to wet seasonal rainforest type to the south . Situated along the modern Thailand–Malaysian border (see Fig. 1b), these biogeographical barriers has been documented to limit the gene flow between the Indochinese M. fascicularis from their Sundaic conspecifics [15, 23, 39, 44] as well as in other vertebrate faunal species [43, 45,46,47,48,49,50]. Therefore, to measure the extent of the M. mulatta male hybridization into M. fascicularis, further research to screen M. fascicularis from this area are needed.
Sundaic M. fascicularis: continental and insular forms
The Sundaic M. fascicularis lineage, which is the original carrier of Y-DNA in M. fascicularis, was further separated into continental and insular forms and shared the LCA at about 2.00 mya. The continental form comprised of haplotypes from southern Thai Peninsular (south of Surat Thani-Krabi depression), Peninsular Malaysia, and Sumatra, including the introduced Mauritian population, while the insular form consisted of haplotypes from southern Sumatra, Java, Borneo, and the Philippines. Likewise, previous molecular studies have identified both forms using mtDNA and nuclear gene [15, 23, 51,52,53,54,55,56]. Furthermore, the MJN suggested that the insular form was the predecessors of the continental form (Fig. 3). This result would be consistent with the hypothesis that the ancestors of M. fascicularis evolved in insular SEA [31, 57, 58]. Wherever and whenever this species originated, the earliest M. fascicularis-like fossil was discovered in Java, suggesting that they have already inhibited the Sundaic region as early as 0.9–1.0 mya .
After diverging out from the insular form, the continental form diverged into 17 haplotypes, which all shared the LCA at around 1.65 mya (Fig. 2). The dominant Haplotype 1 characterized this form, which is shared with 103 individuals and is distributed as far north as Wat Suwankhuha, of the Thai Peninsular, downwards to entire Peninsular Malaysia until Sumatra of Indonesia. According to the MJN (Fig. 3), all the other continental haplotypes were derived from this dominant haplotype. Moreover, within Peninsular Malaysia, the distribution of haplotypes (Fig. 1a) further showed the segregation of several unique haplotypes that can be found only at the north-western region (Haplotype 4, 6–8, 11–14) as well as the east coast region (Haplotype 15–17). In contrast, the central-southern regions are relatively connected with no unique haplotype defining the region.
Compared to the 17 observed haplotypes within the continental form, only two haplotypes were observed within the insular form, which shared the LCA at around 0.71 mya (Fig. 2). M. fascicularis from southern Sumatra, Borneo, and the Sibuyan Island of northern Philippines were all represented by a single dominant insular haplotype, Haplotype 18. The island of Java, on the other hand, was represented by Haplotype 19. The small number of haplotypes observed shows that the insular form consists of significantly low variability and intact Y-DNA. The highly dispersal nature of males would move and homogenize Y-DNA across populations  and therefore contribute to this highly intact Y-DNA across the insular region.
Prolong and recurring connectivity that existed in the past could have facilitated the dispersal of males between the insular regions (Sumatra, Java, Borneo, and the Philippines). Several episodes of sea-level fluctuations in Sundaland during the Pliocene up until recently during the Holocene around seven thousand years ago (kya) would permit for such connectivity with some as low as − 120 m below the present level [25, 59]. In addition to the connectivity, a possible ancient bottleneck could explain for the intact Y-DNA in the insular region. The ancestors of the insular form could have been forced into a common refugium postulated to have existed in northern Borneo [24, 25, 29, 60, 61]. Over time, only a single dominant haplotype could have persevered. Alternatively, overhunting and exploitation by early human could also severely have reduced the effective Y-DNA gene pool. Remains of M. fascicularis-like species were discovered in cave settlements in Java (~ 0.9 to 1.0 mya) and Borneo (~ 30 to 40 kya)  giving evidence to early exploitations by human. Nevertheless, the low number of samples and available sequences (N = 27) could hinder the possibility of detecting other haplotypes within the insular form. Therefore, to further measure the diversity of the insular form, further research utilizing more samples should be conducted on the insular form.
Uniquely, M. fascicularis from Sumatra exhibited both the continental and insular Y-DNA haplotypes (Fig. 1b; Additional file 2: Table S2). A secondary contact hypothesis could account for this condition in Sumatra with two possible explanation. The first explanation offered by Tosi and Coke  hypothesized that Sumatra was already populated by the insular Sundaic form, prior to the immigrating males of the continental form from the Malay Peninsula. Due to the narrow and shallow waterways on the Sunda Shelf separating the Malay Peninsula and Sumatra [25, 59], a recent dispersal of the dominant continental haplotype from the Malay Peninsula is postulated as being in the process of colonizing the native insular form.
In contrast, a second explanation similar to the first secondary contact hypothesis takes into account the biogeographical history of both regions. Several authors have postulated the formation of several refugia during the Pleistocene, whereas the Malay Peninsula was surrounded by a belt of savannah-like vegetation [24,25,26, 63]. This implies that dispersal of the continental form occurred from a postulated refugium in Sumatra into the Malay Peninsula and the surrounding areas. Dispersal of the insular form from a refugium in western Java into Sumatra, which was already populated by the continental Sundaic haplotypes could explain for the presence of both forms in Sumatra.
Findings in this study could also assist the authorities particularly in Malaysia on the population management of conflict M. fascicularis. Prior information on the genetic structure is imperative before any measures of population control such as translocation programs or culling are to be conducted [64,65,66]. Thus, the haplotype mapping in Fig. 1a can be used by the managing authority to correctly manage the local populations without posing risks of contamination to the existing population or accidentally wiping out a unique haplotype due to culling.
Overall, Fig. 1b summarizes the Y-DNA mapping of M. fascicularis according to the lineage and forms observed. Evidently, M. fascicularis are distributed accordingly to the zoogeographical segregation of SEA. The single nucleotide polymorphism observed at np 159 which distinguishes between the continental and insular forms are important from the perspective of biomedical research. The identification of the genetic background of non-human primate model used in experimental research is important since variations could influence the results and repeatability of experiments [67, 68]. Thus, this finding could assist in the identification of the origin of the captive M. fascicularis used in biomedical research, at least paternally.
Finally, findings of this study suggest that the separation between the continental and insular forms warrants for the treatment as separate management units (MUs). MUs are lower ranked conservation units which recognize populations with significant divergence at nuclear or mtDNA and represent populations connected by such low levels of gene flow . The status of the Indochinese lineage, on the other hand, should be handled with caution and until the level of genetic hybridization of Indochinese M. mulatta into M. fascicularis populations are thoroughly assessed, this lineage should be managed separately from the Sundaic lineage.
Sample collection and GenBank sequences
Blood samples were collected from wild free-ranging conflict M. fascicularis by authorized personnel and veterinarians of DWNP (for samples from Peninsular Malaysia), Sabah Wildlife Department (for samples from Sabah), and Forest Department of Sarawak (for samples from Sarawak). The euthanization of the conflict animals were also performed by authorized and qualified veterinarians of DWNP. In brief, the macaques were sedated intramuscularly using general anaesthesia (combination of ketamine, 5–10 mg/kg and xylazine, 0.2–0.4 mg/kg) before a lethal dosage of pentobarbital (equal or more than 60 mg/kg) were given intraveneously. A total of 164 samples of adult male long-tailed macaques from 62 localities in Peninsular Malaysia, Sarawak, and Sabah were collected and used in this study (Additional file 1: Table S1). Figure 4a and b summarizes the localities of the samples according to the states and major islands in Malaysia.
Sixty-eight sequences of M. fascicularis from other regional populations of were also included in the dataset to make up a total of 232 sequences (Fig. 4c; Additional file 2: Table S2) representing the most comprehensive sampling coverage of the species distribution range. In total, 287 sequences were used for analyses comprising of 20 Macaca species and two outgroup species (Papio hamadryas and Mandrillus sphinx). Both outgroup species were selected as representatives from the major branches within the Family Cercopithecidae.
DNA extraction, PCR amplification, and DNA sequencing
Total genomic DNA was extracted from blood samples using the QIAamp DNeasy® Blood and Tissue Kit using the protocol provided by the manufacturer (Qiagen, Germany). Two pairs of oligonucleotides were used to amplify the and the TSPY segments of the Y-DNA as designed by Rovie-Ryan et al. . In Old World monkeys the SRY segment is found on the short arm of the Y-DNA  while the TSPY segment is typically located on the long arm, close to the centromere region in the Y-DNA [70, 71].
PCR amplifications were conducted in 20 µl reactions using a GeneAmp® PCR System 9700 (Applied Biosystems, USA). Each PCR reaction consists of 1.0 µl of DNA template (~ 15 to 20 ng), PCR mixtures containing 4× Green GoTaq® Flexi Buffer (Promega, USA), 0.875 mM of MgCl2 (25 mM), 0.1 mM of each dNTPs (10 mM), 0.1 mM of each primer (10 mM), 1 unit of Taq Polymerase (5 unit/µl), and later added with ddH20 to a total of 20 µl of total reaction mixtures. The cycling profile for the amplification was as follows: a preliminary denaturation at 98 °C for 2 min followed by 45 cycles of 95 °C for 30 s, 55 °C for 60 s and 72 °C for 60 s. This was followed by a final extension period of 72 °C for 5 min before the samples were cooled to 10 °C. Cycle sequencing was done on an ABI PRISM®377 DNA Sequencer by a sequencing service provider (1st Base Laboratories Sdn. Bhd., Malaysia).
Sequence characterizations and genetic diversity indices
Multiple sequence alignments were done by using the program Geneious ver5.6 . Prior to further sequence analysis, the SRY and TSPY segments were combined as both segments are closely linked on the Y-chromosome and the partition homogeneity tests conducted in PAUP ver4  did not find significant differences in their evolutionary signal (P = 1.0) as also previously reported [9, 16].
Sequence characterizations including conserved sites (CS), VS, and PIS were examined by using MEGA ver6  while DnaSP ver5  were used to calculate the standard genetic diversity indices including the number of haplotypes (NHap), haplotype diversity (Hd), and π . Genetic distances were calculated on MEGA ver6 using the Kimura 2-parameter model . To evaluate the amount of population genetic structure pairwise FST values were calculated using DnaSP ver5.
Phylogenetic analysis and divergence time estimates
To infer phylogenetic relationships, the haplotypes data were used to construct phylogenetic trees using the ML, BI and the MJN methods. The best substitution model to run the ML tree and BI was Kimura 2-parameter model with discrete gamma distribution (K2P + G) as determine using MEGA ver6.l. To assess the robustness of the ML tree, bootstrapping  with 2000 replicates were conducted. In BEAST package ver2.5 , BI was analyzed using two independent runs with 10 million of MCMC chain length each using strict clock settings, Yule model prior , and sub-sampled at every 1000 generations. To set the substitution model to K2P + G, HKY85 model  was selected with all the base frequencies set to equal and all rate parameters were fixed to 1.0. The convergence of all parameters was assessed using TRACER ver1.6  and both independent runs were then combined using the software LogCombiner ver2.5 available within the BEAST ver2.0 package. A consensus tree was later created from the combined tree files after a burn-in of 10% using the TreeAnnotator ver2.5 also available within the BEAST ver2.0 package. Finally, the MJN was constructed using the Network ver4.6 .
To estimate the divergence dates of several important events in the evolutionary and dispersal history of M. fascicularis in SEA, the BI tree was recalibrated using the same BI analysis parameters as described above. Using the proposed divergence date of 5.5 mya for the LCA of Macaca [9, 31, 84], the priors were set to normal distribution with the mean set to 5.5 and sigma at 0.5.
Availability of data and materials
Sequence data produced in this study are available in the NCBI (https://www.ncbi.nlm.nih.gov/) with the Accession No.: KC572634-KC572687; KJ690361-KJ690376; KJ733018-KJ733278.
Department of Wildlife and National Parks
Thousand years ago
Last common ancestor
Million years ago
National Center for Biotechnology Information
Fooden J, Albrecht GH. Latitudinal and insular variation of skull size in crab-eating macaques. Am J Phys Anthropol. 1993;92:521–38.
Fooden J. Systematic review of Southeast Asian long-tail macaques, Macaca fascicularis (Raffles, 1821). Fieldiana Zool. 1995;81:1–206.
Eudey AA. The crab-eating macaque (Macaca fascicularis): Widespread and rapidly declining. Primate Conserv. 2008;23:129–32.
Osman NAW. A study on composition and density of long-tailed macaque (Macaca fascicularis) in Kuala Lumpur, 1998. JWP. 1998;16:80–4.
DWNP. Pelan pengurusan kera (Macaca fascicularis) bermasalah di Semenanjung Malaysia. Department of Wildlife and Parks. 2006. https://www.wildlife.gov.my/images/stories/penerbitan/lain_lain/PelanPengurusanKeraBermasalah.pdf. Accessed 11 Oct 2016.
Payne J, Francis CM. A field guide to the mammals of Borneo. Kota Kinabalu: Sabah Society; 2005.
Pusey AE, Packer C. Dispersal and philopatry. In: Smuts BB, Cheney DL, Seyfarth RM, Wrangham RW, Struhsaker TT, editors. Primate societies. Chicago: University of Chicago Press; 1987. p. 250–66.
Melnick DJ, Hoelzer GA. Differences in male and female macaque dispersal lead to contrasting distributions of nuclear and mitochondrial DNA variation. Int J Primatol. 1992;13:379–93.
Tosi AJ, Morales JC, Melnick DJ. Paternal, maternal, and biparental molecular markers provide unique windows onto the evolutionary history of macaque monkeys. Evolution. 2003;57:1419–35.
Avise JC. Phylogeography: retrospect and prospect. J Biogeogr. 2009;36:3–15.
Jobling MA, Tyler-Smith C. Fathers and sons: the Y chromosome and human evolution. Trends Genet. 1995;11:449–56.
Keane B, Dittus WPJ, Melnick DJ. Paternity assessment in wild groups of toque macaques Macaca sinica at Polonnaruwa, Sri Lanka using molecular markers. Mol Ecol. 1997;6:267–82.
Disotell TR. Human evolution: sex-specific contributions to genome variation. Curr Biol. 1999;9:R29-31.
Chan YC, Roos C, Inoue-Murayama M, Inoue E, Shih CC, Vigilant L. A comparative analysis of Y-chromosome and mtDNA phylogenies of the Hylobates gibbons. BMC Evol Biol. 2012;12:150.
Tosi AJ, Morales JC, Melnick DJ. Y-chromosome and mitochondrial markers in Macaca fascicularis indicate introgression with Indochinese M. mulatta and a biogeographic barrier in the Isthmus of Kra. Int J Primatol. 2002;23:161–77.
Tosi AJ, Morales JC, Melnick DJ. Comparison of Y-chromosome and mtDNA phylogenies leads to unique inferences of macaque evolutionary history. Mol Phylogenet Evol. 2000;17:133–44.
Tosi AJ, Buzzard PJ, Morales JC, Melnick DJ. Y-chromosome data and tribal affiliations of Allenopithecus and Miopithecus. Int J Primatol. 2002;23:1287–99.
Tosi AJ, Disotell TR, Morales JC, Melnick DJ. Cercopithecine Y-chromosome data provide a test of competing morphological evolutionary hypotheses. Mol Phylogenet Evol. 2003;27:510–21.
Tosi AJ, Melnick DJ, Disotell TR. Sex chromosome phylogenetics indicate a single transition to terrestriality in the guenons (tribe Cercopithecini). J Hum Evol. 2004;46:223–37.
Tosi AJ, Detwiler KM, Disotell TR. Y-chromosomal markers suitable for noninvasive studies of guenon hybridization. Int J Primatol. 2005;26:685–96.
Rovie-Ryan JJ, Abdullah MT, Sitam FT, Abidin ZZ, Tan SG. Y-chromosomal gene flow of Macaca fascicularis (Cercopithecidae) between the insular and mainland peninsula of Penang state, Malaysia. JoSTT. 2013;9:113–26.
Bunlungsup S, Imai H, Hamada Y, Gumert MD, San AM, Malaivijitnond S. Morphological characteristics and genetic diversity of Burmese long-tailed Macaques (Macaca fascicularis aurea). Am J Primatol. 2016;78:441–55.
Bunlungsup S, Imai H, Hamada Y, Matsudaira K, Malaivijitnond S. Mitochondrial DNA and two Y-chromosome genes of common long-tailed macaques (Macaca fascicularis fascicularis) throughout Thailand and vicinity. Am J Primatol. 2017;79:e22596.
Abegg C, Thierry B. Macaque evolution and dispersal in insular south-east Asia. Biol J Linn. 2002;75:555–76.
Bird MI, Taylor D, Hunt C. Palaeoenvironments of insular Southeast Asia during the Last Glacial Period: a savanna corridor in Sundaland? Quat Sci Rev. 2005;24:2228–42.
Harrison T, Krigbaum J, Manser J. Primate biogeography and ecology on the Sunda shelf Islands: a paleontological and zooarchaeological perspective. In: Lehman SM, Fleagle JG, editors. Primate biogeography: progress and prospects. New York: Springer Science; 2006. p. 331–72.
Cannon CH, Morley RJ, Bush ABG. The current refugial rainforests of Sundaland are unrepresentative of their biogeographic past and highly vulnerable to disturbance. PNAS. 2009;106:11188–93.
Arora N, Nater A, van Shaik CP, Willems EP, van Noordwijk MA, Goossens B, Morf N, Bastian M, Knott C, Morrogh-Bernard H, Kuze N, Kanamori T, Pamungkas J, Perwitasari- Farajallah D, Verschoor E, Warren K, Krutzen M. Effects of Pleistocene glaciations and rivers on the populations structure of Bornean orangutans (Pongo pygmaeus). PNAS. 2010;107:21376–81.
Wurster CM, Bird MI, Bull ID, Creed F, Bryant C, Dungait JAJ, Paz V. Forest contraction in north equatorial Southeast Asia during the Last Glacial Period. PNAS. 2010;107:15508–11.
Karuppannan KV, Saaban S, Firdaus Ariff AR, Mustapa AR. Non-surgical castration in controlling long tailed macaques (Macaca fascicularis) population by Department of Wildlife and National Parks (DWNP) Peninsular Malaysia. MJVR. 2013;4:33–6.
Delson E. Fossil macaques, phyletic relationships and a scenario of deployment. In: Lindburg DG, editor. The macaques: studies in ecology, behavior, and evolution. New York: Van Nostrand Rheinhold Co; 1980. p. 10–30.
Deinard A, Smith DG. Phylogenetic relationships among the macaques: Evidence from the nuclear locus NRAMP1. J Hum Evol. 2001;41:45–59.
Li J, Han K, Xing J, Kim HS, Rogers J, Ryder OA, Disotell T, Yue B, Batzer MA. Phylogeny of macaques (Cercopithecidae: Macaca) based on Alu elements. Gene. 2009;448:242–9.
Fan Z, Zhou A, Osada A, Yu J, Jiang J, Li P, Du L, Niu L, Deng J, Xu H, Xing J, Yue B, Li J. Ancient hybridization and admixture in macaques (genus Macaca) inferred from whole genome sequences. Mol Phylogenet Evol. 2018;127:376–86.
Gibbard PL, Head MJ, Walker MJ. Formal ratification of the quaternary system/period and the Pleistocene Series/Epoch with a base at 2.58 Ma. J Quat Sci. 2010;25:96–102.
Gibbard PL, Head MJ. The newly-ratified definition of the Quaternary System/Period and redefinition of the Pleistocene Series/Epoch, and comparison of proposals advanced prior to formal ratification. Episodes. 2010;33:152–8.
Cranbrook E, Piper PJ. Paleontology to policy: the Quaternary history of Southeast Asian tapirs (Tapiridae) in relation to large mammal species turnover, with a proposal for conservation of Malayan tapir by reintroduction to Borneo. Integr Zool. 2013;8:95–120.
Osada N, Uno Y, Mineta K, Kameoka Y, Takahashi I, Terao K. Ancient genome-wide admixture extends beyond the current hybrid zone between Macacafascicularis and M. mulatta. Mol Ecol. 2010;19:2884–95.
Bunlungsup S, Kanthaswamy S, Oldt RF, Smith DG, Houghton P, Hamada Y, Malaivijitnond S. Genetic analysis of samples from wild populations opens new perspectives on hybridization between long-tailed (Macaca fascicularis) and rhesus macaques (Macaca mulatta). Am J Primatol. 2017;79:e22726.
Yan G, Zhang G, Fang X, Zhang Y, Li C, Ling F, Cooper DN, Li Q, Li Y, van Gool AJ, Du H, Chen J, Chen R, Zhang P, Huang Z, Thompson JR, Meng Y, Bai Y, Wang J, Zhuo M, Wang T, Huang Y, Wei L, Li J, Wang Z, Hu H, Yang P, Le L, Stenson PD, Li B, Liu X, Ball EV, An N, Huang Q, Zhang Y, Fan W, Zhang X, Li Y, Wang W, Katze MG, Su B, Nielsen R, Yang H, Wang J, Wang X, Wang J. Genome sequencing and comparison of two nonhuman primate animal models, the cynomolgus and Chines rhesus macaques. Nat Biotechnol. 2011;29:1019–23.
Malaivijitnond S, Arsaithamkul V, Tanaka H, Pomchote P, Jaroenporn S, Suryobroto B, Hamada Y. Boundary zone between northern and southern pig-tailed macaques and their morphological differences. Primates. 2012;53:377–89.
Wallace AR. The geographical distribution of animals. London: Macmillan; 1876.
Woodruff DS. Neogene marine transgressions, palaeogeography and biogeographic transitions on the Thai-Malay Peninsula. J Biogeogr. 2003;30:551–67.
Doxiadis GGM, de Groot N, de Groot NG, Rotmans G, de Vos-Rouweler AJM, Bontrop RE. Extensive DRB region diversity in cynomolgus macaques: recombination as a driving force. Immunogenetics. 2010;62:137–47.
Inger RF, Voris HK. The biogeographical relations of the frogs and snakes of Sundaland. J Biogeogr. 2001;28:863–91.
Abdullah MT. Biogeography and variation of Cynopterus brachyotis in Southeast Asia. D. Phil. Thesis. Queensland: University of Queensland; 2003.
Hughes JB, Round PD, Woodruff DS. The Indochinese-Sundaic faunal transition at the Isthmus of Kra: an analysis of resident forest bird species distributions. J Biogeogr. 2003;30:569–80.
Round PD, Hughes JB, Woodruff DS. Latitudinal range limits of resident forest birds in Thailand and the Indochinese-Sundaic zoogeographic transition. Nat Hist Bull Siam Soc. 2003;51:69–96.
Luo SJ, Kim JH, Johnson WE, van der Walt J, Martenson J, Yuhki N, Miquelle DG, Uphyrkina O, Goodrich JM, Quigley HB, Tilson R, Brady G, Martelli P, Subramaniam V, McDougal C, Hean S, Huang SQ, Pan W, Karanth UK, Sunquist M, Smith JLD, O’Brien SJ. Phylogeography and genetic ancestry of tigers (Panthera tigris). PLoS Biol. 2004;2:2275–93.
Woodruff DS, Turner LM. The Indochinese-Sundaic zoogeographic transition: a description and analysis of terrestrial mammal species distributions. J Biogeogr. 2009;36:803–21.
Smith DG, Mcdonough JW, George DA. Mitochondrial DNA variation within and among regional populations of longtail macaques (Macaca fascicularis) in relation to other species of the fascicularis group of macaques. Am J Primatol. 2007;68:182–98.
Blancher A, Bonhomme M, Crouau-Roy B, Terao K, Kitano T, Saitou N. Mitochondrial DNA sequence phylogeny of 4 populations of the widely distributed cynomolgus macaque (Macaca fascicularis fascicularis). J Hered. 2008;99:254–64.
Shiina T, Tanaka K, Katsuyama Y, Otabe K, Sakamoto K, Kurata M, Nomura M, Yamanaka H, Nakagawa H, Inoko H, Ota M. Mitochondrial DNA diversity among three subpopulations of cynomolgus macaques (Macaca fascicularis) originating from the Indochinese region. Exp Anim. 2010;59:567–78.
Abdul-Latiff MAB, Ruslin F, Faiq H, Hairul MS, Rovie-Ryan JJ, Abdul-Patah P, Yaakop S, Md-Zain BM. Continental monophyly and molecular divergence of peninsular Malaysia's Macaca fascicularis fascicularis. BioMed Res Int. 2014;897682.
Abdul-Latiff MAB, Ruslin F, Fui VV, Abu MH, Rovie-Ryan JJ, Abdul-Patah P, Lakim M, Roos C, Yaakop S, Md-Zain BM. Phylogenetic relationships of Malaysia’s long-tailed macaques, Macaca fascicularis, based on cytochrome b sequences. ZooKeys. 2014;407:121–39.
Liedigk R, Kolleck J, Böker KO, Meijaard E, Md-Zain BM, Abdul-Latiff MAB, Ampeng A, Lakim M, Abdul-Patah P, Tosi AJ, Brameier M, Zinner D, Roos C. Mitogenomic phylogeny of the common long-tailed macaques (Macaca fascicularis fascicularis). BMC Genomics. 2015;16:222.
Fa JE. The genus Macaca: a review of taxonomy and evolution. Mammal Rev. 1989;19:45–81.
Fooden J. Comparative review of fascicularis-group species of macaques (Primates: Macaca). Fieldiana Zool. 2006;107:1–43.
Voris HK. Maps of Pleistocene sea levels in Southeast Asia: shorelines, river systems and time durations. J Biogeogr. 2000;27:1153–67.
Meijaard E, Groves CP. Biogeography and phylogeography of Presbytis. Primate Rep. 2004;68:71–90.
Stimpson CM. Local scale, proxy evidence for the presence of closed canopy forest in North-western Borneo in the Late Pleistocene: bones of Strategy I bats from the archaeological record of the Great Cave of Niah. Sarawak Palaeogeogr Palaeoclimatol Palaeoecol. 2012;331–332:136–49.
Tosi AJ, Coke CS. Comparative phylogenetics offer new insights into the biogeographic history of Macaca fascicularis and the origin of the Mauritian macaques. Mol Phylogenet Evol. 2007;42:498–504.
Gathorne-Hardy FJ, Syaukani, Davies RG, Eggleton P, Jones DT. Quaternary rainforest refugia in southeast Asia: using termites (Isoptera) as indicators. Biol J Linn Soc. 2002;75:453–66.
O’Brien SJ, Mayr E. Bureaucratic mischief: recognizing endangered species and subspecies. Science. 1991;251:1187–8.
Moritz C. Defining “evolutionary significant units” for conservation. Tree. 1994;9:373–5.
Kanthaswamy S, Kurushima JD, Smith DG. Inferring Pongo conservation units: a perspective based on microsatellite and mitochondrial DNA analyses. Primates. 2006;47:310–21.
Kanthaswamy S, Satkoski J, George D, Kou A, Erickson BJA, Smith DG. Hybridization and stratification of nuclear genetic variation in Macacamulatta and M. fascicularis. Int J Primatol. 2008;29:1295–311.
Stevison LS, Kohn MH. Determining genetic background in captive stocks of cynomolgus macaques (Macaca fascicularis). J Med Primatol. 2008;37:311–7.
Affara N, Bishop C, Brown W, Cooke H, Davey P, Ellis N, Graves JM, Jones M, Mitchell M, Rappold G, Tyler-Smith C, Yen P, Lau YFC. Report of the second international workshop on Y chromosome mapping 1995. Cytogenet Cell Genet. 1996;73:33–76.
Kim HS, Hirai H, Takenaka O. Molecular features of the TSPY gene of gibbons and Old World monkeys. Chromosome Res. 1996;4:500–6.
Glaser B, Grutzner F, Willman U, Stanyon R, Arnold N, Taylor K, Rietschel W, Zeitler S, Toder R, Schempp W. Simian Y chromosomes: Species-specific rearrangements of DAZ, RBM, and TSPY versus contiguity of PAR and SRY. Mamm Genome. 1998;9:226–31.
Drummond AJ, Ashton B, Buxton S, Cheung M, Cooper A, Duran C, Field M, Heled J, Kearse M, Markowitz S, Moir R, Stones-Havas S, Sturrock S, Thierer T, Wilson A. Geneious (Version 5.6.). 2012. http://www.geneious.com.
Swofford DL. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4. Sunderland: Sinauer Associates; 2003.
Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA: molecular evolutionary genetics analysis (Version 6.0). Mol Biol Evol. 2013;30:2725–9.
Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.
Nei M. Molecular evolutionary genetics. New York: Columbia University Press; 1987.
Kimura M. A simple method for estimating evolutionary rate of base substitutions through comparative studies of nucleotide sequences. J Mol Evol. 1980;16:111–20.
Felsenstein J. Confidence limits on phylogenies: an approach using the bootstrap. Evolution. 1985;39:783–91.
Bouckaert R, Heled J, Kühnert D, Vaughan TG, Wu CH, Xie D, Suchard MA, Rambaut A, Drummond AJ. BEAST2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10:e1003537.
Heled J, Drummond AJ. Calibrated tree priors for relaxed phylogenetics and divergence time estimation. Syst Biol. 2012;61:138–49.
Hasegawa M, Kishino H, Yano T. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985;22:160–74.
Rambaut A, Suchard MA, Xie D, Drummond AJ. Tracer (Version 1.6). 2014. http://beast.bio.ed.ac.uk/Tracer.
Bandelt HJ, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.
Delson E. The oldest monkeys in Asia. In: Takenaka O, editor. Evolution of Asian primates. Inuyama, Aichi, Japan: International symposium conducted at Freude & Kyoto University, Primate Research Institute; 1996. p. 40.
We acknowledge the following local authorities for the permission to collect samples: Department of Wildlife and National Parks (DWNP), Forest Department of Sarawak, and Sabah Wildlife Department. We are thankful to the Outbreak Response Team of DWNP, UNIMAS Primate Team, and EcoHealth Alliance for assistance during sampling. We also thank the UNIMAS PETARY Library for Publication Support Fee.
This project was funded by the DWNP for laboratory equipment and consumables, and partly by the Proboscis Genome Research grant for travel and accommodation, awarded to M. T. Abdullah and colleagues. Funding for sampling work for ZMW coded samples was supported in part by the USAID Emerging Pandemic Threat Program—PREDICT Project and by the Skoll Foundation and Google Incorporated through Global Viral. All funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript. Open access funding provided by Universiti Malaysia Sarawak.
Ethics approval and consent to participate
Sampling at Peninsular Malaysia was permitted by the Research Committee of DWNP, at Sarawak by Forest Department of Sarawak (Permit No: 37/2011), and at Sabah by the Sabah Biodiversity Council (Permit No: JKM/MBS.1000-2/2(103)). The sampling protocol was approved by UC DAVIS Institutional Animal Care and Use Committee (Protocol number: 16048).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Sample identification, localities, haplotypes designation, and the GenBank Accession Numbers for each of the samples used in this study.
: Table S2. Species, sample identification, localities, haplotypes designation, and the GenBank Accession Numbers for downloaded sequences from GenBank used in this study.
About this article
Cite this article
Rovie-Ryan, J.J., Khan, F.A.A. & Abdullah, M.T. Evolutionary pattern of Macaca fascicularis in Southeast Asia inferred using Y-chromosomal gene. BMC Ecol Evo 21, 26 (2021). https://doi.org/10.1186/s12862-021-01757-1
- Continental and insular sundaic
- Secondary contact