Population genomics of sika deer reveals recent speciation and genetic selective signatures during evolution and domestication

Abstract

Background

Population genomic analysis can reconstruct the phylogenetic relationship and demographic history, and identify genomic selective signatures of a species. To date, fundamental aspects of population genomic analyses, such as intraspecies taxonomy, evolutionary history, and adaptive evolution, of sika deer have not been systematically investigated. Furthermore, accumulating lines of evidences have illustrated that incorrect species delimitation will mislead conservation decisions, and even lead to irreversible mistakes in threatened species.

Results

In this study, we resequenced 81 wild and 71 domesticated sika deer representing 10 main geographic populations and two farms to clarify the species delimitation, demographic and divergence histories, and adaptive evolution of this species. First, our analyses of whole genomes, Y chromosomes and mitochondrial genomes revealed substantial genetic differentiation between the continental and Japanese lineages of sika deer, representing two phylogenetically distinct species. Second, sika deer in Japan were inferred to have experienced a “divergence-mixing-isolation” evolutionary scenario. Third, we identified four candidate genes (XKR4, NPAS3, CTNNA3, and CNTNAP5) possibly involved in body size regulation of sika deer by selective sweep analysis. Furthermore, we also detected two candidate genes (NRP2 and EDIL3) that may be associated with an important economic trait (antler weight) were under selection during the process of domestication.

Conclusion

Population genomic analyses revealed that the continental and Japanese lineages represent distinct phylogenetic species. Moreover, our results provide insights into the genetic selection signatures related to body size differences and a valuable genomic resource for future genetic studies and genomics-informed breeding of sika deer.

Introduction

Population genomic analysis is crucial for reconstructing the phylogenetic relationship and demographic history, and identifying genomic signatures underlying selection of a species [1, 2]. Native to East Asia, sika deer (Cervus nippon) were once widely distributed in the eastern Eurasian mainland (Far East of Russia, China, the Korean Peninsula, and Vietnam), Taiwan island and the Japanese archipelago [3, 4]. However, wild sika deer populations have gradually declined due to human activities such as overexploitation via hunting and habitat destruction caused by urbanization [5, 6], and they are only distributed in nature reserves in some areas. To date, population genomics research on this specie is scarce [7], as it is distributed in different countries and has limited wild populations, making sampling difficult.

Based on morphological characteristics (such as antlers, skulls, body size, and coat color pattern), this species has been taxonomically classified into 13 subspecies: six subspecies in eastern Eurasian continent and Taiwan island, six in Japan, and one in Vietnam [8,9,10]. Of the six subspecies that have inhabited eastern Eurasian continent, two (C. n. grassianus and C. n. mandarinus) are considered extinct today. The remaining four subspecies are (1) C. n. hortulorum, which is distributed in northeastern China, the Far East of Russia, and the Korean Peninsula; (2) C. n. sichuanicus, which lives in the Tiebu Nature Reserve in Sichuan Province; (3) C. n. kopschi, which lives in the Qingliangfeng Nature Reserve in Zhejiang and the Taohongling Nature Reserve in Jiangxi; and (4) C. n. taiouanus, which is extinct in the wild but currently lives in Kenting National Park after restoration from captive-bred wild individuals. In Japan, sika deer are classified into six subspecies [8]: (1) C. n. yesoensis (Hokkaido Island), (2) C. n. centralis (Honshu mainland and Tsushima Island), (3) C. n. nippon (Kyushu, Shikoku, and Goto Islands), (4) C. n. mageshimae (Mageshima and Tanegashima Islands), (5) C. n. yakushimae (Yakushima and Kuchinoerabu Islands), and (6) C. n. keramae (Ryukyu Islands). The last subspecies, C. n. pseudaxis is distributed in Vietnam but is extinct in the wild; it is only raised on deer farms and relocated to parks, reserves, and zoos (pp. 541–548 in [10]).

Fundamental aspects of evolutionary and conservation biology, such as intraspecies taxonomy, evolutionary history, and adaptive evolution, of this species have not been investigated by genome-wide single nucleotide polymorphism (SNP) markers but only by mitochondrial gene fragments [5, 11,12,13,14,15] and Y-chromosome genes [16], which have some limitations [17, 18]. Accumulating lines of evidences have illustrated that incorrect species delimitation will mislead conservation decisions, such as protection and management practices, and even lead to irreversible mistakes in threatened species [19,20,21]. SNPs identified by whole-genome sequencing show great potential not only to assist in species delimitation and phylogeographic evolutionary analyses at an appropriate temporal scale [19] but also to help elucidate microevolutionary mechanisms on a population scale [22], for instance, dissecting genomic signatures of selection and local adaptation [21, 23].

In addition, ecological pressures and environmental conditions could influence the body size evolution of vertebrates living on islands. The classic insular pattern characterized by gigantism in small animals and dwarfism in large animals has been summarized as a macroevolutionary or biogeographic concept, the so-called “island rule” [24,25,26]. However, the microevolutionary mechanism driving patterns of insular dwarfism in large animals remains largely elusive. Since body size is intimately linked to many physiological and ecological characteristics of vertebrates, it is predictable that shifts in body size would leave signals of selection and local adaptation in the genome during the evolution of a species. For instance, a recent study revealed that selective signatures in the fatty acid desaturase (FADS) gene cluster, likely related to diet, are associated with stature size differences between an extant pygmy population on Flores Island and modern humans, which provided a paradigm to bridge the gap between selective signatures in the genome and insular dwarfism in humans [27]. Overall, research on screening genes related to body size differences between island and mainland mammals is still scarce. Furthermore, it is also unknown whether the genes regulating body size related to the “island rule” in humans and other mammals are conserved. Paired island-mainland populations of mammals, birds, reptiles and amphibians are frequently used as models for exploring the patterns of body size evolution in terrestrial vertebrates [26]. In this study, we have been able to sample populations from the eastern Eurasian continent and the Japanese archipelago, which provide an excellent model for exploring the microevolutionary mechanism of this phenomenon.

Furthermore, the sika deer has a history of domestication and breeding for more than 300 years in China [28], and antler has been used as a traditional Chinese medicine for more than 2,000 years [29, 30]. Due to their economic and medicinal value, sika deer are commonly raised in captivity in eastern Asia. During the domestication process, antler weight is an important target of artificial selection, resulting in a 30–60% increase in antler weight in artificially bred populations compared to uncultivated ones [28]. However, the molecular mechanism regulating the weight of velvet antler has not been well explained.

Recently, we assembled the first high-quality chromosome-level reference genome of sika deer, which represents a valuable genetic resource for research on and the conservation of this species [4]. In the present study, we resequenced the genomes of 81 wild and 71 domesticated sika deer representing 10 main geographic populations and two farms with an average depth of 26.7 × . We used these data to clarify the species delimitation, demographic and divergence histories, and adaptive evolution of wild sika deer. Furthermore, we identified potential candidate genes that may underlie the observed body size variations among population pairs from the Japanese archipelago. Finally, we characterized the signatures of selection related to antler size during the domestication of sika deer.

Results

Whole-genome sequencing and genetic variation

We resequenced the genomes of 81 wild sika deer from 10 sampling sites, including a nature reserve (RU1, n = 16) and a wild region (RU2, n = 6) in the Far East of Russia, the Tiebu Nature Reserve in Sichuan Province (SC, n = 8), the Qingliangfeng National Nature Reserve in Zhejiang Province (ZJ, n = 8), the Liugongdao National Forest Park in Shandong Province (TW, n = 2), Hokkaido (HK, n = 9), Nikko in Honshu (NK, n = 8), Shimane in Honshu (SM, n = 2), Tsushima Island (TS, n = 14), and Yakushima Island (YAK, n = 8). The last five sampling sites are located in Japan. In addition, a total of 71 domesticated sika deer from farms in Jilin and Heilongjiang Provinces (Northeast China, NEC) were also collected for genome resequencing (Fig. 1A and Supplementary table S1).

Fig.1

Geographic distributions and population genetic structure of the sika deer. A Sampling locations in this study. A total of 81 wild sika deer representing populations from China, Russia, and Japan and 71 domesticated sika deer from farms in Jilin and Heilongjiang Provinces (Northeast China, NEC) were sampled. Triangles and circles were used to distinguish between mainland and island populations, and two domestic populations were labelled with black edges. The black squares denoted non-native populations: the Kazakhstan population was introduced from a nature reserve in the Far East of Russia (RU1) in the 1930s, whereas the population in the Liugongdao National Forest Park, Shandong province, was introduced from Taiwan on April 16, 2011. B NJ tree constructed using p-distances of all 152 individuals. C ADMIXTURE results obtained using all 152 sika deer with K = 2, 3 and 8. D Principal component analysis plot of all 152 samples used in this study. E Phylogenetic tree based on mitochondrial genome haplotypes. F Median-joining network map based on Y-chromosome SNP haplotypes

A total of 13.22 Tb of high-quality Illumina clean reads were generated (Supplementary table S2). Genome mapping showed that approximately 98.46% of the reads covered 98.68% of the reference genome, with an average genome depth of 26.7 × . After SNP calling and quality filtering, we identified a total of 31,320,151 SNPs. Among these SNPs, 242,811 were located within exonic regions, 9,605,504 in intronic regions, 796 in splicing sites, 430,223 in upstream or downstream regions, and the remaining 21,040,817 in intergenic regions. In the coding regions, we annotated 102,530 nonsynonymous, 218 stop-loss, and 1,799 stop-gain SNPs that cause amino acid changes, elongated transcripts or premature termination (Supplementary table S3). In addition, we also detected 2,034,186 insertions and deletions (InDels) with a minor allele frequency (MAF) > 0.05 (Supplementary table S4). Genetic variants of sika deer along the 33 chromosomes indicated they were evenly distributed across the genome (Supplementary Figure S1 and S2).

Since our genome assembly was from a female individual [4], the Y chromosome of another sika deer genome (GCA_040085125.1) in National Center for Biotechnology Information (NCBI) RefSeq database (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_040085125.1/) was used to call Y-chromosome (Y-chr) SNPs. We identified a total of 543,807 Y-chr SNPs in 129 male individuals. In addition, we assembled the mitochondrial genome of each individual from the genome resequencing reads of 81 wild and 71 domesticated sika deer and detected 914 mitochondrial SNPs.

The average SNP densities were 12.53/kb, 16.67/kb, and 19.72/kb in autosomes, X-chromosomes, and Y-chromosomes, respectively. This result indicated that the nucleotide diversities of nonrecombining sex chromosomes were higher than those of autosomes, which contrasted with the observed pattern in some livestock species where sex chromosomes exhibit reduced nucleotide diversity compared to autosomes due to selective breeding [31, 32]. This distinction was possibly due to that most of our samples were wild (n = 81), and although 71 domesticated samples were also included, the domesticating and breeding history of sika deer was much shorter than other livestock.

In addition, we detected a total of 132,644 SVs, including insertions (INS), deletions (DEL), intrachromosomal translocations (ITX), interchromosomal translocations (CTX), and inversions (INV) (Supplementary tables S5, S6). In total, 48,570 copy number variant (CNV) were also identified (Supplementary table S7).

Genomic evidence of two phylogenetic species of sika deer

Based on the autosomal SNPs, we performed phylogenetic tree construction, population structure analysis, and principal component analysis (PCA) (Fig. 1B-D). In the phylogenetic tree, we found two divergent genetic lineages, the continental lineage and the Japanese lineage, in sika deer. The continental lineage included four clades, northern (NEC, RU1 and RU2), SC, ZJ, and TW, which was consistent with their geographic origin, corresponding well to the four subspecies C. n. hortulorum, C. n. sichuanicus, C. n. kopschi, and C. n. taiouanus as classified in traditional taxonomy [9]. On the other hand, populations from HK, NK, SM, TS, and YAK were classified into five clades, with individuals in each population forming a single clade (Fig. 1B). In Japan, six subspecies have been reported for sika deer, but unfortunately, only three of these subspecies were sampled during our efforts. The five clades together formed the Japanese lineage, which displayed a smaller cluster size with a shorter branch length than the continental lineage, indicating that the genetic differentiation within the Japanese lineage was relatively shallow. Moreover, the species tree reconstructed with SNAPP [33, 34] also demonstrated that the continental lineage and Japanese lineage were reciprocally monophyletic and the divergence events within the Japanese lineage occurred within the last 0.1 million years (Ma) (Supplementary Figure S3 and S4).

The major lineage classification was further validated through ADMIXTURE analysis (Fig. 1C and Supplement Figure S5). Although cross-validation (CV) error minimization indicated K = 3 as the optimal value (Supplement Figure S6), examination of alternative K-values provided additional insights into the admixture history of sika deer populations. At K = 2, a clear genetic divergence emerged between the continental and Japanese lineages. At K = 3, this primary division persisted while additionally resolving a secondary divergence between northern (RU1, RU2, NEC) and southern (SC, ZJ, TW) subclades within the continental lineage. At K = 8, finer-scale population structure was resolved: within the continental lineage, the SC and ZJ populations exhibited relatively homogeneous genetic composition, whereas RU1, RU2, NEC and TW displayed admixture. Similarly, in the Japanese lineage, the HK, YAK and TS populations showed limited admixture, while NK and SM were mixed populations (Fig. 1C and Supplementary Figure S5).

The PCA revealed clear genetic structure, with samples from each geographic region clustering together (Fig. 1D). The first component separated the continental lineage from the Japanese lineage. Within the continental lineage, separation was detected between the northern (including NEC, RU1 and RU2) and southern (SC, TW and ZJ) populations in eastern Eurasian continent and Taiwan island along the axis of the second component. Our PCA results is consistent with the pattern revealed by the phylogenetic tree and ADMIXTURE.

Complementing autosomal SNP evidence, both uniparental inheritance systems – maternally transmitted mitochondrial DNA and paternally inherited Y chromosome – consistently resolved the deep evolutionary split between continental and Japanese sika deer lineages (Fig. 1E-F and Supplementary Figure S7). Phylogenetic reconstruction of mitochondrial haplotypes (39 haplotypes defined by 914 SNPs) revealed reciprocally monophyletic clades (Fig. 1E and Supplementary Figure S7A), while Y chromosome analyses integrating phylogenetic topology (52 haplotypes derived from 211,326 SNPs filtered with TrimAl) and median-joining network patterns (35 haplotypes resolved by 1,248 high-confidence SNPs after removing heterozygous sites and sites with missing genotypes) further confirmed complete lineage partitioning, with no haplotypes shared between the two groups (Fig. 1F and Supplementary Figure S7B).

To further evaluate whether the two lineages represent distinct species, we conducted a coalescent model-based species delimitation analysis using BPP v. 4.3.0 [35]. The results unequivocally supported their delimitation as distinct species, with posterior probabilities (PP) of 1.0 for the two-species model and 0 for the one-species model. A PP of 1.0 indicated that the genomic data decisively rejected lineage coalescence under a shared ancestral population, strongly favoring independent evolutionary histories. Under the phylogenetic species concept (PSC) [36], species are defined as irreducible monophyletic lineages diagnosable by fixed characters or coalescent thresholds. The BPP results was consistent with PSC since the two lineages form reciprocally monophyletic clades in the species tree inference.

Population genetic diversity

To assess the genetic diversity of different populations of sika deer, we calculated the genetic differentiation statistic (F_ST), nucleotide diversity (θπ), observed heterozygosity (Het), linkage disequilibrium (LD), and inbreeding coefficient of each population (Fig. 2). The TW and SM populations were excluded from the genetic diversity analysis since only two samples were collected at both sites; thus, we were unable to provide an accurate assessment for these two populations.

Fig. 2

Population differentiation (F_ST), genetic diversity, linkage disequilibrium, and inbreeding coefficient of different sika deer populations. A Genetic differentiation statistics (F_ST) of different sika deer populations. B, C. nucleotide diversity (θπ) (B) and Observed heterozygosity (C) of different sika deer populations. D Genome-wide average linkage disequilibrium decay estimated from different sika deer populations. E Inbreeding coefficient of the 9 sika deer populations. The inbreeding coefficient is calculated as the fraction of runs of homozygosity (ROH). The TW and SM populations were excluded in these analyses since only two samples were collected at both sites; thus, we were unable to provide an accurate assessment for these two populations

F_ST measures based on whole-genome SNPs showed that F_ST between the continental and Japanese lineages ranged from 0.24 to 0.55, which was much higher than that between populations within the same lineage (0.19–0.43 within the Japanese lineage, 0.02–0.26 within the continental lineage) (Fig. 2A).

The average θπ value of the populations from the continental lineage was 1.62–2.89 × 10^–3 (Fig. 2B), higher than that of the Japanese lineage (0.77–1.32 × 10^–3). The θπ value of populations in the Japanese lineage was comparable to that of the threatened artiodactyl Budorcas taxicolor distributed in East Asia (0.45–0.75 × 10^–3) [21], suggesting that the Japanese lineage experienced a severe bottleneck effect and still harbors low genetic diversity. The value of observed heterozygosity (Het) showed large variation among populations in the continental lineage, possibly due to the situation that these populations live in nature reserves or farms with different area sizes and engage in non-random mating (Fig. 2C and Supplementary Fig. 8).

In addition, populations in the Japanese lineage showed an overall slow decay rate and a high level of LD (as measured by r²), whereas the populations in the continental lineage exhibited a rapid decay rate and a low level of LD (Fig. 2D). The difference in the genome-wide LD patterns between populations is consistent with the larger effective population size (Ne) for the continental lineage compared with the Japanese lineage. Moreover, the assessment of inbreeding level showed that populations in the Japanese lineage had higher inbreeding coefficient (F values: 1.83–4.56 × 10^–6) than those in the continental lineage (0.64–2.93 × 10^–6) (Fig. 2E).

Demographic history

To explore the demographic histories of sika deer, we inferred historical population size changes using the pairwise sequential Markovian coalescent (PSMC) model based on the distribution of heterozygous sites across the genome [37]. PSMC analysis showed that the demographic histories of sika deer can be traced back from the early Holocene (~ 10,000 years ago) to the middle Miocene (~ 10 Ma). The historical trends of Ne for different populations (Eurasian mainland vs. the Japanese archipelago) of sika deer showed a distinct pattern, especially during the period from the late Pliocene to the late Pleistocene (Fig. 3A).

Fig. 3

Demographic history of different sika deer populations. A The ancestral population size of each population was inferred using a PSMC model. A generation time (g) = 5 years and neutral mutation rate (µ) = 2.2 × 10^–9 per nucleotide site per generation were used. B Relative cross-coalescence rates (RCCR) over time between populations. Each line depicts the RCCR profile estimated using pairs of high-coverage genomes. The pairwise comparisons were between populations from Japan and RU1. C FastSIMCOAL2 v.2.6 simulation used to reconstruct the divergence, admixture, and demographic history of different sika deer populations. D TreeMix analysis of populations in Japan using RU1 as the outgroup when m = 2. E D-statistics with the form D (A, B; X, Y) for populations in Japan using RU1 as the outgroup. Positive D-statistic values indicate gene flow from X to A, whereas negative values indicate gene flow from X to B

We then used the multiple sequentially Markovian coalescent (MSMC) approach [38] to estimate the divergence time between the continental and Japanese lineages. We detected a decrease in the cross-coalescence rate between populations from the Japanese lineage and RU1 to 0.5 approximately 0.2–0.7 Ma (Fig. 3B). Furthermore, we also performed FastSIMCOAL2 v.2.6 analysis to infer demographic histories and estimate divergence times. The divergence between the continental and Japanese lineages occurred at approximately 0.60 Ma (Fig. 3C). Intriguingly, this divergence time coincided with the occurrence of the Naynayxungla Glaciation (0.72–0.50 Ma, corresponding to Marine Isotopic Stages (MIS) 14–18), which is speculated to be the most extensive glaciation among the four major glaciations during Quaternary climatic history [39]. Moreover, gene flow was detected from the population in southern China (ZJ) to the Japanese population in the late Pleistocene (Fig. 3C).

Next, we investigated potential gene flow among populations in the Japanese archipelago using two different approaches. In the TreeMix analysis, which included populations from Japan and population RU1, we detected ancestral gene flow under the optimal migration parameter (m = 2, as determined by OptM) (Supplementary Figure S9). Specifically, gene flow was inferred from the common ancestor of NK and HK to SM, as well as from TS to the common ancestor of NK and HK (Fig. 3D). Furthermore, D-statistics showed significant values (|Z-score|> 3) indicating gene flow from NK to HK, from NK to SM, from TS to YAK, from TS to SM, and from YAK to SM (Fig. 3E).

Genomic selection signatures associated with body size differences between different populations of wild sika deer

To explore the genetic mechanism of selection and local adaptation, we used two population pairs, YAK vs NK + SM and TS vs NK, to detect selective sweep signatures associated with body size (Supplementary Figure S10). We used these two population pairs for the following reasons. First, comparisons between population pairs that are genetically very similar would help to minimize the effects of population genetic differences that are not related to the selection signatures of interest. By reconstructing the phylogenetic relationship, we found that populations from HK, NK, SM, TS, and YAK were clustered together and each clade displayed a shorter branch length than the continental lineage, indicating that the genetic differentiation within the Japanese lineage was relatively shallow (Fig. 1B, D). Second, different populations in the Japanese lineage showed similar demographic histories (Fig. 3A). Therefore, using population pairs within the Japanese lineage could also minimize the effects of demographic changes in population size, such as bottlenecks or population expansions, on selective sweep analysis. Third, since both the TS and YAK populations are characterized by their smaller body size compared to the other populations in their respective pairs, the candidate genes under selection they share are likely to be linked to the small body size trait. Last, we used a population pair (TS vs NK) with relatively close latitudes, so as to reduce the influence caused by Bergman’s rule [40]. Due to similar body size [41], we did not conduct selective sweep analysis between TS and SM populations.

Furthermore, it is widely accepted that all methods detecting signatures of positive selection rely on different explicit or implicit assumptions and thus have limitations [42]. Therefore, using a range of methods in recombination is a powerful approach to resolving different forms of selection in practice [43]. In this study, we used five selective sweep methods, including the genetic differentiation statistic F_ST, the genetic diversity θπ, Tajima’s D, XP-CLR, and XP-EHH, to identify genomic regions under selection. The selection signature plots displaying all the outlier results across the chromosomes detected by five selective sweep methods were shown in Fig. 4A and B, and identified genes were listed in Supplementary table S8-S17. Those supported by all five methods (Supplementary Figure S11A and S11B) were considered as candidate genes associated with the body size. We identified 131 and 117 candidate genes in YAK vs NK + SM and TS vs NK, respectively, and 12 genes were shared between these two population pairs (Fig. 4C, Supplementary table S18-S22). Among the 12 candidate genes, four (XKR4, NPAS3, CTNNA3, and CNTNAP5) were supported in published papers suggesting a possible correlation with body size regulation (see details in discussion) [44,45,46,47,48,49,50,51,52,53]. Selective signals detected by the five selective sweep methods around these gene regions in YAK and TS populations were shown in Fig. 4D and E. Our data revealed a striking pattern at the XKR4 locus among populations. An enrichment of homozygous variant genotypes (green) was observed across SNPs around XKR4 region in population YAK, contrasting with the predominant heterozygous variants (yellow) in the NK + HK populations, suggesting this gene may be under selection in population YAK (Fig. 4F). Similarly, this gene was also likely to be under selection in population TS as suggested by Fig. 4G.

Fig. 4

Genomic regions with strong selection signals in the Yakushima (YAK) and Tsushima Island (TS) populations. A, B Selection signature plots displaying all the outlier results across the chromosomes detected by five selective sweep methods (FST, XP-CLR, XP-EHH, θπ, and Tajima’s D) in population pairs NK+SM vs YAK (A) and NK vs TS (B). C Candidate genes shared between population pairs NK+SM vs YAK and NK vs TS. D, E Selective signals detected by five selective sweep methods around the XKR4 gene region in YAK vs NK+SM populations (D) and TS vs NK populations (E). F, G The patterns of genotypes of the XKR4 gene region in YAK and NK+SM populations (F), and in TS and NK populations (G) based on 833 SNPs

Genomic selection signatures in domesticated sika deer

Sika deer have experienced strong artificial selection during the process of domestication, as their antler size and reproductive traits have changed substantially compared with those of their wild ancestors. Therefore, in this study, we resequenced 81 wild and 71 domesticated sika deer representing 10 main geographic populations and two farms to explore the genetic mechanism of selection during domestication. Aforementioned five selective sweep methods were used to identify genomic regions under selection (Fig. 5A, Supplementary table S23-S27). Genes supported by more than three methods were considered as candidate genes associated with the domesticated process (Fig. 5B and Supplementary table S28, S29) [54, 55].

Fig. 5

Genomic regions with strong selection signals in the domesticated sika deer population. A Selection signature plots displaying all the outlier results across the chromosomes detected by five selective sweep methods (F_ST, XP-CLR, XP-EHH, θπ, and Tajima’s D) in population pair wild vs domesticated. B Genes supported by more than 3 methods were considered as candidate genes in population pair wild vs domesticated. C Overlapping genes from three datasets, namely, genes detected by selective sweep analysis in this study, genes detected with transcriptomic data from different antlers [56], and antler-specific genes from a previous study [57], were used as candidate genes related to antler size to decrease the false-positive rate. D Selective signals detected by five selective sweep methods around the NRP2 (left panel) and EDIL3 (right panel) gene regions in domesticated populations

One of the most important traits under selection during the domestication of sika deer is antler weight, which increased by 30 ~ 60% during the process of domestication [28]. To screen genes related to antler growth, we selected genes that overlapped among three different datasets: genes detected by selective sweep analysis in this study, genes identified by transcriptomic data from different antlers [56], and antler-specific genes from a previous study [57]. We identified 17 overlapping genes between the selective sweep and antler-specific genes, five between the selective sweep and transcriptomic datasets, and two (NRP2 and EDIL3) among all three datasets (Fig. 5C and Supplementary table S30). The candidate genes NRP2 and EDIL3 were shown to be under selection in the domesticated population by the five selective sweep methods (Fig. 5D).

Discussion

The sika deer has historically been subdivided into 13 subspecies based on morphological characteristics [8,9,10]. As phylogenomic approaches can reveal discordances between molecular systematics and conventional morphological taxonomy [58, 59], we conducted phylogenomic investigations on sika deer. Based on the phylogenetic tree, PCA, and ADMIXTURE analysis of autosomal SNPs, together with mitochondrial and Y-chromosome evidence, we identified two lineages, the continental lineage and the Japanese lineage, with deep genetic divergence. Furthermore, Bayesian species delimitation decisively supported the species-level divergence between these lineages, with complete separation of posterior distributions (two-species model: PP = 1.0 vs single-species model: PP = 0). The divergence time between the two lineages was estimated to be approximately 0.60 Ma, which coincided with the occurrence of the Naynayxungla Glaciation (0.72–0.50 Ma; MIS 14–18), the most extensive of the four major glaciations during Quaternary climatic history [39]. This glacial interval has been proposed as a key driver of evolutionary dynamics in species distributed across the eastern Eurasian mainland [21, 60, 61].

The ancestral origins of sika deer in the Japanese archipelago were considered to derive from Eurasian continental populations that migrated via the land bridges repeatedly formed between Japan and Eurasian continent during the middle and late Pleistocene [12, 62]. Although previous paleontological investigation had posited a potential northern China origin based on fossil record [63], the precise continental progenitor of the Japanese lineage remains elusive due to their deep divergence time (~ 0.60 Ma). Our phylogenetic analysis demonstrated that populations from southern China and Taiwan (TW, ZJ and SC) exhibited a closer genetic relationship to the Japanese lineage compared to populations from northern China and Russia (RU1, RU2 and NEC) (Fig. 1B). In addition, ADMIXTURE results at K = 8 revealed that the genetic composition of the TW and ZJ populations were partially shared with populations from Japan (HK, NK, SM and YAK) (Fig. 1C). This admixed genetic structure likely reflected secondary contact events in late Pleistocene, as revealed by our demographic modeling demonstrating post-divergence gene flow from population ZJ to the Japanese lineages (Fig. 3C). These findings suggested complex evolutionary dynamics involving both early divergence and subsequent limited genetic exchange between the continental and Japanese lineages, complicating the ancestral origins inference of the Japanese lineage through contemporary genomic data alone.

Prior phylogeographic studies suggest that following their initial migration to Japan, sika deer ancestral populations likely inhabited in Honshu-Shikoku-Kyushu region [62, 63]. Nevertheless, some uncertainties remain regarding their demographic history, particularly the marked temporal discrepancy between mitochondrial and nuclear genomic data. Our whole mitochondrial genome analysis was consistent with previous D-loop and cytochrome b findings [12, 13], revealing two deeply divergent clades (0.53 Ma) (Supplementary Figure S7A). By contrast, the divergence times among populations in Japan estimated by genomic data were less than 0.10 Ma (Supplementary Figure S3 and S4). This inconsistency can be explained by a three-stage scenario of “divergence-mixing-isolation”.

Stage I, Divergence. Our mitochondrial data supported the divergence of the two clades in Japan, yet it remains unclear whether this divergence resulted from dual colonization events from the continent [12] or occurred post-colonization in Japan [64]. Notably, while a previous study detected no differentiation in Y-chromosomal markers among populations in Japan [16], our analysis of SNPs across the entire Y-chromosome clearly showed that populations in Japan were highly structured, with individuals from the same population forming distinct clusters (Supplementary Figure S7B). This discrepancy can be attributed to our more comprehensive dataset, which encompasses SNPs from the entire Y-chromosome, thus providing higher resolution. The differences in topological structures observed among populations based on mitochondrial and Y-chromosomal data were likely due to female philopatry and male-biased dispersal patterns in sika deer [65,66,67]. The substantial genetic divergence between the northern and southern sika deer clades in Japan renders the hypothesis of simultaneous continental divergence followed by concurrent migration biologically unlikely. While this dual-origin scenario could account for the existence of two distinct lineages, it fails to provide a parsimonious explanation for their pronounced latitudinal segregation. The invocation of female philopatry to reconcile this geographic pattern with the dual migration hypothesis is further weakened by evidence of Y chromosome-based genetic differentiation among populations (Supplementary Figures S7B). This genetic structure suggests that even male-mediated gene flow has been insufficient to homogenize these lineages, thereby challenging any simplistic model predicated solely on continental divergence and synchronous colonization.

Stage II, Mixing. Secondary contact and subsequent genetic admixture between the two clades occurred in Honshu-Shikoku-Kyushu [12, 62, 63]. The relatively recent divergence times among populations estimated by genomic data was likely attributed to this historical admixture event. In contrast, mitochondrial DNA-based estimates between the two clades remained unaffected by post-contact gene flow, as maternal inheritance of mitochondria allows the persistence of deeply divergent matrilineal lineages even when nuclear genomic differentiation becomes homogenized through introgression [68]. Moreover, TreeMix and D-statistics analyses further reveled asymmetry gene flow patterns, with significant migration detected from TS to NK, TS to HK, and NK to SM (Fig. 3D and 3E).

Stage III, Isolation. The sika deer in Japan were isolated to different islands leading to the formation of the current genetic structure (Supplementary Figure S3 and S4).

In addition, the NEC population is a domesticated group living in farms in northeast China. Our analysis of mitochondrial haplotypes revealed that this population had a relatively distant genetic relationship with other populations in the continental lineage (Fig. 1E). By contrast, biparental genomic data suggested that this population was more closely related to RU1 and RU2 than to other populations within the continental lineage (Fig. 1B). This discrepancy may be attributed to admixture of this population with RU1 and RU2, thereby reducing their genetic differences in nuclear markers.

Previous studies on patterns in ungulate body size variations emphasized the effect of a particular environmental factor, such as island area or precipitation, which reflects the availability of food resources [24]. This macroevolutionary phenomenon, termed the “island rule”, is pervasive across terrestrial vertebrates [26]. There is geographical variation in the body size of sika deer [8]. Substantial body weight variations among sika deer populations from the Japanese archipelago were reported: 50 kg in Kyushu (KY), 55 kg in Tanegashima (TN), 35 kg in Yakushima Island (YAK), and 30 kg in Kerama (KR) for adult males [8]. Moreover, the body size differences among different populations in the Japanese archipelago were verified in both adult male and female, respectively [41]. Hence, the conclusion that the body size of sika deer populations in Japan adheres to the island rule was neither influenced by disparate life histories nor attributable to sexual dimorphism [40]. By selective sweep analyses, we identified 12 genes shared between two small body size populations (YAK and TS) were under selection. Through a further literature search, we found four genes (XKR4, NPAS3, CTNNA3, and CNTNAP5) might be related to body size regulation as suggested by published papers. Among these genes, phospholipid scramblase XKR4 encoded a member of XK-related (XKR) membrane protein family, catalyzing the rapid, non-specific, and bi-directional translocation of phospholipids between the two leaflets of the plasma membrane [44]. It was reported to be associated with feed intake and growth phenotypes in cattle [45], carcass weight (CWT) and eye muscle area (EMA) in Korean Hanwoo cattle [46], growth and fatness traits in Sujiang pigs [47], and body size traits in Pekin ducks [48], indicating that this gene played a crucial role in regulating body size in various animals. NPAS3 was reported to play a critical role for lung development in mice [49]. In addition, a genome-wide association study suggested that it was candidate gene associated with growth and fatness traits in large white pigs [50]. CTNNA3 was reported to be involved in growth traits of red deer [51] and Qinchuan cattle [52]. CNTNAP5 was speculated to be a candidate gene associated with growth and carcass traits including body weight and body fat deposition in cattle [53].

These genes may contribute to the difference in body size between sika deer among the Japanese archipelago. To our knowledge, only one study has explored genes under selection that are related to body size differences between island and continental populations [27]. In that work, selective signatures in the fatty acid desaturase (FADs) gene cluster were inferred to be associated with stature size differences between an extant pygmy population on Flores Island and modern humans. In addition, other genes (PLAG1, HMGA2, LCORL, CCND2, TNN2, and TCP11) were found to regulate body size in cattle [69], suggesting that the genetic architecture of body size in large mammals is affected by many polymorphisms of small effect, as observed in cattle and humans [69, 70]. Therefore, more studies are needed to investigate whether the genes regulating body size related to the “island rule” in humans and other mammals are conserved.

Understanding the genetic changes underlying phenotypic variation in sika deer during domestic processes could facilitate our efforts toward further improvement. In this study, we detected genomic regions harboring genes associated with antler size. In particular, a strong selective signature located around two gene regions (NRP2 and EDIL3) implied intense artificial selection during domestication, suggesting that they may be involved in antler growth. Neuropilin-2 (NRP2) is a member of the neuropilin protein family. It was identified as a potential biomarker for chondrocyte programmed cell death in a recent study [71]. Epidermal growth factor (EGF)-like repeats and discoidin I-like domains-containing protein 3 (EDIL3) was reported to maintain cartilage extracellular matrix and protect chondrocyte apoptosis [72, 73]. Our data provide a valuable genomic resource for facilitating future genomics-informed breeding and genetic improvement of domestic sika deer.

Conclusions

Our analyses of whole genomes, Y chromosomes, and mitochondrial genomes revealed substantial genetic differentiation between the continental and Japanese lineages of sika deer, indicating that they belong to two distinct phylogenetic species. The divergence time between the two species was estimated to be approximately 0.60 Ma, which coincided with the Naynayxungla Glaciation (0.72–0.50 Ma, MIS 14–18). Sika deer in Japan were inferred to have experienced a “divergence-mixing-isolation” evolutionary scenario. In addition, we also identified four candidate genes (XKR4, NPAS3, CTNNA3, and CNTNAP5) that may be involved in body size regulation. Furthermore, we also detected two candidate genes (NRP2 and EDIL3) that is putatively associated with an important economic trait (antler weight) under selection during the process of domestication. The genomic data generated in this study will serve as a valuable resource for genomics-assisted breeding to improve economic traits in domesticated sika deer.

Limitation of the study

Our work has several limitations. First, sika deer are listed as an endangered animal in China’s Red Book of Endangered Animals (National Forestry and Grassland Administration, National Park Administration, 2021). Therefore, the number of sika deer that could be collected and used for this study was limited. Second, some population genetic parameters may be affected by sample size, for example, genetic diversity, linkage disequilibrium and even selective sweep statistics. Third, body size has been considered one of the most important traits of an organism or a species, as it covaries with a wide range of morphological, physiological, and ecological features. Our work is only a preliminary analysis of this complex trait.

Methods

Sample collection, DNA extraction, and genome sequencing

We collected blood, ear, antler or skin samples from 81 wild and 71 domesticated sika deer representing 10 main geographic populations and two farms, respectively. During the process of sampling, the deer were anesthetized by intramuscular injection of 5% xylazine with a dose of 0.5 mg/kg body weight, which was performed in accordance with the ARRIVE guidelines (https://arriveguidelines.org). The deer were not treated harshly during sampling and were released after sampling. Of the 81 wild individuals, 16 were from a nature reserve in the Far East of Russia (RU1, including 10 individuals collected in Kazakhstan introduced from this nature reserve in the 1930s), 6 were from the wild region in the Far East of Russia (RU2), 8 were from the Tiebu Nature Reserve in Sichuan Province (SC), 8 were from the Qingliangfeng National Nature Reserve in Zhejiang Province (ZJ), 2 were from the Liugongdao National Forest Park in Shandong Province (introduced from Taiwan on April 16, 2011, TW), 9 were from Hokkaido (HK), 8 were from Nikko in Honshu (NK), 2 were from Shimane in Honshu (SM), 14 were from Tsushima Island (TS), and 8 were from Yakushima Island (YAK). The last five sampling sites are located in Japan. In addition, a total of 71 domesticated sika deer from farms in Jilin and Heilongjiang Provinces were also included in this study (Fig. 1a and Supplementary table S1).

Russian and Kazakhstani samples were provided by Professor Nurlan Toktarov. Chinese samples were collected by Huamiao Liu, Yimeng Dong, Yan Ju, Yang Li, Weilin Su, Ranran Zhang, Shiwu Dong, Hongliang Wang, Yongna Zhou, Yanmin Zhu, Lei Wang, Zhengyi Zhang, Pei Zhao, Shuyan Zhang, Rui Guo, E A, Yuwen Zhang, and Xin Liu. DNA of the Japanese samples were provided by Professor Hidetoshi B. Tamate.

Genomic DNA of all samples sequenced in this study was extracted using the DNeasy Blood and Tissue Kit (QiaGen, Shanghai, China), including RNase A treatment. DNA quality and quantity were evaluated using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and electrophoresis on a 0.8% agarose gel, respectively. Individual DNA libraries with an insert size of 350 bp were constructed with at least 1 µg of genomic DNA according to the manufacturer’s protocol and then sequenced using the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA).

Sequencing read filtering and SNP variant calling

In total, we obtained ~ 13.52 terabases (Tb) of raw sequencing data (~ 88.95 Gb per sample). To avoid reads with artificial bias (i.e., low-quality paired reads, which primarily result from base-calling duplicates and adapter contamination), we used fastp v0.23.4 [74] Wheeler aligner (BWA-mem) software version 0.7.84] to remove the following types of reads: (1) reads with ≥ 10% unidentified nucleotides (N); (2) reads with > 10 nt aligned to the adapter, with ≤ 10% mismatches allowed; and (3) reads with > 50% bases having a Phred quality < 5. Consequently, we retained 13.22 Tb (~ 86.94 Gb per sample) of high-quality genomic data (Supplementary table S2).

The remaining high-quality paired-end reads were aligned to the chromosome-level reference genome of sika deer assembled from a domesticated individual [4] and the Y chromosome of GCA_040085125.1 in NCBI RefSeq database (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_040085125.1/) using Burrows-Wheeler aligner (BWA-mem) software version 0.7.8 [75] with the command “mem -t 4 -M -R”. We then converted the mapped reads into BAM files and sorted the files using SAMtools v1.9 [76] with the parameter setting “-bS -t”. PCR duplicates were subsequently removed using “MarkDuplicates” in GATK version 4.1.2.0 [77]. We calculated the whole-genome sequencing coverage and depth of each sample using SAMtools. Only samples with miss < 10% and depth > 5 × were kept to call SNP variations.

SNPs and InDels were called from the BAM files using the HaplotypeCaller module in GATK version 4.1.2.0 with the “-ERC GVCF” option. Raw GVCFs with the samples called individually were merged using CombineGVCFs and called for SNPs using GenotypeGVCFs. We then selected the candidate SNPs and created the selected SNP dataset using the SelectVariants module in GATK. To obtain reliable candidate SNPs and InDels, we implemented “VariantFiltering” in GATK according to the following parameter settings: “QUAL < 30.0 | QD < 2.0 | MQ < 40.0 | FS > 60.0 | SOR > 3.0 | MQRankSum < -12.5 | ReadPosRankSum < -8.0” for SNPs and “QD < 2.0 | QUAL < 30.0 | FS > 200.0 | ReadPosRankSum < -20.0” for InDels.

After filtering, all the identified SNPs and InDels were further annotated using the ANNOVAR software tool v2022-12–28 [78] and subsequently divided into the following groups on the basis of the sika deer genome annotation information: variants occurring in intergenic regions, within 1 kb upstream of the translation start and downstream of the translation stop sites, in coding sequences, and in introns.

Structural variant and copy number variant calling

Only samples with depth > 15 × were kept to identify structural variants (SVs). SVs, including INS, DEL, ITX, CTX, and INV, were detected by LUMPY v.0.2.13 [79] and BreakDancer v1.4.5 [80] with the following criteria: 50 bp ≤ SV ≤ 1 Mb. CNV calling was performed using CNVnator 0.3.3 [81].

Mitochondrial genome assembly and haplotype analysis

Clean reads were mapped to a complete mitochondrial genome (GenBank: GU457433.1) to generate a SAM file for each individual. After filtering potential PCR duplications by SAMtools software, the resulting BAM file was used as input for MITObim version 1.9.1 [82] to assemble the mitochondrial genome of each individual with the parameters “-start 0 -end 30”.

The mitochondrial genome alignment was refined using TrimAl v1.4.1 [83] with automated parameter optimization (-automated1 mode), which removed ambiguously aligned regions while preserving phylogenetically informative positions, resulting in a conserved alignment length of 15,980 bp. From this curated alignment, we identified 914 variable sites. Thirty-nine unique haplotypes were resolved in wild sika deer populations using DnaSP version 5.10.01 [84]. Bayesian phylogenetic reconstruction and divergence time estimation were performed in BEAST version 2.6.6 [85] using an uncorrelated log-normal relaxed clock, with a calibration node prior set at 1.4 Ma [95% Highest Posterior Density (HPD): 0.73–2.20 Ma] as the divergence time between Cervus elaphus and Cervus nippon [86]. Additionally, a maximum likelihood (ML)-based unrooted phylogeny was reconstructed in IQ-TREE v2.2.0 [87].

Y chromosome SNPs and haplotype analysis

Following aforementioned GATK variant calling procedures, we identified 543,807 Y-chromosomal SNPs across 129 male sika deer. Subsequent quality filtering with TrimAl (-automated1 mode) retained 211,326 SNPs, from which 52 distinct haplotypes were resolved in wild populations. A neighbor-joining tree was reconstructed in MEGA11 [88] under the Kimura two-parameter (K2P) model [89], with nodal support assessed through 100 bootstrap replicates.

After removing heterozygous sites and sites with > 5% missing genotypes, we obtained 1,248 high-confidence SNPs that delineated 35 Y-chromosomal haplotypes. A median-joining network was constructed using NETWORK v10.1 (https://www.fluxus-engineering.com/) [90].

Population genetic analysis

To reduce the required computing resources, we implemented LD-based pruning for the genotype data using PLINK software version 1.9 [91] with the option “-indep-pairwise 50 10 0.2” for the autosomal SNP dataset. To clarify phylogenetic relationships from a genome-wide perspective, we constructed an NJ tree according to the p distance using the software TreeBeST v1.9.2 [92] with 1000 bootstrap replications. Furthermore, population genetic structure was examined via an expectation maximization algorithm, as implemented in the program ADMIXTURE v1.23 [93]. The number of assumed genetic clusters K ranged from 2 to 9, with 10,000 iterations for each run. CV error minimization was used to determine the optimal K value [94]. In addition, principal component analysis (PCA) was performed using PLINK software [91].

To infer the species tree and estimate the divergence times among the major clades, we used the SNAPP package version 1.5.2 [33] implemented in BEAST version 2.6.6 [85] to perform phylogenetic analysis. Due to the high computational demand of SNAPP, we performed this analysis using only four individuals for each population and a reduced set of 5,000 SNPs randomly selected from all biallelic SNPs with no missing genotypes. The input files for SNAPP were prepared with the script snapp_prep.rb [95], implementing a strict-clock model and a pure-birth tree model. We used the divergence time of Cervus elaphus and Cervus nippon of 1.4 Ma (95% HPD: 0.73–2.20 Ma) [86] as the calibration point. We carried out three independent SNAPP analyses, each with a Markov chain Monte Carlo (MCMC) chain length of 1,500,000 and parameter values recorded every 250 generations, resulting in ESS values that were all greater than 400. The posterior tree distribution was again summarized as a maximum clade credibility (MCC) consensus tree. To verify the repeatability of SNAPP analysis, we repeated this analysis by re-selecting the reduced set of 5,000 SNPs randomly for another 9 times, and the topological structures of species trees were integrated into one figure (Supplementary Figure S4).

Species delimitation

We used BPP v. 4.3.0 [35] for joint species delimitation and species tree inference with unguided species delimitation [96]. Due to the relatively high computational burden of Bayesian methods such as BPP, we performed inference on reduced datasets consisting of a subset of 1000 of the most variable loci, which were selected by Geneious v9.1.7 [97] based on the percentage of variable sites. Three independent MCMC runs were implemented with a burn-in of 10,000, a sample frequency of 2, and 100,000 total samples.

Genetic diversity

Genetic differentiation statistics (F_ST), observed heterozygosity (Het), and nucleotide diversity (θπ) were calculated using VCFtools v0.1.14 [98], with sliding windows of 40 kb and 20-kb overlap between adjacent windows. To detect the genomic signatures of inbreeding in different populations of sika deer, genome-wide runs of homozygosity (ROH) were identified for each individual using the tool in PLINK software version 1.9 [91] with the parameter setting “-threshold 0.05 -snp 65 -kb 100 -missing 3 -gap 5,000 -density 5,000”.

We estimated the genome-wide LD decay of each population using PopLDdecay version 3.40 [99], setting the maximum distance between two SNPs to 500 kb, MAF ≥ 0.05, and missing ratio ≤ 10%. In addition, PLINK version 1.9 [91] was used to calculate the LD coefficient (r²) between pairs of high-quality SNPs, with the parameters set as “-ld-window-r2 0 -ld-window 99,999 -ld-window-kb 1000”. The average r² value was calculated for pairs of markers in a 10-kb window, and values were averaged across the whole genome.

Demographic history and divergence time analyses

To trace the demographic history of each sika deer population, we employed the PSMC model [37] to reconstruct changes in historical effective population size (N_e) over time using heterozygous sites. Analyses for each population were based on the two to three individuals with the highest mean sequencing depth. The PSMC analysis was run with the parameters -N25 -t15 -r5 -p “4 + 25 *2 + 4 + 6”, together with 100 bootstrap replicates. Referring to Tarim red deer (C. e. yarkandensis) and other species in Cervidae, we adopted a mutation rate of 2.2 × 10^–9 per nucleotide site per generation [100, 101] and a generation time of 5 years to scale the results.

The MSMC method [38] was used to estimate the relative cross-coalescence rate (RCCR) between pairs of high-coverage individuals. The time at which the RCCR decreased below the 0.5 threshold was taken as a point estimate of the divergence time. For the range of divergence times, we used the first time point at which the cross-coalescence rate was at or below 0.25 and 0.75. All parameter estimates were scaled using a mutation rate of 2.2 × 10^–9 per nucleotide site per generation [100, 101] and a generation time of 5 years to scale the results.

Demographic histories and migration rates among RU1, ZJ and Japan (JPN) were estimated by using FastSIMCOAL2 v.2.6 [102], using a mutation rate of 2.2 × 10^–9 per nucleotide site per generation [100, 101] and a generation time of 5 years. We first inferred the best parameters and the likelihoods for each of the demographic models: no gene flow, early gene flow and recent gene flow, since the divergence between the continental populations and JPN by running 100 independent iterations with 300,000 coalescent simulations and 60 optimization cycles. The demographic model with the highest likelihood was selected to run parameter estimation with block-bootstrap of 100 replicates.

To detect potential gene flow among populations in Japan, we implemented two distinct statistical analyses: TreeMix v1.13 [103] and D-statistics [104]. For the TreeMix analysis, we conducted maximum likelihood trees using blocks of 1,000 SNPs and RU1 as the root, setting the number of tested migration events from 1 to 5. Subsequently, the optimal number of migration edges was determined by OptM [105]. Additionally, we calculated D-statistics using the program AdmixTools [104]. The four populations were set to be ((A, B), (X, Y)), where Y represented the outgroup RU1, and A, B, X were selected population combinations according to the phylogenetic tree constructed by autosomal SNPs (Supplementary Figure S3 and S4). The statistical significance of the D value was assessed using a two-tailed Z test, with |Z-score|> 3 to be significant [106].

Genomic signatures of selection for body size and domestication

To identify genome-wide selective sweeps associated with body size, we used two population pairs: YAK (n = 8) vs NK (n = 8) + SM (n = 2), and TS (n = 14) vs NK (n = 8). In addition, the wild population (n = 81) and domesticated population (n = 71) were compared to detect selective signatures related to the domestication process. We used five selective sweep methods to identify genomic regions under selection, including the genetic differentiation statistic F_ST, the genetic diversity θπ, Tajima’s D, XP-CLR, and XP-EHH. F_ST, θπ and Tajima’s D were calculated using VCFtools v0.1.14 [98], with sliding windows of 40 kb and 20-kb overlap between adjacent windows. We also scanned for signals of selection across the genome using XP-EHH scores [107] as implemented in the Python package “ALLEL” with a sliding window of 40 kb in 20 kb step size. The XP-CLR statistic was calculated using the XP-CLR v1.1.2 [108]. According to the simple majority vote strategy [54, 55], genes supported by more than three methods can be considered as candidate genes related to the phenotype. To identify candidate genes related to body size, we adopted a stricter criterion (supported by all five methods) to decrease the false-positive rates. To identify candidate genes related to antler weight, we adopted a common criterion (supported by more than three methods), since we already used two different datasets (genes identified by transcriptomic data from different antlers [56], and antler-specific genes from a previous study [57]).