Genetic Diversity and Identification of Homozygosity-Rich Genomic Regions in Seven Italian Heritage Turkey (Meleagris gallopavo) Breeds

Italian autochthonous turkey breeds are an important reservoir of genetic biodiversity that should be maintained with an in vivo approach. The aim of this study, part of the TuBAvI national project on biodiversity, was to use run of homozygosity (ROH), together with others statistical approaches (e.g., Wright’s F-statistics, principal component analysis, ADMIXTURE analysis), to investigate the genomic diversity in several heritage turkey breeds. We performed a genome-wide characterization of ROH-rich regions in seven autochthonous turkey breeds, i.e., Brianzolo (Brzl), Bronzato Comune Italiano (BrCI), Bronzato dei Colli Euganei (CoEu), Parma e Piacenza (PrPc), Nero d’Italia (NeIt), Ermellinato di Rovigo (ErRo) and Romagnolo (Roma). ROHs were detected based on a 650K SNP genotyping. ROH_islands were identified as homozygous ROH regions shared by at least 75% of birds (within breed). Annotation of genes was performed with DAVID. The admixture analyses revealed that six breeds are unique populations while the Roma breed consists in an admixture of founder populations. Effective population size estimated on genomic data shows a numeric contraction. ROH_islands harbour genes that may be interesting for target selection in commercial populations also. Among them the PTGS2 and PLA2G4A genes on chr10 were related to reproduction efficiency. This is the first study mapping genetic variation in autochthonous turkey populations. Breeds were genetically different among them, with the Roma breed proving to be a mixture of the other breeds. The ROH_islands identified harboured genes peculiar to the selection that occurred in heritage breeds. Finally, this study releases previously undisclosed information on existing genetic variation in the turkey species.


Introduction
According to historical evidence, turkey domestication originated in Mexico and Central America. After the Spanish conquerors brought the turkeys to Europe, the novelty and the appreciated meat characteristics permitted a rapid spread across Europe, starting from the 16th century. Since their diffusion in Europe, turkey populations were then bred divergently in the next centuries [1]. In the last 40 years, the turkey species underwent a targeted selection for meat production in order to obtain a fast-growing, heavy bird; nowadays several hybrids, born out of an intensive selection program, are used in intensive farming to produce turkey meat.
The Italian autochthonous turkey populations show a wide phenotypic variation among breeds, spanning plumage colour, body size and weight [2]. Heritage breeds can be considered an in situ genetic variability reservoir that should be maintained and their value exploited. The heritage populations, in fact, own unique characteristics making them especially valuable for their capability to adapt to harsh environments and to resist

Sampling and Genotyping
In this study, 181 birds' genotypes from the Axiom ® TurkeyHD Genotyping Array of seven Italian turkey breeds (Brianzolo: Brzl n. 31 [8,9]. The two genotype datasets were appended and underwent a joint quality control. Filtering was applied to retain for the analyses SNPs with genome position on autosome (first 30 chrs) and having a call rate >99%.
A total of 346,155 SNP markers passed the quality control, and their location on the genome was in concordance with the Turkey_5.0 genome assembly-GCA_000146605.1.

Population Genetic Diversity Analyses
The genetic diversity within and among breeds was determined using the following approaches:

•
The number of monomorphic SNPs, the expected and observed heterozygosity (Exp Het and Obs Het), the expected and observed number of homozygous SNPs (Exp Hom and Exp Hom) and minor allele frequency for each breed were calculated. • A principal components analysis (PCA) was performed based on allele genotypes using the SNP & Variation Suite (SVS) v8.9 (Golden Helix Inc., Bozeman, MT, USA). The graphical visualization of PCA was obtained by the ggplot2 R package [18].

•
The pairwise fixation index (i.e., Wright's F-statistic F ST ) was estimated using the dedicated module implemented in SVS. The F ST was estimated for all possible pairs of breed combinations. • Identity-by-state estimates of genetic distances were calculated for all pairs of individuals using PLINK 1.9. A neighbour-joining tree (NJ) was constructed-using the distance matrix as input file-by Phylip software implemented in the online WebPHYLIP tools [19]. A NJ tree was then drawn with MEGA X software [20].

•
The ADMIXTURE analysis: The most probable number of ancestral populations was identified in conjunction with the lowest cross-validation error (CV), setting the analysis with optimal number of clusters (K-value) from 2 to 8. PLINK 1.9 software [21] was used to generate the input file to run ADMIXTURE analysis [22]. • Effective population size (Ne) for each breed was predicted using the SNeP_111 software [23] based on linkage disequilibrium (LD).

Identification of Runs of Homozygosity (ROHs)
ROHs were computed using the consecutive run method implemented in SVS 8.9 software. We restricted ROH segments to longer than 1 Mb to avoid identification of ROH resulting from LD. In addition, a minimum of 150 homozygous SNPs was required to be in each run, and no heterozygote or missing SNPs were allowed in the ROH, and the maximum gap between SNPs was set to 1000 Kb.
The mean number, the average length of ROH, the average sum of ROH segments per breed, as well as the ROHs distribution into four classes of length (1-2, 2-4, 4-8 and 8-16 Mb) were estimated per breed.
The genomic regions most commonly associated with ROH (ROH_islands) were defined within breed as clusters of runs found in more than 75% of samples using a specific option implemented in SVS software during the ROH detection: samples within the same cluster had identical ROH (length, position on genome and boundaries).

Gene Annotation of ROH_Islands and Functional Analyses
The completed list of turkey genes annotated according to the most recent reference genome was downloaded from NCBI online Database (NCBI Annotation Release: 103). Genes with an official gene name and "LOC" ID associated both with a protein coding and official gene name were catalogued within the identified ROH_islands using the "inter-sectBed" command of BEDTools software [24]. The Database for Annotation, Visualization and Integrated Discovery (DAVID), v6.8, was interrogated to perform a gene ontology (GO) functional annotation and KEGG pathway analyses.
Additionally, considering that to date no quantitative trait loci (QTLs) database is available for turkeys, the one accessible for the chicken species, i.e., the Chicken QTLdb-based on Galgal 6.0 genome assembly, was used to identify-using the "Search by associated gene" option-QTLs overlapping the genes found in the ROH_islands.

Inbreeding Coefficients
Two genomic inbreeding coefficients for all the birds were estimated based on: • F HOM -the difference between the observed and expected numbers of homozygous genotypes within the population calculated using the routine implemented within SVS; • F ROH -the ratio calculated between the total length (sum) of all ROH of an individual and the total length of the autosomal genome covered by the SNP marker dataset, 902,020,024 bp in this study.

Breed Diversity among All Breeds
The proportion of monomorphic SNPs was variable according to breed, spanning from 21% (NeIt) to 77% in ErRo (Table 1). The principal component analysis (PCA) results are shown in the graphs of Figure 1A,B: (i) the first two principal components (PC1 eigenvalue of 18.94-29% of total variation; PC2 eigenvalue of 8.44-16.2% of total variation) allow to clearly distinguish ErRo and BrCI breeds from the other ones clustered in the adjacent space; (ii) PC3 separates the remaining breeds, in particular NeIt birds and Brzl are well separated from PrPc, CoEu and Roma birds. Birds belonging to CoEu and Roma breeds appear to be two very close groups, showing partial overlapping.
Comparable results were found by NJ tree analysis ( Figure 1C) and by the pairwise breed comparisons (F ST values) ( Figure 1D), confirming that ErRo is the more distant breed from the others with an F ST value >0.47 for all comparisons).
The ADMIXTURE analysis revealed that at K = 6 (where the lowest CV error value was identified), six breeds seem to be mostly unique populations, with an ancestral genetic composition defined by a proportion of 90% in CoEu and equal (PrPc) or larger than 98% (Brzl, BrCI, ErRo, NeIt). The Roma breed appeared to have a composite ancestral genetic contribution from CoEu (39%), PrPc (30%), NeIt (7%), Brzl (21%) and BrCI (3%) ( Figure 1E). Figure 2 shows the effective population size trends calculated for each breed up to 950 past generations.

Run of Homozygosity
ROHs were identified in all birds of each breed for a total of 20,858 homozygous regions (Table S1). The largest average number of ROHs per animal was observed for the

Run of Homozygosity
ROHs were identified in all birds of each breed for a total of 20,858 homozygous regions (Table S1). The largest average number of ROHs per animal was observed for the

Run of Homozygosity
ROHs were identified in all birds of each breed for a total of 20,858 homozygous regions (Table S1). The largest average number of ROHs per animal was observed for the CoEu breed (156.2), while the lower ones were identified in PrPc (62.8) and Roma (69.3) ( Table 2). The longer average ROH lengths were observed for Brzl and CoEu breeds (about 2 Mbps) while the lowest were found for the BrCI and Roma breeds (about 1.7 Mbps). The portion of the genome covered by ROH is more than two times larger in ErRo and CoEu compared to the ones for NeIt, PrPc and Roma. The graphical representations of ROH statistics by breeds are shown in Figure 3. In detail, (i) Figure 3A represents the relationship between the number and the mean total length of ROHs for each individual according to breed; (ii) Figure 3B shows the proportion of ROH classified by taking into account the length. As shown, the ROH < 2Mb class of length is the most represented for all seven breeds, whereas the 4-8 Mb and the 8-16 Mb classes of length are the less frequent ones. The Roma breed did not have ROHs longer than 8 Mb. CoEu breed (156.2), while the lower ones were identified in PrPc (62.8) and Roma (69.3) ( Table 2). The longer average ROH lengths were observed for Brzl and CoEu breeds (about 2 Mbps) while the lowest were found for the BrCI and Roma breeds (about 1.7 Mbps). The portion of the genome covered by ROH is more than two times larger in ErRo and CoEu compared to the ones for NeIt, PrPc and Roma. The graphical representations of ROH statistics by breeds are shown in Figure 3. In detail, (i) Figure 3A represents the relationship between the number and the mean total length of ROHs for each individual according to breed; (ii) Figure 3B shows the proportion of ROH classified by taking into account the length. As shown, the ROH < 2Mb class of length is the most represented for all seven breeds, whereas the 4-8 Mb and the 8-16 Mb classes of length are the less frequent ones. The Roma breed did not have ROHs longer than 8 Mb.

Incidences of Common Runs per SNP and ROH_Island Identification
In order to explore how the mating plans in each breed have affected their genetic structure, the ROH_island regions were considered and analysed. Using the SNP incidences, the ROH_islands were defined as a genomic homozygous region with (i) the same boundaries (and as such the same chromosome length) and (ii) shared in at least 75% of samples of a given breed (threshold values are a function of each breed's sample size).
As shown in Figure 4, the genomic distribution of SNP incidence defining ROH_islands is clearly non-uniform across autosomes in all breeds.

Incidences of Common Runs per SNP and ROH_Island Identification
In order to explore how the mating plans in each breed have affected their genetic structure, the ROH_island regions were considered and analysed. Using the SNP incidences, the ROH_islands were defined as a genomic homozygous region with (i) the same boundaries (and as such the same chromosome length) and (ii) shared in at least 75% of samples of a given breed (threshold values are a function of each breed's sample size).
As shown in Figure 4, the genomic distribution of SNP incidence defining ROH_islands is clearly non-uniform across autosomes in all breeds. Table S2 reports the ROH_islands together with the genes annotated in each of them and the overlapping QTLs that were identified interrogating the Chicken QTLdb. A total of 130, 21, 17, 15, 3 and 1 ROH_islands were found in ErRo, Brzl, BrCl, CoEu, NeIt and Roma, respectively. No ROH_islands were identified for the PrPc breed.
ROH_islands were found in ErRo, Brzl, BrCl, CoEu, NeIt and Roma, respectively. No ROH_islands were identified for the PrPc breed. A different proportion of proper (identified only in one breed) and common ROH_islands (shared by at least two breeds) were observed among the populations ( Figure 5A,B). Two ROH_islands were shared by four breeds (BrCI, Brzl, ErRo, Roma-chr 10; Brzl, CoEu, ErRo, NeIt-chr 19) (Table 3). A total of 25 other ROH_islands were common among breeds with various combinations.  Table 3 reports the list of ROH_islands identified in more than one breed, together with annotated genes and the overlapping QTLs. The functional annotation of genes annotated in the ROH_islands (p-value <0.05) is reported in Table S3.  A different proportion of proper (identified only in one breed) and common ROH_islands (shared by at least two breeds) were observed among the populations (Figure 5a,b). Two ROH_islands were shared by four breeds (BrCI, Brzl, ErRo, Roma-chr 10; Brzl, CoEu, ErRo, NeIt-chr 19) (Table 3). A total of 25 other ROH_islands were common among breeds with various combinations.  A different proportion of proper (identified only in one breed) and common ROH_islands (shared by at least two breeds) were observed among the populations ( Figure 5A,B). Two ROH_islands were shared by four breeds (BrCI, Brzl, ErRo, Roma-chr 10; Brzl, CoEu, ErRo, NeIt-chr 19) (Table 3). A total of 25 other ROH_islands were common among breeds with various combinations.  Table 3 reports the list of ROH_islands identified in more than one breed, together with annotated genes and the overlapping QTLs. The functional annotation of genes annotated in the ROH_islands (p-value <0.05) is reported in Table S3.  Table 3 reports the list of ROH_islands identified in more than one breed, together with annotated genes and the overlapping QTLs. The functional annotation of genes annotated in the ROH_islands (p-value <0.05) is reported in Table S3.   When we considered all birds as a unique population, an ROH_island on chr10 was identified ( Figure 6) in at least 75% of turkeys (>n. 136, out of a total of 181 birds). This analysis had the goal of finding a genomic homozygous region shared among all breeds as a possible result of common selection for traits related to survival of the species, such as reproductive efficiency.

Parma e Piacenza ROH_Islands
Although no ROH_islands for the PrPc breed were found at the threshold considered (i.e., 75% of animals of a breed having a specific ROH), two homozygous regions were shared by a large proportion of individuals, 60%, corresponding to 15 birds. These two regions mapped on chr 1 (at 185,658,695-185,769,482) and on chr 21 (at 5,328,921-5,922,971) ( Table 3) and harboured a total of 18 genes (Table S5).

Inbreeding Coefficients
Two genomic inbreeding coefficients FHOM and FROH have been calculated and are reported in Table 4. Figure 7 shows the linear regression and correlation values estimated in each turkey breed. Slightly negative FHOM mean values were calculated for CoEu, NeIt, and PrPc while positive close to 0 ones were found for the remaining breeds. CoEu birds showed the largest variability in FHOM values.   Table S4 reports the allele counts for each of these markers calculated per breed: (i) in all ErRo birds, the homozygous allele was always the alternative one, compared to homozygous birds in other breeds for all these 18 SNPs, without any heterozygous SNPs: e.g., if ErRo turkeys were homozygous AA, all the homozygous birds of Brzl, BrCI, CoEU, NeIt, and Roma were BB. For all these breeds, the individuals with heterozygous SNPs were negligible as shown in Table S4. In the PrPc samples, the distribution of genotypes for this region was about 25% AA, 50% AB and 25% BB.

Parma e Piacenza ROH_Islands
Although no ROH_islands for the PrPc breed were found at the threshold considered (i.e., 75% of animals of a breed having a specific ROH), two homozygous regions were shared by a large proportion of individuals, 60%, corresponding to 15 birds. These two regions mapped on chr 1 (at 185,658,695-185,769,482) and on chr 21 (at 5,328,921-5,922,971) ( Table 3) and harboured a total of 18 genes (Table S5).

Inbreeding Coefficients
Two genomic inbreeding coefficients F HOM and F ROH have been calculated and are reported in Table 4. Figure 7 shows the linear regression and correlation values estimated in each turkey breed. Slightly negative F HOM mean values were calculated for CoEu, NeIt, and PrPc while positive close to 0 ones were found for the remaining breeds. CoEu birds showed the largest variability in F HOM values. As the FROH values were calculated based on ROH proportion over the genome length, they reflect the ROH distribution in each sample, its average length across the classes and the total genome length covered by ROH ( Figure S2). Differences in FROH were also found across the chromosomes of all breeds ( Figure S3). Linear regressions between the two inbreeding values (FHOM and FROH) obtained for each breed are reported in Figure  7 together with correlation coefficients (i.e., R in each graph of Figure 7). Correlation values were found to have a medium value, i.e., from 0.52 to 0.75 in Brzl, CoEu, ErRo and NeIt, to a strong one, i.e., 0.93 in BrCI and 0.96 Roma. As the F ROH values were calculated based on ROH proportion over the genome length, they reflect the ROH distribution in each sample, its average length across the classes and the total genome length covered by ROH ( Figure S2). Differences in F ROH were also found across the chromosomes of all breeds ( Figure S3). Linear regressions between the two inbreeding values (F HOM and F ROH ) obtained for each breed are reported in Figure 7 together with correlation coefficients (i.e., R in each graph of Figure 7). Correlation values were found to have a medium value, i.e., from 0.52 to 0.75 in Brzl, CoEu, ErRo and NeIt, to a strong one, i.e., 0.93 in BrCI and 0.96 Roma.

Genetic Diversity of Breeds
Autochthonous turkey breeds have undergone to a population size contraction after the diffusion of meat type hybrids, more efficient in growth performance and carcass yield. Nowadays, heritage turkey breeds represent a unique genetic biodiversity reservoir. Their in situ conservation is particularly important as the autochthonous breeds own unique characteristics, including adaptability to extensive farming in harsh environments. Their characterization, as envisaged in the TuBAvI project, allows the identification of the phenotypic and the genetic characteristics related to their intrinsic value, especially if compared to selected hybrids that, even though very efficient in meat production performances, may exhibit suboptimal performance in other traits, e.g., reproduction efficiency [17,25].
This study showed that the genetic composition of each of the seven autochthonous turkey breeds differs from others. Based on the PCA and the NJ tree analyses, birds of each breed clustered together revealing a consistent genetic composition. The ErRo and BrCI are the two more genetically distant breeds from the other ones and, according to the F ST values, the ErRo is the one more differentiated from all the other ones: the F ST values were 0.47 and 0.63 when ErRo was compared with PrPc and CoEu, respectively. A possible explanation is that the feather colour of this breed, white with black streaks, clearly different from plumage colour of other breeds, has been chosen as a selection criterion by the breeders and farmers causing, as such, a divergent selection with respect to the other breeds. The pairs Roma-PrPc (F ST = 0.18) and Roma-CoEu were less differentiated (F ST = 0.21). The relative closeness of Roma and PrPc could be related to the geographic origin of these breeds, the northern area of Italy below the Po river in the Emilia-Romagna region (eastern part i.e., Romagna). The Roma breed standard is not as stringent as the one of the other breeds, in fact it does not include a specific trait such as the feather colour as in the ErRo. This may have facilitated the cross of birds of the breeds geographically closer to the Romagna region, i.e., Veneto (CoEu), Emilia (PrPc) and Lombardia (Brzl), but not the ErRo due to its colour. It is interesting to note that the proportion of ancestors found in the Roma breed using the ADMIXTURE analysis was inversely proportional to the geographical distance of Veneto, Emilia and Lombardia from Romagna. With the exclusion of the Roma breed, the ADMIXTURE analysis revealed for all the other ones the uniqueness of the genetic background for each of the Italian breeds with a proportion of each ancestor >90%.

ROH and ROH Islands
In the last decades, the in situ maintenance, i.e., the reproduction of Italian autochthonous turkeys, was mostly carried on by private breeders. They had two main selection goals: (i) to maintain the morphological characteristics of each breed; (ii) to improve birds' performance in a semi-extensive farming system.
The meat of these breeds is mainly used for traditional cuisine recipes for private consumption; the backyard farming as a consequence is not centred on growth performance but on meat quality (i.e., muscle fibre consistency for long cooking) and resilience of birds to pastoral extensive farming.
Generally, the reproductive scheme used by private breeders is avoiding as much as possible increases in inbreeding, adopting as strategy of toms' exchange between farms. Although the selection in autochthonous breeds is not as focused and intense as the one in other livestock species in pure breeding, it determined consequences in the genome makeup of the turkey breeds, including the formation of ROH_islands.
The population size of each breed is particularly small because its contraction occurred in the last 4 decades: assuming a generation interval of 1 year (a realistic value due to the seasonality of reproduction in autochthonous breeds), the graph in Figure 2 shows the decrease of Ne in the last 45 generations for each breed.
It is interesting to note that the Ne of the ErRo was the smallest among all, with a fairly constant value in the last 30 generations. This is a further evidence supporting the closeness of the nucleus of the ErRo population in the last decades. The Roma breed, in accordance with the ADMIXTURE results, exhibited the largest Ne possibly because of the inclusion of birds from different breeds.
The population size of these breed is very small as reported by the FAO DAD-IS database and by Castillo et al. ranging from 9 (PrPc) to 445 (BrCI) birds across the breeds and all considered at critical risk of extinction [26,27]. Due to the very limited population size, the total number of homozygous SNPs was then found to be quite high in all the breeds, with the higher proportion in the ErRo, where the total number of homozygous SNPs was 318,213 (92% of the total SNP dataset). As discussed hereinbefore, the ErRo breeders aim at a strict maintenance of the morphological standard of the feather colour, limiting the number of reproducers and the variability existing in the breed. This might be one of the reasons determining a higher homozygosis in this breed with respect to the other ones having a proportion of homozygous loci ranging from 72% to 84%. The ErRo also had the highest proportion of SNPs in ROH_islands suggesting, together with other results here presented, that the selection for the specific colour of plumage might have played a role in building a larger proportion of ROH_islands. The ErRo, in particular, is the only breed showing a clear homozygosity state and ROH_island for all individuals (see Supplementary Materials, Figure S4) on chr1: 41,354,408,938 that harbours the KITLG gene, associated with melanin patterning [28] and involved in skin coloration in vertebrates [29]. Additionally, the ErRo had the highest number of ROH_Islands, 130 vs. 21, 17, 15, 3 and 1 found in Brzl, BrCI, CoEu, NeIt and Roma, respectively. However, the proportion of ROH_Islands overlapping with those identified in other breeds was only 17% (mainly with BrCI, n. 12), suggesting that a large proportion of regions is under specific selection pressure in the ErRo.
The Roma breed is the one with the smallest proportion of homozygous SNPs in ROH. Additionally, there are no ROHs longer than 8 Mb. This result agrees with the ADMIXTURE one, indicating a possible introgression in the breed from various ones as hereinbefore supposed.
A similar high level of homozygosity has been identified in Italian local chicken breeds [30], where low values of heterozygosity underline the difficulties of maintaining a fair level of genetic variability in breeds under an in situ conservation plan.
The ROH_Island located on chr10 is shared by the highest number of birds across breeds, i.e., n. 155 (n. 136 turkeys is considered the threshold of 75% of samples), as shown in Figure 5. In this ROH_Island, 200 consecutive SNPs were homozygous in all ErRo and in the major part of turkeys belonging to the other breeds; 18 homozygous SNPs mapped in the intronic position of the PRGS2 and PLA2G4A genes ( Figure 6).
To assess the genetic variability existing in this species, it is important to characterize the heritage breeds' specific genetic biodiversity, especially in comparison to the highly selected heavy turkey lines. Comparing the evidence of the ROH_Islands in this region to existing data on ROH in commercial hybrids, it is interesting to note that the PTGS2 and PLA2G4A genes were found by several authors to be related to reproductive physiology in avian and other livestock species [31][32][33][34][35][36]. The wide occurrence of this ROH in a large number of animals of autochthonous turkey populations under an outbreeding reproduction scheme suggests that the genes included in this region may be under selection because they are important for reproduction and thus for the survival of the species. NeIt and CoEu are two breeds that in the recent past were appreciated and bred for their good egg production. Except the Roma, the PrPc and one bird of CoEu (heterozygous), all other birds were in a homozygous state in and in the proximity of the two genes here mentioned. The ErRo was in a homozygous state but with the alternate allele with respect to the other breeds (i.e., Brzl, BrCI, NeIt, CoEu) as shown in Supplementary Materials, Table S4. Comparing the results of this study with the ROHs mapped in a commercial turkey hybrid, this region does not appear to be in a homozygous state in the commercial selected, fast-growing, heavy turkey line as reported by Strillacci et al. [11]. The intense selection for body weight in commercial turkeys determined a reduction in reproductive performance, as reported by Nestor et al., who discussed the genetics of growth and reproduction using the performance measured in over 40 generations of selection for 16-week body weight in turkey [37].
Interestingly, the PrPc, having the largest proportion of heterozygous individuals in this region (i.e., f(A) = f(B) = 0.50), is also the heaviest of all the Italian breeds [8,9] reaching an adult weight doubling the one of the other breeds (about 12 kg in males vs. 5-7 kg of other breeds; about 6.5 kg in females vs. 3-4 kg in other breeds). At present, the studies on genomic structure in turkeys are very limited as also recently pointed out by Adams et al. [38] and, to best of our knowledge, no other published data are available on ROHs in hybrid turkeys or in selected pure lines.

Comparison with Literature
To the best of our knowledge, only three studies [10,11,38] performed a genomewide ROH detection in local (Mexican), hybrid and pure-line turkey populations. Among the ROHs detected by these authors, only those found by Strillacci et al. are available online: among these published ROHs, only two overlapped with those here identified (in Brzl breed), both on chr 8 at 4,344,134-5,525,301 and at 5,536,601-6,787,697 [11]. In these two regions, 15 genes are annotated, some of which are related to meat juiciness and tenderness (MMRN2 [43,44]), fat cell development and fatty acid metabolism (ADIRF [45]), and circadian rhythms (Opn4 [46]).

Inbreeding Coefficients
The possibility of calculating genomic inbreeding values is useful in the absence of pedigree information, a very common condition in avian species. This is particularly important in local and endangered populations such as autochthonous Italian turkey breeds. In our study, populations showed genomic inbreeding with a large variability of mean values as shown in Table 4: F HOM mean values range from −0.118 to 0.141 while F ROH mean values vary from 0.126 to 0.401. The variability of inbreeding values within breeds is also very different across them: considering F ROH , the SD of the mean values ranges from 0.035 to 0.153 showing that the effect of selection and the management of reproductive practices breeders are influencing in a variable manner the genetic makeup of these populations.
The two breeds with lower value of F ROH are the Roma and the PrPc, reflecting the reproductive structure and management of the two populations. The Roma breed is an admixture of different turkey breeds as hereinbefore commented. The PrPc is, together with the Roma, the breed with largest Ne as reported in the graph of Figure 2. Additionally, the observed heterozygosity in the PrPc is the highest among all. The size reduction of the populations during the last decades is for sure having an impact on the genomic inbreeding value: the number of farms and the frequency of exchange of birds among breeders may also affect inbreeding values. Unfortunately, up to 2014, with the setup of the RAA, there is no track of the reproducers used in the populations, making the genomic approach a unique resource for mapping the existing variability among breeds.
A positive correlation was found between F HOM and F ROH in all breeds, with R values ranging from moderate to strong (>0.90 in BrCI and Roma). These results, attributable to a higher proportion (from 62% to 77%) of short ROHs in all breeds, are consistent with those reported in the literature on livestock [11,[47][48][49] and confirm the helpfulness of ROHs in the estimation of the inbreeding coefficient, mainly in small populations where pedigree information is difficult to collect or not available.

Conclusions
Information on the genetic variability of turkey breeds is very limited and, to the best of our knowledge, this is the first study on autochthonous turkey populations. The wide diffusion of commercial, fast-growing, heavy hybrids in the last 40 years has impacted the farming of autochthonous populations, affecting in the same manner the existing genetic variability in the turkey species. The wide spread of commercial hybrids has in fact reduced the farming of heritage breeds, causing a loss of genetic variation and biodiversity. Local breeds are well adapted and selected over centuries to deal with environmental harsh conditions of semi-extensive farming system. The knowledge of their genomic variability is indeed an important step forward with respect to the state of the art, today mainly based on selected heavy turkeys. This study is in fact part of the nation-funded project TuBAvI that aims to investigate the existing biodiversity in local avian populations using genomic data and to provide insights for their in situ conservation. This study releases previously unavailable information on the genomic variation existing in autochthonous turkey populations, showing regions under selection. The results of this project represent as such a unique resource to compare the genetic variability of autochthonous turkey breeds with the one of selected commercial turkeys. The use of ROHs to exploit genomic variation in autochthonous populations revealed genomic regions under selection (ROH_islands) that harbour genes potentially interesting also for a targeted selection in commercial lines. The results here obtained are, as such, a unique resource in mapping genomic variation in turkey species.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/genes12091342/s1, Figure S1: Box plot of F HOM values calculated for each breed, Figure S2: Box plot of F ROH values for each class of ROH length; Figure S3: F ROH values calculated for each chromosome; Figure S4: ROH_island on chr1 containing KITGL gene: homozygosity genotypes for all ErRo's birds. Table S1: List of ROH identified for all turkeys; Table S2: List of ROH_Islands per breed; Table S3: Annotation of genes performed using DAVID database; Table S4: SNPs mapped in PTGS2 and PLA2G4A genes and turkey genotypes; Table S5: ROH_Islands identified in Parma e Piacenza (PrPc) breed.