Evolutionary Conservation Genomics Reveals Recent Speciation and Local Adaptation in Threatened Takins

Abstract Incorrect species delimitation will lead to inappropriate conservation decisions, especially for threatened species. The takin (Budorcas taxicolor) is a large artiodactyl endemic to the Himalayan–Hengduan–Qinling Mountains and is well known for its threatened status and peculiar appearance. However, the speciation, intraspecies taxonomy, evolutionary history, and adaptive evolution of this species still remain unclear, which greatly hampers its scientific conservation. Here, we de novo assembled a high-quality chromosome-level genome of takin and resequenced the genomes of 75 wild takins. Phylogenomics revealed that takin was positioned at the root of Caprinae. Population genomics based on the autosome, X chromosome, and Y chromosome SNPs and mitochondrial genomes consistently revealed the existence of two phylogenetic species and recent speciation in takins: the Himalayan takin (B. taxicolor) and the Chinese takin (B. tibetana), with the support of morphological evidence. Two genetically divergent subspecies were identified in both takin species, rejecting three previously proposed taxonomical viewpoints. Furthermore, their distribution boundaries were determined, suggesting that large rivers play important roles in shaping the genetic partition. Compared with the other subspecies, the Qinling subspecies presented the lowest genomic diversity, higher linkage disequilibrium, inbreeding, and genetic load, thus is in urgent need of genetic management and protection. Moreover, coat color gene (PMEL) variation may be responsible for the adaptive coat color difference between the two species following Gloger’s rule. Our findings provide novel insights into the recent speciation, local adaptation, scientific conservation of takins, and biogeography of the Himalaya–Hengduan biodiversity hotspot.


Introduction
The fundamental goals of species delimitation are not only to identify species, but also to clarify the phylogenetic relationships among and within the species. For threatened species, incorrect taxonomic classification will lead to inappropriate conservation decisions, and thus, improper protection and management effects (Gutierrez and Helgen 2013;Zachos 2013). The takin (Budorcas taxicolor) is a largely threatened artiodactyl that was once widely distributed in East Asia. However, as a forest-dwelling species, it currently only exists in the eastern Himalayas and southwest China at elevations of 1500-4500 m, with its distribution range including China, Myanmar, India, and Bhutan. The takin populations have gradually decreased due to human activities such as hunting and deforestation, and resultant habitat loss and fragmentation (Song et al. 2008).

Chromosome-level Genome Assembly and Annotation
We produced the first high-quality chromosome-level genome of male takin, using a combined strategy of Illumina, Nanopore, and Hi-C technologies (supplementary table S2, Supplementary Material online). Through the pure Nanopore long-read assembly strategy, a takin genome of 2.6 Gb was assembled, which consisted of 230 contigs with a contig N50 of 47.9 Mb (supplementary table S3, Supplementary Material online). These contigs were further assembled into 27 pseudochromosomes (PCHRs) at the chromosome scale, with an assembly mount ratio of 99.99% ( fig. 1a, and supplementary fig. S1 and table S4, Supplementary Material online). Benchmarking Universal Single-Copy Orthologs (BUSCO) (Simão et al. 2015) and Core Eukaryotic Genes (CEG) (Parra et al. 2007) assessments showed that 3,880 of 4,104 BUSCO genes (94.6%) were complete and that 234 of 248 CEG genes (94.35%) were complete and partial, indicating high completeness of the assembled genome ( S9, Supplementary Material online). The genome synteny evaluation between the PCHRs of takin and the chromosomes of goat showed a synteny ratio of 97.23%, indicating a high assembly quality of the takin genome ( fig. 1b). PCHRs LG05 and LG27 corresponded to the X and Y chromosomes of takin, respectively ( fig. 1b).

Phylogenomics and Adaptive Evolution of Takin
To reveal the phylogenetic status of takin within Bovidae, we constructed a genome-wide phylogenetic tree including the genomes of 20 other Ruminantia species (supplementary table S10, Supplementary Material online). The phylogenetic tree based on a total of 3.1 million syntenic 4-fold degenerate (4Dtv) sites showed that takin was positioned at the root of Caprinae ( fig. 1c).
To detect the potential genetic basis of large body size in takins as compared with the other Caprinae species, we performed positive selection and rapid evolution analyses using the branch-site and branch models, respectively, implemented in PAML (version 4.9) (Yang 2007). We identified 130 positively selected genes (PSGs) and 333 rapidly evolving genes (REGs) (supplementary tables S11 and S12, Supplementary Material online). The Gene Ontology (GO) enrichment analyses of both PSGs and REGs revealed significant terms related to body size development, including developmental growth (GO:0048589, P = 4.71E−3), head development (GO:0060322, P = 1.53E−3), muscle tissue development (GO:0060537, P = 5.53E−3), and skeletal muscle tissue regeneration (GO:0043403, P = 8.37E−3) (supplementary tables S13 and S14, Supplementary Material online). Among these genes, four genes (BMP3, WNT7A, TGFBR3, and UNC45B) are involved in skeletal development, bone osteogenesis, and muscle structure development (supplementary table S15, Supplementary Material online). The loss of WNT7A and TGFBR3 in humans or mice results in severe skeletal defects and limb malformations (Woods et al. 2006;Kirkbride et al. 2008;Al-Qattan et al. 2013;Hill et al. 2015). These genes may contribute to the difference in body size between takin and the other Caprinae species. To detect the speciesspecific selection signals in takin, we performed selection analyses using all of the other Bovidae species as the background branches, and detected 154 PSGs and 379 REGs (supplementary tables S16 and S17, Supplementary Material online). The GO enrichment analysis showed significant terms related to cartilage development (GO:0051216, P = 4.80E−3), muscle organ development (GO:0007517, P = 7.93E−3), lipid biosynthetic process Yang et al. · https://doi.org/10.1093/molbev/msac111 MBE FIG. 1. Genomic landscape, genome synteny, and genome-wide phylogenetic tree of takin. (a) Circos plot of a chromosome-level takin genome. From outer to inner: (I) sizes of 27 PCHRs (LG01-LG27); (II) gene density; (III) repeat sequence distribution; (IV) GC content. The photo in the circle shows takin from the Sichuan subspecies. (b) Genome synteny between takin (right) and goat (left). Collinear blocks between the two species are linked by lines with the same colors. (c) The genome-wide phylogenetic tree of 21 species from 9 subfamilies within the Bovidae based on 4Dtv sites. Two calibration points (red node) were used for estimating the divergence time. The takin marked with red was at the root of Caprinae.

Genomic Evidence of Two Phylogenetic Species in Takins
We resequenced the genomes of 75 wild takins from seven mountains or regions, including Qinling (QIN, n = 20), Minshan (MIN, n = 17), Qionglai (QIO, n = 15), Xiaoxiangling (XXL, n = 1), Gaoligong (GLG, n = 4), Southeast Tibet to the east of the Yalu Zangbu River (SETE, n = 11), and Southeast Tibet to the west of the Yalu Zangbu River (SETW, n = 7) ( fig. 2a  Because of the small sample size, this individual was excluded from any population-level analyses. Takin individual MIN01, collected in the junction region (Qingmuchuan) between the Qinling and Minshan Mountains, was genetically clustered with MIN rather than QIN, suggesting that the takins from this region belong to the MIN population and the Jialing River, rather than the traditionally thought Bailong River, is most likely the boundary of B. t. bedfordi and B. t. tibetana ( fig. 2a).
In addition to the evidence from autosomal SNPs, X-chr SNPs, paternally inherited Y-chr SNPs, and maternally inherited mtgenome data also revealed the existence of two divergent genetic lineages (fig. 2e-f; Supplementary fig. S4 and S5, Supplementary Material online). The Network maps and phylogenetic trees of the Y-chr SNPs and mtgenome haplotypes and the phylogenetic tree of X-chr SNPs consistently showed that GLG-SETE-SETW clustered together (the Himalayan lineage), separately from the QIN-MIN-QIO-XXL cluster (the Chinese lineage), highlighting the substantial genetic divergence between the two lineages. Furthermore, genetic differentiation was found between SETW and GLG-SETE and between QIN and MIN-QIO-XXL, similar to the results obtained for autosomal SNPs. Particularly, all the haplotypes of Y-chr SNPs and mtgenome for four hybrid SETE individuals (SETE8/9/10/11) clustered with those of GLG-SETE (supplementary fig. S5, Supplementary Material online). So, combining the above results with the ADMIXTURE result and geographical distribution, we classified the four hybrid SETE individuals into the GLG-SETE population (i.e., the Mishmi subspecies).
Furthermore, F ST measure based on whole-genome SNPs found that F ST between the Himalayan and Chinese lineages was 0.2925, much higher than those between sublineages of the same lineage (GLG-SETE vs. SETW, F ST = 0.2045; QIN vs. MIN-QIO-XXL, F ST = 0.1613). F ST between two genetic populations of Sichuan subspecies, MIN and QIO-XXL, was 0.069, much lower than the above measures (supplementary table S23, Supplementary Material online). To further test whether the takins were classified as two, three, four, or five phylogenetic species, we performed species delimitation analyses using the SNAPP method (Grummer et al. 2014). SNAPP results showed that the most significantly supported species delimitation was two species, namely GLG-SETE-SETW and QIN-MIN-QIO-XXL, with a marginal likelihood estimate being the highest (supplementary table S24, Supplementary Material online).
In summary, all of the results based on autosomal SNPs, X-chr SNPs, Y-chr SNPs, and the mtgenomes consistently revealed two evolutionarily significant units, showing substantial genetic divergence between them (F ST = 0.2925). SNAPP species delimitation results supported the classification of two species. Following the phylogenetic species concept, we propose the existence of two phylogenetic species in takins: the Himalayan takin B. taxicolor Hodgson, 1850, and the Chinese takin B. tibetana  . 3b). Therefore, takins can be classified into two species and four subspecies, rejecting the three previously proposed taxonomic viewpoints based only on morphological evidence. The molecular taxonomy was also supported by morphological differences, especially the coat color and body size (supplementary table S1, Supplementary Material online). Specifically, the Himalayan takins have darker coat color (brown-dark or black) than the Chinese takins (white, golden yellow, or light brown-gray); and the Himalayan takins have smaller body size (lighter and shorter body) than the Chinese takins. These two species have a difference in altitude distribution, with the Himalayan takins living at a Yang et al. · https://doi.org/10.1093/molbev/msac111 MBE much higher starting elevation than the Chinese takins (supplementary table S1, Supplementary Material online). The distribution boundary of the two species is the Three Parallel Rivers (Nujiang, Lancangjiang, and Jinshajiang), which stemmed from the collision of the Indian and Eurasian Plates and formed in the middle Pleistocene or earlier (Li et al. 2008;Nie et al. 2018); however, which river being the most important boundary will require further study with more samples across this region (figs. 2a and 3b). The distribution boundary for the Mishmi and Bhutan subspecies is most likely the Yalu Zangbu River, which also stemmed from the collision of the Indian and Eurasian Plates and formed in the early Miocene (Wang et al. 2002), and that for the Sichuan and Qinling subspecies is most likely the Jialing River that was built in the middle Pleistocene (Zhang et al. 2008), as suggested by the genomic findings (figs. 2a and 3b).

Demography, Divergence Histories, and Gene Flow
Pairwise sequential Markovian coalescent (PSMC) analysis (Li and Durbin 2011) showed that the demographic histories of takins could be traced back from the early Pleistocene ( 2 Ma) to the early Holocene ( 10,000 years ago [ka]). On the whole, the demographic trajectories of the two takin species were different, although both had experienced two population bottlenecks and one population expansion ( fig. 3a). The two species underwent an obvious population decline 0.70 Ma, which coincided with the occurrence of the Naynayxungla Glaciation (0.78-0.50 Ma) (Zheng et al. 2002). The population decline resulted in the first bottleneck 0.25 Ma, most likely caused by penultimate glaciation (0.30-0.13 Ma). After the glaciations, the populations started to expand and reached a maximum 80 ka, which coincided with the occurrence of the last interglacial period. Then, the occurrence of the last glaciations resulted in another rapid population decline, and the second bottleneck occurred during the Last Glacial Maximum ( 20 ka). A split of the effective population size (Ne) curves between the two species occurred 0.25 Ma, possibly suggesting the initial divergence of the two species at that time. In addition, the size of population expansion was obviously different in the two species. The Ne of the Himalayan takin was 1.6 times that of the Chinese takin at the peak of population expansion, suggesting that the Himalayan takin had a larger population size than the Chinese takin during the last interglacial period.
To reveal the detailed species/population divergence history of takins, we further performed a Fastsimcoal2 simulation (Excoffier and Foll 2011;Excoffier et al. 2013 . 3b). These two divergence events most likely resulted from the last glaciations and geographical isolation imposed by the Jialing River and Minjiang River. Ne estimation showed that after the divergence of the two species, the Ne of the Himalayan takin decreased more slowly than that of the Chinese takin. Furthermore, the Ne of QIN and QIO continued to decline after divergence, whereas the Ne of SETW, GLG-SETE, and MIN did not change much (

Genomic Signatures of Selection and Local Adaptation in the Two Species
Considering the differences in the morphology (such as coat color and body size), geographical distribution of the two species (supplementary table S1, Supplementary Material online), and their genetic divergence ( 0.23 Ma), we identified signatures of selection and local adaptation in the Chinese takin and the Himalayan takin using the cross-population composite-likelihood ratio (XP-CLR) test (Chen et al. 2010a), in which genomic regions with extreme allele frequency differentiation were scanned based on the allele frequency spectrum (AFS). Employing the top 0.5% XP-CLR values, we identified 270 and 314 selected genes in the Chinese and Himalayan takins, respectively (supplementary tables S27 and S28, Supplementary Material online). These selected genes were significantly enriched in the TGF-beta signaling pathway (ko04350, P = 1.53E−04) and biological processes involved in hair follicle development (GO:0001942, P = 7.20E−3), skin epidermis development (GO:0098773, P = 8.53E−3), ossification (GO:0001503, P = 1.40E−4) Recent Speciation and Local Adaptation in Threatened Takins · https://doi.org/10.1093/molbev/msac111 MBE (supplementary tables S29 and S30, Supplementary Material online), which may contribute to the differences in coat color and body size between the two species.
Coat color is the most significant morphological difference between the two species ( fig. 3b). The Himalayan takin (brown-dark or black) shows a darker coat coloration than the Chinese takin (white, golden yellow, or light brown-gray) ( fig. 3b; supplementary table S1, Supplementary Material online). The color pattern of takins is consistent with Gloger's rule (Delhey 2017) that darker coat colors are often found in humid and warm environments (i.e., the Himalayan takin), whereas lighter individuals are more prevalent in dry and cold areas (i.e., the Chinese takin). Interestingly, one coat color gene, PMEL (encoding a melanocyte-specific type I transmembrane glycoprotein), was shown to be under selection in the Himalayan takin by both the XP-CLR and  (Kuehn and Weikard 2006;Gutierrez-Gil et al. 2007;Schmutz and Dreger 2013), horses (Brunberg et al. 2006), mice (Hellstrom et al. 2011), dogs (Clark et al. 2006), and chickens (Kerje et al. 2004), suggesting that genetic convergence played an important role in the coat color dilution of different species. This finding suggested that the PMEL variant might contribute to the adaptive lighter coat color of the Chinese takin, providing an important candidate gene for future functional verification.
Genomic Diversity, Linkage Disequilibrium, Inbreeding, and Genetic Load To determine the genetic evolutionary potential of the two takin species and four subspecies, we analyzed genomic variation, linkage disequilibrium (LD), inbreeding, and genetic load levels on the basis of whole-genome SNPs. Population-level genomic variation analysis showed that genetic diversity (θ π and θ w ) of the Himalayan takin was higher than that of the Chinese takin on the basis of autosome, X-chr, and Y-chr SNPs, and mitogenome data ( fig. 5a- . Compared with the other ungulates or threatened mammals, the genomic diversity of the two takin species was lower than that of wild yak, golden snub-nosed monkey, and giant panda, similar to that of the Chinese red panda and higher than that of the Himalayan red panda (supplementary table S35, Supplementary Material online including related references). A single-genome level analysis demonstrated that genomic heterozygosity of the Himalayan takin was higher than that of the Chinese takin. By contrast, genomic heterozygosities of the two takin species were lower than those of Tibetan antelope, wild Bactrian camel, and Ursidae species, similar to those of wild horse, Guizhou snub-nosed monkey, tiger, and lion, and higher than those of snow leopard, baiji, and other snub-nosed monkey species (supplementary table S36, Supplementary Material online including related references).
Genomic LD (measured as R 2 ) analysis indicated that the Qinling subspecies (QIN) presented the highest level of LD and slower LD decay than other subspecies, with R 2 becoming stable at a distance of 300 kb, whereas faster LD decay was found in the other three subspecies, with a stable R 2 value found at a distance of 200 kb ( fig. 5e). The assessment of inbreeding levels based on the runs of homozygosity (ROH) demonstrated that the Qinling subspecies presented the highest inbreeding level among the four subspecies ( fig. 5f). Because deleterious mutations should be removed more efficiently in larger populations (Charlesworth 2009;Xue et al. 2015), we further evaluated the genetic load of deleterious variations in takins according to the methods of Hu et al. (2020). We investigated the distributions of three types of variations (deleterious, loss of function [LoF], and synonymous mutations) in proteincoding genes. We found that the ratios of deleterious variants to synonymous variants at both homozygous-derived sites and nonhomozygous-derived sites were higher in the Qinling subspecies than in the other three subspecies ( fig. 5g). The ratios of nonhomozygous-derived LoF variants to nonhomozygous-derived synonymous variants were the highest in the Qinling subspecies; however, the ratios of homozygous-derived LoF variants to homozygous-derived synonymous variants were comparable in the four subspecies ( fig. 5h).
In summary, the above results consistently showed that compared with the other three subspecies, the Qinling subspecies (B. t. bedfordi) presented the lowest genomic diversity, the highest LD level, and relatively higher inbreeding and genetic load levels, highlighting that this subspecies is in the most urgent need of conservation and genetic management.

Discussion
We produced the first high-quality chromosome-level genome of takin using Illumina, Nanopore, and Hi-C technologies, which provides an important genetic resource for the research and conservation of takins.
Yang et al. · https://doi.org/10.1093/molbev/msac111 MBE The phylogenomic tree construction including 21 species from nine subfamilies in Bovidae confirmed that takins belong to Caprinae. Takins are well known for their threatened status and unique morphology when compared with the other Caprinae species. In this study, we detected adaptive divergence signatures in takin through comparisons with the other Caprinae species. Four genes (BMP3, WNT7A, TGFBR3, and UNC45B) associated with body size development were found to be under positive selection or rapid evolution in takin, suggesting that Due to the low population size and elusive habit, the intraspecies taxonomy of takins has long remained unresolved. Here, based on the combined evidence of genetic (i.e., multi-marker genetic structure, F ST , species delimitation results, and different demographic histories), morphological (i.e., coat color and body size), and biogeographic differences, we defined two phylogenetic species in takins: the Himalayan takin and the Chinese takin. Furthermore, we identified two subspecies of the Himalayan takin, the Mishmi and Bhutan subspecies, and two subspecies of the Chinese takin, the Sichuan and Qinling subspecies. The distribution boundary between the Sichuan subspecies and the Qinling subspecies has traditionally been considered to be the Bailong River. In the present study, individual MIN01, collected in Qingmuchuan to the east of the Bailong River, was genetically clustered with Minshan (MIN) rather than Qinling (QIN) ( fig. 2a). In addition, a previous mitochondrial-DNA-based study showed that one takin individual from Maozhai to the east of the Bailong River was also genetically clustered with MIN rather than QIN (Li et al. 2003). These results suggested that the Jialing River is the distribution boundary between the Sichuan subspecies and the Qinling subspecies, rather than the Bailong River. The distribution boundary of the Mishmi subspecies and Bhutan subspecies is the Yalu Zangbu River, as suggested by the genomic evidence. Interestingly, the Yalu Zangbu River is also the boundary between the Himalayan red panda and the Chinese red panda (Hu et al. 2020;Joshi et al. 2021). The distribution boundary of the Himalayan takin and the Chinese takin is most likely the Three Parallel Rivers (Nujiang, Lancangjiang, and Jinshajiang). Because no recent occurrence of takin was reported from the east of Nujiang to the west of Jinshajiang, which river being the most important dispersal barrier requires further study with more samples across these regions (figs. 2a and 3b).
The Hengduan Mountains, situated at the southeastern edge of the Qinghai-Tibet Plateau, underwent dramatic geological and climatic changes during the Pleistocene epoch (Chen et al. 2010b;Li et al. 2016), which could have driven range shifts of allopatric populations, created new microhabitats, and thus accelerated diverse speciation (Meng et al. 2015;Li et al. 2016;Xing and Ree 2017;Dong et al. 2020;Rana et al. 2020). The complex topography provided a network of refugia and maintained high levels of genetic diversity among local animals and plants during the Pleistocene glaciations (Chen et al. 2010b;Sun et al. 2017;Hu et al. 2020Hu et al. , 2021. The divergence of the two takin species occurred at 0.23 Ma, corresponding well to the penultimate glaciation (0.30-0.13 Ma). This suggests that the population bottlenecks resulting from glaciations and subsequent geographical isolation may have played key roles in this species divergence. Genomic variation analysis showed that the GLG-SETE population has the highest genetic diversity ( fig.  5a-d; supplementary table S34, Supplementary Material online), suggesting that this region might be the glacial refugia of takins. PSMC analysis showed that the maximum Ne of the Himalayan takin was 1.6 times that of the Chinese takin, which partially explained why the Himalayan takin shows a higher genetic diversity than the Chinese takin.

Conclusions
We clearly elucidated the taxonomic relationships between and within species in takins, delineating two takin species and two subspecies of each species, which have important conservation implications. First, conservation management measures should be implemented separately for the two species to maintain their evolutionary uniqueness. Second, based on the genetic structure and distribution boundary, at least two management units (SETW and GLG-SETE) and at least three management units (QIO, MIN, and QIN) could be classified for the Himalayan and Chinese takins, respectively. Third, for the Chinese takin, the Qinling subspecies shows a much lower genomic diversity and higher LD, inbreeding, and genetic load levels than the Sichuan subspecies, and thus, it is in more urgent need of genetic management and conservation. For the Himalayan takin, the Bhutan subspecies with lower genomic variation and higher inbreeding needs more conservation attention. Last, to avoid captive interbreeding between the two species and preserve their unique genetic diversity, it is important to build a pedigree management system for captive takins, and independent captive populations should be maintained for the two species separately. The delimitation of takin species and subspecies and the determination of their boundaries have important implications for the biogeography of the Qinghai-Tibet Plateau and Hengduan Mountains and provide insights into the roles of Pleistocene glaciations and large rivers in mammal diversity patterns of this region.

Sample Collection and Sequencing
For takin genome sequencing, we collected fresh blood from a wild adult male takin from Louguantai, Shaanxi Province. We applied three sequencing technologies for de novo genome sequencing: Illumina HiSeq X10, Nanopore PromethION, and Hi-C. For Illumina HiSeq X10 sequencing, we constructed two paired-end libraries with insert sizes of 305-500 bp and generated clean data with a depth of 129× (335.58 Gb). For Nanopore PromethION sequencing, we constructed four long-read libraries and generated clean data with a depth of 100× (259.29 Gb). For Hi-C sequencing, we constructed four Hi-C libraries and generated clean data with a depth of 229× (592.62 Gb) (supplementary table S2

Comparative Genomic Analysis
The genomes of 20 Ruminantia species (Bos taurus, B. mutus, Aepyceros melampus, Cephalophus harveyi, Litocranius walleri, Eudorcas thomsonii, Saiga tatarica, Kobus ellipsiprymnus, Connochaetes taurinus, Oryx gazella, Pantholops hodgsonii, Oreamnos americanus, Hemitragus hylocrius, Ovis aries, Octoknema orientalis, Ammotragus lervia, Pseudois nayaur, Capra sibirica, C. hircus, and Tragulus kanchil) were downloaded from NCBI (supplementary table S10, Supplementary Material online). First, a pairwise alignment for each pair of genomes was generated by LAST (version 802) (Kielbasa et al. 2011) with the cattle genome as the reference. Each genome was aligned to the reference using the "lastal" command with the parameter -E 0.05. Then, we used the "maf-swap" command to change the order of the sequences in the MAF-format alignments and obtained the best pairwise aligned blocks. Finally, we used MULTIZ (version 11.2) (Blanchette et al. 2004) to merge the pairwise alignments into multiple genome alignments using the cattle genome as the reference.
Based on the generic feature format (GFF) file of cattle, we identified a total of 3.1 million syntenic 4-fold degenerate (4Dtv) sites from multiple genome alignments. All 4Dtv sites were combined into one sequence to construct a phylogenetic tree using RAxML (version 8.2.12) (Stamatakis 2014) with 1,000 bootstrap replicates in raxmlHPC-PTHREADS, followed by a search of the best-scoring maximum likelihood (ML) tree. Here, we used kanchil as an outgroup as it was positioned at the root of Ruminantia ). Based on this tree, we used r8s (Sanderson 2003) to estimate the divergence time with two fossil calibration points (supplementary table S40, Supplementary Material online).
To detect the potential genetic bases of the takin's unique morphological characteristics, we performed two groups of selection analyses based on the 4Dtv phylogenetic tree using both the branch and the branch-site models implemented in PAML: (1) takin was used as the foreground branch, and the other Caprinae species-including Oreamnos americanus, Hemitragus hylocrius, Ovis aries, Octoknema orientalis, Ammotragus lervia, Pseudois nayaur, Capra sibirica, and C. hircus-were used as the background branches; (2) takin was used as the foreground branch, and other Bovidae species were used as the background branches. First, we extracted orthologous genes from multiple genome alignments based on the GFF file of cattle using the scripts of Chen et al. (2019) (Supplementary Material online). After removing low-quality genes, we obtained 11,175 high-confidence orthologous genes, which were used in downstream analysis. For PSGs, the branchsite model was used. The likelihood ratio test (LRT) was conducted to compare a model that allowed sites to be under positive selection on the takin branch with the null model in which sites could evolve either neutrally or under purifying selection. P-values were computed based on χ 2 statistics, and genes with P , 0.05 were treated as candidate PSGs. For REGs, the branch model was used, with the null model (model = 0) assuming that all branches have been evolving at the same rate and the alternative model (model = 2) allowing the foreground and background branches to evolve under different rates. LRT with df = 1 was used to discriminate between alternative models for each ortholog in the gene set. Genes with P , 0.05 and a higher ω value on the foreground than the background branches were considered REGs. Both candidate PSGs and REGs were functionally enriched by the online tool Metascape (v1.0) , with a P-value of ,0.01.
Identification of Autosomal SNPs, X-chr SNPs, and Y-chr SNPs To obtain the corresponding contig sequences of autosomes, X and Y chromosomes in the takin genome assembly, we aligned the takin assembly (including 230 contigs) to the sheep chromosome-level genome (Li et al. 2021) (GCA_011170295) using MUMmer and identified 94 contigs from autosomes, 120 contigs from the X chromosome, and 16 contigs from the Y chromosome. We then extracted the autosomal SNPs, X-chr SNPs, and Y-chr SNPs according to the corresponding contig sequences. As our reference genome was from a male individual, we identified 37 males among the 75 individuals.

Mitochondrial Genome Assembly
First, we randomly extracted 40 million reads from raw paired-end data of 75 individuals. The Assembly by Reduced Complexity (ARC) (Hunter et al. 2015) approach, an assembly method based on a reference sequence, was then used to assemble the mitochondrial genome with the published takin mitochondrial genome as a reference (GenBank accession: KU361169). After successful assembly, we manually excluded the highly repetitive sequences within the D-loop region of each mitochondrial genome. The mitochondrial genome (MK748332) of Mishmi takin (B. t. taxicolor) downloaded from NCBI was also included in our mtgenome-based analyses.
Population Genetic Structure Based on Autosomal SNPs, X-chr SNPs, Y-chr SNPs, and Mitochondrial Genome To reduce required computing resources, we first performed LD-based pruning for the genotype data using PLINK (version 1.9) (Purcell et al. 2007) with the "-indep-pairwise 50 5 0.5" option for the autosome-SNP data set. For PCA, we used the Genome-wide Complex Trait Analysis (GCTA) tool (version 1.93.0) (Yang et al. 2011) to estimate the eigenvalues and eigenvectors. For admixture (genetic component sharing) analysis, we ran ADMIXTURE (version 1.3.0) (Alexander et al. 2009) with the default parameters, increasing K from 2 to 6, where K is the assumed number of genetic clusters. For phylogenetic analysis, neighboring joining trees of the autosomal SNPs, X-chr SNPs, Y-chr SNPs, and mitochondrial genomes were constructed by using MEGA (Kumar et al. 2018) with 1,000 bootstrap replicates using cattle/sheep as the outgroup. The NWK format tree files were further visualized using FigTree (version 1.4.4) (http://tree.bio.ed.ac.uk/ software/figtree/) and the online tree illustration website (https://itol.embl.de/). We also used NETWORK (version 10.1) (https://www.fluxus-engineering.com/) to construct a median-joining network map for the Y-chr SNP haplotypes and mitochondrial genome haplotypes.

Population Genetic Divergence (F ST )
We used a sliding-window method (50 kb window with 25 kb increment) to calculate the genome-wide distribution of F ST values between takin species, subspecies, and populations, implemented in VCFtools.

Species Delimitation
We used the SNAPP package (version 1.5.2) in BEAST (version 2.6.6) (Bouckaert et al. 2019) to run species delimitation analysis in takins. Five possible species delimitation models (supplementary table S24, Supplementary Material online) were carried out with 24 steps and each step was set to 100,000 generations with a pre-burnin of 10,000 generations. Then, we obtained the best model by comparing the marginal likelihood estimate values of five models. All Bayes factor (BF) calculations were made against the first taxonomy model (runA). Positive BF values indicate support for the first taxonomy model (runA), and BF . 10 is decisive.

Demographic, Divergence Histories, and Gene Flow of Species and Populations
To trace the demographic history of each takin species and subspecies/population, we employed PSMC to estimate changes in effective population size using heterozygous sites, with the following set of parameters: −N 30 -t 15 -r 5 −p 4 + 25*2 + 4 + 6. We selected two to three individuals sequenced at a high depth ( 30×) from each genetic population for PSMC analysis (supplementary  table S21, Supplementary Material online). The sex chromosome sequences of the takin genome were excluded. The generation time is often set as twice the sexual maturity age (Zhao et al. 2013). Here, the generation time (g) was set to 9 years as the sexual maturity age of takin was 4.5 years old (Wu et al. 1997). The nucleotide mutation rate (μ) of takin was estimated to be 1.11×10 −8 mutations per site per generation with cattle as the reference species for comparison using the following formula: μ = D × g/2T (Kondrashov and Crow 1993), where D is the observed frequency of pairwise differences between the two species, T is the estimated divergence time, and g is the estimated generation time for the takin.
To further infer species/population divergence histories, Fastsimcoal26, a flexible and robust simulation-based composite-likelihood approach, was implemented with the following parameters: −u −n 100000 −d −M 0.001 −l 100 −L 100 −q −C 10 −B 60 −c 10. To avoid the effects of four SETE individuals with potential hybrid signals, we removed the four individuals from the GLG-SETE population. We first converted the VCF file to a site frequency spectrum (SFS) file. Five alternative population divergence models were then simulated (supplementary fig. S6, Supplementary Material online). Three models (a, b, and c) Recent Speciation and Local Adaptation in Threatened Takins · https://doi.org/10.1093/molbev/msac111 MBE were tested to determine the time order between the ancestor node of SETW-GLG-SETE and the other ancestor nodes. The other two models (d and e) were tested to explore the divergence order (supplementary fig. S6, Supplementary Material online). Each model was run 100 times to ensure convergence, and the run with the highest likelihood was retained among the 100 runs. The best model was selected according to the ML value, and the parametric bootstrap estimates were obtained on the basis of 100 simulated data sets (supplementary table S41, Supplementary Material online).
We used D-statistics/ABBA-BABA statistics to test the gene flow among two takin species and four subspecies, implemented in ADMIXTOOLS (version 7.0) (Patterson et al. 2012 (Durand et al. 2011). We also used the Treemix method to detect the gene flow and its directions among the takin species and subspecies (Pickrell and Pritchard 2012).

Genomic Signatures of Selection and Local Adaptation between the Two Species
It is crucial to choose the corresponding selection methods for populations with different divergence time scales in a selection analysis (Sabeti et al. 2006). To detect selective signals between the Chinese and Himalayan takins, we chose two methods (XP-CLR [version 1.0] and F ST − θ π ) according to the divergence time ( 0.23 Ma) of the two species. We first scanned the genome for target regions under natural selection using the XP-CLR approach, which uses the AFS to calculate local deviations across populations. A sliding window of 0.5 cM with a spacing of 2 kb across the whole genome was scanned, and 200 SNPs were assayed in each window with the parameters -w1 0.005 200 2,000 chrN -p0 0.95. The genomic regions in the top 0.5% of the XP-CLR values were considered to represent selective sweeps.
In addition, we used a combined approach involving both the fixation (F ST ) and genetic diversity (θ π ) indexes to scan the whole genome to detect candidate selected regions because positive selection generally gives rise to lower genetic diversity within populations and higher genetic differentiation between populations. We used a slidingwindow method (50 kb windows with 25 kb increments) to calculate the genome-wide distribution of F ST values and θ π ratios for the two species, implemented in VCFtools. We applied Z transformation for F ST values and log 2 transformation for θ π ratios, and considered windows with the top 5% of Z(F ST ) and log 2 (θ π ratio) values, simultaneously, as the candidate outliers under strong selection. All the outlier windows were assigned to corresponding SNPs and genes. The functional enrichment analyses of candidate genes were implemented with the online tool Metascape (v1.0), with a P-value of ,0.01.

Estimation of Genetic Variations, LD, inbreeding, and genetic load
To calculate nucleotide diversity (θ π ) and Watterson's estimator (θ w ) of the whole-genome data in each subspecies, we used VCFtools (version 0.1.17) with the parameters "window-pi 50000 -window-pi-step 25000" to generate the suffix pi file, which was further used to calculate the final θ π and θ w, based on the official formula. We estimated the genetic diversity for mtgenome data using DNASP (version 5.10.01) (Librado and Rozas 2009). For genomic heterozygosity estimation, we calculated the ratio of the number of heterozygous SNPs to the genome size of the takin genome using an individual reference genome and a high sequencing-depth individual (SETW07) as representatives of the Chinese takin and the Himalayan takin. The genome-level LD of each subspecies was assessed using PopLDdecay ) with default settings. We used the genomic inbreeding coefficients (Froh) (Zhang et al. 2015) to evaluate the inbreeding level of each individual of each population: Froh = k L ROH k /L, where ROH is runs of homozygosity, L ROH k is the length of the kth ROH, and L is the genome size. The ROH value of each individual was calculated in PLINK with the parameters "-homozyg-kb 200 -homozyg-snp 50." For genetic load analysis, we first used SnpEff to annotate our filtered VCF files. Then, the coding sequence variants were classified as LoF (variants with the gain of a stop codon), missense, or synonymous variants. The missense variants were further classified as deleterious or tolerated missense mutations by SIFT 4G (Vaser et al. 2016). For the determination of ancestral and derived alleles, we first used LAST to compare the takin genome to the genomes of cattle and sheep, and obtained two derived site sets. Then, we obtained a total of 4.2 million derived allele sites from the intersection of the cattle and sheep data. To detect the genetic load of each subspecies, for each individual, we counted the numbers of nonhomozygous-and homozygous-derived alleles for LoF, deleterious and synonymous variants. Furthermore, we calculated the ratio of derived LoF variants (or deleterious variants) to that of derived synonymous variants at both the homozygous and nonhomozygous sites for each individual (Hu et al. 2020).

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.