Mitochondrial DNA variability in Gyimesi Racka and Turcana sheep breeds

The current knowledge and documentation on the origins and relationship between Gyimesi Racka reared in Hungary and the Romanian Turcana is rather controversial. Lack of information and scientific reliable proofs for the divergent theories found in the two countries motivated us to implement a trial using molecular methods to assess the genetic distance and diversity in the two breeds. Hair follicles were collected from Gyimesi Racka (2 phenotypes) and from Turcana (6 ecotypes). The 599 bp segment of the D-loop region of the mitochondrial DNA was sequenced. Altogether, 42 haplotypes were identified, while 23 were found in both populations. Populations were highly diverse according to the haplotype and nucleotide diversity indices. AMOVA analysis showed that most of the variation was observed within populations (98%), indicating a weak genetic structure between the two breeds. Animals were grouped into seven groups based on their phenotype; however genetic distances among them were also low. Tajima’s D, Fu’s Fs, goodness-of-fit statistics, mismatch distribution and network analysis suggested recent demographic expansion. Current comprehensive mtDNA study indicates that there is very low level of genetic differentiation between the Gyimesi Racka and Turcana populations therefore they are de facto one trans-boundary breed.


INTRODUCTION
Knowledge on the evolutionary history of sheep breeds and information on genetic variation among and within breeds is vital for the conservation and breedspecific management efforts (Rege & Gibson, 2003;Lancioni et al,. 2013).
Zackel sheep are widely dispersed throughout Central, Eastern and Southern Europe and together with the Tsigai group are considered as the main indigenous breed groups found in these regions (Kusza et al., 2008;Savic et al., 2013).Moreover, Lawson-Handley et al. (2007) have found the Tsigai group to be strongly influenced by the Zackel and merged the two groups in their studies regarding the genetic structure of European sheep breeds.
Little is known regarding the origin of the Zackel sheep, although it is commonly accepted that this breed may have originated in the Carpathian basin, and that it was widely spread to Center (as far as Germany), Eastern (up to the Crimean peninsula) and Southern (up to continental Greece) by the nomadic Vlach sheep breeders as early as the 14th century.
Historical evidence suggests the presence of similar types of sheep as the Zackel in ancient Egypt, from where the migrations in different periods took these breeds to the Middle East and Europe (Draganescu & Grosu, 2010).Currently, Zackel sheep can be found in 14 countries from Central-, Eastern-and Southern Europe, being regarded as triple-purpose breeds (meat, milk and wool), commonly managed under extensive lowinput production systems.Phenotypically, breeds from the Zackel group are small or medium sized, with typical long coarse wool and spiraled horns.Production levels of Zackel sheep are in general modest, typical to unimproved/mountain type breeds, with milk yields varying greatly, from 40 to 150 kg of milk/lactation (Padeanu, 2010;Kawecka & Krupinski 2014).
Romanian indigenous Turcana, accounting for over 10 million breeding ewes (Ilisiu et al., 2010), is one of the most representative breed belonging to the Eastern European Zackel group.High numbers of animals, diverse geographical conditions and divergent selection have lead to a high within-breed phenotypical heterogeneity, with adult body weight of 30 to 60 kg in ewes and 50 to 120 kg in rams, growth rates in un-weaned lambs of 110 to 275 g/day and prolificacy rates ranging between 103 and 140% (Padeanu, 2010).Turcana breed has five recognized varieties/ecotypes (Draganescu & Grosu, 2010), out of which one is listed as endangered (Creata de Caransebes, last census 195 breeding ewes) and included in a conservation program.
The Gyimesi racka breed originates from the Carpathian basin (Arnyasi et al., 2013) as the Turcana, accounting for 1664 breeding ewes, being listed as an endangered in Hungary.Body weight of ewes is between 45 to 50 kg and that of rams ranges from 80 to 90 kg.The lambs' growth rates are 240 to 320 g/day, with an average prolificacy up to 110% and milk yields 40-50 kg per lactation.Given the low numbers in Hungary, and thus the higher expected inbreeding rate, the breed has considerably lower within-breed phenotypic heterogeneity compared to Turcana.Debates regarding the origin of the breed and the level of admixture with the Romanian Turcana are frequent.
We hypothesized based on both the common history that Hungary and Romania have and on the recent literature survey that Turcana and Gyimesi racka populations are de facto one trans-boundary breed.To the best of our knowledge, no other comparative study concerning the Zackel strains between European countries exists up to this moment.Furthermore, this is the first attempt to evaluate genetic variation and diversity within the Turcana population, arguably the most numerous sheep breed found in Europe.
The aim of this study was to acquire initial information about the genetic characterization of the studied breeds and to analyze, through the use of the mitochondrial DNA molecular marker, the among-and withinbreed genetic diversity and extent of admixture among the Romanian Turcana and Hungarian Gyimesi Racka sheep breeds.

MATERIAL AND METHODS
The number of animals per population, breeding site, country and phenotypical characteristics are given in Table 1.Sampling was done in private commercial farms from 50 Gyimesi Racka in Hungary and 40 Turcana in Romania, respectively.
Hair samples with follicles were collected from genetically unrelated individuals by plucking out.The genomic DNA was extracted using standard protocoll (FAO/ IAEA, 2005) and then kept at -20°C.
A 599-bp fragment of the mitochondrial control region between positions 15957 and 16555 of the sheep mitochondrial genome (access no.AF010406) was amplified by the polymerase chain reaction (PCR) with primers according to protocol elaborated by Niemi et al. (2013), as follows: forward: 5'-GTTTCACTGAAGCAT-GTAGGG-3' and reverse: 5'-GTATTGAGGGCGG-GATAAAT-3'.PCR was performed on a PTC-200 thermocycler (MJ Research, Inc.) in a total volume of 10 mL of the following mixture: 50-70 ng DNA, 0.1 μM each primer, 0.2 mM dNTP (Promega, USA), 2.5 mM MgCl 2 (Promega, USA), 1.5 U of GoTaq DNA Polymerase (Promega, USA) and 5X Colorless GoTaq Reaction Buffer (Promega, USA).The cycling conditions were as follow: initial denaturation step for 5 min at 94°C, followed by 35 cycles at 94°C for 1 min, 60°C annealing for 1 min and 72°C extension step for 1 min.A final step was also performed at 72°C for 10 min.Sequencing reactions were performed by Macrogen Inc. (The Netherlands).
Ten samples from the Gyimesi Racka and eight from Turcana were not successfully sequenced therefore they were omitted from further analysis.The remaining 72 sequences were aligned using CLUSTALW (Larkin et al., 2007) and trimmed to the shortest one in order to get the same length of segment (526 bp, between positions 16016 and 16541).
DnaSP 5.10 software (Librado & Rozas, 2009) was used to calculate haplotype diversity (Hd), nucleotide diversity (π), number of segregating sites (S) and average number of nucleotide differences (K) (Table 1).ARLE-QUIN 3.1.(Excoffier et al., 2005) was used for calculating Tajima's D test, Fu'Fs test, goodness-of-fit statistics, and reflect the shape of mismatch distribution.Nonsignificant Harpending's raggedness index and sum of square deviations (SSD) indicates a good fit and support of expansion.Smooth, unimodal distributions can be regarded as being representative of a population expansion, while the ragged and bi-or multimodal mismatch distributions indicate a stationary population (Rogers & Harpending, 1992;Excoffier & Schneider, 1999).Population genetic structure was also calculated among and within populations by an analysis of molecular variance (AMOVA) using ARLEQUIN 3.1.
To get the appropriate DNA substitution model for analysing our sequences FindModel test and MEGA version 6 program (Tamura et al., 2013) were used and both supported Kimura two-parameter model (K2, Kimura, 1980).Gyimesi Racka and Turcana individuals were grouped into seven groups based on their phenotypes (Speckled-faced, Black-faced, Breaza of Petrosani, Freckled-faced, Transhumant of Sibiu, White-faced and Curly of Caransebes).These groups are parallel with Romanian Turcana ecotypes.Genetic diversity values were calculated between Gyimesi Racka and Turcana and among the above mentioned groups with MEGA v.6 program (Tamura et al., 2013).Phylogenetic analysis of the studied sheep haplotypes and one outgroup sequence from Capra hircus (GenBank: KJ940969) was conducted using Neighbour joining (NJ) algorithm with 1000 bootstrap replicates by MEGA v.6.program (Fig. 1).In addition, NETWORK 4.111 (Bandelt et al., 1999) was used to construct median-joining networks among our haplotypes (Fig. 3) and between our and previously published 1005 wild and domestic sheep sequences (Fig. 4).Networks fit more to depict intraspecific phylogenies than tree algorithms because they allow the co-existence of ancestral and descendant alleles in a sample (Posada & Crandall, 2001).
The research activities were performed in accordance with the European Union's Directive for animal experimentation (Directive 2010/63/EU).
The results obtained by the AMOVA analysis showed that most of the variation was observed within populations (98%), indicating a weak genetic structure.Pairwise F ST value of 0.019 for Gyimesi Racka and Turcana was not significantly different from zero at the level P < 0.05.
Current results are in agreement with those published by Kawecka & Piorkowska (2011) on genetic distance and structure of two Polish native Zackel breeds, with high genetic similarity being found for the Podhale Zackel and Coloured Mountain, respectively.
Genetic distance was calculated using the Kimura 2-parameter model and low value (0.004) was found between Gyimesi Racka and Turcana populations.Current findings are in accordance with those reported by Tapio and coworkers (2010), which following a genetic cluster analysis for a number of 52 sheep breeds found a sub-cluster anchored by the Zackel populations, and concluded that this reflects the common ancestry for the majority of breeds within the sub-cluster.Furthermore, results by Tapio and coworkers (2010) confirmed the assumption of Lawson-Handley and coworkers (2007) that the Tsigai group was strongly influenced by Zackel.
Neutrality tests were implemented to get information about recent historical demographic events.Significant negative values indicate a sudden expansion in population size, while significant positive values indicate a population subdivision or recent population bottlenecks.In our populations, neutrality tests values Tajima's D (-1.734 and -1.814) and Fu's Fs (-18.549 and -2.721) were both negative but Tajima's D values were significant for both populations indicating expansion, while Fu's Fs was significant only for Gyimesi Racka.
The mismatch distribution for the Gyimesi Racka was unimodal, bell-shaped and smooth suggesting that it had undergone a demographic expansion.It was further supported by a low and statistically non-significant Harpending's raggedness index (0.056), and a low SSD value (0.008).Multimodal mismatch distribution and the low and non-significant raggedness and SSD index (r: 0.032 P=0.65; SSD: 0.022 P=0.15) obtained for the Turcana are not fully consistent with demographic expansion.
Regrettably, information regarding the census of the two populations in Romania and Hungary could not be traced-back or followed throughout the last decades (such as in the case of purebred horses from stud-farms) in order to test our results regarding the expansion of the two populations, due to the lack of data and herd-books for .the breeds.Draganescu (2005) estimates that in 1970, 5.8 million Turcana were reared in Romania, representing roughly 55% of the today's census for the breed.
Neighbor-Joining tree was constructed and used to measure differences within and between the observed distinct haplotype groups (Fig. 1).The low internal resolution, weak bootstrap values and haplotypes without grouping correlation to their sampling sites in the NJ tree of the mitochondrial CR haplotypes also suggested that there was no or very little geographical or breeding structure of the two populations.
Populations were classified into seven ecotypes: Transhumant of Sibiu, Curly of Caransebes, Breaza of Petrosani, Black-faced, White-faced, Speckled-faced and Freckled-faced, according to their phenotypes (Fig. 2).The first three groups are found only in Turcana while the Freckled-faced variant is found only in Gyimesi Racka breed.
Network tree did not show any clear pattern for distribution (Fig. 3).All ecotypes were closely related (Table 3).Highest genetic distance was observed for the Curly of Caransebes Turcana (0.009-0.012), which sets the phenotype apart from the others.This ecotype of Turcana is being reared in the South-Western Romania, in the Banat's highlands region, being well adapted to its environment.The Curly of Caransebes its found in the same area as the Racka sheep.As a result, to some extent, it is possible and plausible that gene flow occurred between the two populations.This might explain the "curly" aspect of the fleece staple, found in both populations.As mentioned, the Curly of Caransebes is the only one of the Turcana ecotypes to be under genetic preservation in Romania.Our findings suggest that this ecotype can be regarded as an emerging-breed, and efforts to better preserve this population should be taken.
Opposite to our results, subdivision within the Turcana breed was expected due to large population size (over 10 million breeding ewes), isolation in terms of geography and different breed selection traits and management.
Conversely to current findings, Draganescu (2005)  Distances were calculated using the Kimura two-parameter model (Kimura, 1980) with 1000 replicates in the bootstrap test.Sequence from Capra hircus (access.no: KJ940969) is used as outgroup.

Gyimesi Racka and Turcana considered as one breed
A total of 42 haplotypes (GenBank accession numbers: KR054136-KR054177) were detected in the 72 animals.Haplotype diversity, nucleotide diversity and average number of nucleotide differences were high, 0.953, 0.004, and 1.424, respectively (Table 1).
The mismatch distribution of complete dataset showed a multimodel shaped distribution indicating stationary population.Fu's Fs neutral test was non-significant and negative (Fs=-17.778,P > 0.05), while Tajima D was significant and negative (D=-1.912,P ≤ 0.05).
Due to the shorter sequences from GenBank, the geographic distribution and origin of our haplotypes and that of contemporary other 1005 wild and domestic sheep, our study was based on a 240 bp fragment of the mtDNA control region (Table S1 at www.actabp.pl).
The haplotype group distribution did not indicate either a clear phylogeographic patterns or geographical structure (Fig. 4).Our haplotypes were mixed with all other haplotypes from different geographical and genetical origin of different sheep breeds.These are in accordance with the historical records, origin theory and geographical localizations of the studied populations (Juler, 2014;Kawecka & Krupinski, 2014) and also confirm Peter and coworkers (2007), who consider the South-East European Zackel breed as a 'genetic hot-spot'.Aspect on genetic diversity documented in this study indicating high levels of gene flow that have occurred between the two populations.
Current results suggest that the Gyimesi Racka and Turcana populations from Hungary and Romania, respectively, are in fact a single trans-boundary breed, thus our hypothesis has been confirmed.Conversely, reports of Cinkulov and coworkers (2008) on genetic diversity and differentiation between two Serbian Tsigai strains revealed substantial differentiation between the two populations, up to the point that authors suggested the two populations to be considered as distinct breeds.Authors' results are contradictory to other reports and to current findings, considering that the Tsigai represents a minority breed in Serbia, with an estimated census of 86 000 breeding animals (Petrovic et al., 2013).Furthermore, results of Cinkulov and coworkers (2008) on the population structure of the Pramenka strains (Svrljig, Dubska and Recka) showed elevated within-population diversity and lack of signatures of admixture between the three strains.
In the Turcana populations, the lack of phylogeographic patterns or geographical structure could be attributed  to two of the ecotypes, namely Breaza of Petrosani and Transhumant of Sibiu, wich are the most productive of the variants and the most widespread, respectivelly.As already mentioned, production levels within the Turcana populations vary greatly (e.g.milk yields from 60 kg up to 150 kg) as a result, rams from the Breaza of Petrosani ecotype are widely used for the upgrading the other Turcana populations.While Transhumat of Sibiu as the most widespread of the ecotypes (40% of the Turcana breed), it currently being used for vertical transhuman of 200 up to 300 km/year in the Carpathian mountains, thus, contributing to the level of ad-mixture among the Turcana ecotypes.A similar situation was described by Ligda and coworkers (2009) in their studies concerning genetic analysis of the Greek sheep breeds, where Karagouniko breed was found to strongly influence the other Zackel breeds due to its widely use for upgrading production levels.
Current results could prove important for the conservation and breed management efforts of the Gyimesi Racka, with a census for the breed in Hungary of less than 1700 heads and estimates for the number of breeding rams of 55 to 60.Thus, high inbreeding rates can be avoided and breed upgrading achieved by importing Turcana rams.
Further large-scale studies should be implemented for the Zackel strains located in Central, Eastern and Southern Europe, in order to assess their census, inbreeding rates and genetic distances among these breeds to help with the conservation efforts and avoid genetic diversity erosion or loss.A very important fact to have in mind is that in most of the 14 countries which rear Zackel sheep in Europe, these breeds are either listed as endangered or represent a minority breed.As such, a further largescale study to investigate the molecular genetic diversity of the Zackel sheep group in a regional or continental context would increase our knowledge of the development of gene pools of the European sheep breeds and sheep biodiversity.
Figure 1.Neighbour Joining tree of mtDNA-CR haplotype sequences.Distances were calculated using the Kimura two-parameter model(Kimura, 1980) with 1000 replicates in the bootstrap test.Sequence from Capra hircus (access.no: KJ940969) is used as outgroup.

Figure 3 .
Figure 3. Median-joining network constructed for seven ecotypes: Transhumant of Sibiu: blue; Curly of Caransebes: yellow; Breaza of Petrosani: pink; Black faced: black; White faced: white; Speckled faced: grey and Freckled faced: green.Numbers on the lines indicate the number of mutations (no number indicates a single mutation).The size of the circle is proportional to the number of animals represented.

Figure 4 .
Figure 4. Median-joining network of 42 Gyimesi Racka, 30 Turcana and 1005 sheep sequences from GenBank.The size of the circle is proportional to the number of animals represented.Numbers on the lines indicates the number of mutations while no number indicates single mutation.Grid represents Gyimesi Racka, horizontal lines show Turcana from this study.Colours represent regions of sequence origin.Italian Peninsula: light green; Iberian Peninsula: dark green; Northern Europe: grey; Baltic region: white; Central-Europe: light blue; East Europe: purple; Balcan Peninsula: dark blue; British Isles: yellow; Russia: light pink; Caucasian: cyclamen; Asia: fuchsia.

Table 1 . Basic parameters of mtDNA CR variability of Gyimesi Racka, Turcana and total dataset
*P < 0.05