High-density genotyping reveals signatures of selection related to acclimation and economically important traits in 15 local sheep breeds from Russia

Domestication and centuries of selective breeding have changed genomes of sheep breeds to respond to environmental challenges and human needs. The genomes of local breeds, therefore, are valuable sources of genomic variants to be used to understand mechanisms of response to adaptation and artificial selection. As a step toward this we performed a high-density genotyping and comprehensive scans for signatures of selection in the genomes from 15 local sheep breeds reared across Russia. Results demonstrated that the genomes of Russian sheep breeds contain multiple regions under putative selection. More than 50% of these regions matched with intervals identified in previous scans for selective sweeps in sheep genomes. These regions contain well-known candidate genes related to morphology, adaptation, and domestication (e.g., KITLG, KIT, MITF, and MC1R), wool quality and quantity (e.g., DSG@, DSC@, and KRT@), growth and feed intake (e.g., HOXA@, HOXC@, LCORL, NCAPG, LAP3, and CCSER1), reproduction (e.g., CMTM6, HTRA1, GNAQ, UBQLN1, and IFT88), and milk-related traits (e.g., ABCG2, SPP1, ACSS1, and ACSS2). In addition, multiple genes that are putatively related to environmental adaptations were top-ranked in selected intervals (e.g., EGFR, HSPH1, NMUR1, EDNRB, PRL, TSHR, and ADAMTS5). Moreover, we observed that multiple key genes involved in human hereditary sensory and autonomic neuropathies, and genetic disorders accompanied with an inability to feel pain and environmental temperatures, were top-ranked in multiple or individual sheep breeds from Russia pointing to a possible mechanism of adaptation to harsh climatic conditions. Our work represents the first comprehensive scan for signatures of selection in genomes of local sheep breeds from the Russian Federation of both European and Asian origins. We confirmed that the genomes of Russian sheep contain previously identified signatures of selection, demonstrating the robustness of our integrative approach. Multiple novel signatures of selection were found near genes which could be related to adaptation to the harsh environments of Russia. Our study forms a basis for future work on using Russian sheep genomes to spot specific genetic variants or haplotypes to be used in efforts on developing next-generation highly productive breeds, better suited to diverse Eurasian environments.


(Continued from previous page)
Conclusions: Our work represents the first comprehensive scan for signatures of selection in genomes of local sheep breeds from the Russian Federation of both European and Asian origins. We confirmed that the genomes of Russian sheep contain previously identified signatures of selection, demonstrating the robustness of our integrative approach. Multiple novel signatures of selection were found near genes which could be related to adaptation to the harsh environments of Russia. Our study forms a basis for future work on using Russian sheep genomes to spot specific genetic variants or haplotypes to be used in efforts on developing next-generation highly productive breeds, better suited to diverse Eurasian environments.

Background
The domestication and selective breeding of livestock species has led to changes in their genomes to meet the needs of humans and adapting livestock to various environments [1]. Many centuries of artificial selection shaped the genomes of contemporary sheep breeds by selective sweeps, which are associated with a variety of economically important traits, such as quality and quantity of wool, milk and meat. On the other hand, the history and development of these breeds is associated with the history of human migrations [2,3]. After domestication about 9000-11,000 years ago on the territory of modern Iran [4], sheep have been spread worldwide, accompanying human migrations, e.g., the nomad's expansion [2][3][4]. Thus, both artificial selection and acclimatization to various environments have contributed to the diversity of local breeds, expressing distinct and sometimes contrasting phenotypes (e.g., polled and horned breeds, black and white fleece, breeds adapted to the hot temperatures in Africa and the cold climates of Siberia). Therefore, domestic sheep should be considered as a valuable model species to study genome reaction to both environmental adaptation and artificial selection.
Searches for genome intervals, genes and polymorphisms which determine performance for economically important traits in sheep have resulted in revealing signatures of selection likely associated with coat colour [5][6][7], tail fat deposition [8], muscle growth [9], milk yield [10], reproductive traits [6,11], wool [12], resistance to parasites [13], presence/absence of horns [6,11], etc. These chromosome regions contain markers (genes or regulatory sequences) which should be the focus of future efforts on the improvement of local and multinational breeds using the rapidly-developing plethora of contemporary genetics tools.
Sheep have been successfully bred in different environments including some harsh ones, suggesting that the process of extreme acclimation could also lead to selective sweeps in the genomes of local breeds. Environmentinfluenced adaptations are complex, effecting multiple biochemical processes; therefore signatures of selection would be expected in genes from different pathways [14]. Thus, adaptation to high altitudes is an essential feature for the sheep breeding industry in countries with a predominating mountain terrain. Genomic studies of Tibetan sheep identified several candidate genes, associated with tolerance to hypoxia (e.g., EPAS1, CRYAA, LONP1, NF1, DPP4, SOD1, PPARG, and SOCS2). It was proposed that a key gene for adaptation to hypoxia is the EPAS1, which effects the mean corpuscular haemoglobin concentration and mean corpuscular volume [15]. A study of sheep breeds from a different region, the Himalayas, suggests that the FGF-7 is a candidate for protection against pulmonary injuries and affects efficiency of lungs in sheep that inhabit high-altitude areas [16].
Temperature regime is one of the key climate factors that influence survival and successful breeding of livestock species. Multiple genes and gene networks are affected in different environments. For instance, Kim et al. (2016) shows that in the genomes of Egyptian indigenous sheep, there are changes in multiple genes that influence adaptation to hot arid environments including genes involved in melanogenesis, body size and development, energy and digestive metabolism, as well as in nervous and autoimmune response [17]. Furthermore, Lv and co-workers (2014) identified 17 genes putatively associated with climate-driven selection when they looked at a set of native sheep populations from a worldwide range of geographic areas [18]. Nine of the genes were directly involved in energy and regulation activities (e.g., TBC1D12 and FBXO8) or encoded enzyme activators (e.g., THY1), while eight additional genes were involved in endocrine and autoimmune processes (e.g., EDNRB, NMUR1, PRL, IL12RB1, and ACVR2A). Signatures of adaptations, caused by temperature and sunlight, were linked to the TBC1D12 [18]. Comparative genomics provides additional evidence for links between genes and acclimation. For instance, TRPM8 (transient receptor potential cation channel subfamily M member 8) plays a role in thermal sensation in mice [19]. The same gene was linked to cold tolerance in sheep [6,7,11]. These examples demonstrate that different genes could be affected by selection to similar environments, suggesting that studies of multiple breeds adapted to a variety of conditions could help revealing a full set of genomic regions affected during the process of adaptation.
As a step toward this goal, we present the scan for signatures of selection in a set of 15 Russian local sheep breeds. For our current analysis we chose a subset of genetically different breeds and breeds adapted to cold climates of Siberia based on results of our previous study focused on the phylogenetic history of 25 Russian sheep breeds, performed using a medium-density ovine SNP array [20]. We added one additional breed from Siberia (putatively adapted to cold climate), to extend the dataset of cold-adapted breeds. Due to its large area and geographical position, Russia comprises a variety of diverse climatic zones present in Eurasia. Nevertheless, about three quarters of Russia is represented by territories in a continental climate with occasionally extremely low winter temperatures (-30°C and below). Because of adaptation to the diverse environments of Russia, we expect that the analysis of Russian native breeds will: a) confirm the previously reported signatures of selection, and b) point to new regions not found in other studies, which could be important for the local adaptations especially in the northern regions. We used a combination of two powerful approaches to identify signatures of selection, one based on differences in frequencies of haplotypes in sets of related breeds (hapFLK), and the second one combining various breed-specific selection statistics into a single integrated framework (DCMS). Together these methods provide a complimentary way of finding signatures of selection in the set of 15 Russian sheep breeds, which we genotyped on a high-density SNP array. We analysed a minimum of 40 autosomes per breed (20 individuals) to ensure detection of the most common signatures of selection by both methods. Our data represent the first comprehensive analysis of selective sweeps in the genomes of Russian local sheep breeds. Our results could be used by the international community to better understand genetic mechanisms of adaptation to various environments, and by breeders looking to develop new breeds that are better adapted to local conditions or to improve adaptation in extant breeds.

Breed groups
The Admixture and Principal Component Analysis (PCA; PC1) suggested the presence of two welldifferentiated clusters of breeds in our dataset ( Fig. 1): Buubei, Lezgin, Karachaev, Karakul, Tuva, Edilbai, Romanov (GROUP1) and the Russian Longhaired, Altai Mountain, Groznensk, Salsk, Volgograd, Krasnoyarsk, Baikal, and Kulundin (GROUP2). Although the Romanov breed demonstrated a substantial differentiation from the rest of the breeds we assigned it to GROUP1 because it had a higher GROUP1 component in the AD-MIXTURE (K = 2) analysis, clustered with other GROUP1 breeds based on the PCA PC1 results, and additionally represents a coarse wool breed (as all other breeds from the GROUP1). In total, 312 animals from 15 Russian local sheep breeds with a mean number of 21 individuals per breed were used in these analyses ( Table 1).
The de-correlated composite of multiple signals (DCMS) and hapFLK statistics for individual breeds and groups of breeds overlapped to some extent: 546 DCMS-detected regions, covering 100.8 Mbp of the sheep genome sequence were found in shared intervals, providing independent support for selected regions (Fig. 2, see Additional file 1). However, because the hapFLK statistic detects signatures of selection within groups of breeds and the DCMS in our study was used to combine statistics within a breed, the hapFLK results could not be added to the DCMS framework. The hapFLK revealed additional selective sweeps within groups of breeds missed by the DCMS, while the DCMS was efficient in detecting narrower selective sweeps often attributed to individual breeds.
Composite measure of selection DCMS statistics were calculated for each single nucleotide polymorphism (SNP) for each breed. After fitting for normal distribution, calculation of p-values and correction for multiple testing, we obtained from 134 to 238 genomic intervals under putative selection per breed (q-value < 0.01) with a total of 3069 regions detected across all breeds (with some overlaps between breeds; Fig. 2, see Additional file 2). The size of the genomic regions putatively under selection varied from 1 bp to 3,510,819 bp, with the average size equal to 120,924 bp. The total number of genes across all selected regions per breed ranged from 146 to 366.

HapFLK
The total number of selected regions identified by the hapFLK analysis (122; Fig. 3, see Additional file 2) was lower than found by the DCMS method. The largest number of hapFLK-detected regions was observed in the GROUP1 set of breeds (62), followed by the all-breed (33) and the GROUP2 (27) sets. The GROUP1 and GROUP2 shared five common hapFLK intervals. One region was found on ovine chromosome (OAR) 6 (34.6-37.5 Mbp) containing multiple genes with known effects on economically important traits such as milk production and growth. Three overlapping regions were detected on this chromosome: one near the UFM1 gene involved in brain development and abnormalities in humans [21]; the second, near neurobeachin (NBEA), a gene previously associated with autism in humans [22] and wool production traits in sheep [23]; and the last, near RXFP2, involved in formation of horns in sheep [24] and cattle [25]. The last overlapping region was found on OAR13 and contained BMP2, the only common gene between the groups in this region. BMP2 is associated with the body size and developmental traits in sheep [17]. Sizes of the hapFLK putatively selected regions ranged from 70 Kbp to 5.7 Mbp with an average size of 667.7 Kbp.
Candidate genes for adaptation of the Russian sheep breeds to environmental and climate challenges In the regions under putative selection we looked for genes that could be related to adaption of Russian sheep breeds to local environments ( Table 2). We noticed that five breeds: Karachaev, Krasnoyarsk, Lezgin, Salsk, and Altai Mountain had a common signature of selection reported by the DCMS method near the ZFHX2 (zinc finger homeobox 2), a gene which was shown to be responsible for an inability to feel pain in humans (including low and high temperatures) [26]. The same authors report that ZFHX2 knock-out mice had significantly higher acute thermal pain thresholds [26]. Interestingly, several other genes directly involved in hereditary sensory and autonomic neuropathies [27] often accompanied with inability to feel pain and cold temperatures in humans were found top-ranked in selected regions in the Russian sheep breeds. Among them were the SCN9A [28], reported for Buubei, WNK1 [29] reported for Salsk, HARS [30] reported for Karakul and Kulundin breeds, TECPR2 [31] reported for Tuva sheep, GAN [27] for Karachaev, and EGR2 [27] for GROUP1.
We also found several genes that were previously proposed to be related to thermal adaptations [14]. Among them was the EGFR, a membrane receptor for epidermal growth factors and a top DCMS candidate in a 33 Kbp region on OAR19 reported for the Altai Mountain sheep breed from Siberia. According to Wollenberg Valero and co-workers (2014) EGFR could represent a functional hub for the relay of thermal signalling [14]. Another gene involved in adaptation to hot/cold environment was the heatshock protein H1 (HSPH1) [14], the top-ranked gene found in a DCMS-reported region on OAR10 for the Baikal sheep. Several genes linked to climate adaptation through their interaction with the POMC (pro-opiomelanocorin receptor), shown to be involved in energy homeostasis, melanocyte stimulation and immune response [14], were top-ranked in selected regions of Russian sheep breeds. Among them are melanocortin receptors: MC1R, MC2R, and MC5R. The MC1R was the top gene in a hapFLK reported region on OAR14 for GROUP1 and in a DCMS reported region for Altai Mountain sheep. The MC2R was the top gene for a DCMS-reported region on OAR23 found in Baikal breed. The MC5R was the top gene in region found by the DCMS method in the Salsk sheep. The melanocortin receptors are involved in the movement and positioning of melanocytes, suggesting their possible contribution to adaptation to different light regimes [14], however these genes could also be involved in formation of coat colours, suggesting their contribution to economically important traits [32]. Another gene, functionally linked to POMC [14] was the MMRN1 with SNPs associated with winter duration in 54 human populations [33]. It was the top ranked gene in a hapFLK reported region on OAR6 for the all-breed analysis, and in the regions found by the DCMS method for Volgograd and Krasnoyarsk sheep.
We further found genes previously identified in a genome-wide scan for climate-mediated selective pressures in sheep [18] in the DCMS or hapFLK-reported regions. These included the NMUR1, EDNRB, and PRL all being the top genes in the DCMS/hapFLK reported intervals for the Romanov, GROUP1, and Volgograd breeds, respectively. The hapFLK method reported another region under putative selection for GROUP1 with PDGFRA found near the most significant SNP. Under cold stress, proliferation of endothelial cells and interstitial cells expressing PDGFRA increases 3-4 fold in brown adipose tissue; it is a key gene involved in cold-induced adipogenesis in mice [34]. Also, the adipocyte determination and differentiation-dependent factor 1 (ADD1) known to be involved in cold adaptation through brown adipocytes [35] was top-ranked for a region on OAR6 reported in GROUP2 by hapFLK. Consistent with the expected role of brown adipose tissue in adaptation to cold climate sheep breeds from Siberia: Kulundin, Altai Mountain, and Baikal all had a strong signature of selection (q-value ranges from 10 − 5 to10 − 7 ; DCMS) near the ADAMTS5, shown to be involved in the adiposity and metabolic health. Enhanced thermogenesis through the browning of white adipose tissue was reported in ADAMTS5 knock out mice upon cold exposure [36,37]. In accordance with the previous findings in Chinese native sheep [11] and indigenous Sunite sheep [38] we detected a signature of selection by DCMS near TSHR in Karakul and Lezgin breeds. TSHR was previously associated with metabolic regulation and photoperiod control or reproduction in vertebrates [11].

Fleece related traits
Quantity and quality of wool is one of the most economically important traits in sheep. Consistent with this, the hapFLK analysis shows a strong signature of selection (q-value < 10 − 5 ) in the group of all Russian breeds on OAR23 in the region containing a cluster of seven desmosomal genes (DSG@ and DSC@). This selection signal becomes much stronger for GROUP1 (q-value < 10 − 23 ), but the region is narrower and contains only the DSG@ genes. DSMC splits this interval in two regions for the Russian Longhaired breed, suggesting separate selection sweeps near the DSG and DSC gene clusters. The Altai Mountain breed, however, exhibits the wider selection signal covering both clusters. The DSG4 (desmoglein 4) is a known candidate for wool length and crimp in sheep [40], likely to be related to white and black coat colour in goats [41] and is a known cause of the recessive hairless phenotype in rats [42]. The DSG1 and DSG3 are associated with hair growth and follicle structure and the DSC2 with woolly/straight hair phenotype in sheep [43,44]. The DSMC method reports another strong selection signal (q-value < 10 − 7 ) near the keratin gene cluster (KTR@) on OAR3 for the fine fleeced breeds: Baikal, Kulundin, Salsk, Krasnoyarsk, Groznensk, and Volgograd. For five out of the six breeds the selection intervals are relatively wide (> 100 Kbp) and contain multiple keratin genes while the Volgograd breed contains a narrower interval (~50 Kbp) with the KRT72 gene only. Multiple members of the keratin cluster were shown to be related to fleece development. KRT5 is related to fleece development and function [45], while KRT71 is associated with the curly hair phenotype in sheep [46].

Coat colour
We identified a large group of known candidate genes related to coat colour, a trait of economic importance, also related to domestication and historical breed formation [47]. As expected, strong signatures of selection were found in the regions containing the genes KIT on OAR6, and KITLG on OAR3. The KIT proto-oncogene receptor tyrosine kinase (KIT) is associated with coat colour in various sheep breeds [5] and other species [48]. In our study, hapFLK identified the chromosomal interval containing this gene to be under strong selection in the set of all Russian breeds and GROUP1. KITLG is associated with pigmentation in sheep [49], roan coat phenotype in goats [50], and cattle [51]. The DCMS method reported this gene to be the top-ranked candidate in the selected regions of the Buubei and roan Edilbai breeds. Another gene found in a region of putative positive selection, the melanocortin 1 receptor (MC1R), was top-ranked in the hapFLK-reported region in GROUP1 and in the DCMS interval identified in the Altai Mountain breed. MC1R has a pleiotropic effect, known, among other traits, to influence coat colour in sheep [52] and cattle [53]. The microphthalmia transcription factor (MITF) is a regulator of melanocyte development, and is associated with coat colour in mouse Fig. 2 Circus plot of the relative density of selected regions along the ovine genome. DCMS statistics corresponding to GROUP1 is in dark red, to GROUP2 in dark blue, and hapFLK statistics for all breeds shown in green [54], dog [55], and other species. A strong signature of selection (q-value < 10 − 6 ) was reported by hapFLK for GROUP1 in a~350 Kbp interval on OAR19 containing MITF, and by DCMS for the Karakul breed. We identified an additional 20 genes in selected regions of Russian sheep breeds, reported by either hapFLK or DCMS (or both) methods that could be related to coat colour, fleece and other phenotypes in Russian sheep breeds (see Additional file 3).

Milk and lactation-related traits
The hapFLK method reported a large,~5.  [39]. A major milk-related trait gene, the ATPbinding cassette, sub-family G (white), member 2 (ABCG2 [56]) was found in intervals reported by DCMS for multiple Russian breeds but was not top-ranked in any of them. On the other hand, the region containing three genes ABCG2, SPP1, and PKD2 was under selection in the Baikal sheep with SPP1 (osteopontin) reported as the top-ranked gene by DCMS. SPP1 is associated with milk protein percentage, milk yield, milk protein yield, and lactation regulation in dairy cattle [57]. Signatures of selection near the SNCA, potentially related to milk protein and fat traits in dairy cattle [58] were reported by DSMC for the Kulundin, Volgograd, Krasnoyarsk, and Baikal breeds. Two family members, the ACSS1 and ACSS2 previously associated with the mammary gland function and milk fatty acid composition in sheep [59,60], cattle [61], and yaks [62] were the top-ranked genes in positively selected regions reported by DCMS. The ACSS1 gene was reported for the Volgograd, Groznensk, and Altai Mountain, and ACSS2 for the Romanov breeds. Another gene related to the mammary gland function was the ATM on OAR15, found in a DCMS-reported region for the Salsk, Baikal, and Groznensk breeds. The ATM contributes to mammary gland homeostasis and its knockout leads to a progressive lactation defect in mice [63]. In addition, we found VPS13B, a known candidate gene for milk-related traits in buffalo [64], in the interval reported by hapFLK for GROUP1. The DCMS analysis further reported this      [65,66], with the HOXA@ also being associated with body composition and structure in pigs [67] and HOXC@ being associated with tail fat deposition in sheep [68]. The region containing LAP3, NCAPG, LCORL, and HERC6 genes on OAR6 is known to be related to growth traits, carcass composition, body size, weight and height in sheep [1], horses [69], and cattle [70,71]. The results obtained using hapFLK and DCMS showed varying patterns of selective sweeps near these genes in Russian sheep breeds. Thus, according to hapFLK, LCORL was top-ranked for the all-breed analysis. This was further supported by DCMS for the Lezgin, Buubei, and Edilbai breeds while the NCAPG was top-ranked by hapFLK for the GROUP1 only. The LAP3 from the same area was top-ranked for the Baikal and Altai Mountain breeds and the HERC6 for Volgograd, Karakul, and Altai Mountain implying multiple signatures of selection in this region. The hapFLK approach reveals a~1 Mbp interval containing only one gene, the CCSER1, on OAR6 in the all-breed set as well as in GROUP2, while DCMS further confirms this signature of selection in Volgograd, Krasnoyarsk, and Baikal breeds. Variants in the CCSER1 are associated with the feed efficiency in beef cattle [72].

Reproduction
Both hapFLK and DCMS reported putative positive selection in multiple breeds for genomic regions containing genes with known effects on reproduction in mammals, including fertility: DNAJB14 [77], gonad development and sperm maturation: GNAQ [78], spermatogenesis: UBQLN1 [79], IFT88 [80], and PPP1CC [81]. The DCMS method detected selection sweeps near CMTM6 in Karakul and Edilbai sheep and near HTRA1 in Karakul, Tuva, and Buubei sheep. These genes are associated with off-season reproduction traits, such as year-around oestrous behavior in sheep (HTRA1 [82] and the evolution of sperm and the circadian rhythm systems in mammals (CMTM6 [83]). The DCMS method further reported narrow selection sweeps for multiple Russian sheep breeds in the areas of B4GALNT2 which was proposed as a strong candidate gene for ewe fertility [84].

Discussion
Here we present the first comprehensive study of signatures of selection in the genomes of 15 native sheep breeds from the Russian Federation, for most of which we recently revealed the phylogenetic and population history in the context of related breeds from other countries [20]. In contrast to our phylogenetic study, based on moderate-density SNP array genotypes, genotyping for the present study was performed on a high-density array required to reveal the majority of selective sweeps. The analysis of the data was performed using complimentary approaches: the hapFLK and DCMS, which allowed us to detect signatures of selection that are putatively related to adaptation of breeds to local environments and human needs.
More than 50% of the putatively selected regions detected in the present study overlap with previously reported signatures of selection in various sheep breeds [1,7,11,39], suggesting that our approach is robust enough to detect the expected signals of selection and that our dataset is different enough from the previously published ones to reveal new strong selective sweeps. The overlap observed between the hapFLK and the DCMS results point to putatively selected regions that are detectable using different models. The overlapping regions mainly contained known targets of selection in sheep. However, in general, hapFLK detects fewer but longer selected regions, while DSMC identifies a larger number of shorter intervals. This indicates that hapFLK could be more efficient in detecting regions with long haplotypes under selection, while DCMS can further dissect these intervals and point to specific genes under putative selection in individual breeds. DCMS is also capable of detecting shorter intervals, not found by hapFLK. This is confirmed by the fact that the region on OAR6, containing multiple genes related to milk production and growth traits was reported as a single interval by hapFLK but DCMS reported different regions within this interval containing different genes to be under selection in individual breeds. We detected the strongest signatures of selection shared between two groups of Russian sheep breeds in the regions of genes related to brain development, growth, milk production, and horned/polled phenotypes. These groups of genes are likely to be related to the process of domestication and historical breed formation.
In contrast to our recent study on the signatures of selection in Russian native cattle breeds [85] and high-density analysis of popular commercial breeds from Zhao and co-workers (2015) [86] we did not observe breeds where major milk-production related genes (e.g., DGAT1 and ABCG2) would be reported as top ranked. This could be related to the fact that the sheep breeds used for milk production in the Russian Federation (e.g., Lezgin and Karachaev) cannot be considered as strictly dairy breeds. They are also breed for wool and meat, and these traits are often more important.
Consistent with this and our previous publication [20], the PCA and ADMIXTURE analysis separated the Russian sheep breeds into two major clusters which follow wool quality type [20]. The DCMS method further suggested that the cluster of KRT genes on OAR3 is under positive selection in fine wool breeds, confirming that the keratin genes could be related to the quality of wool, as it has been shown in other studies [46]. Interestingly, the region containing a cluster of desmosomal genes was under strong selection in coarse wool breeds (GROUP1). The DCMS approach, however, points to two breeds from GROUP2, which also had this region under putative selection. GROUP2 includes breeds with the fine and semi-fine wool. According to DCMS the desmosomal genes were under selection only in the semi-fine Russian Longhaired and Altai Mountain breeds. Combined with the hapFLK results, this may indicate that selection in desmosomal genes could be related to 'coarseness' of the sheep wool but this hypothesis needs further verification using sequenced genomes from breeds of all the three types.
Unlike the Yakut cattle [87], domestic sheep cannot be normally found above the Polar Circle in Russia. However, the harsh environments of Siberia and other parts of the Russian Federation still impose significant selective pressure on local sheep breeds. Consistent with this, multiple genes related to acclimatization were found in putatively selected regions in Russian sheep breeds. In another study [85] we found one of the strongest signals of selection in the Yakut cattle, adapted to survive at − 50°C, in the area of the gene RETREG1, one of the key genes involved in the hereditary sensory and autonomic neuropathy, type II in humans [27]. This disorder in humans is accompanied by an inability to feel pain and low temperatures. In the present study, we identified multiple genes with known key contributions to human hereditary sensory and autonomic neuropathies (types I-IV) in selected regions of individual or multiple sheep breeds. It is tempting to hypothesise that changes within these genes could be related to adaptation to the harsh environments of the Russian Federation. For the breeds from Siberia, changes in these genes could be beneficial for surviving cold winters, but in other breeds they could be advantageous for survival in cold mountain climates. However, a definite answer to this question could only be obtained after the actual sequences of these genes are obtained for breeds dwelling in different environments and the presence and frequencies of missense or regulatory variants are compared.
Similar to Russian native cattle breeds, sheep breeds express signatures of selection in regions of genes related to brown adipogenesis. Brown adipose tissue is an important organ involved in in non-shivering thermogenesis indicating that these genes could be related to adaptation to cold climates. However, the adipose tissue contributes to phenotypic differences between breeds (fat-and thin-tailed sheep) and to meat quality, which is an economically important trait. This suggests that more studies need to be done to distinguish potential adaptive effects of the genes involved into adipogenesis, from their contribution to economically important traits in sheep and cattle.

Conclusions
In conclusion, we identified signatures of selection in Russian local sheep breeds of European and Asian origin. These signatures point to known regions related to economically important traits, domestication and breed formation as well as to intervals of the sheep genome that could contribute to adaptation of breeds to their corresponding local environments. Our results indicate that a detailed study(ies) of Russian local breeds involving whole-genome sequencing should focus on identifying causative genetic variants or haplotypes. These polymorphisms should be in turn the focus of future efforts on the improvement of local breeds, or for selecting multinational commercial breeds which would be better suited for the environments of the Russian Federation and Northern Eurasia.

Sample collection
Tissue samples for the Lezgin, Karachaev, Karakul, Edilbai, Romanov, Russian Longhaired, Groznensk, Salsk, and Volgograd breeds and blood samples for the Buubei, Tuva, Altai Mountain, Krasnoyarsk, Baikal, and Kulundin breeds were collected from farms and breeding centres across Russia. All samples were collected by trained personnel following strict veterinary regulations. We studied the pedigrees of animals to avoid sampling of close relatives (siblings, parents, and offspring). Tissues and blood were stored at -80°C until use.

Genotyping of Russian sheep breeds
DNA from tissue samples of nine sheep breeds was extracted using Nexttec columns (Nexttec Biotechnology GmbH, Germany) following the manufacturer's instructions. DNA from blood samples of six additional sheep breeds was extracted using cell lysation followed by phenol-chloroform extraction [88]. DNA samples of all breeds were genotyped using the Ovine Infinium® HD SNP BeadChip (600 K SNPs) to produce dense genome coverage. Genotypes were called using GenomeStudio 2 software (Illumina, San Diego, USA) and samples with overall calling rate < 95% were removed from further analyses. The produced files with genotypes (.ped) and chromosomal positions (.map) were processed using PLINK v.1.90 whole genome analysis toolkit [89].

Analysis of groups of populations
Divergence between populations can influence methods which assume a certain population structure and low level of genetic differentiation [7]. To reduce the effect of population structure on the hapFLK analysis, we performed the Principal Component Analysis (PCA [90]) and ADMIX-TURE genetic clustering on all the studied breeds. Prior to the analyses, to reduce effects of linkage-disequilibrium between loci we pruned the dataset using the PLINK function --indep-pairwise 50 10 0.1 using only autosomal genotypes resulting in 95,809 variants for 312 animals. Then we ran ADMIXTURE for K = 1-20 and calculated cross-validation error for each run (-cv). The results of the ADMIXTURE analysis were visualized with PONG software [91] the results of PCA analysis were plotted within the R environment.

Identification of signatures of selection with hapFLK statistics
We performed a genome scan for selective sweeps within each group of breeds and for all breeds simultaneously to infer regions which were under selection during the group/breed formation and in the ancestral population of all studied breeds using a haplotype-based statistics hapFLK [92]. Due to the hapFLK model assuming that selection acts on shared ancestral SNP allele frequencies, we excluded rare SNPs with low minor allele frequencies (MAFs) from each of the breed groups (MAF < 0.05). We also excluded poorly genotyped individuals (< 95% of SNPs with genotypes), loci genotyped in < 99% of samples, SNPs without chromosomal assignments, and SNPs on sex chromosomes in PLINK, using the commands: --maf 0.05, --mind 0.05, --geno 0.01, and --chr 1-26 prior to performing the genome selection scans. This resulted in 506,343 autosomal SNPs for the all-breed analysis, 492,607 and 499,219 for GROUP1 and GROUP2, respectively.
The hapFLK method takes the haplotype structure of the population into account. What was important for our dataset is that this method can account for population bottlenecks and migration. Reynolds distances and a kinship matrix were calculated by the hapFLK program v.1.4 [92]. For the hapFLK analysis, the number of haplotype clusters for each breed group were estimated with fastPhase [93] and were set as -K 40, 25, 35 for the all-breed set, GROUP1, and GROUP2, respectively. The expected maximum number of iterations was set to 30 for three groups. We applied midpoint rooting to all sets of breeds.

P-value calculation
For hapFLK, the calculation of raw p-values was performed assuming that the selected regions represent only a small fraction of the genome [7]. The genome-wide distribution of hapFLK statistics could be modelled relatively well with a normal distribution except for a small fraction of outliers from potentially selected regions [7]. Robust estimations of the mean and variance of the hapFLK statistic were obtained using the R MASS package rlm function to eliminate influence of outlying regions following Biotard and co-workers (2016) [94]. This has been done for each group (all breeds, GROUP1, and GROUP2).
The hapFLK values were Z-transformed using these parameter estimates, and p-values were calculated from the normal distribution in R. The R qvalue package was used to correct p-values for multiple testing [95].

Composite measure of selection (DCMS statistics)
Recent studies demonstrated the high efficiency of composite measures of selection over the single-statistic tests or their simple meta-analysis [96,97]. Composite measures of selection such as de-correlated composite of multiple signals DCMS [96] allow more pricisely locate the selection signal and filter out spurious results specific for some methods. For this study we combined five well-established genome-wide statististics into a single DCMS framework [96]. The DCMS works by combining p-values from different statistics at each locus and correcting for the overall correlation between the statistics based on the covariance matrix. We aggregated the following statistics in the present work: haplotype homozygosity (H1 [98]), modified haplotype homozygosity statistics (H12 [98], fixation index (F ST [99], Tajima's D index [100] and nucleotide diversity (Pi [101]).

Haplotype-based statistics
Autosomal genotypes were phased using SHAPEIT2 software [102] with 400 conditioning states (-states 400) and effective population size parameter equal to 3000 (-effective-size 3000) as a safe estimate of genetic variation within our diverse dataset. The recombination rate along the chromosomes was corrected with a high-resolution ovine genetic map [103].
To calculate the haplotype-based H1 and H12 statistics the phased VCF file was converted to the format required by the H12_H2H1.py script (https://github.com/ ngarud/SelectionHapStats) from Garud and co-workers (2015) [98]. The statistics were calculated for each autosome of each breed using overlapping windows of 25 single nucleotide variants (SNVs, -w 25) with the step size equal to 1 (-j 1) and allowing zero false-positive SNVs per window (-d 0).

Tajima's D statistics
To calculate Tajima's D statistics, we first formed chromosome intervals based on the output of the H1 statistics and then passed them to the bcftools (view) software [104] along with the breed-specific gzipped VCF file, before being piped to the vcftools -TajimaD function. The work was performed in a parallel mode with assistance of GNU PARALLEL [105] to reduce calculation time.

Fixation index (F ST )
To quantify the population differentiation for every SNV we calculated F ST index for each breed against the combined pooled sample of all other breeds using the plink --fst function. Negative F ST values were converted into zeros and the statistic was smoothed for each chromosome using R runmed function in windows of 31 SNPs (k = 31, endrule = "constant") to reduce noise.

Nucleotide diversity (pi)
Nucleotide diversity was calculated using the vcftools -site-pi option for each position and breed separately. To reduce the overall noise the statistic was smoothed with the R runmed function with a window size of 31 SNPs (k = 31, endrule = "constant").

De-correlated composite of multiple signals (DCMS)
At the first step, we combined all the calculated statistics (H12, H1, Pi, Tajima's D, F ST ) into a single spreadsheet based on the SNV name. We then calculated genome-wide rank-based p-values for each statistic (stat_to_pvalue MINOTAUR function) using one-tailed tests (Pi and Tajima Dleft-tailed; H1, H12, and F STright-tailed) of the R MINOTAUR package [106]. To adjust for the correlation among the statistics we calculated the covariance matrix based on 300,000 randomly sampled SNPs using the Cov-NAMcd function with alpha = 0.75. The matrix was then used to calculate the DCMS statistic using DCMS function of the MINOTAUR package. We fitted the resulting DCMS statistics for each breed into a normal distribution using the robust fitting of the linear model method implemented in the rlm R function of the MASS package [94]. The fitted DCMS statistics were then converted into p-values using the pnorm function (lower.tail = FALSE, log.p = FALSE) and the p-values were finally converted to the corresponding q-values using the qvalue R function [95].

Identification of chromosome intervals under selection and candidate genes
We downloaded the ovine gene annotations from the Biomart [107] which correspond to the Oar_v3.1 genome assembly [108]. Next, we considered chromosome intervals with SNPs with adjusted p-values < 0.01 to determine putative regions under selection. The boundaries of each interval were defined by the locations of the first flanking SNPs exhibiting adjusted p-values > 0.1. Within the selected intervals, genes were identified within 1σ value from the most significant SNP based on statistical value (DCMS or hapFLK) distribution similar to Fariello and co-workers [7]. This approach helps to balance the number of candidate genes reported between the "sharp" selection peaks, and intervals with many SNPs exhibiting similar statistics values where larger numbers of genes were reported. Finally, the genes were ranked based on their distance from the SNP with the highest statistics value in each region with larger ranks assigned to more distant genes.