Genome-Wide Identification, Annotation of SNPs, and Selective Sweeps in Indian Horse and Pony Breeds

1. Introduction

Horses have been used for sporting, hunting, and armed warfare for many generations. Therefore, certain athletic performance-related genes were positively selected in response to the demand created by the environment, resulting in their unique physical morphology and physiological system. Identifying the genetic variants and positive selection regions will help in understanding their role, which may improve selection based on environmental adaptation and performance. India is bestowed with rich equine germplasm. However, these animals have not been studied in detail on the genomics diversity scale. India possesses about 0.54 million equines as per the 2019 livestock census [1]. The major population of equids comprises donkeys, mules, horses, and ponies. These provide livelihood means to the rural and semi-urban societies living in arid, semi-arid, and hilly regions, especially within the foothills of Himalayas, through transport and drought whereas the remaining small population of equines is employed in the army, police, border security force, racing industry and equestrian sports [2,3]. The population of horses and ponies were reduced by 1.81% per year from 1956 to 2019. A negative annual growth rate (-3.28%) of equids for draught work over the years 1956 to 2019 (2.64 to 0.54 million) indicates farmers’ preference for a faster mode of mechanised transport [1]. India has three horses (Kathiawari, Kachchhi-Sindhi, and Marwari) and four ponies (Manipuri, Spiti, Zanskari, and Bhutia) inhabiting different agro-climatic regions [4,5,6,7,8,9]. Over time, these breeds have adapted with certain unique traits like endurance, relative disease tolerance, sturdiness, and sure-footedness to sustain and work in harsh environmental conditions. Among horse breeds, Marwari is the best-adapted breed in hot, dry, and desert conditions of the Marwar area of Rajasthan. The animals of this breed are selected for their conformation, speed, and stamina, endurance capacity with good export potential, while Kathiawari, the native breed of Kathiawar area in Gujarat, is also known for its speed, stamina, endurance, and capability of fast movement. In India, seven breeds of horses have been identified based on their geographical localisation [4,5,6,7,8]. Different attempts have been made for characterization of horses, ponies and donkey breeds at the phenotypic, biochemical and molecular levels [4,5,6,7,8,9,10,11,12,13,14], but still, it is difficult to assign a breed to an unknown sample In this study we will characterize the genetic diversity of unique and valuable indigenous horses and quantify their relationship with other breeds, aiming to identify breed-specific variants underlying unique phenotypic traits such as inward turning of ears, endurance and racing performance.

Since the domestication of the horse approximately 5,000 years ago, selective breeding has been directed mainly towards the use of the horse in agriculture, transportation, and warfare. Within the past 400 years, the founding of formal breed registries and continued breed specialisation have focused more upon preserving and improving traits related to aesthetics and performance. As a result, most horse breeds today are closed populations with high phenotypic and genetic uniformity of individuals within the breed, but with a great deal of variation among breeds. The use of genetic markers now allows for the identification of genomic variation specific to breeds to be uncovered [15], and to trace genomic relationships and pedigrees [16,17,18]. Genetic markers (i.e. single-nucleotide polymorphisms; SNPs) are an additional source of information that could be incorporated into the breeding evaluation to enable accurate early selection of horses for coat colour, health, and performance traits. High-throughput, whole-genome SNP arrays are now available for these purposes beyond that possible from conventional microsatellite markers. Other newer approaches include the use of genotyping-by-sequencing (GBS), RAD-seq, and whole-genome sequencing to characterise genomic diversity among individuals and within breeds. Once genomic regions targeted by selection are detected, the variants and processes that have contributed to desired phenotypes within breeds and across performance groups can be more readily identified. Population-based approaches to identify signals of selection using the loss of heterozygosity and/or other diversity indices have been successful in several domestic species [19,20]. GBS offers great potential for characterising molecular variants and linking these variants with basic phenotypes and other desired characteristics in livestock populations because it makes it possible to cover large fractions of the genome and to vary the sequence read depth per individual. The ddRAD-Seq, a popular tool for MAS in plants, is a simple, reproducible, highly multiplexed approach based on the Illumina sequencing platform [21]. The method has been proven suitable for population studies, germplasm characterisation, genetic improvement, and trait mapping in a variety of diverse organisms. The major advantages over other protocols are both technical simplicity and that informatics pipelines are publicly available. Very few works with MAS have been reported in livestock such as cattle [22], Yak [23], chicken [24] and ducks [25] and horses [26].

2. Materials and Methods

Location of the Study and Source of Samples

The blood samples of adult, healthy and true-to-breed equines of each breed were selected based on their morphological features and breeders' information, from their breeding tracts in different parts of India (Figure 1).

These selected animals either belonged to government organizations and/or private small or big equine breeders maintaining these animals for breeding purposes. About 19 Kathiawari horses were sampled from state equine farms at Junagarh and Inaz, Police horses from Junagarh, Rajkot, Surender Nagar as well as some private breeders in other parts of Kathiawar (Gujarat), Marwari horses (N=26) from private breeders of Jodhpur, Pali Marwar, Udaipur, Dundlod, Nawalgarh, Jalore, Nagore, Bikaner areas (Rajasthan), the Kachchhi-Sindhi (n=24) were selected based on their morphological features in their home tracts in Bhuj (Gujarat, India) and Jaisalmer (Rajasthan, India). The selected animals belonged to the private horse breeders maintaining the animals for breeding purposes. Efforts were made to select animals from different parts of their home tract to include adequate representatives of the breed. The Manipuri ponies (N=5) and Zanskari ponies (N=5) from Equine Production Centre (EPC), ICAR-NRCE, Bikaner (Rajasthan) (the animals being procured earlier from their native breeding tracts and now maintained at Bikaner). Besides this, Thoroughbred horses (N=16) bred in India and available in different parts of Haryana and Rajasthan states were also included as an outgroup in this study. Efforts were made to cover the entire breeding tract of each horse and pony breed in India. The blood sampling of animals was performed following the relevant guidelines and regulations as approved by the Institutional Animal Ethics Committee (IAEC).

Blood Collection and DNA Isolation

About 8-10 ml of blood from the jugular vein of true-to-breed individuals of the horse/pony breeds were collected aseptically into vacutainers coated with EDTA (0.5 mM, pH 8.0). DNA was extracted from whole blood samples according to the kit protocol (ReliaPrep™ Blood gDNA Miniprep System, Promega).

Quality of Genomic DNA

DNA quality was assessed in 0.8 % Agarose w/v suspension in 0.5X TBE buffer (pH 8) as described by [27]. DNA samples devoid of smearing were considered as good. The Qubit dsDNA HS (high sensitivity) assay kit (Thermo Fisher Scientific) was used to quantify DNA in the Qubit 3.0 fluorometer (Thermo Fisher Scientific) following the manufacturer's instructions.

ddRAD-Based Genotyping Library Preparation

After the initial quality and quantity check of DNA, standard RAD protocol as described by [15] was used for further sample processing. Double digestion of DNA with Sph I and MluCI restriction enzymes was carried out followed by combinatorial barcoding including both Illumina index and an inline barcode for library preparation. Ligation of P1 and P2 adaptors was carried out using T4 DNA ligase, followed by pooling and clean-up of the ligated product. Size selection of the digested DNA product was done after running 2% agarose gel electrophoresis. Then, PCR amplification was done to enrich and assign the Illumina-specific adapters and flow cell annealing sequences were carried out. Further, all the samples were pooled after adapter ligation and size selection was carried out. Finally, the samples were sequenced on Illumina Hi-seq 2000 and bioinformatic analysis was carried out (Figure 2).

SNP Identification and Quality Control for Diversity Analysis

Paired-end fastq files were trimmed for barcodes using PRINSEQ [28]. Low-quality reads with an average PHRED score of < 15 per read were discarded using STACKS 2.2 [29]. After the initial quality control, the reads were aligned with the Equus caballus reference genome (EquCab-3.0) using Bowtie 2 [30]. SNPs and Indels were filtered using VCFtools [31]. SNPs with good quality (Q = 30) at read depth 10 were retained for further analysis. Also, the SNPs with missing genotypes, minor allele frequency (MAF) <0.05, and those significantly (p <0.0001) deviating from the Hardy Weinberg Equilibrium were further removed by VCFtools.

Selection Sweep Identification

All breeds combined SNPs annotated vcf file at RD10 was used as input in SweeD [32] to find common selection signatures among indigenous and thoroughbred horse breeds. Later chromosome-wise selected regions were screened for genes using UCSC Genome Table Browser [33] after figuring out genes present in selected regions were further subjected to functional annotation using Panther GO [34]. Gene to Gene network interaction carried out using Cytoscape [35].

3. Results

Around 252.60 million raw reads from 95 samples of several horse breeds were processed, and of this, 248.08 million reads were of good quality after trimming out adaptors and removal of bad quality sequences. About 122.25 million reads (i.e. 49.20 %) of total processed reads were aligned with the reference genome EquCab-3.0 with an overall alignment rate of 92.32 %. A total of 84.07 million reads were found to be uniquely aligned among the total aligned reads, accounting for 68.11%. SNP and Indels were filtered both as combined and breed-wise horses, respectively, at many read depths (RD @ 2, 5 and 10) (Table 1), and chromosome-wise distribution of genetic variants was found for a combined dataset of Indian horses and ponies (Figure 3).

3.1. Genome-Wide Annotation of SNPs

Filtered SNPs of combined horses and ponies at a read depth of 10 with quality above 30 were subjected to annotation using a database created using the EquCab 3.0 reference genome. About 38.30 % and 37.95 % of total variants account for transcript and intron regions, followed by intergenic regions with 11.28% of variants (Figure 4). The substitution ratio of transitions (Ts) is 68.35%, and transversions (Tv) is 31.75 %, with a ratio of 2.16 transitions per transversion.

3.2. Gene-Wise Mapping of SNPs and Indels

A total of 460 high-quality SNPs and 328 Indels were mapped with 30 already reported candidate genes that have many functions, and their role is discussed in Table 2.

3.3. Selection Sweeps

A total of 1631 selection signatures were enumerated at the top 1 percentile with a greater than 1.72 likelihood score from 1,77,726 filtered SNPs. 700 selected regions were found incorporated in 306 genes; of these, 125 genes have greater than one selected region, and the top 10 genes with their corresponding selected regions were listed in Table 3.

On functional annotation, the majority of genes, about 28 %, accounted for the cellular process (GO: 0009987) followed by biological regulation (GO: 0022610) and signalling (GO: 0023052), with 15.8% & 15.3% respectively of total annotated genes (Figure 5, panel A). At the molecular level, nearly half of the genes aid in binding function (GO: 0005488) and act as function regulators (GO: 0098772) (Figure 5, panel B). Annotated genes in greater proportion assist in integrity, cellular anatomical entity (GO: 0110165), and intracellular (GO: 0005622) cellular components (Figure 5, Panel C).

Pathway analysis of selected region genes associated in several functional pathways revealed their roles in many vital functions related to growth, regulation of normal cell function, immunity, cardiac conduction and its regulation, neuron development and conduction of signal, hormonal regulation, transportation of mineral ion, and many signalling pathways (Table 4).

3.4. Gene-Gene Interaction Network Analysis

On exploring gene interaction among selected region genes, many genes were found as hub genes involved in multiple biological roles that also may account for the physiological adaptation of the horse towards its high performance. Notably, the ANK2, BCL2, and RYR2 were involved in the transportation of calcium ions and cardiac regulation (Figure 6).

4. Discussion

This study aimed to understand the genetic architecture of Indian horses and ponies through a genome-wide selection signature scan. Earlier, to understand the genetic variants in horses, CNVs from six Indian horse breeds, namely, Manipuri, Zanskari, Bhutia, Spiti, Kathiawari and Marwari, were discovered using a genotyping array [36]. Genome-wide single-nucleotide polymorphism-based genomic diversity and runs of homozygosity (ROH) were also elucidated for selection signatures in Indian equine Breeds [37]. The results of the presented study gave an idea of how the putative positive selected genes involved in the physiological regulation of horses in terms of cardiac, neuronal and hormonal. CADM1 and SLC4A10 genes have a greater number of selected regions. Cell adhesion molecule-1 (CADM1), an immune gene that assists in the regulation of cell-to-cell adhesion in diverse human epithelial cells and consequently provides defence as an anti-malignant agent, which inactivates a variety of cancers such as small cell lung cancer and breast cancer in humans [38,39]. Numerically high selected regions in CADM1 genes gave a glimpse that its selection might be a reason for normal lung function regulation without malignancy in horses, even after their long-term high stress involves performance. The solute carrier 4 family (SLC4A10) gene is involved in the transportation of sodium-coupled bicarbonate ions, mainly in the central nervous system and a minor proportion in the ocular receptors [40]. In a study involving positive selected regions over a variety of horse breeds, this gene is found to be common among all relating to the function of regulation of synapse assembly along with SPARC (Osteonectin), Cwcv, and Kazal Like Domains Proteoglycan 1 (SPOCK1) which resembles our findings and it was hypothesized that these genes might be involved in controlling ocular and hearing processes in horses [41].

Ryanodine Receptor type (RYR2) is a candidate gene for athletic performance, predominantly involved in the Ca++ release channel in cardiac muscle [42]. In the mice, the phosphorylated protein encoded by RYR2 increases the inflow of intracellular calcium in myocytes, resulting in cardiac contractility and ventricular arrhythmia during exercise [43]. Researchers reported that the RYR2 gene is closely associated with the QTL region of jumping performance in sport horses [44]. ANK2 (Ankyrin 2) is a vital gene that acts in the regulation of ion channels and transporters to cardiac cell membranes and is also involved in the organisation of cellular receptors such as Na+/ K+ ATPase, the Na+/Ca++ exchanger, and inositol-1,4,5-trisphosphate receptors [45]. ANK2 causes a longer QT interval in the cardiac cycle via an increase in contractility by the influx of calcium in the sarcoplasmic reticulum, resulting in bradycardia in humans. It was observed that a positive selection signature in this gene might be a reason for the lower heartbeat to compensate for its high performance [46]. B-cell lymphoma 2 (Bcl-2) is an anti-apoptotic survival protein in several cells, including cardiac myocytes, against external stress factors such as heavy exercise, training, and toxic drugs, causing cardiac ischemia resulting in cardiac failure [47,48]. This gene has a closer association in the regulation of integrity and cell survival of myocytes through influencing mitochondrial membrane permeability and also by activation of caspases [49]. GRM8 (glutamate receptor, metabotropic 8) and GRIK3 (glutamate receptor, ionotropic, kainate 2) are both genes encoding for glutamate receptor regulation and activation of L-glutamate, a major excitatory neurotransmitter in the central nervous system [50]. Earlier works [51,52,53] idealised the glutamate involved in nervous development, synaptic assembly, behavioural and learning processes, which might influence athletic success and the usefulness of the horse. The Vitamin D receptor (VDR) gene is involved in vitamin D3 metabolism and regulation of a variety of other metabolic pathways, such as those involved in the immune response and muscle strength [54]. A study [55] found that a positive association between vitamin D and muscle strength also improves athletic performance. Tenascin C (TNC) is a glycoprotein abundant in developing tendons, bone, and cartilage, and also commands the interaction of cell-matrix [56]. [57] found tenascin coding genes (i.e., TNC, TNN, TNR) have positive-selected regions which are associated with performance and tendon pathology in thoroughbred horses.

The study provided a genome-wide scan of selection signatures in Indian horses and ponies, significantly supplementing our prior work on CNVs and ROH-based diversity in these breeds. Our identification of 1,631 positive selection signatures—700 mapped to 306 genes—reveals a genetic architecture fine-tuned for performance adaptation. Three key functional clusters emerge, with profound evolutionary and translational implications. These pertain to cardiovascular-exercise Adaptations, neural and sensory Specialization and structural integrity.

5. Conclusions

These findings contextualize Indian equine breeds as reservoirs of adaptive genetic variants shaped by diverse environmental pressures. The convergence of signatures in cardiac, neural, and structural genes across breeds suggests universal selection pressures for performance. Beyond breeding applications (e.g., selecting for TNC variants to reduce tendon injuries), our data hold cross-species relevance: CADM1’s potential role in cancer resistance warrants exploration in comparative oncology, while ANK2/BCL2 insights could inform human cardiovascular research. Future work should validate these signatures in functional assays (e.g., RYR2 calcium flux studies), assess genotype-phenotype links in athletic cohorts, and expand the sampling to underrepresented breeds to decode local adaptations. Such efforts will maximize the translational potential of this genetic blueprint for enhancing equine welfare and performance globally.