Biodiversity of Russian Local Sheep Breeds Based on Pattern of Runs of Homozygosity †

: Russian sheep breeds traditionally raised in speciﬁc environments are valuable parts of sociocultural heritage and economic component of the regions. However, the import of commercial breeds negatively inﬂuences the population sizes of local sheep populations and might lead to biodiversity loss. Estimation of the runs of homozygosity (ROH) in local sheep genomes is an informative tool to address their current genetic state. In this work, we aimed to address the ROH distribution and to estimate genome inbreeding based on SNP data to evaluate genetic diversity in Russian local sheep breeds. Materials for this study included SNP-genotypes from twenty-seven Russian local sheep breeds which were generated using the Illumina OvineSNP50 BeadChip ( n = 391) or the Illumina Ovine Inﬁnium HD BeadChip ( n = 315). A consecutive runs method was used to calculate ROH which were estimated for each animal and then categorized in the ROH length classes. The ROH were found in all breeds. The mean ROH length varied from 86 to 280 Mb, while the ROH number ranged from 37 to 123. The genomic inbreeding coefﬁcient varied from 0.033 to 0.106. Our ﬁndings provide evidence of low to moderate genomic inbreeding in major local sheep populations.


Introduction
The introduction of high throughput arrays for single nucleotide polymorphisms genotyping has led to the development of new bioinformatic approaches, which allow evaluation of genetic diversity more fully and address demographic history of the mammalian species. For example, estimation of genomic inbreeding and the analysis of patterns of distribution of runs of homozygosity regions in the genome are gaining popularity among the geneticists and are used in addition to classical methods to assess genetic processes in the populations.
Runs of homozygosity are contiguous stretches of homozygous loci that inbred offspring inherit from both parents originated from a common ancestor [1,2]. The number and length of ROH reflect individual demographic history and evaluate the homozygosity burden [3,4]. The length of ROH indicates whether inbreeding was recent or ancient in a population [1,2].
Nonetheless, livestock breeding practices often use selection schemes involving inbreeding as a tool to stabilize useful traits in farm animals. Thus, to trace the inbreeding events, the genome scanning for ROH segments was performed in various livestock species including cattle [5] (for example, Angus, Charolais, Hereford, Holstein, Simmental [6], Russian Kholmogory and Yaroslavl breeds [7]) and small ruminant species (for example, German White-headed Mutton sheep [8] and African native goats [9]).
The patterns of ROH distribution were analyzed in local and commercial sheep breeds which were selected for various purposes, inhabit different environments, and are kept under diverse production systems. Al-Mamun et al. [10] performed a search for ROH segments in Australian populations of Border Leicester, Merino, Poll Dorset, and their crosses and found that Border Leicester sheep were characterized by a higher genome coverage by ROH. In addition, analysis of ROH distribution was used to elucidate the demographic history of the six commercial meat breeds including Belclare, Beltex, Charollais, Suffolk, Texel, and Vendeen [11]. Based on estimation of the genomic inbreeding coefficient based on ROH (FROH), He et al. [12] suggested that Chinese Merino had the lowest levels of inbreeding.
However, more precise attention was paid to estimation of genomic inbreeding based on ROH in populations of local sheep. Such populations often lack reliable pedigree information, and according to Purfield et al. [11], ROH might be recommended as a predictor of the pedigree inbreeding coefficient (correlation 0.62). Mastrangelo et al. [13] investigated the occurrence of ROH in 21 Italian sheep breeds using medium-density SNP genotypes and found that Barbaresca, Leccese and Valle del Belice breeds have been affected by recent inbreeding events. Signer-Hasler et al. [14] found a high correlation (0.95) between genomic inbreeding coefficients based on ROH (FROH) estimates from mediumdensity data and HD data in eight local Swiss sheep breeds. Low genomic inbreeding was observed in the Kyrgyz local sheep breeds including Alai, Aykol, Gissar, and Tien-Shan [15]. Abied et al. [16] suggested that some animals have experienced recent inbreeding events in Chinese indigenous sheep breeds. Using Illumina OvineSNP50 BeadChip, Dzomba et al. [17] analyzed ROH distribution in 400 animals from South African sheep populations representing mutton, pelt and mutton and wool dual-purpose breeds, as well as indigenous non-descript breeds and contributed to the better understanding of the genomic landscape of African sheep breeds.
After the severe crisis caused by the USSR collapse, Russian sheep breeding, which has been focused mainly on wool production for several decades, became unprofitable due to weak demand and low prices [18]. Thus, over the eighteen-year period, the share of fine wool sheep breeds decreased by 20.9%, the population number of semi-fine wool breeds reduced by 2.3 times, and the share of coarse wool breeds increased by 5.4 times [19].
Besides an increase in the population number of unproductive sheep may lead to drastic consequences in sheep farming by financial ruining of farmers and smallholders because sale prices for wool and mutton do not cover the costs for keeping sheep [18]. Considering that majority of sheep rearing farms are in the geographical areas of underdeveloped or risky farming [20], the recessions in the regional sheep industry have notable negative consequences for local people. In this regard, along with the economic aspect, sheep farming is of significant social and cultural importance in Russia.
Summarizing, contemporary Russian sheep breeding should be focused on increasing the meat production and using low-cost technologies. The import of commercial meat breeds was not beneficial because the specific harsh feeding and keeping conditions did not allow these sheep to realize full genetic potential. In addition, development of new hybrid breeds was expected to have a positive effect on the sheep rearing industry; however, this prediction did not come true [21].
Sustainable use and management of existing local breeds is a topic priority for the rising national sheep industry. Genetic resources of Russian local sheep include breeds, which were specifically selected for wool and dual-purpose production (wool and meat), and autochthonous breeds, which are adapted to extreme environments and from which all types of sheep products are used by local smallholders [22,23].
However, the levels of genetic diversity of local breeds, including the addressing of inbreeding events, should precede the changes in the selection direction to design scientificbased breeding programs. The evaluation of genetic diversity by calculating heterozygosity and effective population sizes in the most popular Russian local breeds was not a very informative tool [24].
In this regard, the aim of our present study is to address the distribution of the runs of homozygosity and to estimate genome inbreeding in Russian local sheep breeds based on SNP-genotyping data for better understanding of current levels of genetic diversity in these valuable livestock resources.

Samples and Genotyping
Samples for this study included SNP-genotypes of 706 individuals from twentyseven Russian local sheep breeds which were genotyped using the OvineSNP50 BeadChip (Illumina, San Diego, CA, USA) or the Ovine Infinium HD BeadChip (Illumina, San Diego, CA, USA) [20]. The 50k SNP genotypes for 391 samples were generated in our previous research [24]. Additional 315 samples were genotyped using the Ovine Infinium HD BeadChip [25]. The details on the used data collection are given in Table 1.

Quality Control
Genotype quality control (QC) procedures were performed using PLINK v1.90 [26]. To consider the accuracy and efficiency of SNP genotyping, valid genotypes for each SNP were determined by setting a cut-off of 0.5 for the GenCall (GC) and GenTrain (GT) scores [27]. Samples that did not meet the quality criteria (missing genotype call rate 0.1) were eliminated from the analysis.
After merging the genotypic data from the 600K and 50K arrays, a total of 42,230 autosomal SNPs that overlapped between the two DNA chips were left in the analysis. SNPs with a call rate below 0.90, a minor allele frequency (MAF) lower than 0.05, or those which were located on sex chromosomes were eliminated from the analysis.

Principal Component Analysis (PCA)
Principal component analysis (PCA) was performed in PLINK v1.9 and visualized with the R package "ggplot2" [28]. The PCA was performed before the analysis of the runs of homozygosity distribution and showed that the individuals from the Buubei breed were divided into two groups (Buubei_1 and Buubei_2). The final sample numbers are shown in Table 2.

Runs of Homozygosity Estimation
For ROH calculation, we used a window-free method for consecutive SNP-based detection [29] implemented in the R package "detectRUNS" [30]. One SNP with a missing genotype and up to one possible heterozygous genotype was allowed in the run. The minimum ROH length was 1000 kb.
ROH were estimated for each animal and then categorized in the corresponding ROH length classes:

Estimation of Genomic Inbreeding (FROH)
The genomic inbreeding coefficient based on ROH (FROH) was estimated as the sum of the length of all ROH per sheep as a proportion of the total autosomal SNP coverage (2.44 Gb).

Assessment of Genetic Links between Sheep Breeds within the Wool-Type Groups
Principal Component Analysis performed for 15 coarse wool sheep populations ( Figure 1A) showed that the first Principal Component (PC1) accounting for 13.34% of genetic variability clearly separated the Romanov and the Kuchugur breeds from the other populations which were clustered together. The Buubei (2) population was differentiated from the Buubei_1 group as well as the other breeds by the second Principal Component (PC2). Besides the PC2 pierced the joined cluster into two slightly traceable subgroups (Buubei (1) + Mongolian + Edilbai + Kalmyk + Karakul + Tuva and Andean + Karachaev + Tushin + Ossetin + Lezgin).  Table 1.  Table 1.
Based on PCA results for semi-fine wool breeds, the Russian longhaired breed was the most distant from the other breeds and located in the right down quadrant formed by the PC1 and PC2 ( Figure 1B). The Kuibyshev and North Caucasian breeds occupied the left upper quadrant, while the Altai Mountain and Tsigai breeds were placed within the left down quadrant.
Within the fine wool group, the Volgograd breed was separated from the other breeds be the PC1, which accounts for 3.02% of genetic variability, while the Dagestan Mountain breed was differentiated by the PC2 (Figure 1C).

Pattern of Distribution of Runs of Homozygosity in Populations of Russian Local Sheep Breeds
The ROH segments were identified in all studied breeds on all autosomes. In all studied breeds, the highest genome coverage by ROH was found on Oar1 (10.58-12.31% in coarse wool, 9.49-12.24% in semi-fine wool, and 9.69-11.68% in fine wool group), Oar2 (8.88-12.17% in coarse wool, 10.53-12.52% in semi-fine wool, and 10.16-12.83% in fine wool group), and Oar3 (8.33-9.74% in coarse wool, 5.42-10.09% in semi-fine wool, and 8.98-10.55% in fine wool group). The Oar 26 was characterized by the lowest coverage by ROH (≤2.40% in all breeds).
Mean ROH lengths varied significantly in different sheep breeds ( The greatest mean ROH number was found in the Romanov breed (123.14) while the lowest one was detected in the Buubei_2 breed (56.58) within the group of coarse wool breeds. The mean ROH number varied from 37.64 in the Altai Mountain breed to 84.31 in the Russian longhaired breed within the group of semi-fine wool breeds. The maximum ROH number was estimated in the Volgograd breed (77.38), and the minimum was found in the Baikal fine-fleeced breed (52.53).
Among all studied breeds the maximum individual ROH length was found in the Kuchugur breed (872.75 Mb), and the minimum was identified in the Lezgin breed (48.78 Mb). Considering the individual ROH numbers, the greatest number was displayed in the Romanov breed (146) and the lowest number was detected in the Altai Mountain breed (26) (Table 2, Figure 2).

Ranging the Runs of Homozygosity by the Length Classes in Russian Local Sheep Breeds
A genetic pattern of predominance of the shortest ROH segments (1-2 Mb) was found in all studied sheep populations ( Figure 3). Thus, the frequencies of the 1-2 Mb ROH segments were higher 80% (with maximums in 91% in the Lezgin and 93.59% in the Mongolian breeds ( Figure 3A) within the coarse wool breeds, higher 66.82% (with maximums in 77.79% in the Baikal fine-fleeced and 78.39% in the Groznensk breeds ( Figure 3B) within fine wool breeds, and higher 50.70% within semi-fine wool breeds with maximums in 68.69% in the Altai Mountain and 68.90% in the Tsigai breeds ( Figure 3C). The Romanov (71.96%) and Kuchugur (67.73%) breeds had lower frequencies of the shortest ROH segments in comparison with other coarse wool populations.    The 2-4 Mb ROH segments were more distributed in semi-fine wool (from 19.73% in the Altai Mountain to 29.76% in the Russian longhaired breeds) and in fine wool breeds (from 16.64% in the Groznensk to 22.25% in the Volgograd breeds) than in coarse wool breeds (from 5.30% in the Mongolian breed to 21.46% in the Kuchugur breeds). A similar pattern was found in the 4-8 Mb ROH length class (7.66-13.60% in semi-fine wool group and 2.54-10.59% in fine wool group versus 0.62-7.20% in coarse wool group).
The classes of the longer ROH segments (8)(9)(10)(11)(12)(13)(14)(15)(16) Mb and >16 Mb) were less frequent in all studied breeds. The highest share of 8-16 Mb ROH segments was found in semi-fine wool breeds (2.55% in the Tsigai to 4.86% in the Russian longhaired breeds) while the maximums in coarse wool and fine wool breeds did not exceed 3.95% and 2.24%, respectively. The frequencies of the longest ROH segments (>16 Mb) varied from 0.12% in the Tuva to 2.10% in the Kuchugur breeds in coarse wool group, from 0.38% in the Altai Mountain to 1.13% in the Kuibyshev breeds in semi-fine wool group, and from 0.05% in the Groznensk to 1.36% in the Stavropol breeds in the fine-wool group.   Among all studied animals, the maximum individual FROH value was calculated in the Kuchugur breed (FROH = 0.33) from the coarse wool group. The highest individual FROH values were found in animals from the Stavropol breed from fine wool group (FROH = 0.16), from the Russian longhaired breed from semi-fine wool group (FROH = 0.16) and from the Romanov breed from coarse wool group (FROH = 0.17).

Discussion
Estimation of genetic diversity in local livestock species is of special priority to prevent the steady inbreeding increase which leads to negative consequences of inbreeding depression and to endangered status. There are several approaches to address biodiversity and its dynamics in the populations of livestock species: effective population size, heterozygosity and runs of homozygosity [10].
In our previous study, we calculated and analyzed effective population sizes and heterozygosity to unlock the current state of genetic diversity in Russian local sheep breeds [24]. Nonetheless, estimation of runs of homozygosity is a useful tool to reveal the presence of long-term inbreeding in livestock populations [2,5].
A strong primary subdivision of the Russian local sheep populations according to their wool type (fine wool, semi-fine wool, and coarse wool) was reported based on using the medium density DNA arrays [24]. Therefore, in the present study, we divided Russian sheep populations corresponding to their wool type to analyze the specific patterns of distribution of the runs of homozygosity in their genomes.
The breeds in the fine wool group demonstrated a high consistency in the ROH distribution. Thus, most individuals from these breeds fit into «90 ROH number and 200 Mb sum of ROH length» pattern. These findings might be occurred because of similarities in the developmental history and of long-term underling of strong positive selection for wool production [22]. Besides Dzomba et al. [17] reported that the Merino-type breeds had similarities in ROH distribution. These results are agreed with our findings on Russian fine wool breeds which also belong to the group of Merino-derived breeds.
A common trend was not established in semi-fine wool breeds which might be divided into three groups based on ROH distribution. The first group included Kuibyshev, Altai Mountain, and Tsigai breeds. Most individuals of these breeds had 55-60 ROH numbers and 150 Mb sum ROH length. The second group represented by sheep from the North Caucasian breed was characterized by 60-80 ROH numbers and 150-220 (250) Mb sum ROH length per individual. In addition, the individuals from the Russian longhaired had the greatest ROH numbers (≥75) and ROH length (200-350 Mb) as well as the highest genomic inbreeding coefficients.
The coarse wool group included fourteen populations collected from thirteen breeds. Based on the PCA results, animals from the Buubei breed clearly separated into two groups. This pattern was not explained by the sampling locations. The history of the Buubei breed provides a tragic lesson for future generations because the valuable gene pool of this ancient native breed was lost in the Republic of Buryatia. The contemporary Buubei breed was re-introduced into the territory of Russia from a small group of animals that had been previously imported to China and which escaped the extinction in the homeland habitat [31]. However, it might be hypothesized that some re-introduced sheep were of admixed origin which resulted in establishment of a few genetic strains within the contemporary gene pool. Nonetheless, there were no significant differences between two Buubei populations in the ROH distribution.
Most sheep in the coarse wool group had similar ROH numbers and ROH lengths which were up 105 and 250 Mb, which corresponded to previously detected patterns based on high-density genotypes [32].
However, three breeds did not fit into the genetic patterns, which were characteristic for coarse wool (Romanov, and Kuchugur) and semi-fine wool groups (Russian longhaired) and displayed the highest estimates of mean ROH length (282 Mb in Romanov; 257 Mb in Russian longhaired; 223 Mb in Kuchugur).
The contemporary gene pools of studied coarse wool breeds were formed by folk selection (somehow or other) and by required adaptions to survive in severe natural environments. Nonetheless, the Romanov breed was underling a stronger selective pressure by selection individuals, which had the best pelt traits and the highest prolificacy [22]. This could result in fixing definite genome regions that related to desirable traits and might overlap with the ROH segment. Nevertheless, this assumption should be addressed more fully with a larger sample.
However, the severity of the consequences of the recent autozygosity's events on the gene pool of these three breeds is different. Due to higher resilience to feeding and keeping conditions, the Romanov breed is reared in 26 regions. Besides, the pedigree base for the Romanov breed includes twenty breeding enterprises and multipliers [33]. Thus, a rising of the genomic inbreeding might be prevented in the future by smart choices for the unrelated rams and rotation of the founder's lines within the breed.
In comparison with the Romanov breed, the state of genetic resources of the Kuchugur and Russian longhaired breeds are more unstable. The Kuchugur breed was created by intense folk selection in the Voronezh region in the second half of the 19th century [22]. An identification of highly inbred animals (FROH = 0.33) in the Kuchugur breed was expected because this breed is in endangered status (no official census recordings are available) and most likely only a few sires are used to multiply the last existing flocks which are kept by the smallholder farmers in the Voronezh and Kursk regions. Considering Russian longhaired breeds, a rising demand for mutton has contributed to revived interest in raising breeds to produce meat (Kuibyshev, Altai Mountain, North Caucasian breed and Tsigai). However, the Russian longhaired breed has long crossbred wool (in Lincoln type), which is not in high demand currently. Thus, higher inbreeding level in this breed might correspond to the small population size (1400 heads at the end of 2019, Supplementary Material Table S1) [33] and to using of a limited number of rams. Thus, the Kuchugur and Russian longhaired breeds without proper management might be extinct in the nearest future.
Nevertheless, common genetic patterns were found in all studied Russian breeds as well. A prevalence of short ROH segments detected in all Russian local sheep populations are compatible with the relevant patterns identified in other local and cosmopolitan sheep breeds. Thus, Border Leicester, and Poll Dorset breeds predominantly had short ROH segments (1 to 5 Mb) [10] as well as Italian local sheep breeds were characterized by the highest numbers of short ROH segments (<10 Mb) [13]. Analyzing the ROH distribution in genomes of South African sheep breeds, Dzomba et al. [17] showed that 88.2% of identified ROH were in the short (1-6 Mb) category [17]. A similar pattern was reported in five Chinese sheep breeds [16] and six commercial meat breeds including Suffolk, and Texel [11].
Nonetheless, chromosome coverage in ROH varied in different sheep populations. Thus, Abied et al. [16] showed the highest coverage rate on OAR2 in Chinese sheep populations, which corresponded with our findings. Purfield et al. [11] reported the highest and the lowest percentage of the autosome residing in a ROH on OAR15 and on OAR24 in the Charollais and Suffolk populations. Dzomba et al. [17] found that South African sheep breeds were characterized by even ROH distribution amongst chromosomes.
Considering genome coverage in ROH, Russian sheep populations displayed greater mean ROH length (86.77-282.15 Mb) and higher mean ROH number (37.64-123.14) in comparison with those estimated in Italian (3.85-5.51 Mb and 10.58-44.54) [13] and in eight Swiss sheep breeds (1.88-103.25 Mb and 6.58-29.14) [14]. In addition, mean ROH length calculated in our study was a bit higher than those obtained in commercial sheep breeds (92.61-128.31 Mb [11] and 94.88-126.06 Mb [10]). However, larger variation was observed in Chinese sheep breeds for which ROH number ranged from 259 to 796, and mean ROH length varied from 15.23 Mb to 46.8 Mb with individuals values up to 273 and 984 Mb [16].
Variation of ROH lengths within the studied breeds was from moderate to high. In addition, animals with large ROH coverage were found in several breeds including Kuchugur (872.75 Mb), Romanov (457.36 Mb), Ossetin (463.41 Mb), and Buubei (2) (483.11 Mb). However, in general, maximum individual ROH length values estimated in our study were close to those obtained in Australian populations of Border Leicester, Merino, and Poll Dorset breeds (427.2, 410.5 and 396.45 Mb) [10]. Besides, the presence of several individuals, which had experienced recent autozygosity events, is typical in livestock species [13].
Although standard deviation values revealed high variability in autozygosity levels within each population, genomic inbreeding coefficients estimated per breed predominantly demonstrated a pattern of low to moderate inbreeding in Russian sheep populations (FROH from 0.033 to 0.106). Comparable results were observed in Italian (FROH from 0.016 to 0.099) [13] and Swiss local sheep breeds (FROH from 0.021 to 0.102) [14], while most South African breeds exhibited much more high inbreeding levels (FROH from 0.10 to 0.31) [17].

Conclusions
Here, we presented a detailed analysis of the pattern of the runs of homozygosity distribution in twenty-seven Russian local sheep breeds based on SNP profiles. The results corresponded to breed history and used production system under which populations are reared. The calculated levels of ROH reflect the inbreeding history of the studied sheep populations. Our findings provide evidence of a low to moderate genomic inbreeding in major local sheep populations. The results suggest that several animals from Kuchugur, Romanov, Ossetin, and Buubei breeds have experienced recent autozygosity events. The study results provide useful information and might contribute to designing conservation programs for local genetic resources of sheep in Russia.