Peripheral blood samples from six horses (2H, 7H, 16H, 25H, 30H, and 57H) were collected in 1 mL volumes at “Akzhar Ondiris” horse farm (51°32′07.4″N 77°27′16.9″E) in Pavlodar region of Kazakhstan (Figure 1A). All samples were anticoagulated with EDTA and refrigerated at 4°C. The phenotypic characteristics of these horses are detailed in Supplementary Table S1. Genomic DNA was extracted from the samples using Illustra Blood Kit (Cytiva, United State) and Gentra Puregene Blood Kit (Qiagen, Germany) following the manufacturers’ protocols. The concentration and quality of the extracted DNA were checked using a Qubit fluorometer (Invitrogen, United State), a Nanodrop 2000 spectrophotometer (Thermo Scientific, United State), and 1% agarose gel electrophoresis. This high-molecular-weight DNA was then used for library construction and subsequent Nanopore sequencing.

To generate Oxford Nanopore long reads, 3 µg of genomic DNA was randomly sheared to obtain a target size of 20 kbp using g-TUBE (Covaris, United State) and processed according to the Ligation Sequencing Kit (SQK-LSK110) protocol (Oxford Nanopore Technologies, United Kingdom). For genome sequencing, at least 1 µg of sheared DNA from each sample was utilized for library construction. DNA fragments were repaired using NEBNext FFPE Repair Mix (New England Biolabs, United State). End repair and A-tailing were performed using the NEBNext End Repair/dA-Tailing Module kit (New England Biolabs, United State), followed by ligation of Oxford Nanopore sequencing adapters with the NEBNext Quick Ligation Module (E6056) (New England Biolabs, United State). The constructed libraries were sequenced on R9.4.1 flow cells of PromethION sequencer (Oxford Nanopore Technologies, United Kingdom) for 72 h. Basecalling of the raw signal data was performed using Guppy v.5.1.13, which also trimmed adapters and removed low-quality sequencing reads with a Q-score below 9.0. All DNA samples were sequenced with an average coverage of 26X. A summary of the sequenced reads is provided in Supplementary Table S2.

Draft assemblies were produced using one round of Flye v.2.9.2 (Kolmogorov et al., 2019), followed by a polishing round with Oxford Nanopore Technologies (ONT) reads using Medaka v.1.11.1 (https://github.com/nanoporetech/medaka). To evaluate the quality of the final assemblies, we aligned the ONT contigs to EquCab3.0 reference genome assembly (NCBI Accession No. GCF_002863925.1) and assessed them with QUAST v.5.2.0 (Gurevich et al., 2013). Considering the advanced sequencing ability of ONT, the longest contig among the assembled genomes was 92.32 Mb, and the largest contig N50 was 28.26 Mb. The completeness of the genome assemblies was further assessed using BUSCO v.5.4.6 (Simão et al., 2015), which compared the genome against the laurasiatheria_odb10 database containing 12,234 orthologous genes. BUSCO assessment scores ranged from 93% to 95% (Figure 1B; Table 1), indicating high completeness for the obtained assemblies.

To identify SNVs and indels, the wf-human-variation Epi2Me Labs pipeline from ONT (https://github.com/epi2me-labs/wf-human-variation) was used. Samples were analyzed using Clair3 v.1.0.4, which identified small variants in ONT reads (Supplementary Table S3). The number of identified SNVs ranged from 6,336,129 to 7,101,556, while the number of identified indels ranged between 549,718 and 820,662 across samples.

Phylogenetic analysis and tree construction were performed using VCF-kit v.0.2.6 (Cook and Andersen, 2017) and MEGA software v.11.0.13 (Tamura et al., 2021). The neighbor-joining tree was constructed using 1,331,674 mutation points from a merged VCF file (Figure 1C) containing data from Kazakh horses and 88 additional horse samples (Jagannathan et al., 2019) deposited in the European Nucleotide Archive (ENA) database (https://www.ebi.ac.uk/ena/). At the autosomal genetic level, Kazakh Zhabe horses formed a distinct cluster and a separate group compared to other horse breeds. Additionally, a multiple sequence alignment of all mitochondrial D-loop sequences was performed in MEGA using the built-in MUSCLE (Edgar, 2004) alignment option to construct a consensus tree. The analysis included 71 samples of the control region and mtDNA from our assemblies, as well as 25 different horse breeds deposited in the National Center for Biotechnology Information (NCBI) GenBank database (http://www.ncbi.nlm.nih.gov/). All sequences were processed using blastn v.2.12.0+ (Camacho et al., 2009) to extract an early part of the control region (400 bp in the position between 15,469 and 15,868). The consensus tree (Figure 1D) was built using the Neighbor-Joining method (Saitou and Nei, 1987) with 1,000 bootstrap iterations. The D-loop region sequence of the donkey (Equus asinus, NCBI Accession No. NC001788) was used as an outgroup. While phylogeny reconstruction showed Kazakh horse mtDNA sequences are widespread and distributed across many different clusters in the tree, two samples (Kazakh1 and Kazakh5) from the assembled genomes formed a distinct clade with Cheju and Akhal-Teke horses. These results are consistent with previous studies (Gemingguli et al., 2016) reporting tightly linked mtDNA genetic relationships between these breeds. It can be suggested that the Kazakh horse breed has a mixed origin in the maternal lineage, likely due to the use of horse populations in trade and military campaigns, which moved them to distant locations, where they interbred with indigenous populations. The observed lack of Kazakh horse samples clustering in the phylogenetic tree constructed from mtDNA control region sequences may indicate high levels of variability, which, in turn, means that the Kazakh breed may serve as an important reservoir of genetic biodiversity. It is of particular significance for horses, as a species, because its wild ancestors are now extinct and sources of biodiversity that could be used to maintain their functions in certain environments are limited.