G C A T
T A C G genes
G C A T
Article
Capturing Genetic Diversity and Selection Signatures of the
Endangered Kosovar Balusha Sheep Breed
Olusegun O. Adeniyi 1,† , Rebecca Simon 1,† , Hysen Bytyqi 2, Waltraud Kugler 3, Hajrip Mehmeti 2,
Kaltrina Berisha 2 , Mojca Simčič 4, Mohamed Magdy 1 and Gesine Lühken 1,*
1 Institute of Animal Breeding and Genetics, Justus Liebig University, 35390 Giessen, Germany;
olusegun.o.adeniyi@agrar.uni-giessen.de (O.O.A.); rebecca.simon@agrar.uni-giessen.de (R.S.);
mohamed.ismail-magdy-saadeldin-sabri@agrar.uni-giessen.de (M.M.)
2 Faculty of Agriculture and Veterinary, University of Pristina, 10000 Pristina, Kosovo;
hysen.bytyqi@uni-pr.edu (H.B.); hajrip.mehmeti@uni-pr.edu (H.M.); kaltrina.berisha@uni-pr.edu (K.B.)
3 SAVE Foundation, 9000 St. Gallen, Switzerland; waltraud.kugler@save-foundation.net
4 Department of Animal Science, Biotechnical Faculty, University of Ljubljana, SI-1000 Ljubljana, Slovenia;
mojca.simcic@bf.uni-lj.si
* Correspondence: gesine.luehken@agrar.uni-giessen.de
† These authors contributed equally to this work.
Abstract: There is a growing concern about the loss of animal genetic resources. The aim of this study
was to analyze the genetic diversity and potential peculiarity of the endangered Kosovar sheep breed
Balusha. For this purpose, a dataset consisting of medium-density SNP chip genotypes (39,879 SNPs)
from 45 Balusha sheep was generated and compared with SNP chip genotypes from 29 individuals
of a second Kosovar breed, Bardhoka. Publicly available SNP genotypes from 39 individuals of
the relatively closely located sheep breeds Istrian Pramenka and Ruda were additionally included
in the analyses. Analysis of heterozygosity, allelic richness and effective population size was used
Citation: Adeniyi, O.O.; Simon, R.; to assess the genetic diversity. Inbreeding was evaluated using two different methods (FIS, FROH).
Bytyqi, H.; Kugler, W.; Mehmeti, H.; The standardized FST (di) and cross-population extended haplotype homozygosity (XPEHH) methods
Berisha, K.; Simčič, M.; Magdy, M.;
were used to detect signatures of selection. We observed the lowest heterozygosity (HO = 0.351) andLühken, G. Capturing Genetic
effective population size (Ne5 = 25, Ne50 = 228) for the Balusha breed. The mean allelic richness levelsDiversity and Selection Signatures of
the Endangered Kosovar Balusha (1.780–1.876) across all analyzed breeds were similar and also comparable with those in worldwide
Sheep Breed. Genes 2022, 13, 866. breeds. FROH estimates (0.023–0.077) were highest for the Balusha population, although evidence of
https://doi.org/10.3390/ decreased inbreeding was observed in FIS results for the Balusha breed. Two Gene Ontology (GO)
genes13050866 TERMs were strongly enriched for Balusha, and involved genes belonging to the melanogenesis
and T cell receptor signaling pathways, respectively. This could result from selection for the special
Academic Editor: Giovanni Amori
coat color pattern of Balusha (black head) and resistance to certain infectious diseases. The analyzed
Received: 24 March 2022 diversity parameters highlight the urgency to preserve the local Kosovar Balusha sheep as it is clearly
Accepted: 9 May 2022 distinguished from other sheep of Southeastern Europe, has the lowest diversity level and may
Published: 12 May 2022 harbor valuable genetic variants, e.g., for resistance to infectious diseases.
Publisher’s Note: MDPI stays neutral
with regard to jurisdictional claims in Keywords: ovis aries; conservation/preservation; population structure; ROH; GO TERMS; Pramenka;
published maps and institutional affil- Balkan
iations.
1. Introduction
Copyright: © 2022 by the authors. Since the domestication of sheep at around 10.000 BC [1], many breeds with high di-
Licensee MDPI, Basel, Switzerland. versity in appearance, reproductive traits and performance characteristics have developed
This article is an open access article worldwide through breeding and selection. The Balkans are considered to be the main entry
distributed under the terms and point in the immigration of agriculture in general, and sheep especially, to present-day
conditions of the Creative Commons
Europe during the Neolithic period [2–5]. There is evidence that in the 20th century, cross-
Attribution (CC BY) license (https://
breeding of Merinoland sheep with some of the local Pramenka breeds took place in the
creativecommons.org/licenses/by/
4.0/). Balkans [6]. However, due to changes in agriculture and the loss of economic importance of
Genes 2022, 13, 866. https://doi.org/10.3390/genes13050866 https://www.mdpi.com/journal/genes
Genes 2022, 13, x FOR PEER REVIEW 2 of 19 
 
Genes 2022, 13, 866 century, crossbreeding of Merinoland sheep with some of the local Pramenka b2roefe1d8 s took 
place in the Balkans [6]. However, due to changes in agriculture and the loss of economic 
importance of sheep in Europe [7], numbers of flocks are decreasing and many traditional 
bshreeeepdsin tEhuatr oapree[ 7n],ont uumsbeder sino fcfloomckms aerrecidael crheuassbinagnadnrdy maaren ythtrraedaitteionneadl  bwreietdhs ethxatitnacrteion [8] 
(nborteuesde dexinamcopmlemse r[c9i]a)l. hTuhsibsa nisd raylsaore tthhere acatesnee dinw Kithoseoxvtion.c tFioonur[8 ]P(rbarmeeednekxaa mbrpeleesd[s9 ]()r.ough-
wThoiosliesda ltsyoptehse),c washeiicnhK aoreso mvoa.inFloyu ur sPerdam foern ktraipblree epdusr(prousegsh, -mwiolokl–emd etyapt–ews),owolh [i1ch0]a, raere con-
smidaeinrleydu tsoe dbfeo lrotcriapll einp Kuropsoosveso,:m Sihlka–rmri,e Kato–wsoovoal ,[ 1B0a],rdarheockoan s(iBdAerRed, atolsboe kloncoawl in Kaso sWovhoi:te Me-
tSohhairarni, KProasmovean, kBaa)r dahnodk Baa(BluAsRh,aa (lBsoAkLn)o. wBAnRas aWndh iBteAML esthoeheiapn (PFrigamureen 1k)a )aaren dkeBpatl uinsh aareas in 
t(hBeA sLo).uBthAwReasntedrnB ApLarsth oefe tph(eF cigouurnetr1y)  a(rDeukkeapgt jiinnia prelaasnein), twhehesoreu tthhweye swteernrep marationflyth seelected 
fcooru tnrtaryits(D ausskoacgijaintei dpl waniet)h, wmhilekre ptrhoedyuwcetiroenm [a1i1n]l.y selected for traits associated with milk
production [11].
 
FFiigguurree 11.. GGeeooggrraapphhicicd disitsrtibriubtuiotnioonf ofof ufor uanr aalnyzaelydzsehde espheberpee bdrseiendsso iunth seoausttehrenaEstuerronp Ee.uTrohepeb.l aTchke-  black-
hheeaaddeeddB Baalulushsaha(B (ABLA),Lt)h, ethweh iwtehBitaer dBhaorkdah(oBkAaR ()B, oAriRg)in, aotrinigginfraotmintgh efrsoamm etharee asainmwe easrteeran inK owsoevsote, rn Ko-
sthoevoA,l bthane iaAnlbRaundiaan(R RUuDd; ap h(RotUo:DA; .pHhoodtao): aAn.d HthoedIas)tr aianndP trhame eIsntkraia(nIS PT;rapmhoetno:kVa. (RISezTa; rp),hthoetola: tVte.r Rezar), 
tohnee lsaatmtepr loednein saSmlovpelneida .inP lSelaosveennoitae.t Phaletathsee nporotpe otrhtaiotn tshean pdrodpisotarnticoensso af nthde dmisatpaanrceens ootfd tihsep lmayaepd are not 
dreiasplislatiycaeldly r, ethaeliystaircealolnyl,y thfoeryo arrieen otantliyo nf.or orientation. 
The medium-size Balusha is coarse-wooled with a black head and neck while the
fleeceThane dmleegdsiuamre-swizhei tBe.alOusnhlya risa mcosaarrsee-hwoornoeleddw witihths car ebwla-schka pheedadh oarnnds ,nwechke rweahsile the 
felweeecsea arnedp olellgesd .arIen wcohnittrea. sOt,nthlye rmamedsi uamre- shizoerdnecdo awrsieth-w socorelewd-Bshaardpheodk haoirsncso,m wphleetreelays ewes 
awrhe itpeo. lTlehde. mIna jocorintytroafstr,a tmhse hmaveedisupmira-sl ihzoerdn scobaurtsseo-mweooalreedp oBlalerdd,hwokhaic his iscocmonpslideteerleyd white. 
Tuhnety mpiacajol.ritByo othf rbarmeesd hs aavre sapdiarpatle hdofronrs mbuilkt spormodeu acrteio pnobllaesde,d wohnicghr aizsi ncognosfidnearteudra ul ntypi-
cpaals. tuBroetsh[ 1b1r,e1e2d].sM areed iaudma-psitzeedd fRoUr Dmilska pbrroededucwtiiothn cbleaasredM eornin gorianztirnog roesf snioant,uwrahl icphastures 
[w1a1s,1a2l]r.e aMdeydciounmfir-msizededb yRUgeDne itsi caa bnraeleydsi sw[2it]h.  RclUeDar iMs uesreidnofo irntriopglerepsusripoons,e ws,hwicithh wahsit aelready 
cmoendfiiurme-qdu balyit ygeflneeectiec, laengaslaynsdish [e2a]d. ,RbUutDr airse ulys,eindd fiovird turaiplslec apnubrepcoosmesp, lweteitlhy bwlahcikt.e Tmheedium-
qraumalsitcya rfrlyeebcige,s lpeigrasl ahnodrn hs,ewadh,i lebufetm raalreeslya,r einpdoilvleiddu[1a2l,s1 3c]a.nT hbee IcSoTmispwleidteelsyp rbelaadckin. Tthhee rams 
cwaersrtye rbnigre sgpioirnalo hf othrensB,a wlkhanileP feenminasulelsa .aTreh eploalrlgeed I[S1T2,,w13h]i.c Thhies uISsTed isf owr imdeilskpprreoaddu icnti tohne,  west-
has a white coarse fleece with a different amount of black spots, to the extent that some
ern region of the Balkan Peninsula. The large IST, which is used for milk production, has 
animals are nearly black. The rams are always horned, while ewes can be both polled or
ah owrnheidte. IcSoTarasree kflneoewcen wfoirthth ae dhiigffhe-rfeant ct oanmteonutnint tohfe birlamciklk s[p6o,1t2s], .to the extent that some ani-
malsN aorwe andeaayrsly, e bsplaecckia. llTyhteh erabmlacsk a-hreea adlewdaByAs Lhsohreneepdh, awvehibleec eowmeesr acraen.  Ibt ies basostuhm peodlled or 
hthoartn7e0d%. IoSfTt haerew khnoolewpno fpourl athtieo nh,iwghh-ifcaht ccoonnstiesntst oinf  otnhleyir3 m00ilinkd [i6v,1id2u].a ls [14], is kept by
only tNwoowbaredeadyesr,s .eIsnpaedcidailtliyo nt,hteh ebrleaicskt-heepadosesdib BiliAtyLo sfhueneppla hnanveed bcreocsosm-bree eradrine.g Iot fisB AasLsumed 
twhiatht 7B0A%R ,osfi nthcee bwothhoblere pedospaurleatrieoanre, dwthogicehth ceorninsitshtes soafm oenalyre 3a0(0D uinkdaigvjiindiuraeglsi o[n1)4u],n ids ekrept by 
oanleyx ttewnosi bvereseydsteerms. oInf  faadrmdiitnigon[1, 1th].eTrhee irse tfhoree p, tohsesibbrieleidtyi sofa tuanvpelarynnhiegdh crriosksso-fberxeteindcitniogn of BAL 
wanidthl BosAsRo,f sdinivceer bsiotyth.  bPrreeseedrsv ianrge roeraigrienda ltobgreeethdesra ilns othme esaanmsep areresear (vDinugkcahgajirnaic rteergiisotincs) under 
an extensive system of farming [11]. Therefore, the breed is at a very high risk of extinction 
 
Genes 2022, 13, 866 3 of 18
and traits that may become important with changing requirements, e.g., by the farmers or
environmental conditions, offering future breeding options [8,15,16]. The importance of
regional breeds as a cultural heritage, which should be preserved, is also a factor in favor
of preserving old breeds e.g., [7,17,18].
The Illumina Ovine 50K SNP chip [19] is a powerful tool to analyze the genetic
diversity within a breed and across breeds and to determine the status quo of populations.
This was confirmed in many studies on sheep diversity e.g., [20–24]. Various approaches
are available to determine the degree of inbreeding and thus to assess the diversity of
a breed as a whole and of individual animals. Classically, the expected and observed
heterozygosity, the effective population size (Ne) and the Wright’s F statistical indices (Fis)
can be used [25]. Moreover, with genome-wide SNP data, it is possible to quantify the extent
of inbreeding, determined by the runs of homozygosity (ROH) [26], even with missing
pedigree information [27]. The different lengths of the segments containing homozygous
loci can be used to define whether the inbreeding occurred recently or can be rated as an
ancient event [28].
Regions with high ROH are regions that also carry signatures of selection [29]. Differ-
ent methods are available to detect those signatures of selection [30,31] that indicate patterns
of reduced diversity close to genes that have been under strong selection, artificially or
naturally caused, in a certain population [32]. The di statistics is an established method used
to calculate the locus-specific divergence of a specific breed against other breeds [30]. The
cross-population extended haplotype homozygosity (XPEHH) is a haplotype-related analy-
sis that detects differences between two populations in which a selected allele is fixed or
approaching fixation in one population while still polymorphic in the other population [31].
The diversity of some sheep breeds from Eastern Europe was studied and eval-
uated, mainly by using microsatellite markers, but the BAL breed has not yet been
considered [2,33,34]. With our study, we aimed to determine the genetic diversity and
population structure of the Kosovar sheep breed BAL for the first time. Special attention is
paid to the genetic differentiation of this breed from the BAR breed kept in the same region.
Additionally, published Ovine 50K SNP chip data from the two relatively geographically
close sheep breeds Istrian Pramenka (IST) and Ruda (RUD) (Figure 1) were included for
comparison. ROH were analyzed and signatures of selection based on di statistics as well as
on XPEHH were compared. The results should provide a valid basis for the justification and
establishment of conservation programs for the endangered BAL breed in a further step.
2. Materials and Methods
2.1. Animal Samples and Data
For this study, blood samples of BAL (n = 45) and BAR (n = 29) sheep breeds were
collected from two farms in Kosovo (Figure 1). All animals were females, except for
a single male in the BAL sample set. Unrelated animals were sampled where possible
to avoid genetic relation. DNA was extracted due to manufacturer’s instructions from
the blood samples using the NucleoSpin Blood Kit (Macherey-Nagel, Düren, Germany)
and genotyped with the Ovine 50K SNP BeadChip (Illumina, San Diego, CA, USA). In
addition, published Illumina Ovine 50K SNP BeadChip data from the sheep breeds Istrian
(IST, n = 23) and Ruda (RUD, n = 16) were available from a recent study on Balkan sheep
breeds (Figure 1) [2]. These data were merged and a total of 39,879 autosomal SNPs were
available for the analysis. Data were filtered to exclude loci with minor allele frequency
(<5%), SNPs with a low call rate (<95%) and animals with more than 5% of missing
genotypes. To overcome the effect of closely related animals, the data were checked for
relatedness which was estimated as a proportion of identity-by-descent (IBD) between
animal pairs using –related option in KING v2.2.4 [35] and one of each animal pair with an
IBD proportion > 0.35 was removed. After quality control (QC) procedures using PLINK
v1.90 [36], a total of 39,214 autosomal SNPs and 92 samples (Table S1) remained in the
dataset for genomic-based diversity analysis.
Genes 2022, 13, 866 4 of 18
2.2. Analysis of Population Structure and Genetic Diversity
The above dataset was linkage disequilibrium (LD)-pruned using the –indep-pairwise
function in PLINK v1.90 (this version of PLINK was also used for other analyses with that
program; © 2022 Shaun Purcell, GNU General Public License). One of a pair of SNPs in
high LD (r2 > 0.5) within a window size of 50 SNPs and with a window slide of 5 SNPs
was removed. Pairwise identity-by-state (IBS) distances and principal component analysis
(PCA)-based multidimensional scaling analysis were performed in PLINK, and a plot show-
ing the first two principal components was constructed with the package “ggplot2” [37] in
R [38]. The observed heterozygosity (HO) and expected heterozygosity (HE) were calculated
using HIERFSTAT R package [39]. Additionally, rarified allelic richness and private allelic
richness were estimated with the ADZE v1.0 software [40] using a standardized sample size
for each breed. The Wright inbreeding coefficient was estimated as F HE−HOIS = H , where HE O
and HE are observed and expected heterozygosities, respectively [41]. Effective population
size (Ne) for all breeds was calculated based on LD using the SNeP package [42]. Default
parameters were adjusted for sample size, recombination rate according to [43], and occur-
rence of mutation as suggested by [44]. The Ne five (Ne5) and fifty (Ne50) generations ago
were estimated. The ADMIXTURE program (Version 1.3.0; [45] was employed to determine
the ancestry and results were visualized using “pophelper” R package version 2.3.0 [46].
Clustering was calculated for K values from 2 to 4 and optimal number of clusters was
determined by adding the –cv flag [45].
2.3. Estimation of Runs of Homozygosity
Runs of homozygosity (ROH) were estimated using the QC data and –homozyg function
in PLINK. One missing SNP, one heterozygous SNP, and ROH longer than 1Mb were
allowed per run. A scanning window threshold of 0.05, a sliding window of 50 SNPs,
a maximum gap between consecutive SNPs set to 250 kb, and a minimum average SNP
density of more than one SNP/100 kb were used. The minimum number of SNPs (l) per
log αe
run was calculated as suggested by [43]. l nsn= ilog 1−het , where α is the percentage ofe( )
false-positive ROH (set to 0.05 in this study), ns is the number of SNPs per individual, ni
is the number of individuals, het is the mean heterozygosity across all SNPs. Calculated l
was 39 for BAL, 36 for BAR and 34 for IST and RUD. The total and mean ROH length were
estimated for each individual and breed. In addition, the total and mean ROH numbers for
each individual and breed were determined. To compare ROH frequency between breeds,
ROH segments were categorized into 3 length classes (1–5, >5–10 and >10 Mb).
The inbreeding coefficient of ROH (FROH) was determined according to [47]:
F = ∑ LROHROH L , where LAUTO ROH is the sum of the length of all ROH per animal and LAUTO is
the total autosomal SNP coverage (2.648 Gb). The average FROH per breed was determined
for LROH above 1, 5 and 10 Mb, which is indicative of inbreeding up to 50, 10 and 5 gen-
erations ago, respectively [48–50]. Furthermore, box plots were constructed to show the
within-breed distribution of FROH per breed for LROH above 1 Mb.
2.4. Detection of Signatures of Selection
To determine selection signatures of BAL, the FST and the cross-population extended
haplotype-based homozygosity score test (XPEHH) [31] was employed using the QC data.
Pairwise FST [51] values of each SNP were estimated for BAL against a combination of
all other sheep breeds in this study using the VCFtools v0.1.15 software. After this, the
FST values obta[ined] were sta[ndar]dized (d[i) ac]cording to the method published by [30]:Fij −E Fij
d = ∑ ST Si j 6 i [ ]T= ij where E FijST and sd FijST denote the expected value and standardsd FST
deviation of FST between breeds i and j calculated from all 39,214 SNPs. The locus-specific
divergence calculated using the di statistics was averaged for SNPs in non-overlapping
windows of 1 Mb. A total of 2562 informative windows with an average of 16.07 SNPs
per window (SD = 3.6) were left after removal of windows with less than five SNPs. We
Genes 2022, 13, 866 5 of 18
plotted the window-averaged di values across the chromosome, and the top 1% windows
of the empirical distribution were defined as putative selection regions. On the other hand,
pairwise comparison between BAL and each of the other breeds for XPEHH scores were
computed separately with the R package rehh [52] using haplotype information. Haplotypes
were estimated with fastphase 1.4 [53] by using population label information to estimate
phased haplotype background and applying the following options for each chromosome:
-T10 -C25 -Ku40 -Kl10 -Ki5 -Km1000 -H100, where T is the number of the start of EM
algorithm, C is the number of iterations of EM algorithm, Km is the number of SNP loci
used for cross-validation, Ki is the interval between values of number of clusters, Ku
and Ki are the upper and lower limit of number of clusters, respectively. For compar-
ison of selection signatures and visualization, |XPEHH| scores were transformed into
pXPEHH = −log[1−2|φXPEHH− 0.5|] in which φ XPEHH denotes the Gaussian cumula-
tive distribution function, which was calculated by the pnorm R function. Since the XPEHH
scores were approximately normally distributed, pXPEHH can be interpreted as log10(1/P).
Therefore, pXPEHH was the two-sided p-value for testing the null hypothesis that no selec-
tion occurred [48], and we defined putative selection signals as pXPEHH ≥ 5 equivalent to
a p-value of 0.00001. Common candidate selection regions from all three XPEHH analyses
with a minimum of three SNPs were chosen as selection signatures specific to BAL. The
selected regions were defined with a maximum base-pair length of approximately 1 Mb.
2.5. Gene Annotation and Enrichment Analysis
Genes located within the candidate selection regions for BAL were identified using
the UCSC genome browser [54] with the selection of the Oar_v4.0 assembly. Furthermore,
functional enrichment analysis was performed on genes identified by both FST (di) and
XPEHH methods separately using the Database for Annotation, Visualization and Inte-
grated Discovery (DAVID) software (https://david.ncifcrf.gov, accessed on 14 February
2022) [55,56]. The gene list was analyzed by selecting the human gene annotations with
their background due to low number of recognized genes when sheep gene annotations
were selected. The gene ontology annotation category, including biological process (BP) and
molecular function (MF) terms, was investigated in this study. The GO TERMS of enriched
genes with a p-value threshold of ≤0.05 were considered significant, and correction for
multiple testing was carried out using the FDR procedure [57]. Overlapping regions of
signatures of selection were defined as those regions above the threshold for both methods
and within the same chromosomal location [58].
3. Results
3.1. Population Structure
Results of genetic diversity indices for the Balkan breeds in this study are given in
Table 1. The observed and expected heterozygosity ranged from 0.351 (BAL) to 0.384 (IST)
and 0.341 (BAL) to 0.402 (RUD), respectively. Allelic richness was comparable among the
sheep breeds with the highest observed in the RUD (1.876) and the lowest in the BAL
(1.780). Similarly, the highest private allelic richness was found in RUD (0.012), and the
lowest in BAL (0.009). The FIS for the Balusha (−0.029) and Istrian breeds (−0.009) were
observed to be negative with the highest value observed for Ruda (0.054) and the lowest
for Istrian. The recent effective population sizes (Ne5) were estimated for all breeds and
ranged from 25 (BAL) to 35 (BAR), while the effective population sizes 50 generations ago
(Ne50) varied between 228 in BAL and 324 in BAR. The PCA plot (Figure 2) showed that
the first component explained 55.7% of the genetic variability between the breeds. It clearly
distinguished the BAL breed from the others. In addition, the second principal component
accounted for 30.4% of the variability and clearly separated the other breeds from each
other (BAR, IST and RUD). The ADMIXTURE analysis revealed a clear pattern of ancestry
for the breeds (Figure 3). Based on cross-validation error, an optimal value of K = 3 clusters
was identified. At K = 2, the other breeds were clearly separated from the BAL breed.
Genes 2022, 13, x FOR PEER REVIEW 6 of 19 
 
ranged from 25 (BAL) to 35 (BAR), while the effective population sizes 50 generations ago 
(Ne50) varied between 228 in BAL and 324 in BAR. The PCA plot (Figure 2) showed that 
the first component explained 55.7% of the genetic variability between the breeds. It 
clearly distinguished the BAL breed from the others. In addition, the second principal 
component accounted for 30.4% of the variability and clearly separated the other breeds 
from each other (BAR, IST and RUD). The ADMIXTURE analysis revealed a clear pattern 
of ancestry for the breeds (Figure 3). Based on cross-validation error, an optimal value of 
Genes 2022, 13, 866 K = 3 clusters was identified. At K = 2, the other breeds were clearly s6eopf 1a8rated from the 
BAL breed. However, a genetic distinction between the IST breed and the other two 
breeds (BAR and RUD) was observed at K = 3. 
However, a genetic distinction between the IST breed and the other two breeds (BAR and
TRaUbDle)  w1a. sGoebnseetrivce d iavteKrs=ity3. (LD = r2 > 0.5) and effective population size of breeds analyzed in this 
study. 
Table 1. Genetic diversity (LD = r2 > 0.5) and effective population size of breeds analyzed in
this stBudrey.ed n Ho He AR pAR FIS Ne5 Ne50 
BaBlrueesdha 3n0 0H.3o51 He0.341 AR 1.7p8A0R F0IS.009 Ne5−0.02N9e 50 25 228 
BaBradluhshoaka 2306 0.3.35172 0.3401.3781 .780 1.80.40069 −00.0.20910 25 0.0182 28 35 324 
BIasrtdrhioakna 2261 00.3.37284 0.3708.3801 .846 1.80.40140 0.001.8011 35 −0.00392 4 31 277 
Istrian 21 0.384 0.380 1.844 0.011 −0.009 31 277
Rudaa 15  0.3.38080 0.4002.4021 .876 1.80.70162 0.05.4012 29 0.0542 66 29 266 
 
Figure 2. Principal component analysis plot for the first (PC1) and second (PC2) components of
Figure 2. Principal component analysis plot for the first (PC1) and second (PC2) components of 
Balusha, Bardhoka, Istrian and Ruda sheep breeds. The different breeds are color-coded as mentioned
Binatlhueslhegae, nBda. rNdohteotkhaa,t tIhsetrbialune Baanldus hRaucdluas tsehr iesecple abrrlyeesedpsa. rTathede fdroimffetrheenret sbt oref ethdesb arereed sc.olor-coded as men-
tioned in the legend. Note that the blue Balusha cluster is clearly separated from the rest of the 
breeds. 
 
Genes 20G2en2e,s1 230,228,6 163, x FOR PEER REVIEW 7 of 19 7 of 18
 
 
FFiigguurer e3. 3A.dAmdixmturiex tpulorte foprl aoltl sfhoerepa lblrseehdese (pBablursehead, Bsa(rBdahlouksah, Iast,rBiaanr adnhdo Rkuad,aI) satnrailaynzeadn ind thReu da) analyzed in the
sttuuddyy. .AA clecalre apratptearntt eofr nanocefsatrnyc ies svtirsyibiles fvoirs tihbel ebrfeoerdst.h Pelebarse endoste. tPhlaet aKs e= n3 owtaes tihdaenttKifie=d 3asw as identified as the
tohpe toipmtiamlanl nuummbbeerr oof fclculsutesrtse. rs.
33..22. .RRunusn osf HofomHoozymgooszityyg osity
A total of 2562 ROH across all breeds were identified. Out of the 92 animals analyzed, 
90 hadA att loetaastl oonfe2 R5O6H2 .R TOheH nuamcrboesr sofa RllObHr eveadrisedw weirtheinid aenndt bifietewde.enO burteeodfs,t hweith9 2a animals analyzed,
h9i0ghear dfreaqtuleenacsyt oof nsheoRrt OROHH.  Tsehgemnenutms inb earll obfreRedOs H(Fivguarrei e4d). Twhiet hmienana nudmbbetrw ofe en breeds, with a
RhOigHh beertwfreeeqn ubreenedcsy raonfgsehd ofrrotmR 1O2.H67 (sISeTg)m toe 4n2t.9s0 i(nBAaLll). bLrikeeewdisse(, Fthieg suhroerte4s)t. mTehane mean number of
lRenOgHth obf eRtOwHe ewnasb orbeseedrvsedra inn gISeTd (5f8r.o05m M1b2),. a6n7d( tIhSeT lo)ntgoe4st2 i.n9 B0A(LB A(20L4).5. 1L Mikbe)w (Tiasbel,e the shortest mean
2). The frequencies of ROH segments in the different length classes (1–5, >5–10, >10 Mb) 
wleenreg tchomopf aRraObHle bwetawseoenb sberrevedesd. TinheI SfrTeq(u5e8n.0ci5esM ofb )R,OaHn dint hthee loshnogrteesstt ilennBgtAh Lcla(2ss0 4.51 Mb) (Table 2).
rTahngeefdr efrqoume 6n5c.6ie2%s  o(BfARRO) Hto 7s2e.5g6m%e (nIStTs),i nwitth evadluifefse oref n5.t21l%en (gBtAhLc) laansds 7e.s35(%1– (5R,U>D5)– 10, >10 Mb) were
ocbosmervpeadr ainb tlheeb lentwgthe ecnlasbs r>e1e0d Ms.bT (Fhiegufree q5)u. Tehnec imeseaonf FROH H(1 Minbt) hvearsiehdo bretetwsteelne ngth class ranged
bfrreoemds,6 w5.i6th2 %the( hBiAghRes)t tfoou7n2d.5 in6 %the( IBSATL) ,bwreietdh (0v.a07lu7)e, sanodf t5h.e2 1lo%we(sBt Ain Lth)ea InSdT b7r.e3e5d% (RUD) observed
(i0n.0t2h3)e. Tlehne gwthithcinla-bsrsee>d1 0disMtribbu(tFioign uorf eFR5O)H. (T1h Mebm) ies avnisuFalized (i1n ROH M
Figbu)rve a5r, iweidthb tehtew een breeds, with
highest median observed in the BAL breed. In addition, mean FROH at a minimum ROH 
ltehnegthh iogf h5 eMsbt afnodu 1n0d Mibn wtehree cBaAlcuLlabterde, eadnd( t0h.e0 h7i7g)h,eastn FdROHth vealuloesw oef s0.t03in9 (t5h MebI)S aTndb reed (0.023). The
0w.0i1t1h (i1n0- bMrbe)e wdedrei sfotruinbdu itni oBnALo f(TFaRbOleH 2)(. 1 Mb) is visualized in Figure 5, with the highest median
Genes 2022, 13, x FOR PEER REVIEW 
 observed in the BAL breed. In addition, mean FROH at a minimu
8 mof 1R9 OH length of 5 Mb
and 10 Mb were calculated, and the highest FROH values of 0.039 (5 Mb) and 0.011 (10 Mb)
were found in BAL (Table 2).
 
Figure 4. Frequency distribution of number of ROH in  varying class lengths and for each breed.
Figure 4. Frequency distribution of number of ROH in varying class lengths and for each breed. 
Table 2. Descriptive statistics for ROH and FROH of breeds analyzed in this study. 
 Breed 
 Balusha Bardhoka Istrian Ruda 
Number of samples 30 26 21 15 
ROH length (Mb)     
Mean (SD) 204.51 (88.30) 124.85 (125.81) 58.05 (48.49) 113.16 (156.06) 
Range 64.28–361.00 0–516.50 10.77–163.54 0–442.61 
ROH Number     
Mean (S.D) 42.90 (16.31) 25.73 (24.09) 12.67 (9.19) 22.67 (30.35) 
Range 15–73 0–101 2–35 0–88 
FROH (1 Mb)     
Mean 0.077 0.046 0.023 0.042 
Range 0.02–0.14 0–0.20 0–0.06 0–0.17 
F     ROH (5 Mb) 
Mean 0.039 0.026 0.011 0.024 
Range 0–0.08 0–0.12 0–0.03 0–0.10 
F     ROH (10 Mb) 
Mean 0.011 0.007 0.003 0.009 
Range 0–0.03 0–0.05 0–0.02 0–0.04 
 
Genes 2022, 13, 866 8 of 18
Table 2. Descriptive statistics for ROH and FROH of breeds analyzed in this study.
Breed
Balusha Bardhoka Istrian Ruda
Number of samples 30 26 21 15
ROH length (Mb)
Mean (SD) 204.51 (88.30) 124.85 (125.81) 58.05 (48.49) 113.16 (156.06)
Range 64.28–361.00 0–516.50 10.77–163.54 0–442.61
ROH Number
Mean (S.D) 42.90 (16.31) 25.73 (24.09) 12.67 (9.19) 22.67 (30.35)
Range 15–73 0–101 2–35 0–88
FROH (1 Mb)
Mean 0.077 0.046 0.023 0.042
Range 0.02–0.14 0–0.20 0–0.06 0–0.17
FROH (5 Mb)
Mean 0.039 0.026 0.011 0.024
Genes 2022, 13, x FOR PEER REVIEW Range 0–0.08 0–0.12 0–0.03 0–0.10 9 of 19 
 FROH (10 Mb)
Mean 0.011 0.007 0.003 0.009
Range 0–0.03 0–0.05 0–0.02 0–0.04
 
FiFgiugruer e55. .BBooxx pplloottss ooff wwiitthhiinn-b-brereededd idstirsitbruibtiuotnioonf roufn rsuonfsh oofm hoozmygoozsyitgyoisnibtyre eindbinrgeecdoienffigc cieonetfffoicrient for 
eaecahc hshsheeeepp bbrreeeedd.. 
3.3. Signatures of Selection Based on Locus-Specific Divergence (FST)
3.3. SigTnhaetudresst aotfi sSteiclescwtioenre Bcaaslecdu loante Ldofcours-BSApLecaifgica iDnsivt earlgl eontchee r(Fi bSrTe)e ds combined in this
studTyh.eT dhie sdtaisttirsitbiucsti owneoref  dciaslctautliastiecds  wfoinr dBoAwLv agluaeisnsist dailsl polathyedr binreFeidgus rceo6m. bAintoetda lin this 
study. The distribution of di statistics window values is displayed in Figure 6. A total of 
26 informative windows (top 1% of the total number of informative windows) were iden-
tified as putative selective sweeps using the di statistics, with significant regions observed 
on sheep chromosomes 2, 4, 6, 10, 13, 19, 21, 23 and 24 (Table S2a). The informative win-
dow with the highest di statistics value (2.97) was found on chromosome 2. Based on the 
di statistics, 194 genes were identified in the regions under selection (Table S3). Regarding 
enrichment analysis of genes, 177 of the genes identified from selection signals using the 
di statistics were recognized by the DAVID software (Table S2b). Furthermore, the result 
shows that 27 GO TERMS were enriched (p-value < 0.05) (Table 3). 
 
 
Genes 2022, 13, x FOR PEER REVIEW 9 of 19 
 
 
Figure 5. Box plots of within-breed distribution of runs of homozygosity inbreeding coefficient for 
each sheep breed. 
3.3. Signatures of Selection Based on Locus-Specific Divergence (FST) 
Genes 2022, 13, 866 9 of 18
The di statistics were calculated for BAL against all other breeds combined in this 
study. The distribution of di statistics window values is displayed in Figure 6. A total of 
26 informative windows (top 1% of the total number of informative windows) were iden-
toiffie2d6 iansf pourmtaatitvivee sewleinctdivoew ssw(etoepps1 u%sinofg tthhee dtoit satlantiustmicbs,e wr iotfh isnifgonrimficaatnivte rewgiinondso owbss)ewrveerde 
oidne snhteifieepd charsopmuotsaotimveess 2el, e4c,t i6v, e10s,w 1e3e, p1s9, u2s1i,n 2g3 tahnedd 2i4s (tTaatibslteic Ss,2aw).i tThhsei ginnfiofircmanattivreeg wioinn-s
doobwse rwveitdho tnhes hheiegphecshtr odmi sotsaotimsteicss2 v, 4a,lu6,e1 (02,.9173), 1w9a, 2s 1f,o2u3nadn odn2 4ch(TroambleosSo2ma)e.  T2h. eBainsfeodr mona ttihvee 
dwii sntdatoiwsticws,i t1h94th geenheigs hweesrted iidsetnattiifsiteidcs inv athlue ere(g2i.9o7n)s wunads efro usenldecotinonc h(Troabmleo sSo3m). eRe2g. aBrdasinegd 
eonnritchhemdeinstt aantiastlyicssi,s 1o9f4 gegnenese,s 1w77e roef tidhee ngteifineeds iidnenthtiefieredg firoonms  usenldecetriosnel esicgtinoanls( uTasibnlge tSh3e) .
dRie sgtaartdisitnicgs ewnreirceh rmeceongtnainzaeldy sbisy othf eg eDnAesV, I1D77 sooffttwheargee (nTeasbildee Sn2tbifi).e Fdufrrtohmersmeloercet,i othnes irgensuallst 
using the di statistics were recognized by the DAVID software (Table S2b). Furthermore,
shows that 27 GO TERMS were enriched (p-value < 0.05) (Table 3). 
the result shows that 27 GO TERMS were enriched (p-value < 0.05) (Table 3).
 
Figure 6. Genome-wide distribution of di statistic for all 1 Mb windows across all autosomes. The di
values calculated from comparison of BAL population against other Balkan breeds (BAR, IST and
RUD) in this study. The dashed gray line indicates the top one percent of total informative windows.
 Note that on chromosomes 2, 4, 6, 10, 13, 19, 21, and 23 signals exceed the threshold.
Table 3. Enriched GO TERMS determined by DAVID from genes identified in selection regions
specific to the BAL breed using FST (di). Please note that GOTERM_MF refers to molecular function
and GOTERM_BP refers to biological process.
Category Term Genes p-Value FDR
GO:0016616: oxidoreductase activity, acting on the
GOTERM_MF_DIRECT CH-OH group of donors, NAD or NADP LDHC, UEVLD, 0.000 0.095
as acceptor LDHA, MDH1B
GOTERM_BP_DIRECT GO:0060789: hair follicle placode formation CDC42, GNAS, CTNNB1 0.001 0.425
GOTERM_BP_DIRECT GO:0030182: neuron differentiation EDN3, FZD5, SNPH,FZD7, ID3, WNT4 0.001 0.425
KLF7, NRP2, CREB1,
GOTERM_BP_DIRECT GO:0007411: axon guidance DCC, EPHA8, 0.002 0.425
EPHB2, FEZF2
GOTERM_BP_DIRECT GO:0060231: mesenchymal to epithelial transition TCF15, FZD7, WNT4 0.002 0.425
CDC42, TUBB1, GNAS,
GOTERM_MF_DIRECT GO:0005525: GTP binding GIMAP2, SUCLG2,GIMAP4, GIMAP5, 0.004 0.602
PCK1, GIMAP7, RAB22A
GOTERM_BP_DIRECT GO:0060070: canonical Wnt signaling pathway CDC42, FZD5, FZD7,CTNNB1, WNT4 0.004 0.769
GOTERM_BP_DIRECT GO:0019752: carboxylic acid metabolic process LDHC, UEVLD, LDHA 0.005 0.769
BMPR2, CREB1, GTF2H1,
GOTERM_BP_DIRECT GO:0006366: transcription from RNA MYBL2, ZNF831,polymerase II promoter NELFCD, SOX12, FEZF2, 0.007 1.000
POU4F1, CARF, CTCFL
GOTERM_BP_DIRECT GO:0045669: positive regulation of BMPR2, GNAS,osteoblast differentiation CTNNB1, WNT4 0.012 1.000
GOTERM_BP_DIRECT GO:0072033: renal vesicle formation CTNNB1, WNT4 0.016 1.000
GOTERM_BP_DIRECT GO:0048145: regulation of fibroblast proliferation CREB1, CTNNB1 0.016 1.000
GOTERM_BP_DIRECT GO:0051726: regulation of cell cycle L3MBTL1, ID3, E2F2,MYBL2, CDK15 0.017 1.000
Genes 2022, 13, 866 10 of 18
3.4. Selection of Signatures Based on XPEHH
The pXPEHH values for pairwise comparison of BAL with other breeds in this study
are displayed in Figure 7. To define selection signals specific to BAL, common significant
selection regions with at least three SNPs were selected from the pairwise XPEHH analyses
(Table S4a–c). From this selection, we identified 15 candidate selection regions (Table S4d)
on OAR 2, 4, 7, 13 and 17. A total of 82 genes were found in selection signals specific to
Genes 2022, 13, x FOR PEER REVIEW BAL with the XPEHH method (Table S3). Of these, 76 genes were recognized by th1e1 DofA 1V9I D
 software (Table S4e). The result shows that 13 GO TERMS were significantly enriched
(p-value < 0.05) (Table 4).
 
Figure 7. Manhattan plots for cross-population extended haplotype homozygosity (XPEHH) analyses
Figoufr(ea )7B. aMluasnhhaavttsa.nB aprldohtso kfoar,  (bcr)oBssa-lpuoshpaulvasti.oIns treixatnenadnedd( ch)aBpalloutyshpae vhso. mRouzdyag.oTsihtye d(XasPhEeHdHli)n es
analyses of (a) Balusha vs. Bardhoka, (b) Balusha vs. Istrian an≥d (c) Balusha vs. Ruda. The dashed represent the threshold level for significance at pXPEHH 5 equivalent to a p-value < 0.00001.
lines represent the threshold level for significance at pXPEHH ≥ 5 equivalent to a p-value < 0.00001. 
Common selection signals were considered as signatures of selection specific for Balusha. Overall
Common selection signals were considered as signatures of selection specific for Balusha. Overall 
thetyh eayrea mreomsto osbt voibovuios uins cinhrcohmroomsoomsoems e2s, 42 ,a4ndan 1d3.1 3.
Table 4. Enriched GO TERMS determined by DAVID from genes identified in selection regions 
specific to the BAL breed using XPEHH. Please note that GOTERM_MF refers to molecular function 
and GOTERM_BP refers to biological process. 
Category Term Genes p Value FDR 
GOTERM_MF_
GO:0005212: structural constituent of eye lens CRYGC, CRYGD, CRYGA, CRYGB 0.000 0.005 
DIRECT 
GOTERM_MF_
GO:0004181: metallocarboxypeptidase activity CPA2, CPA1, CPO, CPA5 0.000 0.005 
DIRECT 
GOTERM_MF_
GO:0004415: hyalurononglucosaminidase activity HYAL4, SPAM1 0.023 0.666 
DIRECT 
 
Genes 2022, 13, 866 11 of 18
Table 4. Enriched GO TERMS determined by DAVID from genes identified in selection regions
specific to the BAL breed using XPEHH. Please note that GOTERM_MF refers to molecular function
and GOTERM_BP refers to biological process.
Category Term Genes p Value FDR
GOTERM_MF_DIRECT GO:0005212: structural constituent of eye lens CRYGC, CRYGD,CRYGA, CRYGB 0.000 0.005
GOTERM_MF_DIRECT GO:0004181: metallocarboxypeptidase activity CPA2, CPA1, CPO, CPA5 0.000 0.005
GOTERM_MF_DIRECT GO:0004415: hyalurononglucosaminidase activity HYAL4, SPAM1 0.023 0.666
GOTERM_MF_DIRECT GO:0016790: thiolester hydrolase activity ACOT2, ACOT1 0.023 0.666
GOTERM_BP_DIRECT GO:2000601: positive regulation of Arp2/3complex-mediated actin nucleation ABI2, WASL 0.024 1.000
GOTERM_BP_DIRECT GO:0007015: actin filament organization FSCN3, WASL, LMOD2 0.026 1.000
GOTERM_BP_DIRECT GO:0006508: proteolysis CPA2, CPA1, CPO, CTSZ,ADAM23, CPA5 0.029 1.000
GOTERM_BP_DIRECT GO:0031295: T cell costimulation CD28, CTLA4, ICOS 0.030 1.000
GOTERM_BP_DIRECT GO:0007601: visual perception CRYGC, CRYGD,CRYGA, CRYGB 0.032 1.000
GOTERM_MF_DIRECT GO:0016290: palmitoyl-CoA hydrolase activity ACOT2, ACOT1 0.035 0.832
GOTERM_BP_DIRECT GO:0008154: actin polymerizationor depolymerization ABI2, WASL 0.041 1.000
GOTERM_BP_DIRECT GO:0000038: very long-chain fatty acidmetabolic process ACOT2, ACOT1 0.047 1.000
GOTERM_MF_DIRECT GO:0047617: acyl-CoA hydrolase activity ACOT2, ACOT1 0.048 0.909
Some overlapping signatures of selections were observed. Regions on sheep chromo-
somes 2 and 13 were jointly identified by both di statistics and XPEHH. A total of 34 genes
were identified in these regions (Table S3).
4. Discussion
The Kosovar BAL breed is regarded as an endangered breed [14]. Already, the ge-
ographical distribution of the BAL breed is a major negative factor in the risk status of
this breed, as the concentration of remaining animals in only two flocks makes them more
vulnerable in the unexpected event of, for example, an outbreak of infectious animal disease.
Another factor is the small number of breeding individuals remaining and the missing
exchange of animals between flocks in concordance with a missing breeding program. Even
though farmers try to keep BAL and BAR separate, this proves difficult, especially during
the extensive grazing period, which also includes the mating season. This raises difficulties
for maintaining the two breeds as purebreds.
The mean expected heterozygosities of the breeds in this study were in concordance
with the range reported for other sheep breeds in earlier studies [20–22,25]. As expected,
the BAL breed had the lowest genetic diversity of the four breeds analyzed in our study.
The low He (0.341) observed in BAL is in accordance with the geographical isolation and
points to the lack of gene flow between BAL and other breeds in southeastern Europe.
This is corroborated by the PCA (Figure 2) and the admixture analysis (Figure 3), clearly
differentiating BAL from the other analyzed Balkan breeds. For all breeds analyzed, the
effective population sizes at 5 and 50 generations ago were observed to be lower than
those reported for sheep breeds from other European or Eurasian regions, e.g., Russian [59]
and Kyrgyzstan [20] sheep breeds. The highest effective population size values (Ne5 = 35;
Ne50 = 324) were observed in the BAR breed. This fits very well with a report that this
breed is unlike the others widespread across different countries in eastern Europe [60],
but still, it points towards the risk of extinction for all analyzed breeds. The low genetic
diversity in the BAL breed is in accordance with its low effective population size values
(Ne5 = 25; Ne50 = 228). Moreover, Ne50 for the breeds in this study were within the range of
ten European and three non-European breeds reported by [21,25]. Furthermore, the mean
allelic richness, ranging from 1.780 to 1.876, is slightly lower than values observed in sheep
Genes 2022, 13, 866 12 of 18
breeds from Russia and Kyrgyzstan [20,59], but within the range observed in worldwide
breeds [25]. Even though results from microsatellite analysis are not directly comparable
with SNP data, allelic richness and private allelic richness for BAR and IST breeds show
similar patterns compared with [33].
Regarding inbreeding, various methods can be employed including FIS and FROH
whereby the latter is considered the more precise approach [28]. The BAL population had
a negative FIS, which suggests increased heterozygosity and decreased inbreeding. It is
possible that unsupervised crossbreeding with BAR may be responsible. However, the
observed FROH values and the mean ROH length indicated the presence of inbreeding,
which was expected at least for BAL with the highest FROH (0.077). The distribution
of ROH classes (Figure 4) supports the assumption that mainly historical inbreeding
events are relevant for local sheep breeds [21] and only a small proportion of inbreeding
(BAL = 14.29%, BAR = 15.22%, IST = 13.04%, RUD = 21.43%) is due to recent events within
the last five generations.
The inbreeding status is an important factor in qualifying a breed as endangered,
therefore the result for BAL indicates the urgent need to set up a conservation program to
reduce inbreeding and increase genetic diversity. The Bela Krajina Pramenka is an example
of a local breed in a southeastern European country, Slovenia, for which the setup of a
conservation program showed success, enlarging the number of purebred animals from
250 to around 900 within 17 years [61,62].
In BAL, several signatures of selection were detected on different chromosomes,
which may result from selection by humans and/or the adaptation to the environment.
Chromosomes 2, 13 and 19 showed the strongest signals based on the FST (di) method,
while ovine chromosomes 2, 4 and 13 were prominent based on the XPEHH method.
Interestingly, there is no concordance between the genes located in these regions with
previously published results of common genes within signatures of selection in ruminants
(summarized by [21]). However, four genes on ovine chromosome 2 (ALS2, AOX1) and
13 (PPP1R3D, SYCP2) were also found by Manzari et al. [63]. Based on an analysis of
Asian and European sheep breeds, several genes with the same chromosomes such as
those in our study were identified (Chr. 2: KIAA2012, SUMO1, Chr. 13: GNAS, Chr. 19:
CTNNB1) [64,65].
GO TERMS analysis identified the groups T cell co-stimulation (GO:0031295) and
Wnt signaling pathway (GO:0060071). Potential candidate genes observed within selection
signals on chromosome 2 (WNT4, FZD7, FZD5) are implicated in Wnt signaling [66,67],
which is an essential process for melanogenesis [68]. Likewise, the CREB1 gene located
within selection signatures on chromosome 2 is involved in the melanin production process
by transcriptional activation of the microphthalmia-associated transcription factor (MITF)
gene [69]. As possible candidate genes in the region of selection signals on chromosome 13,
we observed the GNAS and GNAS-AS1 genes, which are involved in the production and
regulation, respectively, of G proteins such as the G protein α-subunit (Gs), functioning
in melanogenesis as well [70–72]. Additionally, within this region of chromosome 13,
the endothelin-3 (EDN3) gene is located, which is an isoform of the endothelin-1 (EDN1)
gene [73]. Studies have implicated both genes in the melanogenesis pathway, with equal
affinity for endothelin receptor type B (ENDRB) [73–75]. EDN3 is also involved in the
production of embryonic melanocytes during fetal development [71]. Additionally, the
CTNNB1 (ß-catenin) gene observed in a selection signal region on chromosome 19 is
important for the enhancement of TCF/LEF (T cell factor/lymphoid enhancer factor) and
subsequent transcriptional activation of the MITF gene downstream of the Wnt signaling
pathway [67,68]. Our findings reveal that the agouti signaling protein (ASIP) and MITF
genes associated with coat color in animals are located at a distance of 3.1 Mb and 2.2 Mb
from the nearest selection signal on chromosome 13 and chromosome 19, respectively.
Therefore, due to their participation in the melanogenesis pathway (Figure 8), some or
all of the mentioned putative candidate genes on chromosomes 2, 13 and 19 may have
an influence on coat color differences between BAL and the other analyzed breeds. This
Genes 2022, 13, 866 13 of 18
Genes 2022, 13, x FOR PEER REVIEW 14 of 19 
 
fits well with the predominant difference between BAL and the other breeds, which is the
distinct black-colored head of BAL that makes this breed easily recognizable (Figure 1).
 
Figure 8. Melanogenesis pathway (KEGG database [76], map04916, truncated to the part relevant for
Figure 8. Melanogenesis pathway (KEGG database [76], map04916, truncated to the part relevant 
fothr itshsitsu sdtyu)d. yT)h. eTgheen geesnbeesl obnegloinnggitnogt htoe pthroet epirnostehiingsh hligighhtelidghinterde dinw reerde wideernet iifideedntinifiseedle icnti osenlescigtinoanl s
sisgpneaclifis cspfoerciBfiacl ufoshr aB.aPlulesahsae. nPolteeatshea ntontaem theaotf nthame gee onfe sthaen gdetnheesc aonrrde stphoen cdoirnrgespproontediinnsga prerontoetinalsw aarye s
ntohte aslwamayes( eth.ge. ,sFamrizez (le.dg.=, FFrZizDzl7e,dF =Z DFZ5D). 7A, lFlZinDv5o)l.v Aeldl ipnrvootelvinesd aprreodteisinpsla ayred diinsplliagyhetdg rinee lnigbhot xgerse.en 
boxes. 
The activation and survival of T cells are associated with three genes (CD28, CTLA4,
ICOTSh) elo accattievdatiinona arengdi osnurwviivthals oefl eTc tcieolnlss aigren aalsssoocniactherdo mwiotsho tmhree2e g[7e7n]e(sF (iCguDr2e89, )C. TTLcAe4ll,s ,
ICasOpSa) rltoocfattehde iimn am ruengeiosny swteimth, saeildectthieonp rsoigpnearlfsu onnc tciohnroomf oimsommuen 2e [r7e7sp] o(Fnisgeusr,eh 9o)m. Teo csetlalssi,s 
aasn pdarmt oemf tohrey im[7m8,7u9n]e.  Ssytustdeimes, haiadv ethseh porwonpetrh afut nCcDti2o8n porfo ivmidmeusnceo srteismpuonlasteosr,y hsoimgneaolsst,aasnisd 
apnrdo moemtesorTy c[7e8ll,7p9r]o. lSifteurdaiteios nh,acvyet oshkoinwinn pthraotd CuDct2io8n paronvdidceesll csousrtvimivualla[t7o7ry,8 0s]ig. nLailkse, wanisde ,
pIrCoOmSotisesa Tk ncoewll npcroo-lisfteimrautiloanto, rcoyftoTkcinelilna pctriovdauticotnio, np raonlidfe rcaetlilo nsuarnvdivdailf f[e7r7e,n80ti]a. tiLoinke[w77i,s8e1, ].
ICFuOrSth ise ram konroew, ena croli-esrtirmepuolarttosrr eovf eTa cleedll athctaitvtahtieoCn,T pLrAol4ifiesraathioonm aonldo gdiofffeCreDn2t8iaatinodn a[7c7ts,8a1s].a 
Fsuurpthperremssoorre,o efatrhleielra rtetepro[r8t0s ,r8e2v].eaCleodn stehqaut ethnetl yC,TinLcAr4ea isse ad hCoTmLoAlo4ge oxfp CreDss2i8o nanldea adcstst oast hae 
stueprpmriensastoiro nofo tfhTe lcaetltlear c[t8iv0,a8t2io].n C. oInntseerqeusteinntglyly, ,itnhcirseastsuedd yCsThLoAw4s etxhpartefsascitoonr sleoafdasd taop tthivee 
teimrmminuantiitoyns hoof wT dceiflfle arecnticveastiionns. igInntaetruersetsinogflyse, ltehcitsi osntubdeytw seheonwBsA tLhaatn fdacBtoArRs roefa areddapintivthee 
imsammuenliotcya stihoonw.  differences in signatures of selection between BAL and BAR reared in the 
same lFoicnaatliloy,nw. e identified, as a possible candidate gene, the CTSZ gene (located within se-
lection signals on chromosome 13), which was associated with tuberculosis (Mycobacterium
Finally, we identified, as a possible candidate gene, the CTSZ gene (located within 
tuberculosis) susceptibility in humans [83–85]. Investigation of a possible association of
selection signals on chromosome 13), which was associated with tuberculosis (Mycobac-
this potential candidate gene with paratuberculosis (Mycobacterium avium subspecies
terium tuberculosis) susceptibility in humans [83–85]. Investigation of a possible associa-
paratuberculosis) susceptibility in BAL sheep would be a promising task, as there is an in-
tion of this potential candidate gene with paratuberculosis (Mycobacterium avium sub-
creasing interest in breeding for resistance to infectious diseases in farm animals, including
species paratuberculosis) susceptibility in BAL sheep would be a promising task, as there 
paratuberculosis in ruminants [86]. As it was shown for small ruminant lentiviruses, pro-
ist eacnti vinecvreaarisainngts iangtearinesstt iinnf becrteieoduisndg ifsoera sreessicsatnanscheo two ainbferecetido-ussp edciisfiecasdeiss tirnib fuatrimon aen.igm.,a[l8s7, ].
inTchleudfiinndgi pnagrsaotuf btheriscusltousdisy isnu rgugmesint athnatst [a8t6te].n Atiso nit wshaosu slhdobwenp faoird stmoatlhl eruBmAiLnasnhte leepntwiviit-h
ruthseesp, opsrsoibteiclittiyveo fvuarsiianngtsth aegmainassta insofeucrtcieoufos rdmisoeadseersn cbarne esdhionwg gao barlesesdu-cshpeacsifriecs disitsatnricbeut-o
tiionnfe ec.tgio.,u [s87d]i.s Teahsee fsi,nwdhinicghs iosf athfiusr sthtuedr yre sausgognefsotr tthhaet caottnesnetrivoant isohnouolfdl obcea lpbaridee tdos t.he BAL 
sheep with the possibility of using them as a source for modern breeding goals such as 
resistance to infectious diseases, which is a further reason for the conservation of local 
breeds. 
 
Genes 2022, 13, x FOR PEER REVIEW 15 of 19 
Gene s 2022, 13, 866 14 of 18
 
FFiigguurree 99. .HHumuamna Tn cTellc reelclerpetcoer psitgonrasliinggn aplaitnhgwpaya t(hKwEGaGy  (dKatEaGbaGse d[7a6t]a, bhassae04[67660]),. hThsae 0g4e6n6e0s )b.eT- he genes
bloenloginnggi ntog tthoe tphreotperinost ehinigshhligighhtelidg hint erdedi nwreered iwdeenrteifiedde ninti fiseeldectiinonse sliegcntaiolsn sspiegcnifaicl sfosrp Becailfiuschfao. r Balusha.
Allll iinnvvoollvveedd pprortoetienisn asrae rdeisdpilsapyleady iend ai nligahtl iggrheteng rbeoexn. box.
5. Conncclulussioionns s
Reelilaiabblel eknkonwowledlegde goef tohfet hdievedrisviteyr sbiettywbeetnw aenedn wainthdinw birteheidns bisr,e aemdsonisg, oatmheorns,g a others, a
rreeqquuiirreemeennt tfofro arna enffiecfifiencti ecnont sceorvnasteiorvna sttiroantegsytr. aTtheegrye.foTreh, eoruerf sotrued,yo ounr tshteu gdeyneotnic dthi-e genetic
dveivrseirtysi toyf tohfet hKeosKovoasro vBaArLB bAreLedb rpereodvipdreosv ai dbeassica wbaosrikc fworo trhke fporrestehrevaptrieosne orvf athtiiso nraoref this rare
llooccaall bbrreeeedd aanndd ininddiciactaetse sopotpiotinosn fsofro irtsi tfsutfuurteu rdeevdeelvoeplmopenmt.e Tnht.e Tchome cboinmatbioinna otifo tnheo fret-he results
fsruoltms frsoigmn saitgunraetsuroefs soef lseecleticotinona nanddt thhee RROOHH aannalaylsyesse psrpovroidveisd ae sbeattbere uttnedreursntadnedrsintagn ding of
of how the BAL breed is genetically distinguishable from the other considered sheep 
how the BAL breed is genetically distinguishable from the other considered sheep breeds
breeds originating from the Balkan Peninsula. GO TERMS analysis revealed two interest-
oinrgig ginroautipnsg off rgoemnest hinevBolavlekda nin Pmeenliannsouglean.esGisO anTdE TR cMelSl caon-satliymsuislatrieovne, awlehdichtw moigihntt eresting
ginrtoeurfpesreo wf gitehn BeAs Lin’sv uonlviqeudei ncomlore lpaantotegrenn aensdis tahne dimTmcuenllec soy-sstteimm.u Tlahtei olant,tewr hcoicuhldm piogihntt interfere
wtoiwtharBdA a Lc’hsaunngeiqdu deisceoalsoer spusacteteprtinbialintyd, wthheicihm wmouunlde msyasktee mBA. LT ha evlaaluttaebrlec oguenldetipco rien-t toward
asocuhracen gfoerd brdeiesdeiansge fsour sacneimptailb hileiatylt,hw. hich would make BAL a valuable genetic resource for
breedSincge fourra anniamlyaslesh eshaoltwh.ed that BAL sheep differ strongly from the other compared 
PramSeninkcae sohueerpa bnraeleydsse, sanshdo twhee bdretheda tisB nAoLt ksnhoewepn tdoi fofcecrusrt rino nogthlyer fcroumnttrhiees,o ctohnesrerc-ompared
Pvaratimone onfk tahes hBeAeLp sbhreeeepd ssh,oaunldd bthee gbivreene dpriisornitoyt okvneor wotnhetro roegcciounrailn broetehdesr. countries, conserva-
tion of the BAL sheep should be given priority over other regional breeds.
Supplementary Materials: The following supporting information can be downloaded at: 
www.mdpi.com/xxx/s1, Table S1: Information on sample numbers and SNP data. Table S2: (a–b) 
SSiugpnpifilceamnte rnetgairoynsM ofa steelreicatliso:n Tfrhoemf Foll oawnaiST lnygsiss. uTpabploer Sti3n: gAninnfootartmeda tgieonnesc aidnenbteifidedo winn sleoleacdtieodn at: https:
/si/gwnawtuwre.ms sdppeci.icfiocm to/ tahret icBlAe/L 1b0r.e3e3d9 0d/etgeecnteeds 1b3y0 F5S0T8 (6d6i)/ asn1d,  TXaPbElHe HS1. :TIanbfloer Sm4a: t(iao–ne) oSnigsnaimficpanlet numbers
arengdioSnNs Pof dsealteac.tiToanb flreomS2 X: (PaE–Hb)HS aignnaliyfiscias.n t regions of selection from FST analysis. Table S3: Annotated
gAeunthesori dCeonnttirfiibeudtiionnsse: lCeocnticoenptsuiaglnizaattuiorne,s Os.pOe.Aci.fi, Rc .tSo., tHh.eB.B aAndL Gb.rLe.e; Fdodrmetaelc atneadlybsyis,F OST.O(.dAi). aanndd XPEHH.
TMa.bMle.; SF4u:n(dai–neg) aSciqguniisfiitcioann,t Wre.gKi.o; Pnrsoojefcst ealdemctiinoinstfrraotimon,X GP.ELH.; RHesaonuarlcyessi, sH. .B., H.M. and K.B.; Su-
ApeurtvhisoiornC, Gon.Lt.r;i Vbiusutiaolnizsa:tioCno, nOc.Oep.Atu. aanlidz aRt.iSo.n; W, Orit.iOng.A—.,orRig.Sin.,alH d.rBa.fta, Ond.OG.A..L a.n; dF oRr.Sm.; aWl raitninagly—sis, O.O.A.
review and editing, O.O.A., R.S., H.B., M.S. and G.L. All authors will be informed about each step 
and M.M.; Funding acquisition, W.K.; Project administration, G.L.; Resources, H.B., H.M. and
K.B.; Supervision, G.L.; Visualization, O.O.A. and R.S.; Writing—original draft, O.O.A. and R.S.;
Writing—review and editing, O.O.A., R.S., H.B., M.S. and G.L. All authors will be informed about
 each step of manuscript processing including submission, revision, revision reminder, etc. via emails
from our system or assigned Assistant Editor. All authors have read and agreed to the published
version of the manuscript.
Genes 2022, 13, 866 15 of 18
Funding: The genotyping of the Balusha and Bardhoka samples was financed by the SAVE Founda-
tion e.V. CH-9000 St. Gallen.
Institutional Review Board Statement: All animals were sampled according to the local regulations,
performed, and inspected by the Veterinary and Food Agency of Kosovo (AUV, permit; 173829; date
25.09.2019), with the consent of animal owners.
Informed Consent Statement: Not applicable.
Data Availability Statement: Data supporting reported results are provided in Supplementary Tables S1–S3.
Acknowledgments: We thank the owners of the sheep for providing samples and data.
Conflicts of Interest: The funders had no role in the design of the study; in the collection, analyses,
or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
G.L. is one of the Special Issue editors of MDPI Genes but have not in any way been involved in or
interacted with the journal’s review process or editorial decision-making. The guest editors were
blinded to the review process.
References
1. Chessa, B.; Pereira, F.; Arnaud, F.; Amorim, A.; Goyache, F.; Mainland, I.; Kao, R.R.; Pemberton, J.M.; Beraldi, D.; Stear, M.J.; et al.
Revealing the history of sheep domestication using retrovirus integrations. Science 2009, 324, 532–536. [CrossRef] [PubMed]
2. Ciani, E.; Mastrangelo, S.; Da Silva, A.; Marroni, F.; Ferenčaković, M.; Ajmone-Marsan, P.; Baird, H.; Barbato, M.; Colli, L.;
Delvento, C.; et al. On the origin of European sheep as revealed by the diversity of the Balkan breeds and by optimizing
population-genetic analysis tools. Genet. Sel. Evol. 2020, 52, 25. [CrossRef] [PubMed]
3. Tresset, A.; Vigne, J.D. Substitution of Substitution of Species, Techniques and Symbols at the Mesolithic/Neolithic Transition in Western
Europe; British Academy: London, UK, 2007; pp. 189–210.
4. Bocquet-Appel, J.-P.; Naji, S.; Linden, M.V.; Kozlowski, J.K. Detection of diffusion and contact zones of early farming in Europe
from the space-time distribution of 14C dates. J. Archaeol. Sci. 2009, 36, 807–820. [CrossRef]
5. De Vareilles, A.; Bouby, L.; Jesus, A.; Martin, L.; Rottoli, M.; Vander Linden, M.; Antolín, F. One sea but many routes to Sail.
The early maritime dispersal of Neolithic crops from the Aegean to the western Mediterranean. J. Archaeol. Sci. Rep. 2020,
29, 102140. [CrossRef]
6. Marković, B.; Dovč, P.; Marković, M.; Radonjić, D.; Adakalić, M.; Simčič, M. Differentiation of some Pramenka sheep breeds based
on morphometric characteristics. Arch. Anim. Breed. 2019, 62, 393–402. [CrossRef]
7. FAO. The Second Report on the State of the World’s Animal Genetic Resources for Food and Agriculture, Rome. 2015. Available
online: www.fao.org/3/a-i4787e/index.html (accessed on 24 February 2022).
8. Taberlet, P.; Valentini, A.; Rezaei, H.R.; Naderi, S.; Pompanon, F.; Negrini, R.; Ajmone-Marsan, P. Are cattle, sheep, and goats
endangered species? Mol. Ecol. 2008, 17, 275–284. [CrossRef]
9. Federal Office for Agriculture and Food. Red List: List of Native Livestock Breeds: Sheep. Available online: https:
//tgrdeu.genres.de/en/red-list/list-of-native-livestock-breeds/list-of-genetics/?tx_sttgrdeu_roteliste%5Baction%5D=
listGenetik&tx_sttgrdeu_roteliste%5Ba_id%5D=5&tx_sttgrdeu_roteliste%5Bcontroller%5D=Roteliste&cHash=356115b92cb1
ab653752bd5316ce286f (accessed on 12 April 2022).
10. Bytyqi, H.; Mehmeti, H. Bardhoka Strain. In Catalogue of West Balkan Pramenka Sheep Breed Types: Identification and Conservation
of Animal Genetic Resources in South Eastern Europe; Faculty of Agricultural Sciences and Food, University of Skopje: Skopje,
Macedonia, 2006.
11. Bytyqi, H.; Baumung, R.; Mehmeti, H.; Fuerst-Waltl, B. Phenotypic characterization and description of production systems of
autochthonous sheep breeds in Kosovo. Anim. Genet. Resour. 2014, 54, 163–170. [CrossRef]
12. Porter, V.; Alderson, L.; Hall, S.J.G.; Sponenberg, D.P. Mason’s World Encyclopedia of Livestock Breeds and Breeding; CABI: Wallingford,
UK, 2016; ISBN 978-1-84593-466-8.
13. Hoda, A.; Ajmone, P. Genetic Characterization of Albanian Sheep Breeds by Microsatellite Markers. In Analysis of Genetic Variation
in Animals; Caliskan, M., Ed.; InTech: London, UK, 2012; ISBN 978-953-51-0093-5.
14. SAVE Foundation. Schafhaltung und gelebtes kulturelles Erbe im Kosovo. In SAVE e-New—Safeguard for Agricultural Varieties
in Europe: Der Vierteljährliche Informationsdienst der Europäischen SAVE Foundation; SAVE Foundation: St. Gallen, Switzerland,
2017; pp. 6–9.
15. Notter, D.R. The importance of genetic diversity in livestock populations of the future. J. Anim. Sci. 1999, 77, 61–69. [CrossRef]
16. Groeneveld, L.F.; Lenstra, J.A.; Eding, H.; Toro, M.A.; Scherf, B.; Pilling, D.; Negrini, R.; Finlay, E.K.; Jianlin, H.;
Groeneveld, E.; et al. Genetic diversity in farm animals–a review. Anim. Genet. 2010, 41 (Suppl. 1), 6–31. [CrossRef]
17. Gandini, G.C.; Villa, E. Analysis of the cultural value of local livestock breeds: A methodology. J. Anim. Breed. Genet. 2003,
120, 1–11. [CrossRef]
18. Solti, L.; Crichton, E.G.; Loskutoff, N.M.; Cseh, S. Economical and ecological importance of indigenous livestock and the
application of assisted reroduction to their preservation. Theriogenology 2000, 53, 149–162. [CrossRef]
Genes 2022, 13, 866 16 of 18
19. Kijas, J.W.; Townley, D.; Dalrymple, B.P.; Heaton, M.P.; Maddox, J.F.; McGrath, A.; Wilson, P.; Ingersoll, R.G.; McCulloch, R.;
McWilliam, S.; et al. A genome wide survey of SNP variation reveals the genetic structure of sheep breeds. PLoS ONE 2009,
4, e4668. [CrossRef]
20. Deniskova, T.; Dotsev, A.; Lushihina, E.; Shakhin, A.; Kunz, E.; Medugorac, I.; Reyer, H.; Wimmers, K.; Khayatzadeh, N.; Sölkner,
J.; et al. Population Structure and Genetic Diversity of Sheep Breeds in the Kyrgyzstan. Front. Genet. 2019, 10, 1311. [CrossRef]
21. Signer-Hasler, H.; Burren, A.; Ammann, P.; Drögemüller, C.; Flury, C. Runs of homozygosity and signatures of selection:
A comparison among eight local Swiss sheep breeds. Anim. Genet. 2019, 50, 512–525. [CrossRef]
22. Meyermans, R.; Gorssen, W.; Wijnrocx, K.; Lenstra, J.A.; Vellema, P.; Buys, N.; Janssens, S. Unraveling the genetic diversity of
Belgian Milk Sheep using medium-density SNP genotypes. Anim. Genet. 2020, 51, 258–265. [CrossRef]
23. Al-Mamun, H.A.; Clark, S.A.; Kwan, P.; Gondro, C. Genome-wide linkage disequilibrium and genetic diversity in five populations
of Australian domestic sheep. Genet. Sel. Evol. 2015, 47, 90. [CrossRef]
24. Edea, Z.; Dessie, T.; Dadi, H.; Do, K.-T.; Kim, K.-S. Genetic Diversity and Population Structure of Ethiopian Sheep Populations
Revealed by High-Density SNP Markers. Front. Genet. 2017, 8, 218. [CrossRef]
25. Kijas, J.W.; Lenstra, J.A.; Hayes, B.; Boitard, S.; Porto Neto, L.R.; San Cristobal, M.; Servin, B.; McCulloch, R.; Whan, V.; Gietzen,
K.; et al. Genome-wide analysis of the world’s sheep breeds reveals high levels of historic mixture and strong recent selection.
PLoS Biol. 2012, 10, e1001258. [CrossRef]
26. Curik, I.; Ferenčaković, M.; Sölkner, J. Inbreeding and runs of homozygosity: A possible solution to an old problem. Livest. Sci.
2014, 166, 26–34. [CrossRef]
27. Burren, A.; Neuditschko, M.; Signer-Hasler, H.; Frischknecht, M.; Reber, I.; Menzi, F.; Drögemüller, C.; Flury, C. Genetic diversity
analyses reveal first insights into breed-specific selection signatures within Swiss goat breeds. Anim. Genet. 2016, 47, 727–739.
[CrossRef]
28. Keller, M.C.; Visscher, P.M.; Goddard, M.E. Quantification of inbreeding due to distant ancestors and its detection using dense
single nucleotide polymorphism data. Genetics 2011, 189, 237–249. [CrossRef]
29. Purfield, D.C.; McParland, S.; Wall, E.; Berry, D.P. The distribution of runs of homozygosity and selection signatures in six
commercial meat sheep breeds. PLoS ONE 2017, 12, e0176780. [CrossRef]
30. Akey, J.M.; Ruhe, A.L.; Akey, D.T.; Wong, A.K.; Connelly, C.F.; Madeoy, J.; Nicholas, T.J.; Neff, M.W. Tracking footprints of
artificial selection in the dog genome. Proc. Natl. Acad. Sci. USA 2010, 107, 1160–1165. [CrossRef]
31. Sabeti, P.C.; Varilly, P.; Fry, B.; Lohmueller, J.; Hostetter, E.; Cotsapas, C.; Xie, X.; Byrne, E.H.; McCarroll, S.A.; Gaudet, R.; et al.
Genome-wide detection and characterization of positive selection in human populations. Nature 2007, 449, 913–918. [CrossRef]
32. Qanbari, S.; Simianer, H. Mapping signatures of positive selection in the genome of livestock. Livest. Sci. 2014, 166, 133–143.
[CrossRef]
33. Cinkulov, M.; Popovski, Z.; Porcu, K.; Tanaskovska, B.; Hodzić, A.; Bytyqi, H.; Mehmeti, H.; Margeta, V.; Djedović, R.;
Hoda, A.; et al. Genetic diversity and structure of the West Balkan Pramenka sheep types as revealed by microsatellite and
mitochondrial DNA analysis. J. Anim. Breed. Genet. 2008, 125, 417–426. [CrossRef]
34. Peter, C.; Bruford, M.; Perez, T.; Dalamitra, S.; Hewitt, G.; Erhardt, G. Genetic diversity and subdivision of 57 European and
Middle-Eastern sheep breeds. Anim. Genet. 2007, 38, 37–44. [CrossRef]
35. Manichaikul, A.; Mychaleckyj, J.C.; Rich, S.S.; Daly, K.; Sale, M.; Chen, W.-M. Robust relationship inference in genome-wide
association studies. Bioinformatics 2010, 26, 2867–2873. [CrossRef]
36. Chang, C.C.; Chow, C.C.; Tellier, L.C.; Vattikuti, S.; Purcell, S.M.; Lee, J.J. Second-generation PLINK: Rising to the challenge of
larger and richer datasets. Gigascience 2015, 4, 7. [CrossRef]
37. Wickham, H. ggplot2: Elegant Graphics for Data Analysis; Springer International Publishing: New York, NY, USA, 2016;
ISBN 9783319242774.
38. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria,
2021. Available online: https://www.R-project.org/ (accessed on 22 November 2021).
39. Goudet, J. HIERFSTAT, a package for R to compute and test hierarchical F-statistics. Mol. Ecol. Notes 2005, 5, 184–186. [CrossRef]
40. Szpiech, Z.A.; Jakobsson, M.; Rosenberg, N.A. ADZE: A rarefaction approach for counting alleles private to combinations of
populations. Bioinformatics 2008, 24, 2498–2504. [CrossRef]
41. Nei, M.; Chesser, R.K. Estimation of fixation indices and gene diversities. Ann. Human Genet. 1983, 47, 253–259. [CrossRef]
[PubMed]
42. Barbato, M.; Orozco-terWengel, P.; Tapio, M.; Bruford, M.W. SNeP: A tool to estimate trends in recent effective population size
trajectories using genome-wide SNP data. Front. Genet. 2015, 6, 109. [CrossRef] [PubMed]
43. Sved, J.A.; Feldman, M.W. Correlation and probability methods for one and two loci. Theor. Popul. Biol. 1973, 4, 129–132.
[CrossRef]
44. Corbin, L.J.; Liu, A.Y.H.; Bishop, S.C.; Woolliams, J.A. Estimation of historical effective population size using linkage disequilibria
with marker data. J. Anim. Breed. Genet. 2012, 129, 257–270. [CrossRef] [PubMed]
45. Alexander, D.H.; Novembre, J.; Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009,
19, 1655–1664. [CrossRef]
46. Francis, R.M. pophelper: An R package and web app to analyse and visualize population structure. Mol. Ecol. Resour. 2017,
17, 27–32. [CrossRef]
Genes 2022, 13, 866 17 of 18
47. McQuillan, R.; Leutenegger, A.-L.; Abdel-Rahman, R.; Franklin, C.S.; Pericic, M.; Barac-Lauc, L.; Smolej-Narancic, N.; Janicijevic,
B.; Polasek, O.; Tenesa, A.; et al. Runs of homozygosity in European populations. Am. J. Hum. Genet. 2008, 83, 359–372. [CrossRef]
48. Mastrangelo, S.; Ciani, E.; Sardina, M.T.; Sottile, G.; Pilla, F.; Portolano, B. Runs of homozygosity reveal genome-wide autozygosity
in Italian sheep breeds. Anim. Genet. 2018, 49, 71–81. [CrossRef]
49. Meyermans, R.; Gorssen, W.; Buys, N.; Janssens, S. How to study runs of homozygosity using PLINK? A guide for analyzing
medium density SNP data in livestock and pet species. BMC Genom. 2020, 21, 94. [CrossRef]
50. Mastrangelo, S.; Portolano, B.; Di Gerlando, R.; Ciampolini, R.; Tolone, M.; Sardina, M.T. Genome-wide analysis in endangered
populations: A case study in Barbaresca sheep. Animal 2017, 11, 1107–1116. [CrossRef] [PubMed]
51. Weir, B.S.; Cockerham, C.C. Estimating F-statistics for the analysis of population structure. Evolution 1984, 38, 1358–1370.
[CrossRef] [PubMed]
52. Gautier, M.; Naves, M. Footprints of selection in the ancestral admixture of a New World Creole cattle breed. Mol. Ecol. 2011,
20, 3128–3143. [CrossRef] [PubMed]
53. Scheet, P.; Stephens, M. A fast and flexible statistical model for large-scale population genotype data: Applications to inferring
missing genotypes and haplotypic phase. Am. J. Hum. Genet. 2006, 78, 629–644. [CrossRef] [PubMed]
54. Karolchik, D.; Hinrichs, A.S.; Kent, W.J. The UCSC Genome Browser. Curr. Protoc. Hum. Genet. 2011, 71, 18.6.1–18.6.33. [CrossRef]
55. Huang, W.D.; Sherman, B.T.; Lempicki, R.A. Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis
of large gene lists. Nucleic Acids Res. 2009, 37, 1–13. [CrossRef]
56. Huang, W.D.; Sherman, B.T.; Lempicki, R.A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics
resources. Nat. Protoc. 2009, 4, 44–57. [CrossRef]
57. Benjamini, Y.; Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R.
Stat. Soc. Ser. B 1995, 57, 289–300. [CrossRef]
58. Voight, B.F.; Kudaravalli, S.; Wen, X.; Pritchard, J.K. A map of recent positive selection in the human genome. PLoS Biol. 2006,
4, e72. [CrossRef]
59. Deniskova, T.E.; Dotsev, A.V.; Selionova, M.I.; Kunz, E.; Medugorac, I.; Reyer, H.; Wimmers, K.; Barbato, M.; Traspov, A.A.; Brem,
G.; et al. Population structure and genetic diversity of 25 Russian sheep breeds based on whole-genome genotyping. Genet. Sel.
Evol. 2018, 50, 29. [CrossRef]
60. Kusza, S.; Ivankovic, A.; Ramljak, J.; Nagy, I.; Jávor, A.; Kukovics, S. Genetic structure of Tsigai, Ruda, Pramenka and other local
sheep in Southern and Eastern Europe. Small Rumin. Res. 2011, 99, 130–134. [CrossRef]
61. Bojkovsky, D. Bela Krajina Pramenka Sheep: The Best Lamb Meat in Slovenia. Available online: https://www.fao.org/dad-is/
success-story/detail/en/c/1193631/ (accessed on 17 February 2022).
62. Tomazic, D. The Bela Krajina Pramenka Sheep—Conservation Strategies. Agric. Conspec. Sci. 2003, 68, 323–328.
63. Manzari, Z.; Mehrabani-Yeganeh, H.; Nejati-Javaremi, A.; Moradi, M.H.; Gholizadeh, M. Detecting selection signatures in three
Iranian sheep breeds. Anim. Genet. 2019, 50, 298–302. [CrossRef] [PubMed]
64. Eydivandi, S.; Roudbar, M.A.; Ardestani, S.S.; Momen, M.; Sahana, G. A selection signatures study among Middle Eastern and
European sheep breeds. J. Anim. Breed. Genet. 2021, 138, 574–588. [CrossRef] [PubMed]
65. Eydivandi, S.; Roudbar, M.A.; Karimi, M.O.; Sahana, G. Genomic scans for selective sweeps through haplotype homozygosity
and allelic fixation in 14 indigenous sheep breeds from Middle East and South Asia. Sci. Rep. 2021, 11, 2834. [CrossRef] [PubMed]
66. MacDonald, B.T.; Tamai, K.; He, X. Wnt/β-catenin signaling: Components, mechanisms, and diseases. Dev. Cell 2009, 17, 9–26.
[CrossRef]
67. Van Amerongen, R.; Nusse, R. Towards an integrated view of Wnt signaling in development. Development 2009, 136, 3205–3214.
[CrossRef]
68. D’Mello, S.A.N.; Finlay, G.J.; Baguley, B.C.; Askarian-Amiri, M.E. Signaling Pathways in Melanogenesis. Int. J. Mol. Sci. 2016,
17, 1144. [CrossRef]
69. Arora, N.; Siddiqui, E.M.; Mehan, S. Involvement of adenylate cyclase/cAMP/CREB and SOX9/MITF in melanogenesis to
prevent vitiligo. Mol. Cell. Biochem. 2021, 476, 1401–1409. [CrossRef]
70. Plagge, A.; Kelsey, G.; Germain-Lee, E.L. Physiological functions of the imprinted Gnas locus and its protein variants Galpha(s)
and XLalpha(s) in human and mouse. J. Endocrinol. 2008, 196, 193–214. [CrossRef]
71. Tang, L.; Su, M.; Zhang, Y.; Ip, W.; Martinka, M.; Huang, C.; Zhou, Y. Endothelin-3 Is Produced by Metastatic Melanoma Cells
and Promotes Melanoma Cell Survival. J. Cutan. Med. Surg. 2008, 12, 64–70. [CrossRef] [PubMed]
72. Weinstein, L.S.; Liu, J.; Sakamoto, A.; Xie, T.; Chen, M. Minireview: GNAS: Normal and abnormal functions. Endocrinology 2004,
145, 5459–5464. [CrossRef] [PubMed]
73. Sakurai, T.; Yanagisawa, M.; Masaki, T. Molecular characterization of endothelin receptors. Trends Pharmacol. Sci. 1992, 13, 103–108.
[CrossRef] [PubMed]
74. Saldana-Caboverde, A.; Kos, L. Roles of endothelin signaling in melanocyte development and melanoma. Pigment. Cell Melanoma
Res. 2010, 23, 160–170. [CrossRef]
75. Imokawa, G.; Yada, Y.; Kimura, M. Signalling mechanisms of endothelin-induced mitogenesis and melanogenesis in human
melanocytes. Biochem. J. 1996, 314 Pt 1, 305–312. [CrossRef]
76. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 28, 27–30. [CrossRef]
77. Smith-Garvin, J.E.; Koretzky, G.A.; Jordan, M.S. T cell activation. Annu. Rev. Immunol. 2009, 27, 591–619. [CrossRef]
Genes 2022, 13, 866 18 of 18
78. Bonilla, F.A.; Oettgen, H.C. Adaptive immunity. J. Allergy Clin. Immunol. 2010, 125, S33–S40. [CrossRef]
79. Kumar, B.V.; Connors, T.J.; Farber, D.L. Human T Cell Development, Localization, and Function throughout Life. Immunity 2018,
48, 202–213. [CrossRef]
80. Beyersdorf, N.; Kerkau, T.; Hünig, T. CD28 co-stimulation in T-cell homeostasis: A recent perspective. Immunotargets Ther. 2015,
4, 111–122. [CrossRef]
81. Wikenheiser, D.J.; Stumhofer, J.S. ICOS Co-Stimulation: Friend or Foe? Front. Immunol. 2016, 7, 304. [CrossRef] [PubMed]
82. Walker, L.S.K.; Sansom, D.M. The emerging role of CTLA4 as a cell-extrinsic regulator of T cell responses. Nat. Rev. Immunol.
2011, 11, 852–863. [CrossRef] [PubMed]
83. Adams, L.A.; Möller, M.; Nebel, A.; Schreiber, S.; van der Merwe, L.; van Helden, P.D.; Hoal, E.G. Polymorphisms in MC3R
promoter and CTSZ 3’UTR are associated with tuberculosis susceptibility. Eur. J. Hum. Genet. 2011, 19, 676–681. [CrossRef]
[PubMed]
84. Cooke, G.S.; Campbell, S.J.; Bennett, S.; Lienhardt, C.; McAdam, K.P.W.J.; Sirugo, G.; Sow, O.; Gustafson, P.; Mwangulu, F.; van
Helden, P.; et al. Mapping of a novel susceptibility locus suggests a role for MC3R and CTSZ in human tuberculosis. Am. J. Respir.
Crit. Care Med. 2008, 178, 203–207. [CrossRef]
85. Hashemi, M.; Eskandari-Nasab, E.; Moazeni-Roodi, A.; Naderi, M.; Sharifi-Mood, B.; Taheri, M. Association of CTSZ rs34069356
and MC3R rs6127698 gene polymorphisms with pulmonary tuberculosis. Int. J. Tuberc. Lung Dis. 2013, 17, 1224–1228. [CrossRef]
86. Purdie, A.C.; Plain, K.M.; Begg, D.J.; de Silva, K.; Whittington, R.J. Candidate gene and genome-wide association studies of
Mycobacterium avium subsp. paratuberculosis infection in cattle and sheep: A review. Comp. Immunol. Microbiol. Infect. Dis.
2011, 34, 197–208. [CrossRef]
87. Molaee, V.; Otarod, V.; Abdollahi, D.; Lühken, G. Lentivirus Susceptibility in Iranian and German Sheep Assessed by Determina-
tion of TMEM154 E35K. Animals 2019, 9, 685. [CrossRef]