The Genomic Variation in the Aosta Cattle Breeds Raised in an Extensive Alpine Farming System

Simple Summary Genetic variability among native cattle breeds can disclose the important features that make a population adapted to harsh environments. The Aosta cattle breeds have been raised and selected for centuries to be farmed in a mountain environment, characterized by a semi-intensive system, i.e., summer pasture with winter recovery on the farms. To disclose the genomic variation and its association with known genes, it is important to genetically characterize these breeds. Abstract The Aosta Red Pied (Valdostana Pezzata Rossa (VRP)), the Aosta Black Pied (Valdostana Pezzata Nera (VBP)) and the Aosta Chestnut (Valdostana Castana (CAS)) are dual-purpose cattle breeds (meat and milk), very well adapted to the harsh environmental conditions of alpine territories: their farming is in fact characterized by summer pasture at very high altitude. A total of 728 individuals were genotyped with the GeenSeek Genomic Profiler® (GGP) Bovine 150K Illumina SNP chip as a part of the DUALBREEDING-PSRN Italian-funded research project. The genetic diversity among populations showed that the three breeds are distinct populations based on the FST values, ADMIXTURE and Principal Component Analysis (PCA) results. Runs of Homozygosity (ROH) were obtained for the three populations to disclose recent autozygosity. The genomic inbreeding based on the ROH was calculated and coupled with information derived from the F (inbreeding coefficient) and FST parameters. The mean FROH values were low: CAS = 0.06, VBP = 0.05 and VRP = 0.07, while the average F values were −0.003, −0.01 and −0.003, respectively. The annotation and enrichment analysis, performed in the identified most frequent ROH (TOP_ROH), showed genes that can be linked to the resilience capacity of these populations to harsh environmental farming conditions, and to the peculiar characteristics searched for by farmers in each breed.


Introduction
Animal genetic resources play an important role in local economies and in maintenance of territories and landscapes [1]. Among the large number of autochthonous cattle populations in Italy, the Aosta breeds play an important role for the Aosta valley, located in the northwest Alpine territories of Italy [2]. In addition to their milk and meat production, their economic value is also related to the farming activity itself, closely linked to the use of local territories: farming activity, in fact, allows to Analysis and ADMIXTURE analysis. Finally, we aimed to annotate the genes mapped in the ROH in order to disclose the common and proprietary regions under selection in the three Aosta breeds.

Ethics Statement
This study did not require approval from the Animal Care and Use Committee.

Sampling and Genotyping
The Associazione Nazionale Allevatori Bovini di Razza Valdostana (A.N.A.Bo.Ra.Va.) provided genotypes of 728 individuals-male and females (CAS, 297; VBP, 153; and VRP, 278)-obtained with the GGP Bovine 150K Illumina SNP chip. The sampling of the individuals was homogeneous across breeds regarding their population structure and reflects their actual population size. The genotyping is part of the national project DUALBREEDING, a pluriannual effort funded by the EU EAFRD and by the Ministry of Agriculture of Italy aimed at maintaining the biodiversity in cattle populations. The SNP genotypes were mapped on the ARS-UCD1.2 bovine reference genome. Out of the genotypes available on the SNP chip, those with a sample and marker call rate ≤0.90, without chromosomal position and on non-autosomal chromosomes, were deleted, leaving a total of 128,180 SNP markers for subsequent analyses.

Aosta Breeds Diversity Performed by PCA, F ST and ADMIXTURE
Genetic diversity within and among breeds was determined through a Principal Component Analysis (PCA) and estimating the pairwise Fixation Index through the pipelines implemented in Golden Helix (SVS) 8.8.4 software (Golden Helix Inc., Bozeman, MT, USA). The Fixation Index-equivalent to Wright's F-statistic F ST -is a measure of the total genetic variance that can be explained by the population structure; in other words, F ST estimates the genetic divergence among subpopulations. The F ST was estimated for three pairs of breed combinations (CAS vs. VBP, CAS vs. VRP, and VBP vs. VRP). The expected (He) and observed (Ho) heterozygosity values were also calculated for each breed.
Determination of the most probable number of ancestral populations was obtained with ADMIXTURE v. 1.3.0 software [15]. ADMIXTURE was run from K = 2 to K = 6 and the optimal number of clusters (K-value) was determined as the one having the lowest cross-validation error. The R script suggested by the ADMIXTURE procedure was used to perform a graphical representation of the ADMIXTURE results. PCA, ADMIXTURE and Wright's F-statistic were performed on the pruned SNP dataset (83,776 markers, resulted by using r 2 > 0.5 in a 50-SNP sliding windows in SVS software) in order to reduce the impact of the SNP ascertainment bias from linkage between loci.

F ST Analysis at Marker Level
In order to identify the genome-wide patterns under selection, the outlier loci approach, based on the F ST estimated for each SNP (marker-based F ST ), was also used and obtained with SVS on the same three breed pair combinations. Genomic regions can be considered being under different (positive) selection in a pair comparison if they contained a high proportion of highly differentiated SNPs based on the F ST values (deviation from neutral loci with an F ST threshold of 0.5) [16]. Differentiated SNPs with an F ST > 0.5 were considered to be under different selection.

Runs of Homozygosity Detection
ROH analyses were performed separately for each population, using SVS software. No linkage disequilibrium (LD)-based pruning was performed and, as in [17], the minimum ROH length was set to 1 Mb to avoid the detection of short and common ROH across the genome due to LD. The ROH were defined setting a minimum of 1000 Kb in length and 60 homozygous consecutive SNPs. In addition, no heterozygote SNPs and no missing SNPs were allowed in the ROH, and a maximum gap between Animals 2020, 10, 2385 4 of 18 the SNPs of 1000 Kb was set in order to assure that the SNP density did not affect the ROH. ROH were also grouped into 5 classes of length in each set of ROH: <2 Mb, 2-4 Mb, 4-8 Mb, [8][9][10][11][12][13][14][15][16] Mb and >16 Mb. Descriptive statistics of the ROH were calculated across individuals in each Aosta breed. The genomic regions with the highest frequency of ROH (TOP_ROH) were identified by selecting the top 1% SNPs with the largest occurrence (TOP_SNPs).

Gene Functional Analysis
The full gene set (Bos taurus: Annotation Release 105) was downloaded from NCBI [18] and genes were catalogued within the TOP_ROH using the intersectBed command of BEDTools [19]. In addition, the position of all TOP_SNPs with respect to the annotated genes in the TOP_ROH was identified. The SNPchiMp online database [20] was used to convert the Illumina SNP name to the SNP rsID, the unique SNPs code recognized by the Ensembl Variant Effect Predictor (VEP) [21], which was employed to annotate all the TOP_SNPs. Only genes with official "gene name ID" and LOC genes associated with a protein-coding gene name (excluding uncharacterized ones) were considered. A gene ontology (GO) functional annotation and KEGG pathway analyses were performed using DAVID Bioinformatics Resources, version 6.8 [22].
In addition, bovine QTL, available from the "AnimalQTLdb" database (Cattle-ARS-UCD1.2 genome assembly), were catalogued into the TOP_ROH by overlapping [23]. The same SNP annotation approach was applied for differentiated SNPs with an F ST > 0.5.

Inbreeding Coefficients
Identity-by-descent (IBD) was obtained in order to assess the quality of the dataset, to identify the potential sample replicates and to estimate pair-wise comparison of relatedness (comprising first-degree relatives). IBD was calculated using SVS software on pruned SNPs. SVS software also provided the individual's inbreeding coefficient F [24], which was calculated for each breed. The F values ranged from −1 to +1, representing an excess of heterozygosity and an excess of homozygosity, respectively, where an F value equal to 0 denotes the Hardy-Weinberg equilibrium across all markers. To compare inbreeding across breeds, the average F value was calculated for each breed.
In addition, the inbreeding coefficient based on the ROH (F ROH ) for each sample and for the three breeds was calculated separately, considering the following formula [8]: where L ROH is the total length of all the proper ROH of an individual, and L AUT is the specified length of the autosomal genome covered by SNPs (2,487,082,459 bp).

Aosta Breeds' Diversity
The effective number of polymorphic SNPs (number of SNPs in which at least one heterozygous individual was identified) represented up the 94% of the total SNPs of all three breeds. The He and Ho values were similarly low among breeds: CAS (0.344 and 0.342), VBP (0.345 and 0.341) and VRP (0.339 and 0.338). The Principal Component Analysis (PCA) result is displayed in Figure 1A. On the first principal component (PC_1) (eigenvalue of 16.46), a clear separation of VRP from CAS and VBP is observable. On the other hand, the PC_2 principal component (eigenvalue of 5.12) allowed to separate the CAS and VBP breeds. Samples belonging to these two breeds appear to be two very close groups, showing partial overlapping.
These results are confirmed by the pairwise breed comparisons using F ST , showing that although there is a very low level of genetic differentiation, there is a clear and similar distinction between  To investigate the ancestry composition of the Aosta breeds, ADMIXTURE analysis was run for values of possible ancestors (K) ranging from 2 to 6. The lowest CV value was obtained at K = 3 ( Figure 1B). At this K value, all three breeds appear to be mostly unique populations, represented by an ancestral genetic group in a proportion larger than 90% ( Figure 1C). At K = 2, CAS and VBP belonged to a unique ancestral genetic group, while VRP to a different one.

F ST at Marker Level, Gene Annotation and Gene Functional Analyses
The distribution of F ST values, calculated for each marker across all chromosomes for each pair of comparison, is shown in Figure S1. We consider the SNPs' allelic frequency as differentiated, and then under a possible different positive selection, when the F ST value is >0.5 (above the redline threshold in Figure S1). A total of 53 differentiated SNPs was from the VRP_VBP breeds' comparison. Forty-one SNPs were mapped in the intergenic regions and 12 in the intronic regions (CORIN, FIP1L1, LNX1, CLOCK, PIEZO1, CBFA2T3, TMEM104 and DNAJC12; Table S1, Sheet 1). If we consider the VRP_CAS breeds' comparison, 34 SNPs were differentiated markers, of which 29 were annotated in intergenic regions, one in the 3_prime_UTR_variant region of PDGFRA gene and four in the intronic positions of the KIT (n = 1) and CLOCK (n = 3) genes (Table S1, Sheet 2). No differentiated SNPs were identified for the CAS_VBP breeds' comparison. There were no significant GO terms and KEGG pathways from the DAVID database for the genes mapped in the differentiated regions.

Runs of Homozygosity Detection
The ROH were identified in all 728 individuals of the three Aosta breeds, for a total of 36,400 homozygous regions. At the individual level, the average numbers of ROH per animal were 50, 45 and 53 for CAS, VBP and VRP, respectively (Table 1), with a total mean ROH length of 2.87, 2.93 and 3.15 Mb, as shown in Figure 2A, which graphically represents the relationship between the ROH counts and the averaged total length of the ROH for each individual.  Differences among animals were also found considering the total length of the genome covered by the ROH (sum of all ROH per animal): 6.10-381.55 Mb for CAS, 19.02-403.27 Mb for VBP and 7.31-557. 46 Mb for VRP samples (data not shown).
The ROH were identified for all classes of length and, generally, the lower coverages were in concordance with a low number of regions per samples. Although shorter regions (<2 Mb) are the most frequent classes of length (about 50%; Figure 2B), the proportion of the genome covered by them was relatively small (averaged total length per samples around 1.49 Mb) in all three breeds. Contrariwise, accounting for a small number of ROH per sample, an ROH larger than 16 Mb (CAS and VBP, from 1 to 6; VRP from 1 to 9) covered a wider region of the genome: an ROH >16 Mb was identified in a total of 26%, 30% and 34% of the CAS, VBP and VRP samples, covering up to 5.8%, 6.25%, and 11.3% of their autosome genome.
An ROH was found on all chromosomes and no evident relationship between the chromosomes' length and mean ROH length was observed: the graphical representation of the ROH frequencies on autosomes together with the mean ROH coverage length calculated for each chromosome is shown in Figure S2. The mean length values of the ROH mapped on autosomes in CAS were more uniform with respect to the ones calculated for VBP and VRP.
In order to explore the effect of selection on the Aosta breeds' genome, the TOP_ROH regions were considered. The TOP_ROH regions are those located above the redline threshold, as shown in Figure 3 for each breed. The SNP occurrences defining the thresholds were 48, 24 and 54 for CAS, VBP and VRP, respectively. As shown in Figure 3, the genomic distribution of TOP_ROH is clearly non-uniform across autosomes for all three cattle breeds.
TOP_ROH and the annotated genes are reported in Tables 2-4 for CAS, VBP and VRP, respectively. A total of 18 TOP_ ROH on 9 chromosomes (BTA) were identified in CAS (above the redline in Figure 3), as reported in Table 2. The higher chromosomal peaks were identified on BTA5, BTA6, BTA19 and BTA28.  TOP_ROH and the annotated genes are reported in Tables 2-4 for CAS, VBP and VRP, respectively.  One thousand three hundred and fifty-five SNPs are involved in the definition of TOP_ROH and the summary indication of their annotated position on genome are shown in Figure S3. According to the Ensembl VEP tool, the SNPs are mapped in the intergenic (57.5%) and intragenic (42.5%) positions ( Figure S3). If we consider the intragenic SNPs, the highest number of homozygote SNPs (≥25) was annotated within the KHDRBS2 (n = 33), CTNNA3 (n = 32), ADGRL3 (n = 27), CCSER1 (n = 26) and CACNA2D1 (n = 25) genes. In addition, three homozygote SNPs were annotated in missense positions of the PCDHA13 (Hapmap42803-BTA-110000-BTA7), KRBA2 (BovineHD1900008383-BTA19) and DNA2 (BTB-00981633-BTA28) genes.
The proper TOP_ROH of VRP were 11, identified on 8 chromosomes (Table 4 and above the redline in Figure 3). SNPs delineating the TOP_ROH were 1304, annotated both in the intergenic (69.3%) and intragenic (30.8%) positions ( Figure S3). Among the latter, the three missense SNP positions were in the MAP3K19 (BovineHD0200018075-BTA2), CFAP221 (ARS-BFGL-NGS-16745-BTA2) and GALNT6 (ARS-BFGL-NGS-110943-BTA5) genes. The highest number of homozygous SNPs were harbored within the THSD7B (n = 48) and the KHDRBS2 (n = 33) genes.  Within the CAS, VBP and VRP's TOP_ROH, a total of 312, 212 and 162 genes were annotated, respectively. Among them, 6 genes (C5H12orf50, C5H12orf29, CEP290, TMTC3, KITLG and KHDRBS2) lied within the three regions in common among the three breeds (CAS-VBP-VRP). In addition, other 12 TOP_ROH were shared by two breeds: n = 5, CAS-VBP; n = 4, VBP-VRP; and n = 3, CAS-VRP.  Figure 4 represents the Venn diagram of the proprietary and shared genes annotated in the TOP_ROH; for the common regions, the lists of genes are also reported. For CAS and VBP, a larger proportion of the shared genes was identified (n = 80, including the genes in common with VRP) with respect to CAS_VRP and VBP_VRP, while the lowest one was between CAS and VRP (n = 11). Not all common TOP_ROH harbored genes. The annotation, performed with the DAVID database, is presented in Supplementary Table S2.
The Table 5 reports details of the QTL associated with bovine traits (within the TOP_ROH) according to the nomenclature available in the AnimalQTLdb database [23].

Inbreeding Coefficients
No sample duplication (pairwise IBD > 0.95) was identified for each breed. However, some samples belonging to CAS, VBP and VRP, with a genetic similarity greater than 50%, were found, suggesting a first-degree relationship between pairs of individuals ( Figure 5A). The average genomic relationship within breed was similar in all populations, with a slightly higher value in VRP ( Figure 5B). The average genomic relationship among breeds was higher between CAS and VBP. The VRP samples had a very low relationship with CAS and VBP. Figure 4 represents the Venn diagram of the proprietary and shared genes annotated in the TOP_ROH; for the common regions, the lists of genes are also reported. For CAS and VBP, a larger proportion of the shared genes was identified (n = 80, including the genes in common with VRP) with respect to CAS_VRP and VBP_VRP, while the lowest one was between CAS and VRP (n = 11). Not all common TOP_ROH harbored genes. The annotation, performed with the DAVID database, is presented in Supplementary Tables S2. The Table 5 reports details of the QTL associated with bovine traits (within the TOP_ROH) according to the nomenclature available in the AnimalQTLdb database [23].   The inbreeding coefficient, estimated using the proportion of homozygous SNPs distributed overall the autosomes SNP markers (F), were slightly negative (close to 0) in each population (min, max and mean values in Table 6). The average observed homozygote genotypes (67% in each breed), the expected homozygote ones and F are similar in CAS and VBP, and higher in VRP. Individual FROH values varied from 0.002 (CAS) to 0.224 (VRP), with the highest average value for FROH (0.067) in the VRP ( Table 6).
As shown in Figure 6, the FROH values, calculated within each class of ROH length, differ among them: VRP samples showed clearly higher average FROH values in the three shortest classes of ROH length. FROH values reflect the ROH distribution and its average length across the classes. According to the same principle, differences in FROH were also found along all chromosomes ( Figure S4). The inbreeding coefficient, estimated using the proportion of homozygous SNPs distributed overall the autosomes SNP markers (F), were slightly negative (close to 0) in each population (min max and mean values in Table 6). The average observed homozygote genotypes (67% in each breed), the expected homozygote ones and F are similar in CAS and VBP, and higher in VRP.
Individual F ROH values varied from 0.002 (CAS) to 0.224 (VRP), with the highest average value for F ROH (0.067) in the VRP ( Table 6).
As shown in Figure 6, the F ROH values, calculated within each class of ROH length, differ among them: VRP samples showed clearly higher average F ROH values in the three shortest classes of ROH length. F ROH values reflect the ROH distribution and its average length across the classes. According to the same principle, differences in F ROH were also found along all chromosomes ( Figure S4).

Discussion
The three Aosta breeds have been sharing the same environment and farming practices for centuries. Nowadays, after the structuring of the breeding activities in herd books, even if they are part of the same herd book association, they are managed as three different populations with different selection indexes. A common denominator is that their milk is used for the Fontina cheese production, a DPO product that, in its manufacturing specifications, includes the rule that only the milk from Aosta cattle breeds can be used.
In term of selection, CAS and VBP are sharing a similar selection goal accounting for milk, meat and fighting ability. Differently, the dual-purpose selection of VRP is more oriented toward milk production [3]. The differentiation in selection goal occurred for decades and this may reflect the findings in term of genomic regions under selection that may be identified as the ROH.
Even if the three breeds share the same environment and farming structure, they appear to be genetically different. The genetic relationship between CAS and VBP, previously reported using microsatellite markers [25], was confirmed here with the use of SNP markers, as per the FST (both at population and at marker levels), PCA and ADMIXTURE results.
The FST statistic at the population level showed a value of genetic differentiation among the three breeds, which was around 0.05 when the comparison (both) involved VRP, but a lower one (FST = 0.019) between CAS and VBP. At the single-marker level, no differentiated SNPs (FST > 0.5) were identified in the CAS_VBP breed comparison, but several ones were found for the comparison of VRP with the other two breeds. These different values of FST can be affected by the origin of the three Aosta breeds: as known by historical evidences, the origin of VRP is independent from VBP and CAS. In fact, historical findings indicate that VRP was introduced in the area by the Burgundians in the 5th century AC.
It is interesting to mention that these genomic regions are harboring genes involved in the mammalian circadian rhythms regulation (CLOCK) [26], feed efficiency and growth traits (CORIN) [27], marbling (LNX1) [28] and immune response to mammary gland inflammation (CBFA2T3) [29]. It is worth noting that the VRP_CAS comparison identified the KIT and PDGFRA genes involved in melanogenesis and in coat color (spotting) (KIT) [30], as well as in intramuscular adipocyte development and marbling fat deposition (PDGFRA) [31], respectively. The coat color and shape

Discussion
The three Aosta breeds have been sharing the same environment and farming practices for centuries. Nowadays, after the structuring of the breeding activities in herd books, even if they are part of the same herd book association, they are managed as three different populations with different selection indexes. A common denominator is that their milk is used for the Fontina cheese production, a DPO product that, in its manufacturing specifications, includes the rule that only the milk from Aosta cattle breeds can be used.
In term of selection, CAS and VBP are sharing a similar selection goal accounting for milk, meat and fighting ability. Differently, the dual-purpose selection of VRP is more oriented toward milk production [3]. The differentiation in selection goal occurred for decades and this may reflect the findings in term of genomic regions under selection that may be identified as the ROH.
Even if the three breeds share the same environment and farming structure, they appear to be genetically different. The genetic relationship between CAS and VBP, previously reported using microsatellite markers [25], was confirmed here with the use of SNP markers, as per the F ST (both at population and at marker levels), PCA and ADMIXTURE results.
The F ST statistic at the population level showed a value of genetic differentiation among the three breeds, which was around 0.05 when the comparison (both) involved VRP, but a lower one (F ST = 0.019) between CAS and VBP. At the single-marker level, no differentiated SNPs (F ST > 0.5) were identified in the CAS_VBP breed comparison, but several ones were found for the comparison of VRP with the other two breeds. These different values of F ST can be affected by the origin of the three Aosta breeds: as known by historical evidences, the origin of VRP is independent from VBP and CAS. In fact, historical findings indicate that VRP was introduced in the area by the Burgundians in the 5th century AC.
It is interesting to mention that these genomic regions are harboring genes involved in the mammalian circadian rhythms regulation (CLOCK) [26], feed efficiency and growth traits (CORIN) [27], marbling (LNX1) [28] and immune response to mammary gland inflammation (CBFA2T3) [29]. It is worth noting that the VRP_CAS comparison identified the KIT and PDGFRA genes involved in melanogenesis and in coat color (spotting) (KIT) [30], as well as in intramuscular adipocyte development and marbling fat deposition (PDGFRA) [31], respectively. The coat color and shape (uniform and chestnut in CAS and red and white in VRP) is among the distinctive criteria of the breeds, whereas the differences in fat deposition can be ascribed to the peculiar constitution of the CAS breed. The vitality and attitude for dominance of this breed are a part indeed of a pretty "masculine" phenotype of CAS cows, showing curly hair and a large neck with developed anterior muscle masses, like the neighboring Hérens breed [32,33].
The genetic distinctness of the three populations is supported also by the PCA, showing clearly that the three breeds cluster separately: VBP and CAS are both placed on the same spatial position according to PC_1 (Figure 1A), explaining 41% of the total variance, and both separated from VRP. PC_2 is on the other hand separating the two breeds in two differentiated groups. The ADMIXTURE results also clearly support the evidence that the three breeds have distinct genetic origins.
The ROH analysis is particularly interesting in these three populations as it allows to disclose recent inbreeding caused by the managing of the populations, for example, by performing artificial insemination and with structured breeding plans, even if the selection scheme was carefully evaluated and planned to minimize the increasing of inbreeding [4].
The recent inbreeding's loop produces small numbers of long ROH, influencing the sum of ROH much more than the total number of ROH itself. In the VRP, the regions linked to recent inbreeding are more frequent than the ones found in the other two breeds. In VRP, in fact, an ROH > 16 Mb has been found in 34% of samples (of which 28 animals had an ROH longer than 30 Mb), while in CAS and VBP, it was found in 26% and 30% of animals, respectively. Among the latter, 15 CAS and 12 VBP individuals had an ROH longer than 30 Mb.
Although a similar proportion of homozygous SNPs was identified among the three breeds (around 67%- Table 6), a higher number of homozygous SNPs concentrated in the ROH was detected in VRP (10.12%-calculated as the number of SNPs defining ROH/observed homozygotes). In CAS and in VBP, the proportion of homozygous SNPs mapping within the ROH is lower: 8.8% and 8.04%, respectively.
In VRP, a higher number of ROH and a longer size with respect to CAS and VBP were identified, mainly as a consequence of no introduction of animals from other regions. This is not the case for CAS, who has been recently recognized as genetically similar to the Hérens breed. The possibility to enroll progeny of the two breeds in the studbook of any of them was recently approved by the two breeder's associations. This recent advance in reciprocal recognition occurred after some generation of known exchange of reproducers. In fact, already in 1929 [34], the closeness of CAS with Hérens was highlighted and the author already raised the hypothesis that the two breeds could be recognized genetically as one population. A previous study on microsatellites, carried out on the most important cattle breeds of the Alpine arc [25], also support the strong relationship between these two breeds. Before the approval of using Hérens bulls for breeding, crossbred mating was, however, already occurring between the CAS and VBP breeds, contributing to maintain the average level of inbreeding lower than in VRP (both F and F ROH ).
The ROH are not randomly distributed across the genome and there are regions with a high prevalence of ROH. In the CAS breed, on the BTA5 located around 79-80 Mb, the highest number of samples (n = 94-31%) that shared the same TOP_ROH was identified. This TOP_ROH, defined by 18 homozygous SNPs (of which 5 are intronic), lies within the TECRL gene encoding for an enzyme involved in chemical reactions and pathways involving lipids, also reported to possibly play a role in puberty and female fertility in cattle [35]. In CAS, within the TOP_ROH on BTA4 11, the genes belong to the Homeobox family genes, of which HOXA13, HOXA11, HOXA10, HOXA9, HOXA7, HOXA5, HOXA3 and HOXA4 are involved in the reproductive tract and in development and fertility in males and females [36][37][38]. It interesting to note that a negative genetic relationships occurs between fighting ability and fertility, as observed in CAS [33], and a reduction in fertility has been found in Hérens cows, also empirically, for fighting ability across time [39]. The breed has no problems with fertility, like all Aosta breeds, well known for their hardiness, but a strong selection towards fighting ability; disregarding fitness characteristics could, in the long term, have a detrimental effect on fertility.
The second highest peak (TOP_ROH) found in CAS is located on chromosome 5 (max number of samples = 78, 25.4%), where the KITLG maps. This gene was included in regions under selection in cattle breeds [40,41] and is responsible for the coat color phenotype in different species [42,43]. In VBP, two higher peaks were identified on BTA6 and BTA23, where the maximum number of samples (n = 59, 38.5%) and (n = 45, 29.4%), respectively, shared intergenic and intragenic (KHDRBS2 gene) homozygous regions. The KHDRBS2 gene has been associated with fertility traits in goats [44] and in Brahman cows [45], and more recently with adaptability traits in Colombian cattle breeds [46]. This finding may relate to the VBP as well as other Aosta breeds, which are well adapted to harsh environmental conditions, such as the alpine pastures farming. Lastly, the regions shared by a higher number of VRP samples (more than half) are those found on chromosomes 2 (n = 150, 53.6%), 5 (n = 157, 56%) and 6 (n = 180, 64.3%). Except for the most representative regions on BTA2, for which an SNP maps in the intronic position of the DARS1 gene, the other ones include SNPs lying closely to the KITLG (BTA5) and KIT (BTA6) genes, respectively.
We would like to underline the impact of selection on highly homozygous regions: several interesting genes in TOP_ROH include a large number of SNPs annotated within the gene itself. Although this information may be affected by the density of SNPs on the chip and by the length of the genes, it would be worth mentioning them. The THSD7B gene of the VRP's TOP_ROH has 48 intronic SNPs annotated (Table S3), with an average distance among them of 17.7 kb, which is consistent with the average marker distance for this SNP chip spacing of approximately 19 kb (https://genomics.neogen.com/en/ ggp-hd150k-dairy). For this gene, classified as an integral component of membranes (GO:0016021), no association study, at the best of our knowledge, is available. Another example is represented by the KHDRBS2 gene on BTA23 (described above in the text), in which 33 homozygous intronic SNPs are annotated for all the three Aosta breeds at an average distance of 20.5 kb. Lastly, CTNNA3 (n = 32 intronic SNPs spacing 20 kb), both in CAS and VBP, and ADGRL3 (n = 27 intronic SNPs spacing 20.2 kb), CCSER1 (n = 26 intronic SNPs spacing 17.9 kb) and CACNA2D1 (n = 25 intronic SNPs plus a synonymous variant spacing 18.8 kb) in CAS are genes mostly annotated with homozygous SNPs. These genes were previously shown to be under positive selection and associated with marbling score in Korean cattle (CTNNA3) [47], with protein yield and percentage (ADGRL3) [48], with aggressiveness during gestation (CCSER1) [49] and with carcass and meat quality traits in the cattle (CACNA2D1) [50]. This is in line with the selection performed in the three Aosta breeds.
In VRP and VBP, a homozygous cluster on BTA6, involving GABA-A receptor subunit genes (GABRA2, GABRA4, GABRB1 and GABRG1), was identified. These genes being the major inhibitory neurotransmitters in the mammalian brain, mediating anxiolytic activity and playing a key role in emotional and behavioral control in humans [51,52]. Even if speculative, we may comment that the fact that this region does not appear to be in autozygosity in CAS may lead to variability in expression regarding the behavior of cows and for the specific selection operated by the farmers on this breed for the "Battailes des Reines". We may assume that, during the non-ferocious match, the reaction to the visual and physical view of the opponent may be mediated by these groups of genes. The behaviors exhibited during the matches is the same that cows show in pastures when unfamiliar individuals from different herds meet at the beginning of the summer season.

Conclusions
This study, performed on the three Aosta breeds (Aosta Red Pied, Aosta Black Pied and Aosta Chestnut) on a large number of individuals and on a medium density SNP chip (150K), is indicating that farmers in their mating decisions still prioritized adaptation to the environment and thus the farming system envisaging summer pasture, a practice that have characterized these breeds for centuries. All three breeds were shown to have a large number of regions in autozygosity, harboring genes that appear to be linked to efficiency and functional traits, i.e., characteristics for adapting to the environment. As such, with respect to other populations, e.g., Holstein, where breeding plans have strongly acted on genomic regions linked to milk production, the genomic status of these populations shows that, for decades, breeders pursued among the selection objectives the adaptation to the environment, making these breeds resilient and efficient producers in a changing environment. The low extent of recent autozygosity found in this study is a tangible measure of the success of the reproductive scheme to implement selection in the Aosta breeds, as pursued in the last 30 years, aimed also at the maintenance of the genetic variability in these three breeds.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-2615/10/12/2385/s1, Figure S1: Manhattan plot of marker-based F ST values obtained in each breed comparison, Figure S2: Frequencies of ROH per chromosomes and chromosomes mean length (Mb) calculated for the three breeds, Figure S3: Annotation of SNPs defining the three breed TOP_ROH, Figure S4: F ROH along all chromosomes of the three breeds, Table S1: F ST by marker values for CAS vs. VRP (sheet1) and for VBP vs. VRP (sheet 2) comparisons, Table S2: Annotation of genes mapped in CAS TOP_ROH (sheet 1), in VBP TOP_ROH (sheet 2), and VRP TOP_ROH (sheet 3) according to DAVID online Database, Table S3: SNPs annotation in the TOP_ROH regions for the CAS, VBP and VRP breeds.