Complex Evolutionary History With Extensive Ancestral Gene Flow in an African Primate Radiation Axel Jensen ,1,* Frances Swift,2 Dorien de Vries,3 Robin M.D. Beck,3 Lukas F.K. Kuderna,4 Sascha Knauf,5 Idrissa S. Chuma,6 Julius D. Keyyu,7 Andrew C. Kitchener ,8,9 Kyle Farh,4 Jeffrey Rogers ,10 Tomas Marques-Bonet,11,12,13,14,15 Kate M. Detwiler,16 Christian Roos,17 and Katerina Guschanski 1,2,* 1 Department of Ecology and Genetics, Animal Ecology, Uppsala University, Uppsala SE-75236, Sweden 2 School of Biological Sciences, Institute of Ecology and Evolution, University of Edinburgh, Edinburgh, UK 3 School of Science, Engineering & Environment, University of Salford, Salford M5 4WT, UK 4 Illumina Artificial Intelligence Laboratory, Illumina Inc., Foster City, CA 94404, USA Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 5 Institute of International Animal Health/One Health, Friedrich-Loeffler-Institut, Federal Research Institute for Animal Health, Greifswald – Insel Riems 17493, Germany 6 Tanzania National Parks, Arusha, Tanzania 7 Tanzania Wildlife Research Institute (TAWIRI), Arusha, Tanzania 8 Department of Natural Sciences, National Museums Scotland, Edinburgh EH1 1JF, UK 9 School of Geosciences, University of Edinburgh, Edinburgh EH8 9XP, UK 10 Human Genome Sequencing Center and Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, TX 77030, USA 11 Institute of Evolutionary Biology (UPF-CSIC), PRBB, Barcelona 08003, Spain 12 Institut Català de Paleontologia Miquel Crusafont, Universitat Autò noma de Barcelona, Barcelona, Spain 13 Catalan Institution of Research and Advanced Studies (ICREA), Barcelona, Spain 14 CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Barcelona 08028, Spain 15 Institució Catalana de Recerca i Estudis Avanç ats (ICREA) and Universitat Pompeu Fabra, Barcelona 08010, Spain 16 Department of Biological Sciences, Florida Atlantic University, Boca Raton, FL, USA 17 Gene Bank of Primates and Primate Genetics Laboratory, German Primate Center, Leibniz Institute for Primate Research, Gö ttingen 37077, Germany *Corresponding authors: E-mails: axel.jensen@ebc.uu.se; katerina.guschanski@ed.ac.uk. Associate editor: Mary O’Connell Article Abstract Understanding the drivers of speciation is fundamental in evolutionary biology, and recent studies highlight hybrid­ ization as an important evolutionary force. Using whole-genome sequencing data from 22 species of guenons (tribe Cercopithecini), one of the world's largest primate radiations, we show that rampant gene flow characterizes their evolutionary history and identify ancient hybridization across deeply divergent lineages that differ in ecology, morphology, and karyotypes. Some hybridization events resulted in mitochondrial introgression between distant lineages, likely facilitated by cointrogression of coadapted nuclear variants. Although the genomic landscapes of introgression were largely lineage specific, we found that genes with immune functions were overrepresented in in­ trogressing regions, in line with adaptive introgression, whereas genes involved in pigmentation and morphology may contribute to reproductive isolation. In line with reports from other systems that hybridization might facilitate diversification, we find that some of the most species-rich guenon clades are of admixed origin. This study provides important insights into the prevalence, role, and outcomes of ancestral hybridization in a large mammalian radiation. Key words: phylogenomics, speciation, ancient hybridization, evolutionary genomics, mammalian radiation, guenon. Received: July 05, 2023. Revised: October 17, 2023. Accepted: November 09, 2023 © The Author(s) 2023. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https:// creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Open Access Mol. Biol. Evol. 40(12):msad247 https://doi.org/10.1093/molbev/msad247 Advance Access publication November 21, 2023 1 Jensen et al. · https://doi.org/10.1093/molbev/msad247 MBE Introduction Results Ancient hybridization has been reported in many organ­ Dataset, Sequencing and Genotyping isms, including mammals (Taylor and Larson 2019). We compiled a dataset of whole-genome sequences from 37 However, owing to their large genomes, studies of entire samples belonging to 24 primate species (supplementary mammalian radiations, particularly among species-rich table S1, Supplementary Material online; van der Valk et al. groups, are underrepresented (but see Gopalakrishnan 2020; Ayoola et al. 2021; Kuderna et al. 2023). These included et al. 2018; Chavez et al. 2022), and most cases focus on 22 guenon species from all 6 genera, Allenopithecus, pairs of species (Taylor and Larson 2019). Yet, large radia­ Allochrocebus, Chlorocebus, Erythrocebus, Miopithecus, and tions with lineages of different ages offer a unique oppor­ Cercopithecus (Fig. 1A), which collectively show a high degree tunity to understand how increasing lineage divergence of sympatry throughout sub-Saharan Africa (Fig. 1B). The along a speciation continuum may impact the ability to most species-rich genus, Cercopithecus, is commonly divided hybridize and the outcome of hybridization events. As re­ Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 into 6 evolutionary distinct species groups: cephus, mitis, productive isolation accumulates with genetic distance mona, neglectus, diana, and hamlyni (Grubb et al. 2003; (Coyne and Orr 2004), the impact and consequences Lo Bianco et al. 2017), all represented in our dataset. We of introgression may vary, as well as the mechanisms also included 2 outgroup species: Macaca mulatta (rhesus that allow for genomic incompatibilities to be overcome. macaque) and Cercocebus torquatus (red-capped manga­ Several studies have reported that introgression between bey). Average read mapping depth varied between 15.7 divergent lineages may spark rapid speciation and adaptive and 57.7 (median = 29.8). We called genotypes against radiations by introducing novel genetic variation (Meier the rhesus macaque reference genome (Mmul_10, et al. 2017; Marques et al. 2019; Svardal et al. 2020). GenBank: GCA_014858485.1) and after stringent filtering Comprehensive genomic datasets from species-rich radia­ (Materials and Methods), we obtained 1.18 billion geno­ tions offer unique possibilities to study the interplay be­ typed sites across the autosomes, 52 million sites on the tween hybridization, adaptation, and speciation. X-chromosome, and 63,565 sites on the Y-chromosome. Here, we focus on guenons (tribe Cercopithecini), a Of these, 140 and 4.2 million sites were biallelic single nu­ species-rich group of African primates that radiated over cleotide polymorphisms (SNPs) on the autosomes and the the last ca. 10 to 15 million years (MY) (Guschanski et al. X-chromosome, respectively. Additionally, we assembled 2013; Kuderna et al. 2023). With 89 taxa and over 30 dis­ and annotated the mitochondrial genomes (mtDNA) of tinct species (IUCN 2022), guenons represent one of the all samples. world's largest primate radiations, spanning a broad range of divergence times. Guenons are renowned for their eco­ Species Tree Inferences logical, morphological, and karyotypic diversity and have We constructed a multispecies coalescence tree in attracted the attention of evolutionary biologists and ecol­ ASTRAL (Zhang et al. 2018), using 3,346 independent ogists for decades (Dutrillaux et al. 1988; Glenn and Cords autosomal gene trees, each representing a 25 kb genomic 2002; Grubb et al. 2003; Enstam and Isbell 2007; Moulin region sampled every 500 kb along the genome (Fig. 1C; et al. 2008). Yet, despite possessing multiple characteristics supplementary fig. S1, Supplementary Material online; that act as reproductive barriers in other study systems, Materials and Methods). All genera, species groups and guenons readily hybridize, even across deep evolutionary species were monophyletic with maximal local posterior distances (Detwiler et al. 2005; de Jong and Butynski probability (lpp) support (lpp = 1, supplementary fig. S1, 2010; Detwiler 2019). For example, Detwiler (2019) reports Supplementary Material online). The tree topology was viable and at least partially fertile hybrids between largely consistent with the recent primate phylogeny based Cercopithecus mitis and Cercopithecus ascanius, 2 species on ultraconserved elements (UCE; Kuderna et al. 2023), ex­ that diverged ca. 5 million years ago (MYA) (Kuderna cept for the placement of the genus Erythrocebus and the et al. 2023) and differ in chromosome numbers. Genomic diana species group. Both these discrepancies involve rapid studies have also identified ancient gene flow in several diversifications with short internal branches and strong in­ guenon lineages (Svardal et al. 2017; van der Valk et al. complete lineage sorting (ILS), as illustrated by a high degree 2020; Ayoola et al. 2021), but the extent and role of ances­ of gene tree discordance around these nodes (normalized tral hybridization throughout the clade is unknown. quartet score = 38% and 44%, respectively), which may be As such, guenons provide a highly informative system to difficult to resolve using UCE loci due to their slow evolu­ study the evolution of reproductive isolation and speci­ tion (Bejerano et al. 2004). ation in the context of a large radiation. Specifically, we We estimated divergence times on the ASTRAL topology aim to (i) study genomic patterns of introgression and me­ by combining 4 independent MCMCTree runs (Yang 2007), chanisms that allow barriers to gene flow to be overcome which converged at highly similar node age estimates and at different evolutionary distances, (ii) investigate the re­ jointly resulted in effective sample sizes (ESS) ≥ 285 for all peatability of introgression landscapes along the speci­ nodes (supplementary fig. S2, Supplementary Material on­ ation continuum and (iii) explore the functional role of line). These analyses suggest that the guenons split from introgressed regions and identify loci contributing to re­ the Papionini tribe ca. 14.6 MYA (Fig. 1C; supplementary productive isolation. table S2, Supplementary Material online), followed by the 2 Complex Primate Evolutionary History With Extensive Ancestral Gene Flow · https://doi.org/10.1093/molbev/msad247 MBE A B Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 C. hamlyni C. diana C C. diana C. roloway 88 C. roloway C. neglectus 44 C. neglectus C. pogonias 45 C. pogonias C. lowei 76 C. lowei C. mona C. mona E. patas C. albogularis C. albogularis C. mitis C. mitis C. nictitans C. nictitans C. petuarista A. preussi C. cephus A. lhoesti C. ascanius C. petuarista M. ogouensis C. cephus Allen. nigroviridis C. ascanius E. patas Chl. pygerythrus 64 38 Chl. sabaeus Chl. dryas Chl. pygerythrus Chl. sabaeus Chl. dryas A. solatus A. preussi C. hamlyni A. lhoesti Allen. nigroviridis A. solatus M. ogouensis M.mulatta M. mulatta Cercoc.torquatus Cercoc. torquatus MYA 14 12 10 8 6 4 2 2 4 6 8 10 12 14 MYA Fig. 1. Taxonomy, species richness and mitonuclear discordances among guenons. (A) Taxonomic overview of the 22 species included in this study, shown according to genus and species groups. The color of the vertical bars corresponds to branch colors in (C). (B) Species richness/ degree of sympatry for species included in this study, based on species distributions from IUCN (2022). (C) Astral species tree obtained from 3,346 autosomal gene trees (left) and maximum likelihood tree constructed with RaxML from complete mitochondrial genomes (right), with connectors highlighting phylogenetic discordances. Node annotations show branch quartet support (left) or bootstrap support (right), where this was <95%. Branches are colored based on genus/species groups as in (A). Tree topologies were estimated from all available samples (supplementary figs. S1 and S7, Supplementary Material online) and were subsequently pruned to a single sample per species prior to divergence date estimates with MCMCTree, applying fossil calibrations to the nodes annotated with asterisks. Nodes were rotated to aid visualization of mitonuclear discordances. Photo credit for (A): hamlyni: NRowe/alltheworldsprimates.org; diana, Chlorocebus: K. Guschanski; neglectus: M. D’haen; mona: S. Knauf; mitis: M. Mpongo & K. Detwiler; cephus: S. Crawford & K. Detwiler; Allochrocebus: T. Ukizintambara; Erythrocebus: T. Valkenburg; Miopithecus: P. Paixão; Allenopithecus: D. Sutherland. 3 Jensen et al. · https://doi.org/10.1093/molbev/msad247 MBE first split within the guenons ca. 13.4 MYA, dividing them (supplementary fig. S9, Supplementary Material online). into 2 major clades. One of these clades contains the genera For example, while Erythrocebus patas shifted position Miopithecus (M) and Cercopithecus (C), which diverged across a single node in 19% of 10,000 simulations, only a from each other ca. 12.4 MYA. The mainly arboreal species single simulation resulted in a shift across 2 nodes, which groups belonging to the genus Cercopithecus radiated ca. 8 would still be insufficient to explain its placement in to 11 MYA, referred to as the “Arboreal clade” hereafter. the real data (supplementary fig. S9B, Supplementary The other major clade contains the genus Allenopithecus Material online). Thus, ILS remains a plausible explanation (Allen), as sister to the terrestrial genera Chlorocebus (Chl), for lineages shifting positions across single nodes, e.g. Erythrocebus (E), and Allochrocebus (A), hereafter termed Miopithecus and Cercopithecus diana, whereas the deeper the “Terrestrial clade.” The Terrestrial clade split from mitonuclear discordances are more likely explained by hy­ Allenopithecus ca. 12.5 MYA and underwent rapid radiation bridization and mitochondrial introgression. ca. 8 MYA. Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 We repeated the same process with 118 X-chromosomal loci, resulting in an identical topology with highly Ancient Gene Flow is Prevalent Among Guenons similar divergence times (supplementary figs. S3 to S5, To further investigate ancient hybridization, we used Supplementary Material online). The placement of Dsuite (Malinsky et al. 2021) to calculate D-statistics and Erythrocebus as a sister to Chlorocebus received low sup­ f4-ratios (Durand et al. 2011) on autosomal data for all port (lpp = 0.53), likely due to rapid radiations and a high possible trios of guenon species, using the rhesus macaque degree of ILS. as outgroup. Out of 1,540 tested trios, 865 produced sig­ Samples from male individuals were available for 17 spe­ nificant D-statistics, indicative of gene flow. However, cies and they were used to construct a Y-chromosome these are not independent as shared ancestry may create phylogeny with RaxML (Stamatakis 2014). Compared signals of gene flow in close relatives or descendants of with the autosomal phylogeny, this tree disagreed only the involved lineages. To disentangle this complex pattern, in the placement of Erythrocebus, in line with the incom­ we used different combinations of D-statistics similar to plete support obtained for this node in other datasets the DFOIL-method (Pease and Hahn 2015). Additionally, (supplementary fig. S6, Supplementary Material online). we applied an approach adapted from the partitioned Considering the high topological concordance among nu­ D-statistic (Eaton and Ree 2013), where we excluded clear markers of different inheritance modes, we hereafter shared variants to detect “carryover effects,” a phenom­ refer to the autosomal phylogeny as the “species tree” enon that produces positive D-values in close relatives (Fig. 1C), although some uncertainty remains regarding to the actual source, as a consequence of their shared an­ the exact placement of Erythrocebus. cestry. Briefly, if directional gene flow occurs from species A to species B, this will also generate a signal of excess allele sharing between any species more closely related to A than Mitonuclear Discordances Suggest Ancient to B, whereas this is not the case if gene flow is restricted to Hybridization the opposite direction (more details in Materials and We assembled the mitochondrial genomes of all samples, Methods; Fig. 2). Identifying these effects allowed us to constructed a maximum likelihood tree in RaxML, and es­ disentangle the primary source and recipient lineages timated divergence times with MCMCTree (Fig. 1C; (i.e. infer directionality) for most gene flow events. supplementary figs. S7 and S8, Supplementary Material on­ We detected the strongest signal of excess allele sharing, line; Materials and Methods). The mitochondrial phyl­ indicative of ancestral hybridization, between the cephus ogeny was consistent with previous reports (Guschanski and mona species groups (D = 0.15 to 0.29, Fig. 2, events et al. 2013) and showed several hard incongruences A1 and A2 in Fig. 3). The f4-ratios, which provide an with the species tree (Fig. 1C; supplementary fig. S7, estimate of the proportion of introgression across the gen­ Supplementary Material online). Most strikingly, the ter­ ome, suggested that at least 10% to 20% of the genome restrial genus Erythrocebus is nested within the Arboreal was introgressed (simulations have shown that f4-ratios clade, sister to the cephus and mitis species groups, and may substantially underestimate the actual proportion 2 terrestrial Allochrocebus species—Allochrocebus lhoesti of gene flow [Svardal et al. 2020]). This gene flow event and Allochrocebus preussi—cluster together with the ar­ also produced significant values of D between the cephus boreal cephus species group, whereas their congener group and Cercopithecus neglectus, the latter being sister Allochrocebus solatus is sister to the terrestrial genus to the mona group (Fig. 2). However, this was likely driven Chlorocebus. by shared ancestry between C. neglectus and the mona Topological differences between nuclear and mitochon­ group, as these D-values approached zero when only drial phylogenies can be caused by ILS or hybridization. We private alleles were considered (Fig. 2). This suggests direc­ explored the probability of ILS causing the observed mito­ tional gene flow from the mona group, with carryover- nuclear discordances using simulations and found that effects from C. neglectus, into the cephus group. topological conflicts where a lineage changes position Similarly, we found excess allele sharing between the from being nested in one clade to being nested within an­ mona and mitis groups, also primarily driven by carryover- other were highly unlikely to occur under ILS alone effects due to the shared ancestry between the cephus and 4 Complex Primate Evolutionary History With Extensive Ancestral Gene Flow · https://doi.org/10.1093/molbev/msad247 MBE A B C. neglectus P1: C. neglecus P1: mitis group P1: mitis group P2: mona group P2: mona group P2: cephus group P2: cephus group P3: cephus group P3: mitis group P3: mona group P3: C. neglectus 0.4 P2: C. pogonias 0.3 P2: C. petaurista Excess P2: C. pogonias allele sharing 0.2 P2-P3 Dstatistic P2: C. petaurista 0.1 P2: C. mona/ lowei C. ne M. m mitis ceph mon 0.0 No Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 P2: C. mona excess ulatta us gr a gro grou glect /lowei P2: C. cephus/ P2: C. cephus/ allele ascanius ascanius sharing -0.1 p us oup up Negative D due to the "mirror effect" of partitioned D-statistics (not gene-flow P3-P1) Fig. 2. D-statistics and partitioned D-statistics. (A) Schematic overview of gene flow event A (Fig. 3), with branch colors corresponding to the partitions used to estimate private allele sharing. Solid black arrows show the inferred gene flow, and dashed gray lines show allele sharing inferred as “carryover” effects. (B) Four tests for excess allele sharing between the mona group, cephus group, mitis group, and C. neglectus, similar to the DFI, DOL, DIL, and DFO of the DFOIL test (Pease and Hahn 2015), respectively. Gray points and distributions show the D-statistics for all combinations of taxa from the re­ spective groups as depicted in the panel header, with M. mulatta as outgroup, using all SNPs. Colored points and violin distributions show the D-statistics for the same trios (connected by gray lines) after removing sites with shared alleles between the cephus and mitis groups (orange, two leftmost panels) or Cercopithecus neglecus and the mona group (blue, two rightmost panels). Note that some D values turn negative after removing shared ancestral variation, which is an expected consequence of this method and should not be interpreted as excess allele sharing between P1 and P3 (Pease and Hahn 2015). Error bars show the lowest and highest D-statistic ±3 standard errors. Positive D-statistics that approaches zero or become negative after shared alleles are removed are interpreted as carryover effects, indicating directionality. mitis groups (Fig. 2). Hence, we conclude that there has B1; supplementary fig. S11, Supplementary Material online). been bidirectional gene flow between the ancestors of However, the terrestrial genus Allochrocebus shares more al­ the cephus and mona groups. The D-values varied substan­ leles with the Arboreal clade than its sister genus tially among cephus and mona group species (Fig. 2), sug­ Chlorocebus (supplementary fig. S12A, Supplementary gesting more than a single ancestral pulse of gene flow. The Material online), suggesting that gene flow continued simplest scenario that could explain the observed D-values or reoccurred along the ancestral Allochrocebus branch involves 2 bidirectional pulses of gene flow: one between (event B2; Fig. 3). This is also supported by the presence Cercopithecus petaurista and the mona group ancestor, of a cephus group-like mtDNA in A. lhoesti/preussi and a second between Cercopithecus pogonias and the (Fig. 1C; discussed below) and a trend towards excess al­ ancestor of Cercopithecus cephus/ascanius (A1 and A2 in lele sharing between specifically the cephus group and Fig. 3). Notably, the split between C. cephus and C. ascanius Allochrocebus (supplementary fig. S12B, Supplementary pre-dates that between C. pogonias and Cercopithecus Material online). The reason that this trend is largely non­ mona/lowei in our divergence time estimates, contradicting significant is likely due to the extensive gene flow within the proposed second gene flow pulse into the C. cephus/ the Arboreal clade (Fig. 3). ascanius ancestor. However, divergence time estimates Furthermore, we found evidence of gene flow between from genomic regions with low introgression between the the ancestor of the mitis group and Cercopithecus hamlyni cephus and mona groups produced a substantially (Fig. 3; event C; D ∼ 0.11, f4-ratio = 0.04), the cephus and younger split time between C. cephus and C. ascanius, post­ mitis groups (event D, 2 pulses of gene flow; D ∼ 0.07, dating the divergence between C. pogonias and C. mona/ f4-ratio = 0.02 to 0.04; Fig. 3; supplementary fig. S13, lowei by ca. 1 MY (supplementary fig. S10, Supplementary Supplementary Material online), C. pogonias and the eastern Material online). These estimates are broadly consistent C. mona population (event E; D = 0.04, f4-ratio = 0.01; also with the divergence times proposed by Kuderna et al. reported by Ayoola et al. [2021], Chlorocebus dryas and (2023) and suggest that introgression from the mona group Chlorocebus pygerythrus (event F; D = 0.03, f4-ratio = 0.01; likely inflated our divergence time estimates among cephus as reported by van der Valk et al. [2020]) and A. solatus group species. and A. lhoesti (event G; D = 0.02, f4-ratio = 0.01). Using The second strongest gene flow signal was detected be­ the polarization approach, we could infer that event D tween the Arboreal and Terrestrial clades (D ∼0.12, f4-ratio was bidirectional. The most parsimonious explanation for ∼0.03 to 0.04; event B1 and B2 in Fig. 3). We applied the the increasing D-statistics with more distant P1-lineages ob­ same polarization approach and found that most D-values served in event C and E is directionality from the mitis group can be attributed to a bidirectional gene flow event be­ into C. hamlyni and from the eastern C. mona populations tween the ancestors of the 2 clades, ca. 11 to 12 MYA (event into C. pogonias, respectively (Fig. 3). We cannot test for the 5 Jensen et al. · https://doi.org/10.1093/molbev/msad247 MBE A B Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 Fig. 3. Excess allele sharing and gene flow events among guenons. (A) D-statistic (circles) and f4-ratios (diamonds) as estimated by Dsuite illustrating excess of allele sharing caused by the 7 identified gene flow events. The combinations of tested taxa are shown on the left (with M. mulatta as outgroup). Positive D-statistic values indicate gene flow between P3 and P2 taxa and are connected by lines to the f4-ratios that reflect the proportion of gene flow. The values were obtained for all possible members of the respective species groups, for which more than a single representative was available in our dataset. In events C and E, we observed decreasing D-statistics with more closely related P1 taxa (annotated as dashed shapes), indicating directionality from P2 into P3. (B) Schematic overview of the identified gene flow events. Arrows are shown where directionality could be inferred. Mitochondrial introgressions, accompanying nuclear gene flow between internal branches, are shown as dashed circles with an arrow. Dashed circles on tips highlight species that carry an introgressed mitochondrion, with circle color corresponding to the most likely introgression event. The donor in B2 was not fully resolved, as illustrated by the question mark, and the dashed line of the C. cephus/ascanius divergence illustrates the likely overestimate of this node's age (see main text). Where available, annotations to the right of species labels show the chromosome number of extant species. (C) Schematic overview of the minimum time since divergence at the time of introgression, based on the divergence time estimates. The divergence time estimate between the eastern and western populations of C. mona was retrieved from Ayoola et al. (2021). opposite direction in these 2 events since our dataset lacks subsets of the tree to test each independent event and found sister lineages to C. hamlyni and C. pogonias. strong support for all of them (supplementary figs. S14 to S21, In total we pinpoint 7 major ancestral gene flow events (A Supplementary Material online). However, we fully acknow­ to G, Fig. 3), some of which occurred in several pulses (A1 and ledge that the presented scenarios are likely an oversimplifi­ A2, B1 and B2, D1 and D2, treated as separate pulses of the cation, and additional gene flow events may have occurred. same event as they are strongly interdependent). The differ­ Among the identified gene flow events, events D and E ence in age between the nodes separating P2 from its closest most likely involved lineages with different chromosome sister and that separating P3 from P2, provides a minimum numbers based on the known karyotypes of the extant spe­ time of divergence at the time of gene flow. These times ran­ cies (Fig. 3; Sineo et al. 1986; Lo Bianco et al. 2017). To inves­ ged from 0.7 to over 5 MY (Fig. 3C), suggesting that repro­ tigate if this was the case in older events, we reconstructed ductive isolation may remain incomplete for long periods. the ancestral karyotypes using ChromEvol (Hohna et al. We also inferred hybridization events with an independent, 2016; Freyman and Höhna 2018). The 4 independent runs re­ model-based approach based on phylogenetic networks turned unique ancestral states for internal nodes and con­ using PhyloNetworks (Solís-Lemus and Ané 2016; Solís- verged poorly (low support [lpp ≤ 25]; supplementary figs. Lemus et al. 2017). Since PhyloNetworks does not allow mul­ S22 and S23, Supplementary Material online). However, all tiple reticulations along the same branch, we analyzed reconstructions suggest that the more ancient events A 6 Complex Primate Evolutionary History With Extensive Ancestral Gene Flow · https://doi.org/10.1093/molbev/msad247 MBE (primarily A2), B, and C occurred across different karyotypes. oxidative phosphorylation (OXPHOS), and several reports Thus, 5 out of 7 major gene flow events are most likely exam­ attributed hybrid inviability to malfunctioning mitonuc­ ples of cross-karyotypic gene flow. The remaining 2 events lear complexes (Burton 2022). Therefore, we investigated occurred within the karyotypically stable Allochrocebus and whether mitochondrial introgressions were facilitated by Chlorocebus genera (F and G), although the karyotype of cointrogression of mitonuclear interacting genes. We Chl. dryas is not known and karyotypic differences may in­ tested this in one of the most extreme cases, the introgres­ volve more than differences in chromosome numbers. sion of a cephus group-like mitochondrion into the ances­ tral Allochrocebus lineage, which involves lineages A Combination of Introgression and ILS Best Explains separated by at least 5 MY of independent evolution Mitonuclear Discordances (Fig. 3C). While A. lhoesti and A. preussi are fixed for an in­ The detected gene flow events directly explain 2 of the 5 hard trogressed mitochondrial genome, their sister, A. solatus, mitonuclear discordances (Figs. 1C and 3; events F and E), retains the ancestral mtDNA. To test if alleles in genes Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 whereas the remaining cases are more complex. Our results with known mitonuclear interacting function (N-mt genes, suggest that these discordances are best explained by a com­ supplementary table S3, Supplementary Material online) bination of introgression and ILS: Introgressed mitochondria were differentially sorted alongside the introgressed segregate in the ancestral population and become fixed only mtDNA, we investigated if these genes were more cephus- in a subset of descendant lineages. For example, the hamlyni like in A. lhoesti/preussi compared with genes without such group is a member of the Arboreal clade in the species tree, interactions. However, N-mt and control genes showed but sister to the Terrestrial clade in the mitochondrial phyl­ similar levels of cephus group ancestry in A. lhoesti/preussi, ogeny. Contradicting the direct introgression of mitochon­ and there was no difference in absolute divergence to dria into hamlyni from the Terrestrial clade, we found no A. solatus, suggesting a lack of broad-scale cointrogression signal of nuclear gene flow into hamlyni after it diverged of N-mt genes (supplementary fig. S26, Supplementary from the rest of the Arboreal clade (supplementary fig. S24, Material online). Next, we assessed SNP patterns on a Supplementary Material online). A scenario in which a gene-by-gene basis to identify signatures of fine-scale coin­ Terrestrial-like mtDNA introgressed into the Arboreal clade trogression of N-mt genes (Fig. 4A). Using the same gene ancestor during event B1 and was subsequently retained sets (N-mt and control), we counted cephus group-like al­ only in the hamlyni group is thus most compatible with leles retained only in either A. lhoesti/preussi or A. solatus. our results (Fig. 3B). Similarly, despite indistinguishable levels We used Chlorocebus to polarize the ancestral Terrestrial of nuclear gene flow between all members of the genus clade allele and counted 2 SNP patterns: category Allochrocebus and the Arboreal clade (supplementary fig. 1—grouping A. lhoesti/preussi with the cephus group and S25, Supplementary Material online), only the A. lhoesti/pre­ A. solatus with Chlorocebus (in agreement with the ussi lineage carries the introgressed cephus group-like mtDNA topology) and category 2—grouping A. solatus mtDNA. This haplotype was most likely transferred into with the cephus group and A. lhoesti/preussi with the Allochrocebus ancestor in event B2 but differentially re­ Chlorocebus (opposing mtDNA topology, Fig. 4A). If only tained during subsequent divergence. The Arboreal-like ILS is involved, we expect equal frequencies of both mtDNA of the genus Erythrocebus is also likely the result of categories. However, if a non-neutral process favors the ancestral introgression and subsequent differential sorting, retention of coadapted mitonuclear variants, we expect as this lineage shows lower levels of nuclear gene flow with an excess of category 1 SNPs. In line with cointrogression the Arboreal clade than other members of the Terrestrial of N-mt alleles, there was a clear excess of category 1 clade (supplementary fig. S11, Supplementary Material on­ SNPs in N-mt compared with control genes (Fig. 4B; line). A complicating factor is that this scenario requires supplementary table S4, Supplementary Material online). the retention of mitochondrial polymorphism over 2 speci­ In N-mt genes, we found 196 category 1 SNPs and 96 ation events (nodes), according to our species tree. category 2 SNPs, whereas similar numbers of both SNP However, as mentioned above, uncertainty remains regard­ categories were found among control genes (a maximum ing the phylogenetic position of Erythrocebus, as alternative of 133 and 131, respectively). More than half of the topologies place it as sister to a clade containing the category 1 N-mt SNPs were located in the genes Chlorocebus and Allochrocebus genera (Kuderna et al. NDUFA10 and LRPPRC (70 and 31, respectively) and 2023). Given this latter placement, the mtDNA polymorph­ included several nonsynonymous variants (3 and 1, ism only needs to be retained over a single speciation event. respectively, supplementary table S3, Supplementary Material online). No nonsynonymous cephus-like variants were found in A. solatus in these genes, and across all con­ Mitochondrial Introgression Across Deeply Divergent sidered N-mt genes, we detected a greater number of Lineages was Facilitated by the Co-sorting of Alleles in cephus-like nonsynonymous variants in lhoesti/preussi Nuclear Genes With Mitonuclear Interacting than in solatus (11 vs. 1). Mutations in NDUFA10 and Functions LRPPRC have been linked to dysfunction of mitochondrial We identified several mitochondrial introgressions across complexes I (Hoefs et al. 2011) and IV (Oláhová et al. 2015), deeply divergent lineages. Nuclear and mitochondrial both of which are involved in OXPHOS. Our results thus genes interact during the essential production of ATP via suggest that retaining cephus-like variants in these genes 7 Jensen et al. · https://doi.org/10.1093/molbev/msad247 MBE A Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 B C Fig. 4. Prevalence of privately retained introgressed cephus group alleles in A. lhoesti/preussi and A. solatus in nuclear genes involved in mito­ chondrial functions (N-mt) compared with the genomic background. (A) Schematic overview of the 2 considered SNP categories. The dashed line in the species tree illustrates the mitochondrial introgression from the cephus group to the ancestral Allochrocebus branch, subsequently retained only in the A. lhoesti/preussi lineage. The circles illustrate the 2 categories of SNPs differentially fixed between the cephus group (light circles) and Chlorocebus (black circles): Category 1 groups A. lhoesti/preussi with cephus and can arise through ILS and mitonuclear cosegrega­ tion, whereas category 2 groups A. solatus with cephus and is expected only from ILS. (B) Number of SNPs per 1,000 base pairs of categories 1 and category 2 across 199 N-mt genes (red diamonds) and 100 samples of 199 other nuclear genes (gray dots). (C) Proportion of category 1 SNPs contributed by 33 N-mt genes with at least one such site, highlighting the identity of the 2 most important genes. may have been essential to overcome mitonuclear incom­ the genome and assigned them to 3 main topologies: patibilities and to facilitate mitochondrial introgression. Tree 1 corresponds to the species tree, Tree 2 shows the introgressed topology, and Tree 3 is consistent with ILS. The species tree—the most frequent topology in all in­ Admixed Genomes are a Mosaic of Different stances—was interspersed with short regions of intro­ Ancestries With Elevated Heterozygosity in gressed and ILS-derived ancestry in all events (Fig. 5A; Introgressed Regions supplementary figs. S27 to S32, Supplementary Material Next, we investigated the spatial patterns of gene flow to online). Consecutive windows of introgressed ancestry understand the genomic architecture of introgression. were rare, in line with ancient timing of gene flow For simplicity, we only tested the trio setup ((P1, P2) P3) (supplementary fig. S33, Supplementary Material online). with the highest D-values for events occurring in several Notably, Tree 2 may also arise through ILS and the propor­ pulses (A1, B1, D1, Fig. 3). We constructed neighbor- tion of introgressed ancestry must be assessed relative to joining trees in non-overlapping 25 kb windows along Tree 3, which is expected to result only from ILS. Tree 2 8 Complex Primate Evolutionary History With Extensive Ancestral Gene Flow · https://doi.org/10.1093/molbev/msad247 MBE A 1 C. neglectus 2 C. petaurista 3 C. neglectus C. neglectus mona group C. petaurista mona group (1) mona group mona group (2) C. petaurista C. neglectus mona group C. petaurista M. mulatta M. mulatta M. mulatta M. mulatta B chr1 chr2 chr3 chr4 Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 chr5 chr6 chr7 chr8 chr9 chr10 2.82 % chr11 C 2.88 % 3.16 % 4.76 % chr12 X-chromosome 9.57 % 9.77 % chr13 Autosomes 18.72 % chr14 33.13 % chr15 chr16 chr17 65.54 % chr18 49.65 % chr19 chr20 chrX 0 50,000,000 100,000,000 150,000,000 200,000,000 Genomic position (Basepair) D Allenopitecus Allochrocebus Erythrocebus Chlorocebus mitis cephus mona neglectus diana hamlyni Miopithecus Papionini 0.006 Mean Heterozygosity 0.004 0.002 Topology 1 2 3 3 11 22 33 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 nigroviridis lhoesti patas dryas albogularis ascanius pogonias neglectus diana hamlyni ogouensis torquatus preussi pygerythrus mitis cephus lowei roloway mulatta solatus sabaeus nictitans petaurista mona Fig. 5. Distribution, prevalence and heterozygosity of introgressed genomic segments. (A) Five possible tree topologies obtained from 25 kb genomic windows for the introgression event A1 (Fig. 3), involving C. neglectus, mona, and cephus groups (here represented by C. petaurista), rooted with M. mulatta as outgroup. Tree 1 groups the lineages according to the species tree. Tree 2 groups the mona group species and C. petaurista monophyletically to the exclusion of C. neglectus and can be caused by ancestral introgression or ILS. Tree 3 places the C. petaurista with C. neglectus and is expected to be caused only by ILS. The remaining 2 trees (light blue and gray) represent more complex topologies, which could be caused by, e.g. more recent introgression or ILS. (B) The genomic location of the different tree topologies along the M. mulatta chro­ mosomes, and (C) their relative abundance on the autosomes and the X-chromosome. Black blocks in (B) correspond to regions of the genome with insufficient information for inferences, frequently located around centromeres and telomeres. (D) Heterozygosity in windows correspond­ ing to Tree 1, Tree 2, and Tree 3 topologies (shown as colored circles along the x axis) calculated for all samples and species. Connectors show significant within-sample differences in heterozygosity between tree topologies, as assessed by ANOVAs followed by post hoc Tukey's tests cor­ rected for multiple testing. was always more frequent than Tree 3 and, in line with our local recombination rate, mainly due to reduced linked D-statistics analyses, we found the highest relative propor­ selection (Nachman 2001; Cutter and Payseur 2013). The tions of Tree 2 in events A to D. same mechanism is predicted to create a positive correl­ To assess how genetic diversity relates to local ancestry, ation between introgression and recombination rate: we calculated heterozygosity in the same 25 kb autosomal Introgressed haplotype blocks, which are predominantly windows as above. Heterozygosity correlates positively with deleterious, are efficiently purged in regions of low 9 Jensen et al. · https://doi.org/10.1093/molbev/msad247 MBE recombination, whereas high recombination rate may de­ (n = 188 to 345 per event, supplementary table S5, couple neutral or beneficial alleles from their deleterious Supplementary Material online). Most genes were present neighbors, leading to a higher local introgression rate in a single event (n = 1,412), in line with low repeatability (Figueiró et al. 2017; Veller et al. 2023). Since most of our of introgression. No functional enrichment was detected study species lack recombination maps, we used heterozy­ when considering the gene lists from each event separately gosity as a rough proxy for recombination rate, hypothesiz­ (Fischer's exact test corrected with false discovery rate ing that it would be higher in regions with introgressed [FDR] > 0.05), but the total set of introgressing genes ancestry. Indeed, we observe this pattern in 3 of the 4 stron­ was significantly enriched for functions involved in “regu­ gest gene flow events (A, B, and D), with Tree 2 windows lation of lymphocyte mediated immunity (GO:0002706)” showing higher heterozygosity than Tree 1 and Tree 3 win­ (FDR = 0.0493; supplementary table S6, Supplementary dows (Fig. 5D; supplementary figs. 27D and 29D, Material online), suggesting that sympatric species may Supplementary Material online). Remarkably, this pattern benefit from exchanging alleles adapted to circulating Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 was not restricted to taxa involved in the particular gene pathogens. For example, the gene CCR2 was repeatedly in­ flow event (cephus and mona in the event A; Fig. 5D), argu­ trogressed in 3 separate events (B, C, and D). This gene is ing against higher heterozygosity being driven by segregat­ involved in the immune response to the simian immuno­ ing introgressed alleles, which would only affect the deficiency virus (SIV; Bissa et al. 2023), which is known to recipient lineages. Instead, it is consistent with a globally be a strong selective force in guenons (Svardal et al. 2017). higher recombination rate in regions prone to retaining in­ We also note 11 introgressed genes involved in meiotic trogressed material. We do not find elevated heterozygosity chromosome segregation (supplementary table S5, in introgressed regions in some of the weaker gene flow Supplementary Material online), as they may play a role events (C, E, F, G; supplementary figs. S28D, S30 to S32D, in cross-karyotypic gene flow events. Nine of them intro­ Supplementary Material online). It is possible that the intro­ gressed between lineages that most likely differed in karyo­ gressed regions are too short and too few in these events or types. Among them are SMC2, which is critical for bivalent that the positive correlation between introgression and lo­ chromosome formation during meiosis (Lee et al. 2011), cal recombination rate is not a universal pattern. and the well-studied PRDM9, which is vital for shaping the recombination landscape (Parvanov et al. 2010). Both introgressed in event B, most likely from the Arboreal to Introgressed Regions Show Low Repeatability Across the Terrestrial clade ancestor (supplementary fig. S35, Independent Gene Flow Events and are Enriched for Supplementary Material online). Immune Genes Next, we quantified the introgression proportion in regions along the genome by calculating FD in non- Introgression Deserts Contain Genes Encoding overlapping 25 kb windows (Martin et al. 2015). If genomic Morphological Traits architecture and local recombination rate are the main de­ Regions with strong signals of introgression were significantly terminants of the genomic location of introgressed regions depleted for several gene ontology (GO)-terms related and if synteny is conserved in different guenon species, we to anatomy and development (e.g. cell differentiation expect to repeatedly find the same genomic regions intro­ [GO:0030154], animal organ development [GO:0048513] gressing in separate gene flow events (Langdon et al. 2022). and anatomical structure morphogenesis [GO:0009653], We found weak, albeit significant, correlations in the intro­ FDR < 0.05, supplementary table S6, Supplementary gression landscape across all 21 comparisons (Pearson's R2 Material online). This suggests that genes related to morpho­ range 0.008 to 0.100 P < 0.05, supplementary fig. S34, logical traits may reduce gene flow, maintaining species in­ Supplementary Material online), suggesting low levels of tegrity. To better understand the processes involved in parallelism. reproductive isolation, we also specifically investigated the We investigated the functional relevance of introgres­ gene content of long genomic regions (≥500 kb) devoid of sion for each event by identifying genes overlapping introgression, here termed “introgression deserts” (see genomic regions with exceptionally high signals of intro­ Materials and Methods). Such regions were likely purged gression (FD ≥ 99th percentile). To reduce the risk of false of introgressed variants rapidly after hybridization and may positives, we also required that P2 showed higher sequence thus contain genes acting as barriers to gene flow. We found similarity with P3 than with P1, and excluded windows 59 to 135 such regions across all introgression events, con­ where sequence divergence between P1 and P3 was below taining 3,354 protein-coding genes (n = 391 to 995 per the 10th percentile (as FD tend to be biased upwards in re­ event, supplementary table S7, Supplementary Material on­ gions of low divergence [Martin et al. 2015]). In principle, line), without any significant GO enrichments. In events A, B, differential sorting of ancestral polymorphisms could gen­ and G, the proportion of the X-chromosome contained in erate similar patterns (Guerrero and Hahn 2017), but the introgression deserts was substantially larger than for any regions would be on average shorter due to their old age autosome, suggesting that the X-chromosome plays an im­ and less likely to pass our outlier criteria. Hence, we portant role in reproductive isolation in these events call these windows putatively introgressed. We found (supplementary fig. S36, Supplementary Material online). 1,588 genes in these putatively introgressed regions In the remaining events, the introgression deserts were 10 Complex Primate Evolutionary History With Extensive Ancestral Gene Flow · https://doi.org/10.1093/molbev/msad247 MBE approximately equally common across the autosomes and Petrogale spp. (Potter et al. 2015). Furthermore, cross- the X-chromosome, spanning ca. 0% to 5% of the chromo­ karyotypic hybridization may generate novel, intermediate somal length. The majority of genes located in introgression karyotypes (Lukhtanov et al. 2020). Considering their ex­ deserts were unique to one event (n =2,647). However, 2 treme karyotypic diversity and highly reticulate evolution, genomic regions were identified as introgression deserts in guenons provide an excellent system to better understand 5 independent events. These contained 10 genes, including the influence of hybridization on chromosomal reorgan­ HPS6 which is involved in melanogenesis and has been asso­ ization and, in turn, the role of such events in reproductive ciated with pelage coloration in the African wild dog Lycaon isolation and speciation. We found several genes with pictus (Chavez et al. 2019). Guenons have distinct face col­ reported functions in chromosomal segregation that orations, suggested to be important for species recognition show strong signatures of introgression across lineages and mate choice (Allen et al. 2014). It is thus likely that genes with different karyotypes (Fig. 3; supplementary table S5, with functions related to pelage coloration are important for Supplementary Material online). Introgression of such Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 maintaining species identity in sympatry. Among genes re­ genes has been suggested to be an important mechanism peatedly found in introgression deserts are also several genes to retain stable meiotic chromosome segregation, follow­ involved in reproductive traits. For example, APPBP2 and ing whole-genome duplication in plants (Marburger et USP32 (events A, B, and C) are associated with reproductive al. 2019; Seear et al. 2020). Hence, they may be involved seasonality in sheep (Martinez-Royo et al. 2017), a trait that in facilitating cross-karyotypic gene flow by promoting is also highly variable in guenons (Enstam and Isbell 2007). pairing and recombination across homologous chromo­ somes despite karyotypic differences. Particularly intri­ guing is the ancestral introgression of PRDM9, which is Discussion known for shaping the recombination landscape in mam­ Using a comprehensive whole-genome dataset, we reveal a mals (Parvanov et al. 2010) and rapidly causing sterility in complex evolutionary history of guenons that was shaped mouse hybrids (Forejt et al. 2021). Homogenization of by multiple ancient hybridization events despite the pres­ PRDM9 during hybridization between the Arboreal and ence of ecological, morphological and karyotypic differ­ the Terrestrial clades may have promoted hybrid fertility, ences. In addition, we present a well-resolved guenon thus facilitating gene flow. species tree that provides a foundation for future studies Mitochondrial introgression has been found in many of guenon evolution. study systems, including primates (Zinner et al. 2018; Supporting the idea that ancestral introgression may Bailey and Stevison 2021), and is more likely to occur spark subsequent radiations through the introduction of across recently diverged lineages as mitonuclear incompat­ novel genetic variation (Han et al. 2017; Meier et al. ibilities accumulate rapidly and can act as early contribu­ 2017; Grant and Grant 2019; Svardal et al. 2020), we find tors to reproductive isolation (Tobler et al. 2019; Burton that several species-rich guenon lineages are descendants 2022). In line with previous studies (Guschanski et al. of hybridization events. For example, we detect strong sig­ 2013), we identified multiple cases of ancestral mitochon­ nals of gene flow between the ancestors of the Arboreal drial introgressions in guenons, several of which occurred and Terrestrial clades, both of which subsequently under­ across deeply divergent lineages separated by at least 4 went extensive radiations. However, their nonadmixed to 5 MY of independent evolution. Biogeographical ana­ sister genera, Miopithecus and Allenopithecus, remained lyses suggest that mitochondrial introgression commonly largely monotypic (2 and 1 species, respectively [IUCN occurs from the resident into the dispersing taxon 2022]). Similarly, the ancestors of the diverse mona and ce­ (Toews and Brelsford 2012; Mastrantonio et al. 2016). phus species groups showed the most pronounced ances­ Female guenons tend to be philopatric, whereas males tral gene flow and exchanged up to 20% of their genomes. disperse at sexual maturity (Enstam and Isbell 2007). Similar admixture proportions have been reported as Mitochondrial introgression is thus likely to happen “hybrid speciation” in other studies of primates (Zhang when males at the front of an expanding population hy­ et al. 2023) and may have provided genetic fuel for bridize with resident females of a different species, a pro­ subsequent radiations in these species-rich groups. cess sometimes referred to as “nuclear swamping” Explicitly testing for the contribution of hybridization to (Zinner et al. 2018; Detwiler 2019; Sørensen et al. 2023). diversification requires broader sampling across the pri­ While the nuclear genome remains homogenized across mate phylogeny, but it is possible that rampant ancient the population, mainly through male-mediated gene hybridization events have contributed to the remarkable flow, the introgressed mitochondria can persist at high fre­ diversity of guenons. quency at the expansion front. If such edge populations Our results suggest frequent gene flow across guenon subsequently become isolated, the pattern of differential lineages with different chromosome numbers. While kar­ retention of mitochondria in the presence of equal levels yotypic differences can act as reproductive barriers of nuclear gene flow may arise. Such a scenario fits well (Giménez et al. 2013; Garagna et al. 2014; Augustijnen with current and inferred historical distribution ranges of et al. 2023; Mackintosh et al. 2023), gene flow across differ­ guenon species with introgressed mitochondria. Ancient ent karyotypes has been demonstrated in the common mitochondrial introgression events tend to have happened shrew Sorex araneus (Horn et al. 2012) and rock wallabies in lineages that likely diverged from the source population 11 Jensen et al. · https://doi.org/10.1093/molbev/msad247 MBE through range expansion leading to secondary contact with immunodeficiency virus (HIV). Supporting the importance other species (Guschanski et al. 2013). It is intriguing to of adaptively introgressed loci with immune function, we speculate whether introgression of mitochondrial and coa­ detect CCR2 in multiple introgression events (B, C, and dapted nuclear variation, as we found e.g. in the A. lhoesti/ D). This gene plays an important role in the inflammatory preussi lineage, may have contributed to reproductive isola­ response following SIV infection, mainly through the mo­ tion from the ancestral population. A similar idea was raised bilization of monocytes (Tsou et al. 2007; Bissa et al. 2023). by Moran et al. (2022), who reported a cointrogression of We also find that elevated heterozygosity in intro­ mitochondrial and coadapted nuclear genes in swordtail gressed genomic regions, a predicted outcome given a fish, Xiphophorus spp., and proposed reproductive isolation positive correlation between recombination rate and arising from mitonuclear incompatibility between the re­ introgression, is not a universal pattern. Although the kar­ sulting lineage and its closest sister species. If this is a general yotypic diversification of guenons may have altered their pattern, mitochondrial introgressions, particularly between recombination landscape, heterozygosity is still expected Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 deeply divergent lineages, might be an important mechan­ to be elevated in introgressed regions in the admixed lin­ ism in the early stages of speciation. eage (but not necessarily globally). However, a recent the­ An intriguing question in speciation research is: To what oretical study suggested that the relationship between extent is the outcome of introgression predictable? Studies recombination and introgression rates may be absent or analyzing multiple independent hybridizing populations negative in the early stages of species divergence and between the same pairs of taxa report highly similar ad­ only later become positive (Dagilis and Matute 2023). mixed genome compositions (Runemark et al. 2018; Two out of 4 events lacking this relationship occurred be­ Langdon et al. 2022; Nouhaud et al. 2022). Guenons consti­ tween closely related lineages (F and G, diverged < 2 My tute an excellent model to study the predictability of before hybridization) and may thus not have been di­ introgression between different taxa in a speciation con­ verged enough to develop the positive relationship be­ tinuum, which has not been previously done. In contrast tween recombination and introgression. However, we to the aforementioned studies, our results reveal introgres­ caution that accurately identifying introgressed haplotype sion landscapes primarily unique to each admixed lineage. blocks from phylogenomic data is challenging, particularly This difference could be due to the level of evolutionary when hybridization is ancient and these blocks are short divergence among hybridizing lineages. High repeatability and contain few informative sites. There may also be other and predictability may be present in closely related popu­ factors obscuring these inferences, e.g. large proportions of lations or species, whereas low repeatability may be the ILS. hallmark of independent hybridizations between different, In summary, our study untangles a complex evolution­ deeply divergent lineages. The weak repeatability of the ary history of guenons shaped by extensive ancient intro­ introgression landscapes in guenons suggests that genomic gression. We found gene flow across deeply divergent incompatibilities accumulated in different genomic re­ lineages, suggesting that reproductive isolation in this pri­ gions and genomic barriers to gene flow are predominantly mate clade may remain incomplete for extended periods lineage specific. It is possible that the karyotypic diversifi­ despite ecological, morphological, and karyotype differ­ cation in guenons has altered their recombination land­ ences. We note that while hybridization frequently resulted scape, which is expected to be a critical factor in driving in interspecific transfers of mitochondrial genomes, not a introgression (Edelman et al. 2019; Martin et al. 2019; single case of introgression of the paternally inherited Langdon et al. 2022). As chromosome-level genome as­ Y-chromosome was identified in our dataset. This strongly semblies from nonmodel organisms become increasingly suggests that male hybrids experience lower fitness, in line available, this is an emerging future research topic. with Haldane's rule (Haldane 1922), which predicts that Nevertheless, we find an underrepresentation of genes in­ postzygotic barriers will likely be stronger in the hetero­ volved in morphological structures and development gametic sex. Our study offers insights into the role and among introgressed regions (supplementary table S6, prevalence of ancient gene flow in mammalian radiations Supplementary Material online) and a repeated occur­ and provides an important framework for future studies rence of genes involved in pelage coloration and repro­ of this diverse primate group. ductive characteristics among introgression deserts (supplementary table S7, Supplementary Material online), pointing to their functional importance in maintaining Materials and Methods species boundaries. Supporting the notion of repeatability on the level of Mapping and Variant Calling gene function, we detect an enrichment for immune genes Our dataset consisted of short-read whole-genome sequen­ among introgressed loci. Sympatric species are likely to be cing data from 37 individuals (22 guenon species and 2 out­ exposed to similar pathogens. Viruses in particular are group species) from different sources (van der Valk et al. known to be a strong selective force in guenons (Svardal 2020; Ayoola et al. 2021; Kuderna et al. 2023; see et al. 2017), and exchanging genes adapted to circulating supplementary table S1, Supplementary Material online for de­ pathogens may thus be highly beneficial. Guenons are nat­ tails). We followed the Genome Analysis Toolkit best practices ural hosts to the SIV, which is closely related to the human for mapping and variant calling (Van der Auwera et al. 2013). 12 Complex Primate Evolutionary History With Extensive Ancestral Gene Flow · https://doi.org/10.1093/molbev/msad247 MBE Briefly, we marked adapters with MarkIlluminaAdapters (6.5 to 15 MY), and (iii) the oldest split within (Picard v. 2.23.4, http://broadinstitute.github.io/picard) and Cercopithecini (6.5 to 12.51 MYA). We used the HKY85 mapped reads to the rhesus macaque (M. mulatta) reference substitution model and ran 4 independent runs of genome (Mmul_10, GenBank: GCA_014858485.1) using 2,000,000 iterations, sampled every 100, and discarded bwa-mem (v 0.7.17). Duplicate reads were marked using the first 2,000 samples as burn-in. After confirming that MarkDuplicates (v. 2.23.4), and we used HaplotypeCaller (v. the runs converged at similar age estimates, the 4 runs 4.1.4.1) and GenotypeGVCFs (v. 4.1.1.0) to call genotypes, were combined and analyzed jointly. Convergence and including invariable sites. Next, we excluded insertions and ESS were assessed in R (R Core Team 2017). deletions (indels) and filtered the remaining calls using VariantFiltration (v. 4.1.1.0) with the recommended exclu­ Mitochondrial Genome Assembly and mtDNA Phylogenic sion criteria (QD < 2.0, QUAL < 30, SOR > 3, FS > 60, MQ Tree < 40, MQRankSum < −12.5, ReadPosRankSum < −8.0). The mitochondrial genomes were assembled using the Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 Next, we excluded multiallelic sites and repetitive regions an­ MitoFinder pipeline (Allio et al. 2020). We trimmed adapt­ notated in the repeat-masker track for the Mmul_10 refer­ er sequences from input reads using Trimmomatic (Bolger ence with BCFtools (Li et al. 2009) and set heterozygous et al. 2014) and then ran MitoFinder with the metaspades genotypes with minor allele support to <0.25 to no call. assembler. Next, assemblies were rotated to the same Last, we used BCFtools to remove sites with coverage below starting position using a custom Python script, and aligned and above one standard deviation of the genome-wide aver­ with MAFFT (Katoh et al. 2019). The alignment was visu­ age or with a fraction of missing genotypes higher than 0.1. ally inspected, and assembly errors from incomplete circu­ larization leading to a short, duplicated sequence were corrected. The cleaned assemblies were then annotated Phylogenetic Inferences with MitoFinder, using the annotation of Chlorocebus Autosomal and X-chromosomal Trees sabaeus (NC_008066.1) as a reference, and divided into We used IQTree (v 2.0, Minh et al. 2020) to construct gene 42 blocks as follows: Each of the 13 protein-coding genes trees across autosomes and the X-chromosome, using the were divided into 1st, 2nd, and 3rd codon position result­ GTR model with 1,000 rapid bootstraps. We used a window ing in 39 blocks, the tRNAs were concatenated into a single size of 25 kb and sampled one region every 500 kb, as the block, and the 2 rRNAs were kept as 2 separate blocks. We multispecies coalescence model (see below) assumes inde­ used PartitionFinder2 (Lanfear et al. 2017) to find the op­ pendent gene trees without internal recombination. Sites timal partitioning scheme. First, retaining all samples, we with missing data were filtered out and windows with a fil­ built a maximum likelihood tree with RaxML (Stamatakis tered length below 10 kb were removed. The filtered gene 2014), with the GTRGAMMA substitution model and trees were then used to independently estimate phylogenies 1,000 rapid bootstraps. The obtained topology was then for autosomes and the X-chromosome under the multispe­ pruned to keep one sample per species, and divergence cies coalescence model with ASTRAL (Zhang et al. 2018). times were estimated with MCMCTree, using the likeli­ ASTRAL was run on all samples and, after confirming the hood function (–usedata 1). As this function is more com­ monophyly of all species, we pruned the tree to include a putationally intense than the approximate likelihood, we single sample per species when estimating divergence times. ran 4 independent runs of 1,000,000 iterations, sampling We followed the IUCN taxonomy (IUCN 2022), with the ex­ every 50th and discarded the initial 2,000 samples as ception that we elevated C. (mitis) albogularis to species sta­ burn-in. The 4 runs converged to highly similar age esti­ tus, following several recent publications (Lo Bianco et al. mates, and after merging we obtained good ESS (≥219) 2017; Detwiler 2019; Shao et al. 2023). for 20/23 nodes, and acceptable ESS (≥181) for all nodes. We used MCMCTree (Yang 2007) for divergence time estimates. To make the analyses computationally feasible, Phylogenetic Simulations to Explore the Role of ILS we generated sequence alignments of 20 randomly To explore the probability that ILS could cause the ob­ sampled genomic windows from the filtered ASTRAL served topological conflicts between the nuclear and mito­ input for the autosomal analyses and 10 windows for the chondrial trees, we used msprime (Baumdicker et al. 2022) X-chromosome. Using a custom python script, genotypes to simulate a demography mimicking the guenon radi­ were converted into haploid alignments by selecting a ation. We simulated a single lineage per species group/ random allele at heterozygous sites. The windows were genus, using the divergence date estimates from the pre­ treated as partitions to allow for rate heterogeneity sent study. The ancestral effective population sizes (Ne) and we used the approximate likelihood calculation of guenons were estimated in the approximate range of (usedata = 3). We applied fossil data calibrations with 100,000 to 250,000 by Kuderna et al. (2023), and since uniform prior distributions, using hard minimum and the mitochondrial genome has an Ne of 1/4th of the auto­ soft maximum bounds (2.5% probability of exceeding somal genome, we used a fixed Ne of 50,000, which is likely upper limit), following de Vries and Beck (2023), on 3 to be in the upper range of the true mitochondrial Ne nodes: (i) the split between Cercocebus and Macaca (5.33 for most lineages. We simulated a haploid sequence of to 12.51 MYA), (ii) the split between the tribes Papionini 15,000 bp in 10,000 replicates. Next, we identified all un­ (including Cercocebus and Macaca) and Cercopithecini ique topologies and quantified their relative frequencies 13 Jensen et al. · https://doi.org/10.1093/molbev/msad247 MBE using a custom Python script. Additionally, we measured 2015). The partitioned D-statistics (Eaton and Ree 2013) the node distance between Erythrocebus (one of the also make use of the shared ancestry between the respect­ lineages involved in a deep mitonuclear discordance) ive sister lineages, such that if a positive D-value between D and its sister genus Chlorocebus in all replicates, to quantify and B is solely driven by sites introgressing from C into B the frequency of deep ILS (supplementary fig. S9, (here referred to as “carryover”-effects), there should be Supplementary Material online). no introgression of alleles emerged on the D lineage after the split from C. Based on both these methods, we calcu­ Y-chromosome Phylogenetic Tree lated D-statistics in a 2-step fashion for all significant quar­ We identified males in our study dataset as samples with tets for which a symmetrical 4-taxon topology was mapped read coverage on X/Y at approximately half of available. In step 1, the baseline D-statistics were calculated that across the autosomes. The Y-chromosomal genotypes for the 4 possible combinations (similar to the DFO, DIL, DFI, were then filtered to obtain high-confidence consensus se­ and DOL of DFOIL, except that we only considered sites Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 quences for these samples by excluding sites where any fe­ where P3 had the derived allele, following Svardal et al. male had a called genotype or any male was heterozygous. [2020]). In step 2, we used a custom Python script to ex­ The filtered sites were concatenated into a single align­ clude SNPs with shared alleles between either A and B ment, and RaxML was run using the same parameters as (DFI and DOL) or C and D (DIL and DFO), thus removing for the mitochondrial data. any signal caused by shared ancestry, and then recalcu­ lated the 4 D-statistics with Dsuite. For example, we can consider the case of bidirectional gene flow between the Gene Flow mona and cephus group ancestors (Fig. 2). In addition to D-statistics the positive D-statistics between these lineages (DFI and We used Dsuite (Malinsky et al. 2021) to calculate the DIL) that are the result of direct gene flow, mona will intro­ D-statistics (Patterson et al. 2012) for autosomal biallelic duce some alleles into cephus that it shares with its sister SNPs for all possible trios of species in agreement with C. neglectus, returning a positive DFI, and cephus will, in our inferred species tree (see Results), using the macaque turn, transfer some alleles to mona that it shares with as outgroup. For this analysis, we grouped the samples by the mitis group (positive DOL). The positive D-values result­ species and focused only on interspecific gene flow. The ing from these indirect effects disappear after step 2, as the only exception was C. mona, for which we kept the eastern shared alleles are removed. and western lineages separate, as they were shown to have experienced differential gene flow and mitochondrial Phylogenetic Networks introgression with an unknown lineage (Ayoola et al. We used PhyloNetworks (Solís-Lemus et al. 2017), a model- 2021). A significant D-value indicates gene flow between based approach, to independently test for the gene flow the mid-group (P3) and either of the in-group sisters P1 events inferred from D-statistics. We used IQTree to con­ (negative D) or P2 (positive D). Since a single hybridization struct ML trees in 25 kb windows, separated by 75 kb, to event can cause multiple significant D-values, additional reduce linkage between windows. We applied the GTR steps were taken to identify hybridization events. First, model of evolution and 1,000 rapid bootstraps. Windows we used a parsimonious approach to explain all significant with <5,000 called sites across all samples were discarded. D-values with as few gene flow events as possible: if all spe­ As network approaches quickly become computationally cies of a particular clade produced similar D-values, we in­ unfeasible with increasing taxa, and only first-level net­ terpreted this as gene flow in their common ancestor. works can be inferred with this approach (Solís-Lemus Second, we used an approach combining the methods of and Ané 2016), this method does not allow for testing DFOIL (Pease and Hahn 2015) and the partitioned all the identified events (see Results) in a single run. D-statistic (Eaton and Ree 2013), which uses a symmetrical Therefore, we ran PhyloNetworks on subsections of the 4-taxon topology plus an outgroup, e.g. (((A, B), (C, D)), species tree (subsets of species), both to give an overview Outgroup), and offers a more detailed understanding of of the most pronounced gene flow events and test the the timing and directionality of gene flow. We used the fol­ support for each (supplementary figs. S14 to S21, lowing rationale: If gene flow occurs from C into B in this Supplementary Material online). Ten runs were performed example topology, the shared ancestry between C and D with hmax 0 to 4 to identify the optimal upper limit of hy­ will lead to significant D-values also between D and B, al­ bridization events. Following this, we ran 30 additional beit lower than between C and B (i.e. when using A as runs with the optimal hmax and selected the network P1, B as P2 and either D or C as P3, respectively). with the best likelihood score with a topology concordant However, if gene flow occurred in the opposite direction, with our inferred species tree. from B into C, there would be no excess allele sharing be­ tween B and D, but instead between A and C (due to the shared ancestry between A and B). DFOIL uses 4 different Reconstructing Ancestral Chromosome Numbers D-statistics (DFO, DIL, DFI, DOL), where the directionality To identify which of the detected gene flow events occurred of gene flow will result in various combinations of positive, across karyotypes, we used ChromEvol, a Bayesian approach negative, or nonsignificant D-values (Pease and Hahn implemented in RevBayes (Höhna et al. 2016), to 14 Complex Primate Evolutionary History With Extensive Ancestral Gene Flow · https://doi.org/10.1093/molbev/msad247 MBE reconstruct the ancestral karyotypic states along the gue­ investigate gene content of introgressed regions and to re­ non species tree. All lineages/species with known chromo­ duce the risk of false positives, we required a number of cri­ some numbers were used in this analysis (Sineo et al. 1986; teria to be fulfilled: FD ≥ 99th percentile of the genome-wide Lo Bianco et al. 2017). First, we used equal rates of chromo­ FD, topology classified as category 2 (“introgressed”), DXY P2– somal gains and losses (delta ∼ dnExponential(10)) and ran P3 < DXY P2–P1, and, since FD tend to be higher in regions of 100,000 iterations, discarding 25% as burn-in. Second, since low absolute divergence (Martin et al. 2015), DXY P1–P3 > the karyotypic diversity of guenons has been attributed 10th percentile of genome-wide divergence. Consecutive mainly to fissions (Moulin et al. 2008), we also performed windows were merged, and protein-coding genes overlap­ 2 runs with higher probability of fissions than fusions (delta ping and contained within 50 kb up and downstream such ∼ dnExponential(11 and 12)). regions were identified using the Mmul_10 annotation. The lists of genes were analyzed for GO enrichments using the online tool PANTHER (Ashburner et al. 2000; Mi et al. Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 Characterizing and Comparing the Introgression 2021). We used the human annotation database, included Landscapes all annotated protein-coding genes in the macaque To map the introgression landscape for all inferred gene annotation as a reference list, and restricted the analyses of flow events, we conducted a number of statistical analyses GO-terms to biological processes. in nonoverlapping windows along the reference genome. We also specifically investigated long genomic regions We chose a window size of 25 kb for these analyses, which with low levels of introgression, as this indicates that intro­ we considered adequately short given the ancient timing gressed alleles were purged rapidly after hybridization and of gene flow but long enough to retain informative sites. suggests a potential role in reproductive isolation. We call We restricted the analyses to windows with at least such regions “introgression deserts,” defined as a region of 1,000 called genotypes and constructed neighbor-joining at least 500 kb (≥20 consecutive windows) with FD below trees for each window using PhyML (Guindon et al. the genome-wide average, allowing for a maximum of 5% 2010). For each introgression event, we then used a custom of windows with missing data. Gene content was assessed Python script to classify the topology of each window into and analyzed in the same way as for highly introgressed 1 of 5 categories: (i) “((P1/P2),P3)”; “species tree”, (ii) regions. “((P2/P3),P1);” “ancestral introgression,” (iii) “((P1/P3), P2)”; “ILS control,” (iv) “((P1/P3_1),(P2/P3_2))”; “differential The Effect of Introgression on Species Divergence introgression” or (v) “complicated/other.” Note that the Time Estimates latter 2 trees only occur in cases where P3 or P2 contain To investigate if introgression affected the species diver­ more than one sample (either multiple samples per species gence time estimates, we recalculated divergence times or multiple species per group in case of ancestral introgres­ using sequence data from 20 windows, 25 kb in length, sion events). The relative abundance of these 5 topologies windows, randomly sampled from 3 different sets of win­ and their distribution along the genome was then visua­ dows depending on the levels of local introgression (FD) lized in R to infer the spatial distribution of introgressed between cephus and mona: (i) low introgression, FD < regions. 25th percentile, (ii) random, and (iii) high introgression, We also investigated how local genomic ancestry relates FD > 75th percentile. MCMCTree was used to estimate di­ to heterozygosity (as a rough proxy for recombination, see vergence times with the same settings as for the species Results), by testing if topology category 1, 2, and 3 was as­ tree. sociated with different levels of heterozygosity. Window heterozygosity for each sample was calculated using a cus­ tom Python script, the effect of topology was tested with Coevolution/introgression of Mitonuclear Genes analysis of variance (ANOVA) in R and post hoc Tukey's To investigate if the coevolution of nuclear genes facili­ tests were applied where significant P-values (P < 0.05) tated mitochondrial introgression, we focused on the were retrieved. P-values from the Tukey's test were cor­ mitochondrial introgression event from the cephus group rected for multiple testing using the Benjamini and ancestor into A. lhoesti/preussi. This event stands out in 2 Hochberg method (Benjamini and Hochberg 1995). ways: (i) it occurred across branches separated by at least 5 Next, we used the python scripts ABBABABAwindows.py MY (according to the mitochondrial divergence dates clo­ and popgenWindows.py (https://github.com/simonhmartin/ ser to 9 MY) and (ii) A. solatus, the sister of A. lhoesti/pre­ genomics_general) to calculate FD and nucleotide divergence ussi, retained the ancestral mtDNA genome and thus (DXY), respectively, along the reference genome, for each provides an opportunity for comparative analyses. We com­ event. FD is related to the D-statistic and estimates discordant piled a list of 199 nuclear-encoded genes with reported allele-sharing patterns in a 4-taxon topology but is more ro­ mitochondrial function (Signes and Fernandez-Vizarra bust in short genomic windows where data is sparse 2018; Bailey and Stevison 2021; supplementary table S3, (Martin et al. 2015). To assess whether windows harbored Supplementary Material online). Next, we constructed con­ similar levels of introgression among independent gene trol gene sets by sampling one gene from the remaining flow events, we used R to calculate Pearson's correlations protein-coding genes for each N-mt gene, requiring that of FD in all pairwise combinations of gene flow events. To they have the same gene length ±10%. This process was 15 Jensen et al. · https://doi.org/10.1093/molbev/msad247 MBE repeated 100 times with replacement, retrieving 100 Author Contributions control gene sets of 199 genes each. We calculated FD and Conceptualization: A.J., K.G., Methodology and analyzes: DXY for each gene set with popgenWindows.py and A.J., K.G., F.S., D.d.V., R.B.; Sample acquisition: K.G., ABBABABAwindows.py (https://github.com/simonhmartin/ L.F.K.K., S.K., I.S.C., J.D.K., A.C.K., K.F., J.R., T.M.-B., C.R.; genomics_general). Fd was calculated for both A. lhoesti/pre­ Writing—original draft: A.J., K.G.; Writing—review and ussi and A. solatus as P2, using all members of the cephus editing: all authors. group as P3, Allenopithecus nigroviridis as P1 and M. mulatta as outgroup. DXY was calculated across 3 pairs: A. lhoesti/pre­ Conflict of interest statement. Employees of Illumina, Inc. ussi versus cephus group, A. lhoesti/preussi versus A. solatus, are indicated in the list of author affiliations. All other and A solatus versus cephus group. If N-mt genes coevolved authors declare no conflict of interest. with the introgressed mtDNA genome, we expect to find a stronger signal of cephus group introgression in these genes Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 compared with control genes in A. lhoesti/preussi, but no Data Availability such pattern in A. solatus. This would also lead to lower The sequencing data used in this project are available on DXY between A. lhoesti/preussi and the cephus group but the European nucleotide archive (https://www.ebi.ac.uk/ increased DXY between A. lhoesti/preussi and their sister ena), under accession numbers as listed in supplementary A. solatus in the N-mt genes. table S1, Supplementary Material online. Custom scripts We also counted private cephus group SNPs in A. lhoesti/ used for data analyses are available at (https://github. preussi in relation to A. solatus. This was done to assess com/axeljen/guenon_phylogenomics). fine-scale variation that may have facilitated the mitochon­ drial introgression. A custom python script was used to identify and count SNPs in 2 categories: Category 1: References A. lhoesti/preussi fixed for the same allele as the cephus Allen WL, Stevens M, Higham JP. Character displacement of group, whereas A. solatus and Chlorocebus spp. are fixed Cercopithecini primate visual signals. Nat Commun. 2014:5(1): for the alternative allele (a nuclear pattern that is concord­ 4266. https://doi.org/10.1038/ncomms5266. ant with mtDNA genome ancestry), Category 2: A. solatus Allio R, Schomaker-Bastos A, Romiguier J, Prosdocimi F, Nabholz B, fixed for the cephus group allele, A. lhoesti/preussi and Delsuc F. MitoFinder: efficient automated large-scale extraction Chlorocebus spp. fixed for the alternative allele (discordant of mitogenomic data in target enrichment phylogenomics. Mol Ecol Resour. 2020:20(4):892–905. https://doi.org/10.1111/1755- with mtDNA genome ancestry). The number of SNPs in 0998.13160. each category was then compared between N-mt and con­ Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, trol genes, and the predicted impact of these mutations was Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: estimated using the Ensembl Variant Effect Predictor tool for the unification of biology. Nat Genet. 2000:25(1): (McLaren et al. 2016). 25–29. https://doi.org/10.1038/75556. Augustijnen H, Baetscher L, Cesanek M, Chkhartishvili T, Dinca V, Iankoshvili G, Ogawa K, Vila R, Klopfstein S, de Vos J, et al. Supplementary Material 2023. A macroevolutionary role for chromosomal fusion and fis­ sion in Erebia butterflies. bioRxiv 2023.01.16.524200. https://doi. Supplementary material is available at Molecular Biology org/10.1101/2023.01.16.524200, 18 January 2023, preprint: not and Evolution online. peer reviewed. Ayoola AO, Zhang B-L, Meisel RP, Nneji LM, Shao Y, Morenikeji OB, Adeola AC, Ng’ang’a SI, Ogunjemite BG, Okeyoyin AO, et al. Acknowledgments Population genomics reveals incipient speciation, introgression, and adaptation in the African Mona Monkey (Cercopithecus We thank Christophe Escudé, Laurianne Cacheux, and mona). Mol Biol Evol. 2021:38(3):876–890. https://doi.org/10. Bertrand Bed’Homme at the Muséum National 1093/molbev/msaa248. d’Histoire Naturelle, Paris, for providing guenon cell cul­ Bailey NP, Stevison LS. Mitonuclear conflict in a macaque species ex­ hibiting phylogenomic discordance. J Evol Biol. 2021:34(10): ture samples, and Jean-Pierre Gautier for tissue samples 1568–1579. https://doi.org/10.1111/jeb.13914. that form the bulk of data used in this project, as well Baumdicker F, Bisschop G, Goldstein D, Gower G, Ragsdale AP, as Mareike Janiak, Tom van der Valk, Simon Martin, Tsambos G, Zhu S, Eldon B, Ellerman EC, Galloway JG, et al. and Konrad Lohse for bioinformatic support and helpful Efficient ancestry and mutation simulation with msprime 1.0. discussions. The computations were enabled by re­ Genetics. 2022:220(3):iyab229. https://doi.org/10.1093/genetics/ iyab229. sources in projects SNIC 2022/6-325 and SNIC 2022/ Bejerano G, Pheasant M, Makunin I, Stephen S, Kent WJ, Mattick JS, 5-561, provided by the Swedish National Infrastructure Haussler D. Ultraconserved elements in the human genome. for Computing (SNIC) at Uppsala University Science. 2004:304(5675):1321–1325. https://doi.org/10.1126/ (UPPMAX), partially funded by the Swedish Research science.1098119. Council through grant agreement no. 2018-05973. The Benjamini Y, Hochberg Y. Controlling the false discovery rate: a prac­ project was supported by the Swedish Research tical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. 1995:57(1):289–300. https://doi.org/10.1111/j.2517- Council VR (2020-03398) to K.G., Zoologiska Stiftelse 6161.1995.tb02031.x. grants to A.J., and a UKRI NERC Standard grant (NE/ Bissa M, Kim S, Galli V, Fourati S, Sarkis S, Arakelyan A, de Castro IS, T000341/1) to D.d.V and R.M.D.B. Rahman MA, Fujiwara S, Vaccari M, et al. HIV vaccine candidate 16 Complex Primate Evolutionary History With Extensive Ancestral Gene Flow · https://doi.org/10.1093/molbev/msad247 MBE efficacy in female macaques mediated by cAMP-dependent ef­ Forejt J, Jansa P, Parvanov E. Hybrid sterility genes in mice (Mus mus­ ferocytosis and V2-specific ADCC. Nat Commun. 2023:14(1): culus): a peculiar case of PRDM9 incompatibility. Trends Genet. 575. https://doi.org/10.1038/s41467-023-36109-8. 2021:37(12):1095–1108. https://doi.org/10.1016/j.tig.2021.06.008. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Freyman WA, Höhna S. Cladogenetic and anagenetic models of Illumina sequence data. Bioinformatics. 2014:30(15):2114–2120. chromosome number evolution: a Bayesian model averaging ap­ https://doi.org/10.1093/bioinformatics/btu170. proach. Syst Biol. 2018:67(2):195–215. https://doi.org/10.1093/ Burton RS. The role of mitonuclear incompatibilities in allopatric sysbio/syx065. speciation. Cell Mol Life Sci. 2022:79(2):1–18. https://doi.org/10. Garagna S, Page J, Fernandez-Donoso R, Zuccotti M, Searle JB. The 1007/s00018-021-04059-3. Robertsonian phenomenon in the house mouse: mutation, mei­ Chavez DE, Gronau I, Hains T, Dikow RB, Frandsen PB, Figueiró HV, osis and speciation. Chromosoma. 2014:123(6):529–544. https:// Garcez FS, Tchaicka L, de Paula RC, Rodrigues FHG, et al. doi.org/10.1007/s00412-014-0477-6. Comparative genomics uncovers the evolutionary history, dem­ Giménez MD, White TA, Hauffe HC, Panithanarak T, Searle JB. ography, and molecular adaptations of South American canids. Understanding the basis of diminished gene flow between hy­ Proc Natl Acad Sci U S A. 2022:119(34):e2205986119. https:// bridizing chromosome races of the house mouse. Evolution. Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 doi.org/10.1073/pnas.2205986119. 2013:67(5):1446–1462. https://doi.org/10.1111/evo.12054. Chavez DE, Gronau I, Hains T, Kliver S, Koepfli K-P, Wayne RK. Glenn ME, Cords M, editors. 2002. The guenons: diversity and adap­ Comparative genomics provides new insights into the remarkable tation in African monkeys. Boston (MA): Springer. adaptations of the African wild dog (Lycaon pictus). Sci Rep. Gopalakrishnan S, Sinding MHS, Ramos-Madrigal J, Niemann J, 2019:9(1):8329. https://doi.org/10.1038/s41598-019-44772-5. Samaniego Castruita JA, Vieira FG, Carøe C, Montero MM, Coyne JA, Orr HA. Speciation. Sunderland (MA): Sinauer Associates; Kuderna L, Serres A, et al. Interspecific gene flow shaped the evo­ 2004. lution of the genus Canis. Curr Biol. 2018:28(21):3441–3449.e5. Cutter AD, Payseur BA. Genomic signatures of selection at linked https://doi.org/10.1016/j.cub.2018.08.041. sites: unifying the disparity among species. Nat Rev Genet. Grant PR, Grant BR. Hybridization increases population variation dur­ 2013:14(4):262–274. https://doi.org/10.1038/nrg3425. ing adaptive radiation. Proc Natl Acad Sci U S A. 2019:116(46): Dagilis AJ, Matute DR. The fitness of an introgressing haplotype 23216–23224. https://doi.org/10.1073/pnas.1913534116. changes over the course of divergence and depends on its size Grubb P, Butynski TM, Oates JF, Bearder SK, Disotell TR, Groves CP, and genomic location. PLoS Biol. 2023:21(7):e3002185. https:// Struhsaker TT. Assessment of the diversity of African primates. doi.org/10.1371/journal.pbio.3002185. Int J Primatol. 2003:24(6):1301–1357. https://doi.org/10.1023/B: de Jong YAD, Butynski TM. Three Sykes's monkey Cercopithecus mitis IJOP.0000005994.86792.b9. × vervet monkey Chlorocebus pygerythrus hybrids in Kenya. Guerrero RF, Hahn MW. Speciation as a sieve for ancestral poly­ Primate Conserv. 2010:25(1):43–56. https://doi.org/10.1896/052. morphism. Mol Ecol. 2017:26(20):5362–5368. https://doi.org/10. 025.0109. 1111/mec.14290. Detwiler KM. Mitochondrial DNA analyses of Cercopithecus mon­ Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. keys reveal a localized hybrid origin for C. mitis doggetti in New algorithms and methods to estimate maximum-likelihood Gombe National Park, Tanzania. Int J Primatol. 2019:40(1): phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 28–52. https://doi.org/10.1007/s10764-018-0029-7. 2010:59(3):307–321. https://doi.org/10.1093/sysbio/syq010. Detwiler KM, Burrell AS, Jolly CJ. Conservation implications of hy­ Guschanski K, Krause J, Sawyer S, Valente LM, Bailey S, Finstermeier K, bridization in African cercopithecine monkeys. Int J Primatol. Sabin R, Gilissen E, Sonet G, Nagy ZT, et al. Next-generation mu­ 2005:26(3):661–684. https://doi.org/10.1007/s10764-005-4372-0. seomics disentangles one of the largest primate radiations. Syst de Vries D, Beck R. 2023. Twenty-five well-justified fossil calibrations Biol. 2013:62(4):539–554. https://doi.org/10.1093/sysbio/syt018. for primate divergences. Palaeontol Electron. 26:1–52. https://doi. Haldane JBS. Sex ratio and unisexual sterility in hybrid animals. J org/10.26879/1249. Genet. 1922:12(2):101–109. https://doi.org/10.1007/BF02983075. Durand EY, Patterson N, Reich D, Slatkin M. Testing for ancient ad­ Han F, Lamichhaney S, Rosemary Grant B, Grant PR, Andersson L, mixture between closely related populations. Mol Biol Evol. Webster MT. Gene flow, ancient polymorphism, and ecological 2011:28(8):2239–2252. https://doi.org/10.1093/molbev/msr048. adaptation shape the genomic landscape of divergence among Dutrillaux B, Muleris M, Couturier J. Chromosomal evolution of Darwin's finches. Genome Res. 2017:27(6):1004–1015. https:// Cercopithecinae. In: Gautier-Hion A, Bourliére F, Gautier JP, doi.org/10.1101/gr.212522.116. Kingdon J, editors. A primate radiation: evolutionary biology of Hoefs SJG, Van Spronsen FJ, Lenssen EWH, Nijtmans LG, Rodenburg the African guenons. New York: Cambridge University Press; RJ, Smeitink JAM, Van Den Heuvel LP. NDUFA10 mutations 1988. p. 150–159. cause complex I deficiency in a patient with Leigh disease. Eur J Eaton DAR, Ree RH. Inferring phylogeny and introgression using Hum Genet. 2011:19(3):270–274. https://doi.org/10.1038/ejhg. RADseq data: an example from flowering plants (Pedicularis: 2010.204. Orobanchaceae). Syst Biol. 2013:62(5):689–706. https://doi.org/ Hohna S, Landis MJ, Heath TA, Boussau B, Lartillot N, Moore BR, 10.1093/sysbio/syt032. Huelsenbeck JP, Ronquist F. RevBayes: Bayesian phylogenetic in­ Edelman NB, Frandsen PB, Miyagi M, Clavijo B, Davey J, Dikow RB, ference using graphical models and an interactive model- García-Accinelli G, Van Belleghem SM, Patterson N, Neafsey specification language. Syst Biol. 2016:65(4):726–736. https:// DE, et al. Genomic architecture and introgression shape a butter­ doi.org/10.1093/sysbio/syw021. fly radiation. Science. 2019:366(6465):594–599. https://doi.org/ Horn A, Basset P, Yannic G, Banaszek A, Borodin PM, Bulatova NS, 10.1126/science.aaw2090. Jadwiszczak K, Jones RM, Polyakov AV, Ratkiewicz M, et al. Enstam Karin L., Isbell Lynne A. 2007. The guenons (genus Chromosomal rearrangements do not seem to affect the gene Cercopithecus) and their ailies – behavioral ecology of polyspeci­ flow in hybrid zones between karyotypic races of the common fic associations. In: Campbell Christina J. Mackinnon Katherine shrew (Sorex araneus). Evol Int J Org Evol. 2012:66(3):882–889. C. Panger Melissa Bearder Simon K, editors. Primates in https://doi.org/10.1111/j.1558-5646.2011.01478.x. Perspective. New York: Oxford University Press. p. 252–274. IUCN. 2022. The IUCN red list of threatened species. Version 2022-1. Figueiró HV, Li G, Trindade FJ, Assis J, Pais F, Fernandes G, SantosSHD, https://www.iucnredlist.org. Hughes GM, Komissarov A, Antunes A, et al. Genome-wide sig­ Katoh K, Rozewicki J, Yamada KD. MAFFT online service: multiple se­ natures of complex introgression and adaptive evolution in the quence alignment, interactive sequence choice and visualization. big cats. Sci Adv. 2017:3(7):1–13. https://doi.org/10.1126/sciadv. Brief Bioinform. 2019:20(4):1160–1166. https://doi.org/10.1093/ 1700299. bib/bbx108. 17 Jensen et al. · https://doi.org/10.1093/molbev/msad247 MBE Kuderna LFK, Gao H, Janiak MC, Kuhlwilm M, Orkin JD, Bataillon T, API. Nucleic Acids Res. 2021:49(D1):D394–D403. https://doi. Manu S, Valenzuela A, Bergman J, Rousselle M, et al. A global catalog org/10.1093/nar/gkaa1106. of whole-genome diversity from 233 primate species. Science. Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, Von 2023:380(6648):906–913. https://doi.org/10.1126/science.abn7829. Haeseler A, Lanfear R. IQ-TREE 2: new models and efficient methods Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. for phylogenetic inference in the genomic era. Mol Biol Evol. PartitionFinder 2: new methods for selecting partitioned models 2020:37(5):1530–1534. https://doi.org/10.1093/molbev/msaa015. of evolution for molecular and morphological phylogenetic ana­ Moran BM, Payne CY, Powell DL, Iverson ENK, Banerjee SM, Madero lyses. Mol Biol Evol. 2017:34(3):772–773. https://doi.org/10.1093/ A, Gunn TR, Langdon QK, Liu F, Matney R, et al. 2022. A lethal molbev/msw260. genetic incompatibility between naturally hybridizing species Langdon QK, Powell DL, Kim B, Banerjee SM, Payne C, Dodge TO, in mitochondrial complex I. bioRxiv 2021.07.13.452279. https:// Moran B, Fascinetto-Zago P, Schumer M. Predictability and paral­ doi.org/10.1101/2021.07.13.452279, 24 August 2023, preprint: lelism in the contemporary evolution of hybrid genomes. PLoS not peer reviewed. Genet. 2022:18(1):e1009914. https://doi.org/10.1371/journal.pgen. Moulin S, Gerbault-Seureau M, Dutrillaux B, Richard FA. 1009914. Phylogenomics of African guenons. Chromosome Res. 2008:16(5): Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 Lee J, Ogushi S, Saitou M, Hirano T. Condensins I and II are essential 783–799. https://doi.org/10.1007/s10577-008-1226-6. for construction of bivalent chromosomes in mouse oocytes. Nachman MW. Single nucleotide polymorphisms and recombin­ Mol Biol Cell. 2011:22(18):3465–3477. https://doi.org/10.1091/ ation rate in humans. Trends Genet. 2001:17(9):481–485. mbc.e11-05-0423. https://doi.org/10.1016/S0168-9525(01)02409-X. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Nouhaud P, Martin SH, Portinha B, Sousa VC, Kulmuni J. Rapid and Abecasis G, Durbin R. The sequence alignment/map format and predictable genome evolution across three hybrid ant popula­ SAMtools. Bioinformatics. 2009:25(16):2078–2079. https://doi. tions. PLoS Biol. 2022:20(12):e3001914. https://doi.org/10.1371/ org/10.1093/bioinformatics/btp352. journal.pbio.3001914. Lo Bianco S, Masters JC, Sineo L. The evolution of the Cercopithecini: Oláhová M, Hardy SA, Hall J, Yarham JW, Haack TB, Wilson WC, a (post)modern synthesis. Evol Anthropol. 2017:26(6):336–349. Alston CL, He L, Aznauryan E, Brown RM, et al. LRPPRC muta­ https://doi.org/10.1002/evan.21567. tions cause early-onset multisystem mitochondrial disease out­ Lukhtanov VA, Dincă V, Friberg M, Vila R, Wiklund C. Incomplete side of the French-Canadian population. Brain. 2015:138(12): sterility of chromosomal hybrids: implications for karyotype evo­ 3503–3519. https://doi.org/10.1093/brain/awv291. lution and homoploid hybrid speciation. Front Genet. 2020:11: Parvanov ED, Petkov PM, Paigen K. Prdm9 controls activation of 583827. https://doi.org/10.3389/fgene.2020.583827. mammalian recombination hotspots. Science. 2010:327(5967): Mackintosh A, Vila R, Laetsch DR, Hayward A, Martin SH, Lohse K. 835. https://doi.org/10.1126/science.1181495. Chromosome fissions and fusions act as barriers to gene flow be­ Patterson N, Moorjani P, Luo Y, Mallick S, Rohland N, Zhan Y, tween brenthis fritillary butterflies. Mol Biol Evol. 2023:40(3): Genschoreck T, Webster T, Reich D. Ancient admixture in hu­ msad043. https://doi.org/10.1093/molbev/msad043. man history. Genetics. 2012:192(3):1065–1093. https://doi.org/ Malinsky M, Matschiner M, Svardal H. Dsuite—fast D-statistics and 10.1534/genetics.112.145037. related admixture evidence from VCF files. Mol Ecol Resour. Pease JB, Hahn MW. Detection and polarization of introgression in a 2021:21(2):584–595. https://doi.org/10.1111/1755-0998.13265. five-taxon phylogeny. Syst Biol. 2015:64(4):651–662. https://doi. Marburger S, Monnahan P, Seear PJ, Martin SH, Koch J, Paajanen P, org/10.1093/sysbio/syv023. Bohutínská M, Higgins JD, Schmickl R, Yant L. Interspecific intro­ Potter S, Moritz C, Eldridge MDB. Gene flow despite complex gression mediates adaptation to whole genome duplication. Nat Robertsonian fusions among rock-wallaby (Petrogale) species. Commun. 2019:10(1):5218. https://doi.org/10.1038/s41467-019- Biol Lett. 2015:11(10):20150731. https://doi.org/10.1098/rsbl.2015. 13159-5. 0731. Marques DA, Meier JI, Seehausen O. A combinatorial view on speci­ R Core Team. R: a language and environment for statistical ation and adaptive radiation. Trends Ecol Evol. 2019:34(6): computing. R Foundation for Statistical Computing, Vienna, 531–544. https://doi.org/10.1016/j.tree.2019.02.008. Austria. https://www.R-project.org/2017. Martin SH, Davey JW, Jiggins CD. Evaluating the use of ABBA–BABA Runemark A, Trier CN, Eroukhmanoff F, Hermansen JS, Matschiner statistics to locate introgressed loci. Mol Biol Evol. 2015:32(1): M, Ravinet M, Elgvin TO, Sætre G-P. Variation and constraints 244–257. https://doi.org/10.1093/molbev/msu269. in hybrid genome formation. Nat Ecol Evol. 2018:2(3):549–556. Martin SH, Davey JW, Salazar C, Jiggins CD. Recombination rate vari­ https://doi.org/10.1038/s41559-017-0437-7. ation shapes barriers to introgression across butterfly genomes. Seear PJ, France MG, Gregory CL, Heavens D, Schmickl R, Yant L, Higgins PLoS Biol. 2019:17(2):e2006288. https://doi.org/10.1371/journal. JD. A novel allele of ASY3 is associated with greater meiotic stability pbio.2006288. in autotetraploid Arabidopsis lyrata. PLoS Genet. 2020:16(7): Martinez-Royo A, Alabart JL, Sarto P, Serrano M, Lahoz B, Folch J, Calvo e1008900. https://doi.org/10.1371/journal.pgen.1008900. JH. Genome-wide association studies for reproductive seasonality Shao Y, Zhou L, Li F, Zhao L, Zhang B-L, Shao F, Chen J-W, Chen C-Y, traits in Rasa Aragonesa sheep breed. Theriogenology. 2017:99: Bi X, Zhuang X-L, et al. Phylogenomic analyses provide insights 21–29. https://doi.org/10.1016/j.theriogenology.2017.05.011. into primate evolution. Science. 2023:380(6648):913–924. Mastrantonio V, Porretta D, Urbanelli S, Crasta G, Nascetti G. https://doi.org/10.1126/science.abn6919. Dynamics of mtDNA introgression during species range expan­ Signes A, Fernandez-Vizarra E. Assembly of mammalian oxidative sion: insights from an experimental longitudinal study. Sci Rep. phosphorylation complexes I–V and supercomplexes. Essays 2016:6:300355. https://doi.org/10.1038/srep30355. Biochem. 2018:62(3):255–270. https://doi.org/10.1042/EBC201 McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, Flicek 70098. P, Cunningham F. The Ensembl variant effect predictor. Genome Sineo L, Stanyon R, Chiarelli B. Chromosomes of the Cercopithecus Biol. 2016:17(1):122. https://doi.org/10.1186/s13059-016-0974-4. aethiops species group: C. aethiops (Linnaeus, 1758), C. cynosurus Meier JI, Marques DA, Mwaiko S, Wagner CE, Excoffier L, Seehausen (Scopoli, 1786), C. pygerythrus (Cuvier, 1821), and C. sabaeus O. Ancient hybridization fuels rapid cichlid fish adaptive radia­ (Linnaeus, 1766). Int J Primatol. 1986:7(6):569–582. https://doi. tions. Nat Commun. 2017:8:14363. https://doi.org/10.1038/ org/10.1007/BF02736662. ncomms14363. Solís-Lemus C, Ané C. Inferring phylogenetic networks with max­ Mi H, Ebert D, Muruganujan A, Mills C, Albou LP, Mushayamaha T, imum pseudolikelihood under incomplete lineage sorting. PLoS Thomas PD. PANTHER version 16: a revised family classification, Genet. 2016:12(3):e1005896. https://doi.org/10.1371/journal. tree-based classification tool, enhancer regions and extensive pgen.1005896. 18 Complex Primate Evolutionary History With Extensive Ancestral Gene Flow · https://doi.org/10.1093/molbev/msad247 MBE Solís-Lemus C, Bastide P, Ané C. PhyloNetworks: a package for phylo­ monocyte mobilization from bone marrow and recruitment to genetic networks. Mol Biol Evol. 2017:34(12):3292–3298. https:// inflammatory sites. J Clin Invest. 2007:117(4):902–909. https:// doi.org/10.1093/molbev/msx235. doi.org/10.1172/JCI29919. Sørensen EF, Harris RA, Zhang L, Raveendran M, Kuderna LFK, van der Auwera GA, Carneiro MO, Hartl C, Poplin R, Del Angel G, Walker JA, Storer JM, Kuhlwilm M, Fontsere C, Seshadri L, et al. Levy-Moonshine A, Jordan T, Shakir K, Roazen D, Thibault J, Genome-wide coancestry reveals details of ancient and recent et al. From FastQ data to high confidence variant calls: the male-driven reticulation in baboons. Science. 2023:380(6648): Genome Analysis Toolkit best practices pipeline. Curr Protoc eabn8153. https://doi.org/10.1126/science.abn8153. Bioinformatics. 2013:43(1110):11.10.1–11.10.33. https://doi.org/ Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and 10.1002/0471250953.bi1110s43. post-analysis of large phylogenies. Bioinformatics. 2014:30(9): van der Valk T, Gonda CM, Silegowa H, Almanza S, Sifuentes- 1312–1313. https://doi.org/10.1093/bioinformatics/btu033. Romero I, Hart TB, Hart JA, Detwiler KM, Guschanski K. Svardal H, Jasinska AJ, Apetrei C, Coppola G, Huang Y, Schmitt CA, The genome of the endangered dryas monkey provides new Jacquelin B, Ramensky V, Müller-Trutwin M, Antonio M, et al. insights into the evolutionary history of the vervets. Mol Ancient hybridization and strong adaptation to viruses across Biol Evol. 2020:37(1):183–194. https://doi.org/10.1093/molbev/ Downloaded from https://academic.oup.com/mbe/article/40/12/msad247/7439455 by guest on 20 December 2023 African vervet monkey populations. Nat Genet. 2017:49(12): msz213. 1705–1713. https://doi.org/10.1038/ng.3980. Veller C, Edelman NB, Muralidhar P, Nowak MA. Recombination and Svardal H, Quah FX, Malinsky M, Ngatunga BP, Miska EA, Salzburger W, selection against introgressed DNA. Evolution. 2023:77(4): Genner MJ, Turner GF, Durbin R. Ancestral hybridization facilitated 1131–1144. https://doi.org/10.1093/evolut/qpad021. species diversification in the Lake Malawi cichlid fish adaptive radi­ Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. ation. Mol Biol Evol. 2020:37(4):1100–1113. https://doi.org/10.1093/ Mol Biol Evol. 2007:24(8):1586–1591. https://doi.org/10.1093/ molbev/msz294. molbev/msm088. Taylor SA, Larson EL. Insights from genomes into the evolutionary im­ Zhang B-L, Chen W, Wang Z, Pang W, Luo M-T, Wang S, Shao Y, portance and prevalence of hybridization in nature. Nat Ecol Evol. HeW-Q, Deng Y, Zhou L, et al. Comparative genomics reveals 2019:3(2):170–177. https://doi.org/10.1038/s41559-018-0777-y. the hybrid origin of a macaque group. Sci Adv. 2023:9(22): Tobler M, Barts N, Greenway R. Mitochondria and the origin of spe­ eadd3580. https://doi.org/10.1126/sciadv.add3580. cies: bridging genetic and ecological perspectives on speciation Zhang C, Rabiee M, Sayyari E, Mirarab S. ASTRAL-III: polynomial time processes. Integr Comp Biol. 2019:59(4):900–911. https://doi. species tree reconstruction from partially resolved gene trees. org/10.1093/icb/icz025. BMC Bioinformatics. 2018:19(S6):153. https://doi.org/10.1186/ Toews DPL, Brelsford A. The biogeography of mitochondrial and nu­ s12859-018-2129-y. clear discordance in animals. Mol Ecol. 2012:21(16):3907–3930. Zinner D, Chuma IS, Knauf S, Roos C. Inverted intergeneric introgres­ https://doi.org/10.1111/j.1365-294X.2012.05664.x. sion between critically endangered kipunjis and yellow baboons Tsou C-L, Peters W, Si Y, Slaymaker S, Aslanian AM, Weisberg SP, in two disjunct populations. Biol Lett. 2018:14(1):20170729. Mack M, Charo IF. Critical roles for CCR2 and MCP-3 in https://doi.org/10.1098/rsbl.2017.0729. 19