Pattern of Runs of Homozygosity and Selection Signatures in Nili Ravi Buffalo

Atiq Ahmad; Hamid Mustafa; Waqas Ahmad Khan; Asif Rafiq; Ishtiaque Ahmad; Haiba Kaul; Seema Baloch; Abdul Manan; Waqas Akram

doi:10.20944/preprints202602.0650.v1

Submitted:

08 February 2026

Posted:

09 February 2026

You are already at the latest version

Abstract

The present study analyzed the genomic landscape of Nili Ravi buffalo using the Axiom 90K SNP chip to assess SNP density, Runs of Homozygosity (ROH), and genomic differentiation. SNP density analysis revealed substantial variation across chromosomes, with Chromosome 1 exhibiting the highest density (>50 SNPs per 1Mb window), while Chr7, Chr13, and Chr20 had lower densities. High-density SNP regions were identified as potential genomic hotspots under selection, whereas regions with sparse SNP coverage indicated gaps in genetic diversity. ROH analysis identified four ROH length categories, with medium-length ROH (2–4 Mb and 4–8 Mb) constituting 68.9% of total ROHs, indicating a moderate inbreeding history. Chromosome 5 showed the highest short ROH occurrences (1–2 Mb), while Chromosome 3 had the longest ROHs (>16 Mb). The highest ROH coverage was observed on Chromosomes 5 (10%), 3 (9.9%), and 19 (8.8%). Eight genomic regions with high ROH frequencies on chromosomes 5, 11, and 19 contained candidate genes related to adaptability, productivity, and disease resistance. Key genes identified included COL8A1, ACACA, SREBF1, AKT3, and WNT11, linked to milk production, fat metabolism, and fertility traits. Genomic differentiation analysis using FST values (mean = 0.682) between Nili Ravi and Kundi buffalo populations identified 24 genomic regions under selection, with Chromosome 11 harboring the highest-ranked SNP (FST = 0.740787) associated with ARPP19. Other key candidate genes included CSMD1, ROR1, DNAJC15, and FBXW7, involved in disease resistance, metabolism, and reproductive traits. These findings provide crucial insights for genetic selection programs aimed at improving economically significant traits in Nili Ravi buffalo.

Keywords:

Nili Ravi buffalo

;

SNP density

;

runs of homozygosity

;

FST

;

genetic diversity

;

genomic selection

Subject:

Biology and Life Sciences - Animal Science, Veterinary Science and Zoology

Introduction

Buffalo (Bubalus bubalis) is an important livestock species that significantly contributes to the dairy and meat industries worldwide. It possesses unique physiological and genetic traits that make it highly adaptable to diverse environmental conditions, particularly in tropical and subtropical regions. The superior milk quality of buffalo, characterized by higher fat (6–9%), protein (4%), and lactose (5%) content compared to cattle, makes it a preferred dairy species in several countries, including Pakistan, India, Italy, and Brazil [1]. Additionally, buffalo meat is considered a healthier alternative to beef due to its lower cholesterol (46 mg per 100g) and saturated fatty acid content (0.46 g per 100g) compared to cattle (59 mg and 4.33 g, respectively) [2]. In Pakistan, buffaloes account for a significant proportion of the agricultural GDP, the Nili Ravi breed is regarded as the backbone of the dairy sector, often referred to as the "Black Gold" due to its high milk productivity and adaptability [3]. However, genetic challenges such as inbreeding and loss of genetic diversity pose significant threats to the long-term sustainability of Nili Ravi buffalo farming.

Inbreeding has been a critical concern in animal breeding programs, as it can lead to inbreeding depression, reduced genetic variation, and decreased fitness of the population [4]. The presence of excessive inbreeding within a breed results in reduced reproductive performance, lower disease resistance, and decreased milk yield, ultimately affecting economic sustainability. Studies indicate that the Nili Ravi buffalo population has an inbreeding coefficient of approximately 15%, with an estimated inbreeding rate of 1% per generation [5,6]. The prevalence of inbreeding in Pakistani buffalo populations necessitates the development of strategies for effective genetic management to maintain the breed’s productivity and adaptability. With advancements in genomic tools, particularly high-density single nucleotide polymorphism (SNP) arrays, inbreeding was assessed through Runs of Homozygosity (ROH) analysis [7]. This method was applied to directly estimate autozygosity, enabling breeders to implement genomic selection strategies to control inbreeding and improve desirable traits.

The application of genomic technologies has been instrumental in revolutionizing animal breeding by providing insights into the genetic architecture of economically significant traits [8]. The availability of high-density SNP panels, such as the Axiom 90K SNP chip, has facilitated the identification of genetic markers associated with milk yield, fertility, disease resistance, and environmental adaptability in buffalo [9]. In addition to ROH analysis, genome-wide selection signature studies have identified genomic regions under natural or artificial selection, helping to pinpoint genes linked to productivity and resilience traits [10]. Selection signatures in buffalo have been associated with key metabolic, reproductive, and immune functions, which are crucial for optimizing breeding programs [11]. By leveraging these genomic resources, breeders will able to develop data-driven selection strategies to enhance the genetic potential of Nili Ravi buffalo while ensuring the preservation of its genetic diversity.

This study was aimed at exploring the genetic architecture of the Nili Ravi buffalo by analyzing SNP density, ROH patterns, and selection signatures using Axiom 90K SNP chip genotyping. The findings have provided valuable insights into the genetic structure of this breed, highlighting the distribution of SNPs, the extent of inbreeding, and genomic regions under selection. The results will contribute to the sustainable genetic improvement of Nili Ravi buffalo, ensuring its productivity and resilience. By integrating genomic tools into breeding programs, this research is concluded to provide a foundation for data-driven selection strategies, ultimately ensuring the long-term viability of this economically significant buffalo breed in Pakistan.

Materials and Methods

Ethical Statement

This study used previously generated genotype data obtained from published sources and a public repository. No new biological sampling, animal handling, or experimental interventions were conducted for the present analyses. Ethical approvals for the original sample collection were secured in the respective primary studies and/or institutional frameworks under which the datasets were generated. Public data were accessed and used in accordance with repository terms and conditions; therefore, additional ethical approval was not required for this secondary analysis.

Genotype Datasets and Populations

Genome-wide SNP genotypes generated using the Axiom 90K buffalo array were assembled from three independent sources. First, genotypes for 35 buffalo were obtained from the HEC-NRPU (No. 16844) funded project. Second, genotypes for 35 animals were incorporated from a previously reported dataset (Irem, 2022). Third, publicly available genotypes for 15 Nili Ravi and 10 Kundi buffalo were downloaded from Harvard Dataverse (2024; https://doi.org/10.7910/DVN/UGA2QX).

Before downstream analyses, datasets were harmonized for SNP identifiers, chromosome labels, genomic coordinates, and allele coding. Variants with unresolved allele/strand inconsistencies across datasets were excluded to avoid spurious ROH calls or inflated differentiation signals.

Quality Control and Marker Filtering

Quality control (QC) was applied at the marker level to retain high-confidence SNPs for ROH and selection analyses. SNPs were removed if they met any of the following criteria: (i) minor allele frequency (MAF) < 0.05, (ii) genotype call rate < 0.90, or (iii) significant deviation from Hardy–Weinberg equilibrium (p < 0.001). Only autosomal SNPs passing QC were retained, and sex chromosomes were excluded to prevent sex-linked marker behavior and hemizygosity from biasing homozygosity-based metrics. After QC, 43,389 autosomal SNPs were available for analyses.

SNP Density Profiling

To evaluate marker coverage across the genome, SNP density was summarized using non-overlapping 1-Mb autosomal windows. The number of retained SNPs per window was calculated and visualized as a genome-wide heatmap. This step provided a diagnostic view of marker distribution and identified chromosomal segments with sparse or dense SNP coverage that could influence ROH detection power and regional signal interpretation.

Runs of Homozygosity Detection and Genomic Inbreeding (FROH)

Runs of homozygosity (ROH) were used to quantify autozygosity and to infer the genomic footprint of inbreeding. ROHs were detected with PLINK v1.9 using autosomal SNPs and the following thresholds: minimum ROH length of 1 Mb, at least 15 SNPs per ROH, maximum inter-SNP gap of 1,000 kb, allowance of one heterozygous genotype within a segment, and a minimum marker density of one SNP per 50 kb.

ROHs were classified into four length categories—1–4 Mb (short), 4–8 Mb (medium), 8–16 Mb (long), and >16 Mb (very long)—to reflect different time scales of shared ancestry (shorter ROH generally representing older haplotypes and longer ROH reflecting more recent inbreeding). For each animal, ROH burden was summarized by the number of ROHs, total ROH length, and chromosomal distribution.

ROH Islands and Region Prioritization

To highlight genomic segments repeatedly affected by homozygosity, ROH “islands” were defined as regions with elevated ROH occurrence across individuals. Autosomes were scanned to identify segments with high ROH frequency, and these regions were retained for downstream gene annotation and biological interpretation.

Population Differentiation and Selection Signatures (FST)

To detect genomic regions potentially shaped by divergent selection between Nili Ravi and Kundi buffalo, pairwise population differentiation was quantified using Wright’s fixation index (FST), estimated with the Weir and Cockerham approach. SNPs exceeding a genome-wide cutoff of FST > 0.65 were considered highly differentiated candidates. Genome-wide FST patterns were visualized using Manhattan plots, and chromosome-level plots were produced for chromosomes carrying peak signals.

Candidate Gene Mapping and Functional Interpretation

Genomic regions highlighted by ROH islands and elevated FST values were mapped to annotated genes using the UOA_WB_1 buffalo reference genome (GCF_003121395.1). Genes overlapping or located closest to prioritized regions were compiled as positional candidates. Functional interpretation was supported using AnimalQTLdb and targeted literature, focusing on traits relevant to buffalo productivity and adaptation, including milk yield and composition, fertility, metabolic efficiency, immune response, and environmental resilience.

Parameter Justification and Rationale

ROH and selection analyses are sensitive to marker density, genotyping quality, and threshold choices; therefore, parameters were selected to balance sensitivity and specificity for a medium-density SNP array. A minimum ROH length of 1 Mb was used to reduce the likelihood of calling short ROH that can arise by chance under limited marker density, while still retaining informative segments for inbreeding inference. Requiring ≥15 SNPs per ROH and enforcing a minimum density of one SNP per 50 kb helped ensure that ROH segments were supported by sufficient marker evidence, minimizing false positives in regions of low SNP coverage. Allowing one heterozygous genotype per ROH accommodated occasional genotyping errors while preserving long homozygous tracts. The maximum inter-SNP gap of 1,000 kb prevented merging of distant homozygous stretches across poorly covered intervals.

For population differentiation, the Weir and Cockerham estimator was used because it is widely applied in genome-wide SNP datasets and performs well under unequal sample sizes. The threshold of FST > 0.65 was applied to prioritize the most strongly differentiated loci as candidate selection signals, and visualization with Manhattan plots facilitated identification of chromosomal clusters consistent with selection rather than isolated outliers.

Software and Computational Environment

All analyses were conducted on a Linux-based high-performance computing (HPC) system. PLINK v1.9 was used for genotype QC, ROH detection, and FST estimation. Data summaries and plots were generated in R v4.2.0 using ggplot2 and qqman. SNP density plots were created using SRplot. Genome annotation and QTL cross-referencing were supported using the UOA_WB_1 reference genome resources and AnimalQTLdb.

Results

SNP Density Analysis of Nili Ravi Buffalo Genome

Figure 1 shows the SNP density analysis of the Nili Ravi buffalo genome was conducted using the Axiom 90K SNP chip to evaluate the distribution of SNPs across all chromosomes within 1Mb window sizes. The results demonstrated substantial variation in SNP density among chromosomes, with certain regions exhibiting high concentrations of SNPs while others displayed sparse coverage. Chromosome 1 showed the highest SNP density, with several windows containing over 50 SNPs, whereas chromosomes such as Chr7, Chr13, and Chr20 exhibited lower SNP densities. Regions of high SNP density, depicted in red on the heatmap, were particularly evident on chromosomes like Chr1, Chr4, and Chr14, highlighting areas with up to 50 SNPs per 1Mb window. These regions may represent genomic hotspots under selection or with elevated mutation rates.

In contrast, chromosomes such as Chr6 and Chr19 revealed regions with sparse SNP coverage, containing only 0 to 5 SNPs per 1Mb window. Such low-density regions may indicate areas with lower genetic diversity or gaps in SNP coverage. Despite these variations, the Axiom 90K SNP chip provided a reasonably uniform SNP distribution across most chromosomes, with a significant proportion of the genome covered by windows containing 10–20 SNPs, as represented in green and yellow on the heatmap. High SNP density regions identified in this analysis are crucial for genome-wide association studies (GWAS) and for detecting quantitative trait loci (QTLs) associated with economically significant traits. However, regions with low SNP densities highlight the necessity for further marker development or sequencing to improve genome coverage.

Distribution of Runs of Homozygosity

The different classes of ROH are present across all chromosomes of the Nili Ravi buffalo genome (Figure 1). Chromosome 5 exhibits the highest number of short ROHs (1–2 Mb), while Chromosome 3 displays the highest number of long ROHs (>16 Mb). The entire ROH length for Nili Ravi buffalo predominantly consists of medium-length segments (2–4 Mb and 4–8 Mb), collectively accounting for approximately 68.9% of all detected ROHs (Table 1).

In contrast, the larger ROHs, greater than 16 Mb, are less frequent, comprising only 1.99% of all detected ROHs. However, the variations in ROH proportion across the classes indicate that medium-length ROHs are more commonly distributed across the genome. This distribution suggests that while Nili Ravi buffalo show a moderate occurrence of inbreeding events, longer ROHs remain relatively rare, emphasizing the genetic diversity within the population.

The average percentage of chromosome coverage by Runs of Homozygosity (ROH) is shown in Figure 1. The highest average coverage was observed on Chromosome 5 (approximately 10%), followed by Chromosome 3 (9.9%) and Chromosome 19 (8.8%). In contrast, Chromosome 11 represented the lowest percentage of ROH coverage, with only 0.9%. These results highlight the variation in ROH distribution across chromosomes, which may reflect differences in chromosomal size, genetic composition, or selective pressures affecting certain regions.

Genomic Regions with Highest ROH Frequency

The Runs of Homozygosity (ROH) analysis across all autosomes of Nili Ravi buffalo identified genomic regions with a high frequency of ROH, reflecting areas potentially influenced by selection pressures or associated with specific traits. Regions on chromosomes 1, 3, 4, 5, 11, 12, 16, 19, and 21 exhibited the highest incidence of ROH, with some spanning over 10 Mb. These extended ROH segments suggest historical inbreeding or selective breeding practices.

Eight prominent regions with higher SNP density were observed on chromosomes 5, 11, and 19, indicating genetic richness and potential functional significance. These regions may harbor genes involved in critical traits such as adaptability, productivity, or disease resistance. The findings highlight key genomic areas that could be targeted in breeding programs to enhance desirable traits or preserve genetic diversity in Nili Ravi buffalo.

Figure 2. illustrates the distribution of SNPs (Single Nucleotide Polymorphisms) within Runs of Homozygosity (ROH) segments, providing insights into the genetic structure of the population. Most ROH segments contain fewer than 100 SNPs, as indicated by the high frequency of bars on the left-hand side of the graph. This predominance of low-SNP ROHs suggests that shorter ROHs or regions with sparse SNP density are more common, reflecting historical recombination events or shared ancestry.

As the number of SNPs per ROH increases, the frequency of segments sharply declines, indicating that longer or denser ROHs are rare within this population. However, a few outliers with SNP counts exceeding 400 are observed, which may represent extended homozygous regions caused by recent inbreeding or selective pressures. These regions could be indicative of conserved genomic segments associated with desirable traits or recent bottleneck events.

Figure 3 show the analysis of SNP distribution within Runs of Homozygosity (ROH) segments in the Nili Ravi buffalo genome revealed a concentration of ROH segments with a low number of SNPs. The majority of ROH segments contained fewer than 100 SNPs, as indicated by the high frequency in the leftmost bins of the histogram. A smaller proportion of ROH segments contained between 100 and 300 SNPs, with a few outliers extending up to 700 SNPs per ROH.

This distribution pattern suggests that most homozygous regions in the genome are relatively short or sparsely populated with SNPs, while longer or SNP-rich ROH segments are rare. These findings provide insights into the genetic structure of the Nili Ravi buffalo, including the levels of inbreeding and the genomic landscape of homozygosity, which are critical for understanding genetic diversity and designing effective breeding strategies.

Table 2 and Table 3 The identified ROH (Runs of Homozygosity) regions on various chromosomes are linked to several milk production candidate genes, which play critical roles in bovine lactation and overall milk yield traits. On chromosome 1, a significant region spanning 35.62 Mb contains COL8A1, GAP43, and MAP3K7CL, genes associated with structural and metabolic functions relevant to milk synthesis. Chromosome 3 harbors genes such as AADAT, ACACA, and SREBF1 in a region of 33.13 Mb, which are vital for fatty acid metabolism and milk fat production. Similarly, chromosome 4 includes ISX and SOX5 in a region of 6.97 Mb, both implicated in regulatory pathways affecting milk yield. On chromosome 5, two critical regions spanning 8.48 Mb and 8.82 Mb contain AKT3, B3GALT6, PRDX6, and CACYBP, genes essential for cellular signaling, oxidative stress response, and protein synthesis. Chromosome 11 features a 26.69 Mb region that includes MAP2K5 and ABCD4, involved in cellular transport and metabolic pathways.

Table 2. Genomic Regions with SNP Density and Gene Count in Nili Ravi Buffalo Genome.

Chromosome	Start Position	End Position	Length (Mb)	Number of SNPs	Number of Genes
1	15019865	152683005	35.62	809	404
1	65180343	73416245	8.24	183	91
3	25972800	99186056	33.13	796	398
3	40855142	49734555	8.88	190	95
3	69721967	86828369	17.11	394	197
4	28698913	70666413	6.97	156	78
4	10983916	21919119	10.94	232	116
5	44849610	84921927	8.48	193	96
5	24221296	33036520	8.82	181	90
5	26297959	59116943	43.34	931	465
11	13266082	93884820	26.69	628	314
12	214630	79549295	23.81	528	264
12	45100371	54904361	9.80	207	103
16	10633369	53866062	21.32	474	237
16	49469455	55213600	5.74	129	64
16	38486559	49280329	10.79	231	115
19	411773	54239343	32.15	724	362
19	11116210	40154446	13.66	322	161
21	37582611	52347729	5.86	131	65

Table 3. Chromosome, Position and annotated candidate genes of ROH region in Nili Ravi buffalo.

Chromosome	Start Position	End Position	Length (Mb)	Candidate Genes
1	15019865	152683005	35.62	COL8A1, GAP43, LOC102402668, BDH1, LOC112587334, LOC102397246, URB1, MRAP, LOC102391318, MAP3K7CL
3	25972800	99186056	33.13	AADAT, ACACA, SREBF1
4	28698913	70666413	6.97	ISX, SOX5, TRNAC-GCA
5	44849610	84921927	8.48	LOC102413112, AKT3
5	24221296	33036520	8.82	B3GALT6, PRDX6, CACYBP
11	13266082	93884820	26.69	DPH6, LOC112587918, LOC112577867, LOC112577895, SKOR1, MAP2K5, LOC112587919, LIN52, ABCD4, VSX2, LOC102409375,
12	214630	79549295	23.81	TRNAS-GGA, TRNAC-GCA
16	10633369	53866062	21.32	LOC112579629, TRNAP-UGG, WNT11, PARVA
19	411773	54239343	32.15	LOC112580680, CDH12, LOC102411033, FBXW11
19	11116210	40154446	13.66	ADAMTS12
21	37582611	52347729	5.86	WDR82, ERC2, ITIH3, ITIH1, LOC112581178, LOC112581179, NEK4, SPCS1, LOC112581322, GLT8D1, GNL3, LOC112581321, LOC112581323, LOC112581320, PBRM1, SMIM4

Table 4. Traits associated with candidate genes found on ROH in Nili Ravi.

Trait	Chr	Candidate Genes	Ref.
Milk Production^*	1	COL8A1	[12]
	1	URB1, LOC112587334, LOC102397246, URB1, MRAP, LOC102391318, MAP3K7CL	[13]
	3	AADAT, ACACA, SREBF1	[14,15]
	4	TRNAC-GCA	[13]
	5	LOC102413112, AKT3, PRDX6	[13,16,17]
	11	LOC112587918, LOC112577867, LOC112577895, SKOR1, MAP2K5, LOC112587919, LIN52, ABCD4, VSX2, LOC102409375	[13]
	12	TRNAS-GGA, TRNAC-GCA	[13]
	16	LOC112579629, WNT11, TRNAP-UGG, PARVA	[13]
	19	LOC112580680, CDH12, LOC102411033, FBXW11, ADAMTS12	[13]
	21	ERC2, ITIH3, ITIH1, LOC112581178, LOC112581179, NEK4, SPCS1, LOC112581322, GLT8D1, GNL3, LOC112581321, LOC112581323, LOC112581320, PBRM1, SMIM4	[13]
Feed Efficiency	3	AADAT	[19]
Fat Deposition	3	ACACA	[20]
Sex Determination	4	SOX5	[21]
Immunity	5	AKT3	[16]
Fertility	11	DPH6	[22]
^*Milk Production includes all traits (Milk Yield (MY), Protein%, Fat%, Milk Synthesis etc.,)

Furthermore, chromosome 12 contains TRNAS-GGA and TRNAC-GCA in a 23.81 Mb region, which may be linked to translational efficiency in protein production for milk. Chromosome 16 harbors WNT11 in a 21.32 Mb region, a gene linked to developmental processes influencing lactation. On chromosome 19, two regions spanning 32.15 Mb and 13.66 Mb include CDH12 and ADAMTS12, respectively, which are associated with tissue remodeling and lactation sustainability. Finally, chromosome 21 features a 5.86 Mb region containing genes such as ITIH1 and GLT8D1, contributing to milk composition and yield. These findings emphasize the genetic basis of milk production traits and highlight critical loci for genetic improvement in dairy animals.

Fixation Index (FST)

The pairwise Wright’s FST values were analyzed to detect signals of diversifying selection between Nili Ravi and Kundi buffalo populations. The mean FST value across the selected loci was 0.682, indicating a moderate to high level of genetic differentiation between the two populations. A total of 24 significant genomic regions were identified as being under divergent selection, based on FST values above 0.65. The highest-ranked SNP, with a Wright’s FST value of 0.740787, was located on Chromosome 11 at position 45,506,934 and is associated with the ARPP19 gene, a candidate involved in regulatory processes.

The SNPs above the 0.65 threshold were distributed across multiple chromosomes, with notable regions on Chromosomes 1, 3, 6, 11, and 19. Candidate genes identified in these regions include CSMD1, ROR1, ARPP19, FAM214A, KCNK12, DNAJC15, SVIL, CADM1, FBXW7, GPT2, and CDH18. Among these, genes such as CSMD1, located on Chromosome 1 (FST = 0.69978), were linked to developmental and immune pathways. SNPs near ROR1 on Chromosome 6 (FST = 0.716468) were associated with multiple QTLs related to milk production, fertility, and growth, underlining their importance in buffalo productivity. Chromosome 11 exhibited a cluster of highly differentiated SNPs (FST = 0.740787) near ARPP19 and FAM214A, which may influence regulatory and developmental functions.

On Chromosome 12, SNPs near KCNK12 (FST = 0.65841) suggested potential adaptation related to metabolic or stress tolerance traits, while Chromosome 16 highlighted SNPs near CADM1 and MAML2, with MAML2 being linked to a growth-related QTL. Other significant SNPs included those near FBXW7 on Chromosome 17, which may regulate metabolic efficiency and reproduction, and DNAJC15 on Chromosome 13, associated with stress responses and energy metabolism.

Additionally, 12 key genes (CSMD1, ROR1, ARPP19, FAM214A, KCNK12, DNAJC15, SVIL, CADM1, FBXW7, GPT2, CDH18, and others) were identified as being associated with traits such as disease resistance, milk production, and body conformation. These findings emphasize the potential of these loci to improve buffalo breeding strategies.

Figure 4, Figure 5 and Figure 6 highlighting the genome-wide distribution of pairwise FST values provide a detailed view of regions under divergent selection. On Chromosome 11, a highly differentiated cluster of SNPs near ARPP19 suggests strong selection in this region. Similarly, Chromosome 6 showed two regions of high differentiation, including one near ROR1, linked to multiple QTLs for milk production and fertility. These findings collectively reveal genomic regions and candidate genes under selection in Nili Ravi and Kundi buffalo populations, offering valuable information for genomic selection and breeding programs. Targeting regions with high differentiation can help improve traits such as milk production, reproduction, and adaptability, while maintaining genetic diversity to support sustainable genetic improvement and conservation efforts in these populations.

Table 5. Genomic regions under selection identified by high FST values with associated candidate genes and QTLs.

CHR	SNP		POS	FST		Candidate Genes		QTLs
1	Affx-79545940_AX-85063870		28538253	0.668151		--		--
1	Affx-38149651_AX-85130207		41198544	0.69978		CSMD1		--
3	Affx-79606245_AX-85124552		71378491	0.668151		LOC112583840		--
3	Affx-79537395_AX-85055280		75328696	0.681419		--		--
4	Affx-79588822_AX-85107011		24307759	0.669166		--		--
6	Affx-79591522_AX-85109728		29202200	0.659603		--		--
6	Affx-79564503_AX-85082549		81014300	0.716468		ROR1		QTL:39967, QTL:39991; QTL:39968; QTL:39977; QTL:39986, QTL:39985; QTL:39974; QTL:39990, QTL:161879, QTL:39973; QTL:39989; QTL:39984, QTL:39971, QTL:39988; QTL:39983; QTL:39969, QTL:39992; QTL:39978; QTL:39987, QTL:39980, QTL:39970; QTL:39979, QTL:39981, QTL:39972, QTL:39982, QTL:39975
11	Affx-79586621_AX-85104799		45506934	0.740787		ARPP19		--
11	Affx-79590053_AX-85108255		45529532	0.740787		FAM214A		--
11	Affx-79562201_AX-85080232		45553437	0.740787		FAM214A		--
12	Affx-79556226_AX-85074216		29453692	0.65841		KCNK12		--
12	Affx-79556227_AX-85074217		29479420	0.65841		KCNK12		--
12	Affx-79540743_AX-85058645		55901962	0.669166		--		--
13	Affx-79535435_AX-85053315		76250277	0.681419		DNAJC15		--
14	Affx-79533482_AX-85051353		49246825	0.681419		SVIL		--
16	Affx-79588482_AX-85106669		59661175	0.683124		CADM1		--
16	Affx-79563844_AX-85081887		71460237	0.669166		MAML2		QTL 181820
17	Affx-79529129_AX-85046979		67044810	0.683124		FBXW7		--
CHR		SNP			POS		FST		Candidate Genes	QTLs
18		Affx-79562665_AX-85181992			14966032		0.664361		GPT2	QTL 25267, QTL 34630, QTL 34630, QTL 25268, QTL 25269
18		Affx-79585107_AX-85103273			53011744		0.693		--	--
19		Affx-79546221_AX-85064152			52241748		0.695334		--	--
19		Affx-79551876_AX-85069842			52458736		0.654996		--	--
19		Affx-79585129_AX-85103295			53279700		0.713347		CDH18	--
24		Affx-79550460_AX-85068420			61183215		0.721758		--	--

Discussion

SNP Density Analysis

The SNP density analysis revealed significant variability in the genomic distribution of SNPs across the Nili Ravi buffalo genome, highlighting regions of both high and low genetic diversity. Chromosome 1 exhibited the highest SNP density, suggesting areas under strong selective pressures or historical mutations. Similar findings have been observed in dairy cattle, where high-density SNP regions are often linked to quantitative trait loci (QTLs) associated with economically important traits such as milk production and fertility [23]. In contrast, chromosomes such as Chr6 and Chr19, characterized by sparse SNP coverage, underscore the necessity for enhanced marker development or sequencing to ensure more uniform genome representation.

The presence of high-density SNP regions provides a valuable resource for genome-wide association studies (GWAS) and detection of QTLs associated with production and adaptability traits [24]. However, low-density regions may impede the detection of such associations, suggesting the need for complementary genomic tools like whole-genome sequencing [12].

Distribution and Implications of Runs of Homozygosity (ROH)

The ROH analysis highlighted a predominance of medium-length segments (2–4 Mb), accounting for 68.9% of all detected ROHs, suggesting moderate historical inbreeding. This aligns with studies in other livestock species where medium-length ROHs often indicate balanced selection pressures maintaining genetic diversity [25]. The scarcity of long ROHs (>16 Mb), comprising only 1.99%, points to effective management practices in minimizing excessive inbreeding.

Chromosomes 5 and 3 showed the highest ROH coverage, indicating regions under possible selection pressures. Genes such as AKT3 and CACYBP, found in these regions, have been linked to cellular signaling and oxidative stress responses critical for dairy productivity and health [16]. In cattle, similar ROH patterns have been associated with loci regulating milk yield and reproductive traits, emphasizing the genetic potential of these regions [26].

Genomic Regions with High ROH Frequency

Eight genomic regions with high ROH frequency were identified, particularly on Chromosomes 1, 5, and 19, emphasizing their significance in selection and adaptability. These regions harbor candidate genes such as AADAT, ACACA, and SREBF1, which are associated with fatty acid metabolism and milk synthesis. Similar findings have been reported in dairy cattle, where genes involved in fat deposition and milk production were localized within high-ROH regions [20].

The presence of key genes in high-ROH regions highlights their potential role in breeding programs to improve milk yield, fat content, and disease resistance in Nili Ravi buffalo. For example, COL8A1 and WNT11, identified on Chromosomes 1 and 16, respectively, are integral to tissue remodeling and developmental processes influencing lactation efficiency [13].

Candidate Genes and Functional Insights

The functional annotation of genes within ROH regions provided insights into their association with economically important traits. For instance, COL8A1 is implicated in structural integrity, essential for maintaining tissue function during lactation [12]. Similarly, AKT3 and PRDX6 are involved in oxidative stress responses and protein synthesis, critical for sustaining milk production under varying environmental conditions [17,27].

Genes associated with fat metabolism, such as ACACA and AADAT, play pivotal roles in enhancing milk fat content, a trait highly desirable in dairy animals. These findings align with studies in Holstein cattle, where fatty acid metabolism genes have been linked to improved milk quality and production [23,24].

Fixation Index (FST) and Divergent Selection

The pairwise FST analysis between Nili Ravi and Kundi buffalo populations revealed moderate genetic differentiation (mean FST = 0.682). Significant loci under divergent selection were identified on Chromosomes 1, 6, and 11, associated with genes such as CSMD1 and ROR1, which are implicated in immune response and reproductive traits [28,29]. These findings suggest population-specific adaptations driven by environmental pressures or selective breeding.

Chromosome 11 exhibited the highest differentiation, with genes like ARPP19 and FAM214A implicated in regulatory and developmental functions. These loci may play a role in enhancing adaptability and productivity in Nili Ravi buffalo, similar to findings in Korean cattle where highly differentiated SNPs were linked to carcass traits and fat deposition [28].

Regions with high FST values, such as those near KCNK12 on Chromosome 12, may be associated with metabolic or stress tolerance traits. These findings are consistent with research in Holstein cattle, where loci under selection were linked to improved stress resilience and energy efficiency [26].

Implications for Breeding and Conservation

The genomic insights gained from this study provide a roadmap for breeding programs aimed at enhancing Nili Ravi buffalo productivity while maintaining genetic diversity. High-ROH regions identified in this study can be targeted for selecting animals with superior milk yield, fat content, and disease resistance. Moreover, the presence of moderate levels of inbreeding inferred from ROH patterns highlights the importance of managing breeding programs to minimize genetic bottlenecks.

The identification of highly differentiated loci (FST) between Nili Ravi and Kundi populations offers opportunities to exploit population-specific adaptations. Targeting these loci through genomic selection can improve traits such as milk production, reproduction, and disease resistance while preserving genetic diversity.

Data Availability Statement

The genotype datasets analyzed in this study were obtained from previously published sources and publicly available repositories. Public data were accessed from Harvard Dataverse (https://doi.org/10.7910/DVN/UGA2QX). Additional datasets were derived from earlier studies and an HEC-NRPU funded project. All processed data supporting the findings of this study are available from the corresponding author upon reasonable request, subject to institutional data-sharing policies.

Author Contributions

H.M. conceived and designed the study. A.A. coordinated data acquisition and assisted with dataset preparation. H.M. performed data curation, quality control, ROH, and FST analyses and led the bioinformatics workflow. W.A.K. provided technical support for genomic analyses and contributed to data interpretation. A.R. and A.M. facilitated sample logistics and contributed domain expertise related to buffalo breeding systems. S.B. contributed to population-level context and breed interpretation, while M.Z.M. supported institutional coordination and data validation. W.A. contributed to sample facilitation and provided logistical support through Maxim Agri (Pvt.) Ltd. H.M. carried out candidate gene annotation, interpreted the results, and drafted the manuscript. All authors critically reviewed the manuscript, provided intellectual input, and approved the final version for submission.

Funding

This work utilized genotype data partially generated under an HEC-NRPU funded project (16844). No additional external funding was received specifically for the present secondary data analysis.

Acknowledgments

The authors acknowledge the Higher Education Commission (HEC) of Pakistan for facilitating access to genotype data. We thank the Livestock and Dairy Development Department (L&DD), Punjab, and other collaborating institutions for their support in sample collection and data generation. The authors also appreciate the technical support provided by colleagues at the Department of Animal Breeding and Genetics, University of Veterinary and Animal Sciences (UVAS), Lahore. Maxim Agri (Pvt.) Ltd. is gratefully acknowledged for providing samples and logistical support.

Conflicts of Interest

The author declares no commercial or financial relationships that could be construed as a potential conflict of interest.

References

Becskei, Z; Mila, S; Dragan, C; Mladen, R; Nikola, P; Marija, P; Sonja, D; Snezana, P. Assessment of water buffalo milk and traditional milk products in a sustainable production system. Sustain. 2020, 12, 12(16):6616. [Google Scholar] [CrossRef]
Naveena, BM; Kiran, M. Buffalo meat quality, composition, and processing characteristics: Contribution to the global economy and nutritional security. Ani. Front. 2014, 4(4), 18–24. [Google Scholar] [CrossRef]
Bilal, MQ; Suleman, M; Raziq, A. Buffalo: black gold of Pakistan. Livestock research for rural development. Pakistan. Pak. Vet. J. 2006, 18(9), 140–151. [Google Scholar]
Falconer, DS; Mackay, TFC. Introduction to Quantitative Genetics, 4th Edition ed; Addison Wesley Longman: Harlow, 1996. [Google Scholar]
Bashir, MK. Genetic and phenotypic aspects of some performance traits of Nili-Ravi buffaloes in Pakistan. Doctoral dissertation, PhD Thesis, Univ. Agri., Faisalabad, Pakistan, 2006. [Google Scholar]
Bashir, MK; Khan, MS; Lateef, M; Mustafa, MI; Khalid, MF; Shahid-ur-Rehman, FU. Environmental factors affecting productive traits and their trends in Nili-Ravi buffaloes. Pak J Life S Sci. 2015, 13(3), 137–44. [Google Scholar]
Purfield, DC; Berry, DP; McParland, S; Bradley, DG. Runs of homozygosity and population history in cattle. BMC Genet. 2012, 13, 70. [Google Scholar] [CrossRef]
Borquis, RR; Baldi, F; Camargo, GMF; Cardoso, DF; Santos, DJA; Lugo, NH; Sargolzaei, M; Schenkel, FS; Albuquerque, LG; Tonhati, H. Water buffalo genome characterization by the Illumina Bovine HD BeadChip. Genet. Mol. Res. 2014, 1(1), 4202–4215. [Google Scholar] [CrossRef]
Peripolli, E; Munari, DP; Silva, MVGB; Lima, ALF; Irgang, R; Baldi, F. Runs of homozygosity: current knowledge and applications in livestock. Anim. Genet. Available from. 2017, 48, 255–271. [Google Scholar] [CrossRef]
Li, J; Liu, J; Liu, S; Campanile, G; Salzano, A; Gasparrini, B; Plastow, G; Zhang, C; Wang, Z; Liang, A; Yang, L. 2020 Genome-wide association study for buffalo mammary gland morphology. J. Dairy Sci. 87(1), 27–31.
Ferenčaković, M; Sölkner, J; Curik, I. Estimating autozygosity from high-throughput information: effects of SNP density and genotyping errors. Genet. Sel. Evo 2013, 45, 1–9. [Google Scholar] [CrossRef]
Lu, XR; Duan, AQ; Li, WQ; Abdel-Shafy, H; Rushdi, HE; Liang, SS; Ma, XY; Liang, XW; Deng, TX. Genome-wide analysis reveals genetic diversity, linkage disequilibrium, and selection for milk production traits in Chinese buffalo breeds. J. Dairy Sci. 2020, 03(5), 4545–56. [Google Scholar] [CrossRef] [PubMed]
Lázaro, S. F.; Tonhati, H.; Oliveira, H. R.; Silva, A. A.; Nascimento, A. V.; Santos, D. J. A.; Stefani, G.; Brito, L. F. Genomic studies of milk-related traits in water buffalo (Bubalus bubalis) based on single-step genomic best linear unbiased prediction and random regression models. Journal of Dairy Science 2021, 104(5), 5768–5793. [Google Scholar] [CrossRef] [PubMed]
Macciotta, NP; Colli, L; Cesarani, A; Ajmone-Marsan, P; Low, WY; Tearle, R; Williams, JL. The distribution of runs of homozygosity in the genome of river and swamp buffaloes reveals a history of adaptation, migration and crossbred events. Genet. Sel. Evol. 2021, 53(1), 1–21. [Google Scholar] [CrossRef] [PubMed]
Stachowiak, M.; Nowacka-Woszuk, J.; Szydlowski, M.; Switonski, M. The ACACA and SREBF1 genes are promising markers for pig carcass and performance traits, but not for fatty acid content in the longissimus dorsi muscle and adipose tissue. Meat Science 2013, 95(1), 64–71. [Google Scholar] [CrossRef] [PubMed]
Farmanullah, F.; Gouda, M.; Min, Z.; Sutong, X.; Kakar, M. U.; Khan, S. U.; Salim, M.; Khan, M.; Rehman, Z. U.; Talpur, H. S.; Khan, F. A.; Pandupuspitasari, N. S.; Shujun, Z. The variation in promoter sequences of the Akt3 gene between cow and buffalo revealed different responses against mastitis. Journal of Genetic Engineering and Biotechnology 2021, 19(1), 164. [Google Scholar] [CrossRef]
Gu, M; Cosenza, G; Iannaccone, M; Macciotta, NPP; Guo, Y; Di Stasio, L; Pauciullo, A. The single nucleotide polymorphism g.133A&C in the stearoyl CoA desaturase gene (SCD) promoter affects gene expression and quali-quantitative properties of river buffalo milk. J. Dairy Sci. 2019, 102, 442–51. [Google Scholar] [CrossRef]
Anwer, I. Characterization of Nili Ravi Buffalo using HD SNP genotyping array. M.Phil thesis, University of Veterianry and Animal Scioiences, Lahore-Pakistan, 2022. [Google Scholar]
Seabury, CM; Oldeschulte, DL; Saatchi, M; Beever, JE; Decker, JE; Halley, YA; Bhattarai, EK; Molaei, M; Freetly, HC; Hansen, SL; Yampara-Iquise, H; Johnson, KA; Kerley, MS; Kim, J; Loy, DD; Marques, E; Neibergs, HL; Schnabel, RD; Shike, DW; Spangler, ML; Weaber, RL; Garrick, DJ; Taylor, JF. Genome-wide association study for feed efficiency and growth traits in U.S. beef cattle. BMC Genomics 2017, 18(1), 386. [Google Scholar] [CrossRef] [PubMed] [PubMed Central]
Du, L.; Li, K.; Chang, T.; An, B.; Liang, M.; Deng, T.; Cao, S.; Du, Y.; Cai, W.; Gao, X.; Xu, L.; Zhang, L.; Li, J.; Gao, H. Integrating genomics and transcriptomics to identify candidate genes for subcutaneous fat deposition in beef cattle. Genomics 2022, 114(4), 110406. [Google Scholar] [CrossRef]
Abdullah, M; Rehman, MS; Rehman, MSN; AlKahtane, AA; Al-Hazani, TM; Hassan, FU; Rehman, SU. Genome-Wide Identification, Evolutionary and Mutational Analysis of the Buffalo Sox Gene Family. Animals (Basel) 2023, 13(14), 2246. [Google Scholar] [CrossRef] [PubMed] [PubMed Central]
Chen, S. Y.; Schenkel, F. S.; Melo, A. L. P.; Oliveira, H. R.; Pedrosa, V. B.; Araujo, A. C.; Melka, M. G.; Brito, L. F. Identifying pleiotropic variants and candidate genes for fertility and reproduction traits in Holstein cattle via association studies based on imputed whole-genome sequence genotypes. BMC Genomics 2022, 23(1), 331. [Google Scholar] [CrossRef]
Pedrosa, V. B.; Schenkel, F. S.; Chen, S. Y.; Oliveira, H. R.; Casey, T. M.; Melka, M. G.; Brito, L. F. Genomewide association analyses of lactation persistency and milk production traits in Holstein cattle based on imputed whole-genome sequence data. Genes 2021, 12(11), 1830. [Google Scholar] [CrossRef]
Matsumoto, H.; Inada, S.; Kobayashi, E.; Abe, T.; Hasebe, H.; Sasazaki, S.; Oyama, K.; Mannen, H. Identification of SNPs in the FASN gene and their effect on fatty acid milk composition in Holstein cattle. Livestock Science 2012, 144(3), 281–284. [Google Scholar] [CrossRef]
Guo, J.; Jorjani, H.; Carlborg, Ö. A genome-wide association study using international breeding-evaluation data identifies major loci affecting production traits and stature in the Brown Swiss cattle breed. BMC Genetics 2012, 13, 82. [Google Scholar] [CrossRef]
Chen, Y.; Steeneveld, W.; Frankena, K.; Leemans, I.; Aardema, H.; Vos, P. L. A. M.; Nielen, M.; Hostens, M. Association between days post-conception and lactation persistency in dairy cattle. Journal of Dairy Science 2024, 107(8), 5794–5804. [Google Scholar] [CrossRef]
Qiu, J.; Ma, Z.; Hong, Z.; Yin, X.; Chen, Y.; Ahmed, H. Q.; Zan, L.; Li, A. Comparative analysis of the whole transcriptome landscapes of muscle and adipose tissue in Qinchuan beef cattle. BMC Genomics 2025, 26, 32. [Google Scholar] [CrossRef]
Lee, S.H.; Van Der Werf, J.H.J; Kim, N.K.; Lee, S.H.; Gondro, C.; Park, E.W.; Oh, S.J.; Gibson, J.P.; Thompson, J.M. QTL and gene expression analyses identify genes affecting carcass weight and marbling on BTA14 in Hanwoo (Korean Cattle). Mamm. Genome 2011, 22, 589–601. [Google Scholar] [CrossRef]
Hoff, J.L.; Decker, J.E.; Schnabel, R.D.; Seabury, C.M.; Neibergs, H.L.; Taylor, J.F. QTL-mapping and genomic prediction for bovine respiratory disease in U.S. Holsteins using sequence imputation and feature selection. BMC Genom. 2019, 20, 1–15. [Google Scholar] [CrossRef]

Figure 1. SNP density distribution across the genome of Nili Ravi buffalo based on the Axiom 90K SNP chip. The heatmap represents the number of SNPs within 1Mb windows across chromosomes, with SNP density color-coded from low (green) to high (red). High-density regions indicate areas of greater genetic variation, while low-density regions highlight sparse SNP coverage.

Figure 2. Average percent coverage by runs of homozygosity in Nili Ravi Buffalo.

Figure 3. Distribution of SNPs per Runs of Homozygosity (ROH) in the Nili Ravi buffalo genome.

Figure 4. Genome-wide distribution of FST values between Nili Ravi and Kundi buffalo populations. The Manhattan plot illustrates FST values, representing genetic differentiation, across all chromosomes. Each dot corresponds to an individual SNP, with its position on the x-axis indicating its chromosomal location and its FST value on the y-axis. The red horizontal line marks the threshold of significant FST values (e.g., 0.6), indicating regions of the genome with high differentiation between the two buffalo populations. These regions may harbor loci under selection or associated with population-specific adaptations.

Figure 5. Genome-wide FST analysis for Chromosome 11 between Nili Ravi and Kundi buffalo populations. The scatter plot illustrates the FST values for SNPs along Chromosome 11, with positions on the x-axis and FST values on the y-axis. The red horizontal line represents the threshold (FST = 0.6) for significant genetic differentiation. A single SNP (green dot) exceeds this threshold, indicating a region of high differentiation that may be under selection or linked to population-specific traits.

Figure 6. Genome-wide FST analysis for Chromosome 6 between Nili Ravi and Kundi buffalo populations. The red horizontal line marks the threshold (FST = 0.65) for significant genetic differentiation. Two SNPs (green dots) exceed this threshold, highlighting regions potentially influenced by selective pressures or associated with unique traits in the populations.

Table 1. Here is the summary table of Runs of Homozygosity (ROH) based on different length categories.

Class	Number of ROH	Percent (%)	ROH length mean (Kb)	Standard Deviation (Kb)
ROH_1-4 Mb	117	77.48	2764.18	517.75
ROH_4-8 Mb	21	13.91	4920.44	849.84
ROH_8-16 Mb	10	6.62	10698.10	1816.42
ROH_16+ Mb	3	1.99	22275.55	9131.46

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.