Genomic analysis of the rhesus macaque (Macaca mulatta) and the cynomolgus macaque (Macaca fascicularis) uncover polygenic signatures of reinforcement speciation

Abstract

Speciation can involve phases of divergent adaptation in allopatry and ecological/reproductive character displacement in sympatry or parapatry. Reproductive character displacement can result as a means of preventing hybridization, a process known as reinforcement speciation. In this study, we use whole-genome sequencing (WGS) of two closely related primate species that have experienced introgression in their history, the rhesus (Macaca mulatta) and cynomolgus (M. fascicularis) macaques, to identify genes exhibiting reproductive character displacement and other patterns consistent with reinforcement speciation. Using windowed scans of various population genetic statistics to identify signatures of reinforcement, we find 184 candidate genes associated with a variety of functions, including an overrepresentation of multiple neurological functions and several genes involved in sexual development and gametogenesis. These results are consistent with a variety of genes acting in a reinforcement process between these species. We also find signatures of introgression of the Y-chromosome that confirm previous studies suggesting male-driven introgression of M. mulatta into M. fascicularis populations. This study uses WGS to find evidence of the process of reinforcement in primates that have medical and conservation relevance.

Keywords: genomics; macaques; reinforcement; speciation.

FIGURE 1

Range map of Macaca mulatta and Macaca fascicularis with samples in this study. Ranges (shaded in yellow for M. mulatta and red for M. fascicularis) were traced from IUCN Red List range maps last evaluated on November 20, 2015 (M. mulatta) and March 7, 2022 (M. fascicularis). Squares represent samples considered allopatric in this study, and circles represent parapatric samples. Numbers indicate approximate geographical origins of samples in Table 1. Blue shapes indicate sequences introduced in this study, and black shapes indicate previously sequenced samples. The sample numbers correspond to the order in which they appear in Table 1. The map was downloaded from https://d‐maps.com/carte.php?num_car=32142&lang=en. The top photo is of an M. mulatta individual from the National Zoological Park in Delhi, India, taken by Rajiv Bajaj. The bottom photo is of an M. fascicularis individual from Tanjung Puting National Park in Indonesia, taken by Dimitry B. Both were available through Unsplash.

FIGURE 2

Conceptual figure of analyses applied in this study. From left to right, macaques represent allopatric Macaca mulatta, parapatric M. mulatta, parapatric Macaca fascicularis, and allopatric M. fascicularis. In a reinforcement scenario, we expect that some regions of the genome should exhibit increased divergence between parapatric populations (D _p ), increased selection in the population receiving introgression (TD _pf ), and reduced admixture in the population receiving introgression (ADMIX). Variable names in the figure represent variables described for Equations 1 and 2, which are given here as well. The image was created using BioRender.

FIGURE 3

Principal component analysis and ADMIXTURE plots for samples in this study. (a) PCA and ADMIXTURE plot for all autosomes combined. (b) PCA and ADMIXTURE for Y‐chromosomes. Population labeling for all plots is as follows: AlloFas, Allopatric M. fascicularis; AlloMul, Allopatric M. mulatta; ParaFas, Parapatric M. fascicularis; ParaMul, Parapatric M. mulatta. In ADMIXTURE plots, the color red represents the proportion of M. mulatta ancestry, whereas the color cyan represents the proportion of M. fascicularis ancestry.

FIGURE 4

gProfiler results for the candidate gene list. Dot plot showing results of category overrepresentation tests of candidate genes compared to total genes annotated in the rheMac10 genome. Each dot represents a GO, KEGG, or HP category, and the y‐axis represents a negative log transform of the adjusted p‐value, such that higher dots are more statistically significant. Significant categories are opaque, and non‐significant categories are transparent. Categories outlined in black with numbers correspond to the neurological and behavioral categories depicted in the table below.

FIGURE 5

Summary statistics of the gene RARA. RARA, one of the candidate genes resulting from our analysis, is shown here as a representation of signatures in all candidate genes. From top to bottom, different panels represent D _XY Difference, F _st Difference (both defined in Equation 1), Tajima's D Difference (defined in Equation 2), and ADMIXTURE values. X‐axis represents coordinates along chromosome 16. Orange dashed lines signify total length of gene as annotated by Ensembl. Black dashed lines signify the window of the gene that was considered a significant candidate window in our analysis. This window has an unusually high D _XY Difference and F _st Difference as well as an unusually low Tajima's D Difference and ADMIXTURE, whereas other windows in this figure may fulfill one or more of these conditions but not all simultaneously.