Gene conservation of six Hungarian local chicken breeds maintained in small populations over time

Investigation into the genetic diversity of certain endangered native breeds of domestic animals has been in common practice for several decades. The primary objective of these investigations has been to reveal the exceptional genetic value of such breeds, both for their conservation and also to gain insight into their current genetic status, as they have been undergoing a progressive decrease in population size and general diversity; this has been compounded by the general lack of an optimal breeding scheme. In this study, we have investigated changes in the genetic diversity of six Hungarian local chicken breeds based on 29 microsatellite loci over a period of 15 years. In terms of the basic diversity measures, populations sampled in 2017 generally exhibited a lower heterozygosity and mean number of alleles and thus, experienced a higher degree of inbreeding. Although the effective population size increased, the estimates of populations sampled over different periods indicated comparatively low values, suggesting overall lower genetic variance. Pairwise FST estimates were higher in the populations sampled in 2017, showing a larger genetic distance between them. Considerable differences exist between the populations of the same breeds, which can most likely be attributed to genetic drift. STRUCTURE results have shown a clear separation between the Hungarian populations, which is in agreement with the principal coordinate analysis. The most likely clustering was found at K = 6, classifying the populations of the same breed as one group. No considerable allele loss was found in the Hungarian indigenous chicken breeds after 15 years of conservation. In general terms, after 15 years, the level of inbreeding within the populations was, in fact, higher, although this could be effectively reduced through the use of an improved mating system. Consequently, the breed management applied in the case of Hungarian local chicken breeds was found to be effective at adequately conserving their genetic variability.


Introduction
Over the last several decades, exploration has been under way regarding the genetic diversity of endangered domestic animal breeds. The loss of genetic variability, together with the limited number of commercial hybrids, both suggest that genetic diversity in poultry is one of the most endangered genetic resources [1][2][3]. There are numerous grounds supporting gene conservation in indigenous breeds of chicken. From the point of view of conservation, it is necessary for native chickens to be preserved as a genetic resource, revealing their particular and distinctive genetic value for the purpose of future breeding. They are capable of providing a source for unique alleles and facilitate enrichment with genes associated with health and quality traits [4,5]. Another justification for studies in this area is to carefully observe the genetic status of the populations themselves, as native chicken breeds in conservation breeding programmes are often maintained in small populations, and therefore more affected by loss of genetic variance and inbreeding depression. Local chicken flocks typically have no pedigree data, suffer from changing population sizes (bottlenecks), and lack properly planned breeding programmes [6]. It is clear that those unique combinations of genes which make up a specific genotype are at the risk of disappearance, in spite of the fact that they may constitute potentially beneficial traits [7].
Maintaining an effective population size (Ne) plays an important role in conservation genetics, as it is a determining factor in the development of the population in regard to inbreeding and genetic drift. Minimising loss of genetic variability is one of the main concerns of conservation breeding programmes and can be achieved through the maximisation of the effective population size, and hence, decrease the occurrence of genetic drift. Estimation of the effective size of a population and the monitoring of its changes over time are important to accurately assess the effectiveness of the breed management, particularly in small populations, where the allele frequencies change more easily over time. In the absence of selection, any change in allele frequency could solely be the result of drift and thus indicates the effective size of the population [8].
Although high-throughput SNP techniques have been more preferred for a while [9], microsatellite markers have similarly proven to be efficient molecular tools particularly in temporal studies [10] for assessing changes in genetic diversity as well as for setting priorities for conservation [11]; this has been especially applicable in the case of local chicken breeds [12][13][14][15][16][17][18]. Nevertheless, despite the high level of thoroughness and attention to detail in the investigations, any results obtained can easily become obsolete in the case of small populations where lifespans are more limited, and thus generations tend to change more frequently. Genetic diversity within these small populations is slowly, but steadily, eroded by the limited number of offspring (retained from the parents), especially in the absence of a specialised breeding scheme [19]. This is mostly the case, since it can be generally observed that the costs of animal breeding are constantly increasing and its income potential is decreasing, owing to the aggressive development and competition from intensive production practices [20].
Having a long history in this area, one of the primary activities of the National Centre for Biodiversity and Gene Conservation (Godollo, Hungary) has been the preservation of local Hungarian poultry species. The aim of this study was to investigate the extent of divergence which has developed in the genetic diversity of six Hungarian indigenous chicken breeds over a period of 15 years. We trust that our results will provide useful information regarding small populations which have been managed for the purpose of maintaining their genetic variation in gene banks for gene conservation as well as future breeding programmes.

Materials and methods
Animals were kept and maintained according to general animal welfare prescriptions of the Hungarian Animal Protection Law (1998; XXVIII). Permission to undertake experimental animal research at the National Centre for Biodiversity and Gene Conservation was granted by the National Food Chain Safety Office, Animal Health and Animal Welfare Directorate, Budapest (Permission no. PE/EA/2485-6/2016). All experimental methods described herein were approved by the Institutional Ethics Review Board.

Sample collection and DNA extraction
Six native Hungarian chicken breeds were included in this study: the Yellow (HY), the Speckled (HS) and the White Hungarian (HW), as well as the Speckled (TNNS), the White (TNNW) and the Black Transylvanian Naked Neck (TNNB). A major breeding program has resulted in good dual-purpose, medium-sized Hungarian chicken breeds which have been propagated not only for their sizable egg production (140-150 pieces per year per hen, egg weight: 50-55 g) but also for their excellent meat quality (fine-fibred, highly-palatable). In the case of the HY, HS and HW breeds, the body weight of the hens is 2.0-2.3 kg, while that of the roosters is typically 2.5-3.0 kg. The Transylvanian Naked Neck chickens are characterized by their featherless necks and only slight plumage on the top of their heads. These breeds are extraordinarily resistant and known for their good winter-laying ability; however, their body weight (hens: 1.7-2.2 kg, roosters: 2.2-2.4 kg) and egg production (80-100 pieces per year per hen, 48-50 g) are slightly lower when compared to the feathered Hungarian breeds [21]. 1 ml of individual blood samples were collected from each breed following the 15 year period, from thirty birds per population (20 females and 10 males), and sodium citrate solution was used as anticoagulant. Overall, 180 DNA samples were extracted using the traditional salting-out method [22], modified for chickens [23]. DNA concentration of the samples was measured using a Nano-Drop spectrophotometer (Thermo Fisher Scientific, USA), then equalized to 5 ng/μl. The samples were stored at -20˚C until use.

Microsatellite genotyping
Individuals were genotyped using the same set of 29 microsatellite loci recommended by the FAO [24], as in our previous study [23], to obtain information as to how the genetic diversity of Hungarian indigenous chicken breeds has changed over a period of 15 years. To reduce the costs of the study, tailed primers were used for genotyping; therefore, instead of applying directly-labelled primers, a universal sequence (5'-CAGGACCAGGCTACCGTG-3') was placed before the 5' end of the forward primers [25]. This universal primer (tail) was labelled with varying fluorescent WELL-RED dyes (D2, D3, D4); and thus-with proper optimizationthe markers could be grouped into 10 multiplex PCR reactions, taking into consideration the fragment sizes, and the PCR products were pooled according to the fluorescent colours for the electrophoretic analysis (S1 Table).
The final volume of 15 μL master mix contained 5 μM primer, 10x Dream Taq Buffer with 20 mM MgCl 2 (Thermo Scientific), 25 mM dNTP mix (Thermo Scientific), 5 U/μL Taq DNA polymerase (Dream Taq DNA polymerase, Thermo Scientific) and 25 ng genomic DNA. The ratio of the tailed forward and reverse primers and its fluorescently-labelled tail sequence was 1:2:1. PCR profile: 95˚C for 4 min denaturation followed by 30 cycles of amplification: 95˚C for 15 sec, 60˚C for 30 sec and 72˚C for 1 min, with a final extension at 72˚C for 9 min (Kyratec Trinity Supercycler). The annealing temperature and the number of cycles were optimized for each multiplex reaction respectively (S1 Table).
PCR products were detected by capillary gel electrophoresis using an automatic DNA sequencer (GenomeLab ™ GeXP Genetic Analysis System, USA) according to the manufacturer's instructions. A 400 bp allele ladder (400 bp size standard) was used to determine the allele sizes. For each population, 5 reference samples were selected from our earlier study [23] and genotyped with the tailed primer sets. The results were evaluated per sample using the GenomeLab Genetic Analysis System software (Beckman Coulter). Data correction was performed on the basis of the previously-obtained allele frequencies [23] and the allele sizes of the reference samples.

Genetic diversity estimates within and between populations
Allele frequencies, mean number of alleles (MNA), and expected and observed heterozygosity (H E and H O ) were estimated across all loci in each population by the use of an Excel Microsatellite Toolkit [26]. Effective population sizes (Ne) were calculated using the linkage disequilibrium (LD) method (minimum allele frequency: 0.01) [27] for each population as well as the Jordan/Ryman temporal method [28] for generations as they are implemented in NeEstimator v2.1 [29]. Estimating the effective population size based on gametic disequilibrium, has the advantage that it only requires a single sample from the population. The LD method has greater significance when the effective size of the population is small and may therefore be of considerable use to conservation biologists, who are often faced with low or potentially low Ne [30]. In general, small population sizes are more expected to give rise to relatively-high levels of LD, and conversely, large population sizes to low LD levels. A noteworthy feature of this means of estimation is that a measurement taking at a single, given point in time is suitable for providing a reliable assessment of effective size. In the field of conservation, it is usually the most-recent population sizes that are of greatest interest, and for which unlinked loci are the most relevant [31]. The temporal method utilizes data on the rate of change in allele frequencies between samples taken at different points in time and calculates an average estimate of genetic drift over the loci [28]. The inbreeding coefficients (F IS ) of each Hungarian population, including the test for significance and Weir & Cockerham [32] estimates of Wright's fixation indices (F IT , F ST , F IS ), were calculated with FSTAT software [33]. Values of standard error were generated by jack-knifing over both loci and the populations. A molecular analysis of variance (AMOVA) was calculated using GenAlEx software [34,35], which partitions genetic variability into separate components. A principal component analysis (PCA) by the varimax rotation method was performed on the basis of diversity estimates (MNA, H E , H O , F IS ) using SPSS 11.5 software [36]. Pairwise F ST values were computed with GENEPOP software [37].

Cluster analysis
Clustering of individuals from multi-locus genotypes was performed using a Bayesian-based approach with the STRUCTURE software [38]. The analysis was performed using an admixture model with correlated allele frequencies and run with 20,000 burn-in periods, followed by 50,000 iterations for each K number ranging from 2 to 12. For each K value, 100 independent runs were performed. Pairwise comparisons of the 100 solutions for each K value were made in a greedy algorithm, and clusters with the highest average similarity index were visualized using the CLUMPAK software package, which provided a summary and graphical representation of the STRUCTURE results [39]. The most likely clustering (ΔK) was estimated by STRUCTURE HARVESTER web v0.6.94 [40] using Evanno's method [41]. In contrast, a principal coordinate analysis (PCoA) via a covariance genetic distance matrix with data standardization was conducted with GenAlEx software 6.5 [34,35] and a two-step cluster analysis was performed on the individual diversity scores (originating from the PCA) using the SPSS 11.5 software [36]. PCA creates a clustering of the individuals based on their diversity scores and the direction of the changes observed can be easily visualized as the two axes of the plot represent different meanings.

Genetic diversity of Hungarian chicken populations after 15 years of conservation breeding
The diversity of six Hungarian chicken breeds was analysed in relation to a new set of the same breeds sampled following a period of 15 years.
A principal component analysis (PCA) of the summary diversity statistics was performed to collapse several parameters into two principal components (Fig 1).     Table).  (Table 3).
The results of the STRUCTURE analysis are given in Fig 2. Due to the fact that this software places individuals into clusters without any previous information, all the 12 populations were analysed together for the purpose of acquiring reliable information on how the populations differentiated after 15 years of conservation. At the lowest K value (K = 2), populations of the HW and TNNW breeds diverged from the others and clustered together until K = 5. At K = 3, the HY 2002 and HY 2017 stocks formed a cluster and remained together until K = 12. The Transylvanian Naked Neck breeds (except for the TNNW) grouped with the HS stocks, which are slightly overlapping with the HY breed. At K = 4, the two TNNB populations split to form a separate cluster. The HS stocks are still overlapping with the HY and TNNB breeds, but at K = 5 they showed a clear population structure. The clustering of highest probability was obtained at K = 6 (N = 100) based on Evanno's method, where the HW and the TNNW breeds separated from each other and this state remained until K = 12. From K = 6 to K = 9 the clustering was the same, however, the HW 2002 showed a slightly-mixed population structure from K = 7. At K = 10 the two populations of the TNNS breed formed separate clusters, while at K = 11 the TNNB 2002 and TNNB 2017 stocks split off as individual groups. K = 11 and K = 12 showed identical solutions, indicating that a maximum of eleven clusters illustrate the data.

Discussion
In the current study, we investigated the changes in the genetic diversity of six Hungarian local chicken populations based on genotyping individuals at 29 microsatellite loci over a period of 15 years. The observed heterozygosity values were decidedly similar to the expected ones, and F IS estimates were not significantly different from zero, suggesting that the populations sampled in 2002 were close to the Hardy-Weinberg equilibrium state. After 15 years of conservation, heterozygosity estimates differed significantly in most of the breeds, consequently, the inbreeding coefficient (F IS ) increased in all breeds with the exception of the HS. This is in agreement with the principal component analysis based on diversity estimates (Fig 3). Our basic diversity measurements, obtained after a period of 15 years, maintained their similarity to those estimated for 65 chicken populations at the same 29 microsatellite loci [42], and comparable to those reported by Zanetti et al. [15] for chickens native to Italy; however, authors of the latter study found a decrease of inbreeding coefficient for three out of five Italian populations after only 4 years of conservation [10]. Although the effective population size increased over the 15 year period, the Ne estimates of populations sampled at different points in time indicated quite low values, suggesting less genetic variance in those populations which were sampled 15 years later, due to the higher level of inbreeding ( Table 1). The two methods used in our study for the purpose of estimating effective population size make calculations on different platforms: on the basis of linkage disequilibrium (LD method) and based on genetic drift (temporal method). The LD value is known to be less affected by drift compared to other factors. Accordingly, the value obtained by the LD method was higher, due to the fact that the effect of the factors influencing the LD value within a stock (e.g. population fragmentation, exchange of individuals within population, inbreeding) was smaller, owing to the appropriate breeding system. However, genetic drift did also occur over the 15 year period, as the population sizes were limited due to financial constraints. Thus, all things considered, the effective population size is lower in the case of the temporal method. These results are consistent with the findings in 37 European local chicken populations based on Wright's Ne calculation [19].
In that study, the HW, TNNW and TNNB breeds were found to be the most inbred among the Hungarian populations, which was shown in our principal component analysis (Fig 1) as well.
The overall population inbreeding coefficient (F IT ) indicated that heterozygote deficiency was lower in the Hungarian populations sampled in 2002 and it was almost fully explained by between-population differentiation (F ST , Table 2). In contrast, the populations, after 15 years of conservation, showed a higher F IT estimate, suggesting the increased ratio of homozygote animals, this figure is still not as high as in the Swedish, Polish and Italian local chicken breeds [15,43,44]. The high F ST values showed a clear population structure of both Hungarian sets, which is in agreement with the STRUCTURE results and notably similar to the findings in French native chicken breeds [45]. Compared to those Hungarian chicken populations which were sampled in 2002, those sampled in 2017 exhibited a higher degree of inbreeding (F IS ); however, they were well managed in terms of maintaining the genetic variation found 15 years prior. The F IS value of the Hungarian populations sampled in 2017 is even coherent with that of the Polish local breeds investigated [44]. After 15 years of conservation, the inbreeding level increased from 1% to 6% (F IS , Table 2), and the between-population differentiation changed from 23% to 25% (F ST , Table 2), which is consistent with the molecular analysis of variance (AMOVA, S2 Table). The reason for this might be attributed to the small population size itself (approx. 200 birds), as, in the cases of these populations the decreasing effects of genetic drift in the diversity may occur even more readily than allele losses.
The pairwise F ST estimates observed between populations sampled in 2002 were generally a bit smaller than that 15 years later, indicating lower genetic similarity between the stocks sampled in 2017 (Table 3). The population differentiation was still not as high as local chicken breeds in Britain [16]. The difference had almost doubled between the HW and TNNW breeds by 2017 (from 0.15 to 0.27); however, in other cases, the similarity actually increased, mostly between the TNNS and HY breeds (from 0.27 to 0.21). The highest observed similarity between the two populations of the same breed was 0.08 (TNNS 2002 -TNNS 2017 and TNNB 2002 -TNNB 2017 ), which is remarkably close to the pairwise F ST estimate (0.15) of two distinct breeds (HW 2002 -TNNW 2002 ) sampled in the same year. This suggests that the diversity of TNNB and TNNS populations has changed more than that of the others. This concurs with the STRUCTURE results.
The STRUCTURE analysis (Fig 2) showed a clear clustering of the Hungarian local chickens. A similar population structure was observed for Italian native chicken populations even following a mere 4 years of conservation [10]. The results suggest that the two populations of the HY breed are indistinguishable, which is in concordance with the findings of the PCoA analysis based on genetic distances between individuals (Fig 3). Furthermore, the HW and TNNW breeds clustered together, similar to our previous study [23], but from K = 6 (most likely clustering), they were clearly separated, and their populations sampled in 2002 and 2017 remained together until K = 12, respectively, which was confirmed by the PCoA analysis ( Fig  3). Individuals of the TNNB populations clustered with the HS and TNNS breeds at low K value, but constituted their own cluster at K = 5. Similar phylogenetic relationships were displayed by the PCoA analysis (Fig 3). The two stocks of this breed remained together up to a moderate level of resolution in the STRUCTURE analysis, but appeared as clearly distinguishable populations at K = 11. The two populations of the TNNS breed already appeared as separate groups at K = 10, suggesting a higher genetic differentiation between the stocks sampled at different points in time. This was confirmed by their relatively high pairwise F ST estimate ( Table 2).
The findings of the principal coordinate analysis are displayed in Fig 3 and show a strong concurrence with the STRUCTURE results. Individuals of the HY 2002 and HY 2017 populations formed a group distinct from the other breeds. The populations of the HS and TNNS breeds clustered with the TNNB flocks. Some birds of the HW 2002 population bore similarity to the TNNS and TNNB breeds. This is also indicated by the STRUCTURE results at K levels higher than the most likely clustering; a slightly complicated population structure was discovered for the HW 2002 , which is comparable with the results obtained for an Italian local chicken breed [15]. Nevertheless, the majority of individuals from the HW 2002 population were isolated from the other breeds on the left side of the graph; populations of the HW and TNNW breeds were almost completely differentiated.
In summary, no considerable allele loss was found in the Hungarian local chicken breeds after 15 years of conservation. However, one allele of the TNNB breed, having a frequency of 30%, disappeared; there was no significant change for the rest of the microsatellite loci. This may indicate a natural or artificial selection effect. It appears that there is a notable distance between the populations of the same breeds sampled at different points in time, especially between those of the TNNS 2002 -TNNS 2017 and TNNB 2002 -TNNB 2017 stocks, respectively, most likely due to genetic drift. This was also found in our previous study, which investigated subpopulations of the same breed maintained separately for 30-40 years as closed populations [23]. As for the HS and TNNS breeds, there were not enough eggs available in 2015 for the reproduction in our gene bank; therefore 300 eggs (nearly half of the hatched eggs) were taken from an elite breeder who had been maintaining these breeds for almost 20 years. This replacement had a remedial effect, especially for the HS breed, as its inbreeding coefficient not only did not deteriorate further, but rather displayed an improvement. The aforementioned facts may also have resulted in a slightly greater change in their genetic makeup.
The change of alleles between generations may be attributed to inbreeding, genetic drift, or even the sampling itself, due to the fact that these breeds are maintained in the gene bank having a relatively small population size (approx. 200 birds). Generally, the level of inbreeding appears to be higher after 15 years of conservation, but it is capable of being reduced through the use of a sire-rotation scheme different from that being used at present. The allele frequency of small populations may be subject to fluctuations by leaps, therefore, the allele set itself has proven to be a more stable indicator, as it changes only when the gene pool of the populations changes significantly, and this is not the case here.
The breed management which has been applied to these populations, i.e. maximising the genetic variance based on phenotypic traits and sire rotation used between families, has sufficiently demonstrated its suitability for adequately conserving genetic variability with the addition of some improvements. Thus, we hereby suggest the implementation of a controlled random male rotation scheme in place of the currently constant pattern of rooster exchange between families. Introducing new individuals into the breeding system from time to time may also help reduce the inbreeding level. Increasing the effective population size with either the use of a higher number of roosters in smaller families, or through the application of rooster exchange during periods of egg collection for pedigree, is recommended to conserve the genetic variance more effectively, however this practice does have its financial constraints.
Supporting information S1