GENETIC CHARACTERIZATION OF THE YUGOSLAVIAN SHEPHERD DOG – SHARPLANINA, A LIVESTOCK GUARD DOG FROM THE WESTERN BALKANS

Department of Animal Breeding and Genetics, Faculty of Veterinary Medicine, University of Belgrade, Bul. oslobođenja 18, PO Box 310, 11000 Belgrade, Serbia; Department of Reproduction, Fertility and Artificial Insemination, Faculty of Veterinary Medicine, University of Belgrade, Bul. oslobođenja 18, PO Box 310, 11000 Belgrade, Serbia; Department of Biology, Faculty of Veterinary Medicine, University of Belgrade, Belgrade, Serbia; Ministry of Agriculture, Forestry and Water Management of the Republic of Serbia, Plant Protection Directorate, Department for Plant Variety Registration, Omladinskih brigada 1, SIV 3, 11070 Belgrade, Serbia; Institute of Molecular Genetics and Genetic


IntroductIon
Domestic dog, Canis lupus familiaris L., is the first domesticated animal species characterized today by unprecedented conformational and functional variability [1,2]. It comprises more than 400 distinct breeds, most of which have been developed in the mid-19 th century when reproductive isolation between breeds has been formalized and breed standards defined [3,4]. The breeding practice in dogs, however, is rather specific, because of breeding in closed populations, selective pressures towards desired traits and extensive use of popular sires. That, along with the historical variations in population sizes, have resulted in the considerable inbreeding and loss of genetic variability in certain purebred dog populations [5][6][7][8], as well as in the increased incidence of inherited diseases [9][10][11]. Although genetic make-up of numerous dog breeds has been studied to date [4,[11][12][13][14][15][16][17][18][19][20][21][22][23], the knowledge regarding genetic structure of non-cosmopolitan dog breeds and local canine populations is still grossly lacking [21,24].
Yugoslavian Shepherd Dog -Sharplanina (YSD), is a livestock guard dog from the Western Balkans. This large pastoral dog breed has a long historical presence in the Sharplanina Mountain region. It had been registered by the Fédération Cynologique Internationale (FCI) since 1939 under the designation "Ilirski Ovcar" (Illyrian Shepherd Dog). In 1957, the FCI accepted a motion proposed by the Yugoslavian Federation of Cynology to change the name of the breed to "Yugoslavian Shepherd Dog -Sharplanina". The breed is currently classified to the FCI group 2, section 2.2., along with other two livestock guard dogs from the Western Balkans, Tornjak (Bosnian and Herzegovinian-Croatian Shepherd Dog) and Karst Shepherd Dog. While these three breeds share many morphological and behavioral similarities, they are currently genetically well differentiated [25]. According to Ceh and Dovc [25], the official separation of these closely related dog breeds mainly accounts for their genetic differentiation which, however, is affected also by the stochastic events in small populations. The information regarding the origin and the size of the YSD foundation stock is still lacking. In the study of YSD blood proteins, Dimitrijevic et al. [26] hypothesized that the possible ancestors of this breed are Tibetan Mastiff, ancient Greek moloscian and Roman dogs, and Turkey livestock guard dogs. On the other hand, Ceh and Dovc [25] found that YSD and Caucasian Shepherd could not be distinguished at the nuclear DNA level with the used set of molecular markers. Furthermore, the authors observed that mitochondrial (mtDNA) haplotypes found in 12 individuals of the Caucasian Shepherd represented a subset of those recorded in 129 YSD dogs. Therefore, the authors hypothesized periodic past connection between Balkan breeds and Caucasian livestock guard dogs or, alternatively, their recent admixture. In the case of the YSD, such admixture was likely especially after the Second World War when this breed became rare. Furthermore, the study reported that YSD, in comparison to the Karst Shepherd and Tornjak, is characterized by the highest levels of genetic diversity, allelic richness and the number of private alleles at the nuclear DNA level, and also by the highest number of mtDNA haplotypes. The observed haplotypic diversity, along with a finding of a unique mtDNA haplotype which has not been reported previously in other dog breeds [25], pointed towards a rather diverse maternal ancestry of the YSD breed. The findings of Ceh and Dovc [25] on rather high levels of gene diversity in YSD assessed with nuclear microsatellites (H E =0.770, 109 studied individuals) are concordant with those reported previously by Dimitrijevic [27]. Dimitrijevic [27] analyzed 88 unrelated YSD individuals with a different set of nuclear microsatellites in comparison to that used by Ceh and Dovc [25] and found H E of 0.702 in this breed. In addition, the author analyzed 15 parent-offspring dogs in order to assess whether these molecular markers may be used for parentage verification and individual identification in this breed [28]. However, while possible substructure of the analysed YSD sample was not addressed in the study of Dimitrijevic [27], Ceh and Dovc [25] noticed that nuclear genetic profiles of a few individuals deviate from those typical for the majority of dogs of this breed, implying possible breed subdivision. To the author's best knowledge, a substructure of the YSD breed has not been reported in other available studies. The aim of our work is to contribute towards better understanding of the genetic make-up of the YSD breed by assessing whether implications on its substructure and possible subdivision of its foundation stock are supported. To this end, we took advantage of the dataset of Dimitrijevic [27] comprising nuclear microsatellite genetic profiles of 94 YSD breed individuals, sampled at three dog shows and at a military training centre for dogs in Serbia in 2006 and 2007, to study genetic diversity and genetic differentiation in this population.

Study sample and nuclear microsatellite genetic profiles
We took advantage of the dataset of Dimitrijevic [27] with genetic profiles of 88 unrelated YSD individuals assessed with 10 nuclear microsatellites, namely PEZ01, PEZ03, PEZ05, PEZ06, PEZ08, PEZ12, PEZ20, FHC 2010, FHC 2054 and FHC 2079, included in the StockMarks® for Dogs Canine Genotyping Kit (Applied Biosystems, Foster City, CA, USA). The locus PEZ03 was excluded from further analyses in our study, because Dimitrijevic [27] noticed an excess of homozygotes only at this locus, which may be due to the amplification of multiple loci. Furthermore, Völkel [29] reported inconsistent and unreliable PCR amplification of this locus in 375 individuals belonging to 14 breeds. We included into analyses six additional individuals, which had been genotyped at the same time with the same panel, but were excluded from the study of Dimitrijevic [27] due to the missing data at certain microsatellite loci. Therefore, the overall sample used in our study comprised 94 registered YSD individuals, 63 males and 31 females, with age ranging from one to 12 years (4.7 years in average) ( Table 1). It is worth mentioning that animals were considered unrelated if they did not share a common ancestor for at least two generations. Two pairs of fullsibs were sampled as well, because they represented an offspring of popular parents. Native locally adapted YSD dogs outside the registry system, present in individual households in the rural parts of the mountainous Western Balkans, were not used.
As explained by Dimitrijevic [27], samples (buccal swabs) were taken in 2006 and 2007 from dogs that fit breed standards at dog shows specialized for the YSD breed  Table 1). Procedures for the extraction of the genomic DNA, PCR amplification of microsatellite loci, separation of PCR products and allele scoring are given in Dimitrijevic [27] and Dimitrijevic et al. [28].   All owners gave consent to use samples for our population genetics study. The samples were collected in accordance with institutional and European directives for research animal use. The Ethics Commission of the Faculty of Veterinary Medicine, University of Belgrade, waived the need for ethics approval and the need to obtain consent for the collection, analysis and publication of the data for this non interventional study which used a noninvasive sampling technique.

data analyses
The number of individuals used for calculating standard genetic diversity parameters was 92, because we excluded from these calculations one individual from each of the two full-sib pairs, while all individuals were used for all other analyses. Standard genetic diversity parameters, namely, the number of alleles (A), the effective number of alleles (Ae), observed heterozygosity (H O ), expected heterozygosity (H E ) and fixation index (F), were assessed using GenAlEx 6.5 [30]. Polymorphic Information Content (PIC) for individual microsatellite loci was calculated according to Botstein et al. [31] using MolKin v2.0 [32]. Deviations from Hardy-Weinberg expectations and pairwise linkage disequilibrium among loci were assessed with Arlequin 3.5 [33].
We used Bayesian clustering algorithm implemented in STRUCTURE 2.3.4. [39] to study possible genetic differentiation of the investigated population. Admixture model and allele frequency correlated model were employed, and burn-in and the analysis set to 500,000 and 100,000 iterations, respectively. Ten independent runs for the number of genetic groups (K) set from 1 to 5 were performed. The most likely K value was determined using the ΔK method of Evanno et al. [40]. Then, we used the matrix of pairwise Nei's Da genetic distance [34] based on allele frequencies for constructing a Neighbour-Joining (NJ) tree using the Populations software (Olivier Langella: Populations 1.2.32, https://bioinformatics.org/populations/, accessed 7 th December 2019). The NJ tree was visualised with FigTree 1.4.4 [35], and the statistical support for branches was assessed based on 100 bootstrap (BS) replicates, with BS values >75% considered as good support, and those >50% and <74% as moderate support. To further study the relationship among individuals, we computed pairwise genotypic distances among individuals and used them for principal coordinates analysis (PCoA) with GenAlEx. The obtained genotypic distances were also summarized by two-dimensional scaling (MDS) using PAST ver. 3.0 [36], and minimal spanning tree, which is the shortest possible set of lines connecting all points, was constructed. Minimal spanning tree analysis [37,38] investigates the spatial distribution of a point pattern with focus on small scales, and it is comparable to nearest neighbour analysis but with somewhat different properties [36].
We also calculated mean relatedness (r) according to Wang [41] for the entire population and two groups obtained in STRUCTURE analysis at K=2 as a measure of co-ancestry among individuals that may also serve as a surrogate to pedigree data [42][43][44]. Calculations were performed with COANCESTRY software [45].

Genetic diversity
The percentage of missing data in overall dataset was 1.06%. It is worth mentioning that alleles outside of the range of lengths defined for microsatellites comprising the StockMarks® for Dogs Canine Genotyping Kit were observed at two loci, PEZ05 (two alleles, lengths of 125 and 133 base pairs) and PEZ06 (length of 222 base pairs).
The number of alleles per locus ranged from five (FHC2010 and FHC2079) to 11 (PEZ12) with a mean value of 7±0.667 alleles per locus (  Table 2). Fixation index was negative only at locus FHC2054, but deviations from Hardy-Weinberg expectations at individual loci were not statistically significant at the 95% level. However, low but significant excess of homozygotes was observed in the overall sample (F=0.041±0.010). Linkage disequilibrium was significant at the 95% level among nine out of 36 tested pairs of loci (25%), which is well above the 5% expected to occur by chance.

Genetic differentiation
Bayesian clustering analysis revealed the presence of two gene pools, named green gene pool (GGP) and red gene pool (RGP) (Figure 2a, b). At K=2, 40 individuals were assigned to the GGP (assignment probability (q) >0.50), and 54 individuals were assigned to the RGP with q >0.50. Interestingly, at K=3, 28 individuals were assigned to the GGP with q>0.50, and the remaining individuals were characterized by similar admixed genetic profiles comprising the green gene pool in a minor proportion. Although individuals assigned to the second and the third gene pool both comprised red and a novel (blue) gene pool, the red gene pool was slightly more represented in the second gene pool, and the blue gene pool in the third gene pool. This finding indicates not only substructure of the studied population, but also different integrity of two observed gene pools. Both gene pools were almost equally represented at all four sampling sites (Table 1). In order to assess whether genetic differentiation of the YSD gene pool obtained with the Bayesian clustering may be detected with the Wright's F statistics, we assembled two groups, each comprising individuals assigned to one of the two gene pools with q>0.70 as inferred from the analysis at K=2 Nei's Da genetic distance among pairs of individuals ranged from 0.1436 (SPP_88 and SP_501, one of the two full-sibs pairs) to 0.9444 (SPP_79 and SP_147) (data not shown). In the NJ tree constructed using the matrix of pairwise Nei's Da, only eight moderately supported tip clades, each comprising two individuals, were observed, while deeper nodes were not supported (Figure 3). Two of these clades comprised two pairs of full-sibs. Nonetheless, it is noteworthy that all individuals were grouped into two not supported clades, one predominantly comprising individuals strongly assigned to the GGP obtained in STRUCTURE analysis, and the other comprising mainly individuals strongly assigned to the RGP. The scores of individuals in the PCoA graph were dispersed within the two-dimensional space defined by principal coordinates 1 and 2 that explained 9.95% and 7.84% of the total variability, respectively ( Figure 4). Individuals strongly assigned to the GGP were mainly grouped within the upper left quadrant, and those belonging to the RGP were mainly found in the lower right and left quadrant (Figure 4). The scores of individuals in the MDS graph were also dispersed, but, in this analysis, individuals belonging to the two gene pools were mainly separated along the Coordinate 2 ( Figure 5). Several lineages comprising groups of more closely related individuals were observed in the MDS graph in which minimal spanning tree was presented. Individuals from each of the two full-sib pars were assigned to the same GGP in analysis at K=2. However, the proportion of the GGP in individuals from these pairs was different, i.e., it was high in one pair (q>0.80, SPP_88 and SP_501) and lower in the second pair (q<0.80, SP_353 and SP_354). Increased proportion of the admixed RGP in individuals from the second pair may account for the positioning of scores of these individuals in the PCoA and MDS graph, which were more distant that scores of individuals from the first pair. The distribution of pairwise relatedness variables was centred near 0 ( Figure 6), and the mean relatedness in the studied YSD population was -0.066 (variance 0.044). The mean relatedness in two YSD subpopulations corresponding to the GGP and RGP were 0.052 (variance 0.055) and -0.0001 (variance 0.049), respectively. The relatedness of individuals from the two pairs of full-sibs was slightly higher (0.65) in comparison  to the expected value of 0.5 [45], and furthermore, r of c. 0.5 was recorded in several other pairwise comparisons suggesting that some other closely related individuals were present in our sample despite our attempt to sample unrelated individuals.

dIscussIon
Although recognized as a distinct group 2 dog breed by the FCI only 60 years ago, YSD is a well-known dog breed in the Western Balkans present in this region for a long time, and highly valued by both local people and professional dog breeders. It is genetically well distinguished from other livestock guard dogs from the Western Balkans, Karst Shepherd and Tornjak [25]. We report levels of genetic diversity in YSD which are in the range of those reported by Dimitrijević [2008] and Ceh and Dovc [25], and are comparable and in some cases even higher than those found in numerous dog breeds studied to date by means of nuclear microsatellites [4,[11][12][13][14][15][16][17][18][19][20][21][22][23]. Thus, despite the decline of the YSD population after both 20 th century World Wars, this breed remains genetically highly diverse. The finding of three alleles at two nuclear microsatellite loci in YSD which are outside the expected length range was not fully unexpected because a similar phenomenon has also been observed by Völkel [29]. Assuming that these alleles were not due to the genotyping errors and have not emerged due to de novo mutations, their occurrence in the YSD may be explained by past hybridization with wolves which has been hypothesized previously [26]. Admixture between dogs and wolves is not uncommon [46][47][48][49][50], it is more frequently observed in ancient breeds [51], and it has also been recorded between shepherd dogs and wolves in Georgia, Caucasus [52]. It is worth mentioning that although transferable to closely related species, microsatellite loci are likely to display differences in polymorphism level and/or to harbour different sets Figure 6. Distribution of frequencies of pairwise relatedness according to Wang (2011) of alleles in different species [53]. Such alleles as well those that are rare but expected in dog breeds [14] and those whose recent emergence has been documented (e.g., somatic mutation at locus PEZ20 in one Rottweiler individual, Pádár et al. [54]) are particularly valuable in forensics [54]. Given the occurrence of such alleles in YSD individuals strongly assigned to each of the two observed gene pools, it is very likely that they were introduced into the foundation stock of this breed rather early. On the other hand, past hybridization with wolves may also account for the occurrence of a unique mtDNA haplotype in YSD which has not been reported previously in other dog breeds [25]. The occurrence of wolf-specific alleles in YSD due to past hybridisation, however, may be confirmed by further comparative analyses of both YSD and wolfs.
It has been demonstrated previously that inbreeding and loss of genetic diversity are common in certain purebred dog populations due to the specific breeding practice in dogs [5][6][7][8], and that they are associated with an increased incidence of a number of inherited diseases [9][10][11]. In the studied YSD population, fixation index was low, positive and significant, but a significant excess of homozygotes was not recorded at individual microsatellite loci. On the other hand, linkage disequilibrium among microsatellite loci was observed in more pairwise comparisons than expected simply by chance (25% vs. 5%). This outcome was unexpected because loci used in our study, which comprise a standard kit for canine genotyping, are not supposed to be linked in unrelated individuals, as observed in studies in which this kit was used [13,14]. Altogether, an excess of homozygotes and detection of linkage disequilibrium among non-linked loci are indicative of population substructure (i.e., Wahlund effect), which has been recorded in STRUCTURE analysis and implied also in the NJ tree, PCoA and MDS analyses. A possible YSD population subdivision has also been observed in the study of Ceh and Dovc [25].
An almost equal number of individuals sampled in 2006 and 2007 at four different sites was assigned to each of the two gene pools, and furthermore, both gene pools were almost equally represented at each sampling site (Table 1), suggesting the lack of preferences of dog owners towards individuals belonging to a particular gene pool. This is possibly due to the lack of marked morphological and/or other differences among them. However, in the study of Ceh and Dovc [25], only a handful out of 109 studied YSD individuals were distinct at the nuclear DNA level. Unfortunately, since different sets of nuclear microsatellites were used in our and the study of Ceh and Dovc [25], we were unable to infer which out of the two gene pools recorded in our work was prevalent in the sample studied by Ceh and Dovc [25]. Nonetheless, since the two pairs of full-sibs used in our study are the offspring of popular parents assigned to the GGP with relatively high assignment probability, it is tempting to hypothesize that the GGP is expanding in the YSD population due to the specific breeding practice and usage of popular individuals for mating. In this case, genetically distinct YSD individuals observed by Ceh and Dovc [25] may correspond to our GGP.
We found that two groups assembled from individuals assigned to each of the two gene pools with q>0.70 were genetically differentiated (F ST =0.134, P<0.05). In addition, GGP comprised individuals with somewhat lower genetic diversity (H E =0.599±0.064), more homogeneous genetic profiles (as inferred from a more limited dispersal of scores of individuals strongly assigned to this gene pool in PCoA and MDS graphs, Figures 4  and 5), but with pronounced genetic integrity that was kept in STRUCTURE analysis at K=3 (Figure 2). On the other hand, RGP was characterized by somewhat higher levels of genetic diversity (H E =0.749±0.019) and more heterogeneous (admixed) ancestry as inferred from the wider dispersal of scores of individuals strongly assigned to this gene pool in PCoA and MDS graphs (Figures 4 and 5). This finding may be explained by a large and possibly non-uniform foundation stock of the YSD, comprising not only individuals that fit morphologically and in other features to the breed standards but also those whose conformational and functional characteristics were improved through mating with dogs from native locally adapted populations or with other breeds such as, for instance, Caucasian Shepherd. This is possible, because unexplored native breeding stock of the YSD outside the registry system is still present, and because mating with the Caucasian Shepherd was indeed inferred in the study of Ceh and Dovc [25]. As already mentioned, past hybridization between YSD and wolves has been observed as well [26]. Thus, along with the previously reported diverse maternal ancestry of this breed [25], our data support rather complex and diverse ancestry of the YSD in general, and very likely genetically heterogeneous foundation stock of this breed.
In conclusion, given the lack of knowledge on the origin and the size of the YSD foundation stock, our work represents an important contribution in this respect. We found that the studied YSD population of 94 individuals sampled in 2006 and 2007 in Serbia is genetically differentiated despite the fact that all sampled dogs fit well to the conformational and functional characteristics of the breed. Thus, genetically heterogeneous foundation stock of this breed has been assumed, as well as increased incidence of hybridization of individuals belonging to one of the two observed gene pools, possibly with native landrace breeding stock which is still present in rural parts of the Western Balkans or with other breeds such as the Caucasian Shepherd. In addition, past admixture of the YSD with wolves cannot be excluded. Overall, we point to the rather complex and diverse ancestry of the YSD in general, which is in accordance with previously hypothesized diverse maternal ancestry of this breed [25]. We also support previous findings on rather high levels of nuclear genetic diversity in YSD [27,25] and demonstrate that, contrary to some purebred dog populations that also experienced population decline in the past, this breed does not suffer from inbreeding.