S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E GENETICS Copyright © 2023 The Authors, some Comparative genomics reveals the hybrid origin of a rights reserved; exclusive licensee macaque group American Association for the Advancement of Science. No claim to Bao-Lin Zhang1†, Wu Chen2†, Zefu Wang3,4†, Wei Pang5, Meng-Ting Luo5, Sheng Wang1, original U.S. Government Yong Shao1, Wen-Qiang He5, Yuan Deng6,7, Long Zhou8, Jiawei Chen6, Min-Min Yang1, Works. Distributed Yajiang Wu2, Lu Wang9, Hugo Fernández-Bellon10, Sandra Molloy11, Hélène Meunier12,13, under a Creative Fanélie Wanert14, Lukas Kuderna15, Tomas Marques-Bonet16,17,18,19, Christian Roos20,21, Commons Attribution Xiao-Guang Qi9, Ming Li22, Zhijin Liu23, Mikkel Heide Schierup24, David N. Cooper25, NonCommercial License 4.0 (CC BY-NC). Jianquan Liu3,26, Yong-Tang Zheng5,27*, Guojie Zhang1,8,28,30*, Dong-Dong Wu1,27,29,31* Although species can arise through hybridization, compelling evidence for hybrid speciation has been reported only rarely in animals. Here, we present phylogenomic analyses on genomes from 12 macaque species and show Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 that the fascicularis group originated from an ancient hybridization between the sinica and silenus groups ~3.45 to 3.56 million years ago. The X chromosomes and low-recombination regions exhibited equal contributions from each parental lineage, suggesting that they were less affected by subsequent backcrossing and hence could have played an important role in maintaining hybrid integrity. We identified many reproduction-associ- ated genes that could have contributed to the development of the mixed sexual phenotypes characteristic of the fascicularis group. The phylogeny within the silenus group was also resolved, and functional experimenta- tion confirmed that all extant Western silenus species are susceptible to HIV-1 infection. Our study provides novel insights into macaque evolution and reveals a hybrid speciation event that has occurred only very rarely in primates. INTRODUCTION framework to unequivocally demonstrate hybrid speciation has Interspecific hybridization can facilitate species adaptation by intro- yet to be established. ducing new genetic material and novel allelic combinations (1). In The genus of macaques (Macaca) represents an excellent model some circumstances, hybridization may result in the near-instanta- in which to study the interplay between interspecific hybridization neous formation of new species, i.e., hybrid speciation, thereby pro- and speciation. As one of the most successful primate lineages, ma- moting biodiversity (2). Hybrid speciation has been an important caques now include 23 species that are widely distributed across mode of speciation during plant evolution (3). However, hybrid South, East, and Southeast Asia, the only exception being M. sylva- speciation has been considered to occur only rarely in animals nus, which is confined to the Atlas mountains in North Africa (15). owing to the fact that the derived species from interspecific hybrid- Paleontological and molecular data suggest that macaques originat- ization typically tend to be evolutionarily less fit (e.g., Dobzhansky- ed in North Africa ~7 million years (Ma) ago and then experienced Muller incompatibility) or are only weakly reproductively isolated a burst of speciation in Asia during the past 5 Ma (16, 17). Some from their parental species and hence could be swamped by one macaque species are widely used as laboratory models for studying of their parents (4). human disease and for vaccine development (18). Although Both theoretical and empirical studies suggest that those species members of this genus are morphologically and behaviorally dis- groups that have experienced rapid adaptive radiation are prone to tinct, interspecific hybridization is potentially possible between hybridization because the incompatibilities between these species any geographically overlapping pair of species because their repro- are weak (1). Hybrid lineages are more readily established when ductive isolation is likely to be incomplete. This notion is supported hybrids mate assortatively with each other (5) and/or when they by both field observations (19, 20) and molecular studies (21, 22). became adapted to a new environment (6). An increasing number While some macaque species (e.g., the fascicularis group of ma- of empirical studies have documented the establishment of hybrid caques and the stump-tailed macaque) exhibit distinctive mixed lineages in rapidly radiating groups, including butterflies, canids, phenotypes, which could, in principle, have arisen by genomic ad- baboon, birds, and bears (7–11), suggesting that hybrid speciation mixture (23, 24), the broader role of hybridization and its specific may be more common and hence more important than previously outcomes remain largely unknown. Here, we perform multiple thought. Nevertheless, concrete evidence of hybrid speciation in genome analyses to ascertain species phylogeny in the macaque animals remains scarce and putative examples are often highly con- genus, assess the role of interspecific hybridization in speciation, tentious (12, 13), impeding our understanding of the genomic and explore the genetic basis of mixed phenotypes in mechanisms underlying hybrid speciation. In particular, it is macaque species. often challenging to reliably differentiate between hybrid speciation and post-speciation genetic introgression (4), which may be further complicated by the stochastic sorting of ancestral polymorphisms RESULTS among descendant lineages, a process known as incomplete We generated 10 high-quality macaque genome assemblies (tables lineage sorting (ILS) (14). Thus, a comparative phylogenomic S1 to S3) using the long-read sequencing strategy of Nanopore and Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 1 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E long fragment read technologies (stLFR and 10X Genomics), re- One supported a sister relationship between the fascicularis and spectively. Scaffold N50 sizes of these newly assembled genomes sinica groups (T1), whereas the other (T2) favored the fascicularis ranged from 17.7 to 33.2 Mb, while the Benchmarking Universal group as a sister lineage to the silenus group. These two tree topol- Single-Copy Orthologs (BUSCO) completeness scores ranged ogies were supported by a total of 78% of windows interleavedly dis- from ~93.3 to 94.5%. The final annotated gene numbers ranged tributed across the genome, with 47 and 31% referring to T1 and T2, from 20,662 to 21,811 between different species using a combina- respectively (Fig. 2, A and B). Specifically, we found that the X chro- tion of ab initio and homology-based gene prediction approaches. mosome, where genomic incompatibilities normally first develop We also included two previously published genomes of M. mulatta during speciation (29), exhibited a nearly equal proportion of T1 (25) and M. nemestrina, and one outgroup species Papio hamadryas and T2 topologies versus the autosomes irrespective of window (26). Thus, our genome assemblies cover all known macaque species size (fig. S10). Such widespread mixed ancestry of the autosomes groups (Fig. 1A) (17). and the X chromosome led us to speculate that the fascicularis We first performed pairwise whole-genome alignments against group may have originated from an ancient hybridization the Chinese rhesus macaque (M. mulatta, rheMacS) genome using between the progenitors of the sinica group and those of the the LASTZ program (27) and then merged all of them into multiple silenus group. genome alignments with MULTIZ (28). These genomes showed If the hybrid origination hypothesis is correct, we would expect Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 high collinearity, with the aligned base pairs spanning more than to observe an equal level of sequence divergence between the hybrid 2.7 Gb or 90% of the reference assembly (fig. S1 and table S5). To species and its two parental lineages. To test this postulate, we com- deduce the correct branching pattern of the sequenced macaque pared the sequence divergence (DXY) and relative divergence time species and to characterize the heterogeneous phylogenetic signals between fascicularis versus sinica and fascicularis versus silenus for across genomes, we partitioned the genome alignments into non- the autosomes and X chromosome separately to allow for their dif- overlapping windows of 50 kb and then performed maximum like- ferent rates of evolution (30). We found no obvious difference in lihood (ML) analyses for each window sequence. We initially these two statistics calculated using only the sex chromosome data obtained the multispecies coalescent-based species tree (ASTRAL (both P values > 0.05, Wilcoxon’s test; Fig. 2, C and D). With the and STAR; fig. S3) from windows with high bootstrap values (e.g., autosomes, although both parameters were statistically significant mean value of >80%). Although the topologies obtained from these (P < 0.001, Wilcoxon’s test; Fig. 2, C and D), the mean DXY value two methods were mutually consistent and largely in accord with (0.1616 versus 0.1631) and the mean estimated age (3.45 Ma versus recent studies (17, 21), we observed substantial genealogical dis- 3.56 Ma) between the fascicularis group and its two progenitors cordance within and between species groups among all window were only slightly different. The observed pattern therefore trees (n = 7392 alternative topologies; Fig. 1, B and C, and fig. concurs with that expected under the hybrid speciation scenario. S5), mostly pertaining to the phylogenetic position of the fascicula- However, extensive post-speciation gene flow may also have pro- ris group with respect to the sinica and silenus groups. The discord- duced a similar signature (31). To test this possibility, we performed ance of this position was not a trivial consequence of the choice of D-statistic analyses in a sliding window (32). If either T1 or T2 rep- window size because analyses using either smaller (20 kb) or larger resent the true species phylogeny, the windows least affected by gene (100 kb) window sizes yielded broadly similar results (figs. S6 to S9). flow (absolute D-statistic values close to zero) would overwhelming- To better understand the interrelationship of macaque groups with ly support one of the tested trees, whereas hybridization would incongruent phylogenies, we pruned the trees so as to include five support both (33). Our results supported the latter prediction ingroup species, each representing a particular lineage proposed by showing that the 1% of windows with the lowest absolute D-statistic previous studies (16, 17). Two topologies were found to dominate: values supported T1 and T2, respectively, when using T1 and T2 as 1 State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming 650223, China. 2Guangzhou Zoo and Guangzhou Wildlife Research Center, Guangzhou 510070, China. 3Key Laboratory for Bio-resource and Eco-environment of Ministry of Education, College of Life Sciences, Sichuan University, Chengdu 610065, China. 4Co-Innovation Center for Sustainable Forestry in Southern China, College of Biology and the Environment, Nanjing Forestry University, Nanjing 210037, China. 5Key Laboratory of Animal Models and Human Disease Mechanisms of the Chinese Academy of Sciences, KIZ-CUHK Joint Laboratory of Bioresources and Molecular Research in Common Diseases, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming 650223, China. 6BGI-Shenzhen, Shen- zhen 518083, China. 7Section for Ecology and Evolution, Department of Biology, University of Copenhagen, Copenhagen DK-2100, Denmark. 8Center for Evolutionary and Organismal Biology and Women’s Hospital at Zhejiang University School of Medicine, Hangzhou 310058, China. 9Shaanxi Key Laboratory for Animal Conservation, College of Life Sciences, Northwest University, Xi’an, China. 10Barcelona Zoo, Parc de La Ciutadella, Barcelona 08003, Spain. 11Dublin Zoo, Dublin 8, Ireland. 12Centre de Primatologie, de l’Université de Strasbourg, Niederhausbergen, France. 13Laboratoire de Neurosciences Cognitives et Adaptatives, UMR 7364, Université de Stras- bourg, Strasbourg, France. 14Plateforme SILABE, Université de Strasbourg, Niederhausbergen, France. 15Genome Interpretation Department, Illumina Inc., Foster City, CA, USA. 16Institute of Evolutionary Biology (UPF-CSIC), PRBB, Dr. Aiguader 88, Barcelona 08003, Spain. 17Catalan Institution of Research and Advanced Studies (ICREA), Passeig de Lluís Companys, 23, Barcelona 08010, Spain. 18CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Baldiri i Reixac 4, Barcelona 08028, Spain. 19Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, c/Columnes s/n, 08193 Cer- danyola del Vallès, Barcelona, Spain. 20Primate Genetics Laboratory, German Primate Center, Göttingen, Germany. 21Gene Bank of Primates, German Primate Center, Göttingen, Germany. 22CAS Key Laboratory of Animal Ecology and Conservation Biology, Institute of Zoology, Chinese Academy of Sciences, Beijing 100101, China. 23 College of Life Sciences, Capital Normal University, Beijing 100048, China. 24Bioinformatics Research Centre, Aarhus University, Aarhus C DK-8000, Denmark. 25 Institute of Medical Genetics, School of Medicine, Cardiff University, Cardiff CF14 4XN, UK. 26State Key Laboratory of Grassland Agro-ecosystem, Institute of Innovation Ecology and College of Life Sciences, Lanzhou University, Lanzhou 730000, China. 27National Resource Center for Non-Human Primates, Kunming Primate Research Center and National Research Facility for Phenotypic and Genetic Analysis of Model Animals (Primate Facility), Kunming Institute of Zoology, Chinese Academy of Sci- ences, Kunming, Yunnan 650107, China. 28Liangzhu Laboratory, Zhejiang University Medical Center, 1369 West Wenyi Road, Hangzhou 311121, China. 29Center for Ex- cellence in Animal Evolution and Genetics, Chinese Academy of Sciences, Kunming 650223, China. 30Villum Center for Biodiversity Genomics, Section for Ecology and Evolution, Department of Biology, University of Copenhagen, Copenhagen 2100, Denmark. 31Kunming Natural History Museum of Zoology, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming, Yunnan 650223, China. *Corresponding author. Email: wudongdong@mail.kiz.ac.cn (D.-D.W.); guojiezhang@zju.edu.cn (G.Z.); zhengyt@mail.kiz.ac.cn (Y.-T.Z.) †These authors contributed equally to this work. Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 2 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 Fig. 1. Distribution map and the discordance of phylogeny. (A) Distribution map of macaque species used in this study. Nomenclature follows Delson (16) in defining the species group. We further split the silenus group into Western (W) and Eastern (E) forms on the basis of their similar divergence to other species groups. Macaque drawings are copyright, 2013, Stephen D. Nash, International Union for Conservation of Nature Species Survival Commission Primate Specialist Group and are used with permission. All species except M. mulatta and M. nemestrina were newly sequenced in this study. (B) DensiTree plot for 50-kb window trees. Blue, red, and green colors represent the first, second, and third most common topologies, respectively, whereas gray represents other topologies. (C) Majority-rule consensus tree of all 50-kb window trees (n = 46,003). The numbers above the branches indicate the absolute number of topologies supporting the splits, whereas the numbers below represent the percentage values. the tested topologies (Fig. 2, E and F). Furthermore, on the basis that than post-speciation gene flow, was responsible for the pattern of regions of lower recombination are “cold spots” for introgression genomic mosaicism. (29), we classified the trees by recombination rate and found that Another underlying assumption of hybrid speciation is that all T1 and T2 were still dominant in the low recombination rate the individuals from the hybrid species are uniformly admixed. regions of 0 to 0.1 cM/Mb (fig. S11). This observation was even However, if ancient directional gene flow and ILS were to have oc- more evident with a smaller window size (20 kb), where the ratio curred, it is likely that not all individuals in the hybrid group would of T1 and T2 was almost identical (T1:T2 = 32%:31%; fig. S13), have contained the same amount and length distribution of intro- probably due to the averaging effect with the larger window size. gressed alleles because subsequent recombination and negative se- These findings corroborate the view that hybrid speciation, rather lection would have acted so as to purge the deleterious intruder Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 3 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 Fig. 2. Phylogenetic relationship of the five major macaque lineages across the genome. (A) Distribution of the three most common topologies by reference to the Chinese rhesus macaque (M. mulatta) genome. Colored bands represent tree topologies of each 50-kb window. White interval regions denote missing data. (B) Three most common trees (T1 to T3) recovered by ML analysis and their frequency in consecutive 50-kb windows. Values below the tree refer to the percentage of windows recovering that topology. The outgroup P. hamadryas is not shown. (C) Genetic divergence (DXY) and (D) estimated divergence times from MCMCTREE between the fascicularis group and its two putative parental lineages (sinica and silenus) based on autosomal 50-kb window sequences (left) and the X chromosome (right), respec- tively. P values were estimated by the Wilcoxon rank sum test. The substitution rates estimated from MCMCTREE for the autosomes and X chromosome were 1.2 × 10−9 per site per year and 0.9 × 10−9 per site per year, respectively. (E) DensiTree plot of ML trees derived from the 1% windows with the lowest absolute D-statistic values using T1, and (F) T2, as the tested relationship. The values below show the percentage of windows that recovered the topology. alleles from the acceptor genome during a long period of divergence fuscata) or widely distributed (M. mulatta and M. fascicularis). (34). To test this prediction, we used HyDe (35) and LOTER (36) to This result was further corroborated by PhyloNet-MPL analyses, quantify the genomic contributions of the two parental lineages an alternative method that allows for ILS and hybridization simul- (sinica and silenus) to each of the hybrid species in the fascicularis taneously, based on the maximum pseudo-likelihood method (37), group. The results were consistent with our expectation under the where the section of the fascicularis group was invariably identified hybridization hypothesis, showing that all three species from the as a reticulate node in the scenarios allowing one and two past hy- fascicularis group were uniformly mixed and exhibited a general ex- bridization events (Fig. 4A and fig. S17). ponential decay of the tract size (Fig. 3, A and B, and table S7), ir- Together, these phylogenomic analyses concur in terms of pro- respective of whether the species are geographically isolated (M. viding consistent support for the hybrid origin hypothesis of the Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 4 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 Fig. 3. Genomic ancestry inherent in two parental lineages. (A) Density plot of estimated γ values across 500 bootstrap replicates in HyDe analyses. γ represents the estimated probability of inheritance from the ancestor of the sinica group, whereas 1-γ represents the probability of inheritance from the silenus group. The mean γ values of the three hybrid species are given top right. Mmul, M. mulatta; Mfas, M. fascicularis; Mfus, M. fuscata. (B) Distribution of tract sizes (in 5-kb bins) of sinica group ancestry in three hybrid species. fascicularis group but were not informative regarding its direction (M. nemestrina) are sister species (39). Although these two in the initial hybridization. Because the mitochondrial genome and species are close ecologically and geographically and are also the Y chromosome are representative of the maternal and paternal similar in phenotypic appearance (Fig. 1A), all our phylogenomic lineages, respectively, we next performed phylogenetic analyses analyses with different genomic data types support the postulate based on mitochondrial and Y chromosome sequences (~480 kb that the lion-tailed macaque (M. silenus) from the South Indian in length) to determine the hybridization pattern. We found a Western Ghats and M. leonina are sister species (Fig. 4, A to C, clear topology discrepancy between the mitochondrial and the Y and figs. S14 and S15). Our previous understanding of the phylog- chromosomal data, with the fascicularis group having a mitochon- eny of the silenus group was also confused by the notion that M. drial genome similar to that of the sinica group and a Y chromo- leonina and M. nemestrina were the only two Old World primates some similar to that of the silenus group (Fig. 4, B and C), known to be susceptible to HIV-1 (human immunodeficiency virus suggesting that the ancient hybridization occurred predominantly type 1) infection due to the highly unusual retrotranspositional in- or exclusively between proto-sinica group females and proto- sertion of a cyclophilin A2 (CypA2) gene into the 3′ untranslated silenus group males. region of the TRIM5 locus (40, 41). On the basis of the orthologous With this resolved phylogeny, we further inferred directional sequence, we identified the same TRIM5-CypA2 fusion gene in M. gene flow events between macaque groups after the hybrid origin silenus (fig. S25), suggesting that this retrotranspositional event of the fascicularis group. The results obtained revealed a complex evolved in the common ancestor of the Western silenus group network of ancestral admixture in Macaca (Fig. 4A and fig. S20). (W) at least 2.17 Ma (Fig. 4A). In vitro infection experiments It is interesting to note that after the hybrid origin, various using peripheral blood mononuclear cells (PBMCs) confirmed species of the fascicularis group maintained gene flow with the de- that M. silenus can be infected with HIV-1 (Fig. 4D). scendant species of the two parental groups. Among them, the most The success or otherwise of a given hybrid speciation depends notable direction of gene flow occurred between the stump-tailed upon whether the hybridization events give rise to an established, macaque (M. arctoides) and the common ancestor of mulatta/ persistent, morphologically and ecologically distinct hybrid fuscata. On the basis of discrepancies between mitochondrial and lineage (42). Today, the species of the fascicularis group are well seg- Y chromosomal topologies, it has been previously proposed that regated from their two parental groups by either ecogeographic bar- the stump-tailed macaque originated from hybridization between riers or behavioral differences (43). The extant members of the sinica and mulatta/fuscata macaques (22). However, our analyses fascicularis group, however, display a distinctive mixture of their pa- are consistent with another phylogenomic study (38) that indicated rental species’ characteristics as a result of the ancient hybridization. the stump-tailed macaque to be a member of the sinica group with a For example, the unique bluntly bilobed and narrow penile mor- low level of introgression from mulatta/fuscata (0.05 ≤ γ ≤ 0.09; phology in males, and the bright red sexual skin with little or no table S7). In the phylogenetic trees based on low recombination evident swelling in females, are approximately intermediate rate regions and the X chromosome (figs. S14 and S15), the between the sinica and silenus groups (Fig. 5D) (23, 43). As these stump-tailed macaque also clustered with the sinica group, suggest- morphological traits are mostly associated with the reproductive ing that mitochondrial introgression from mulatta/fuscata ma- system and often have a polygenic basis, it is reasonable to expect caques occurred after the initial speciation. that genes related to these traits should also exhibit a mosaic Our phylogeny and introgression network also revealed a pattern with respect to the putative parental lineages. branching pattern within the silenus group that contradicts the pre- To identify the genetic mechanisms underlying these mixed vious morphological hypothesis, which suggested that the north- phenotypes, we used a recently developed method (44) to detect ern-tailed macaque (M. leonina) and southern-tailed macaque genes that have been subject to positive selection in the fascicularis Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 5 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 Fig. 4. Hybrid origin of the fascicularis group. (A) Schema of the reticulated evolutionary relationships between macaque species, illustrating the hybridization and major admixture events and approximate divergence times (Ma ± 1 SD) inferred from autosomal windows. The solid red line denotes the interspecific gene flow obtained from PhyloNet analyses, whereas the gray dashed line denotes the gene flow inferred from Dsuite. For the Dsuite results, we show the intersection gene flow events modeled from the two most common trees, i.e., one supporting the sister relationship between the fascicularis and sinica groups (fig. S5A), the other supporting the sister relationship between the fascicularis and silenus groups (fig. S5B). (B) Mitochondrial and (C) Y chromosomal trees. Numbers at nodes refer to bootstrap values. M. ton- keana2 and M. sylvanus2 are two additional sequenced male samples. Detailed sample information is shown in table S1. (D) Plasma viral load in PBMC aliquots from two lion-tailed macaques (Msil-001 and Msil-002) after 48 and 72 hours of infection with HIV-1NL4-3 virus. Two northern pig-tailed macaques (NPM-15216 and NPM-17214) were used as positive controls. group and each of its parental species. In addition, we applied three have become fixed in known functional domain regions (Fig. 5B criteria to filter those positively selected genes (PSGs) for the and fig. S23), whereas another two (WDR48 and PLPP1) have purpose of minimizing the potential bias of the small sample size sinica-derived mutations in the 1-kb upstream region (table S14). in each group (see Materials and Methods for details). We identified Of special interest is a gene involved in flavonoid glucuronidation 216 PSGs in the fascicularis group that were inherited from the (UGT1A9); all four sinica-derived nonsynonymous mutations were sinica lineage (table S14). Functional annotation identified 22 found to be located in the UDPGT domain region (Fig. 5B). The genes as being overrepresented in several gene regulatory networks UGT1A9 gene encodes UDP glucuronosyltransferase, which cata- related to reproductive function, such as male gamete generation lyzes the glucuronidation of endogenous estrogen hormones into [Gene Ontology (GO):0048232] and sexual reproduction water-soluble excretable metabolites (45). Because estrogen is (GO:0019953) (P < 0.001; table S15). Of these genes, seven crucial for the development of sexual swellings (46), it is certainly (UGT1A9, ADAM20, WFDC2, FANCF, DMRTC2, SIAH1, and conceivable that the shared substitutions in this gene between fas- POC1A) exhibit sinica-derived nonsynonymous substitutions that cicularis and sinica may have contributed to the reduced sexual Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 6 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 Fig. 5. The genetic basis underlying the mixture of phenotypes in the fascicularis group macaques. (A) Genomic distribution of PSGs identified in the fascicularis group. Only PSGs associated with reproductive functions are depicted. (B and C) are two examples of reproduction-related PSGs with the largest number of fixed non- synonymous mutations that were inherited from the sinica and silenus groups, respectively. The position of the mutation and predicted functional domain are both shown. (D) Diagram showing the three-pole morphocline in macaque (23, 24). Species of the fascicularis group generally exhibit a mixture of glans penis and sexual skin morphology from their parental species groups (silenus and sinica), although there is variation within the group. Examples of species are M. nemestrina (silenus group), M. mulatta ( fascicularis group), and M. thibetana (sinica group). Macaque drawings by J. Shi. swelling in females of both groups. In similar vein, the comparison DISCUSSION involving fascicularis-silenus versus the sinica group identified 195 The origin of a new species from hybridization between two preex- PSGs in the fascicularis group that were inherited from the silenus isting species is one of the most spectacular modes of speciation and lineage (table S16). Functional annotation revealed 10 to be related has consequently attracted the attention of evolutionary biologists to human reproductive phenotypes (HP0000134: female hypogo- for decades (47, 48). Although traditionally considered to be nadism; HP:0000137: abnormality of the ovaries, ovarian disease, more prevalent in plants, hybrid speciation has in recent years P < 0.005; table S17). Four of them (PRLR, PTPN11, KCNU1, and been recognized as being more common in animals than previously MSH6) harbor silenus-derived nonsynonymous mutations, fixed in thought, including examples from both invertebrates and verte- functional domain regions (Fig. 5C and fig. S24), whereas three brates (49). However, convincing evidence of hybrid speciation in (IZUMO4, BAZ1B, and SPIDR) have silenus-derived mutations all animals remains scarce, and our understanding of the underlying occurring in 1-kb upstream regions (table S16). Together, these genomic mechanisms is still quite limited. Here, our phylogenomic results provide additional genetic evidence to support the mosaic analyses using different data types and methods as well as morpho- model of reproductive morphologies observed in the fascicula- logical evidence, have provided consistent support for the hypoth- ris group. esis that the fascicularis group of macaques originated from an ancient hybridization between the sinica and silenus groups, thereby providing us with an unparalleled opportunity to Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 7 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E investigate the factors driving the origin and maintenance of hybrid Last, we have resolved several long-standing evolutionary co- species in higher animals. nundrums during the rapid speciation of macaques, involving con- On the basis of the age estimates derived from the putatively troversies surrounding phylogeny and complex ancestral neutral sequences of autosomal windows (Fig. 4A), our results hybridizations. For example, our analyses confirm that the stump- suggest that the ancient hybrid formation of the fascicularis group tailed macaque is a member of the sinica group with a low level of occurred ~3.45 to 3.56 Ma, soon after the initial separation of the introgression from mulatta/fuscata, and that the lion-tailed two parental lineages (proto-sinica and proto-silenus) ~3.86 Ma. macaque (M. silenus) and M. leonina are sister species, both of The divergence time of the two parental lineages and subsequent which can be infected by HIV. Considering the high species diver- hybridization coincided with the rapid glacio-eustatic fluctuations sity of macaques and the strong lineage structure within some in the early-middle Pliocene (50). As a very dynamic geographical species, genome sequencing data from all macaque species and dif- region in Asia, the land bridge (Isthmus of Kra) between the Malay ferent subspecies will undoubtedly help to further elucidate the peninsula and the Sunda region was repeatedly affected by glacio- complex evolutionary history and biogeography of this genus. We eustatic fluctuations during the early-middle Pliocene. We therefore are, however, confident that our conclusions regarding the hybrid speculate that the early sea level highstand (when the sea level was speciation of the fascicularis group will not change with the acqui- above the edge of the continental shelf ) may have led to the initial sition of new sequence data. Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 separation of proto-sinica and proto-silenus, while the subsequent In summary, our study reports the occurrence of an unusual lowering of the sea level facilitated the secondary contact required ancient hybridization event in primates and illustrates how specia- for hybridization. However, because no obvious postzygotic isola- tion through natural hybridization can arise via the reshuffling of tion mechanisms were observed among macaque species, the most standing genetic variation and how hybrid species may maintain insurmountable difficulties would be how hybrid integrity was their genetic integrity through the action of selection. Our study maintained and reinforced after the initial hybridization. Prezygotic provides both a strategy and a pipeline of genome analyses to iden- barriers, such as ecological divergence or/and geographic isolation, tify hybrid speciation, which should pave the way for the identifica- are thus likely to have been of greater importance in the establish- tion and exploration of further such events in the future. ment of hybrid taxa. We propose that geographic isolation may not be a necessary requirement for hybrid speciation because of the continual overlapping distribution between members of the fascicu- MATERIALS AND METHODS laris group of macaques (e.g., M. mulatta and M. fascicularis) and Sample information and ethics statement their parental group species (Fig. 1A). Nevertheless, the fascicularis A total of 12 macaque blood or DNA samples were obtained for this group macaques are now well segregated from their two parental project. Ten were used for full-genome assembly, whereas samples species by ecological barriers and/or behavioral differences (43, from the other two males, from M. tonkeana and M. sylvanus, re- 51). In particular, the fascicularis group macaques show distinctive spectively, were used for whole-genome shotgun resequencing. mixed sexual phenotypes, and the footprints of natural selection These samples were collected from multiple sources in China and have been detected in these characteristics, suggesting that assorta- Germany. Detailed information is given in table S1. All samples tive mating behavior may have played an important role in the were collected legally and in accordance with the policy of the origin and maintenance of the fascicularis group as a hybrid Animal Care and Use Ethics of Kunming Institute of Zoology (ap- species. In primates, primary and secondary sexual characteristics proval ID: SMKX-20180701-01 and SMKX-2021-01-002), which are often the target of male-female mate choice/sexual selection conforms to the regulatory standards for the human care and treat- (52, 53). Such novel combinations of sexually selected traits in the ment of animals in research. fascicularis group may have led to the emergence of novel mate pref- erences, which could have promoted the establishment of a dis- Genome sequencing, assembly, and annotation tinct lineage. The 10 newly assembled macaque genomes were sequenced on dif- We further examine the genomic makeups from the origination ferent platforms due to the protracted process of collecting tissue of a new species as the consequence of hybridization and how the samples (see table S1 for details). We used Supernova (version subsequent genomic changes in the hybrid lineage facilitated the es- 2.0.0, 10X Genomics Inc., Pleasanton, CA, USA) with default pa- tablishment of reproductive barriers toward both parental species. rameters to assemble the reads generated from the stLFR and 10X We found that the sex chromosomes (X) and regions of low recom- Genomics Chromium reads platforms. For the Nanopore long bination display more conspicuous patterns of mosaicism. General- reads, we first performed self-error correction for all Nanopore ly, these regions are considered to be less permeable to introgression long reads using NextDenovo software (v2.4.0) and then assembled as they are less permissive of foreign genes due to lower hybrid them into contigs using wtdbg-1.2.8 (57). The raw assemblies were fitness (54). However, species that have experienced rapid adaptive further polished by Illumina short reads in Pilon v 1.22 (58) three radiation provide just such an opportunity because these incompat- times under the default settings. The completeness of the new ible regions in the descendant lineages are likely to be less genetical- macaque genomes was assessed by BUSCO (v3.0.2) (59) based on ly divergent from each other, allowing them to hybridize and the mammal-specific set of 4104 single-copy orthologs produce viable offspring (29). In similar vein, these less permeable (mammalia_odb9). regions could also serve as a strong barrier against further parental Repeat elements were predicted by RepeatMasker v4.0.6 (60). introgression in hybrids and hence would be expected to exhibit For protein-coding gene annotation, we first obtained protein se- more mosaicism than other autosomal regions once hybridization quences from five well-annotated mammalian species, namely, occurred (55, 56). human (Homo sapiens, GCA_000001405.28), chimpanzee (Pan troglodytes, GCA_002880755.3), gorilla (Gorilla gorilla, Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 8 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E GCA_900006655.3), orangutan (Pongo abelii, GCA_002880775.3), set to 33, and the mitogenome of M. sylvanus (AJ309865) down- and mouse (Mus musculus, GCA_000001635.8). These protein se- loaded from the National Center for Biotechnology Information quences were mapped to each de novo genome using TBLASTN (NCBI) database was used as a starting reference. The reliability v2.2.26 (61) with an E value cutoff of 1 × 10−5. Proteins with mul- of mitochondrial contig assemblies was further validated via tiple adjacent hits were connected to each other using genBlastA BLAST searches against the reference mitogenome (AJ309865). v1.0.4 (62). We filtered out those candidate loci with homologous The D-loop and all transfer RNA (tRNA) genes were not used block lengths shorter than 30% of the length of the query protein. owing to their high mutation rate and high rate of loss. Protein- For the ab initio prediction, we used Augustus v3.0.3 (63) with op- coding genes were translated into amino acid sequences to ensure timized parameters trained from 1000 randomly selected homolo- open reading frames and to avoid NUMTs (nuclear mitochondrial gous genes. Last, we integrated all gene sets to form a DNA segments). An additional 18 macaque mitochondrial genome comprehensive and nonredundant gene set using in-house sequences from NCBI were added to expand our datasets and to Perl scripts. further validate our sample identification. The aligned genomes were partitioned into protein-coding genes and noncoding frag- Whole-genome alignments ments, and the protein-coding genes were further partitioned into We first performed pairwise whole-genome alignments against first, second, and third codon positions. PartitionFinder v2.1.6 (73) Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 Chinese rhesus macaque (M. mulatta, rheMacS) using the LASTZ was used to evaluate the best partitioning scheme under the Baye- program (27) (v 1.04.03) under the following parameters: K = 4500, sian information criterion. The ML analyses were carried out using l = 3000, Y = 15,000, E = 150, H = 2000, O = 600, T = 2. The original RAxML v8.1.15 with 1000 bootstrap replications under the best par- alignments were then processed by the chainNet package (64) to tition scheme and the GTRGAMMA model. generate the reciprocal best net alignment (default parameters, Y chromosome phylogeny except for axtChain where we used “-minScore = 5,000 -linearGap We confined our analyses to five Y-linked genes (TSPY10, SRY, ZFY, = medium”). Subsequently, we used the “maf-swap” (65) command USP9Y, and RPS4Y1) owing to the substantial technical challenges (LAST software) to sort the alignment results and obtain the presented by the sequence alignment of repetitive sequence regions optimal pairwise synteny blocks between genomes. Last, we used of the Y chromosome. Seven of our genome assemblies were from MULTIZ (v 11.2) (28) to merge all pairwise alignments into multi- males. To expand our dataset, we generated whole-genome shotgun ple genome alignments using rheMacS as the reference. sequences of two male samples from M. tonkeana and M. sylvanus, and further downloaded one male sample from M. nemestrina Phylogenetic analysis (SRR5947292). The Illumina short reads of each species were Sliding window tree and signal of conflict mapped to the Chinese rhesus macaque genome (rheMacS) using On the basis of the multiple genome alignment, we used scripts pro- BWA-MEM v0.7.12 (74) with the default settings. After obtaining vided in (66) to partition the alignment (20 autosomes and one X the bam files, we used ANGSD (75) to obtain the consensus se- chromosome) into nonoverlapping windows of 50 kb without con- quence of these genes with the following filtering parameters: -b sidering the protein-coding content. Meanwhile, we excluded bamfiles.txt -minQ 20 -minMapQ 20 -remove_bads -uniqueOnly windows that contained more than 10% gaps or hard-masked -rf region.txt -dohaplocall 1 -doCounts 1. The concatenation of repeat sequences. For each window of the alignment, an ML tree the abovementioned five genes comprised a total of 480,333 base was constructed using RAxML v8.1.15 (67) with the GTRGAMMA pairs (bp). Again, we performed the phylogenetic analyses using model and rapid bootstrapping for 100 replicates and specifying P. RAxML v.8.1.15 under the GTRGAMMA model with 1000 boot- hamadryas as the outgroup. The number of alternative tree topolo- strap replications. gies and their relative frequency were categorized by PhyBin v0.3 program (68). Conflicts between the gene trees and their frequen- Phylogenetic signal relative to the local recombination rate cies were summarized using the CONSENSE program in PHYLIP No recombination map is currently available for the reference v3.697 package (69). The coalescent species-tree was estimated from genome of the Chinese rhesus macaque (rheMacS). To construct the aforementioned sliding window trees by two summary coales- this, we used the liftover tool (76) to convert the latest fine-scale cent-based methods, ASTRAL (70) and STAR (71). Because the linkage map of Indian rhesus macaque (rheMac8) (77). First, the major disagreement among all topologies concerned the relative genome of rheMac8 was aligned to that of rheMacS using the placement of the fascicularis group (Fig. 1, B and C), for the sake LASTZ program v1.04.03 (27) with the same parameters mentioned of simplicity, we pruned the trees so as to include only five for the genome alignments. After obtaining the chain file, which ingroup species: M. sylvanus, M. thibetana, M. silenus, M. nigra, records the links of the reciprocal best orthologous regions of the and M. mulatta. These five species represented the sylvanus genome, we transferred the coordinates of the recombination map group, sinica group, silenus group (E), silenus group (W), and fas- of rheMac8 to rheMacS coordinates using the program liftover (76). cicularis group, respectively, and were selected on the basis of their Only the successfully derived liftover positions (~96.3%) were used longest scaffold N50 size. Similarly, we used PhyBin v0.3 (68) to es- in downstream analyses. timate the categories of alternative topologies and their relative fre- quencies. We also evaluated smaller (20 kb) and larger (100 kb) Hybridization and introgression analyses window sizes, but this did not change our results materially. We Although ILS certainly underpins some aspects of phylogenetic dis- therefore used the 50-kb block size for all subsequent analyses. cordance within and between clades, it is unlikely to be the only ex- Mitochondrial genome tree planation for the notable differences we observed across all window The mitochondrial genome sequence of each species was obtained trees. Therefore, we used several different methods to test for hy- from Illumina short reads using NOVOplasty 2.4 (72). K-mer was bridization and introgression events in the presence of ILS. We Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 9 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E first used the InferNetwork_MPL program in PhyloNet (37) to re- was visualized using the dtools.py script, which is provided within construct optimal phylogenetic networks between all macaque the Dsuite package. species. Owing to the extensive computational requirements, we Following the logic of Zhang et al. (33), we computed the D sta- performed this analysis on 2000 topologies that were randomly se- tistics in sliding windows to test phylogenetic hypotheses, as the lected and had a mean bootstrap support threshold of more than genomic region least affected by introgression should reflect the 80%. Three independent analyses were performed assuming reticu- true relationship. For example, if the phylogenetic relationship of lation scenarios ranging from 1 to 3, and each analysis executed 100 T1 were true, then the genomic region least affected by gene flow runs. The branch lengths and inheritance probabilities of the re- should be dominated by T1 topology. In this scenario, we set P1 turned species networks were optimized under full-likelihood by as the sinica group, P2 as the fascicularis group, and P3 as the specifying the “-po” option [full parameters: InferNetwork_MPL silenus group. The sylvanus group was treated as the outgroup (O) (all) h –b 75 -n 5 -di -po -x 100 -pl 10]. because it was always located outside the other three Asian species HyDe is another powerful and computationally efficient method groups (Fig. 1B). The same was true for T2. We used the Python to detect the level of interspecies hybridization based on site pattern script ABBABABAwindows.py from (32) to compute the D statistics frequencies of single-nucleotide polymorphism (SNP) data (35). To with window size set to 50 kb, and at least 100 biallelic SNPs were acquire high-quality SNPs for our analyses, we used the “BWA- allowed per window (-w 50,000 -m 100). Windows with absolute D- Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 GATK-SAMtools” pipeline (78) and further filtered the SNPs statistic values close to zero (e.g., low 1% windows) were regarded as using the following criteria: (i) removal of high probability miscalls the genomic regions least affected by gene flow, whereas the highest around 6 bp from predicted indels, (ii) removal of sites with consen- D-statistic values (e.g., top 1% windows) were regarded as those sus quality of <40, (iii) removal of sites with triallelic alleles and most affected by gene flow. indels, and (iv) removal of sites present in <90% of individuals. We first assessed whether ancient hybridization events had oc- Hybridization and divergence time estimates curred between species groups (sylvanus group, silenus group, To minimize the negative impact of post-speciation gene flow upon sinica group, and fascicularis group) using the “run_hyde.py” divergence time estimation, we used the genomic window sequenc- script. P. hamadryas was used as an outgroup, and sites with es that supported the most common (fig. S5A) and second most missing/ambiguous bases were ignored (--ignore_amb_sites). common (fig. S5B) trees to infer the timing of species divergence Next, an analysis at the individual level was performed using the and hybridization because they only differed in the placement of “individual_hyde.py” script to detect hybridization in individuals the fascicularis group. According to the annotation of the within species groups that had significant levels of admixture. rheMacS assembly, we filtered any genomic window that contained Last, bootstrap resampling (500 replicates in “bootstrap_hyde.py” exonic sequences, and the 10-kb regions flanking them on either script) was performed on individuals within the hybrid group to side, to reduce the potential bias introduced by selection at linked obtain a distribution of gamma values to assess heterogeneity in sites (e.g., background selection and hitchhiking) into our analysis levels of gene flow. The tract size of two parental ancestry was in- (81). Furthermore, to ensure that each genomic window was phylo- ferred by LOTER (36). The genomic distributions of each ancestry genetically independent, we only retained the windows that had a were calculated using custom R and Perl scripts. distance of at least 100 kb from each other. After applying these fil- Because HyDe is a powerful method with which to identify tering strategies, a total of 1663 windows were obtained on the au- hybrids when the parental contribution is symmetrical, we further tosomes and 103 on the X chromosome, respectively. For each applied D-statistic analysis (also known as the ABBA-BABA test) window sequence, we used the MCMCTREE program in PAML (79) to explore broader introgression patterns between lineages v.4.9e package (82) to estimate the relative divergence time for the and species. We performed these analyses in Dsuite package (80), tree imputed from that window. Two well-justified and widely used which can estimate the D statistics across all possible triplets at fossil record–based calibration points were used to calibrate evolu- one time based on the user-specified topology and can also assign tionary time (16, 17): (i) the split between African (M. sylvanus) and gene flow to specific, possibly internal, branches. Our former steps Asian macaques: 4.5 to 6.5 Ma and (ii) the split between Macacina suggested that the fascicularis group was of hybrid origin; thus, two and Papionina: 5.8 to 8.0 Ma. The prior of the overall substitution bifurcating topologies could be used: one supported the sister rela- rate was estimated using the BASEML program in PAML (82), as- tionship between the fascicularis and sinica groups (fig. S5A), suming a mean split of the most recent common ancestor of Mac- whereas the other supported the sister relationship between the fas- acina and Papionina at 7 Ma, based on the GTR + Γ substitution cicularis and silenus groups (fig. S5B). We used the program Dtrios model (rgene gamma = 1, 7.27, 1). The gamma-Dirichlet prior in Dsuite to calculate the D statistics of all possible combinations of for the rate-drift parameter (sigma2 gamma) was set to G (1, 4.5, species trios as well as species group trios. P. hamadryas was used as 1). Fossil record constraints were scaled to units of 100 Ma, and the outgroup (O) in all analyses. The resulting P values were further the constraints of minimum and maximum bounds were soft, adjusted by Benjamini and Hochberg (BH) correction for multiple with the default 2.5% probability that allowed bounds to be violated. testing bias via the p.adjust function in R. Using a significance A Markov chain Monte Carlo (MCMC) chain was first run for threshold (α = 0.05) that was BH-corrected for a total of 221 com- 1,000,000 generations as burn-in, then sampled every 500 genera- parisons (the maximum number of trios involving a given pair), z tions, until 10,000 samples had been collected. Node ages for each scores of >2.1 were considered significant. We conservatively topology were obtained using ape v5.1 package (83) in R. display only z scores of >3 as a significant signal of introgression. Results from Dtrios were further processed using the Fbranch func- Identification of PSGs under hybridization tion to generate a matrix of z scores, and a heatmap of the matrix To identify candidate genes that may have fulfilled key functions in hybrids, we used a recently developed method (44) to identify PSGs Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 10 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E in the fascicularis group and each of its parental groups (sinica and 72 hpi. The levels of viral RNA in supernatant were quantified by silenus groups). Under the null hypothesis of a neutral model, poly- a real-time polymerase chain reaction (PCR) method based on am- morphism within a lineage and divergence observed between two plification of an HIV-1NL4-3–derived Gag coding sequence de- lineages at homologous loci would be highly correlated, whereas scribed previously with slight modification (41). For each sample, positive selection would reduce levels of polymorphism relative to total cell-associated RNA was extracted by RNAiso Plus (TaKaRa) that associated with divergence. We performed two independent reagent and dissolved in 50 μl of diethyl pyrocarbonate H2O. Then, analyses. First, we grouped fascicularis and sinica together (group the cell-associated viral RNA was determined using the RNA-Direct 1) and compared them to the silenus group (group 2) to identify Real-time PCR Master Mix Kit (Toyobo). The primers were 6F (5′- the PSGs that originated from the sinica group. Similarly, to identify CATGTTTTCAGCATTATCAGAAGGA-3′) and 84R (5′- PSGs originating from the silenus group, we grouped fascicularis TGCTTGATGTCCCCCCACT-3′), and the probe was 5′-FAM- and silenus together (group 1) and compared them to the sinica CCACCCCACAAGA-TTTAAACACCATGCTAA-TAMRA-3′. group (group 2). Both coding regions and 1-kb upstream sequences PCRs were performed on the ABI Vii7 Sequence Detection System were considered in our analyses because these mutations were under conditions of 1 cycle of 95°C for 10 min, followed by 45 cycles deemed more likely to have a major effect on gene function. For of 95°C for 15 s and 60°C for 1 min. HIV-1–specific p24 antigen was each gene, we counted the number of polymorphic sites (SNPs) measured by enzyme-linked immunosorbent assay (ELISA) using Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 in group 1 (supposition A) and the number of fixed differences the HIV-1 p24 ELISA Kit (ZeptoMetrix, USA). [the SNPs with population genetic diffierentiation (FST) value > The M. silenus TRIM5 and TRIMCyp genes were amplified by 0.95] between group 1 and group 2 (supposition B) based on the TRIM5/Cyp-F: 5′-ATGGCTTCTGGAATCCTGCTTAATGTA-3′, filtered SNP file for all species. Then, a test was performed by com- TRIM5-R: 5′-TCAAGAGCTTGGTGAGCACAGAGTCA-3′, and paring the ratio of A/B to the genome-wide average A/B, which was TRIMCyp-R: 5′-TTATTCGAGTTGTCCACAGTCAGCA-3′. The calculated as the sum of A and B values across all genes analyzed. A PCR conditions were 94°C for a 3-min hot start, followed by 30 Pearson’s chi-square test on the 2 × 2 contingency table was used to cycles at 94°C for 30 s, 55°C for 30 s, 72°C for 120 s, and a final obtain the significance score that rejected the null hypothesis extension of 10 min at 72°C. The PCR products were analyzed on A(gene)/B(gene) = A(genome-wide)/B(genome-wide) (44). We an- a 1% agarose gel with ethidium bromide, purified using a DNA gel notated all SNP variants using SnpEff software (version 5.0c) (84) extraction kit (Generay Biotech, Shanghai, China), cloned into and recorded the number of fixed nonsynonymous mutations and pMD19-T simple vector (Takara, Dalian, China), and finally fixed 1-kb upstream mutations in each gene using in-house Perl sequenced. scripts. PSGs were further filtered using the following criteria: (i) the Yates’ corrected P value should be significant (<0.01) in the hybrid group and one of the parental groups, but not significant Supplementary Materials (>0.05) in another parental group; (ii) the number of fixed nonsy- This PDF file includes: Figs. S1 to S27 nonymous mutations together with 1-kb upstream mutations was Tables S1 to S17 within the top 2.5% of gene; and (iii) the topology of the gene tree should be consistent with the relationship being tested (built with the following parameters: iqtree -s SNP_data.fasta -st DNA REFERENCES AND NOTES -nt AUTO -ntmax 5 -mem 5G -bb 1000 -bnni -o Pham -m 1. O. Seehausen, Hybridization and adaptive radiation. Trends Ecol. Evol. 19, 198–207 (2004). MFP+ASC -quiet -redo). Functional enrichment analyses among 2. R. Abbott, D. Albach, S. Ansell, J. W. Arntzen, S. J. E. Baird, N. Bierne, J. Boughman, the PSGs were performed using the web-based toolkit WebGestalt A. Brelsford, C. A. Buerkle, R. Buggs, R. K. Butlin, U. Dieckmann, F. Eroukhmanoff, A. Grill, (85). The top 40 most significant function categories from the over- S. H. Cahan, J. S. Hermansen, G. Hewitt, A. G. Hudson, C. Jiggins, J. Jones, B. Keller, T. Marczewski, J. Mallet, P. Martinez-Rodriguez, M. Most, S. Mullen, R. Nichols, A. W. Nolte, representation analysis were reported (tables S15 and S17). C. Parisod, K. Pfennig, A. M. Rice, M. G. Ritchie, B. Seifert, C. M. Smadja, R. Stelkens, J. M. Szymura, R. Vainola, J. B. W. Wolf, D. Zinner, Hybridization and speciation. J. Evol. Biol. HIV infection experiment 26, 229–246 (2013). Blood samples from two lion-tailed macaques were provided by 3. L. H. Rieseberg, J. H. Willis, Plant speciation. Science 317, 910–914 (2007). G. Zoo (female Msil-001, male Msil-002), China, and used to 4. M. Schumer, G. G. Rosenthal, P. Andolfatto, How common is homoploid hybrid speciation? perform the HIV infection experiment. As a positive control, two Evolution 68, 1553–1560 (2014). 5. M. C. Melo, C. Salazar, C. D. Jiggins, M. Linares, Assortative mating preferences among northern pig-tailed macaques (NPM-15216 and NPM-17216) hybrids offers a route to hybrid speciation. Evolution 63, 1660–1665 (2009). from Kunming Primate Research Center, Kunming Institute of 6. B. L. Gross, L. H. Rieseberg, The ecological genetics of homoploid hybrid speciation. Zoology, Chinese Academy of Science were used. J. Hered. 96, 241–252 (2005). Whole blood was collected into EDTA vacutainer tubes by veni- 7. J. Mavarez, C. A. Salazar, E. Bermingham, C. Salcedo, C. D. Jiggins, M. Linares, Speciation by puncture, and PBMCs were separated by Ficoll (GE Healthcare) hybridization in Heliconius butterflies. Nature 441, 868–871 (2006). gradient centrifugation. PBMCs were activated with concanavalin 8. B. M. vonHoldt, J. A. Cahill, Z. Fan, I. Gronau, J. Robinson, J. P. Pollinger, B. Shapiro, J. Wall, R. K. Wayne, Whole-genome sequence analysis shows that two endemic species of North A (Sigma-Aldrich, USA) and interleukin-2 (IL-2) (5 U/ml). After American wolf are admixtures of the coyote and gray wolf. Sci. Adv. 2, e1501714 (2016). 72 hours of activation, the PBMCs were grown in RPMI 1640 9. S. Lamichhaney, F. Han, M. T. Webster, L. Andersson, B. R. Grant, P. R. Grant, Rapid hybrid medium supplemented with 10% fetal bovine serum (FBS) and speciation in Darwin’s finches. Science 359, 224–228 (2018). IL-2 (10 U/ml). Freshly activated PBMCs were infected with HIV- 10. J. Rogers, M. Raveendran, R. A. Harris, T. Mailund, K. Leppala, G. Athanasiadis, 1NL4-3 particles at a multiplicity of infection (MOI) of 0.05. At 4 M. H. Schierup, J. Cheng, K. Munch, J. A. Walker, M. K. Konkel, V. Jordan, C. J. Steely, hours post-infection (hpi), cells were washed four times with phos- T. O. Beckstrom, C. Bergey, A. Burrell, D. Schrempf, A. Noll, M. Kothe, G. H. Kopp, Y. Liu, S. Murali, K. Billis, F. J. Martin, M. Muffato, L. Cox, J. Else, T. Disotell, D. M. Muzny, J. Phillips- phate-buffered saline and resuspended in fresh RPMI–10% FBS Conroy, B. Aken, E. E. Eichler, T. Marques-Bonet, C. Kosiol, M. A. Batzer, M. W. Hahn, J. Tung, supplemented with IL-2. Supernatants were collected after 48 and Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 11 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E D. Zinner, C. Roos, C. J. Jolly, R. A. Gibbs, K. C. Worley; Baboon Genome Analysis Consortium, 32. S. H. Martin, J. W. Davey, C. D. Jiggins, Evaluating the use of ABBA-BABA statistics to locate The comparative genomics and complex population history of Papio baboons. Sci. Adv. 5, introgressed loci. Mol. Biol. Evol. 32, 244–257 (2015). eaau6947 (2019). 33. D. Zhang, F. E. Rheindt, H. She, Y. Cheng, G. Song, C. Jia, Y. Qu, P. Alstrom, F. Lei, Most 11. T. Zou, W. Kuang, T. Yin, L. Frantz, C. Zhang, J. Liu, H. Wu, L. Yu, Uncovering the enigmatic genomic loci misrepresent the phylogeny of an avian radiation because of ancient gene evolution of bears in greater depth: The hybrid origin of the Asiatic black bear. Proc. Natl. flow. Syst. Biol. 70, 961–975 (2021). Acad. Sci. U.S.A. 119, e2120307119 (2022). 34. M. W. Nachman, B. A. Payseur, Recombination rate variation and speciation: Theoretical 12. P. A. Hohenlohe, L. Y. Rutledge, L. P. Waits, K. R. Andrews, J. R. Adams, J. W. Hinton, predictions and empirical results from rabbits and mice. Philos. Trans. R. Soc. Lond. B Biol. R. M. Nowak, B. R. Patterson, A. P. Wydeven, P. A. Wilson, B. N. White, Comment on “Whole- Sci. 367, 409–421 (2012). genome sequence analysis shows two endemic species of North American wolf are ad- 35. P. D. Blischak, J. Chifman, A. D. Wolfe, L. S. Kubatko, HyDe: A Python package for genome- mixtures of the coyote and gray wolf”. Sci. Adv. 3, e1602250 (2017). scale hybridization detection. Syst. Biol. 67, 821–829 (2018). 13. G. G. Rosenthal, M. Schumer, P. Andolfatto, How the manakin got its crown: A novel trait 36. T. Dias-Alves, J. Mairal, M. G. B. Blum, Loter: A software package to infer local ancestry for a that is unlikely to cause speciation. Proc. Natl. Acad. Sci. U.S.A. 115, E4144–E4145 (2018). wide range of species. Mol. Biol. Evol. 35, 2318–2326 (2018). 14. C. Scornavacca, N. Galtier, Incomplete lineage sorting in mammalian phylogenomics. Syst. 37. Y. Yu, L. Nakhleh, A maximum pseudo-likelihood approach for phylogenetic networks. Biol. 66, 112–120 (2017). BMC Genomics 16 (suppl. 10), S10 (2015). 15. P. Fan, Y. Liu, Z. Zhang, C. Zhao, C. Li, W. Liu, Z. Liu, M. Li, Phylogenetic position of the white- 38. Z. Fan, A. Zhou, N. Osada, J. Yu, J. Jiang, P. Li, L. Du, L. Niu, J. Deng, H. Xu, J. Xing, B. Yue, J. Li, cheeked macaque (Macaca leucogenys), a newly described primate from southeastern Ancient hybridization and admixture in macaques (genus Macaca) inferred from whole Tibet. Mol. Phylogenet. Evol. 107, 80–89 (2017). genome sequences. Mol. Phylogenet. Evol. 127, 376–386 (2018). 16. E. Delson, Fossil macaques, phyletic relationships and a scenario of deployment, in The 39. C. Abegg, B. Thierry, Macaque evolution and dispersal in insular South-East Asia. Biol. Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 Macaques: Studies in Ecology, Behavior, and Evolution, D. G. Lindburg, Ed. (van Nostran- J. Linn. Soc. 75, 555–576 (2002). d–Reinhold, 1980), pp. 10–30. 40. C. H. Liao, Y. Q. Kuang, H. L. Liu, Y. T. Zheng, B. Su, A novel fusion gene, TRIM5-Cyclophilin A 17. C. Roos, M. Kothe, D. M. Alba, E. Delson, D. Zinner, The radiation of macaques out of Africa: in the pig-tailed macaque determines its susceptibility to HIV-1 infection. AIDS 21 (Suppl Evidence from mitogenome divergence times and the fossil record. J. Hum. Evol. 133, 8), S19–S26 (2007). 114–132 (2019). 41. W. Pang, G.-H. Zhang, J. Jiang, H.-Y. Zheng, L.-T. Zhang, X.-L. Zhang, J.-H. Song, M.-X. Zhang, 18. M. B. Gardner, P. A. Luciw, Macaque models of human infectious disease. ILAR J. 49, J.-W. Zhu, A.-H. Lei, R.-R. Tian, X.-M. Liu, L. Zhang, G. Gao, L. Su, Y.-T. Zheng, HIV-1 can infect 220–255 (2008). northern pig-tailed macaques (Macaca leonina) and form viral reservoirs in vivo. Sci. Bull. 19. E. L. Bynum, D. Z. Bynum, J. Supriatna, Confirmation and location of the hybrid zone (Beijing) 62, 1315–1324 (2017). between wild populations of Macaca tonkeana and Macaca hecki in central Sulawesi, In- 42. G. N. Feliner, I. Alvarez, J. Fuertes-Aguilar, M. Heuertz, I. Marques, F. Moharrek, R. Pineiro, donesia. Am. J. Primatol. 43, 181–209 (1997). R. Riina, J. A. Rossello, P. S. Soltis, I. Villa-Machio, Is homoploid hybrid speciation that rare? 20. Y. Hamada, N. Urasopon, I. Hadi, S. Malaivijitnond, Body size and proportions and pelage An empiricist’s view. Heredity (Edinb.) 118, 513–516 (2017). color of free-ranging Macaca mulatta from a zone of hybridization in Northeastern Thai- 43. J. Fooden, Comparative review of fascicularis-group species of macaques (primates: land. Int. J. Primatol. 27, 497–513 (2006). Macaca). Fieldiana Zool. 2006, 1–43 (2006). 21. Y. Song, C. Jiang, K.-H. Li, J. Li, H. Qiu, M. Price, Z.-X. Fan, J. Li, Genome-wide analysis reveals 44. Z. Wang, Y. Jiang, H. Bi, Z. Lu, Y. Ma, X. Yang, N. Chen, B. Tian, B. Liu, X. Mao, T. Ma, signatures of complex introgressive gene flow in macaques (genus Macaca). Zool. Res. 42, S. P. DiFazio, Q. Hu, R. J. Abbott, J. Q. Liu, Hybrid speciation via inheritance of alternate 433–449 (2021). alleles of parental isolating genes. Mol. Plant 14, 208–222 (2021). 22. A. J. Tosi, J. C. Morales, D. J. Melnick, Paternal, maternal, and biparental molecular markers 45. J. Lepine, O. Bernard, M. Plante, B. Tetu, G. Pelletier, F. Labrie, A. Belanger, C. Guillemette, provide unique windows onto the evolutionary history of macaque monkeys. Evolution 57, Specificity and regioselectivity of the conjugation of estradiol, estrone, and their cate- 1419–1435 (2003). cholestrogen and methoxyestrogen metabolites by human uridine diphospho-glucuro- 23. J. E. Fa, The genus Macaca: A review of taxonomy and evolution. Mamm. Rev. 19, nosyltransferases expressed in endometrium. J. Clin. Endocrinol. Metab. 89, 45–81 (1989). 5222–5232 (2004). 24. J. Fooden, Provisional classification and key to living species of macaques (primates: 46. L. R. Gesquiere, E. O. Wango, S. C. Alberts, J. Altmann, Mechanisms of sexual selection: Macaca). Folia Primatol. (Basel) 25, 225–236 (1976). Sexual swellings and estrogen concentrations as fertility indicators and cues for male 25. Y. He, X. Luo, B. Zhou, T. Hu, X. Meng, P. A. Audano, Z. N. Kronenberg, E. E. Eichler, J. Jin, consort decisions in wild baboons. Horm. Behav. 51, 114–125 (2007). Y. Guo, Y. Yang, X. Qi, B. Su, Long-read assembly of the Chinese rhesus macaque genome 47. R. J. Abbott, M. J. Hegarty, S. J. Hiscock, A. C. Brennan, Homoploid hybrid speciation in and identification of ape-specific structural variants. Nat. Commun. 10, 4233 (2019). action. Taxon 59, 1375–1386 (2010). 26. Y. Shao, L. Zhou, F. Li, L. Zhao, B.-L. Zhang, F. Shao, J.-W. Chen, C.-Y. Chen, X.-P. Bi, X.- 48. J. Mallet, Hybrid speciation. Nature 446, 279–283 (2007). L. Zhuang, H.-L. Zhu, J. Hu, Z. Sun, X. Li, D. Wang, I. Rivas-González, S. Wang, Y.-M. Wang, 49. J. Mavarez, M. Linares, Homoploid hybrid speciation in animals. Mol. Ecol. 17, W. Chen, G. Li, H.-M. Lu, Y. Liu, L. Kuderna, K. Farh, P.-F. Fan, L. Yu, M. Li, Z.-J. Liu, G. P. Tiley, 4181–4185 (2008). A. D. Yoder, C. Roos, T. Hayakawa, T. Marques-Bonet, J. Rogers, D. N. Cooper, H. Li, 50. D. S. Woodruff, Neogene marine transgressions, palaeogeography and biogeographic M. H. Schierup, Y.-G. Yao, Y.-P. Zhang, W. Wang, X.-G. Qi, G. Zhang, D.-D. Wu, Phylogenomic transitions on the Thai–Malay Peninsula. J. Biogeogr. 30, 551–567 (2003). analyses provide insights into primate genomic and phenotypic evolution. Science, (2021). 51. J. Fooden, Ecogeographic segregation of macaque species. Primates 23, 574–579 (1982). 27. R. S. Harris, “Improved pairwise alignment of genomic DNA,” thesis, The Pennsylvania State 52. T. Caro, K. Brockelsby, A. Ferrari, M. Koneru, K. Ono, E. Touche, T. Stankowich, The evolution University (2007). of primate coloration revisited. Behav. Ecol. 32, 555–567 (2021). 28. M. Blanchette, W. J. Kent, C. Riemer, L. Elnitski, A. F. Smit, K. M. Roskin, R. Baertsch, 53. C. Dubuc, L. J. N. Brent, A. K. Accamando, M. S. Gerald, A. MacLarnon, S. Semple, K. Rosenbloom, H. Clawson, E. D. Green, D. Haussler, W. Miller, Aligning multiple genomic M. Heistermann, A. Engelhardt, Sexual skin color contains information about the timing of sequences with the threaded blockset aligner. Genome Res. 14, 708–715 (2004). the fertile phase in free-ranging Macaca mulatta. Int. J. Primatol. 30, 777–789 (2009). 29. N. B. Edelman, P. B. Frandsen, M. Miyagi, B. Clavijo, J. Davey, R. B. Dikow, G. Garcia-Accinelli, 54. O. Seehausen, R. K. Butlin, I. Keller, C. E. Wagner, J. W. Boughman, P. A. Hohenlohe, S. M. Van Belleghem, N. Patterson, D. E. Neafsey, R. Challis, S. Kumar, G. R. P. Moreira, C. L. Peichel, G. P. Saetre; C. Bank, A. Brannstrom, A. Brelsford, C. S. Clarkson, C. Salazar, M. Chouteau, B. A. Counterman, R. Papa, M. Blaxter, R. D. Reed, F. Eroukhmanoff, J. L. Feder, M. C. Fischer, A. D. Foote, P. Franchini, C. D. Jiggins, F. C. Jones, K. K. Dasmahapatra, M. Kronforst, M. Joron, C. D. Jiggins, W. O. McMillan, F. Di Palma, A. K. Lindholm, K. Lucek, M. E. Maan, D. A. Marques, S. H. Martin, B. Matthews, J. I. Meier, A. J. Blumberg, J. Wakeley, D. Jaffe, J. Mallet, Genomic architecture and introgression shape M. Most, M. W. Nachman, E. Nonaka, D. J. Rennison, J. Schwarzer, E. T. Watson, a butterfly radiation. Science 366, 594–599 (2019). A. M. Westram, A. Widmer, Genomics and the origin of species. Nat. Rev. Genet. 15, 30. B. Charlesworth, J. A. Coyne, N. H. Barton, The relative rates of evolution of Sex chromo- 176–192 (2014). somes and autosomes. Am. Nat. 130, 113–146 (1987). 55. T. O. Elgvin, C. N. Trier, O. K. Torresen, I. J. Hagen, S. Lien, A. J. Nederbragt, M. Ravinet, 31. M. C. Fontaine, J. B. Pease, A. Steele, R. M. Waterhouse, D. E. Neafsey, I. V. Sharakhov, H. Jensen, G. P. Saetre, The genomic mosaicism of hybrid speciation. Sci. Adv. 3, X. F. Jiang, A. B. Hall, F. Catteruccia, E. Kakani, S. N. Mitchell, Y. C. Wu, H. A. Smith, R. R. Love, e1602996 (2017). M. K. Lawniczak, M. A. Slotman, S. J. Emrich, M. W. Hahn, N. J. Besansky, Extensive intro- 56. K. Kunte, C. Shea, M. L. Aardema, J. M. Scriber, T. E. Juenger, L. E. Gilbert, M. R. Kronforst, Sex gression in a malaria vector species complex revealed by phylogenomics. Science 347, chromosome mosaicism and hybrid speciation among tiger swallowtail butterflies. PLOS 1258524 (2015). Genet. 7, e1002274 (2011). Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 12 of 13 S C I E N C E A D VA N C E S | R E S E A R C H A R T I C L E 57. J. Ruan, H. Li, Fast and accurate long-read assembly with wtdbg2. Nat. Methods 17, 78. A. O. Ayoola, B.-L. Zhang, R. P. Meisel, L. M. Nneji, Y. Shao, O. B. Morenikeji, A. C. Adeola, 155–158 (2020). S. I. Ng’ang’a, B. G. Ogunjemite, A. O. Okeyoyin, C. Roos, D.-D. Wu, Population genomics 58. B. J. Walker, T. Abeel, T. Shea, M. Priest, A. Abouelliel, S. Sakthikumar, C. A. Cuomo, Q. Zeng, reveals incipient speciation, introgression, and adaptation in the African mona monkey J. Wortman, S. K. Young, A. M. Earl, Pilon: An integrated tool for comprehensive microbial (Cercopithecus mona). Mol. Biol. Evol. 38, 876–890 (2021). variant detection and genome assembly improvement. PLOS ONE 9, e112963 (2014). 79. N. Patterson, P. Moorjani, Y. Luo, S. Mallick, N. Rohland, Y. Zhan, T. Genschoreck, T. Webster, 59. F. A. Simao, R. M. Waterhouse, P. Ioannidis, E. V. Kriventseva, E. M. Zdobnov, BUSCO: As- D. Reich, Ancient admixture in human history. Genetics 192, 1065–1093 (2012). sessing genome assembly and annotation completeness with single-copy orthologs. Bi- 80. M. Malinsky, M. Matschiner, H. Svardal, Dsuite—Fast D-statistics and related admixture oinformatics 31, 3210–3212 (2015). evidence from VCF files. Mol. Ecol. Resour. 21, 584–595 (2021). 60. G. Benson, Tandem repeats finder: A program to analyze DNA sequences. Nucleic Acids Res. 81. A. H. Freedman, I. Gronau, R. M. Schweizer, D. Ortega-Del Vecchyo, E. Han, P. M. Silva, 27, 573–580 (1999). M. Galaverni, Z. Fan, P. Marx, B. Lorente-Galdos, H. Beale, O. Ramirez, F. Hormozdiari, 61. S. F. Altschul, W. Gish, W. Miller, E. W. Myers, D. J. Lipman, Basic local alignment search tool. C. Alkan, C. Vila, K. Squire, E. Geffen, J. Kusak, A. R. Boyko, H. G. Parker, C. Lee, V. Tadigotla, J. Mol. Biol. 215, 403–410 (1990). A. Wilton, A. Siepel, C. D. Bustamante, T. T. Harkins, S. F. Nelson, E. A. Ostrander, T. Marques- 62. R. She, J. S.-C. Chu, K. Wang, J. Pei, N. Chen, genBlastA: Enabling BLAST to identify ho- Bonet, R. K. Wayne, J. Novembre, Genome sequencing highlights the dynamic early history mologous gene sequences. Genome Res. 19, 143–149 (2009). of dogs. PLOS Genet. 10, e1004016 (2014). 63. O. Keller, M. Kollmar, M. Stanke, S. Waack, A novel hybrid gene prediction method em- 82. Z. Yang, PAML 4: Phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, ploying protein multiple sequence alignments. Bioinformatics 27, 757–763 (2011). 1586–1591 (2007). 64. W. J. Kent, R. Baertsch, A. Hinrichs, W. Miller, D. Haussler, Evolution’s cauldron: Duplication, 83. E. Paradis, K. Schliep, ape 5.0: An environment for modern phylogenetics and evolutionary deletion, and rearrangement in the mouse and human genomes. Proc. Natl. Acad. Sci. U.S.A. analyses in R. Bioinformatics 35, 526–528 (2019). Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 100, 11484–11489 (2003). 84. P. Cingolani, A. Platts, L. L. Wang, M. Coon, T. Nguyen, L. Wang, S. J. Land, X. Lu, D. M. Ruden, 65. S. M. Kielbasa, R. Wan, K. Sato, P. Horton, M. C. Frith, Adaptive seeds tame genomic se- A program for annotating and predicting the effects of single nucleotide polymorphisms, quence comparison. Genome Res. 21, 487–493 (2011). SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 66. L. Chen, Q. Qin, Y. Jiang, K. Wang, Z. Lin, Z. Li, F. Bibi, Y. Yang, J. Wang, W. Nie, W. Su, G. Liu, (Austin) 6, 80–92 (2012). Q. Li, W. Fu, X. Pan, C. Liu, J. Yang, C. Zhang, Y. Yin, Y. Wang, Y. Zhao, C. Zhang, Z. Wang, 85. Y. Liao, J. Wang, E. J. Jaehnig, Z. Shi, B. Zhang, WebGestalt 2019: Gene set analysis toolkit Y. Qin, W. Liu, B. Wang, Y. Ren, R. Zhang, Y. Zeng, R. R. da Fonseca, B. Wei, R. Li, W. Wan, with revamped UIs and APIs. Nucleic Acids Res. 47, W199–W205 (2019). R. Zhao, W. Zhu, Y. Wang, S. Duan, Y. Gao, Y. Zhang, C. Chen, C. Hvilsom, C. W. Epps, L. G. Chemnick, Y. Doug, S. Mirarab, H. R. Siegismund, O. A. Ryder, M. T. P. Gilbert, Acknowledgments: We thank W. Wang, X. Ren, and X. Li for their early assistance with this H. A. Lewin, G. Zhang, R. Heller, W. Wang, Large-scale ruminant genome sequencing project. Funding: This work was funded by the National Key Research and Development provides insights into their evolution and distinct traits. Science 364, eaav6202 (2019). Program of China (nos. 2021YFF0702700 and 2022YFF0710901), the National Natural Science 67. A. Stamatakis, RAxML version 8: A tool for phylogenetic analysis and post-analysis of large Foundation of China (no. 31822048 to D.-D.W. and no. 32270500 to B.-L.Z.), Yunnan Applied phylogenies. Bioinformatics 30, 1312–1313 (2014). Basic Research Projects (no. 2019FI010 to D.-D.W.), CAS Light of West China Program (xbzg- 68. R. R. Newton, I. L. Newton, PhyBin: Binning trees by topology. PeerJ 1, e187 (2013). zdsys-202213) and the Animal Branch of the Germplasm Bank of Wild Species of Chinese Academy of Science (the Large Research Infrastructure Funding). This work was also partially 69. J. Felsenstein, PHYLIP—Phylogeny inference package (version 3.2). Cladistics 5, supported by a Villum Investigator Grant to G.Z. (no. 25900) and the National Natural Science 164–166 (1989). Foundation of China to Y.-T.Z. (U1802284). T.M.-B. was supported by funding from the European 70. S. Mirarab, R. Reaz, M. S. Bayzid, T. Zimmermann, M. S. Swenson, T. Warnow, ASTRAL: Research Council (ERC) under the European Union’s Horizon 2020 research and innovation Genome-scale coalescent-based species tree estimation. Bioinformatics 30, programme (no. 864203), BFU2017-86471-P (MINECO/FEDER, UE), and Howard Hughes i541–i548 (2014). International Early Career. Author contributions: Conceptualization and project 71. L. Liu, L. Yu, D. K. Pearl, S. V. Edwards, Estimating species phylogenies using coalescence administration: D.-D.W., G.Z., and Y.-T.Z. Data analysis and visualization: B.-L.Z., Z.W., S.W., Y.S., times among sequences. Syst. Biol. 58, 468–477 (2009). Y.D., L.Z., J.C., and L.W. HIV infection experiments: W.P., M.-T.L., and W.-Q.H. Sample collection 72. N. Dierckxsens, P. Mardulyn, G. Smits, NOVOPlasty: De novo assembly of organelle and provider: W.C., M.-M.Y., Y.W., H.F.-B., S.M., H.M., F.W., L.K., T.M.-B., and C.R. Writing—original genomes from whole genome data. Nucleic Acids Res. 45, e18 (2017). draft: B.-L.Z. and D.-D.W. Writing—review and editing: B.-L.Z., D.-D.W., G.Z., Y.-T.Z., J.L., D.N.C., 73. R. Lanfear, B. Calcott, S. Y. W. Ho, S. Guindon, Partitionfinder: Combined selection of par- M.H.S., Z.L., M.L., X.-G.Q., C.R., T.M.-B., H.M., and Z.W. Competing interests: L.K. is an employee titioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 29, of Illumina Inc. The authors declare that they have no competing interests. Data and materials 1695–1701 (2012). availability: All data needed to evaluate the conclusions in the paper are present in the paper 74. H. Li, Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. and/or the Supplementary Materials. The raw sequencing data from this project were deposited arXiv:1303.3997v2 [q-bio.GN] (16 March 2013). in the Genome Warehouse (GWH) database under project accession number PRJCA007326. 75. T. S. Korneliussen, A. Albrechtsen, R. Nielsen, ANGSD: Analysis of next generation se- quencing data. BMC Bioinform. 15, 356 (2014). Submitted 6 June 2022 76. R. M. Kuhn, D. Haussler, W. J. Kent, The UCSC Genome Browser and associated tools. Brief. Accepted 25 January 2023 Bioinform. 14, 144–161 (2013). Published 1 June 2023 77. C. Xue, N. Rustagi, X. Liu, M. Raveendran, R. A. Harris, M. G. Venkata, J. Rogers, F. Yu, 10.1126/sciadv.add3580 Reduced meiotic recombination in rhesus macaques and the origin of the human re- combination landscape. PLOS ONE 15, e0236285 (2020). Zhang et al., Sci. Adv. 9, eadd3580 (2023) 1 June 2023 13 of 13 Comparative genomics reveals the hybrid origin of a macaque group Bao-Lin Zhang, Wu Chen, Zefu Wang, Wei Pang, Meng-Ting Luo, Sheng Wang, Yong Shao, Wen-Qiang He, Yuan Deng, Long Zhou, Jiawei Chen, Min-Min Yang, Yajiang Wu, Lu Wang, Hugo Fernndez-Bellon, Sandra Molloy, Hlne Meunier, Fanlie Wanert, Lukas Kuderna, Tomas Marques-Bonet, Christian Roos, Xiao-Guang Qi, Ming Li, Zhijin Liu, Mikkel Heide Schierup, David N. Cooper, Jianquan Liu, Yong-Tang Zheng, Guojie Zhang, and Dong-Dong Wu Sci. Adv., 9 (22), eadd3580. DOI: 10.1126/sciadv.add3580 Downloaded from https://www.science.org at Universitat Autnoma de Barcelona on June 12, 2023 View the article online https://www.science.org/doi/10.1126/sciadv.add3580 Permissions https://www.science.org/help/reprints-and-permissions Use of this article is subject to the Terms of service Science Advances (ISSN ) is published by the American Association for the Advancement of Science. 1200 New York Avenue NW, Washington, DC 20005. The title Science Advances is a registered trademark of AAAS. Copyright © 2023 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).