Genomic Insights into the Origins, Population Structure, and Local Adaptation of Philippine Visayan Native Cattle

TL;DR

This study explores the genetic origins and adaptations of native cattle in the Visayas region of the Philippines, revealing their indicine ancestry and traits suited to tropical island conditions.

Contribution

The study identifies genomic regions linked to local adaptation in Visayan native cattle, offering insights into their ancestry and traits for conservation and breeding.

Findings

Visayan native cattle show indicine ancestry bridging mainland Southeast Asia and China with additional taurine and South Asian indicine components.

Candidate genes suggest adaptations for small stature, heat tolerance, and reproductive performance in tropical island environments.

Low genetic diversity and long ROH segments in Siquijor cattle indicate recent inbreeding.

Abstract

Philippine farmers utilize native cattle for meat, draft power, and income security, but until now their origins and perceived traits have been poorly understood. This study utilized genome-wide SNP to analyze native cattle from the Panay and Siquijor islands, Visayas, comparing them to global breeds to uncover their origin and the genomic regions that potentially allow them to thrive in a tropical island environment. Visayan native cattle mostly belonged to indicine cattle and occupied a middle position between mainland Southeast Asian and Southern Chinese cattle, reflecting past mixing among these regions. Their genetic history is consistent with a complex sequence of ancestry signals that may reflect historical livestock movements associated with Southeast Asian maritime connectivity and later colonial-era introductions. Finally, several candidate regions were identified, and these could be related to the traits perceived by local farmers.

The introduction of domestic cattle to the Philippines is often attributed to Spanish and Chinese sources, yet the origins and adaptive history of Philippine Visayan native cattle remain unclear. This study examined the ancestry, structure, and putative selection signals of the Visayan native cattle from Panay and Siquijor islands (VNC) in a global context. Using genome-wide SNP data, population structure was assessed by PCA, IBS/Nei/FST trees, and ADMIXTURE; historical relationships were explored with migration, f-statistics, and an admixture graph; and positive selection was scanned using commonly used methods such as ROH, Tajima’s D, iHS/XP-EHH, and FST with cross-validation across methods and functional enrichment of the overlapping regions. VNC exhibited low-to-moderate genetic diversity (Ho and He ≈ 0.21; and FIS = 0.01 to 0.02) with Siquijor enriched for long ROH segments indicating recent inbreeding. Across multiple complementary analyses, VNC showed predominantly indicine ancestry and occupied an intermediate bridge-like position between indicine from mainland Southeast Asia and from Southeastern China, with additional components that were most similar to Iberian taurine cattle and South Asian indicine. Moreover, the current study identified putative selection signatures that would possibly provide insights to better understand the local adaptation of VNC under insular tropical conditions of the Philippines: (1) small stature (HOXC cluster, STAC3, NXPH4, STARD13, RTN1), (2) heat tolerance and immune robustness (NDUFA4L2, SHMT2, ATP5MC2, ATF7, R3HDM2, CALCOCO1); (3) early reproductive and maturity reproductive performance (IGF2BP2, KL, LRP1, PDS5B). Overall, the VNC in Panay and Siquijor showed a predominantly indicine ancestry with putatively island-adapted physiology, emphasizing the need for conservation and island-specific breeding that preserves local adaptation while managing inbreeding.

Figures

Full text

1. Introduction

Cattle are an integral yet relatively understudied and underdeveloped component of the Philippine livestock sector. The national cattle inventory reached about 2.5 million head in 2025, with more than 90% of animals raised by smallholder farmers. The Philippines is composed of three major islands, namely Luzon, Visayas, and Mindanao, with the Visayas at the center of the archipelago. The Central and Western Visayas regions are among the top producers of cattle in the country [1]. In these mixed crop–livestock systems, cattle provide meat, limited milk, and draft power [2], and they serve as a financial safety net for rural households during emergencies [3].

Historical records suggest that cattle were first introduced into the Philippines only in the 1500s, when Chinese traders and Spanish colonizers brought domesticated cattle into the archipelago [4]. These early introduced cattle subsequently adapted to local environments and management systems. During the American occupation (1898–1946), exotic commercial breeds (e.g., Gir, Brahman, Sahiwal, Angus, Holstein, Jersey) were imported to upgrade local genetic stocks. From the 1960s onwards, government-led crossbreeding programs between exotic breeds and Philippine native cattle became widespread, ultimately giving rise to the largely non-descript cattle population now found in many parts of the country [2]. Philippine native animals have traditionally been classified and named according to their geographical origin. Based on microsatellite markers, the native cattle in the Philippines were admixtures of Bos taurus taurus, B. t. indicus, and B. javanicus [5]. Using MHC sequencing, Philippine native cattle are clustered together with zebu cattle [6,7] and can be subdivided into two: a northern group comprising the island of Luzon and a southern group comprising Panay, Bohol, and Leyte [7]. Furthermore, analyses of mtDNA d-loop [8] and COI markers [9] indicate that both indicine and taurine maternal lineages can be found in Philippine native cattle.

Among the several documented native cattle populations in the Visayas, the Panay (also known as Iloilo cattle) and Siquijor cattle are the two most prominent populations currently being developed by government programs. These animals are promoted as dual-purpose cattle for both beef and dairy production. National assessments of animal genetic improvement programs in the Philippines highlight fragmented recording systems, low adoption of structured breeding schemes, and a lack of genomic information to guide selection and conservation, all of which hinder effective genetic improvement in native animals [10,11].

The increasing use of SNP arrays over the past decade has enabled large-scale studies of genetic diversity, relationships, and the origin and gene flow of native or indigenous cattle in neighboring Southeast Asian countries and China [12,13,14,15,16]. Understanding origin and gene flow is critical for reconstructing demographic history, quantifying admixture with commercial breeds, and identifying unique genetic resources that warrant conservation. Moreover, even in the absence of performance records, comparisons, using genomic data, among breeds exposed to different selective pressures can be used to identify genomic regions under positive selection [17]. Positive selection refers to the process by which alleles conferring a fitness or productivity advantage increase in frequency more rapidly than expected under neutrality, leaving characteristic footprints such as reduced local diversity, high-frequency derived alleles, or long-range haplotypes in the genome [18,19]. Genome-wide scans can therefore detect loci underlying adaptation to heat, disease, microclimate, and production traits, and knowledge of these candidate regions provides a blueprint for designing island-specific breeding programs. However, SNP-based studies on Philippine native cattle remain extremely limited.

Accordingly, this study addresses the following questions: (1) how Philippine Visayan native cattle are classified and what their origins are; (2) what the current diversity status of the gene pool in the Visayas islands is; and (3) which traits are being targeted by selection in Visayan native cattle that may underpin adaptation to tropical island environments and smallholder production systems. By providing the first integrated genome-wide view of selection signatures in Philippine native cattle, the results are expected to inform science-based breeding and conservation strategies, support ongoing national genetic improvement programs, and contribute evidence for the sustainable use of these locally adapted genetic resources in future cattle development plans.

2. Materials and Methods

2.1. Ethics Approval

For this study, hair follicles were collected from institutional, private, and farmer’s herds. Protocol for the biological sample collection was reviewed and approved by the Institutional Animal Care and Use Committee of Chungbuk National University (Cheongju, Republic of Korea) (CBNUA-2198-23-01).

2.2. Sample Collection and SNP Genotyping

The FAO practical guide on genomic characterization [20] was used in designing the research and determining the number of samples. Hair samples from 61 animals were collected from Siquijor (n = 32) and Panay (n = 29) cattle. Siquijor cattle were collected from the Ubay Stock Farm on the island of Bohol. Panay cattle were collected from the four provinces on the islands of Panay and nearby Guimaras. Animals were randomly selected from several cattle farmers on each island. Villages were identified based on consultation with local agricultural extension officers to target areas within a strong presence of native cattle populations. Precautions were taken to avoid sampling from related individuals. Tail hair samples (containing 20–30 hairs) were collected, ensuring the cleanliness of the extraction area to eliminate potential contamination. Samples were stored on paper cards or resealable plastic containers and kept at room temperature with silica beads until delivery to the genotyping service providers. These providers were responsible for DNA extraction and sample genotyping. Siquijor samples were sent and genotyped using Illumina BovineHD by GeneSeek, Inc., Lincoln, NJ, USA. Panay samples were sent and genotyped using Illumina bovine SNP50 v.3 with the help of TNT Research Co. Ltd., Sejong, Republic of Korea.

2.3. Dataset Preparation and Quality Control

Since the two populations were genotyped using different platforms, only SNPs common to both 50k and HD were utilized. All data were merged, and SNP quality control was performed using the SNP Variation Suite (SVS; Golden Helix, Inc., Bozeman, MT, USA). The SVS Append function was used to merge datasets generated on different SNP arrays. Prior to merging, marker annotations were harmonized to a common reference genome assembly and a consistent marker naming scheme. The physical positions of the SNPs were mapped according to the ARS-UCD1.2 available from https://www.animalgenome.org/repository/cattle/UMC_bovine_coordinates/ (accessed on 1 March 2025). SVS then matched markers across datasets using the harmonized marker ID and physical position coordinates; markers that could not be unambiguously matched (i.e., non-overlapping or inconsistent entries) were excluded so that only the shared SNP set was retained in the final merged dataset.

Moreover, to encompass breeds that might have contributed to the genetic background of the Panay and Siquijor cattle or were related to them, the public WIDDE cattle database (http://widde.toulouse.inra.fr/widde/, accessed on 1 March 2025) and other databases were utilized. Reference breeds were selected based on historical records. These cattle populations were from Southeast Asia [12,13], China [15], Northeast Asia [12,21], South Asia [12,15], Africa [21], Europe [22,23,24,25], and America [12,22].

Different data processes were conducted to further analyze the main dataset. The strand orientations were checked and flipped using the –flip-scan and flip functions of PLINK v1.9 [26], respectively. SNPs assigned to sex chromosomes, lacking genomic locations, and repeat loci were excluded from the analysis. Individuals with more than 10% missing genotypes were excluded [27,28]. To limit ascertainment bias favoring SNPs from European origin, SNPs that were not polymorphic (MAF < 0.01) were discarded [23,29,30]. SNPs with missing genotypes in more than 10% of individuals were also excluded [23,30]. However, SNP filtering based on the Hardy–Weinberg equilibrium was not performed since the studied populations are expected to deviate from Hardy–Weinberg due to their small and possibly substructured population and genetic drift [27,29,31]. Thus, the main dataset was comprised of 740 individuals and 27,980 SNPs spanning the 29 bovine autosomes (Table 1).

2.4. Genetic Diversity Analyses

The genetic diversity indices of the studied populations were assessed using different parameters. Minor allele frequency (MAF), observed homozygosity (Ho), expected homozygosity (He), and inbreeding coefficients (FIS) per population were calculated using PLINK 1.9 [26].

Runs of homozygosity (ROHs) were detected using the consecutiveRUNS function of the detectRUNS v0.9.6 R package [32]. Consecutive runs were selected because they can accurately identify ROHs compared to the sliding window method [33]. The following parameters to define an ROH were adopted from [34,35]: no missing and heterozygote allowed; minimum number of consecutive SNPs was set to 15; minimum ROH length was 1 Mb; and maximum gap between homozygous SNPs was 1 Mb. ROHs shorter than 1 Mb were excluded to avoid short LD-driven islands [36]. ROHs were grouped into five length classes as previously reported by [35]: 1–2, 2–4, 4–8, 8–16, and >16 Mb. An average number of ROHs (nROH) was computed per class. Moreover, inbreeding coefficients based on ROHs (FROH) were calculated genome-wide for each class as the sum of ROH lengths within the class divided by the total autosomal genome length covered by SNPs (2.49 Gb). As described by [37], ROH length classes can be interpreted as rough indicators of the number of generations to the most recent common ancestor, with longer size indicating more recent shared ancestry: 1–2 Mb ≈ 50 generations ago; 2–4 Mb ≈ 25 generations ago; 4–8 Mb ≈ 12 generations ago; 8–16 Mb ≈ 6 generations ago; >16 Mb ≈ 3 generations ago.

Lastly, linkage disequilibrium (LD) and recent effective population size (Ne) were estimated using PLINK v1.9 and SNeP v1.11 [38] and GONE v1.0 [39,40] software, respectively. Pairwise LD (r^2^) was computed in PLINK using the commands: --r2 command to get the LD value of SNP pairs; and the --ld- window-r2 0 to get reports for all pairs. The Ne curve in SNeP typically begins at 13 generations ago. Past demography was also reconstructed with GONE [39] using its default configuration and following the methods described in [40]. This software can start Ne estimation for the very recent past, 1 generation ago.

2.5. Genomic Relationship and Population Structure

Prior to genomic relationship and population structure analysis, linkage disequilibrium (LD) pruning was performed using GoldenHelix SVS with the following parameters: window size: 50; window increment: 5; and LD threshold: 0.2. The LD pruned data resulted in a dataset of 18,036 SNPs.

To detect genetically distinct clusters, principal component analysis (PCA) was conducted using PLINK and visualized with the R package ggplot2 v4.0.0 [41]. Genome-wide identity-by-state (IBS) pairwise distances were estimated to cluster the samples using PLINK. Additionally, F_ST_ and Nei genetic distances among the populations were calculated using the R packages SNPRelate v1.42.0 [42] and adegenet v2.1.11 [43], respectively. The neighbor-joining (NJ) tree was constructed using R package ape 5.8-1 [44] and visualized using the NeighborNet graph of R package phangorn v2.12.1 [45].

Moreover, the ancestry proportion of K ancestral populations was estimated using an unsupervised model-based clustering approach implemented in ADMIXTURE 1.30 [46]. Ten cross-validations for K values ranging from 2 to 30 were conducted. Individual ancestry patterns were generated using CLUMPP [47] and were visualized using the R package pophelper v2.3.1 [48].

LD-pruned data were trimmed to remove populations containing fewer than 5 animals. These data were used to estimate the f3-statistics and the number of migrations using the threepop and treemix functions of TreeMix v 1.13 software [49]. Investigation of historical introgression events and past migrations among Asian indicine cattle populations were modeled using the maximum likelihood algorithm implemented in TreeMix software. The LD-pruned dataset was further subdivided with B. javanicus as the outgroup: (1) only Asian indicine cattle; and (2) selected indicine and taurine cattle. Phylogenetic trees with migration events from 0 to 15 were constructed using blocks of 1000 SNPs, allowing global rearrangements. Each analysis was run with 1000 bootstrap replicates to verify the consistency of edges and nodes. In addition, the analysis was repeated 10 times for each number of allowed migration events. The optimal number of migration events and the consensus tree of the selected migration events were built from the R package BITE2 v2.1.2 [50]. Furthermore, the f3-statistics were utilized to test whether Panay or Siquijor cattle show evidence of admixture from two source populations. The significance of the f3-statistics was assessed using blocks of 1000 SNPs, and a Bonferroni adjustment of the p-value was applied [51].

2.6. Demographic Inference

The origins of PAN and SIQ were further explored using the R package poolfstat v3.0.0 [52] using the methods implemented in [23]. Preliminary f3-statistics, migration results, and historical considerations guided the selection of putatively unadmixed lineages for the initial scaffold: BLI (B. javanicus), JER and LID (European taurine), HAN (Northeast Asian taurine), GIR (Indian indicine), and THS (Southeast Asian indicine). A rooted neighbor-joining scaffold was constructed with rooted.njtree.builder. Admixture graphs were inferred with graph.builder (default settings), and the candidate populations PAN, SIQ, WEL, and WAN were placed on the scaffold. To mitigate order dependence, all permutations of insertion order were evaluated; at each step, the minimum-BIC graph was retained. For each retained and final graph, compare.fitted.stats was used to contrast observed f-statistics with values predicted by the fitted graph parameters. Goodness of fit was summarized by: BIC and worst-fitted Z-scores.

The possible admixture events in PAN, SIQ, and WEL were estimated (in generations) using MALDER v1.0.1 [53] with the “jackknife: YES”, “checkmap: NO”, and “mindis: 0.005” options. Genetic distances between pairs of SNPs were derived from physical distances assuming a cM to Mb ratio of 1.0 [54]. MALDER two-reference tests were run with mindis = 0.5 cM to reduce short-range LD; dates were taken from single-exponential fits with significant amplitude (|Z| ≥ 3) and concordant across multiple reference pairs, while mixed-exponential solutions were retained only when the secondary amplitude was significant and stable across pairs.

Estimated generation times from GONE and MALDER were converted to calendar years as Year = 2000 − (g × 6), with g denoting generations and a 6-year generation interval for local cattle breeds [55].

2.7. Scans for Selection Signatures

Investigating selective sweeps in Visayan cattle provides insights into the evolutionary processes underlying local adaptation. Selection signatures were evaluated at two levels: (1) intrapopulation, with Panay and Siquijor combined as a single Visayan cattle data pool; and (2) interpopulation, by contrasting the Visayan group with Temperate reference composed of Holstein (dairy-type) and Hanwoo (beef-type). Commonly used selection scans were selected to detect complementary signatures of selection: within-population (ROH, Tajima’s D, and iHS); and between-population (F_ST_, and XP-EHH). To detect positive selection underlying local adaptation, complementary scans were integrated to capture distinct genomic footprints of selection, ranging from recent haplotype-based signals to longer-term shifts in diversity and between-population differentiation.

Runs of homozygosity (ROHs), a reduced local variability-based method, are contiguous stretches of homozygous genotypes that serve as indicators of recent positive selection and demographic history [19]. ROHs of the pooled Visayan cattle data were estimated using the same parameters as previously described in the detectRUNS R package. Candidate genomic regions for selective sweeps were subsequently identified using the tableRuns function, and the threshold was set to 0.50 as reported by [56].

Tajima’s D, a site frequency spectrum-based method, is a statistic that distinguishes between neutral and non-neutral evolution by contrasting two measures of genetic diversity, with extreme negative values indicating positive selection and extreme positive values suggesting balancing selection [19,57]. Tajima’s D value was calculated using VCFtools [58] with a window size of 500 kb. Candidate regions for positive selective sweeps were then defined as those falling within the lowest 1% of empirical Tajima’s D [59].

The integrated haplotype score (iHS), an LD-based method, is a statistic used to detect evidence of recent positive selection at a locus by comparing the extended haplotype homozygosity (EHH) of the ancestral and derived alleles [19,57]. Haplotype phasing was done using the default parameters of BEAGLE v5.2 [60]. Information on ancestral and derived alleles for the investigated SNPs to compute iHS was obtained from [61]. The iHS was calculated for each SNP using the R package rehh v3.2.2 [62,63]. Moreover, standardized iHS values were also calculated. Potential candidate regions were identified by using the calc_candidate_regions function with a window of 100 kb and a 50 kb overlap [64]. Strongly negative results indicate that the derived allele has swept across the population, whereas strong positive values indicate that the ancestral allele has been selected [57]. Windows at the top 1% of the |iHS| were considered as candidate regions for the selective sweeps [59,65].

The XP-EHH scores, a haplotype-based differentiation, were calculated for the comparison between Visayan native cattle and Temperate cattle (Hanwoo and Holstein) using the rehh R package, taking the Temperate cattle as the reference population. To detect positive selection in Panay and Siquijor, average XP-EHH scores were computed for 100 kb regions with a 50 kb overlap. Windows at the top 1% of the absolute XP-EHH scores were considered as putative candidate regions.

Population differentiation was assessed using F_ST_ to detect loci under divergent selection, which are characterized by unusually high allele frequency differences between populations [19]. The calculation was performed across the genome in 100 kb windows with a 50 kb step size [64] using VCFtools. The theoretical range of F_ST_ is from 0 to 1, representing no differentiation and fixed differences between populations, respectively. As negative values are considered statistical artifacts of the estimator, any such values were set to zero. Windows with F_ST_ values in the top 0.1% of the empirical distribution were considered candidate regions under positive selection [59,64,65].

2.8. Functional Annotation of Candidate Genomic Regions

Candidate genomic regions were identified by integrating five complementary genome-scan methods (ROH, Tajima’s D, iHS, F_ST_, and XP-EHH) with the following thresholds: ROHs present in ≥50% of the population; Tajima’s D in the bottom 1%; |iHS| in the top 1%; XP-EHH in the top 1%; F_ST_ in the top 0.1%). Overlaps among significant windows were computed using the findOverlaps and pintersect functions of the GenomicRanges R package [66]. Thresholds were chosen to match common practice in comparable selection scans that were conducted in indicine cattle, and genomic regions that passed the threshold criteria in three or more independent methods were considered high-confidence candidates. This approach reduces method-specific false positives [67]; however, demographic processes can generate signatures similar to selection. Accordingly, overlapped regions were treated as putative candidate regions rather than definitive evidence of adaptation.

The GALLO R package was used to identify putative genes and QTL within the candidate genomic regions. QTL and genes were retrieved from the Animal Genome database [68] and Ensembl Genes database (https://www.ensembl.org, accessed on 1 March 2025), respectively. The QTL enrichment analysis was also performed for all the QTLs within the candidate genomic regions using the same GALLO R package. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for further gene function analysis using the DAVID platform (https://davidbioinformatics.nih.gov/, accessed on 1 March 2025). Significantly enriched QTL, GO terms, and KEGG pathways were defined by a threshold using the false discovery rate (FDR ≤ 0.05). Additionally, a comprehensive literature search was conducted to further explore the biological significances of the identified genes.

3. Results

3.1. Intra-Population Genetic Diversity

The genetic diversity indices are shown in Table 2. Across the two Visayan native cattle (VNC) populations, minor allele frequencies (MAFs) and observed (Ho) and expected (He) heterozygosities were nearly identical between populations (PAN–MAF: 0.194 ± 0.140; Ho: 0.211 ± 0.180; He: 0.216 ± 0.177; SIQ–MAF: 0.184 ± 0.142; Ho: 0.210 ± 0.180; He: 0.213 ± 0.177), and FIS was close to zero and slightly positive in both groups (PAN: 0.021 ± 0.193; SIQ: 0.010 ± 0.177). Comparable MAFs and Ho and He indicate similar standing genetic diversity in PAN and SIQ. The small, positive FIS values suggest only mild heterozygote deficit overall.

The FROH from short-to-medium segments (1–8 Mb) was broadly comparable between populations, and, as expected, FROH decreased with increasing ROH length class. However, SIQ had higher FROH_8–16_ and FROH_>16_ than PAN. Numerically, ~90% of ROH segments in both populations fell within the 1–4 Mb range; however, SIQ was enriched for long ROH segments relative to PAN. Short ROHs largely reflect older, background autozygosity, whereas long ROHs arise from recent common ancestors. The excess of long ROHs in SIQ points to more recent inbreeding and/or stronger recent drift in SIQ than in PAN, despite similar overall heterozygosity.

Pairwise LD (r^2^) was evaluated up to 500 kb (Supplementary Figure S1). The LD curves for PAN and SIQ were nearly indistinguishable across short and longer distance bins: r^2^ was high at short inter-marker distances and decreased monotonically as physical distance increased through the 0–500 kb range. Thus, both populations showed the expected genome-wide LD decay with distance, with no pronounced between-population differences in the shape or level of the decay. Recent Ne estimates consistently ranked PAN higher than SIQ across methods (Supplementary Figures S2 and S3). SNeP estimated Ne at ~13 generations as 148 (PAN) and 109 (SIQ), whereas GONE estimated Ne at ~1 generation as 742.55 (PAN) and 78.49 (SIQ). Although the methods probe different recent time windows and return different absolute magnitudes, both trajectories showed a marked recent decline in Ne for the two populations and a smaller contemporary in SIQ. This ranking is concordant with the ROH results, where SIQ harbored longer ROH tracts, consistent with stronger recent drift and/or more recent inbreeding relative to PAN.

3.2. Population Structure and Genomic Relationships of Visayan Native Cattle with Other Cattle Populations

Principal component analysis (PCA) of all samples (Figure 1A) revealed a clear separation between indicine and taurine cattle along PC1 (29.80% of the variance), whereas PC2 (8.45%) further distinguished domestic cattle from other Bos species (B. javanicus and B. frontalis). VNC were clustered within the broader indicine space. In the analysis restricted to indicine cattle with B. javanicus and B. frontalis (Figure 1B), PC1 captured the divergence between indicine and the other Bos species, while PC2 arranged indicine populations along a geographic cline from South Asia through Southeast Asia (SEA) to China. VNC were spread in between SEA and Chinese cattle populations, forming a bridge-like position in PCA space suggesting connectivity or gene flow from mainland SEA toward Southeast and Southwest China. Notably, some BNA individuals (Far Southwest Chinese cattle) clustered within VNC, reinforcing this continuity. Within VNC, SIQ spread toward SEA and Far Southwest Chinese cattle, whereas PAN clustered between the Southeast and Southwest Chinese cattle.

The identity-by-state (IBS) distance tree (Figure 1C) revealed strong clustering by population, with most samples from the same breed grouping together. The topology recovered the deep split between indicine and taurine cattle, with the other Bos lineages forming external branches. Within indicine, a clear regional substructure was evident: South Asian and SEA indicine formed a sister assemblage, whereas Chinese indicine resolved into distinct regional lineages. Far Southwestern Chinese indicine grouped closer to the SEA cluster than to other Chinese clades. The VNC nested within the Chinese indicine assemblage, remaining distinct from both the Southeast and Central Chinese clusters. PAN and SIQ were separated by short internal branches, indicating shallow divergence between the two Visayan populations.

Hierarchical clustering of pairwise F_ST_ values (Figure 1D) recovered the expected deep split between indicine and taurine cattle, with other Bos lineages forming an outgroup branch. Within the indicine clade, South Asian and Ethiopian indicine clustered together, whereas a cohesive SEA indicine lay adjacent to, but distinct from, Chinese indicine, and cattle breeds that have significant B. javanicus ancestry [1,3] formed their distinct cluster. The two Far Southwest Chinese cattle showed little differentiation from the SEA cattle population. Moreover, the VNC fell within the SEA indicine cluster, indicating closer affinities to SEA breeds rather than to Chinese indicine. PAN and SIQ were very close to each other on the tree, consistent with a low divergence (pairwise F_ST_ = 0.04).

The neighbor-net tree constructed from Nei genetic distance (Figure 1E) recovered the major split between indicine and taurine cattle, with other Bos species on long external branches. Within indicine, a pronounced regional substructure was evident: a South Asian and Ethiopian indicine block (India–Pakistan–Bangladesh, Ethiopian, and Brahman) connected to a SEA block (Thailand–Indonesia) and to several Chinese groupings (Southeast, Southwest, and Central). The network displayed extensive reticulation among these indicine blocks, indicating non-tree-like relationships (cline/admixture). The topology broadly mirrored the clustering reported in Gao et al. [15]. However, in this current study, Far Southwest Chinese indicine grouped consistently adjacent to SEA indicine rather than to other Chinese indicine. Moreover, VNC occupied an intermediate position along the Southeast Asian–Southeast China axis, and the short internal splits separating PAN from SIQ indicated shallow divergence between the two Visayan populations.

The ADMIXTURE analysis provided insights into the possible genetic ancestry of the VNC. Unsupervised ADMIXTURE (Figure 2) showed a clear split between indicine and taurine ancestries at K = 2. The VNC populations were predominantly indicine but, relative to SEA and Southeastern Chinese (SEC) cattle, carried a higher proportion of taurine ancestry. At K = 4, a component consistent with the B. javanicus component became apparent, and the taurine background resolved into Northeast Asian and European components. B. javanicus and both taurine ancestries were present in the VNC. At K = 5, indicine ancestry further partitioned into South Asian and East Asian indicine (EAI); SIQ exhibited a higher proportion of South Asian indicine, whereas PAN showed more EAI, with both retaining notable taurine admixture.

The cross-validation identified K = 15 as the optimal model (Supplementary Figure S4). At this model, five indicine (zebu) ancestry components were differentiated: African zebu, South Asian zebu, SEA zebu, EAI, and a Philippine zebu component. African Zebu was concentrated in Ethiopian cattle, with Bangladeshi zebu, Brebes, and Madura showing a detectable shared signal. South Asian Zebu extended from South Asia into parts of SEA. DEC (Southwest Chinese cattle), Ethiopian Zebu, and American Creole showed minor South Asian contributions. SEA indicine occurred broadly in Thailand, Indonesia, the Philippines, and many Chinese indicine. EAI predominated in Southeastern, Southwestern, and Central Chinese indicine and was also present in Visayan native cattle. Lastly, the Philippine zebu component was restricted to Philippine samples, with SIQ exhibiting a larger proportion than PAN.

The f3-statistics was used to test admixture in PAN and SIQ with triplets of the form f3 (Target; Source1; Source2), considering indicine, taurine proxies, and B. javanicus as potential sources. Most negative Z-scores were consistent with admixture between a Southeast Asian indicine lineage (THS or ACE) and a taurine-related lineage (Supplementary Table S1). For PAN, combinations pairing THS with nearly all taurine proxies (except WAG and CCC) yielded significantly negative f3, indicating admixture between THS-related indicine and a taurine-related source. For SIQ, the same THS–taurine pairings generally produced negative f3, but none were significant. Overall, the f3 pattern was concordant with ADMIXTURE results at K = 2, which resolved a taurine component within Visayan native cattle.

To focus on historical migration patterns and gene flow, TreeMix was run in two settings: (1) Asian indicine only; and (2) selected indicine cattle with taurine groups. In the indicine-only analysis, the no migration model explained 87.08% of the variance (Supplementary Figure S5) and placed VNC on the SEA–China axis, consistent with the Asian indicine-restricted PCA. Based on the Δm criterion, the nine-migration event was the most optimum, with the amount of explained variance reaching 97.53%. A further increase in the number of migration events improved the model fit only marginally. The tree divided the cattle populations into three groups (Figure 3): the Bali group, South Asian indicine, and EAI. Unlike the zero-migration model, the nine-migration model placed the VNC inside the EAI cluster, in concordance with the results in the hierarchical clustering of IBS values. Specifically, the VNC shared a common ancestor with the SEC. This shared common ancestor received a gene flow from THS-like and Bali cattle, consistent with the f3-statistics and K = 5 results, respectively. Moreover, SIQ received a gene flow from GIR-like cattle.

In the subset including selected indicine populations together with taurine groups (Supplementary Figure S6), the no-migration model already explained 98.58% of the covariance. Although the Δm criterion favored two migration edges, the improvement was minimal, and no supported edge originating from the taurine cluster into the Visayan branch was observed.

3.3. Infering the Demographic History of Visayan Native Cattle

To further visualize the shared history among SEA and SEC with Visayan indicine, an admixture graph was constructed (Figure 4). For the scaffold, the analysis used B. javanicus, taurine, and indicine anchors: JER and LID (European taurine), HAN (Northeast Asian taurine), and GIR (South Asian indicine). THS was retained as the SEA indicine baseline because, in this study and in prior work [13], it preserved an ancestral signal with limited admixture. WEL was initially considered for the scaffold, but migration analyses indicated that WEL shares a common ancestor with VNC that received gene flow from THS and BLI; consequently, WEL was excluded from the scaffold to avoid confounding scaffold relationships.

The inferred best fitting admixture graph (Figure 4) was the consolidation of the previous results. It shows that Visayan indicine and SEC both derived from an admixed population, Filipino-Sino indicine (FSI); this was in accordance with the results of the hierarchical clustering of IBS values. This FSI had an admixed origin with an ancestry (51%) predominantly from an SEA indicine lineage, which was also reflected in the TreeMix migration and f3-statistics results, and the remaining ancestry was derived from an unsampled Bos population.

The ancestral Visayan indicine further received taurine ancestry (11%) from an Iberian taurine population, since it clustered with LID, combining with the FSI. This was in accordance with K = 5, where SIQ and PAN retained a greater European taurine ancestry compared to WEL and WAN. SIQ shared the based admixed parents with SIQ, but it received South Asian indicine ancestry from GIR (20%). This gene flow was also in parallel with the optimum K value and the migration analysis of TreeMix.

The VNC experienced three admixtures: $[eqn]$ (formation of the Filipino-Sino indicine); $[eqn]$ (formation of the Visayan indicine); and $[eqn]$ (formation of the SIQ). For $[eqn]$ , WEL was used as the population target, with THS as one of the source populations and a taurine proxy source as either HAN or JER. The THS and HAN as source proxies obtained the highest amplitude, suggesting that these populations were the closest to the actual source population. The estimated timing for the admixture events leading to FSI was found to be $[eqn]$ = 181 ± 28 generations (calendar year 914 ± 168). For estimating $[eqn]$ , PAN was used as the population target, and WEL and LID were used as source population proxies. The estimated timing for the admixture events leading to Visayan indicine was found to be $[eqn]$ = 60 ± 6 generations (calendar year 1640 ± 36). Lastly, the $[eqn]$ was estimated to be 8 ± 1 generations (calendar year 1952 ± 6) using SIQ as the population target and GIR and PAN as source proxies.

3.4. Candidate Regions Under Positive Selection

Runs-of-homozygosity (ROH) analysis identified multiple ROH islands shared by ≥50% of Visayan native cattle (Supplementary Figure S7A), collectively harboring 1857 genes (Figure 5A; Supplementary Table S2). Other genomic scan methods detected additional candidates (Figure 5A; Supplementary Figure S7B–E): the iHS scan identified 62 genes (Supplementary Table S3); Tajima’s D identified 267 genes (Supplementary Table S4); XP-EHH identified 151 genes (Supplementary Table S5); and windowed F_ST_ identified 56 genes (Supplementary Table S6). To increase robustness and minimize method-specific false positives, signals were cross-validated, retaining loci supported by ≥3 methods. This yielded seven putative genomic regions on chromosomes 1, 5, 10, 12, and 14 comprising 21 genes which included IGF2BP2, HOXC6 and HOXC8-13, CALCOCO1, ATP5MC2, ATF7, R3HDM2, STAC3, NDUFA4L2, SHMT2, NXPH4, LRP1, RTN1, STARD13, KL, and PDS5B (Table 3).

QTL annotations mapped 132 traits associated into these candidate regions. The most frequent category was reproduction, whereas milk and meat/carcass showed similar representation; exterior, production, and health were less frequent categories (Figure 5B). SNP-trait enrichment analysis further identified that the candidate genomic regions harbored SNPs putatively related to inhibin level, hematocrit, coat color, body height, and age at first calving (Figure 5C), which are traits that plausibly support the hypothesis of adaptation of VNC to hot, low-input production environments.

GO enrichment analyses revealed putative processes related to anatomical and physiological adaptation of VNC under tropical island conditions. These are “anterior/posterior pattern specification”, “proximal/distal pattern formation”, “regulation of transcription by RNA polymerase II”, “Nucleus”, “DNA-binding transcription factor-activity RNA polymerase II specific”, and “RNA polymerase II cis-regulatory region sequence-specific DNA binding” (Figure 5D).

4. Discussion

In the Philippines, genetic research on native cattle is very limited. In the past, molecular methods were used to assess genetic diversity and relationships using microsatellites [5,69], mitochondrial markers [8,9], and BoLA-DRB3 sequences [7]. The current study is the first time that the whole genome of Philippine native cattle has been studied using medium-density SNP chips and related to worldwide cattle breeds.

Genetic variation is needed for science-based decisions in breeding and conservation. The genetic diversity in VNC was higher than the values previously reported for Thai, Indonesian, and SEC in Gao et al. [15] but lower than those in the re-analysis by Sudrajad et al. [14]. By contrast, Bangladeshi cattle showed comparable diversity to the present study [29]. As expected, diversity remained lower compared to taurine breeds [15]. Between-study discrepancies are expected because SNP-array heterozygosity is sensitive to ascertainment bias and QC filtering parameters which can inflate or deflate Ho or He and shift allele-frequency spectra [16,70,71]. Moreover, the design of an SNP assay can produce lower average MAFs because of ascertainment bias [72,73]. Hence, the observed diversity likely underestimates the full extent of genome-wide variation that would be captured by sequence-based methods or less biased panels.

Most ROH segments in both populations fell in the short-to-medium classes (1–4 Mb), which typically reflect older background autozygosity attributable to distant common ancestors (up to ~50 generations) [20]. In contrast, SIQ showed an enrichment of long ROH segments (8–16 Mb and >16 Mb) compared to PAN, which are characteristics of recent inbreeding/strong drift (~5–8 and ≤ 3 generations). This ROH profile explains how overall heterozygosity can remain moderate while recent autozygosity increases and is consistent with smaller recent Ne.

The declining Ne trajectories observed mirror patterns reported for SEA [14], Bangladeshi [29], and Indian cattle [74]. From a management perspective, effective sizes of 50–100 are often recommended to limit inbreeding and preserve adaptive potential [75,76,77]. The most recent Ne of SIQ (~78 in GONE) sits near the lower end of this range and combined with its excess of long ROHs, this underscores the need to prevent genetic erosion. Although PAN shows higher recent Ne estimates, its downward trend is concerning.

The observed low-to-moderate diversity in the current study can be attributed to the combinations of adaptation to local environment, geographical isolation, and historical population bottleneck [16]. Overall, these findings emphasize the urgent need for structured, island-specific breeding programs that preserve local adaptation and productivity while reducing long ROH segments and lowering extinction risk. Practical actions include widening and rotating sire exchange within islands, minimizing co-ancestry in mate allocation, and instituting a national herdbook database and routine genomic monitoring to track and manage inbreeding over time [78,79,80].

The genomic relationship between the two VNC and other cattle populations was evaluated using PCA, IBS, F_ST_, and Nei’s genetic distance. All methods were supported by the admixture and migration results. Inferred demographic history was reconstructed using admixture graph and time.

Visayan native cattle, PAN and SIQ, were predominantly indicine and consistently occupied an intermediate position between SEA and SEC indicine across PCA and IBS clustering and Nei’s and F_ST_ trees. This pattern, together with supported migration edges and f3-statistics, inferred a history of gene flow among the cattle population in this region. The results also highlight the role of human-mediated (e.g., colonization and upgrading programs) activities in shaping animal genomes. Notably, this study hypothesized the connections among EAI (e.g., WAN and WEL) and Iberian cattle with the Visayan native cattle.

Indicine cattle are thought to have dispersed from the Indian subcontinent to East Asia through the inland routes [81], but the accumulating genomic and geographic evidence supports a different story. The rugged ranges of the Rakhine/Arakan Yoma mountains between Bangladesh and Myanmar can act as a strong barrier that can hinder the migration of indicine cattle [82]. Chen et al. [83] proposed that the migration of indicine cattle eastward was facilitated by a coastal route. Moreover, Sudrajad et al. [14] hypothesized that the migration of indicines to Indonesia might be facilitated by sea routes. Archaeological and historical records documented a sustained maritime exchange linking South Asian kingdoms and Southeast Asian thalassocracies from the late first millennium BCE onward, providing plausible corridors for cattle migration alongside human movement [84,85,86]. Recent whole-genome studies of indicine have postulated a route from the Bay of Bengal coastal corridor though Myanmar littoral into Thailand/Cambodia and onward into northern Vietnam and south China. This scenario is further supported by clear footprints of introgression from regional wild bovines in mainland SEA [83,87].

In this study, Far Southwestern Chinese cattle (DEH and BNA) clustered more closely with South and Southeast Asian cattle than with other Chinese indicine. This pattern suggests that indicine reached East Asia via at least two maritime corridors: (1) a coastal route skirting mainland Southeast Asia and (2) a sea route into the islands of Southeast Asia. This study explored the latter hypothesis.

Thai cattle share ancestry with Indonesian cattle in Sumatra (e.g., ACE and PES), which also carries B. javanicus ancestry [13,14]. Southern Thailand and Sumatra were once part of the Sri Vijaya empire [88,89]. Consequently, the estimated admixture time of the FSI, 181 ± 28 generations (calendar year 914 ± 168), falls within the timeframe of the recorded period of interaction between the Sri Vijaya empire and the Song Dynasty of China with early Philippine kingdoms [90]. Before the arrival of Magellan, the Philippines was considered as one of the centers of trade among China, India, the Spice Islands, and Sri Vijaya [90]. This network could have facilitated the movement of indicine into China. Moreover, the Philippines is inferred not to have been the sole source of indicine in China but rather a connected contact zone within a wider maritime network.

Studies conducted on Chinese indicine inferred BLI-like and unknown or unsampled Bos-like ghost species introgression into EAI [83,87]. Visayan zebu ancestry shares features of other Bos species introgression with EAI (e.g., WAN and WEL) based on migration analysis and admixture graph. The geographic location of admixture (FSI) to between THS-like and this unsampled Bos remains uncertain; however, this study hypothesizes that it might have occurred in the Philippines. Unlike neighboring countries that harbor Bos javanicus, Bos frontalis, and Bos sauveli [91], archeological and historical records of wild Bos species in the Philippines are limited. Even so, parts of the Philippines were connected to mainland Asia via Sundaland during the last glacial maximum [92], thus facilitating migration of various species from Borneo to the Philippines [93]. The Philippines also harbors two endemic Bubalus species, a genus closely related to Bos: (1) the extant Bubalus mindorensis, whose ancestors are believed to have reached Mindoro via land bridges from Sundaland [94]; and (2) B. cebuensis, an extinct species from the Pleistocene or Holocene discovered on Cebu Island within the Negros-Panay Philippine Faunal Region [95]. The earliest archaeological evidence of cattle in the Philippines comes from Metal Age deposits at the Nagsabaran site in northern Luzon, dated to ~500 BCE [96]. Whether this remain represents B. taurus, B. javanicus, or B. gaurus requires further investigation [97]. Accordingly, it is recommended to do whole-genome sequencing using a larger sample size that includes diverse Philippine cattle populations, alongside genetic analyses of the Nagsabaran remain, to clarify the presence and role of wild Bos lineages in the archipelago.

The inferred timing of the Visayan indicine formation (~60 ± 6 generations; ca. 1640 CE) is temporally compatible within the Spanish colonial period (1565–1898), when the Manila–Acapulco galleon network (1573–1815) tightly linked the Philippines and New Spain (Mexico) and large-scale ranching took place in the archipelago [96]. Documentary work on colonial land tenure and stock-raising shows that cattle herds and estancias proliferated under Spanish rule, consistent with Iberian livestock inputs during this era [2,96]. The much more recent GIR-like contribution into SIQ (~8 ± 1 generations; ca. 1950s) is chronologically consistent with post-World War II Philippine government programs that expanded exotic breed dissemination and farmer-level herd build-up through livestock dispersal initiatives (e.g., “Operation Livestock Dispersal,” R.A. 1499 of 1956). SIQ exhibited an excess of very long ROH segments (>16 Mb), a hallmark of recent inbreeding/consanguinity and a small contemporary effective population size (Ne). This pattern, together with the lower recent Ne observed in SIQ, is fully consistent with recent mating among close relatives within a reduced breeding pool, likely exacerbated by the absence of systematic pedigree recording.

Allele sharing between domestic cattle and other bovines can arise from true introgression but also from incomplete lineage sorting/shared ancestral polymorphism, particularly among recently diverged Bos lineages. In addition, inferences from SNP-array data may be influenced by ascertainment bias and by the composition of the reference panel, and model-based methods (ADMIXTURE, TreeMix, admixture graphs) can attribute residual covariance to unsampled or ghost sources when relevant donors are unsampled or the demographic history is more complex than the fitted model [98,99,100]. Accordingly, the BLI-like and unsampled lineage signals were treated as suggestive hypotheses rather than conclusive evidence of multispecies ancestry, pending validation using broader wild bovid sample and whole-genome sequence data, ideally including historically dated or archaeological specimens from the Philippines and other regions. Consequently, links between the inferred ancestry signals and specific historical episodes (maritime trade connectivity, Spanish-era introductions, and modern South Asian indicine introgression) should be viewed as plausible hypotheses rather than definitive reconstructions [49,101,102,103]. LD-based admixture dates are approximate and depend on assumptions (e.g., generation interval, recombination map, and the suitability of reference proxies), so the converted calendar year estimates should be interpreted as broad time windows rather than exact historical dates [104].

The integrated selection scans highlight three putative biological themes among candidate regions in the sampled VNC populations (PAN and SIQ): (1) potential morphological adaptation consistent with insular tropical conditions; (2) putative regulatory and metabolic processes plausibly relevant to heat tolerance and disease resistance; and (3) candidate loci potentially associated with age at sexual maturity and reproductive performance under subsistence farming. Importantly, these themes should be interpreted cautiously because island isolation, bottlenecks, and admixture can generate sweep-like patterns and can enrich functional categories even under neutrality [105,106,107].

The Philippine native cattle are small-framed at maturity (~105 cm withers height at 3 years) [108]. In this study, a putative genomic signal consistent with this phenotype was the enrichment of anterior–posterior pattern specification (GO:0009952) and proximal–distal patterning (GO:0009954), driven largely by candidate selective sweep across the HOXC cluster (notably HOXC6 and HOXC8–13). HOX genes are canonical transcription factors that specify axial and limb identities, including the cervico-thoracic boundary and elements of the fore/hind limb, in vertebrates. The co-occurrence of HOXC signals with enrichment for body height QTL is compatible with selection affecting small stature and limb conformation in VNC [109,110,111,112,113]. This pattern is also consistent with the broader phenomenon of insular or island dwarfism, in which large mammals isolated on islands may evolve smaller body size over relatively short evolutionary timescales; this has been documented for feral cattle on Amsterdam Island, which underwent rapid dwarfing under restricted resources and weak predation [114,115]. Recent syntheses further emphasize that insular dwarfism can reflect integrated shifts in growth rate, age at first reproduction, and longevity under resource-restricted environments [116]. Taken together, the HOXC-region signals and body height QTL enrichment are consistent with a candidate hypothesis of selection on smaller body size in insular, resource-limited systems.

Additional QTL enrichment for coat color, together with neuromuscular and neuronal candidate genes such as STAC3 [117,118,119,120], NXPH4 [121], STARD13 [122,123], and RTN1 [124], suggests the putative involvement of pathways that may contribute to morphology and performance in tropical island environments. While these candidates are plausibly relevant to locomotor function and neuromuscular control that could matter for free-ranging, small-bodied cattle in rugged landscapes, this interpretation remains speculative without direct phenotype or functional validation. Coat color loci are well-known selection targets in cattle and shape breed-specific pigmentation that influences heat load, UV exposure, and ectoparasite burden, especially in zebu and tropical dairy breeds [125]. Although well-known pigment genes (MC1R, ASIP, KIT) were not among the detected candidate genes, the coat color QTL enrichment could reflect linkage disequilibrium with nearby pigmentation loci or pleiotropic regions where developmental and metabolic genes influence both color and other traits.

Beyond skeletal patterning, GO enrichment for regulation of transcription by RNA polymerase II (GO:0006357), nucleus (GO:0005634), and RNA Pol II-specific DNA-binding and cis-regulatory activity (GO:0000981, GO:0000978) is consistent with selection acting on regulatory homeostasis that buffers metabolic stress [126,127,128]. Several identified candidate genes align with metabolic homeostasis and health, including HOXC12 and HOXC13 [111,129], ATF7 [130], PDS5B [131,132,133], and R3HDM2 [134]. Furthermore, CALCOCO1 has roles in energy/protein metabolism affecting meat quality and carcass traits in beef cattle [135,136,137], and ATF7 has been linked with milk, fat, and protein yield [138]. IGF2BP2 is an mRNA-binding protein that recognizes m^6^A-modified transcripts and stabilizes key mRNAs involved in growth and metabolism [139] and is repeatedly associated with metabolic disease and energy-balance traits in humans [140,141,142]. Collectively, these observations support a hypothesis that the detected candidate regions include transcriptional and post-transcriptional regulators relevant to growth, metabolism, and health, consistent with polygenic fine-tuning for tropical smallholder environments.

This putative regulatory theme may extend into blood physiology and cellular stress responses. The enrichment for hematocrit-related (HCT) QTL, together with colocalization with candidate genes NDUFA4L2 [143,144], SHMT2 [145,146,147,148], CALCOCO1, LRP1 [149,150], ATP5MC2 [151], and RTN1 [152,153], is compatible with selection on oxygen transport and tissue-level coping with hypoxia, heat, and oxidative stress. However, these inferences remain provisional because hematological traits and stress-response pathways are influenced by complex polygenic architecture and demographic history. HCT-related QTL suggests selection on blood oxygen transport and systemic health. In beef cattle, complete blood count traits (including RBC, HCT, and hemoglobin) are heritable and associated with performance, health, and feed efficiency [154,155]; in yaks and high-altitude small ruminants, elevated hemoglobin/HCT are central to hypoxia adaptation [156,157]. The colocalization of hematocrit QTL with these candidate genes is putatively related not only on RBC indices but also on mitochondrial and endoplasmic reticulum stress pathways that limit oxidative damage and sustain productivity in hot, parasite-rich island environments.

Philippine native animals are often perceived to be early maturing [3]. Enrichment for inhibin level and age at first calving (AFC), along with candidate genes such as IGF2BP2, KL, LRP1, and PDS5B, suggest putative involvement in polygenic adjustment of early reproductive maturity under subsistence production systems. AFC integrates growth, puberty onset, and management; large-scale GWAS in Holstein, Hanwoo, and Nellore show that AFC is highly polygenic, enriched for genes in folliculogenesis, hormone signaling, and energy balance, and sensitive to environmental conditions [158,159,160,161]. Analyses of beef heifers further show genetic correlations between body size/structure and AFC, reflecting life-history trade-offs between growth and early reproduction [162]. IGF2BP2 variants are associated with litter size in goats [163,164], and IMP2/IGF2BP2 has roles in oocyte/embryo development [165,166]. The KL/Klotho axis modulates oxidative stress and developmental competence in early bovine embryos [167]. KL also interfaces with FGF23–vitamin D and IGF/insulin signaling, and KL-deficient mice show gonadal atrophy and reduced fertility, linking KL to reproductive longevity and metabolic resilience [168,169,170]. The PDS5B, via cohesin regulation, is essential for normal gametogenesis and embryogenesis [131,132,133], while LRP1, a multifunctional lipoprotein receptor, links lipid uptake to thermogenesis. Adipocyte-specific LRP1 knockout in mice disrupts lipid uptake and thermogenesis, indicating that LRP1 is important for systemic energy homeostasis and adaptive thermogenic responses [171,172]. In cattle, efficient energy use is closely tied to the resumption of cyclicity and maintenance of pregnancy. LRP1 coordinates lipid handling with systemic energy homeostasis in adipose tissue, a pathway known to influence the resumption of cyclicity and maintenance of pregnancy under nutritional stress [173].

Finally, it is important to emphasize the key limitations of genome-wide selection scans in island populations. Bottlenecks, isolation, population structure, and admixture can generate reductions in diversity, long haplotypes, skewed site-frequency spectra, and high differentiation that resemble selective sweeps. This is particularly relevant for populations with small effective size and recent inbreeding, where ROH enrichment and extreme values for the methods used may reflect demographic history in addition to selection. Although multiple complementary statistics were integrated, and candidate regions were prioritized by overlapping across methods, this approach reduces but does not eliminate false positives driven by demography and recombination-rate heterogeneity. Moreover, functional annotation and QTL overlap can suggest biologically plausible processes but do not establish causal genotype–phenotype relationships: genes are pleiotropic, QTLs originate from heterogeneous breeds and environments, and enrichment analyses are sensitive to annotation completeness and bias. Accordingly, the three themes presented here should be viewed as hypothesis-generating interpretations of candidate genomic regions rather than confirmed adaptive mechanisms. Therefore, the regions reported here are best considered candidate loci potentially relevant to local adaptation in the sampled VNC populations and that require validation through explicit demographic modeling, environmental association, and phenotype- or function-based studies. Finally, because PAN and SIQ may have island-specific demographic histories, some signals may be population-specific rather than representative of Visayan native cattle across all islands or Philippine native cattle more broadly; additional sampling across other Visayan and/or Philippine populations is needed in future studies to assess generality.

5. Conclusions

This study presents the first genome-wide study of Philippine Visayan native cattle in a global context. PAN and SIQ are predominantly indicine and occupy a bridge-like position between mainland Southeast Asian and China with historical taurine input likely of Iberian origin and a very recent GIR-like contribution. Analyses on genetic diversity further indicate a declining Ne and signals of recent inbreeding underscoring the urgent need for intervention. To safeguard these genetic resources, the establishment of a national herdbook linked to routine genomic monitoring is recommended to track diversity, manage inbreeding, and coordinate selection decisions. Selection signature analyses converge on three putative biological themes that align with tropical, low-input island production systems: (1) body size, (2) thermotolerance and disease resistance, and (3) earlier reproductive maturity. Accordingly, selection objectives should emphasize heat tolerance, health, and reproductive efficiency, rather than attempting to mirror the pursuit of maximal profit typical of intensive commercial breeding programs. Overall, the findings support the development of island-specific breeding programs that maintain local adaptation while actively controlling inbreeding, thereby strengthening the long-term productivity and resilience of Philippine native cattle.