The experimental methods were approved by the Hankyong National University Animal Ethics Committee, Anseong, Republic of Korea (No.2021-1). A total of 269 whole blood samples were obtained from the four populations of the five geographically distinct provinces (Table 1, Figure 1). Sampling sites included i) East Mongolia Khentii Province (KTP, n = 53); ii) West Mongolia Uvs Province (USP, n = 43); iii) Southern Mongolia Omnogovi and Dundgovi provinces (GOP, n = 133); and iv) Northern Mongolia Khovsgol province (KGP, n = 40). Genomic DNA was extracted from blood samples using the methods described for QuickGene 810 (Kurabo, Osaka, Japan). The concentration and purity of the extracted genomic DNA were evaluated using an ND-1000 UV-Vis spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA).

The genetic diversity of Mongolian horse populations was identified using 14 MS markers (AHT4, ASB2, ASB17, ASB23, CA425, HMS1, HMS2, HMS3, HMS6, HMS7, HTG4, HTG6, HTG7, and VHL20), as recommended by the ISAG.

Multiplex polymerase chain reaction (PCR) was conducted using the Equine Genotypes Panel 1.1 Kit (Thermo Fisher Scientific, Waltham, MA, USA) for genotyping of the 14 MS markers. Following the manufacturer’s instructions, the reaction mixture was prepared using 2 μL genomic DNA (1.0 ng/μL), 9 μL master mix, and 9 μL primer mix, obtaining a total volume of 20 μL. The PCR was conducted using a GeneAmp PCR system 9700 (Applied Biosystems, Waltham, MA, USA). The PCR amplification was conducted using the following conditions: pre-denaturation at 98°C for 3 min followed by 30 cycles of 15 s at 98°C, 75 s at 60°C, and 1 cycle of 30 s at 72°C. The final extension step was performed at 72°C for 5 min and the mixture was cooled to 4°C.

Using Hi-Di formamide, amplified PCR products were diluted 50 to 100 times, and the diluted PCR products (1 μL) were further diluted using a mixture of 9 μL Hi-Di formamide and GeneScan 500LIZ size standard (Applied Biosystems, USA). After denaturation at 95°C for 3 min, capillary electrophoresis was conducted using an ABI 3730xl Genetic Analyzer (Applied Biosystems, USA). The size of each MS marker was determined using GeneMapper ver 5.0 (Applied Biosystems, USA). The determined alleles were collated using Microsoft Excel (Microsoft, Redmond, WA, USA) and used for statistical analysis.

The MS Toolkit software [11] was used to calculate the number of alleles, expected and observed heterozygosity (H_Exp and H_Obs, respectively), and polymorphism information content (PIC) values. To identify genetic correlations between populations, allele frequencies for each marker were confirmed using GenAlEx 6.4 [12] and Genetix software [13]. For this, principal coordinate analysis (PCoA) was conducted using GenAlEx 6.4 to display the two-dimensional coordinates. Factorial component analysis (FCA) was conducted using Genetix to display the three-dimensional coordinates. The genetic distance between populations was calculated based on allele frequencies according to Nei’s DA genetic distance [14] using POPTREE2 [15]. A phylogenetic tree was estimated from the genetic distances using the neighbor-joining method [16] with the DISPAN program [17]. Population structure [18,19] was used to estimate the uniformity of the population, and the K value was set to estimate the number of distinct populations (ΔK). To calculate the average estimate and standard deviation of each K value, the length of the burn-in period and the number of Markov chain Monte Carlo (MCMC) repetitions after the burn-in frequency were set, and the optimal K value and genetic uniformity for each cluster were calculated The results were applied to the Structure harvester [20] using the Evanno method [21].

The number of alleles, H_Exp, H_Obs, and PIC values for the Mongolian horse populations are summarized (Table 2). A total of 160 different alleles were detected, ranging from 6 (HTG7) to 22 (ASB2). H_Exp and H_Obs ranged from 0.572 (HTG4) to 0.885 (ASB17) and 0.561 (HTG4) to 0.868 (ASB17), with mean values of 0.772 and 0.757, respectively. PIC values ranged from 0.535 (HTG4) to 0.866 (ASB17), with a mean value of 0.737. Except for the HTG4 marker (H_Exp = 0.572, PIC = 0.535) used in this study, all other MS markers were considered to be highly useful in analyzing polymorphism of Mongolian horse. MS markers were previously used to assess the genotypic diversity of heterozygosity and PIC value during animal breed selection. According to Botstein et al [22], polymorphisms of MS markers are determined using the following criteria: if the sum of H_Exp is ≥0.6 and PIC is ≥0.5, then the marker is determined to be highly polymorphic. The heterozygosity (H_Exp and H_Obs) and PIC values with respect to the four Mongolian horse populations are summarized (Table 3). Among the provinces, the values for mean number of alleles (MNA), H_Exp, H_Obs, and PIC were highest for KTP horses (8.57, 0.787, 0.789, and 0.751, respectively) and lowest for KGP populations (7.07, 0.751, 0.738, and 0.706, respectively).

The genetic divergences among the four Mongolian horse populations based on allele frequencies were calculated according to Nei’s DA genetic distance [16]. The phylogenetic relationships among these Mongolian horse populations were determined using the neighbor-joining tree (Figure 2). Grouping values were determined by 1,000 repetitive “bootstrap” tests to determine the reliability of the neighbor-joining tree and were specified at the branching points of the tree. The smaller the estimated genetic distance is, the closer each population is genetically. In Nei’s DA genetic distance, the KTP and USP populations were the closest (DA = 0.0535), and the largest difference was calculated for the GOP and KGP populations (DA = 0.1797) (Supplementary Table S1). In the phylogenetic tree, the Mongolian horse populations was mainly divided into two clusters. The KTP and GOP populations formed the first cluster, and the KGP population formed the second cluster. The USP breed was located between the two clusters. The division of the populations into two distinct clusters highlighted the presence of clear genetic separation between each province.

As the phylogenetic tree may not take into account the effects of admixture among the four populations, we conducted PCoA and FCA, using allele frequencies of the 14 MS markers, as an alternative approach to understand the genetic relationships among populations. The PCoA contributed to 100% of the variation, including the third ingredient. The first three principal coordinates represented 60.13%, 29.75%, and 10.12% of the total variation, respectively. The KGP and GOP populations were distinct from the other populations (Figure 3). In contrast, the KTP and USP populations were located between the closest genetic distances. The FCA analysis revealed that the three dimensions contributed to a total of 100 with Axis 1 (56.08), Axis 2 (29.86), and Axis 3 (14.06). The KTP and GOP populations were genetically differentiated from other populations and showed a close relationship among certain populations, including KTP and USP (Figure 4). Thus, results of the PCoA analysis were consistent with those of the FCA.

A Bayesian clustering method and population structure were used for clustering algorithms of multif-locus genotypes to identify the population structure and pattern of admixture within the populations. Population structure was used to estimate the number of groups that were classified by the group surveyed. The Bayes analysis set K values from 2 to 4 and examined the formation of colonies by group (Figure 5). The bar plot in Figure 5 confirms that the KGP and GOP populations showed different clusters from the other two Mongolian horse populations in all cases from the set K values. The remaining two populations showed similar patterns. The KGP and GOP populations were separated into different clusters when K3 and the KTP and USP populations were represented as a mixture with among. When the K value reached 4, the KGP and GOP populations were still separated into different clusters, but the KTP and USP populations primarily clustered with each other. Burn-in and MCMC repetitions (200,000 times and 1,000,000 times each, respectively) were used to estimate the optimal number of groups (ΔK values) when classifying groups between 2 and 4, which were repeated 10 times by setting the K values from 2 to 4, and the ΔK value was estimated using a structure harvester. The highest ΔK value (1.375%) was obtained from K value 3, which uses allele frequency data calculated with the 14 MS markers and distinguishes the three groups when the group is divided (Supplementary Table S2). This indicates that it is the most appropriate method.