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.

Keywords: Cervus nippon; Domestication; Genetic selective signature; Population genomics; Sika deer; Species delimitation.

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

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

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

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

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