Distribution of Runs of Homozygosity and Their Relationship with Candidate Genes for Productivity in Kazakh Meat–Wool Sheep Breed

Increasing the fertility of sheep remains one of the crucial issues of modern sheep breeding. The Kazakh meat–wool sheep is an excellent breed with high meat and wool productivity and well adapted to harsh conditions. Nowadays, runs of homozygosity (ROHs) are considered a suitable approach for studying the genetic characteristics of farm animals. The aims of the study were to analyze the distribution of ROHs, describe autozygosity, and detect genomic regions with high ROH islands. In this study, we genotyped a total of 281 Kazakh meat–wool sheep using the Illumina iScan® system (EquipNet, Canton, MA, USA) via Ovine SNP50 BeadChip array. As a results, a total of 15,069 ROHs were found in the three Kazakh meat–wool sheep populations. The mean number of ROH per animal across populations varied from 40.3 (POP1) to 42.2 (POP2) in the category 1+ Mb. Furthermore, the number of ROH per animal in ROH1–2 Mb were much higher than ROH2–4 Mb and ROH8–16 Mb in the three sheep populations. Most of individuals had small number of ROH>16 Mb. The highest and lowest genomic inbreeding coefficient values were observed in POP2 and POP3, respectively. The estimated FROH presented the impact that recent inbreeding has had in all sheep populations. Furthermore, a set of interesting candidate genes (BMP2, BMPR2, BMPRIB, CLOCK, KDM2B, TIAM1, TASP1, MYBPC1, MYOM1, and CACNA2D1), which are related to the productive traits, were found. Collectively, these findings will contribute to the breeding and conservation strategies of the Kazakh meat–wool sheep breed.


Introduction
Sheep breeding is a traditional branch of livestock breeding in Kazakhstan [1].The availability of vast areas of pasture lands and the centuries-old experience of the Kazakh people in sheep breeding determines the priority of development of this sector in the country.Increasing the fertility of sheep remains one of the crucial issues of modern sheep breeding, since the level of productivity, the production of cheap lamb and wool, as well as the profitability of the industry, depend on its positive solution.The Kazakh meat-wool sheep have high meat and wool productivity, are well adapted to the extreme conditions of semi-desert and desert zones of south-eastern Kazakhstan.To increase the fertility of Kazakh meat-wool sheep, the rams of Romanov and Finn landrace sheep were used, and Genes 2023, 14,1988.https://doi.org/10.3390/genes14111988https://www.mdpi.com/journal/genes the obtained offspring were bred inter se.A herd of sheep with high fertility was created as a result of targeted selection and breeding.In 2011, an intrabreed multiparous type of the Kazakh meat-wool breed was approved [2,3].One of the most promising methods of improving productivity livestock is the use of molecular genetic analysis, the significance of which has been proved by numerous studies [4].There is also an active search for promising marker genes in sheep breeding.In this regard, the growth differentiation factor 9 (GDF9), bone morphogenetic protein receptor IB (BMPRIB), and bone morphogenetic protein 15 (BMP15) genes that influence reproductive qualities are of great interest [5].
It is known that close inbreeding is happening during breed formation.It is a natural and inevitable phenomenon in animal populations [6].Meanwhile, due to a high degree of inbreeding some genetic characteristics, such as selection signatures, specific haplotypes, selective sweeps, and selection pressure, are fixed in a population.A decrease in the productivity of inbred animals is called inbred depression [7].In addition to providing a more powerful measure of inbred depression, genomic data offer additional opportunities for understanding the genetic background of inbreeding.Nowadays, runs of homozygosity (ROHs) are considered a suitable approach for studying the genetic characteristics of farm animals [8].ROHs are continuous homozygous regions that are present in an animal's genotype because both parents passed the same haplotypes to their offspring [9].When two identical alleles at a locus come from a common ancestor through non-random mating (inbreeding), the genotype is said to be autozygous, or identical by descent (IBD).While two identical alleles come from different sources, the genotype is called allozygous, or "identity by state" (IBS).ROH methods allow for the more accurate identification of alleles at the IBD locus and are widely used in human and animal studies to accurately estimate the level of autozygosity [10].ROH length can be used to study the effect of inbreeding age on inbred depression in addition to existing pedigree-based methods.Also, calculations using genotype information are more accurate than pedigree information.The length of the ROH pattern can be used to track the ontogeny of the autozygous segments: short ROH segments reflects an older inbreeding origin, as opposed to long ROH segments, reflecting a recent inbreeding origin [11].For instance, calculating the genomic inbreeding levels from ROHs has shown an accurate estimator for cattle population [12].
Selection traits are the result of genotypic changes in populations subjected to some form of selection pressure [13].These changes are characterized by an increase in the frequency of alleles in one or more genes or clusters of genes involved in the processes of population adaptation to specific conditions or in the improvement of productivity.Methods that allow for the identification of such traits may open up opportunities for the identification of genes involved in these processes [14].Given the increasing availability of genomic information, particularly single-nucleotide polymorphism (SNP) data, genome-wide homozygosity can be better predicted and therefore the negative consequences of homozygosity compared to "pedigree" breeding can be better reflected [15].Recently, several studies revealed genes in sheep within the ROH islands [16][17][18].For instance, Li et al. (2022) detected several interesting candidate genes, including BMPRIB in famous high-prolific Hu sheep [18].The IGF1, TGFBR2, etc. genes that affect muscle development have been reported in Chaka sheep based on the ROH analysis [19].Also, ROH outcomes displayed distinctness because of increased inbreeding among the Spanish Merino sheep [20].
Therefore, the objective of this research was to detect ROH features in Kazakh meatwool sheep population and select candidate genes within the ROH islands that are related to breed-specific traits in tested populations.

Sample Collection
A total number of 281 Kazakh meat-wool sheep bred in "Kuatzhan" (n = 95 female (n = 95♀(POP1)) and "Dukeev" (n = 115 female (POP2) and 71 male (POP3)) farms were used in this work.The two named farms have been working on improving sheep fertility for a long time in the animal husbandry sector of our country.Therefore, it is of interest to breeders to study their productivity and fecundity.Animals were selected for the study on the basis of their fecundity.Ewes with twin or more twin births and lambs were selected for this study.The blood samples were collected by the veterinarian into EDTA-containing tubes to prevent blood clotting.All animal experiments were approved by the Local Ethics Committee of the RSE Institute of Genetics and Physiology SC MSHE RK (19 October 2021, Almaty, Kazakhstan).The specimens were transported to the Laboratory of Animal Genetics and Cytogenetics of the RSE Institute of Genetics and Physiology using the blood bank refrigerator and stored in a freezer (at −25 • C) until DNA extraction.

DNA Extraction and SNP Genotyping
Genomic DNA was extracted from the whole blood using a GeneJET Whole Blood Genomic DNA Purification Mini Kit, #K0782 (Thermo Scientific, Waltham, MA, USA).This kit utilizes silica-based membrane technology in the form of a convenient spin column.The standard procedure takes less than 20 min following cell lysis and yields purified DNA more than 30 kb in size.200 µL of whole blood was mixed with 20 µL of Proteinase K Solution by vortexing.After adding 400 µL of Lysis Solution and vortexing, the tubes were incubated at 56 • C for 10 min using a thermomixer.Then, 200 µL of ethanol (96-100%) was added to the tubes and mixed by pipetting.After transferring the prepared mixture to the spin column, the tubes were centrifuged at 6000× g for 1 min.Wash Buffer I and Wash Buffer II in 500 µL were used to clean the purification column from unnecessary cellular debris.It is crucial to centrifuge tubes at high speed (≥20,000) after every step of the washing process.In the final step, 200 µL of Elution Buffer was added to the center of the column membrane to elute genomic DNA.To increase the output DNA concentration, it is necessary to reduce the amount of Elution Buffer added.The presence of DNA was checked by 1% agarose gel electrophoresis and the quality of genomic DNA was assessed by Nanodrop One (Thermo Scientific, Waltham, MA, USA) spectrophotometer.The DNA concentration was diluted to 50 ng/µL before genotyping.Genotyping was performed using the Illumina iScan system (EquipNet, Canton, MA, USA), which supports the rapid, sensitive, and accurate imaging of Illumina BeadChips for genetic analysis results, and the OvineSNP50 BeadChip microarray, consisting of more than 54,000 SNPs to deliver the densest coverage available for the ovine genome.Further, the ped and map format files were generated using the PLINK Input Report Plug-in parameter implemented in GenomeStudio v2.0 software.Here, the ped files contained the information on the family and sample IDs, paternal and maternal IDs, sex of the animals, and SNPs information, while the map files contained the chromosome ID, unique SNP identifier, genomic distance, and SNP position information.

Quality Control of Data
The primary processing of genotyped data was conducted in the GenomeStudio software.Further quality control of the data was carried out using Plink v1.07 software.The SNPs that unmapped contigs, located on sex chromosomes and mtDNA, were excluded.For the quality filtration of genotyped data, the following parameters were applied: missing genotypes >10%, SNPs with call rate <90%, the minor allele frequency p < 0.05, and Hardy-Weinberg equilibrium p < 0.0001.
For the identification and statistical analysis of runs of homozygosity, Plink v1.07 software and the detectRuns package were used.The following parameters and thresholds were applied for the identification of ROHs: the minimum ROH length was set to 1 Mb with 15 SNP, the window threshold to call a ROH was 0.05, one missing call was allowed per window, the sliding window size in SNPs was set to 15, the required minimum SNP density in an ROH was 1 SNP per 100 kb, and the maximum gap between two consecutive homozygous SNP was 250 kb.

Calculation of Inbreeding Coefficient
Runs of homozygosity is one of the common methods to evaluate genomic inbreeding coefficient in animals.The inbreeding coefficient on the basis of ROHs (F ROH ) for each individual was evaluated as follows: where L ROH is the total length of all ROHs, and L aut is the specified length of the autosomal genome covered by SNPs on the chip [21].Here, we used ROHs above 1 Mb in different classes to calculate F ROH for Kazakh meat-wool sheep.The literature review show that ROHS above 1 Mb are the most recurrently used to evaluate of genomic inbreeding coefficient and recent animal relatedness for farm animals.

Detection of ROH Islands
To determine of common ROH islands, which are the regions in which ROHS are most frequent in that population, the percentage of occurrence of SNPs falls inside ROHS was calculated, and it was plotted against SNP positions in all chromosomes.The top 1% of the highest occurrences SNPs were chosen as the threshold to identify the genomic regions most commonly associated with ROHS.For the visualization of high frequency of SNPs fall in the homozygous fragments, the Manhattan plot was generated.

Runs of Homozygosity Analysis
In the present study, Ovine SNP50 BeadChip array was used to characterize runs of homozygosity level in the genome of Kazakh meat-wool sheep breed.After quality control, 281 animals with 46,143 SNPs were retained for further analysis.A total of 15,069 runs of homozygosity fragments were found in the three populations of Kazakh meat-wool sheep.The distribution of ROHs across all chromosomes is displayed in Figure 1.Descriptive statistics on the number and length of ROHs in different classes are described in Table 1 for the three sheep populations.The mean number of ROHs per animal decreased with the growing length category of autozygosity in all populations.The results presented a rapid decrease in the mean number of ROHs per animal starting from 2+ Mb in the studied sheep.The average number of ROHs per individual ranged from 40.3 to 2.1 Mb at categories 1+ and 16+ in the POP1 population, while this value was varied from 42.2 to 1.6 Mb at categories 1+ and 16+ in the POP2 population.Regarding POP3, the mean number of ROHs was between 40.7 and 1.1 at classes 1+ and 16+ Mb, respectively.The estimated mean number of ROHs per animal was higher in POP2 compared to POP1 and POP3 in all categories, except 16+ Mb.The average ROH length varied from 51.4 (1+ Mb) to 17.2 (4+ Mb), values in category 8+ and 16+ Mb were higher than that of 2+ and 4+ for POP1.The highest and lowest of mean length of ROHs was observed in category 1+ and 4+ for POP2, while this value observed in classes 1+ and 4+ for POP3.

Genomic Inbreeding Coefficients
To evaluate of genomic inbreeding level in the three populations, the inbreeding coefficients on the basis of runs of homozygosity (F ROH ) in different ROH length classes was estimated (Table 2).For the estimated inbreeding coefficients above 1 Mb, the highest value was observed in POP2 (0.047), with the lowest value in POP3 (0.036).According to F ROH > 2 Mb, the genomic inbreeding level varied from 0.025 in POP2 to 0.014 in POP3.However, the highest value of F ROH > 16 Mb was observed in POP1, and the lowest was recorded in POP3.POP3 had the lowest genomic inbreeding level in all F ROH categories.In all populations, F ROH was higher in >8 Mb compared to >4 Mb.The inbreeding coefficient values were stable at the last three F ROH classes in POP1.F ROH coefficients varied from 0.047 (>1 Mb) to 0.013 (>16 Mb) for POP2, while these values were between 0.036 (>1 Mb) and 0.008 (>16 Mb) for POP3.F ROH coefficient rates rapidly declined at >2 Mb in all populations.

Gene Annotation
Next, to investigate effect of selection in the sheep populations, the occurrence of ROHs across genomes were explored.The Manhattan plot showed that the highest frequency of ROH occurrence was observed on 1, 2, 6, 10, and 25 autosomes.Further, the top of 1% of SNPs were selected for the studied sheep populations (Figure 2).We annotated all genes within potential ROH islands and found that the genes that were identified are related to reproductive (BMP2, BMPR2, BMPR1B, CLOCK, KDM2B, TASP1), carcass, or meat traits (MYBPC1, MYOM1, CACNA2D1), with the inclusion of one novel TIAM1 gene.Table 3 displays the genes within the ROH islands.

Discussion
Controlling the level of inbreeding and maintaining adequate genetic diversity is crucial for the robustness of sheep population.The estimation of the inbreeding coefficient based on pedigree information was not accurate due to incomplete pedigrees, the small

Discussion
Controlling the level of inbreeding and maintaining adequate genetic diversity is crucial for the robustness of sheep population.The estimation of the inbreeding coefficient based on pedigree information was not accurate due to incomplete pedigrees, the small number of generations, and errors.One of the strengths of ROH analysis is the ability to reliably identify long homozygous segments even at relatively low marker densities.This makes ROHS attractive for studying and understanding the role of autozygosity and as a possible tool for managing inbreeding levels for breeding programs in livestock production [22].Therefore, inbreeding evaluation coefficients based on genome-wide SNP markers by ROHs gives more rigorous results [23].The number and length of ROHs can directly reflect the genetic history of a population [24].
A total of 15,069 ROHs were found in the three Kazakh meat-wool sheep populations, and the high frequency of ROHs in these populations may be due to their restricted geographical location.Furthermore, POP2 had the highest number of short ROH fragments, indicating more ancient inbreeding events.Numerous research studies have also reported the abundance of ROHs in the short-length category in other sheep breeds [25,26].In terms of length, the average length of ROHs across population ranged from 12.5 Mb (POP3) to 54.1 Mb (POP2).The mean length of ROHs per animal (Mb) was higher in ewes compared to rams in all ROH length categories except 1+ Mb class, which is in agreement with the mean number of ROHs per animal and F ROH results.The lowest values of inbreeding coefficients (0.008-0.036) in POP3 may indicate a relatively short breeding history in this population.This is in accordance with reality; after weaning, the rams usually are separated to temporary stock to buy and slaughter.Similar outcomes were found in Spanish sheep breeds [26].Besides, ROH patterns have specificity in various chromosomes in the present study, which is consistent with previous findings [27,28].
In our study, significant homozygous regions were found in Kazakh meat-wool sheep, which may indicate that some homozygous regions are contributing to fertility complications due to inbreeding depression.The inbreeding coefficient for each of the categories (F ROH1Mb , F ROH2Mb , F ROH4Mb , F ROH8Mb , and F ROH16Mb ) was estimated.The number of segments longer than 16 Mb is much smaller, indicating that the level of inbreeding in the studied sheep population is low, and that they are unaffected by recent inbreeding events.In general, these results are consistent with previously observed results in other sheep breeds [18,29].On the other hand, recent inbreeding resulted in a long ROH (>16 Mb), which may indicate the ineffectiveness of recent genetic management attempts.In contrast, the abundance of short F ROH1Mb indicates inbreeding that occurred in ancient times.Ancient and recent inbreeding emphasize the relevance of measures to maintain high fertility and limit inbreeding in Kazakh meat-wool sheep.
Although not all identified genes are proliferation-associated genes, some important genes have been identified in this study.It has been well documented that bone morphometric proteins (BMPs) are multifunctional growth factors belonging to the transforming growth factor β (TGF-β) superfamily, whose role in sheep reproduction has been widely studied in recent years [30].Bone morphogenetic protein 2 (BMP2) plays an important role in bone formation and regeneration.BMP-2 protein is one of the most effective inducers of osteogenesis that activates the expression of osteoblast markers [31].Of interest, polymorphisms within the BMP2 gene were associated with tail length in Chinese Hu and Tibetan sheep [32].Bone morphogenetic protein receptor type 2 (BMPR2) encodes the bone morphogenetic protein receptor type II, which belongs to the transforming growth factor receptor (TGFR) superfamily.The BMPR2 receptor regulates the number of cells in tissues by triggering and controlling their proliferation and apoptosis [33].Importantly, it has been revealed that the expression of BMPR2 in cumulus-oocyte complex and ovary tissues were related to the sheep antral follicle number [34].The BMPR1B gene is well known to be play an essential role in sheep reproductive traits [35].Interestingly, the Clock gene, which is a transcription activator, was highly expressed in the cerebellum, hypothalamus, and uterus in Small-tail Han sheep [36].In addition, Lysine demethylase 2B (KDM2B) has been reported to play a role in spermatogenesis [37].
T lymphoma invasion and metastasis protein (TIAM1) encodes the invasion of T-cell lymphoma, and it is critical for the integrity of adherent junctions and cell-matrix interactions [38].Myosin regulatory light chain-2 (MYBPC1) was associated with lamb meat quality characteristics [39].Voltage-dependent calcium channel subunit α-2/delta-1 (CACNA2D1) encodes the CACNA2D1 protein, which is cleaved into the α 2-and delta-subunits of a potential-dependent calcium channel.These channels play an important role in coupling muscle excitation and contraction.The SNP mutations in this gene was significantly associated with carcass traits in cattle [40,41].Recently, the GWAS revealed that the Threonine aspartase 1 (TASP1) gene was related to the fertility and precocity in Nelore cattle [42] as well as increased fiber muscle diameter in ducks [43].Furthermore, Myomesin-1 (MYOM1) was found to be differentially methylated by gene ontology enrichment analysis in sheep [44].

Conclusions
In this study, we investigated the use of selection signature analysis within the Kazakh meat-wool sheep populations based on homozygosity.To the best of our knowledge, this is the first study identifying the genomic regions of a Kazakh meat-wool sheep breed.Moreover, genomic regions associated with reproductive traits as well as carcass and meat traits have been identified in tested population.The results obtained can be used in practice to improve the genomic characterization of this breed.

Figure 1 .
Figure 1.Distribution of runs of homozygosity fragments for each autosomal chromosome, with the length of OAR (x-axis) in Mb.Red long lines mean ROH (>5 Mb ROH), blue short lines mean (<5 Mb).

Figure 1 .
Figure 1.Distribution of runs of homozygosity fragments for each autosomal chromosome, with the length of OAR (x-axis) in Mb.Red long lines mean ROH (>5 Mb ROH), blue short lines mean (<5 Mb).

Table 1 .
Descriptive statistics of ROHS for three sheep populations.

Table 1 .
Descriptive statistics of ROHS for three sheep populations.

Table 2 .
The genomic inbreeding coefficient values within different classes.

Table 3 .
The identified top of 1% genes within the ROH islands.
Figure 2. The occurrence (%) of SNPs in ROHs across individuals.

Table 3 .
The identified top of 1% genes within the ROH islands.