Tracing India’s Canine Heritage Through SNP-based Haplotype Identification

Dog (Canis lupus familiaris) holds a unique position as one of the earliest animals to be tamed, evolving alongside humans over millennia to become one of the most beloved companions in history. With a domestication history spanning over 14,000 years, dogs have been selectively bred for various roles, including as pets, guardians and hunting partners. This deliberate breeding has resulted in the emergence of a diverse array of dog breeds, with 199 recognized internationally by the American Kennel Club (AKC) out of an estimated 340 breeds worldwide, each exhibiting distinct morpho-physiological and behavioral traits (https://www.akc.org/dog-breeds/). The astonishing variety of dog breeds emerging relatively recently can be attributed to factors such as single nucleotide polymorphisms (SNPs) (Ali et al., 2020).
       
Single nucleotide polymorphisms (SNPs) represent variations at a single base pair in the DNA sequence, constituting the most prevalent form of genetic variation within animal genomes. These variations play a crucial role in understanding the genetic diversity of species and are frequently utilized to assess population structures across different organisms (Liu et al., 2022). In this study, a comprehensive collection of 2,18,433 genome-wide SNPs obtained through the ddRAD-GBS approach serves as a valuable resource for deciphering the genetic makeup and differentiation among dog breeds in India.
       
Haplotypes, a group of linked SNPs inherited together on a chromosome, offer additional insight into genetic diversity (Cheng-Li et al., 2023). Understanding the distribution and diversity of haplotypes is integral to unravelling the genetic architecture and population structure of different dog breeds (Vychodilova et al., 2018). On the other hand, population structure, defined as a systematic differentiation in allele frequencies among sub-populations due to non-random mating, is fundamental in tracing the evolutionary history of dog breeds.
       
To overcome the challenges posed by the high cost and time-consuming nature of high-throughput sequencing, this study adopts the genotyping-by-sequencing (GBS) approach. Genome-wide SNP-based diversity and population structure analysis using ddRAD have been carried out in different livestock species like buffalo, yak, horse and camel (Van Dijk et al., 2014; Tezuka et al., 2018; Kumar et al., 2020; Sivalingam et al., 2020). Notably, the ddRAD-GBS technique, proposed by Poland et al., (2012), utilizes two restriction enzymes (PstI and MspI), effectively reducing genomic complexity and enhancing sequencing library uniformity.
       
The Indian dog breeds have received little research attention. Research reports on the molecular characterization of indigenous dogs are very limited. The number of indigenous canine SNPs is still unknown and nearly all of them are still unidentified. Identification of haplotypes, along with analysis of population structure, remains unknown in the owned dog breeds reared in India. Therefore, this study seeks to fill this gap by using advanced genomic approaches to understand breed relationships and to look for marker haplotypes for key genes linked with these differentiated characteristics. The results are useful in elucidating the molecular description of the native dogs and provide valuable information about the molecular genetic distinctiveness of these dogs, which will be of great value in breed conservation and the improvement of the breeding program.

Sample collection and DNA extraction for ddRAD-GBS
 
The genotyping-by-sequencing (GBS) data, previously generated in the Animal Genomics Laboratory, College of Animal Biotechnology, Guru Angad Dev Veterinary and Animal Sciences University, Ludhiana, Punjab, as part of the DBT-funded project entitled “Parentage Determination and Cytogenetic Profiling in Dogs (DBT-19I),” were available. Table 1 presents a brief description of the experimental samples.
       
The experiment was conducted from February 2023 to August 2023 at the College of Animal Biotechnology, Guru Angad Dev Veterinary and Animal Sciences University, Ludhiana, Punjab. Genomic DNA was extracted from the blood samples using the phenol-chloroform-isoamyl alcohol technique (Green and Sambrook, 2018). DNA purity was assessed using the Thermo Scientific™ NanoDrop 2000 spectrophotometer and molecular quality was confirmed via agarose gel electrophoresis. All samples were initially processed for genotyping-by-sequencing (GBS), with an initial QC report received for 50 samples (results of 46 samples were finally received).
 
ddRAD-genotyping by sequencing library preparation and sequencing
 
The samples were subsequently forwarded to Novogene Lifesciences Pvt. Ltd., Mumbai (https://en.novogene.com/), where the GBS analysis (Peterson et al., 2012) was carried out along with mapping the SNPs to the reference dog genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_014441545.1/).
       
ddRAD-GBS is a genomic approach that combines restriction enzyme digestion and next-generation sequencing to identify genetic markers. The GBS DNA library was prepared using the ddRAD approach for each sample by digesting them with two restriction enzymes (MseI and HaeIII_MspI), to reduce genomic complexity and enhance sequencing library uniformity. The resulting fragments were joined by two barcoded adapters that either had a compatible sticky end with universal Illumina P5 or P7 sequences or with the primary digestion enzyme. This process allows for the multiplexing of samples and enables their identification during sequencing, ensuring accurate data collection and analysis on the Illumina platform.
       
After the sequencing assessment report, the raw data for 46 samples were sequenced. The effective sequencing data was aligned with the reference sequence using BWA software (Li and Durbin, 2009). The mapping rate and coverage were then calculated based on the alignment results. SAMtools (Li et al., 2009) was used to handle the BAM files.
 
Filtration of raw reads and SNP identification
 
SNP-calling is the process of identifying and cataloging single nucleotide polymorphisms (SNPs) in a genome. It was performed using SAMtools (Li et al., 2009) to identify individual SNP variations. To minimize the error rate in SNP identification, the following criteria were utilized: (a) There should be more than 4 support reads for each SNP and (b) Each SNP should have a mapping quality (MQ) greater than 20.
       
The raw data is further processed by validating accession numbers with chromosomal details and discarding items missing this information. The samples were then separated into 13 trios using R-coding for Mendelian error evaluation. Errors like tri-allelic genotypes separated by a tab were eliminated. Trio’s paternal genotypes are tested for Mendelian errors. Non-Mendelian SNPs are filtered based on segregation distortion, departing from Mendelian inheritance. The SNP file is then transformed into genotype files and subsequently used for filtering based on minor allele frequency (MAF<0.05).

Population structure analysis using R-programming
 
Filtered VCF files were examined in R v4.3.1 using the vcfR, poppr, ape and RColorBrewer packages. The “dist()” function generates a pairwise genetic distance matrix, reflecting genetic distances between all individuals. Genetic linkages were analysed using a UPGMA tree based on data from Punjab, Haryana and Karnataka (46 samples, threshold 50) (Tomar et al., 2022). For population-wide linkages and multi-loci genotypes (MLG), a minimum spanning network (MSN) was developed utilizing a ‘genlight’ object and distance matrix. Allelic differences were estimated via bitwise.dist() from Poppr.
       
Principle component analysis (PCA) was performed using the glPCA() method on the genlight object and visualization was performed with ggplot2. Discriminant analysis of principal components (DAPC) also employs a ‘genlight’ object with predefined populations, utilising PCA findings for sample assignment. Here, the purpose is to allocate each sample to a population based on the findings of the PCA. Therefore, the same parameters used in the PCA (n.pca = 3 and n.da = 2) were used to rebuild the DAPC. ggplot2 and ‘tidyr’ transform data were used for scatter plot display.
       
For more straightforward comprehension, the graphs were recreated using ggplot2 and categorized samples by their positions. DAPC-calculated population assignments were extracted and a new data frame was built using original population labels and sample names. Using pivot_longer from ‘tidyr’, the data frame was rebuilt to match ggplot2’s structural requirements and the columns were renamed for familiarity. The reconstructed dapc.results data frame was displayed using ggplot2. Samples are on the X-axis, membership probabilities are on the Y-axis and fill color represents original populations. Each facet in the bar plot reflects the initial population assignment, offering an ordered picture comparing membership probability with the original population.
 
Identification and validation of haplotypes
 
Haplotype-based analysis, compared to individual SNPs, may significantly lessen the frequency of false discoveries owing to the reduced number of association tests run (Hamblin and Jannink, 2011). The haplotype identification from the SNP genotyping data was conducted using an R-programming environment in the Ubuntu Linux terminal. The following criteria were followed for the identification of the haplotypes: (a) distance between SNPs: less than 80kb; (b) haplotype length: up to 150kb; and (c) at least 4 SNPs should be included in one haplotype.
       
Three haplotypes were selected for the validation trials, considering SNP density and MAF values of 0.05 or 5%. NCBI (https://www.ncbi.nlm.nih.gov/) was searched for the relevant RefSeq IDs to find related genes with the chosen haplotypes. Three genes, namely AFAP1, GBGT1 and CELSR1, were identified and the primers for these genes were self-designed using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). The details of the primers are shown in Table 2.


 
Sample collection for haplotype validation
 
The Institutional Academic Ethical Committee (IAEC) allowed blood sample collection from dogs, as per Memo No.: IAEC/2023/134-55. Dated 22/05/2023 and Memo No. IAEC/2023/176-186, Dated 14/09/2023. Samples were obtained from Guru Angad Dev Veterinary and Animal Sciences University (GADVASU) multispecialty clinics. Samples from seven distinct breeds were collected using a 21G needle (DISPO VANTM) from the cephalic vein (Table 3). Healthy canine samples were collected, independent of age and sex, into 15-ml sterile tubes containing anticoagulant (EDTA). Samples were kept at -20°C until analysis.


       
Genomic DNA was obtained from the blood samples using the phenol-chloroform-isoamyl alcohol method (Green and Sambrook, 2018). The purity of the DNA was evaluated using the Thermo Scientific™ NanoDrop 2000 spectrophotometer and the molecular quality was verified through agarose gel electrophoresis.
 
Sequencing of PCR results for the detection of SNPs
 
PCR reactions were performed in a 30µl reaction volume using primers for AFAP1, CELSR1 and GBGT1 genes using GoTaq® Colorless Master Mix (Promega), which includes Taq DNA polymerase, dNTPs, MgCl2, reaction buffers and nuclease-free water. The purified PCR products of the GBGT1, CELSR1 and AFAP1 (AF3A and AF3B-fragments) genes underwent bidirectional sequencing (forward and reverse) at GeneSpec Pvt. Ltd., Kerala. Sequencing data were then visualized using FinchTV software (https://digitalworldbiology.com/FinchTV). The obtained FASTA sequences were aligned to reference chromosomal sequences of relevant genes using Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/). All chromatograms were manually inspected for SNP genotyping at each discovered SNP site.

Population Structure Analysis Using PCA
The GBS data files were submitted to NCBI SRA (BioProject: PRJNA843534; SRA accession number SRR26636791 to SRR26636836). The hierarchical distance tree constructed for the forty-six canine samples revealed distinct clustering patterns that align with geographic origins. Branch lengths and confidence values (100%) indicated genetic distances and sample similarities. Breeds from Punjab and Haryana exhibited closer genetic relationships, whereas Karnataka breeds demonstrated less resemblance, underscoring regional differentiation (Fig 1).


       
The Minimum Spanning Network (MSN), developed from reference-based SNP data, elucidated breed-specific genetic variations. Node Id 76, representing the Labrador retriever from Punjab, occupied a central position, signifying ancestral connections with shared SNPs across breeds. Branches radiating from this node highlighted regional diversity, with elongated branches corresponding to specific breeds such as Labrador retriever (21 and 79), Pug (44) and Belgian Malinois (90), affirming the genetic resemblance between Punjab, Haryana and Karnataka breeds (Fig 2).

​
       
Principle component analysis (PCA) further delineated genetic diversity, with PC2 effectively separating Punjab and Karnataka samples from a cluster predominantly including Punjab and Haryana breeds. The Karnataka samples formed a distinct orange ellipse, while the data of Punjab and Haryana exhibited overlapping clusters (Fig 3).  Moreover, DAPC revealed consistent clustering with PCA results, assigning breeds to geographic states (Fig 4). The identification of genetic differences through discriminant analysis emphasized the contribution of breed-specific SNPs to regional clustering.

​

​
       
Membership probability analysis, depicted through colored vertical lines (Fig 5), quantified the likelihood of each sample belonging to a specific region. Visualizing posterior probabilities enhanced the interpretation of genetic affiliations, confirming shared ancestry between breeds from Punjab, Haryana and Karnataka. This shared lineage was supported by SNP genotype similarity, as reflected in the plot of genetic diversity (Fig 6).

​

​
 
Haplotype identification
 
A systematic approach was employed to identify haplotypes using stringent criteria: SNP distance thresholds of 80 kb or cumulative differences exceeding 150 kb, with a minimum of four SNPs constituting a haplotype. This methodology uncovered 15,552 haplotypes from 2,18,429 SNP locations.
 
Isolation and quality check of genomic DNA
 
High-quality genomic DNA was extracted from 21 canine blood (Table 3) samples using the conventional Phenol: Chloroform: Isoamyl Alcohol (25:24:1) technique (Green and Sambrook, 2018) and verified through spectrophotometric analysis (OD 260/280 nm values: 1.6-1.8). The integrity of DNA was confirmed via gel electrophoresis on a 0.8% Agarose gel (70V for 45-50 min.), displaying unambiguous high molecular weight bands under GelDoc system, suitable for downstream applications.


 
Functional analysis of haplotypes
 
Three of the identified haplotypes were prioritized based on SNP frequency and minor allele frequency (MAF > 5%).
· Haplotype 1: Located on chromosome 3, spanning 187 bp with four SNPs associated with the AFAP1 gene.
· Haplotype 2: Spanning 111 bp on chromosome 9 with five SNPs linked to the GBGT1 gene.
· Haplotype 3: Covering 117 bp on chromosome 10 with six SNPs connected to the CELSR1 gene.
 
Identification of Polymorphisms
 
To identify SNPs, the reference sequence of the corresponding chromosomes of the GBGT1, CELSR1 and AFAP1 genes was aligned with the FASTA sequences acquired after sequencing using Clustal Omega. Eleven SNPs were discovered for GBGT1 and CELSR1, whereas no SNP was detected in AFAP1. In the GBGT1 gene, SNPs such as 503 g>A and 510 g>A were prevalent across multiple samples, whereas 516 a>G and 426 g>C were breed-specific (Table 4). Whereas, for the CELSR1 Gene, Six SNPs were identified, with specific variants (e.g., 414 g<A, 469 g>A) displaying consistent presence in specific samples (Table 5).




       
Dogs (Canis lupus familiaris) have been integral companions to humans for thousands of years, fulfilling roles ranging from hunting and protection to emotional and therapeutic support. Domesticated from wolves over 14,000 years ago, dogs are considered the first domesticated animals and have become indispensable members of society (Morell, 1997). Despite their historical and cultural significance, genetic studies on dog breeds reared in India remain limited, focusing primarily on breeds from Western and European origins. This research aims to bridge this gap by analyzing the genetic structure and diversity of significant dog breeds from North and South India, utilizing high-quality SNP data derived from ddRAD-GBS technology.
       
Genotyping-by-sequencing (GBS) has become a powerful tool for obtaining genome-wide SNPs, facilitating genetic diversity and population structure studies. The study identified 8,13,580 SNPs, of which 2,18,433 high-quality SNPs were used to analyze the population structure of divergent dog breeds. Previous studies, such as Kaur et al., (2023), demonstrated the utility of GBS in analyzing genome-wide SNPs in indigenous Indian breeds like Gaddi dogs, identifying over 75,000 high-quality SNPs. This reinforces the role of GBS as a robust tool for population genetic studies and association mapping.
       
The population structure analysis revealed tight ancestral relationships among Indian dog breeds, with geographical proximity often correlating with genetic similarity. These findings align with Zhang et al., (2018), who highlighted that regional factor strongly influence genetic structure in domestic animals. In line with this, the phylogenetic analysis confirmed breed-specific genetic differences shaped by environmental adaptations and historical breeding practices. Similar observations were reported by Freedman et al., (2014), who documented how selection pressures and geography impact breed evolution in dogs globally.
       
A significant outcome of this study was the haplotype analysis, which revealed 15,552 haplotypes. These genome regions are critical for understanding genetic architecture, as they represent blocks of conserved alleles inherited together due to linkage disequilibrium. Three haplotypes associated with AFAP1, GBGT1 and CELSR1 genes were particularly interesting. The AFAP1 gene on chromosome 3 plays a pivotal role in cellular motility and has been linked to physiological adaptations that might benefit specific dog populations (Cunnick et al., 2015).
       
Similarly, GBGT1 and CELSR1, located on chromosomes 9 and 10, are implicated in glycosylation (Indellicato and Trinchera, 2021) and planar cell polarity pathways (Chen et al., 2022), respectively. These pathways are essential for neurological development and cellular signaling, potentially impacting behavioral traits and environmental adaptability. The functional implications of these genes underscore their importance in understanding breed-specific characteristics and genetic diversity.
       
These findings not only provide insights into the genetic structure of Indian dog breeds but also contribute to broader efforts in breed conservation and functional genomics. Identifying specific haplotypes linked to traits such as disease resistance or behavioural attributes could guide breeding programs and improve breed welfare (Gutierrez-Reinoso et al., 2021). Furthermore, this study complements global efforts to map dog genetic diversity, as highlighted in studies of conserved haplotypes in humans and animals (Guryev et al., 2006).

The application of ddRAD-GBS technology has significantly advanced our understanding of the genetic variation within Indian dog breeds, revealing unique haplotypes associated with critical physiological functions. This study provides valuable insights into the genetic composition, specific genetic variants (SNPs) and distinctive DNA sequences (haplotypes) in prominent canine breeds in North and South India. The extensive SNP variants identified enrich our knowledge of genetic diversity in these breeds and serve as a foundational resource for future research. The findings hold practical implications for breeding programs, enabling breeders to make informed decisions to preserve and enhance indigenous dog populations’ genetic health and diversity. Future research should focus on expanding sampling and incorporating advanced genomic tools, such as transcriptomics, epigenomics and functional genomics, to explore the genetic underpinnings of breed-specific traits further. These efforts will contribute to the conservation and sustainable management of India’s diverse and culturally significant canine breeds, bridging the gap in genetic studies on indigenous populations.

The authors thankfully acknowledge the funding the Department of Biotechnology, Government of India, provided through the collaborative research project “Parentage Determination and Cytogenetic Profiling in Dogs (DBT-19I)”.
 
Authorship
 
1) DS: Data analysis, blood sample collection, wet lab analysis for validation studies and manuscript writing.
2) SKM: Provided with the blood samples for validation studies and checked the manuscript.
3) NK: Manuscript editing.
4) CSM: Designing, supervising and editing the manuscript.

All authors read and approved the manuscript. All authors contributed to the manuscript revision and read and approved the submitted version.

The authors declare no conflict of interest.