Whole-genome sequencing of Tahe red deer (Cervus hanglu yarkandensis) reveals genetic diversity and selection signatures

The Tahe red deer (TRD), domesticated and artificially raised from wild Tarim red deer, is valued for its high-quality antlers and ability to survive tough desert conditions. Nowadays, the decline in the population of TRD has significantly impacted their genetic diversity, posing a serious threat to their conservation and utilization. However, information based on whole-genome sequencing data of TRD is scarce, and the mechanisms underlying adaptive characteristics remain poorly understood. Additionally, research on Tahe red deer holds great importance for elucidating the evolutionary history and adaptability of the genus Cervus. This study aimed to investigate the genetic diversity, population structure, and selection signals of TRD using whole-genome sequencing data. The results revealed that TRD exhibited high inbreeding level and relatively low genetic diversity, and that TRD had a closer relationship with Cervus canadensis. Using three methods, including the fixation index, nucleotide diversity and cross-population extended haplotype homozygosity, there were 573 genes annotated in 2,303 overlapping candidate selection regions such as IL1R1, F13B, ARHGAP15, DCLK3, CACHD1, NDEL1, and UPF1, most of which were associated with adaptation to a hot arid environment. In summary, this study offered genomic markers and candidate genes associated with these traits, providing valuable insights for improving future breeding strategies of TRD.

1 Introduction

Adaptation represents a continuous and long-term evolutionary process, during which beneficial alleles accumulate at gene loci (1). Species within the genus Cervus have evolved distinct survival strategies and physiological adaptations over time, shaped by various environmental factors such as climate and habitat types. Understanding the mechanisms by which populations or species respond to environmental changes is essential for the conservation of key species under the pressure of global climate change (2).

The genus Cervus with its diverse populations across China has long been extensively farmed for a variety of valuable products such as antler velvet and venison. The Tarim red deer found in the Tarim River basin, as an essential component of the local ecosystem, exhibits considerable biological adaptability under harsh desert environment (3). In recent decades, the fragmentation of the wild Tarim red deer’s habitat, caused by increased human activities and climate change, has led to a further decline in the wild population (4). Accordingly, China has classified the wild Tarim red deer as a first-class nationally protected wild animal and implemented a series of conservation measures to protect this endangered species. The Tahe red deer (TRD) was domesticated and artificially raised from the wild Tarim red deer (Cervus hanglu yarkandensis) by the Second Agricultural Division of Xinjiang Construction Corps since the 1950s (5). As an artificially bred population with significant economic and ecological value, TRD is highly regarded for its high antler yield and remarkable ability to cope with arid conditions. Due to a lack of awareness and the impact of market forces, the number of TRD has witnessed a sharp decline, with the current stock being less than 10,000, posing a serious threat to their genetic diversity (6).

Currently, studies on the genetic diversity of TRD mainly focuses on the sequences from Y chromosome and mitochondrion, and the results indicates an imbalance with high haplotype diversity and low nucleotide diversity (6). Analyses on the phylogenetic relationships about the Tarim red deer have primarily relied on mitochondrial genomes and markers (6–8). Mitochondrial control region analyses have identified two distinct lineages: the Western lineage including red deer (Cervus elaphus) and the Tarim red deer, and the Eastern lineage comprising wapiti (Cervus canadensis) and sika deer (Cervus nippon) (7). However, studies based on cytochrome b gene and D-loop region have identified the red deer, the Tarim red deer, and wapiti as three individual species (6–9). Studies based on mitochondrial DNA and microsatellite markers have not reached a concordant conclusion on the phylogeny of Cervus. Compared with mitochondrial DNA and microsatellite markers, nuclear DNA data can more accurately reveal genetic differences between breeds, making it widely used to evaluate genetic diversity and detect selective sweeps (10). Undoubtedly, the assessment of genetic diversity and population structure based on genomic information will be essential for the better breeding management and sustainable use of genetic resources.

TRD has unique adaptability that makes it an ideal subject for studying how animals cope with extreme conditions (11). The antler velvet with high yield serves as a vital ingredient in traditional Chinese medicines and health products (6). Given that the valuable genetic resource is severely depleted, rational conservation and utilization are urgent. In addition, the study of TRD is of great significance for understanding the evolution and adaptability of the genus Cervus. Presently, information on genomic diversity of TRD based on the whole genomic sequences is scarce and the publicly available genomic data remain limited. This study assessed the genetic diversity in TRD and identified candidate genes associated with adaptation to hot arid environments using whole-genome data from 17 TRD individuals, thereby providing valuable insights into important genetic variations and the conservation of genetic resource.

2 Materials and methods

2.1 Sample selection

Antler slices of 17 TRD were purchased as the commercial products from the farm located in Korla City, Xinjiang Uygur Autonomous Region of China (Supplementary Figure S1). Genomic DNA was extracted for WGS using the EasyPure Genomic DNA Kit (TransGen Biotech). DNBSEQ-T7 was used at the Novogene Bioinformatics Institute company (Beijing, China) to generate 2 × 150 bp paired-end read data for each individual (Supplementary Table S1). Additionally, to enhance our understanding of the genomic genetic diversity and selection signals in TRD, WGS data for 27 published individuals including Cervus elaphus hispanicus (ERD, n = 14) and Cervus canadensis nelsoni (CRD, n = 13) were acquired from the Sequence Read Archive.¹

2.2 Reads mapping and quality control

Raw data were filtered using FASTP v 0.20.1 software (12). Burrows-Wheeler Aligner (BWA) software (v0.7.13) was used for mapping all cleaned data to the red deer reference genome (GCF_910594005.1, mCerEla1.1) using “bwa mem” parameters (13). The aligned BAM files were sorted using SAMtools v1.19 (14). Subsequently, the Picard MarkDuplicates tool (v1.115)² was employed to mark duplicate reads from each alignment. After sorting reads, variants were called using the Genome Analysis Toolkit (GATK v4.1.4) and filtered with GATK’s “VariantFiltration” module (15). The hard filters were applied to the raw SNPs according to the criteria as follows: QD < 2.0 || FS > 60.0 || MQ < 40.0 || SOR > 3.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0. Quality control was then performed using PLINK (v1.9) with the parameters “--geno 0.05 --mind 0.1,” ensuring that each locus retained a minimum of 3 alleles, and biallelic loci on autosomes were acquired using BCFtools (v1.8) (14, 16). The remaining SNPs were subsequently used for further analysis. Based on the mCerEla1.1 reference genome annotation file, the variants after quality control in TRD were annotated using SnpEff software (v5.1d) (17).

2.3 Genetic diversity analyses

To assess genetic diversity, expected heterozygosity (H_E), observed heterozygosity (H_O), nucleotide diversity (pi) and the runs of homozygosity (ROH) were estimated. PLINK (v1.9) software was used to calculate H_O and H_E with the option “—hardy” (16). VCFtools software (v0.1.16) was employed to assess pi with the parameters “--window-pi 20,000 --window-pi-step 10,000” (18). To evaluate the inbreeding degree, the number and length of ROH fragments were calculated using VCFtools with the “—LROH” command (18). After that, the runs of homozygosity-based inbreeding coefficient (F_ROH) were calculated as the total length of ROH fragments divided by the length of the autosomes. PopLDdecay (v3.42) was applied to assess the degree of linkage disequilibrium decay (LD) of each population by calculating the pairwise SNP correlation coefficients (r²) with the parameter “-MaxDist 1,000” (19).

2.4 Population structure analysis

The NJ method was applied to construct a tree based on Nei genetic distance using VCF2Dis (v1.47) and visualized with the Splitstree software (20–22). Sites with high linkage disequilibrium (LD) were removed using PLINK with the parameter “-indep-pairwise 50 25 0.2” and the remaining SNPs were used for the following population structure analysis (16). Principal component analysis (PCA) was performed using GCTA software (v1.92.3) with the parameter “grm” and plotted in the R package ggplot2 (23, 24). Admixture analyses were conducted using ADMIXTURE software (v1.3.0) with the parameters “admixture -cv” to infer ancestral populations in our dataset (25). Cross-validation error values for clustering at K = 2 to 4 were also calculated. BEAGLE (v5.4) software was employed to first phase the data file and impute missing genotypes (26, 27). Haplotypes were then inferred with Refine IBD, using a 40-SNP sliding window, a 1.5 cM minimum haplotype length, a 0.15 cM trimming threshold, and a LOD score of 1 (28). The counts of haplotypes shared between populations were visualized with the circlize package in R (29).

2.5 Genome-wide scanning for selection signatures

To identify the candidate regions and genes associated with the unique traits of TRD compared to CRD and ERD, three analytical methods were performed: the fixation index (F_ST), pi and cross-population extended haplotype homozygosity (XP-EHH). Selection signal values were computed based on a 20 kb sliding window. The pi and F_ST values were calculated using VCFtools (v0.1.16) with the parameters “--window-pi 20,000 --window-pi-step 10,000” and “--fst-window-size 20,000 --fst-window-step 10,000,” respectively (18). Additionally, haplotype phasing was conducted using BEAGLE (v5.4), and the phased data were then used to estimate XP-EHH (30). To detect the positive selection signatures between TRD and other populations, XP-EHH was calculated using selscan (v1.1), and the results were subsequently normalized using the “-norm” module of selscan for each 20 kb region (31). To reduce false positives, only the top 5% genomic windows from the three methods were selected for further analysis. The genomic regions identified by all three methods were intersected using BEDtools (v2.30.0) with the “intersect” parameters to determine the potential candidate regions of selection (32). Functional annotations of the identified regions were carried out using SnpEff software (v5.1d). To gain a better understanding of the identified candidate genes, Gene Ontology (GO) and Kyoto Encyclopedia of Gene and Genomes (KEGG) pathway enrichment analyses were performed using KOBAS, with cow (Bos taurus) genome selected as ortholog (33). Only when the p-value was less than 0.05 were the pathways considered to be significantly enriched.

3 Results and discussion

3.1 SNP genotyping and annotation of TRD

In total, 675.51 Gb of raw data from 17 TRD individuals was obtained through the whole-genome resequencing (Supplementary Table S1). After the filters, 671.74 Gb of clean data were retained. The sequencing quality metrics indicated an average Q20 of 99.49%, an average Q30 of 98.08%, a GC content of 43.37%, mapping rates between 99.23 and 99.88%, and an average sequencing depth of 13.44×. Furthermore, 12,086,946 high-quality autosomal bi-allelic SNPs were detected in the SNP dataset of 17 TRD after quality control, of which no individuals were excluded based on the “mind 0.1” criterion. Among the identified variants, there were 8,819,706 transition SNPs (T_S) and 3,267,240 transversion SNPs (T_V), leading to a T_S/T_V ratio of 2.70. The least number of SNPs was observed in chromosome 26, while the highest number of SNPs was found in chromosome 30, with an average variants frequency of 4.47 SNPs/Kb (Supplementary Table S2). A detailed SnpEff annotation showed that most SNPs were found in intergenic regions (47.26%) and intronic regions (38.42%), with only 1.01% located in exonic regions (Supplementary Table S3). Among the 121,553 exonic SNPs identified, there were 66,662 synonymous variants, 54,167 missense variants, and 724 nonsense variants.

3.2 Genetic diversity analysis

The average values of H_O, H_E, ROH and pi were computed to evaluate the genetic diversity of all populations (Supplementary Table S4). The results indicated that the pi value ranged from 0.0013 to 0.0020. Among them, TRD (0.0016) was higher than CRD (0.0013) but lower than ERD (0.0020). Besides, the H_O and H_E of TRD was between CRD and ERD (Figure 1A). Consistent with the research on the genetic diversity of TRD based on mitochondrial sequences and the Y chromosome, the results also indicated a low level of nucleotide diversity (6). The F_ROH value of TRD was 0.2729, which was much higher compared to ERD. As for LD decay, the r² values declined sharply as genomic distance increased in all populations, with the most rapid decline occurring within the first 50 kb (Figure 1B). This indicated that there had been a significant amount of recombination in the recent history of TRD. When the distance between markers exceeded the first 50 kb, the results indicated that TRD had a relatively high LD level, whereas ERD had the lowest. These values suggested that TRD may have relatively low genetic polymorphism and inbreeding within the population, which reminds us that rational conservation measures and breeding plans should be established to conserve the genetic resource.

Figure 1

Figure 1. Genetic diversity of 3 populations. (A) Distribution of expected heterozygosity (H_E) and observed heterozygosity (H_O) in each population. (B) Decay of linkage disequilibrium (LD) on autosomes was estimated from 3 populations.

3.3 Population structure analysis

To investigate the genetic distances among different populations, phylogenetic analysis was performed on 44 individuals representing three populations. It can be observed from the Neighbor-Net network that three populations formed clearly distinct clusters from one another and TRD was more closely related to CRD than to ERD, which was different from the previous research finding on the basis of comparative genomics (Figure 2B). The 2,506,148 SNPs after linkage pruning were used for PCA and admixture analysis to further investigate the cluster patterns between TRD with other populations (Figure 2A). The lowest cross-validation error value was observed at K = 3, as the hypothetical ancestral groups were evaluated for K values between 2 and 3. When K = 2, these different populations can be genetically divided into two groups: one group consisted of ERD, while the other group included CRD and TRD. When K = 3, TRD showed distinct admixture component proportions compared to the CRD. In PCA analysis, the first principal component (PC1) accounted for 19.02% of the variation in the genomic data, while the second principal component (PC2) explained 15.49% (Figure 2C). The PCA result showed that TRD exhibited a closer genetic relationship with CRD compared to ERD, which was consistent with Neighbor-Net network and ADMIXTURE results. The haplotype exchange analysis between populations identifies more pronounced haplotype exchange between CRD and TRD, while revealing reduced IBD sharing between TRD and ERD (Figure 2D). Moreover, TRD was genetically closer to CRD than to ERD. In addition, mitochondrial phylogenetic tree indicated that TRD had a closer genetic relationship with ERD than CRD. Phylogenetic analysis based on comparative genomics and mitogenomes were inconsistent with the nuclear DNA phylogeny, probably due to incomplete lineage sorting or genetic introgression (10, 34). Whole genomes would provide more accurate evidence in this respect because nuclear DNA data contains a greater quantity of genetic markers. One might interpret it as the Tarim red deer differentiated relatively early due to either an extreme cold climate or human activities, evolved independently in an isolated environment, and had a distant genetic relationship with other CRD populations (4). Given the individuals selected for this study were domesticated, the results of population structure and phylogenetic relationships may be influenced by the size and source of the sampled populations. Therefore, samples from a larger number of TRD and wild Tarim red deer should be included to verify whether this pattern remains consistent.

Figure 2

Figure 2. Population structure and relationships of TRD with other populations. (A) Result of admixture (K = 2 and 3). (B) Neighbor network constructed from Nei. (C) Principal component analysis (PCA). (D) Circos plot of the number of shared haplotypes.

3.4 Selection signature analysis

Two wild populations were selected as the reference populations, including the European red deer (ERD), which primarily inhabit mountain forests with a Mediterranean climate, and the North American Elk (CRD), which live in forests dominated by pine trees with distinct seasonal changes. The Tarim red deer have evolved to cope with the heat and arid desert environment, but the mechanisms underlying their adaptive features remain poorly understood. In the selection signature analyses, there were 11,673 regions, 12,091 regions and 6,441 regions ascertained through pi, F_ST and XP-EHH methods in the top 5%, respectively. For example, CACHD1, CPLX1, DIS3L2, MAGI2, MYH10, NCAM2, NDEL1, POC1A, SUPT3H, TWF2, and UPF1, which have high selective sweeps identified by both pi and F_ST methods in TRD, have also been identified as candidate genes in sheep breeds adapting to desert environments (35). Only the top 5% of windows identified by all three methods were used to pinpoint candidate regions and genes. In total, 2,303 candidate selection regions were identified including 573 genes likely subject to strong selection in TRD. The detailed information regarding these regions was depicted in Figure 3 and documented in Supplementary Tables S5–S8. Study indicates that ARHGAP15, which regulates various biological processes such as cytoskeletal dynamics and cell movement, is linked to physiological traits related to tropical adaptation in Zaobei Beef cattle (Table 1) (36). F13B, a gene strongly selected in TRD, has also been pinpointed as a candidate gene for the desert adaptation of Bactrian camels (37). Meanwhile, DCLK3 (doublecortin like kinase 3), characterized by its biased expression in testis and brain, might be involved in the heat stress and essential for improving thermoregulation in sika deer (38). Several candidate genes, including MYO18A, FCER1G, IL10RB, CCR8, and CXCR6, are connected to the immune system and environmental adaptation, which might be useful for the survival of the Tarim red deer in the Tarim Basin (39–42). ANKH regulates bone formation and skeletal development by inhibiting the mineralization process through the encoding of an inorganic pyrophosphate transport regulator (43). Additionally, we noted several genes associated with meat quality (RPL29, ANK1, and GSK3A), production and growth traits (CDK6, ZBTB7C, KIF22, ATP6V1H, UBR2, and CUBN) (44–49).

Figure 3

Figure 3. Manhattan plots of selection signature between TRD and other populations. The dashed line is the 5 percent threshold. (A) The fixation index (F_ST). (B) Nucleotide diversity (pi). (C) Cross-population extended haplotype homozygosity (XP-EHH). (D) Signal intersection analysis of pi, F_ST, and XP-EHH.

Table 1

Table 1. Potential selected genes associated with important traits identified by three methods in TRD.

The functional enrichment analysis for 573 overlapped genes were further analyzed. The significantly enriched GO terms and KEGG pathways are shown in Supplementary Table S9. Among them, a total of 267 significantly enriched GO terms with p-value < 0.05 were observed, such as heat shock protein binding (GO:0031072, p = 0.015360), cellular response to oxidative stress (GO:0034599, p = 0.004995), response to bacterium (GO:0009617, p = 0.021797), regulation of neuron projection development (GO:0010975, p = 0.004311), and positive regulation of phagocytosis (GO:0050766, p = 0.011622), which may play role in desert adaptation. Several genes, including FCER1G, ALCAM, CCR8, CXCR6, CNTFR, TNFRSF13C, and IL1R1, were significantly enriched in the external side of plasma membrane, which is key for environmental signal perception and immune cell communication (GO:0050766, p = 0.037169). In addition, metabolism-related and growth-related biological functions were significant, such as regulation of bone mineralization (GO:0030500, p = 0.026519), Wnt signaling pathway (GO:0016055, p = 0.000777), activation of protein kinase B activity (GO:0032148, p = 0.021715), ATP binding (GO:0005524, p = 0.000015), and microtubule-based movement (GO:0007018, p = 0.048545). ANKH and SGMS2 were both enriched in the regulation of bone mineralization (GO:0030500, p = 0.026519), suggesting their potential involvement in antler formation and skeletal development in red deer. MT3, which is involved in activation of protein kinase B activity (GO:0032148, p = 0.021714), exerts a key influence on osteoclastogenesis and osteoporosis via dual routes involving reactive oxygen species and specificity protein 1 (50). Furthermore, 45 significant enriched pathways were obtained, including MAPK signaling pathway (bta04010, p = 0.001910), cell cycle (bta04110, p = 0.006348), cytokine-cytokine receptor interaction (bta04060, p = 0.027605), and Jak–STAT signaling pathway (bta04630, p = 0.005778), which are related to immune response and skeletal development. The enrichment of interleukin-1 receptor activity (GO:0004908, p = 0.003269) and MAPK signaling pathway (bta04010, p = 0.001910) further supports the role of IL1R1 as a vital mediator in many immune and inflammatory responses induced by cytokines, which has been reported to be involved in the heat stress response (51). Overall, the results suggest that TRD has experienced strong selection on genes related to immunity, skeletal development, and stress resistance, reflecting a genetic adaptation to arid desert environments.

4 Conclusion

The study showed that TRD population has relatively low genetic diversity and a high level of inbreeding. The results also validate the classification of TRD, ERD, and CRD into three distinct branches, with TRD being more closely related to CRD. Additionally, based on whole-genomic data, candidate genes related to adaptation to heat environments were identified and annotated. These findings will be very useful for the future conservation and management of TRD and determining the potential genomic mechanisms in harsh environments.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Ethics statement

All experiments were conducted in accordance with the Institutional Animal Care and Use Committee guidelines (IACUC No. SY202306046) under currently approved protocols at Jilin University. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

TP: Validation, Writing – review & editing, Writing – original draft, Conceptualization. WY: Conceptualization, Writing – review & editing, Validation. ZM: Writing – review & editing, Project administration, Data curation, Resources. ZW: Data curation, Writing – review & editing, Formal analysis. HS: Visualization, Writing – review & editing. SY: Writing – review & editing, Funding acquisition, Project administration, Supervision.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This study was supported by the Science and Technology Development Project of Jilin Province, China (no. 20240303071NC).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Gen AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2025.1642382/full#supplementary-material

Supplementary Figure S1 | The picture of an adult Tahe red deer.

Footnotes

1. ^https://www.ncbi.nlm.nih.gov/sra/

2. ^https://github.com/broadinstitute/picard

References

1. Barrett, RDH, and Schluter, D. Adaptation from standing genetic variation. Trends Ecol Evol. (2008) 23:38–44. doi: 10.1016/j.tree.2007.09.008

PubMed Abstract | Crossref Full Text | Google Scholar

2. Bradshaw, WE, and Holzapfel, CM. Climate change. Evolutionary response to rapid climate change. Science. (2006) 312:1477–8. doi: 10.1126/science.1127000

PubMed Abstract | Crossref Full Text | Google Scholar

3. Ababaikeri, B, Abduriyim, S, Tohetahong, Y, Mamat, T, Ahmat, A, and Halik, M. Whole-genome sequencing of Tarim red deer (Cervus elaphus yarkandensis) reveals demographic history and adaptations to an arid-desert environment. Front Zool. (2020) 17:31. doi: 10.1186/s12983-020-00379-5

PubMed Abstract | Crossref Full Text | Google Scholar

4. Mahmut, H, Ganzorig, S, Onuma, M, Masuda, R, Suzuki, M, and Ohtaishi, N. A preliminary study of the genetic diversity of Xinjiang Tarim red deer (Cervus elaphus yarkandensis) using the microsatellite DNA method. Jpn J Vet Res. (2001) 49:231–7.

Google Scholar

5. Cao, Z. Techniques for harvesting antlers of Tahe Red Deer. Rur Sci Technol. (2022):69–70. doi: 10.19777/j.cnki.issn1002-6193.2022.01.023

Crossref Full Text | Google Scholar

6. Yang, S, Wang, H, Dong, Y, Zhao, X, Dong, S, Wang, X, et al. Analysis of the genetic diversity of Tahe red deer and population structure of Cervus elaphu/hanglu/canadensis. Gene. (2025) 933:148960. doi: 10.1016/j.gene.2024.148960

PubMed Abstract | Crossref Full Text | Google Scholar

7. Doan, K, Mackiewicz, P, Sandoval-Castellanos, E, Stefaniak, K, Ridush, B, Dalén, L, et al. The history of Crimean red deer population and Cervus phylogeography in Eurasia. Zool J Linnean Soc. (2018) 183:208–25. doi: 10.1093/zoolinnean/zlx065

Crossref Full Text | Google Scholar

8. Polziehn, RO, Hamr, J, Mallory, FF, and Strobeck, C. Microsatellite analysis of north American wapiti (Cervus elaphus) populations. Mol Ecol. (2000) 9:1561–76. doi: 10.1046/j.1365-294x.2000.01033.x

PubMed Abstract | Crossref Full Text | Google Scholar

9. Ludt, CJ, Schroeder, W, Rottmann, O, and Kuehn, R. Mitochondrial DNA phylogeography of red deer (Cervus elaphus). Mol Phylogenet Evol. (2004) 31:1064–83. doi: 10.1016/j.ympev.2003.10.003

PubMed Abstract | Crossref Full Text | Google Scholar

10. Mackiewicz, P, Matosiuk, M, Świsłocka, M, Zachos, FE, Hajji, GM, Saveljev, AP, et al. Phylogeny and evolution of the genus Cervus (Cervidae, Mammalia) as revealed by complete mitochondrial genomes. Sci Rep. (2022) 12:16381. doi: 10.1038/s41598-022-20763-x

PubMed Abstract | Crossref Full Text | Google Scholar

11. Skelly, DK, Joseph, LN, Possingham, HP, Freidenburg, LK, Farrugia, TJ, Kinnison, MT, et al. Evolutionary responses to climate change. Conserv Biol. (2007) 21:1353–5. doi: 10.1111/j.1523-1739.2007.00764.x

PubMed Abstract | Crossref Full Text | Google Scholar

12. Chen, S, Zhou, Y, Chen, Y, and Gu, J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. (2018) 34:i884–90. doi: 10.1093/bioinformatics/bty560

PubMed Abstract | Crossref Full Text | Google Scholar

13. Li, H, and Durbin, R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. (2009) 25:1754–60. doi: 10.1093/bioinformatics/btp324

PubMed Abstract | Crossref Full Text | Google Scholar

14. Danecek, P, Bonfield, JK, Liddle, J, Marshall, J, Ohan, V, Pollard, MO, et al. Twelve years of SAMtools and BCFtools. Gigascience. (2021) 10:giab008. doi: 10.1093/gigascience/giab008

PubMed Abstract | Crossref Full Text | Google Scholar

15. McKenna, A, Hanna, M, Banks, E, Sivachenko, A, Cibulskis, K, Kernytsky, A, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. (2010) 20:1297–303. doi: 10.1101/gr.107524.110

PubMed Abstract | Crossref Full Text | Google Scholar

16. Purcell, S, Neale, B, Todd-Brown, K, Thomas, L, Ferreira, MAR, Bender, D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. (2007) 81:559–75. doi: 10.1086/519795

PubMed Abstract | Crossref Full Text | Google Scholar

17. Cingolani, P, Platts, A, Wang, LL, Coon, M, Nguyen, T, Wang, L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly. (2012) 6:80–92. doi: 10.4161/fly.19695

PubMed Abstract | Crossref Full Text | Google Scholar

18. Danecek, P, Auton, A, Abecasis, G, Albers, CA, Banks, E, DePristo, MA, et al. The variant call format and VCFtools. Bioinformatics. (2011) 27:2156–8. doi: 10.1093/bioinformatics/btr330

PubMed Abstract | Crossref Full Text | Google Scholar

19. Zhang, C, Dong, S-S, Xu, J-Y, He, W-M, and Yang, T-L. PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics. (2019) 35:1786–8. doi: 10.1093/bioinformatics/bty875

PubMed Abstract | Crossref Full Text | Google Scholar

20. Huson, DH, and Bryant, D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. (2006) 23:254–67. doi: 10.1093/molbev/msj030

PubMed Abstract | Crossref Full Text | Google Scholar

21. Nei, M. Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics. (1978) 89:583–90. doi: 10.1093/genetics/89.3.583

PubMed Abstract | Crossref Full Text | Google Scholar

22. Xu, L, He, W, Tai, S, Huang, X, Qin, M, Liao, X, et al. VCF2Dis: an ultra-fast and efficient tool to calculate pairwise genetic distance and construct population phylogeny from VCF files. GigaScience. (2025) 14:giaf032. doi: 10.1093/gigascience/giaf032

PubMed Abstract | Crossref Full Text | Google Scholar

23. Ginestet, C. Ggplot2: elegant graphics for data analysis. J R Stat Soc Ser A Stat Soc. (2011) 174:245–6. doi: 10.1111/j.1467-985X.2010.00676_9.x

Crossref Full Text | Google Scholar

24. Yang, J, Lee, SH, Goddard, ME, and Visscher, PM. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet. (2011) 88:76–82. doi: 10.1016/j.ajhg.2010.11.011

PubMed Abstract | Crossref Full Text | Google Scholar

25. Alexander, DH, Novembre, J, and Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. (2009) 19:1655–64. doi: 10.1101/gr.094052.109

PubMed Abstract | Crossref Full Text | Google Scholar

26. Browning, BL, Zhou, Y, and Browning, SR. A one-penny imputed genome from next-generation reference panels. Am J Hum Genet. (2018) 103:338–48. doi: 10.1016/j.ajhg.2018.07.015

PubMed Abstract | Crossref Full Text | Google Scholar

27. Browning, BL, Tian, X, Zhou, Y, and Browning, SR. Fast two-stage phasing of large-scale sequence data. Am J Hum Genet. (2021) 108:1880–90. doi: 10.1016/j.ajhg.2021.08.005

PubMed Abstract | Crossref Full Text | Google Scholar

28. Browning, BL, and Browning, SR. Improving the accuracy and efficiency of identity-by-descent detection in population data. Genetics. (2013) 194:459–71. doi: 10.1534/genetics.113.150029

PubMed Abstract | Crossref Full Text | Google Scholar

29. Gu, Z, Gu, L, Eils, R, Schlesner, M, and Brors, B. Circlize implements and enhances circular visualization in R. Bioinformatics. (2014) 30:2811–2. doi: 10.1093/bioinformatics/btu393

PubMed Abstract | Crossref Full Text | Google Scholar

30. Browning, SR, and Browning, BL. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. Am J Hum Genet. (2007) 81:1084–97. doi: 10.1086/521987

PubMed Abstract | Crossref Full Text | Google Scholar

31. Szpiech, ZA, and Hernandez, RD. Selscan: an efficient multithreaded program to perform EHH-based scans for positive selection. Mol Biol Evol. (2014) 31:2824–7. doi: 10.1093/molbev/msu211

PubMed Abstract | Crossref Full Text | Google Scholar

32. Quinlan, AR. BEDtools: the swiss-army tool for genome feature analysis. Curr Protoc Bioinformatics. (2014) 47:11.12.1–34. doi: 10.1002/0471250953.bi1112s47

Crossref Full Text | Google Scholar

33. Bu, D, Luo, H, Huo, P, Wang, Z, Zhang, S, He, Z, et al. KOBAS-i: intelligent prioritization and exploratory visualization of biological functions for gene enrichment analysis. Nucleic Acids Res. (2021) 49:W317–25. doi: 10.1093/nar/gkab447

PubMed Abstract | Crossref Full Text | Google Scholar

34. Tang, L, Dong, S, and Xing, X. Comparative genomics reveal phylogenetic relationship and chromosomal evolutionary events of eight Cervidae species. Animals. (2024) 14:1063. doi: 10.3390/ani14071063

PubMed Abstract | Crossref Full Text | Google Scholar

35. Yang, J, Li, W-R, Lv, F-H, He, S-G, Tian, S-L, Peng, W-F, et al. Whole-genome sequencing of native sheep provides insights into rapid adaptations to extreme environments. Mol Biol Evol. (2016) 33:2576–92. doi: 10.1093/molbev/msw129

PubMed Abstract | Crossref Full Text | Google Scholar

36. Shi, L, Zhang, P, Liu, Q, Liu, C, Cheng, L, Yu, B, et al. Genome-wide analysis of genetic diversity and selection signatures in Zaobei beef cattle. Animals (Basel). (2024) 14:2447. doi: 10.3390/ani14162447

PubMed Abstract | Crossref Full Text | Google Scholar

37. Wu, H, Guang, X, Al-Fageeh, MB, Cao, J, Pan, S, Zhou, H, et al. Camelid genomes reveal evolution and adaptation to desert environments. Nat Commun. (2014) 5:5188. doi: 10.1038/ncomms6188

PubMed Abstract | Crossref Full Text | Google Scholar

38. Goli, RC, Mahar, K, Manohar, PS, Chishi, KG, Prabhu, IG, Choudhary, S, et al. Insights from homozygous signatures of Cervus nippon revealed genetic architecture for components of fitness. Mamm Genome. (2024) 35:657–72. doi: 10.1007/s00335-024-10064-1

PubMed Abstract | Crossref Full Text | Google Scholar

39. Eydivandi, S, Roudbar, MA, Karimi, MO, and Sahana, G. Genomic scans for selective sweeps through haplotype homozygosity and allelic fixation in 14 indigenous sheep breeds from Middle East and South Asia. Sci Rep. (2021) 11:2834. doi: 10.1038/s41598-021-82625-2

PubMed Abstract | Crossref Full Text | Google Scholar

40. Yang, H, Yang, Y-L, Li, G-Q, Yu, Q, and Yang, J. Identifications of immune-responsive genes for adaptative traits by comparative transcriptome analysis of spleen tissue from Kazakh and Suffolk sheep. Sci Rep. (2021) 11:3157. doi: 10.1038/s41598-021-82878-x

PubMed Abstract | Crossref Full Text | Google Scholar

41. Weldenegodguad, M, Popov, R, Pokharel, K, Ammosov, I, Ming, Y, Ivanova, Z, et al. Whole-genome sequencing of three native cattle breeds originating from the northernmost cattle farming regions. Front Genet. (2018) 9:728. doi: 10.3389/fgene.2018.00728

PubMed Abstract | Crossref Full Text | Google Scholar

42. Lv, F-H, Agha, S, Kantanen, J, Colli, L, Stucki, S, Kijas, JW, et al. Adaptations to climate-mediated selective pressures in sheep. Mol Biol Evol. (2014) 31:3324–43. doi: 10.1093/molbev/msu264

PubMed Abstract | Crossref Full Text | Google Scholar

43. Lee, J, Mun, H, Koo, Y, Park, S, Kim, J, Yu, S, et al. Enhancing genomic prediction accuracy for body conformation traits in Korean Holstein cattle. Animals. (2024) 14:1052. doi: 10.3390/ani14071052

PubMed Abstract | Crossref Full Text | Google Scholar

44. Ahmad, SF, Gangwar, M, Kumar, A, Kumar, A, Dige, MS, Jha, GK, et al. Dissecting genomes of multiple yak populations: unveiling ancestry and high-altitude adaptation through whole-genome resequencing analysis. BMC Genomics. (2025) 26:214. doi: 10.1186/s12864-025-11387-2

PubMed Abstract | Crossref Full Text | Google Scholar

45. Dang, D, Zhang, L, Gao, L, Peng, L, Chen, J, and Yang, L. Analysis of genomic copy number variations through whole-genome scan in Yunling cattle. Front Vet Sci. (2024) 11:1413504. doi: 10.3389/fvets.2024.1413504

PubMed Abstract | Crossref Full Text | Google Scholar

46. Zhang, L, Liu, J, Zhao, F, Ren, H, Xu, L, Lu, J, et al. Genome-wide association studies for growth and meat production traits in sheep. PLoS One. (2013) 8:e66569. doi: 10.1371/journal.pone.0066569

PubMed Abstract | Crossref Full Text | Google Scholar

47. Han, Y, Liu, J, Zhang, C, Sun, M, Li, X, Liu, H, et al. Relationship between the CUBN and the MIA3 gene copy number variation and growth traits in different cattle breeds. Anim Biotechnol. (2025) 36:2450355. doi: 10.1080/10495398.2025.2450355

PubMed Abstract | Crossref Full Text | Google Scholar

48. Zhang, Y, Huang, H, Zhao, G, Yokoyama, T, Vega, H, Huang, Y, et al. ATP6V1H deficiency impairs bone development through activation of MMP9 and MMP13. PLoS Genet. (2017) 13:e1006481. doi: 10.1371/journal.pgen.1006481

PubMed Abstract | Crossref Full Text | Google Scholar

49. Liu, Y, Jiao, Y, Li, P, and Zan, L. MALDI-TOF-MS-based high throughput genotyping of mutations associated with body measurement traits in cattle. Mamm Genome. (2020) 31:228–39. doi: 10.1007/s00335-020-09840-6

PubMed Abstract | Crossref Full Text | Google Scholar

50. Mo, S, Kim, MK, Jang, JS, Lee, SH, Hong, SJ, Jung, S, et al. Unique expression and critical role of metallothionein 3 in the control of osteoclastogenesis and osteoporosis. Exp Mol Med. (2024) 56:1791–806. doi: 10.1038/s12276-024-01290-3

PubMed Abstract | Crossref Full Text | Google Scholar

51. Rowland, K, Saelao, P, Wang, Y, Fulton, JE, Liebe, GN, McCarron, AM, et al. Association of candidate genes with response to heat and Newcastle disease virus. Genes. (2018) 9:560. doi: 10.3390/genes9110560

PubMed Abstract | Crossref Full Text | Google Scholar

52. Lange, S, Ouyang, K, Meyer, G, Cui, L, Cheng, H, Lieber, RL, et al. Obscurin determines the architecture of the longitudinal sarcoplasmic reticulum. J Cell Sci. (2009) 122:2640–50. doi: 10.1242/jcs.046193

PubMed Abstract | Crossref Full Text | Google Scholar

53. Jaquet, V, Wallerich, S, Voegeli, S, Túrós, D, Viloria, EC, and Becskei, A. Determinants of the temperature adaptation of mRNA degradation. Nucleic Acids Res. (2022) 50:1092–110. doi: 10.1093/nar/gkab1261

PubMed Abstract | Crossref Full Text | Google Scholar

54. Xing, S, Gai, K, Li, X, Shao, P, Zeng, Z, Zhao, X, et al. Tcf1 and Lef1 are required for the immunosuppressive function of regulatory T cells. J Exp Med. (2019) 216:847–66. doi: 10.1084/jem.20182010

PubMed Abstract | Crossref Full Text | Google Scholar

55. Ibrahim, MA, Tolone, M, Barbato, M, Alsubaie, FM, Alrefaei, AF, and Almutairi, M. Geographical distribution, genetic diversity, and environmental adaptations of dromedary camel breeds in Saudi Arabia. Front Vet Sci. (2025) 11:1490186. doi: 10.3389/fvets.2024.1490186

PubMed Abstract | Crossref Full Text | Google Scholar

56. Bazile, J, Jaffrezic, F, Dehais, P, Reichstadt, M, Klopp, C, Laloe, D, et al. Molecular signatures of muscle growth and composition deciphered by the meta-analysis of age-related public transcriptomics data. Physiol Genomics. (2020) 52:322–32. doi: 10.1152/physiolgenomics.00020.2020

PubMed Abstract | Crossref Full Text | Google Scholar