Runs of homozygosity analysis and genomic inbreeding estimation in Sumba Ongole cattle (Bos indicus) using a BovineSNP50K BeadChip

Background and Aim: Runs of homozygosity (ROH) is a biocomputational technique for identifying homozygous regions in the genomics of livestock. This study aimed to determine the ROH in Sumba Ongole (SO) bulls (n = 48) using the BovineSNP50K BeadChip. Materials and Methods: GenomeStudio 2.0 software was used to generate the BovineSNP50K BeadChip output. The ROH and ROH-based inbreeding coefficients (FROH) were determined using the detect RUNS R v4.1.0 package. Using the following filtering criteria, PLINK v1.90 software was used to perform genotype quality control: (1) Individuals and single-nucleotide polymorphism (SNPs) had call rates >0.95; (2) more than 0.05 was the minor allele frequency; (3) the list contained only SNPs linked to autosomes; and (4) SNPs that strongly deviated (p < 1e-6) from Hardy–Weinberg equilibrium were removed. Subsequently, 25,252 autosomal SNP markers were included in the ROH and FROH analyses. Results: In general, the number and length of ROH segments in pool animals were 149.77 ± 16.02 Mb and 486.13 ± 156.11 Mb, respectively. Furthermore, the ROH segments in the animals under study can be discriminated into two classes of 1–4 Mb (83.33%) and 4–8 Mb (16.67%). Subsequently, Bos taurus autosomes (BTA) 1, BTA6, and BTA14 had significant homozygous segments comprising 13 genes. Despite this, the average FROH in pool animals was 0.20 ± 0.06. Conclusion: These findings indicate that a recent inbreeding event in SO cattle occurred many generations ago. Furthermore, the candidate genes identified from the ROH analysis indicate phenotypic attributes associated with environmental adaptation and economic traits.


Introduction
Sumba Ongole (SO) cattle, a breed of Bos indicus that is indigenous to Indonesia, have demonstrated strong adaptation capabilities on Sumba Island.The introduction of this cattle breed can be traced back to the year 1900, when the breed was brought from India by the Dutch colonial government for use as a draught animal resource [1].SO bulls can reach an adult weight of 474.08 ± 25.98 kg [2] and a carcass weight of 264.06 ± 14.72 kg at the same adult weight [3].A genome-wide association study (GWAS) was conducted in SO cattle to determine candidate genes and population structure in these cattle [4,5].However, the assessment of GWAS for runs of homozygosity (ROH) in SO cattle has not been reported.ROH represent lengthy uninterrupted segments of homozygous genetic material within the genome, originating from a combination of two identical haplotypes inherited from a common progenitor [6].Moreover, ROHs serve as crucial resources for investigating genome architecture, particularly in relation to alleles that contribute to genetic enhancement in livestock [7].The presence of homozygous genomic segments may be influenced by various factors, such as intense selection pressure, historical population dynamics, and consanguinity levels [8].Recently, ROH analysis has been employed to calculate the ROH-based inbreeding coefficients (F ROH ) of livestock species [9].
In general, the inbreeding coefficient of livestock is calculated from pedigree data through statistical analysis.However, genome sequences comprising many single-nucleotide polymorphism (SNP) loci have been used to estimate the inbreeding coefficient when the pedigree data are not available [10].Theoretically, an individual's inherited allele throughout a genome with hundreds of thousands of loci can be used to directly infer the likelihood of an allele at a specific location [11].Hence, estimating the inbreeding coefficient using genome information yields accurate results, as reported by Nishio et al. [12].
(USA) launched the BovineSNP50K BeadChip, a genome-wide genotyping array for cattle, in partnership with the US Department of Agriculture Agriculture Research Service (ARS), the University of Alberta, and the University of Missouri.With a 49.4-Kb average probe spacing and 54,609 SNP probes, the BeadChip offers more than enough SNP density to support strong genome associations in cattle [13].Hence, this technique, which includes ROH and F ROH analyses, is important for genomic selection in cattle.Recently, ROH and F ROH analyses have been performed in many breeds of cattle, including Fleckvieh [14], Gyr [8], Wagyu [15], Friesian Holstein [16], Creole [17], Kazakh White-headed, Auliekol [18], and Sahiwal [19].Since 1900, SO cattle have been kept by farmers on Sumba Island in Indonesia.Therefore, most SO cattle were managed using an extensive management system without recording.Furthermore, none of Indonesia's numerous artificial insemination centers have generated frozen sperm (straw) from SO bulls.In addition, no scheme is available for sperm distribution.As a result, the inbreeding coefficient of SO cattle can be increased annually, thereby reducing their productivity.
This study aimed to analyze ROH and F ROH levels in SO cattle using the BovineSNP50K BeadChip.The results of this study are important for the genomic selection of SO cattle in the future.

Ethical approval
The study protocol was approved by the Animal Ethics Committee of the Indonesian Agency for Agricultural Research and Development (Permit number: Balitbangtan/Lolitsapi/Rm/08/2018).

Study period and location
The study was conducted from June to August 2018 at the SO cattle feedlot farm and the Bubulak slaughterhouse in Bogor, West Java, Indonesia.

Detection and classification of homozygosity runs
ROH were computed for each animal using a sliding window in R-4.4.1.software(https://cran.r-project.org/bin/windows/base/).The ROH were calculated using the following parameters [21]: (1) One missing SNP was permitted in the ROH and up to one possible heterozygous genotype; (2) 20 was the minimum number of consecutive SNPs that made up a ROH; (3) there was a minimum density of one SNP every 100 kb; and (4) there was a maximum gap of 1 Mb between consecutive SNPs.This study defined an ROH segment with 100 or more consecutive SNPs as a homozygous segment.Zambrano et al. [16] and Liu et al. [22] recommended that all identified ROH segments be classified into two distinct classes of 1-4 Mb and 4-8 Mb.

Common ROH and gene annotation
The number of times each SNP appeared in the ROH was considered and normalized by dividing it by the total number of animals in the analysis to identify genomic areas with high homozygosity.The values were plotted against the SNP location on the chromosome was performed [23].If not, neighboring SNPs that were more than this criterion eventually combined to form genomic regions of ROH segments, which are shared by most people in the population [24].The genes in the ROH region were then annotated using the bovine genome of Bos indicus (GCF_000247795.1)from the NCBI database (https:// ncbi.nlm.nih.gov).Every annotated gene in the ROH region has a biological function deduced from numerous precise literature searches.
Means and standard deviations were calculated using Statistical Package for the Social Sciences 16.0 for Windows (IBM SPSS, NY, USA) packages for the length of ROH segments, total length of ROH segments, genome coverage, number of SNPs per ROH segment, and number of ROH segments.

F ROH
The F ROH was computed in each animal using the formula mentioned by Bjelland et al. [25] as follows: F ROH is the genomic inbreeding coefficient; L ROH is the entire length of ROH in the genome of an animal.
L Auto is the length of the autosomal genome, i.e., 2,671,695,104 bp [26].

ROH
The mean length and number of ROH segments per animal are presented in

Genomic inbreeding and detection of associated genes
The average F ROH of the animals used in this study ranged from 0.10 to 0.40, as shown in Table-2.In general, the F ROH value in the 4-8 Mb class was higher than 1-4 Mb class (0.17 ± 0.03 vs. 0.31 ± 0.06).Therefore, the three loci of BTA had high homozygous SNPs, that is, BTA1, BTA6, and BTA14 (Table-3).However, the total length of ROH segments in BTA6 (687,618 bp) was lower than those in the BTA1 (2,247,563 bp) and BTA14 (1,691,438 bp) regions.Furthermore, a total of 38 candidate genes was detected in the highest frequency of ROH segments (ROH hotspot) and spread to the BTA1 (6 genes), BTA6 (4 genes) and BTA14 ( [23].Subsequently, the number of ROH segments in the studied animals was higher than that in Brown Swiss (94.60 ± 11.60), Marchigiana  [7,23].
The length of the ROH segment can be classified as small (1-4 Mb), moderate (4-8 Mb), or large (>8 Mb), as described by Liu et al. [22].Very long ROH segments (>20 Mb) are believed to have originated from common autozygosity.In contrast, most short elements are thought to have originated from more distant ancestors, which were the result of intensive selection [21].Ferencakovic et al. [14] explained that the ROH >1, >2, >4, >8, and >16 Mb represent the presence of inbreeding from 50, 25, 12.5, 6, and 3 generations ago in cattle, respectively.In this study, most of the animals had low and moderate lengths of ROH segment (2-5 Mb), which explains the recent inbreeding from about 10-50 generations ago.In this study, the short length of the ROH segments indicated that there were no signals of strong recent selection in the SO cattle.In general, the length and number of ROH segments can be influenced by selection, inbreeding, and the type of BeadChip used for analysis [29].
In Modicana cattle (B.taurus), many candidate genes were detected in BTA1 (25 genes/39 homozygous SNPs) and BTA6 (25 genes/112 homozygous SNPs) according to ROH analysis [31].Despite this, Santos et al. [28] reported that BTA6 and BTA14 in two Italian cattle breeds (Curraleiro Pe-Duro and Pantaneiro) had high-length ROH segments, but they did not identify candidate genes in either region.Subsequently, Zambrano et al. [16] identified candidate genes for BTA1 (six genes) and BTA6 (five genes) in the Friesian Holstein cattle of Colombia based on ROH analysis.
Previous studies by Viale et al. [32], Persichilli et al. [33], Passamonti et al. [34] and Hiltpold et al. [35] reported that many candidate genes detected in Table-3 were associated with body conformation   (RAB2A), milk production (DGKG, VPS8), and reproduction (MAP3K13, WDR19) in cattle.Hence, the candidate genes selected from the ROH analysis in the present study were the main functional genes involved in economic traits and environmental adaptation in the tropical climate of Indonesia.Moreover, genetic improvements to adaptability and survivability traits are important when climate changes can influence the economic traits of cattle [36].

Conclusion
The findings of this study indicate that a recent inbreeding event in SO cattle occurred many generations ago.Furthermore, the candidate genes identified through ROH analysis indicate phenotypic attributes associated with environmental adaptation and economic traits.
Several factors may influenced this research, such as, sample size, BovineSNP50K BeadChip resolution, environmental and management factors, historical data that affect the genetic variants of the population and limit the interpretation of the results of ROH analysis and inbreeding estimates.
Assessing the whole-genome sequencing (WGS) or SNP chips with more markers, comparative studies with other Bos indicus populations, long-term studies combining genomic data with environmental and livestock management information to understand inbreeding dynamics and its impact on productivity, conservation strategies to preserve genetic diversity and reducing the risk of inbreeding are some of the scopes of research that can be carried out in the future.

Figure- 2 :
Figure-2: Relationship between the number of runs of homozygosity (ROHs) per individual (dots) and genome length covered by ROHs.

Figure- 3 :
Figure-3: Manhattan plot of the occurrence of single-nucleotide polymorphisms in runs of homozygosity among Sumba Ongole cattle.
Table-1.The length of ROH segments, total length of ROH segments, genome coverage, and number of SNPs per ROH segments in 4-8 Mb class were higher than in 1-4 Mb class.In contrast, the number of ROH segments in 1-4 Mb class was higher than in 4-8 Mb class.Overall, the SO cattle had 3.28 ± 1.10 Mb of ROH segment length, 486.13 ± 156.11 Mb of total ROH segment length, 33.28 ± 8.52 SNPs per ROH segment and 149.77 ± 16.02 ROH segments.In general, the highest number of ROHs per chromosome was observed in B. taurus autosomes (BTA)1 (562 segments), followed by BTA6 (502 segments), as shown in Figure-1.Despite this, the total length of ROH segments in the animals under study ranged from 2 to 8 Mb, and the number of ROH segments ranged from 80 to 180, as illustrated in Figure-2.Therefore, a Manhattan plot with a threshold line of 75% revealed many potential SNPs in the BTA1 and BTA14 regions Figure-3.
3 genes) regions.Despite this, a total of 56 homozygous SNPs was found at the highest frequency of ROH hotspots and distributed at BTA1 (21 SNPs), BTA6 (12 SNPs), and BTA14 (23 SNPs).Mb class) of ROH length based on the same BeadChip [8].In addition, Santos et al. [28] obtained the length of an ROH segment (4-8 Mb class) that was close to the length obtained in Curraleiro Pé-Duro (5.75 Mb with 506 ROH segments) and Pantaneiro (5.94 Mb with 62 ROH segments) based on the same BeadChip.In this study, the total length of ROH segment SO cattle was highest compared with many B. taurus cattle, such as Friesian Holstein (290.60 ± 67.20 Mb), Polish Red (142.80 ± 67.40 Mb), Limousin (180.50 ± 79.90 Mb), and Simmental (201.80 ± 99.40 Mb) based on the same BeadChip

Table - 1
: Detail of ROH segments per animal in Sumba Ongole cattle.
Chromosome Figure-1: Number of runs of homozygosity segments per chromosome in Sumba Ongole cattle.