Two Y‐chromosomal lineages in White‐cheeked macaque (Macaca leucogenys) indicate a possible male‐mediated introgression in ‘sinica’ group macaques

White-cheeked macaque (Macaca leucogenys: henceforth WCM) is a recently discovered macaque that belongs to the ‘sinica’ species group. Ancestral male-mediated introgression is a probable cause of the origin of WCM (M. leucogenys) and Arunachal macaque (M. munzala). Female philopatry and male-biased gene flow caused phylogenetic discordance between maternal and paternal lineages in macaques. Mitochondrial genomes and Y-chromosomal TSPY gene sequences of WCM were obtained from the whole-genome sequencing of two distant populations suggesting the presence of a single maternal lineage but two paternal lineages. Interestingly, one paternal lineage of WCM clustered with a shared lineage of M. assamensis and M. munzala indicating a possible historic male-mediated introgression in the ‘sinica’ group. The other paternal lineage clustered with the sequences of WCM from the type locality in Southeastern Tibet. A similar evolutionary history was also reported in M. munzala that suggests multiple introgression events within ‘sinica’ group macaques. 白颊猕猴（Macaca leucogenys:以下简称WCM）是最近发现的一个新物种，属于猕猴属“sinica”种组。与藏南猕猴（M. munzala）的情况类似，WCM的物种形成原因，可能源于祖先群中由雄性介导的基因渐渗。其自然种群中，雌性白颊猕猴在出生群留守繁殖，而雄性在群间扩散并带来基因流，这使得分别基于线粒体和Y染色体所得的母系和父系的系统发育关系出现不一致的情况。本研究通过全基因组测序，获得两个相距较远的白颊猕猴种群的线粒体基因组和Y染色体TSPY基因序列，并构建谱系关系。结果显示，尽管它们来源于同一个母系谱系，但呈现出2条不同的父系谱系。有趣的是，其中一个父系谱系与熊猴（M. assamensis）和藏南猕猴聚到同一分支，这表明猕猴属“sinica”种组内，不同物种可能在历史上经历了由雄性介导的基因渐渗。此外，另一个父系谱系与西藏东南地区的白颊猕猴种群聚到同一分支。与藏南猕猴中发现的类似结果一致，这些证据表明“sinica”种组在进化过程中发生过多次渐渗事件。 The identified relic populations of White-cheeked macaques from two districts of Arunachal Pradesh may be prioritized for conservation and management. The male-biased introgression provides empirical evidence for seeking drivers of species divergence and managing conservation units. Macaca leucogenys (White-cheeked Macaque; henceforth WCM) is one of the newly discovered primates, which was first found in Modog County, Southeastern Tibet, China (Li et al., 2015a). Afterward, their distributions were also confirmed from the eastern Himalayan highlands of north-eastern India (Ghosh et al., 2022). A previous study confirmed its phylogenetic position as a member of ‘sinica’ group macaque under the family Cercopithecidae (Fan et al., 2017). Morphologically, they are distinct from other species within the ‘sinica’ group, having general dark dorsal pelage, hairy ventral pelage, prominent white side-whiskers, and round face dark facial skin, long and thick hairs on the neck (Li et al., 2015a). M. leucogenys diverged from its closest relatives around 2.5 Mya (transition from late Pliocene to early Pleistocene, characterized as the onset of major Northern Hemisphere glaciation) and is likely to be basal to other members of this group (Fan et al., 2017). However, different gene tree topologies were found between maternal (mitochondrial tRNAGlu–cyt b and D-loop) and paternal (Y-chromosomal TSPY gene) marker sets in M. leucogenys (Fan et al., 2017). A few studies also uncovered phylogenetic incongruence in wide-ranging species of macaques along the geographic breadth (Morales and Melnick, 1998; Tosi et al., 2000). For example, Tosi et al. (2002) found evidence of Y-chromosomal paraphyly among populations of Vietnamese long-tailed macaque (M. fascicularis) with M. mulatta and mitochondrial incomplete lineage sorting in Chinese M. mulatta that shares haplotype with M. fuscata. Burmese long-tailed macaques also showed discordance between maternal and paternal lineage that suggest ancient hybridization (Matsudaira et al., 2018). A recent genomic study detected gene flow within and between species groups (Song et al., 2021; Vanderpool et al., 2020). Whether widespread gene flow among primates is emblematic of their initial radiation (which began 60–75 Mya, Chatterjee et al., 2009; Springer et al., 2012) or it is a consequence of their current distribution overlaps which include higher environmental occupancy and more secondary contact—remains an open question (Tung and Barreiro, 2017). In the present study, we re-sequenced the genomes of two WCMs, that is, one each from central and eastern Arunachal Pradesh, and analyzed whole mitogenome and TSPY gene sequences of Y-chromosome to further investigate the complex evolutionary history and divergence patterns of WCM within ‘sinica’ group macaques. Genomic DNA was extracted from two skins of hunting trophies (possessed by villagers of West Siang district and Dibang district of Arunachal Pradesh: Figure 1) using Qiagen DNeasy Blood and Tissue Kit (Qiagen Paired-end libraries were prepared using Illumina TruSeq PCR-free HT Library Prep Kit (Illumina Inc.) and sequenced on Illumina HiSeq. 2500 platform (2×150 bp). We performed a quality check of the raw reads using FastQC and filtered out low-quality sequences (Q < 30) and short reads (<100 bases), removed homopolymers (>8 bases) and ambiguous bases (N > 0), and trimmed bases from start (up to 10 bases) using fastp (Chen et al., 2018). The qualified reads were then aligned with the reference mitochondrial genome of M. leucogenys (NC_031156.1) (Hou et al., 2016), using the BWA-MEM (Li, 2014) implemented in Galaxy. Mapping quality was visualized by QualiMap BamQC v2.2.1 (Okonechnikov et al., 2016) and the mapping file was filtered to remove sequences with low mapping quality (<50) using BamTools v2.4.2 (Barnett et al., 2011). The filtered mapping (.bam) file was then converted to sequences (.fasta) using Samtools fastx v1.15.1 (Li et al., 2009). Contigs were generated from the filtered mapped sequences using MEGAHIT v1.2.9 (Li et al., 2015b). The same methodology was used to map and extract contigs of TSPY gene sequences using the TSPY gene sequence of reference M. leucogenys (KX066413.1). The reference mitochondrial genomes and TSPY gene sequences of other macaque species were retrieved from NCBI/GenBank and aligned using MAFFT online server (Katoh et al., 2019). Suitable base substitution model identification was carried out based on the Akaike information criterion (AIC) and a Maximum likelihood phylogenetic tree was constructed using RAxML v8 (Stamatakis, 2014) following GTR+G+I substitution model and 1000 bootstrap replicates. The resulting phylogenetic trees were visualized using Interactive Tree Of Life (iTOL) v5 (Letunic & Bork, 2021). Furthermore, to validate the consistency of tree topology, we reconstructed phylogeny using the maximum parsimony method in MEGA-X (Kumar et al., 2018) with the whole mitochondrial genome and TSPY gene sequences. Interestingly, the phylogenetic tree derived from the Y-chromosomal TSPY gene showed two distinct paternal lineages of M. leucogenys-subclade A and B (Figure 2a), while the phylogenetic tree from the mitochondrial genome showed monophyly (Figure 2b). More specifically, WCM sampled from eastern Arunachal Pradesh (this study) clustered with the sequences of M. leucogenys sampled from the type locality (Fan et al., 2017)—subclade A. However, WCM from the West Siang district of central Arunachal Pradesh clustered with the M. munzala and M. assamensis, subclade-B. We observed the monophyletic origin of WCM based on mitogenome sequences from eastern and central Arunachal Pradesh, that clustered with the sequences of M. leucogenys from the type locality (Fan et al., 2017). The observed discordance between the maternally and paternally inherited markers was consistent with the former studies (Chakraborty et al., 2007; Fan et al., 2017; Khanal et al., 2021). Further, a relatively low variance was observed in the TSPY gene in comparison to whole mitochondrial sequences. Positions of the variable sites (n = 63) across the aligned 1336 bp TSPY gene fragments are listed in Supporting Information: Table S1. The occurrence of sub-clade B in the TSPY gene clustered with M. munzala and M. assamensis indicates the inter-specific male-mediated introgression in the past. Similar tree topologies were also observed in maximum parsimony trees (Figure 3a,b) supporting the results obtained with ML trees. Assuming male-mediated dispersal from natal groups and female philopatry in macaque societies (Thierry, 2011), the findings of the present study support the hypothesis of the introgressive origin of M. leucogenys proposed by Fan et al. (2017) similar to M. munzala (Chakraborty et al., 2007). This study also agreed with other previous studies (Fan et al., 2017; Khanal et al., 2021) that all ‘sinica’ group macaques and M. arctoides belong to the same Y-chromosomal clade. The phylogenetic incongruence between maternal and paternal markers was previously reported in M. mulatta, M. fascicularis, and M. arctoides (Fan et al., 2018; Matsudaira et al., 2018; Roos and Zinner, 2015; Rovie-Ryan et al., 2021; Tosi et al., 2002). Interbreeding between species of the same species group and between species from different species groups of macaques has been widely observed in both the wild and in captivity (Bernstein, 1980; Fa, 1989). In this study, M. leucogenys from central Arunachal Pradesh shared a clade with M. assamensis and M. munzala indicating a possible male-mediated hybridization among M. assamensis, M. munzala, and M. leucogenys. We presumed as M. assamensis harbors a relatively large distribution range and population size, possibly M. leucogenys diverged first from other macaque species within the ‘sinica’ group. A similar evolutionary history (shared phylogenetic clade) of paternal lineage observed in M. munzala indicates the potential hybridization event as proposed by Chakraborty et al. (2007). However, further study using a large number of nuclear loci or whole genome is required to understand the complex evolutionary history including multiple introgression events in the sinica group macaque, and this can be studied when the genomes of all species under the genus Macaca are sequenced. Avijit Ghosh: Data curation; formal analysis; investigation; methodology; visualization; writing—original draft. Lalit Kumar Sharma: Conceptualization; funding acquisition; project administration; supervision. Mukesh Thakur: Conceptualization; methodology; project administration; resources; supervision; validation; visualization; writing—original draft; writing—review & editing. The authors thank the State Forest Department of Arunachal Pradesh for granting permission to undertake the field survey and sample collection. The authors also acknowledge the local guides and field assistants who helped during fieldwork. The research is funded by the Department of Science and Technology, New Delhi under INSPIRE FACULTY SCHEME (Grant No. DST/INSPIRE/04/2016/002246) and the NMHS-Large Grant project of the Ministry of Environment, Forest and Climate Change, New Delhi (Grant No. NMHS/2017-18/LG09/02/476). The authors declare no conflict of interest. Field work and sample collection and were carried out after obtaining research permission from the State Forest Department of Arunachal Pradesh (letter no. CWL/Gen/173/2018-18/Pt VII/1728-37/6.09.2018). Avijit Ghosh is a research scholar of the Zoological Survey of India who is enrolled for a PhD at Department of Zoology, University of Calcutta, Kolkata, West Bengal. He has gained about 5 years of research experience in biodiversity assessment. Lalit Kumar Sharma is a senior scientist and in charge of the wildlife section at the Zoological Survey of India, with over 15 years of professional experience in wildlife monitoring and biodiversity assessment. Mukesh Thakur is a senior scientist, in charge of mammal & osteology section at the Zoological Survey of India with over 15 years of professional experience of wildlife monitoring, conservation genetics, and wildlife forensics. He has a special interest in exploring the utility of disregarded genetic diversity in conservation management and policy planning. Sequences of whole mitochondrial genome and TSPY gene generated in this study are available on NCBI/GenBank (accession no. of whole mitochondrial genomes: MT755630 [from Dibang] and MT949692 [from West Siang] and partial TSPY gene fragments: OM909019 [from Dibang] and OM909020 [from West Siang]). Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.