Genome-Wide Analysis of Genetic Diversity and Selection Signatures in the Jiami Donkey

1. introduction

The donkey (Equus asinus) is an important domestic livestock species that has historically contributed to agriculture through transportation, draft work, and the production of meat, milk, and hides. Modern domestic donkeys originated from the African wild ass (Equus africanus), including the Nubian (E. a. africanus) and Somali (E. a. somaliensis) subspecies [1]. Based on genomic evidence, domestic donkeys can be broadly classified into three major genetic groups: Tropical African, North African–Eurasian, and Australian [2]. China has a donkey domestication history spanning more than 4,000 years and possesses abundant indigenous donkey genetic resources [3]. Long-term natural and artificial selection has resulted in the formation of numerous local breeds that exhibit distinct phenotypic and adaptive characteristics across different geographical regions [4].

The Jiami donkey is primarily distributed in Yulin and Yan’an, northern Shaanxi Province, on the Loess Plateau. This region is characterized by a dry climate, pronounced seasonal variation, frequent sandstorms, and large diurnal temperature fluctuations. Consequently, the Jiami donkey has evolved a compact body conformation, well-developed musculature, and a thick coat, enabling it to withstand harsh climatic conditions. Historically, it served as an important draft and riding animal in local agricultural production. However, the widespread mechanization of agriculture has led to a continuous decline in the breeding population, particularly of purebred Jiami donkeys, raising concerns regarding the conservation of their genetic diversity and long-term sustainability. Therefore, understanding the genomic characteristics of this breed is essential for developing effective conservation and breeding strategies.

Whole-genome sequencing (WGS) has become a powerful tool for investigating genetic diversity, population history, and adaptive evolution in domestic animals [5,6]. It has been widely applied to identify genomic regions associated with economically important traits and environmental adaptation in livestock. For instance, Abigar cattle have been studied for their adaptation to local disease pressures [7], Tibetan cattle have been extensively studied for their high-altitude adaptations [8], and Sanjiang cattle for their heat tolerance [9]. In pigs, the HMGA1 gene, associated with body size, has been identified in Yorkshire pigs [10], and the KIT gene is associated with coat color [11]. Numerous genes related to economic traits and environmental adaptability have been found in goats. For example, the DSG3 and FGF5 genes are linked to high-altitude adaptation [12,13], while the ASIP and KITLG genes are associated with coat color [14]. In donkeys, previous WGS studies have elucidated breed origin and evolution, chromosome evolution, body size, coat color, locomotion ability, and environmental adaptation [3,15]. For example, NCAPG-LCORL has been associated with body size in Liangzhou donkeys, whereas CYP4A11 contributes to drought adaptation [16]. Despite these advances, the genomic diversity, population structure, and selection signatures of the Jiami donkey remain largely unexplored.

In this study, we performed whole-genome resequencing of 30 Jiami donkeys and integrated these data with publicly available genomes from other domestic donkey breeds. Using population genetic analyses and multiple complementary selection-scan approaches, we characterized the genetic diversity, population structure, demographic history, and genomic signatures of selection in the Jiami donkey. Our findings provide valuable genomic resources and candidate loci for the conservation and genetic improvement of this indigenous breed.

2. Materials and Methods

2.1. Sample Collection, Dna Extraction, and Next-Generation Sequencing

Blood samples were collected from 30 Jiami donkeys (Equus asinus) in Yulin City, Shaanxi Province, China. Genomic DNA was extracted using the standard phenol–chloroform method. Sequencing libraries with an insert size of approximately 350 bp were prepared for each individual according to the manufacturer’s protocol. Whole-genome sequencing was performed on the Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, USA) using a 150 bp paired-end strategy at Novogene Bioinformatics Institute (Beijing, China).

In addition, 210 previously published whole-genome datasets of domestic donkeys were incorporated into this study, including 33 foreign and 177 Chinese individuals. The Chinese samples represented multiple indigenous breeds, including Sichuan, Hetianqing, Turpan, Guangling, Liangzhou, Yangyuan, Guanzhong, Dezhou, Qinghai, Tibetan, Yunnan, Qingyang, Biyang, and Xinjiang donkeys. Together with the newly generated data, a total of 240 domestic donkey genomes were included for downstream analyses (Table S1). All newly generated and publicly available raw sequencing reads were processed using the same bioinformatics pipeline. Specifically, all FASTQ files were independently aligned to the EquAss-T2T_v2 reference genome, followed by uniform quality control, variant calling, and joint genotyping using the BWA–GATK workflow described below. Additionally, 4 wild-ass genomes were included in specific analyses and were used as the outgroup where appropriate.

2.2. Short-Reads Mapping and Snp Calling

Clean reads were aligned to the donkey reference genome (EquAss-T2T_v2) using BWA with default parameters [17]. The resulting BAM files were processed using SAMtools, Picard Tools, and GATK v4.0 [18]. SNPs were identified following the GATK Best Practices workflow, including the HaplotypeCaller, GenotypeGVCFs, SelectVariants, and VariantFiltration modules.

To obtain high-confidence variants, hard filtering was performed using the following criteria: QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < −12.5, ReadPosRankSum < −8.0, and SOR > 3.0. Variants with a mean sequencing depth lower than one-third or greater than three times the genome-wide average were excluded. Functional annotation of SNPs was performed using ANNOVAR [19].

2.3. Population Genetic Analyses

2.3.1. Population Structure

SNPs were pruned for linkage disequilibrium using PLINK v1.90 [20] with the parameter --indep-pairwise 50 5 0.2. An unrooted neighbor-joining (NJ) tree was constructed from pairwise genetic distances using MEGA v11.0 [21] and visualized with iTOL v7 (Interactive Tree of Life v7; https://itol.embl.de/). Principal component analysis (PCA) was performed using the smartPCA module implemented in EIGENSOFT v5.0 [22]. Population genetic structure was inferred using ADMIXTURE v1.3 [23], with the number of ancestral populations (K) ranging from 2 to 5.

2.3.2. Genetic Diversity

Linkage disequilibrium (LD) decay was estimated using PopLDdecay v3.42 [24]. Nucleotide diversity (θπ) and the inbreeding coefficient (F) were calculated using VCFtools v0.1.16 [25], whereas observed heterozygosity (Ho) and expected heterozygosity (He) were calculated using PLINK.

Runs of homozygosity (ROH) were identified using the --homozyg function in PLINK. The parameters were as follows: (1) a sliding window size of 50 kb (-homozyg-window-snp 50); (2) a required minimum density of 1000 (-homozyg-density 1000); (3) allowance for one heterozygote per window (-homozyg-window-het 1); (4) allowance for up to one missing pair per window (-homozyg-window-missing 1). Only ROHs longer than 500 kb were retained and subsequently classified into four length categories (0.5–1 Mb, 1–2 Mb, 2–4 Mb, and >4 Mb). Figures were generated using R.

Because LD decay, nucleotide diversity, and ROH analyses are sensitive to sample size, random downsampling was performed to minimize potential bias caused by unequal numbers of individuals among breeds.

2.3.3. Demographic History

Pairwise genetic differentiation (FST) between Jiami donkeys and other donkey populations, was calculated using the smartPCA module. Historical gene flow among Chinese indigenous donkey breeds was inferred using TreeMix v1.12 [26], with the number of migration events (m) varying from 1 to 10. The optimal migration model was determined according to the second-order rate of change in likelihood (Δm) [27].

Historical effective population size (Ne) trajectories were reconstructed separately for Jiami donkeys and the representative domestic donkey populations shown in Figure 3A using SMC++ v1.15.5 [28] based on unphased autosomal whole-genome data. A mutation rate of 7.242×10^-9 mutations/site/generation and a generation interval of 8 years were assumed. The same parameter settings were applied to all populations to facilitate comparison of demographic trajectories.

2.4. Genome-Wide Selection Signal Analysis

To identify genomic regions potentially subjected to recent positive selection within the Jiami donkey population, three complementary statistics were applied: nucleotide diversity (θπ), the composite likelihood ratio (CLR), and the integrated haplotype score (iHS). Nucleotide diversity was calculated using VCFtools v0.1.16 in 50-kb sliding windows with a 20-kb step. CLR scores were calculated using SweepFinder2 [29] and summarized across the same genomic windows.The iHS statistic was subsequently calculated using Selscan v1.1 [30] and normalized using the accompanying normalization procedure. Absolute normalized iHS values (|iHS|) were used to identify unusually extended haplotypes irrespective of allele orientation. Windows in the lowest 1% of the empirical θπ distribution or the highest 1% of the CLR and |iHS| distributions were considered candidate regions under selection. Genes located within candidate regions identified by all three methods were considered high-confidence candidates for selection within the Jiami donkey population.

To identify genomic regions showing divergent selection between Jiami and Guanzhong donkeys, three complementary statistics were used: the fixation index (FST), nucleotide-diversity ratio (θπ-ratio), and cross-population extended haplotype homozygosity (XP-EHH). Genome-wide FST was calculated using VCFtools in 50-kb sliding windows with a 20-kb step. The θπ-ratio was defined as θπ_JM/θπ_GZ, such that reduced nucleotide diversity in Jiami donkeys resulted in lower θπ-ratio values. XP-EHH was calculated using Selscan v1.1 based on phased genotypes and subsequently normalized. Positive XP-EHH values indicated longer-range haplotype homozygosity in Jiami donkeys relative to Guanzhong donkeys and were therefore interpreted as evidence consistent with selection in the Jiami population. Candidate windows were defined as those falling within the highest 1% of the empirical FST distribution, the lowest 1% of the log2(θπ_GZ/θπ_JM) distribution, or the highest 1% of the XP-EHH distribution. Genomic regions supported by at least two of the three methods were retained as candidate regions to reduce method-specific false-positive signals.

Tajima’s D was further calculated using VCFtools for selected candidate regions to characterize local allele-frequency spectra. Functional enrichment analyses of candidate genes were performed using KOBAS v3.0, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. Enriched terms with P < 0.05 were considered significant.

2.5. Statistical Analysis

Unless otherwise specified, all statistical analyses were performed using R (v4.3.0). Figures were generated using the ggplot2 package, and default parameters were applied for all software unless otherwise stated.

3. results

3.1. Genome Resequencing and Snp Identification

High-throughput sequencing generated whole-genome data for 30 Jiami donkeys, with an average sequencing depth of 10.71× and an average mapping rate of 99.88% (Table S1). In addition, 210 previously published whole-genome datasets representing other domestic donkey breeds were incorporated for comparative analyses. After variant calling and quality filtering, 23,752,978 high-quality biallelic SNPs were retained for downstream analyses. Functional annotation using ANNOVAR (Table S2) showed that 67.4% of SNPs were located within introns, 16.5% in intergenic regions, and 2.4% in exons. Among the exonic variants, 336,653 were nonsynonymous and 261,262 were synonymous (Figure 1B).

3.2. Population Structure and Phylogenetic Relationships

Population structure was investigated using ADMIXTURE, PCA, and neighbor-joining (NJ) phylogenetic analysis. ADMIXTURE analysis of the 240 individuals (K = 2-4) revealed clear genetic stratification among domestic donkey populations (Figure 1E). At K = 2, the lowest cross-validation(CV) error (0.3167) was obtained, and Chinese domestic donkeys were clearly differentiated from foreign domestic populations. At K = 3, a distinct ancestral component emerged in both Jiami and Guanzhong donkeys. At K = 4, Jiami donkeys exhibited a unique ancestral component, whereas most other Chinese indigenous breeds remained genetically homogeneous.

PCA produced consistent results, with PC1 (1.66%) separating foreign from Chinese donkey breeds and PC2 (1.31%) distinguishing Jiami donkeys from other domestic populations (Figure 1D). Similarly, the NJ tree showed clear geographic clustering and supported the genetic distinctiveness of the Jiami donkey (Figure 1C). Together, these analyses consistently demonstrated that Jiami donkeys represent a genetically distinct domestic donkey population.

3.3. Genetic Diversity

Genetic diversity was evaluated using linkage disequilibrium (LD), nucleotide diversity (θπ), runs of homozygosity (ROH), heterozygosity, and inbreeding coefficients. Genome-wide LD analysis showed that Jiami donkeys exhibited a relatively low LD level among Chinese breeds, comparable to that of Dezhou and Hetianqing donkeys (Figure 2A). Similarly, Jiami donkeys displayed an above-average level of nucleotide diversity compared with other Chinese indigenous breeds (Figure 2B, Table S3).

ROH analysis revealed that most homozygous segments across all breeds were between 0.5 and 1 Mb in length. Jiami donkeys exhibited the highest number of ROHs across all four size categories (Figure 2C). The inbreeding coefficient (F) of the Jiami donkey was intermediate among the investigated breeds (Figure 2D). Consistent with these observations, the observed heterozygosity (Ho = 0.1462) was slightly lower than the expected heterozygosity (He = 0.1598) (Table S4). Overall, these results indicate that Jiami donkeys maintain moderate genetic diversity despite evidence of localized homozygosity.

3.4. Genetic Differentiation, Migration, and Demographic History

Pairwise genetic differentiation between Jiami donkeys and other populations was evaluated using the fixation index (FST). Among domestic donkey breeds, Jiami donkeys showed the lowest genetic differentiation from Dezhou donkeys (FST = 0.032), followed by Qingyang donkeys (FST = 0.034). In contrast, the greatest differentiation among domestic populations was observed between Jiami and foreign donkey breeds (FST = 0.069) (Figure 3E). As expected, substantially higher differentiation was detected between Jiami donkeys and wild asses (FST = 0.269) (Table S5).

Gene flow was inferred using TreeMix, with the wild ass designated as the outgroup. Based on the likelihood curve and the Evanno-like Δm criterion, the optimal number of migration events was determined to be m = 3. Under the optimal m = 3 model, no major migration edge involving Jiami donkeys was detected, suggesting that strong directional gene flow involving the Jiami population was not evident under the fitted TreeMix model. (Figure 3B–D).

Historical effective population size (Ne) was reconstructed using SMC++. The Jiami donkey exhibited a long-term decline in effective population size over the past 100,000 years, despite several temporary fluctuations. Compared with representative domestic breeds, Jiami donkeys showed a relatively small inferred effective population size toward the recent end of the reconstructed trajectory (Figure 3A).

3.5. Genome-Wide Signatures of Selection Within the Jiami Donkey

Genome-wide signatures of recent positive selection were identified using θπ, CLR, and iHS analyses. Candidate selective regions were defined as genomic windows with empirical P < 0.01 supported by all three methods (Figure 4B). A total of 896, 113, and 1301 candidate genes were identified by the iHS, CLR, and θπ analyses, respectively, with 30 genes detected by all three approaches (Figure 4A, Table S6).

Functional annotation indicated that these overlapping genes were primarily associated with growth and body size (LCORL, NCAPG, ASPA, SCARF1, and FAM184B), reproduction (SMG6, SPATA22, and SLC43A2), immune function (CLDN25, FBXO7, and RPH3AL), nervous system development (TMPRSS5 and RAP1GAP2), and disease-related processes (PWP1 and SGSM2).

Functional enrichment analysis identified three significantly enriched KEGG pathways and 19 GO terms (P < 0.05). These pathways were mainly associated with amino acid metabolism, whereas the enriched GO terms included histone methyltransferase activity, telomerase holoenzyme complex, and telomerase RNA binding (Figure 4C).

3.6. Comparative Selection Analysis Between Jiami and Guanzhong Donkeys

To identify genomic regions potentially associated with local adaptation, Jiami and Guanzhong donkeys were compared using FST, θπ-ratio, and XP-EHH analyses. Genomic windows with empirical P < 0.01 supported by at least two independent methods were considered candidate selective regions. A total of 1491, 1039, and 536 genes were detected by FST, θπ-ratio, and XP-EHH, respectively, with 427 candidate genes shared by at least two methods (Figure 5A, Table S7). Functional annotation showed that these candidate genes were associated with nervous system development (ASTN2, KMT2C, UBA6, BST1, and DRD2), lipid metabolism and reproduction (LRP1B), DNA damage repair (HUS1), and heat tolerance (SIRT1) (Figure 5B).

Functional enrichment analysis revealed significant enrichment in three KEGG pathways and 16 GO terms (P < 0.05), including lysine degradation, RNA degradation, and cell adhesion molecules (CAMs), together with GO categories related to neuron cell-cell adhesion, neuron migration, histone methyltransferase activity, and neurexin family protein binding (Figure 5C). Among these candidates, BTG4, ASTN2, and NLGN1 were highlighted because of their consistent population-genetic signals and their reported biological functions in other mammals.

4. Discussion

This study provides a dedicated population-genomic characterization of the Jiami donkey based on newly generated whole-genome resequencing data and publicly available genomes from representative donkey populations. Previous genomic studies have included Jiami donkeys as part of broader surveys of Chinese indigenous breeds and have revealed substantial geographic structuring and selection among Chinese donkey populations. By expanding the sampling of Jiami donkeys and focusing specifically on their genetic diversity, demographic history, and genomic signatures of selection, the present study provides a more detailed assessment of the genomic characteristics and conservation status of this local breed.

Population structure analyses consistently indicated that Jiami donkeys share a broad genetic background with other Chinese indigenous donkey populations while retaining a distinguishable genomic composition. ADMIXTURE, PCA, and phylogenetic analyses all revealed clear differentiation between Chinese and foreign donkey populations, consistent with previous genomic studies showing that most Chinese domestic donkeys share a predominant ancestral component but exhibit geographic genetic differentiation among regional breeds [4,31]. The relatively close relationship between Jiami and Guanzhong donkeys is also consistent with their geographic proximity and likely reflects shared ancestry among neighboring indigenous populations. In contrast, the greater differentiation between Jiami donkeys and foreign domestic donkeys or wild asses reflects their distinct demographic and evolutionary histories.

Genetic diversity is an important component of the long-term evolutionary potential and conservation value of livestock populations [32]. Jiami donkeys exhibited relatively low linkage disequilibrium and comparatively high nucleotide diversity among the Chinese indigenous populations investigated in this study, indicating that substantial genetic variation remains within the breed. At the same time, Jiami donkeys possessed a relatively large number of ROHs, most of which were shorter than 1 Mb. Short ROHs are generally more consistent with older demographic events and historical inbreeding, whereas long ROHs more commonly reflect recent consanguinity [33]. The predominance of short ROHs, together with intermediate inbreeding coefficients and heterozygosity levels, therefore suggests that the observed homozygosity in Jiami donkeys is more strongly associated with historical demographic processes than with intensive recent inbreeding. These observations indicate that the breed retains appreciable genomic diversity, although continued monitoring of inbreeding will remain important for future conservation management.

Historical demographic reconstruction further revealed an overall decline in the effective population size of Jiami donkeys over the inferred long-term demographic trajectory. This pattern primarily reflects historical demographic processes and should not be directly attributed to recent changes in agricultural practices. Nevertheless, the documented decline in the census population of Jiami donkeys following the reduced use of draft animals represents an additional contemporary conservation concern. A persistently small breeding population may increase genetic drift and accelerate the loss of rare alleles, potentially reducing the long-term adaptive capacity of the breed. Maintaining an adequate effective breeding population and avoiding excessive use of a limited number of sires should therefore be considered in future conservation and breeding programs.

Genome-wide selection scans identified candidate regions associated with several biologically important processes, including growth, reproduction, immune function, and environmental adaptation. The involvement of multiple functional categories suggests that the selective history of Jiami donkeys has affected several biological processes rather than a single major pathway. Among the candidate loci identified within the Jiami population, the NCAPG–LCORL region is particularly noteworthy. This genomic region has repeatedly been associated with body size and growth in livestock, and previous studies of Chinese donkeys have also identified NCAPG and LCORL as candidate genes contributing to body-size differentiation [4,16]. Experimental and association studies in cattle and sheep further support roles of NCAPG in skeletal growth, muscle development, and related production traits [34,35]. Its identification in the present analysis therefore suggests that this region may have contributed to the characteristic body conformation of the Jiami donkey, although direct phenotypic association in this breed remains to be established.

To investigate genomic differentiation potentially associated with local adaptation, we further compared Jiami and Guanzhong donkeys. These populations are geographically close but inhabit contrasting ecological environments: northern Shaanxi is relatively cold and dry and is characterized by large seasonal and diurnal climatic variation, whereas the Guanzhong Plain generally has a warmer and more humid environment. Such ecological differences may have imposed different selective pressures on the two populations. By integrating FST, θπ-ratio, and XP-EHH, we identified 427 candidate genes supported by at least two approaches. Functional enrichment highlighted pathways related to lysine degradation, RNA degradation, and cell adhesion, together with GO terms involving neuronal development, cell adhesion, and epigenetic regulation. These results suggest that differentiation between Jiami and Guanzhong donkeys involves multiple cellular and regulatory processes and provides candidate genomic regions for further investigation of local adaptation.

Among the differentiated candidate genes, BTG4, ASTN2, and NLGN1 were of particular interest because they were supported by multiple population-genetic signals and have biologically relevant functions reported in other mammals. BTG4 is a maternal-effect gene that plays an essential role in oocyte maturation and the maternal-to-zygotic transition; loss of Btg4 in mice results in female infertility and early embryonic developmental arrest [36]. Its identification therefore makes BTG4 a plausible candidate locus associated with reproductive processes, although its specific effects in donkeys remain unknown. ASTN2 is involved in neuronal development and synaptic regulation, and experimental studies have demonstrated altered neural and behavioral phenotypes following disruption of ASTN2 in mammalian models [37,38]. NLGN1 encodes a synaptic adhesion protein with an established role in neuronal signaling and synaptic function [39]. The occurrence of ASTN2 and NLGN1 among the differentiated loci is also consistent with the enrichment of neuronal-development and cell-adhesion categories observed in our analysis. Nevertheless, selection scans identify genomic regions associated with population differentiation rather than causal genes or phenotypes, and therefore these loci should be regarded as candidates requiring further functional validation.

The present study has several limitations. First, the biological functions of the candidate genes were inferred from population-genomic signals and knowledge obtained primarily from other species; direct functional evidence in donkeys is currently limited. Second, although the integration of several complementary selection statistics reduces dependence on any single method, demographic history and population structure can still influence genome-wide selection signals. Future studies incorporating larger and geographically broader Jiami donkey populations, together with transcriptomic, epigenomic, and phenotype-associated data, will be important for validating the candidate loci identified here and clarifying their functional contributions. Overall, our results provide a population-genomic framework for understanding the genetic characteristics of the Jiami donkey and offer useful genomic information for its conservation, genetic improvement, and sustainable utilization.

5. Conclusion

Overall, this study provides a comprehensive overview of the genomic diversity, demographic history, and selection signatures of the Jiami donkey. Our findings not only expand current knowledge of the evolutionary history of Chinese indigenous donkeys but also provide valuable genomic resources for the conservation, genetic improvement, and sustainable utilization of the Jiami donkey.