Genetic diversity and substructuring of the Hungarian merino sheep breed using microsatellite markers

Abstract The Hungarian Merino sheep breed (Ovis aries) is the most significant animal resource of the Hungarian sheep sector which, unfortunately, has gone through a huge reduction in number during the last decades and became endangered in 2014. A modern molecular genetic survey is now becoming more than necessary in order to characterize the within-breed genetic diversity and structure. For that reason, six Hungarian Merino flocks were genotyped in 16 microsatellite markers. In total, 144 different alleles were found and the mean values of observed and expected heterozygosity were 0.714 and 0.705, respectively, suggesting a noticeable genetic variability of the breed. The genetic differentiation of the Hungarian flocks was generally low, as reflected by the estimated total FST value (0.036), the extended pattern of admixture in Structure analysis, as well as, by the noticeable level of genetic clustering in UPGMA and FCA analyses. However, two out of the six studied flocks tended to be genetically more distant. The outcome of our study could be a starting point for a planned breeding strategy of the Hungarian Merino breed, by keeping the within-flock genetic variability in priority, as well as, by preserving the potential genetic uniqueness with close monitoring of the inbreeding.


Introduction
The oldest established sheep (Ovis aries) breed in the world, the Merino, was developed in Iberia and its name probably came from the 'Beni Merin' people (Moroccan tribe). From the eighteen century onward, Merinos were spread all over the world and are now the most numerous sheep breed and the major source of the world's wool supply. The Merino and Merinoderived sheep breeds represent around 20% of the world's sheep population [1]. The dispersal of Merino sheep in Eastern Europe (e.g., Hungary, Poland) started in the eighteen century. With the cheap wool being imported from Australia and spreading across Europe starting in the 1950s, the prices of wool and added value have collapsed. That caused the increase in the number of mutton sheep breeds. European Merino sheep have experienced a dramatic numerical decrease to the point that they are now considered endangered breeds and are no longer undergoing active genetic improvement [2].
More than 80% of sheep flocks in Hungary belong to the breed Hungarian Merino [3]. The Hungarian Merino was developed by crossing with many different Merino and Merino-derived breeds through the years. The first 300 Merino sheep were brought to Hungary from Spain for the request of Queen M aria Ter ezia in 1774. The Hungarian Merino was developed by crossing with Rambouillet Merino and German Mutton Merino. After World War II., the Hungarian Merino sheep breed was crossed with Russian Merino breeds (e.g., Groznen, Stavropolian, etc.). After the year 1960, to develop the meat production, the Hungarian Merino ewes were crossed with Merino Precoce rams from France and with German Mutton Merino rams. During the 1970s and 1980s, the breed was crossed for finer wool with Kent, Corriedale and Australian Merino as well as with Booroola Merino to develop its prolificacy. The Hungarian Merino Herd Book was established in 1993. Since that time no other Merino breeds were used for crossing with the Hungarian Merino [4]. This breed is bred for a dual purpose: its meat and wool are the products for the market, but meat production is more significant and profitable. The Hungarian Merino is well-adapted to the poor pasture and to the hot and dry climate. Due to a long history of adaptation to the local environment, production linked to local populations is often the only sustainable activity in areas of harsh environmental conditions [5].
Characterization of animal genetic resources is a prerequisite for successful management programmes and the formation of directives in national livestock development. Data on genetic variation within and among breeds provide information for management at breed level (i.e., control of inbreeding) and help, as well, to identify divergent breeds/populations that may harbor distinct genotypes due to local adaptation and are, therefore, worthy of conservation efforts. The Hungarian Merino breed has gone through a huge reduction in number during the last decades and a modern molecular genetic survey is now becoming more than necessary, in order to characterize the within-breed genetic diversity and structure and to establish as well the appropriate guidelines for preserving these, possibly, unique genetic resources.
The aim of the present study was to obtain information about the levels and patterns of genetic diversity and substructuring in Merino populations located in the central, north and east parts of Hungary, as well as, to examine the relationship between the genetic structure and spatial distribution of the herds, as a precursor toward the implementation of future management and conservation programmes.

Sample collection and DNA extraction
Hair samples were taken from 179 individuals grouped into six Hungarian Merino flocks (populations) ( Table 1). Four (Balmaz ujv aros, Karcag, Dorm and and Tiszacsege) out of the six populations are located in the east part of Hungary, where the other two were sampled from the north and central parts of the country (Herceghalom and Allampuszta, respectively) (Supplemental Figures 1 and 2). All the animals are registered in the Herd Book of the Hungarian Merino breed and, based on that information, sampling of close relatives was avoided. Genomic DNA was extracted from the hair follicles with simple proteolysis using proteinase K (Sigma-Aldrich, St. Louis, Missouri, United States), according to the FAO protocol [23].

Ethical statement
No harmful invasive procedure was performed and there was also no animal experimentation according to the legal definitions in Europe (Subject 5f of Article 1, Chapter I of the Directive 2010/63/UE of the European Parliament and Council).

Genotyping of microsatellite loci
Three multiplex PCR reactions amplifying simultaneously 16 microsatellite loci were carried out ( Table 2). Microsatellite markers were chosen based on previously published studies, as well as, from the panels recommended by the Food and Agriculture Organization [24] and the International Society for Animal Genetics [25]. All the PCR multiplexes were performed in a final volume of 5 ll, containing 2 ll of 1X KAPA 2 G Fast Multiplex PCR Kit (Kapa Biosystems, Wilmington, Massachusetts, USA), 2 ll of primer mix and $20 ng of template DNA.
Cycling conditions for all multiplex amplification consisted of an initial 95 C denaturation step for 3 min followed by 35 cycles of 30 s at 94 C, 90 s at 60 C and 60 s at 72 C, with a final extension at 72 C for 20 min. Fluorescently labeled PCR products were separated on an ABI3500 Genetic Analyzer (Applied Biosystems, Foster City, California, USA). Alleles were sized and individuals genotyped using the STRand v.2.4.59 software [26].  [30]. Estimates and confidence intervals based on resampling schemes were provided for all F-statistics: the Jackknife procedure was applied over loci and a confidence interval of 95% was computed with 1000 bootstraps. Moreover, the significance of deviation from Hardy-Weinberg equilibrium (HWE) for each population (based on F IS values) was assessed with randomization procedures, implemented also in Fstat v.2.9.3.2, followed by strict Bonferroni correction [31]. In order to define the degree of differentiation among the studied populations, we used Factorial Correspondence Analysis (FCA) of individual multilocus genotypes through Genetix v.4.05.2 software. Flock relationships were also explored by estimating the genetic distances according to [32] using the Gendist program in Phylip v.3.6 statistical packages [33]. A tree topology based on the genetic distances was obtained according to UPGMA (Unweighted Pair Group Method with Arithmetic Mean) algorithms using a Neighbor subroutine, implemented also in Phylip v.3.6. Assessment of support for clusters on the tree (through 1000 bootstrap resampling of loci) was performed with Seqboot subroutine. Finally, the genetic structure of the studied flocks was investigated using Structure v.2.3.4 software [34]. The software infers the number of populations into which the analyzed genotypes can be divided, by estimating the natural logarithm of the probability that a given genotype X is part of a given population K: Ln[Pr(X/K)]. In order to choose the appropriate number of inferred clusters (K), we performed ten runs fitting K from 2 to 10. All runs used a burn-in period of 100,000 iterations and a period of data collection of 100,000 iterations. Moreover, Bottleneck v.1.2.02 package [35] was used in order to validate whether the studied populations had undergone a bottleneck effect. The applied models were the SMM (one-step stepwise mutation model) and the TPM (twophase model with 95% single-step mutations and 5% multiple-step mutations). For both models, three different tests were executed: the sign test, the standardized differences test and the Wilcoxon's test.  Table 3). Summary statistics for genetic diversity within the six studied flocks are presented in Table 4. The mean number of alleles over all loci ranged from 4.2 (Herceghalom) to 7 (Dorm and). Unique alleles were observed in every population except for the Herceghalom flock. The H o values ranged from 0.669 (Herceghalom) to 0.733 (Karcag), having a mean of 0.714, whereas the respective values for expected heterozygosity (H e ) were from 0.644 (Herceghalom) to 0.734 ( Allampuszta). Average F IS values over all loci ranged from À0.010 to 0.054 and were positive in five out of the six flocks with a mean of 0.019, indicating a slight heterozygote deficiency. However, the deviation from Hardy-Weinberg equilibrium, based on 3200 randomizations and Bonferroni correction, was not significant (P < 0.05) for all tests, meaning that all F IS values were statistically non-different from zero.

Genetic structure and differentiation
Average F ST over loci and populations was generally low (0.036), showing that most of the genetic variation is explained within populations. According to pairwise multilocus F ST estimates, the highest genetic differentiation was found between the flocks of Herceghalom and Balmaz ujv aros (F ST ¼ 0.065), while the lowest was among Tiszacsege and Dorm and (F ST ¼ 0.013) ( Table 5).
According to factorial correspondence analysis results ( Figure 1) the most differentiated population seems to be the flock of Karcag, followed by that of Balmaz ujv aros. The four remaining populations formed two groups ( Allampuszta-Dorm and and Tiszacsege-Herceghalom) being close to each other. Additionally, the UPGMA tree topology verified with high bootstrap values the differentiation of Karcag and Balmaz ujv aros flocks (1000 and 815 bootstraps, respectively), while the remaining populations were grouped in two subclusters, though with lower bootstrap values (Figure 2). The differentiation patterns observed in the factorial correspondence analysis and the UPGMA dendrogram are, generally, in agreement with the pairwise F ST estimates of the studied populations.  According to Structure analysis, the most likely value of Ln [Pr(X/K)] was obtained for K ¼ 7 (i.e., seven clusters). In the graphical representation of population clustering (Figure 3) each color represents one cluster and the length of the colored segment (vertical-colored lines) shows the individual's estimated proportion of membership in that cluster. The vertical black lines separate one population from each other. Regarding the proportion of membership for each one of the six flocks in the seven clusters ( Table 6), none of the flocks was assigned in an independent cluster. The highest proportions of membership were found in Herceghalom ( ATK), Balmaz ujv aros and Karcag flocks (0.480, 0.415 and 0.441, respectively), suggesting that these populations are less heterogeneous among others. However, the resulting high membership value of the Herceghalom flock should be treated with caution, as it could also be attributed to its lower sample size (n ¼ 15). Finally, Bottleneck analysis did not reveal any statistically significant results for both SMM and TPM models in all three inferred tests, pointing out that none of the Hungarian Merino flocks have experienced a significant bottleneck effect.

Discussion
The Hungarian Merino sheep breed has more than 250 years of breeding history and has derived from the original fine-wool Spanish Merino rams (like most of the European Merino sheep breeds). During these   years many other Merino and local Hungarian breeds have contributed to the development of the breed. But as the cheap wool arrived to Europe from Australia, the fine-wool Merino breeds started to decrease in number [2]. The Hungarian Merino also went through this reduction and the breed became endangered with a number of individuals less than 5.000. Despite the very low population size, the Hungarian Merino breed seems to possess a quite noticeable level of genetic diversity. The mean number of alleles and observed heterozygosity over loci and populations were 6.1 and 0.714, respectively. These values are comparable to those reported in other studies of Merino sheep populations using microsatellite markers. According to [22], the average observed heterozygosity was 0.769 in two Merino-derived breeds [9] reported a mean H o of 0.726 across six Merino populations, whereas [16] showed a mean observed heterozygosity of 0.71 in three Italian Merino-derived sheep breeds and [36] presented an average H o of 0.705 in three Australian Merino flocks. The level of genetic variation found herein is also similar to that reported in other (non-Merino) European and non-European sheep breeds/populations. For example [10], calculated a mean observed heterozygosity of 0.69 in four Romanian local sheep breeds, while in [21] the mean H o value was 0.75 across two Bulgarian indigenous breeds. Moreover [11], reported an average H o value of 0.67 while studying seven North-African ovine breeds and [8] estimated a mean observed heterozygosity value of 0.75 in five Algerian and four Turkish local sheep breeds. However, the mean allelic number in our study was, in some cases, lower compared to the above studies, which could somehow be attributed to the lower sample (flock) sizes used. Finally, according to the estimated inbreeding coefficients (F IS values), the six studied Hungarian Merino flocks exhibited a slight heterozygote deficiency. However, none of them deviated significantly from Hardy-Weinberg equilibrium, suggesting that no inbreeding phenomena were present in our sample set, which might also indicate efficient management of the flocks by the breeders.
The genetic differentiation of the Hungarian flocks understudy was generally low, as reflected by the estimated total F ST value (0.036), the extended pattern of admixture in Structure analysis, as well as by the noticeable level of genetic clustering in UPGMA and FCA analyses, indicating that most of the genetic variation is shared within the Merino populations. In addition, the geographic distribution of the Merino herds doesn't seem to be in concordance with the low level of genetic differentiation found between them. For example, the Allampuszta flock, being located in the central part of Hungary, is genetically close to the Dorm and population that comes from the east part of the country. On the other hand, Kargag and Balmaz ujv aros flocks tend to differentiate genetically from Dorm and and Tiszacsege herds, despite the fact that all four populations are geographically close to each other, mapped in the east part of Hungary.
The low level of genetic substructuring found herein is somewhat expected as many different sheep breeds have been involved in the formation of the Hungarian Merino breed. However, it might be also due to cross-breeding practices followed often by farmers, in order to avoid inbreeding and/or increase productivity, which is common in many other  countries as well. This scenario was also evident in [20] studying the genetic structure of various herds in the Lesvos sheep breed. In addition [15], showed that the genetic differentiation of Hungarian Tsigai sheep populations was, in most cases, negligible and [22] reported a pairwise F ST value of 0.021 between two Merino-derived breeds. Nevertheless, two Hungarian Merino flocks, Karcag and Balmaz ujv aros, tended to differentiate more compared to all others. These two flocks showed the highest mean F ST values (0.048 and 0.041, respectively) and were also quite isolated in both FCA and UPGMA analyses. Despite the fact that none of the six flocks was assigned in an independent cluster in Structure analysis, Karcag and Balmaz ujv aros populations were less heterogeneous showing, relatively, high proportions of membership in one cluster (0.441 in cluster 1 and 0.415 in cluster 6, respectively). Consequently, the results of our study suggest that these two flocks are genetically unique to some extent. Regarding the Karcag flock, this belongs to an experimental station (Research Institute of Karcag, Hungarian University of Agriculture and Life Sciences) and is being intensively selected for many generations now to increase meat productivity and resistance to parasites. To our knowledge, it's considered to be a quite closed population, as only Hungarian Merino rams that meet the phenotypic standards of the breed are accepted for reproduction in this flock, thus justifying its relative genetic homogeneity.
The Hungarian Merino sheep are well-adapted to the semi-arid climate and pastures of the country and utilize in the most economic and effective way the poor grasslands, thus being worthy of conservation efforts. The involvement of various sheep breeds in the development of the Hungarian Merino breed, as well as, a possible gene flow between herds through cross-breeding practices, favor the existence of a high degree of genetic diversity at the subpopulation (within flocks) scale, despite the huge numerical decline of the breed during the last decades. A high level of genetic variability preserves the potential of adaptation to environmental changes; this becomes increasingly important due to recent climate changes and reduces the negative effects of inbreeding that are common in livestock farming [20]. On the other hand, we should also keep in priority the genetic homogeneity that some populations may have to some extent, such as Karcag and Balmaz ujv aros flocks, as these might constitute unique genetic resources.

Conclusions
In brief, we assessed the genetic diversity and structure of the Hungarian Merino sheep breed by genotyping six populations in 16 microsatellite markers. A high level of within-population genetic variation was revealed which may have been the result of the involvement of many different sheep breeds in the formation of the Hungarian Merino breed and of gene flow through cross-breeding. In addition, two flocks tended to be genetically more distant, being also less heterogeneous compared to all others. The outcome of our study could be a starting point for a planned breeding strategy of the Hungarian Merino breed, by keeping the withinflock genetic variability in priority, as well as, by preserving the potential genetic uniqueness with close monitoring of the inbreeding. Additional work, with more flocks and representatives of them, analyzed with a higher number of molecular markers, would be also very useful in this direction.