All animal procedures were reviewed and approved by the Research Ethics Committee at Faculty of Veterinary Medicine, Universitas Gadjah Mada No.0120/EC-FKH/Eks./2022. Genotypes for the worldwide cattle breeds were obtained from previously published data [14,15], and therefore, no ethical approval was required.

This study utilized genomic data of SO cattle in Sumba island, consisting of 79 animals with medium density SNPs. Animals were collected from three SO cattle populations on East Sumba island, which had similar conditions and climate. The genomic data was obtained from DNA which were extracted from the blood samples. DNA extraction was conducted at the animal genetics laboratory of the Indonesian Research Institute for Animal Production in Bogor, following the procedure outlined in the Wizard® Promega isolation kit method. DNA evaluations, i.e. concentration quantification and purity, were performed using Picogreen and NanoDrop devices, respectively. DNA samples were deemed to be of sufficient quality if the concentration was at least 20 ng/μL and the absorbance ratio at 260 and 280 nm was more than 1.8. Appropriate samples were then genotyped using the Illumina BovineSNP50 v3 Beadchip. These DNA evaluations and genotyping process were finished by the service of Macrogen Inc. in Republic of Korea. We added genotypes of 209 animals consisted of 11 cattle populations which were selected from previously published genome datasets by Dixit et al [14] and Decker et al [15] to support the analyses in our study, i.e. Bali, Ongole Grade, Ongole, Nellore, Hariana, Tharparkar, Brahman, Limousine, Simmental, Hereford, and Shorthorn. The last two populations, in particular, were only used in the admixture analysis to determine whether a true Brahman cross (BX) mixture exists in the SO cattle population. The detail of datasets information were summarized in Table 1.

Quality control procedures were performed using PLINK v1.9 software [16] to ensure the reliability and accuracy of the data. SNP variants would be maintained if the Hardy–Weinberg Equilibrium (HWE) value was not <1×10⁻⁴. Any samples with a call rate of less than 90% were removed, as low call rates can indicate poor DNA quality or technical issues during genotyping. We also excluded any SNPs with a call rate of less than 90% or a minor allele frequency of less than 0.01, as these could indicate genotyping errors or non-polymorphic markers. Only autosomal SNPs with known genomic position based on UMD3.1 assembly and found in all populations were used in the subsequent analyses.

This study began by calculating several diversity parameters, including expected heterozygosity (H_e), observed heterozygosity (H_o), inbreeding coefficient (F_IS), pairwise fixation indices between-population (F_ST), and Nei’s genetic distance to assess the genetic diversity of the SO and the comparators cattle populations. PLINK v1.9 [16] and R v4.1.0 [17] software were used to calculate these parameters. A genetic relationship matrix (GRM) of SO cattle was estimated using GCTA v1.93.2 [18]. A diagonal GRM value referring to the inbreeding of the animal itself and off-diagonal value referring to relationship between animals in the population.

This study also investigated the population structure using two different clustering methods: multi-dimensional scaling (MDS) and admixture analyses, which were performed using PLINK v1.9 [16] and fastSTRUCTURE [19], respectively. To determine the optimal number of subpopulations in the studied populations, the admixture analysis was ran using various estimated ancestral populations (referred to as K values). The cross-validation (CV) analysis was performed to select a minimum error value for K. Both clustering results were plotted in R v4.1.0 [17] to get understandable visualization. TreeMix v1.12 [20] was used to conduct extensive analysis of the relationships between cattle populations. A phylogenetic tree model was generated by positioning Bali cattle as root. The threepop command line was used to perform f3 statistical analyses introduced by Reich et al [21] to identify possible traces of introgression. As possible introgression events of populations B and C into population A, three-population statistical models (A, B, and C) with a significant negative value for both the f3 statistic and the Z-score were chosen.

The cattle population's demographic history was reflected in the estimated current and former effective population sizes (N_e). Sved and Hill [22] explained how past N_e was estimated from the linkage disequilibrium (LD) value bases on the following equation: E(r²) = 1/(1+4N_ec). LD was annotated as r² before N_e to measure the correlation of alleles at two loci and c was a recombination fraction in Morgan's unit. To estimate LD and N_e, the default PLINK v1.9 [16] and SNeP v1.0 [23] approaches were used, respectively. The --r2 command was used to get the LD value of SNP pairs, and the --ld- window-r2 command was used to get reports for all pairs. The N_e analysis was performed by using the standard input parameters in GONE approach [24] to infer a population's recent demographic history from SNP data of a small sample of individuals. Following the estimated times in the horizontal ordinate, the historical N_e values were plotted using R v4.1.0 [17].

The genotyping process using the Illumina BovineSNP50 v3 Beadchip yielded a total of 53,218 SNPs on SO cattle. The observed genotyping call rate was 99.28%. Due to the quality control process, 28% of the variants were removed, leaving only 38,310 SNPs to be merged with the other cleaned genotypes of other cattle breeds (Table 1). Combined file consisted of 22,652 variants which have been encountered on 288 animals from 12 cattle breeds. Variability of the samples in SO population was proven by the negative off-diagonal variances in the GRM analysis, i.e. −0.0131. SO cattle showed commensurate level of both H_e and H_o with the other indicines, where the position is lower than taurines and higher than Bos javanicus (Bali cattle). Ongole Grade cattle had the highest H_e and H_o among indicines. Discovered in all populations, the H_o was slightly higher than the H_e. Likewise with the F_IS of SO cattle, the value (−0.21) was almost the same as other indicines. Moreover, the F_IS of cattle populations observed in this study were all negative.

An initial investigation into genetic diversity is summarized in the pairwise F_ST and Nei's genetic distance matrix (Table 2). This matrix clarified the distinctions and differences between SO and other cattle breeds. SO cattle were close to indicines but far from taurines and Bali cattle. Among indicines, Nellore cattle were the closest to Ongole from India (F_ST 0.0025; Nei 0.0061), while SO cattle were the farthest (F_ST 0.0074; Nei 0.0096). However, SO cattle had the shortest distance to Ongole cattle than to the others.

Best illustration of the cattle genetic cluster could be explained by the MDS plot (Figure 1). A three-dimensional plot was used to clarify SO cattle position among indicines, taurines, and Bali cattle (Bos javanicus) coordinate. Interestingly, SO cattle populations tend to assemble in one colony but still near to indicines cluster, while Ongole, Nellore, Ongole Grade, Hariana, Tharparkar, and Brahman were close to one another. The admixture diagram (Figure 2) also shows the results of a more detailed cluster analysis. Based on the minimum CV error calculation (Supplementary Figure S1), K = 8 was chosen as the sensible model in this analysis. SO cattle could clearly be distinguished from other indicines, but the majority of their mixtures were similar to those found in Ongole cattle from India. Furthermore, the admixture with K = 3 was still displayed because it can be used to describe the sub species clusters, namely Bos javanicus (black color), Bos indicus (orange color), and Bos taurus (gray color).

The phylogenetic tree was developed specially for indicines and Bali cattle positioned as a root (Figure 3). Nellore and Ongole cattle from India were directly linked in the tree and SO cattle were closely related to those branch, while Ongole Grade cattle were in line after SO branch. The other indicines, i.e. Brahman, Tharparkar, and Hariana were in the different direction. Moreover, only one migration edge was observed, originating from the direction of Bali to Ongole Grade cattle. The f3 statistic was generated and the most significant results marked by negative value of both f3 statistics and Z-score were summarized in Table 3. This analysis only detect the possible introgression of Bali cattle (Bos javanicus) and indicines into Ongole Grade cattle. Here, SO cattle were the most likely candidate representing indicine ancestry (f3 statistic −0.0055; Z-score −11.4887). We looked more closely at the f3 statistics in which the Population A was SO cattle. This step was completed in order to identify the cattle breeds that have an impact to the SO population. Despite the fact that the values of both f3 statistic and Z-score were not negative, we select the combination of three-population statistical models with the lowest value. As a result, we identified that Ongole or Nellore (indicines originated from India) and Brahman cattle were the populations that suitable in the columns of Population B and C, respectively.

An estimation of LD values based on r² equation were done for all of SNP pairs and depicted following interallelic distances (Supplementary Figure S2). The observed LD values were high over a short allele distance and decayed steadily as the allele distance increased, while the average LD of adjacent markers (~100 kb) was 0.16±0.25. The historical N_e was calculated based on the LD value across the genome and was used to represent demographic changes in SO cattle population (Figure 4; Supplementary Table S1). Results from the SNeP software started the N_e estimation from 13th generation ago, while in GONE started from more recent time, i.e. 1st generation ago. Given in the similar generation time, both analyses were resulted in different values of N_e, e.g. N_e values in 13th generation ago were 922.48 and 342 for GONE and SNeP, respectively. However, both GONE and SNeP were detect similar pattern of dramatic Ne decline at almost in the simultaneously generation time estimates, i.e. GONE at 34 to 26 generations ago and SNeP at 54 to 13 generations ago. Moreover, according to the GONE estimation results, the values of N_e from the 25th generation ago to the current generation were stable with a slight increase.