High-Density Genomic Characterization of Native Croatian Sheep Breeds

A recent comprehensive genomic analysis based on 50K SNP profiles has shown that the regional Balkan sheep populations have considerable genetic overlap but are distinctly different from surrounding breeds. All eight Croatian sheep breeds were represented by a small number of individuals per breed. Here, we genotyped 220 individuals representing the native Croatian sheep breeds (Istrian Sheep, Krk Island Sheep, Cres Island Sheep, Rab Island Sheep, Lika Pramenka, Pag Island Sheep, Dalmatian Pramenka, Dubrovnik Sheep) and mouflon using the Ovine Infinium® HD SNP BeadChip (606,006 SNPs). In addition, we included publicly available Balkan Pramenka and other Mediterranean sheep breeds. Our analyses revealed the complex population structure of Croatian sheep breeds and their origin and geographic barriers (island versus mainland). Migration patterns confirmed the historical establishment of breeds and the pathways of gene flow. Inbreeding coefficients (FROH>2 Mb) between sheep populations ranged from 0.025 to 0.070, with lower inbreeding coefficients observed in Dalmatian Pramenka and Pag Island Sheep and higher inbreeding in Dubrovnik sheep. The estimated effective population size ranged from 61 to 1039 for Krk Island Sheep and Dalmatian Pramenka, respectively. Higher inbreeding levels and lower effective population size indicate the need for improved conservation management to maintain genetic diversity in some breeds. Our results will contribute to breeding and conservation strategies of native Croatian sheep breeds.


INTRODUCTION
Sheep were domesticated about 10,000 years ago in the region of Anatolia and, along with goats, were among the first domesticated farm animals. Sheep were first hunted by humans and over time became managed wild populations, then kept in controlled herds, and finally, humans began to breed sheep (Zeder, 2008;Larson and Burger, 2013). Trought the history breeding for diseriable traits was always present but only about 200 years ago, the organized formation of sheep breeds began, resulting in fragmented populations, and reduced genetic diversity (Taberlet et al., 2008). Genetic diversity is defined as the variety of alleles and genotypes present in a population (Frankham et al., 2002). Each breed has a unique genetic characteristic due to mutations and drift caused by geographic isolation and bottlenecks, artificial selection and adaptation to climate, nutrition and diseases and parasites (Barker, 2001). Many sheep breeds are local breeds that have adapted to specific locations for thousands of years and are closely associated with culture and history (Gandini and Villa, 2003). These small, unique breeds contribute to the overall diversity of the species, but diversity within breeds is low. Highly productive breeds have displaced local breeds and unique gene combinations of some local breeds are being lost (Groeneveld et al., 2010). Erosion of within-breed diversity could be a problem even for breeds whose overall populations remain very large. Monitoring population trends is a prerequisite for rapid and effective action to protect breeds from extinction. Measures to prevent the loss of livestock diversity will be more effective if the factors that drive genetic erosion and extinction risk are well understood. Maintaining genetic diversity is important for rapid adaptation to challenges (Andersson, 2012). Thus, characterizing genetic diversity is an important aspect of developing sustainable strategies for breed improvement (Groeneveld et al., 2010).
Livestock breeding, especially sheep species, was an important part of agriculture in Croatia. The first written records of sheep breeding in Croatia date back to 1781, and two decades later, the number of sheep per citizen in Europe was highest in Dalmatia (Defilipis, 1966), a region of Croatia whose name comes from Illyrian words dalma, meaning sheep. Since then, the number of sheep on Croatian territory has steadily decreased. The erosion of agricultural production combined with wars led to a drastic decrease in the number of some native Croatian breeds, which were threatened by extinction. Maintaining genetic diversity in livestock production is critical for meeting future challenges such as climate change, emerging diseases, and food security for a growing human population (Barker, 2001;Groeneveld et al., 2010). Croatia has a number of unique sheep breeds. The native Croatian sheep breeds belong to the Pramenka type, which is characterized by coarse wool, low production, and high resistance to environmental conditions. Within the Pramenka type, there are many breeds/populations that differ in size, wool quality and colour, as well as in their adaptation to the specific microclimatic conditions in their breeding area. Throughout history, Croatian farmers have been eager to improve production characteristics and therefore imported high-performance breeds. Merino rams (Pag Island Sheep, Rab Island Sheep, Krk Island Sheep, Cres Island Sheep and Dubrovnik sheep) and meat breeds such as Merinolandschaf (Dalmatian Pramenka and Lika Pramenka) or Bergamasca Sheep (Istrian Sheep) were mainly used to improve primitive Pramenka type. With the progress of selection, different breeds were created. Genetic diversity in sheep needs to be estimated to identify unique breeds that may be in danger of extinction.
The genetic study of native Croatian breeds started with βlactoglobulin polymorphisms in Pag Island Sheep (Cubric-Curik et al., 2002). Since then, scientists have studied genetic diversity (Feligini et al., 2005;Ivanković et al., 2005;Cubric-Curik et al., 2009;Kasap et al., 2021;Špehar et al., 2022) and population structure (Lawson Handley et al., 2007;Ćinkulov et al., 2007;Ferencakovic et al., 2012;Salamon et al., 2014) of native Croatian sheep breeds using different approaches. Recently, 50 K SNP markers (Ciani et al., 2020) and whole-genome sequences (Deng et al., 2020;Lv et al., 2021) have been used in population genetic studies presenting native Croatian breeds, but with a very small number of animals per breed. Ciani et al. (2020) showed the separation of Balkan Pramenka from the rest of the European breeds and the belonging of native Croatian sheep breeds to the Pramenka cluster. The authors suggested that Balkan breeds are evolutionary connected with the domestication process and are one of the hub regions from which the migration of sheep spread to the rest of Europe. Recently, Shihabi et al. (2022) identified specific adaptive selection signals on the X chromosome of the individuals used in this analysis.
The main objective of this study was to evaluate the conservation status (diversity, inbreeding, and effective population size), population structure, and admixture of eight native Croatian sheep breeds. Our analysis was based on highdensity autosomal SNP chip genotype information sampled from 20 to 45 individuals per breed. Our study is an extension of the study by Ciani et al. (2020), in which only a few individuals representing several native Croatian breeds were analyzed using 50 K SNP chip information. More specifically, the data used refer to a sample of 201 individuals representing eight native Croatian breeds (Cres Island Sheep, Dalmatian Pramenka, Dubrovnik Ruda, Istrian Sheep, Krk Island Sheep, Lika Pramenka, Pag Island Sheep, and Rab Island Sheep) and 10 mouflon sampled in Croatia. They provide an accurate estimate of genetic diversity, ROH-based inbreeding, current and historical effective population size, and allow in-depth comparison with other Mediterranean breeds and some other Pramenka sheep breeds.

Sample Collection
France) and stored at +4°C until further use. All samples were collected according to national and European ethical protocols and directives.

DNA Extraction and SNP Genotyping
DNA extraction was performed using the Qiagen DNeasy Blood & Tissue Kit (Qiagen, Germany) according to the sample preparation and extraction protocol. DNA quality was checked by 1% agarose gel electrophoresis, while DNA quantity was determined using a NanoPhotometer P330 spectrophotometer (IMPLEN, Germany).
Croatian sheep and mouflon samples were genotyped using the Ovine Infinium HD BeadChip (Illumina, San Diego, CA, United States) by Gene Seek, Neogen Genomics (Neogene Europe Ltd., Scotland, United Kingdom). Genotypes were mapped to the Oar4.0 map. Quality control and filtering of SNPs were obtained using SNP & Variation Suite v8.7.0 (Golden Helix, Inc., Bozeman, MT, www.goldenhelix.com). The accuracy and efficiency of SNP genotyping were assessed by applying a cut-off value of 0.7 for the Illumina GenCall score and 0.4 for the Illumina GenTrain score. All SNPs with more than 10% missing genotypes and individuals with more than 5% missing genotypes were removed. SNPs without position or SNPs that had been assigned to sex chromosomes or mitochondrial genome were also excluded.  Figure 2A.

Population Divergence and Relationship
To estimate genetic diversity within the population, expected heterozygosity, observed heterozygosity and the inbreeding coefficient (F IS ) were calculated using the software PLINK 1.9 (Purcell et al., 2007). In addition, Wright's F ST fixation indices were calculated to determine the degree of genetic differentiation among selected sheep breeds using SNP & Variation Suite v8.7.0 (Golden Helix, Inc., Bozeman, MT, www.goldenhelix.com) and mean F ST was calculated. Principal Component Analysis (PCA) was performed to investigate the genetic relationship among Croatian sheep breeds as well as at the Mediterranean level. PCA was performed for both datasets: Mediterranean and Croatian, using the R package "SNPRelate" (Zheng et al., 2012). Phylogenetic relationships based on Nei's genetic distances (Nei, 1978) were presented with a Neighbor network using SplitsTree4 software (Huson and Bryant, 2006). Nei genetic distances were calculated in the program R using the package "stAMPP" (Pembleton et al., 2013).

Population Structure
Genetic relationships and population structures were assessed with the software STRUCTURE v.2.3.4. (Pritchard et al., 2000) FIGURE 1 | Map of sampling locations for native Croatian sheep breeds. The map illustrates the geographic locations where samples were collected. Each breed is presented with a different colour.
Frontiers in Genetics | www.frontiersin.org July 2022 | Volume 13 | Article 940736 using the Bayesian model based on 83,795 SNPs in eight Croatian sheep breeds and 104,623 SNPs in 28 Mediterranean sheep breeds after pruning based on LD. Pruning of the data was performed in a sliding window of 50 SNPs, moving in increments of five SNPs along each chromosome and removing one SNP pair at a time with a pairwise r 2 > 0.1 using the software PLINK 1.9 (Purcell et al., 2007). After pruning the SNP missing rate was 0.000-0.003 and 0.000-0.031 per individual for Croatian dataset and Mediterannean dataset, respectively, and 0.003-0.011 for Mourerous and Altamurana population, respectively. For the Croatian and Mediterranean datasets, an ancestry model was constructed for a putative number of ancestral populations of 1-10 and 1-30, respectively. For each K value, 10 runs of 100,000 Markov chain Monte Carlo iterations were performed after a burn-in period of 10,000 iterations. The most likely number of clusters was determined using the ΔK method (Evanno et al., 2005) implemented in software Structure Selector (Li and Liu, 2018). Results from STRUCTURE were analyzed and visualized using the software Structure Selector (Li and Liu, 2018).

Gene Flow
Gene flow between Croatian and Mediterranean sheep breeds was examined using TreeMix v. 1.13 (Pickrell and Pritchard, 2012). TreeMix was applied to infer population trees and migration events between divergent populations. To determine the position of the root in the maximum likelihood tree, the European mouflon (n = 10) was defined as the outgroup population. Between 0 and 15 migration events were inferred.

Genomic Inbreeding Based on Runs of Homozygosity
The proportion of the genome that is autozygose is estimated by identification of Runs Of Homozygosity (ROH) with SNP & Variation Suite v8.7.0. ROH was detected for each individual separately using the following criteria: (I) the minimum number of SNPs in ROH was set to 50, which was calculated using the method proposed by Lencz et al. (2007) to minimize the number of false-positive ROH: where α is the percentage of false-positive ROH (set to 0.05 in this study), n s is the number of SNPs per individual, n i is the number of individuals, het is heterozygosity across all SNPs; (II) the maximal gap between adjacent SNPs was set to 250 Kb; (III) the minimum SNP density per ROH was set to 1 SNP every 50 Kb; (IV) the minimum length constituting ROH was set to 2 Mb. To account for genotyping error, ROHs were calculated separately for each of the five categories specified according to the ROH length: >2 Mb, 2-4 Mb, 4-8 Mb, 8-16 Mb, and >16 Mb. Based on ROH length, the number of allowed heterozygotes and missing genotypes were defined according to Ferenčaković et al. (2013). The individual genomic inbreeding coefficients based on ROH (F ROH ) were estimated as the sum of the length of all ROH per individual divided by the total length of the autosomal genome covered by SNPs, as described by McQuillan et al. (2008). The F ROH was calculated for different minimum ROH lengths because the short segments are associated with individual autozygosity derived from a distant common ancestor, whereas long ROH segments are correlated with recent inbreeding (Ferenčaković et al., 2013). For each animal, we calculated the F ROH2-4 Mb , F ROH4-8 Mb , F ROH8-16 Mb , and F ROH>16 Mb based on ROH with different minimum lengths. Based on the estimated individual F ROH>2 Mb , the breed inbreeding coefficient was calculated by averaging the F ROH estimates per breed.

Effective Population Size
To better understand the demographic history of native Croatian sheep breeds, we estimated current and historical effective population size (N e ) form linkage disequilibrium (LD). Trends in effective population size were estimated using an optimization method that implements a genetic algorithm (Mitchell, 1998) to infer the recent demographic history of a population from SNP data of a small sample of contemporary individuals, as implemented in the software GONE (Santiago et al., 2020). Default parameters were applied. The most recent estimate of N e was taken in the current generation. Furthermore, N e estimates from 50 generations ago were used for comparison with results from other authors (Kijas et al., 2012;Beynon et al., 2015;Deniskova et al., 2019).

Population Divergence and Relationships
Genetic diversity parameters are presented in Table 1. Observed heterozygosity ranged from 0.327 ± 0.006 in Rab Island Sheep to 0.346 ± 0.003 in Dalmatian Pramenka. The average expected heterozygosity within the breeds ranged from 0.336 in Istrian Sheep to 0.348 in Dalmatian Pramenka. The inbreeding coefficient (F IS ) varied from −0.017 ± 0.011 in Istrian Sheep to 0.022 ± 0.018 in Rab Island Sheep. The results of observed and expected heterozygosity are similar to those found by Al-Mamun et al. (2015) in five Australian sheep breeds using a 50k SNP chip and Edea et al. (2017) in Ethiopian sheep using a 50 and 600K SNP chip. Similar values of expected heterozygosity were reported for native Italian breeds (Ciani et al., 2014), southern and western European sheep breeds (Kijas et al., 2012) and native Russian sheep breeds (Deniskova et al., 2018).
Genetic relationships between individuals of multiple sheep breeds were assessed by principal component analysis (PCA). The PCA plot of the first and second principal components (PC) is shown in Figure 2B. The breeds were generally grouped according to their geographic origin. Similar results were also presented in Rochus et al. (2020) and Ciani et al. (2020). Total variation explained by the PC components was 94.52%. The largest component (8.6% of the total variation) showed an east-west cline between the Western European sheep and the Balkan Pramenka sheep. The second PC separated Croatian island breeds (ISS, CIS, RIS, KIS, PIS) from the rest of the Pramenka cluster and highlighted a division of native Croatian breeds into mainland breeds (DRS, LPS, and DPS) and island breeds. Two Croatian breed groups have contrasting structures. In the mainland group, breeds were clearly separated and individuals had clear cluster assignments, indicating greater differentiation between breeds. In contrast, the island breeds showed less differentiation with tighter clusters and breed overlapping. A similar differentiation in the structure of the two French sheep groups was observed in Rochus et al. (2018). Two Italian breeds occupy a central position between Balkan Pramenka sheep and French and Spanish sheep. The genetic proximity of Pramenka sheep, Italian and Spanish breeds is interpreted as a possible sheep migration route along the Mediterranean coast from the Balkans to Spain via Italy, as presented in Ciani et al. (2020).
To better estimate the relationships between native Croatian sheep breeds PCA was performed only on the Croatian dataset (Supplementary Figure S1). The first two components together explained 14.9% of the total variation, representing a large portion of the variability. Principal component analysis showed a clear separation of Istrian Sheep and Dubrovnik Sheep, while Lika Pramenka and Dalmatian Pramenka created one close cluster and Rab Island Sheep, Pag Island Sheep, Krk Island Sheep, and Cres Island Sheep formed a second cluster. Mainland and island cluster are not well defined here. The result shows a north-south cline following the geographical breeding area of each breed. A similar pattern following the geographical distribution of breeds in the PCA graph was found in Ciani et al. (2014) and Edea et al. (2017). The largest PC (8.1% of the total variation) positioned Dubrovnik Sheep apart from the rest of the Croatian native sheep breeds. The second PC (6.8% of the total variation) separated the Istrian Sheep from all others. The genetic uniqueness of Istrian Sheep was observed earlier in Ćinkulov et al. (2007), Lawson Handley et al. (2007), and Salamon et al. (2014) using microsatellite markers. Istrian Sheep and Dubrovnik Sheep are also geographically isolated from other native Croatian sheep breeds, but they are also phenotypically the most diverse.
The relationship between the Mediterranean sheep breeds and the mouflon was determined using genetic distances according to Nei (1972), which include a correction for drift and mutations. Higher genetic distances also indicate a longer temporal divergence of the breeds. The European mouflon was set as an outgroup. Neighbor-net analysis of Mediterranean sheep ( Figure 2C) confirms the separation of the Pramenka from the rest of the European sheep population. The genetic uniqueness of the Croatian island breeds is emphasized here by their separation from the Pramenka group. Analysis by Ciani et al. (2020) using a 50K SNP chip shows similar separation of Croatian island sheep breeds from other European breeds and their positioning between Pramenka and Italian sheep. The Lika and Dalmatian Pramenka were in a group with the Serbian Pramenka and the Ukrainian Mountain Carpathian sheep, confirming the relationships identified by Ciani et al. (2020). The separation of Croatian island and continental breeds suggests that geographic barriers have an influence on shaping the current structure of native Croatian sheep breeds. Greater geographic distances with barriers such as the sea and mountains have reduced the exchange of genetic material between these two groups of breeds, resulting in greater separation of these breeds. Earlier gene flows and population mixing, which most likely occurred between the Balkan Pramenka populations, are still visible in the genetic structure of the Croatian mainland breeds and should not be neglected here. In the past, a Pramenka was bred on the territory of today's Republic of Croatia, the socalled Balkan Pramenka, which was divided into different types depending on where it was bred (e.g., Travnik type, Sjenica type . . .). Later, with the disintegration of Yugoslavia, the gene flow between Pramenka types decreased and new breeds were formed, which still show a great genetic connection. The more frequent and stronger introduction of Italian breeds into island populations, especially Pag Island Sheep, is reflected in the proximity of these two sheep groups in the graph, although in this dataset were no representatives of Italian Merino sheep involved in breeding native Croatian sheep breeds.
The level of genetic differentiation based on F ST estimation between breeds were calculated. Genetically similar populations will have lower F ST

Population Structure
The population structure of Mediterranean breeds was identified using a model-based STRUCTURE analysis with an assumed number of populations set to K = 30. The Evanno ΔK method did not reveal a clear number of genetic clusters, peaking at K = 4, K = 6, and K = 11 (Supplementary Figure S2), indicating a hierarchical population structure. The results of model-based population structuring revealed a complex genetic structure of native Croatian sheep breeds when placed in the Mediterranean context ( Figure 2D). At K = 2, we observed the separation of Balkan Pramenka from Western European breeds. A similar split is presented in the genetic diversity results of 57 European and Middle Eastern sheep breeds by Peter et al. (2007) To identify population structure within native Croatian sheep breeds on a finer scale, we performed another STRUCTURE analysis only on the Croatian dataset with an assumed number of populations (K) between one and ten. The most informative number of ancestral populations (K = 3) was estimated using the Evanno ΔK method (Supplementary Figure S3). The results for K = 2-6 are presented in Supplementary Figure S4. The clearest ancestry components when assuming only two ancestral populations showed a clear separation of the Dubrovnik Sheep. At K = 4, the Istrian Sheep are completely separated. The specific genetic structure of the Istrian Sheep was also shown in the research of Ciani et al. (2014), where it stands out as a distinct breed at a very low K = 8. At K = 5, Dalmatian Pramenka and Lika Pramenka form a clearly separate group. These two breeds are not separated even at K = 10, the maximum number of hypothesized groups for this data set that suggest common ancestry. At K = 6, the Krk Island Sheep is divided into three subpopulations, with two subpopulations showing the dominance of a single ancestral genome, while the third population is genetically similar to other island breeds. The cause of such separation within the Krk Island Sheep breed could be due to reduced exchange of genetic material within Krk Island Sheep populations, when under the influence of inbreeding and genetic drift, gene fixation occurs and population structure changes.

Gene Flow
Migration patterns between native Croatian and Mediterranean sheep breeds were studied using population tree models accounting for migration events and implemented in the program TreeMix. The maximum likelihood tree, created without assuming migration events, was calculated and rooted in the European Mouflon (EMC) (Supplementary Figure S5). The ML tree generally confirmed the relationships revealed by the  PCA and STRUCTURE analysis, as in the Mediterranean context the Croatian native sheep breeds were clustered with the other Balkan Pramenka sheep. The Croatian island breeds additionally separated from the other Balkan Pramenka breeds, creating a new branch. As in the PCA graph, the Italian breeds together with the Corse sheep represent the connection between the Balkan and West-European breeds. The West-European sheep breeds show a similar phylogenetic structure as shown in Rochus et al. (2018). In the population tree, the sheep breeds, including the native Croatian breeds, generally had short branch lengths, indicating higher genetic diversity within the breed. By adding migration events, the population tree showed a more detailed view of the genetic relationships between the Mediterranean and Croatian breeds. Estimation of the optimal number of migration events using the ΔM method showed that the population tree with m = 7 migrations best explained the genetic relationships between the breeds (Figure 3). Inference of population trees with m = 7 showed that the migration event with the greatest weight was from the root of the Balkan Pramenka sheep breeds to the root of the two Italian sheep breeds. There are no historical records of the importation of Balkan sheep into Italy, so it is assumed that this is an ancient migration, and it would be interesting to estimate how old it is. The migration edge from the French Merino branch to the Dubrovnik Sheep branch confirms the evolutionary history of the Dubrovnik Sheep, which has a higher presence of the Merino genotype than other native Croatian sheep breeds. The Dubrovnik Sheep is also known by the name Ruda, which was given to the Dubrovnik Sheep in reference to its finer wool. The migration edge with low weight from European mouflon to Corse sheep is presented. This migration is explained by the knowledge that European mouflon living in European countries are imported from populations in Sardinia and Corsica (Guerrini et al., 2015) and that the introduction of mouflon into domestic sheep populations has been recorded in Sardinia and Corsica (Barbato et al., 2017). Migrations between French breeds are the result of breed formation or crossbreeding and are discussed in detail in Rochus et al. (2018).

Genomic Inbreeding
In sheep populations, estimation of inbreeding coefficient based on pedigree information has not been accurate due to incomplete pedigrees, low generations and frequent errors in recording (Ferenčaković et al., 2013), so estimation of inbreeding coefficients based on high-density SNP markers by ROHs gives us more accurate results. The overall genomic inbreeding (F ROH>2 Mb ) for Croatian native sheep breeds and three selected breeds from the Mediterranean dataset (Altamurana, Leccese, Churra and Corse sheep) was estimated as well as inbreeding for each of the defined categories (F ROH2-4 Mb , F ROH4-8 Mb , F ROH8-16 Mb , F ROH>16 Mb , Figure 4).
Inbreeding coefficients (F ROH>2 Mb ) between sheep populations ranged from 0.025 to 0.070 (Table 1), with lower inbreeding coefficients observed in Dalmatian Pramenka and Pag Island Sheep. Similar inbreeding coefficients were reported for native Russian breeds (Deniskova et al., 2019) and most Italian sheep breeds (Mastrangelo et al., 2018), while higher levels of inbreeding were observed in Barbaresca, Delle Langhe, Valle del Belice (Mastrangelo et al., 2017), Leccese (Mastrangelo et al., 2018, this study) and African sheep breeds (Dzomba et al., 2021). The low F ROH values observed for almost all native Croatian sheep breeds suggest that the animals in this study are not highly related. During recombination, long ROH segments are interrupted. Long ROH segments are associated with recent inbreeding, while short ROH segments indicate ancient inbreeding (Ferenčaković et al., 2013;Curik et al., 2014). The largest inbreeding values for F ROH2-4 Mb = 0.010 and F ROH4-8 Mb = 0.013 were observed in Istrian Sheep, indicating the presence of IBD individuals and more ancient  (Mioč et al., 2003), which represent the genetic base of the current Dubrovnik Sheep population. The highest degree of inbreeding in Dubrovnik Sheep (F ROH = 0.070) is a consequence of the low number of animals available for revitalization. The results show an increase in the number of homozygous regions from the 20th generation to the third generation in all breeds. Accordingly, there is an increase in the average inbreeding coefficient from 0.005 to 0.021, respectively from the 25th to the third generation. The increase in the recent inbreeding coefficient could be due to the extensive use of few rams within herds, as suggested by Mastrangelo et al. (2017). In Croatia, as in Sicily, livestock breeders use natural mating where rams are used in herds with closely related individuals for several years or one ram is used in several herds, what leads to increased inbreeding and consequently lower variability.

Effective Population Size
Effective population size (N e ) is an important parameter used in population genetics. It is useful for monitoring genetic diversity and explaining population trends. In addition, effective population size is an indicator of the risk of genetic erosion (Tenesa et al., 2007). It has been known for many years that monitoring effective population size is an important tool for the long-term conservation of endangered populations (Notter, 1999;Gandini et al., 2004;Biscarini et al., 2015). The current effective population size for Croatian sheep breeds estimated using the software GONE ranged from 61 to 1039 for Krk Island Sheep and Dalmatian Pramenka, respectively ( Table 1). Higher N e values were also observed for Pag Island Sheep (1005), Lika pramenka (598)  When considering historical effective population size, the highest N e across all generations was observed in Dalmatian Pramenka. Fifty generations ago, the lowest N e values were 3159 (Cres Island Sheep), 4030 (Istrian Sheep) and 4814 (Krk Island Sheep) and the highest were 51537 (Dalmatian Pramenka) and 35988 (Dubrovnik Sheep). These values are higher than those determined by Kijas et al. (2012) and Ciani et al. (2014) for the European sheep breeds and Deniskova et al. (2018) for the Russian breeds. Very high values of historical effective population size for certain breeds indicate that these breeds were very widespread on the territory of Croatia in the past and confirm the historical records of a large number of sheep (1,105,078 head) bred in the area of present Dalmatia (Defilipis, 1966). Also, a very high linear trend in the effective population size of Pag Island Sheep can be attributed to the systematic activities of breeders since 1870, when the Gregge Modella Society was established for the improvement of sheep breeding (Pavlinić, 1936). The much higher values obtained for the Croatian sheep breeds compared to other European breeds could be due to the application of different software. The presence of population substructuring into herds could also affect the estimate of genomic effective population size. In the simulations, the software GONE showed an accurate estimation of effective population size up to 200 generations. Unlike other programs that are based on LD for N e estimation, GONE very accurately reflects changes in N e over history, even when there were substantial increases or decreases in value (Wang et al., 2016;Santiago et al., 2020). All N e estimates based on LD are sensitive to population mixing because the presence of mixing in populations increases N e values. In addition, migration affects N e estimates, as the exchange of rams between flocks and the natural mating systems in sheep breeding can lead to an overestimation of N e values.

CONCLUSION
The present study is the first detailed analysis of genomic diversity and population structure of eight native Croatian sheep breeds using high-density SNP markers. PCA, Neighbour network, Maximum likelihood trees and STRUCTURE analyses consistently revealed a common clustering pattern showing a clear separation between the island and mainland breeds. The eight Croatian breeds were generally unique. Some degree of admixture was observed in the population structure, confirming common ancestry. For all Croatian breeds except Krk Island Sheep, the estimated effective population size was above 100 -a threshold above which the breed is sustainable. Inbreeding coefficients estimated from runs of homozygosity were low to moderate. We compared the native Croatian breeds with the other Mediterranean and Pramenka breeds genotyped with highdensity markers. In the global context, Croatian sheep clustered with other Pramenka type breeds and differed from West-European breeds. We observed gene flow from the mouflon population to domestic sheep and from Croatian native sheep breeds to Italian breeds. This migration is very interesting and it would be useful to determine the age of the migration as it could be confirmations of the Mediterranean migration route from the domestication center. This study contributes to a better understanding of the genetic background of Croatian native sheep breeds and provides information to support the genomic improvement of these local breeds.

DATA AVAILABILITY STATEMENT
Genotypic data of 201 animals representing eight Croatian sheep breeds are deposited and available at https://doi.org/10.5061/ dryad.pg4f4qrsn.

ETHICS STATEMENT
The animal study was reviewed and approved by the Bioetic Committee for Animal Welfare and Protection of the University of Zagreb, Faculty of Agriculture. Written informed consent for participation was not obtained from the owners because Animal samples were collected within National Gene Bank procedure in Croatia.

AUTHOR CONTRIBUTIONS
ID, VC-C, and IC conceived and designed the study; ID, VC-C, JK, M-HL, SM, EC, and BL organized sample collection, extraction and genotyping; ID and MS performed the analysis with the supervision of IC, VC-C, and JAL; ID wrote the manuscript with the supervision of VC-C and IC. All authors read and approved the final manuscript.