A total of 120 donkeys were selected for whole-genome sequencing in Yangyuan County, Hebei Province, China. All donkeys were adult jennies with unknow pedigree information. The animal feed consisted of grass and hay ad libitum and water. Blood samples were taken from the jugular vein and rapidly stored at −80°C. Genomic DNA was extracted from the blood using a TIANamp genomic DNA kit (TIANGEN, DP304, China). The quality and integrity of the extracted DNA were examined based on the A260/A280 ratio and 1% agarose gel electrophoresis. We also conducted measurements of the body height, body length, chest circumference, and shin girth. EXCEL 2019 was used to generate the descriptive statistics of the four body size traits.

Genomic libraries with an insert size of ∼400 bp were constructed and sequenced using an Illumina NovaSeq instrument (2 × 150 paired-end mode), which was conducted by the Personalbio Biotechnology Co., Ltd. (Shanghai, China). FastQC was used to provide quality control checks on the raw sequence data (http://www.bioinformatics.babraham. ac.uk/projects/fastqc). Fastp (Chen et al., 2018) was used to perform data filtering to provide clean data for downstream analysis with the following criteria: 1) removal of the joint contamination at the 3′end; 2) sliding window quality pruning, where if the average Q value ≤ 20, then the bases in the window would be marked as discarded and the window would be moved forward by one base; and 3) if the length of any of the double-end reads was ≤50 bp, the double-end reads would be removed.

High-quality trimmed read pairs were aligned to the reference donkey genome (GenBank: PRJNA431818) by Burrows-Wheeler Aligner (BWA, versin0.7.12-r1039) with default parameters and SAM files were generated (Li and Durbin, 2009). The SAM files were sequenced using Picard1.107 software. Duplicates were removed using the MarkDuplicates command in the Picard package. Thus, if multiple paired reads had the same chromosome coordinates after comparison, only the paired reads with the highest score would be retained, and the reads near InDel were most prone to mapping errors. To minimize SNPs caused by mapping errors, it was necessary to rematch the reads near InDel to improve the accuracy of SNP calling, using the IndelRealigner command in the Genome Analysis Toolkit (GATK) program to rematch all of the reads near InDel to improve the accuracy of SNP prediction. SNP prediction accuracy was improved by using the IndelRealigner command in the GATK program to re-compare all of the reads near InDel. For all 120 samples, mapping statistics based on high-quality mapped reads for each accession included the average coverage depth and the proportion of the donkey genome covered by different read depths.

High-quality mapped reads were used for variant calling with GATK (version 3.8) (Zhu et al., 2015). To improve the accuracy of SNP calling, two steps were conducted, including using GATK to output files that contained all possible InDels and using IndelRealigner to re-compare all InDels near reads. The SNPs that met the following criteria were filtered out to ensure reliability: 1) Fisher test of strand bias (FS) ≤ 60; 2) HaplotypeScore ≤13.0; 3) Mapping Quality (MQ) ≥ 40; 4) Quality Depth (QD) ≥ 2; 5) ReadPosRankSum ≥ −8.0; 6) MQRankSum > −12.5; 7) Read depth (DP) > 4; and 8) --maf 0.01; 9) --max-missing 0.7. All SNPs were annotated using ANNOVAR (Wang et al., 2010) software with default parameter settings.

PCA of the 120 samples was conducted with GCTA (Yang et al., 2011) software. To evaluate the phylogenetic relationships between individuals, a maximum likelihood (ML) tree was built based on the SNPs using Fasttree software (Price et al., 2009), with n = 1,000 bootstrap replications. We used PopLDdecay (Zhang C. et al., 2019) software with default parameters to calculate linkage disequilibrium (LD) decay.

GWASs for body size traits, including body height, body length, chest circumference, and shin girth, were performed using linear mixed models (LMMs) by the efficient mixed-model association expedited (EMMAX) software (http://genetics.cs.ucla.edu/emmax.) (Kang et al., 2010). The model is y = X α + Z β + W μ + ε where y is a vector of phenotypic records (body height, body length, chest circumference or shin girth); X represents the vector of the covariates (the PCA principal components obtained from the analysis of population structure), including a column of 1; α represents the vector of the corresponding coefficients, including the intercept; Z represents the vector of the marker genotypes, β represents the effect size of the marker; W μ represents the vector of the random polygenic effects, with μ ∼ (0, G σ a 2 ), where σ a 2 is the genetic variance, and G represents the genomic relationship matrix constructed with identity-by-state (IBS); ε is the vector of random residuals with ε ∼N (0, I σ 2 ), where σ 2 is the residual variance.

The Manhattan and quantile-quantile plot (Q-Q) plots were drawn by R packages. The SNPs reached the genome-wide significant level -log10 (p-value) = 7 were determined as the significant loci. The search for candidate genes was extended from 50 kb upstream and downstream from the significant SNPs.

To provide insight into the functional enrichment of the candidate genes, we performed GO enrichment analysis (Ashburner et al., 2000) and the KEGG pathway (Kanehisa et al., 2004). GO enrichment analysis was performed using InterProScan. A hypergeometric p-value was calculated where the background was set as genes in the entire genome. GO terms with p-value <0.05 were considered significantly enriched, and GO enrichment analysis elucidated the molecular function, biological process, and cellular components. KEGG pathway analysis was performed using KEGG Automatic Annotation Server (KAAS) to predict the main biological pathways involved in the candidate genes.

The body height, body length, chest circumference, and shin girth were collected from 120 Yangyuan jennies. The raw measurement data are detailed in Supplementary Table S1. Ranking from high to low according to body height, we divided the population into two groups, the top 60 donkeys are the larger ones, and the bottom 60 donkeys are the smaller ones. The order of the samples in Supplementary Table S1 is the result of ranking according to body height. The general descriptive statistics of the body size traits are shown in Table 1. The average body height was 139.33 cm with a standard deviation of 10.94, the average body length was 138.62 cm with a standard deviation of 10.43, the average chest circumference was 138.62 cm with a standard deviation of 9.68, and the average shin girth was 15.93 cm with a standard deviation of 1.19.

To sequence the Yangyuan donkey genome, libraries with 400 bp inserts were sequenced and 861 Gb of raw sequence data were produced. The quality of the raw data is shown in Supplementary Table S2, with Q20 ≥ 96.7, Q30 ≥ 91.48, and GC content from 43.08 to 44.00, which confirmed the high quality of the sequence data. The construction of the library and sequencing were successful. To enhance the quality of the data, Fastp was used to generate high quality data. In total, 718 Gb of clean data were generated. The filtered clean data were aligned to the reference genome, and the mapping rate was 99.57%–99.79% (Supplementary Table S3), implying that the mapping rate could be used for subsequent variant detection. In general, by using next-generation sequencing technology, and after filtering and mapping, a total of 9,735,956 SNPs were identified using the Illumina NovaSeq sequencing platform with an average sequencing depth of ∼4.1 (Supplementary Table S4). Their distribution on the chromosomes is shown in Figure 1, and the SNP annotation results are shown in Table 2. We found that 65.8% of the SNPs were concentrated in the intergenic region, 31.48% of the SNPs were concentrated in the intronic region, 0.91% of the SNPs were located upstream, 0.8% of the SNPs were concentrated downstream, and only 0.98% of the SNPs were located in the exon region.

We performed principal component analysis (PCA) with 9,735,956 SNPs to analyze the population structure of 120 Yangyuan donkeys. PCA revealed the relationship between the first two principal components, and we found that most of the individuals were clustered together and only a few individuals showed outlier distribution, which demonstrated that the genetic background of the population was relatively uniform (Figure 2A). The phylogenetic tree showed that the samples were mixed together without obvious clustering into separate groups, which was consistent with PCA (Figure 2B).

Based on the detected SNPs, correlation analysis was performed according to the SNP and body height, body length, chest circumference, and shin girth traits to identify the markers and candidate genes that were closely associated with the target traits. In total, 16 SNPs reached a significance threshold of 7, and most of the SNPs were within the intronic, intergenic, and downstream locations of chromosomes 8, 13, and 31. The Manhattan plots for the four body size traits and corresponding Q-Q plots of the p-values against the expected p-values are presented in Figure 3. According to the LD analysis results (Supplementary Figure S1), all positional candidate genes were identified for 50 kb before and after the most significant SNP. Detailed information regarding the significant SNPs identified using GWAS and the putative candidate genes are presented in Table 3, including sphingomyelin phosphodiesterase 4 (SMPD4), ribosomal protein S6 kinase A6 (RPS6KA6), lysophosphatidic acid receptor 4 (LPAR4), glucagon-like peptide 2 receptor (GLP2R), bromodomain, and the WD repeat domain containing 3 (BRWD3), magnesium transporter 1 (MAGT1), zinc finger DHHC-type palmitoyltransferase 15 (ZDHHC15), and cysteinyl leukotriene receptor 1 (CYSLTR1). Other genes located in the 50-kb genome region flanking the significantly associated SNPs are shown in Supplementary Table S5.

The GO function annotation and KEGG pathway enrichment analyses were utilized to explore the biological functions of the candidate genes. The GO enrichment analysis results were classified based on the molecular function, biological process, and cellular components. The top 10 GO terms with the smallest p-value, indicating the most significant enrichment, were selected for display, and the results are shown in Figure 4. According to the KEGG enrichment results, the KEGG pathways with the smallest FDR values, which were the most enriched, are illustrated in Figure 5.