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

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).

Some studies based on mitochondrial markers (e.g., Cytb, D-loop, and mitochondrial genome) have found that the takin is more closely related to Caprinae (Groves and Shields 1997; Feng et al. 2016; Yao et al. 2016; Kumar et al. 2019; Zhou et al. 2019). However, a genome-wide phylogenetic analysis is still lacking to confirm the phylogenetic status of takin. On the other hand, the takin has a much larger body size and a particular horn shape, different from the Caprinae species. Thus, its unique appearance has also aroused extensive interest from evolutionary biologists.

Furthermore, the intraspecies taxonomy of takin has long remained controversial as it has been classified into four subspecies or even four different species: (1) Mishmi takin B. t. taxicolor Hosdgson,1850, which is distributed in the southeastern Tibet and northwestern Yunnan of China, northern Myanmar, and northeastern India; (2) Bhutan takin B. t. whitei Lydekker, 1907, which is distributed in the southwest of the Yalu Zangbu River including the southern Tibet of China, northern Bhutan, and northern India; (3) Sichuan takin B. t. tibetana Milne-Edwards,1874, which lives in the western Sichuan and southern Gansu of China; and (4) Qinling takin B. t. bedfordi Thomas, 1911, which lives in the Qinling Mountains, Shaanxi, China. On the basis of differences in morphology (e.g., coat color and body size) and geographic distribution (supplementary table S1, Supplementary Material online), previous studies have proposed three possible taxonomic classifications within takins: (1) four subspecies, which was the mainstream viewpoint for a long time (Ellerman 1951; Wu 1986; Song et al. 2008); (2) three species (i.e., B. taxicolor, B. tibetana, and B. bedfordi) with two subspecies of B. taxicolor (B. t. taxicolor and B. t. whitei) (Lydekker 1915); and (3) four species, as proposed by Groves and Grubb in 2011, which updated the four subspecies into four independent species (Groves and Grubb 2011). Additionally, the distribution boundaries of these subspecies or species have remained undetermined. One main reason for these controversies is the lack of genetic evidence. To date, no study that includes samples from all four subspecies has been performed to explore the key problem (Li et al. 2003; Yao et al. 2016; Feng et al. 2017).

Here, we first sequenced and de novo assembled a high-quality chromosome-level genome of takin using a combined strategy involving Illumina, Nanopore, and Hi-C technologies. We resequenced the genomes of 75 wild takins from four subspecies covering most of its distribution range, with a mean depth of 18.4×. We aimed to clarify the species delimitation, demographic and divergence histories, adaptive evolution, and local adaptation for scientific conservation decisions and provide insights into the biogeography of the Himalaya–Hengduan biodiversity hotspot.

Results

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 (supplementary tables S5 and S6, Supplementary Material online). Repeat elements occupied 43.44% of the genome (supplementary table S7, Supplementary Material online). We identified 22,614 protein-coding genes (supplementary table S8, Supplementary Material online) and 22,145 (97.93%) of these genes were functionally annotated (supplementary table 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).

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.

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 species-specific 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 (GO:0008610, P = 8.31E−4), and inflammatory response (GO:0006954, P = 9.19E−5) (supplementary tables S18 and S19, Supplementary Material online), which might be related to takin’s unique appearance and local adaptation.

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 and supplementary table S20, Supplementary Material online). A total of ∼4.0 Tb of Illumina clean reads were generated (supplementary table S21, Supplementary Material online). Genome mapping showed an average mapping rate of 96.95% and an average genome depth of 18.43× (supplementary table S21, Supplementary Material online). Single nucleotide polymorphism (SNP) calling and stringent quality filtering obtained 7.27 million whole-genome SNPs (supplementary table S22, Supplementary Material online), including 6.92 million autosomal SNPs, 337.7 thousand X-chromosome (X-chr) SNPs, and 1,834 Y-chromosome (Y-chr) SNPs in 37 male individuals. Moreover, a total of 74 whole mitochondrial genome (mtgenome) sequences were successfully assembled from the genome resequencing reads of 74 individuals (supplementary table S20, Supplementary Material online).

Fig. 2.

Population genetic structure based on autosomal SNPs, Y-chr SNPs, and mitochondrial genomes of takin. (a) The geographic locations of wild takin samples under the background of distribution range from IUCN. QIN, Qinling population; MIN, Minshan population; QIO, Qionglai population; XXL, Xiaoxiangling population; GLG, Gaoligong population; SETE, Southeastern Tibet population to the east of the Yalu Zangbu River; SETW, Southeastern Tibet population to the west of the Yalu Zangbu River. (b) PCA result based on autosome SNPs. (c) ADMIXTURE results with K values 2–5 based on autosome SNPs. (d) NJ tree based on autosome SNPs with cattle as the outgroup. The values on the tree nodes indicate the bootstrap support of ≥60%. (e) Network map based on Y-chr SNP haplotypes. (f) Network map based on mitochondrial genome haplotypes.

On the basis of the autosomal SNPs, principal component analysis (PCA), ADMIXTURE, and phylogenetic tree analyses clearly revealed two divergent genetic lineages, GLG–SETE–SETW (referred to as the Himalayan lineage here) and QIN–MIN–QIO–XXL (referred to as the Chinese lineage here) (fig. 2b–d; supplementary fig. S2a and b, Supplementary Material online). In the Chinese lineage, the QIN population first diverged from the MIN–QIO–XXL populations, corresponding well to the two subspecies B. t. bedfordi and B. t. tibetana. The Himalayan lineage was found to have further diverged into two sublineages, GLG–SETE and SETW, corresponding to the two subspecies B. t. taxicolor and B. t. whitei (fig. 2b–d; supplementary fig. S2a, Supplementary Material online). Four SETE individuals (SETE8/9/10/11) may be the hybrids between SETW and GLG–SETE (fig. 2b–d). On the whole, this four sublineage classification was supported by the lowest cross-validation error rate of ADMIXTURE analysis (K = 4) (supplementary fig. S3, Supplementary Material online). Furthermore, the MIN population diverged from the QIO–XXL population (fig. 2c–d). In total, five genetic clusters were identified (fig. 2c–d). One individual, XXL01 from the XXL population, was genetically assigned to the QIO population (fig. 2b–d; supplementary fig. S2, Supplementary Material online). 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, FST measure based on whole-genome SNPs found that FST between the Himalayan and Chinese lineages was 0.2925, much higher than those between sublineages of the same lineage (GLG–SETE vs. SETW, FST = 0.2045; QIN vs. MIN–QIO–XXL, FST = 0.1613). FST 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 (FST = 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 Milne-Edwards, 1874. The Himalayan takin includes two subspecies: the Mishmi subspecies B. t. taxicolor Hodgson, 1850 (i.e., the GLG–SETE population), and the Bhutan subspecies B. t. whitei Lydekker, 1907 (i.e., SETW). The Chinese takin is also classified into two subspecies: the Sichuan subspecies B. t. tibetana Milne-Edwards, 1874 (i.e., MIN–QIO–XXL), and the Qinling subspecies B. t. bedfordi Thomas, 1911 (i.e., QIN) (fig. 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 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).

Fig. 3.

Demographic, divergence histories, and gene flow of takins. (a) PSMC results showed different demographic histories of the two species with a generation time (g) of 9 years and a mutation rate (μ) of 1.11× 10−8 per site per generation. The time axis is logarithmically transformed. (b) The divergence history of two takin species and their subspecies/populations was reconstructed by Fastsimcoal26. The time on the node indicated the divergence time, and the light blue dotted line and the river name represented the distribution boundary. The photos of four takin subspecies were taken in the wild and the photo credits are listed in the Acknowledgments. (c) Gene flow among the two takin species and four subspecies was estimated by Treemix.

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). A best-supported species/population divergence model was obtained by the comparison of five alternative models (fig. 3b; supplementary fig. S6 and table S25, Supplementary Material online). The divergence between the Himalayan takin (SETW and GLG–SETE) and the Chinese takin (QIN, MIN, and QIO) occurred at ∼0.23 Ma (fig. 3b), coincident with the split time (∼0.25 Ma) of the Ne curves of the two species, and was most likely caused by the penultimate glaciation. SETW diverged from GLG–SETE at ∼0.10 Ma (fig. 3b), possibly caused by geographical isolation imposed by the Yalu Zangbu River. QIN diverged from MIN–QIO at ∼56 ka, after which MIN separated from QIO at∼34 ka (fig. 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 (fig. 3b; supplementary table S25, Supplementary Material online).

To estimate the gene flow among two takin species and four subspecies, we performed the analyses of D-statistics/ABBA–BABA statistics. All the ABBA–BABA combinations of the two species and four subspecies detected significant signals of gene flow between the two species (D = 0.001334–0.009968, |Zscore| >3) (supplementary table S26, Supplementary Material online). Specifically, GLG–SETE (Mishmi subspecies) has more gene flow with MIN–QIO–XXL (Sichuan subspecies) than SETW (Bhutan subspecies) (D = 0.009968, |Zscore| = 24.317) (supplementary table S26, Supplementary Material online). In contrast, MIN–QIO–XXL has more gene flow with GLG–SETE than QIN (Qinling subspecies) (D = 0.001616, |Zscore| = 9.15) (supplementary table S26, Supplementary Material online). In addition, Treemix analyses detected the directions of gene flow. There were significant signals of gene flow from the ancestor of Chinese takin to GLG–SETE, and from GLG–SETE to MIN–QIO–XXL (fig. 3c).

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) (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 FST − θπ methods (fig. 4a–c; supplementary tables S28, S31, and S32, Supplementary Material online). Further examination of this gene identified three SNPs, including one missense variant in exon 6 (Chr03 100927526, G/A, Ala343Thr) (fig. 4d), one intron variant in intron 1, and one synonymous variant in exon 5 (supplementary table S33, Supplementary Material online). In particular, the allele frequencies of the missense variant were obviously different in the two species. The allele frequency of G was fixed in the Chinese takin populations, whereas the allele frequency of G was 0.477 in the Himalayan takin populations (fig. 4e; supplementary table S33, Supplementary Material online). PMEL is responsible for the formation of fibrillar sheets within the melanosome of melanocyte. The fibrillar sheets further serve as a template on which melanins polymerize after they are synthesized (Theos et al. 2005; Leonhardt et al. 2013; Watt et al. 2013; Bissig et al. 2016) (fig. 4f). Mutations in PMEL are associated with coat color dilution in cattle (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.

Fig. 4.

A coat color gene PMEL of takin under natural selection. (a) The XP-CLR values were calculated using XP-CLR with a nonoverlapping 25 kb window around the candidate locus (chr3: 99.55–101.55 Mb) including the PMEL gene (the black dotted square). (b and c) FST and nucleotide diversity (θ) were calculated using VCFtools for each 50 kb window with a 25 kb increment around the candidate locus. (d) The gene structure of PMEL includes 11 exons in takin. The missence mutation (Chr03 100927526 G/A Ala343Thr) was on the sixth exon. (e) Allele frequencies of the missence mutation in each takin population. (f) Role of the PMEL gene in melanosome formation, modified from reference (Watt et al. 2013).

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–d; supplementary table S34, Supplementary Material online). Among the four subspecies, the Qinling subspecies showed the lowest genetic diversity and the Mishmi subspecies harbored the highest genetic diversity (fig. 5a–d; supplementary table S34, Supplementary Material online). 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).

Fig. 5.

Genomic variation, linkage disequilibrium, inbreeding, and genetic load of takins. (a–d) Genetic diversity (θ and θw) of different takin species and subspecies based on the autosome-, X-chr-, Y-chr- SNPs, and mitochondrial genomes. (e) Linkage disequilibrium of the four subspecies. (f) Inbreeding level of the four subspecies measured by runs of homozygosity. (g) Ratios of homozygous- or nonhomozygous-derived deleterious variants to homozygous- or nonhomozygous-derived synonymous variants for different subspecies. (h) Ratios of homozygous- or nonhomozygous-derived LoF variants to homozygous- or nonhomozygous-derived synonymous variants for different subspecies.

Genomic LD (measured as R2) analysis indicated that the Qinling subspecies (QIN) presented the highest level of LD and slower LD decay than other subspecies, with R2 becoming stable at a distance of ∼300 kb, whereas faster LD decay was found in the other three subspecies, with a stable R2 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 protein-coding 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. 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 these genes might contribute to the development of large body size in takins.

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, FST, 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. 2020, 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.

Materials and Methods

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, Supplementary Material online). For the details of the estimation of genome size, genome assembly, assembly evaluation, genome annotation, and chromosome synteny analyses (supplementary tables S3–S9 and S37–S39, Supplementary Material online), see the Supplementary Materials.

For takin whole-genome resequencing, we collected skin and muscle samples of 75 dead wild takins from seven mountains or regions, including the Qinling Mountains (QIN, n = 20), Minshan Mountains (MIN, n = 17), Qionglai Mountains (QIO, n = 15), Xiaoxiangling Mountains (XXL, n = 1), Gaoligong Mountains (GLG, n = 4), the Southeast Tibet to the east of the Yalu Zangbu River (SETE, n = 11), and the Southeast Tibet to the west of the Yalu Zangbu River (SETW, n = 7) (fig. 2a; supplementary table S20, Supplementary Material online). The DNA of all samples was extracted using the DNeasy Blood & Tissue Kit (Qiagen, Valencia, CA, USA). Then, we constructed genomic libraries with insert sizes of 200–500 bp and performed genome resequencing at an average depth of 18.43× for each individual using the Illumina Nova-PE150 sequencing platform, resulting in the generation of ∼4.0 Tb of clean data (supplementary table S21, Supplementary Material online).

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 (Chen et al. 2019). 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 branch-site 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) (Zhou et al. 2019), with a P-value of <0.01.

Variant Calling and Filtering

For the raw whole-genome resequencing data of 75 takins, we first used fastp (Chen et al. 2018) to filter low-quality reads with the default parameters, and high-quality sequence reads were further mapped against the takin reference genome using BWA-mem (version 0.7.17) (Li and Durbin 2009) with the default parameters. We then used SAMtools (version 1.9) (Li et al. 2009) to sort BAM files and create index files. Potential PCR duplications were filtered using “MarkDuplicates” in Picard (version 2.21.6) (http://broadinstitute.github.io/picard). The Genome Analysis Toolkit (GATK version 4.1.2.0) (McKenna et al. 2010) was used for SNP calling with the analysis type of HaplotypeCaller with the “-ERC GVCF” option across 75 individuals. All the individual GVCF files were merged using “CombineGVCFs” in GATK. We called and selected candidate SNPs from the combined GVCF file using “GenotypeGVCFs” and “SelectVariants”, respectively. To obtain reliable candidate SNPs, hard filters were applied to the raw SNPs using GATK and VCFtools (version 0.1.17) (Danecek et al. 2011) according to the following criteria: QD <2.0, MQ <40.0, FS >60.0, SOR >3.0, MQRankSum <−12.5, ReadPosRankSum <−8.0, minDP <4, maf <0.05 and max-missing >0.2. After filtering, the remaining high-quality SNPs were annotated according to the GFF file of the takin reference genome using SnpEff (version 4.3) (Cingolani et al. 2014).

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 (FST)

We used a sliding-window method (50 kb window with 25 kb increment) to calculate the genome-wide distribution of FST 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) 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). For [[[P1, P2], P3], O], we performed multiple tests (Supplementary table S26, Supplementary Material online), with sheep as the outgroup (O). Significant signal of gene flow occurs between P2 and P3 with a D-value of >0 and |Zscore| >3, and between P1 and P3 with D < 0 and |Zscore| >3 (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 FST − θ) 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 (FST) 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 sliding-window method (50 kb windows with 25 kb increments) to calculate the genome-wide distribution of FST values and θ ratios for the two species, implemented in VCFtools. We applied Z transformation for FST values and log2 transformation for θ ratios, and considered windows with the top 5% of Z(FST) and log2(θ 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 (Zhang et al. 2019) 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: , where ROH is runs of homozygosity, 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.

Click here for additional data file.