Genomic approach to manage genetic variability in dairy farms

Abstract In this study we investigated the genetic variability, the inbreeding and allele frequencies of monogenic traits in seven herds of Holstein breed and provided insight to farmers on the value of genomic management of reproduction in their herds. A total of 3,953 Holstein cows were sampled and genotyped with the Neogen GGP Bovine 100K SNP chip within the activities of the Regione Lombardia funded GO-PEI project ‘GENOmic tool for the management of reproduction in dairy cattle and for the control of inbreeding – GENORIP’. Principal component analysis was applied for analysing the genetic variability within and among farms using the SVS software of Golden Helix. Run of Homozygosity (ROH) and the genomic inbreeding were obtained with the detectRUNS package of the R software. Genotype frequencies for mendelian disease, fertility and production traits were also obtained. A total of 458,267 ROH were identified and ROH were distributed on all autosomes with an average length of 2,703,811 bp covering 12.7% of the genome. Several genomic regions appear under selection, while a specific region on BTA4 was identified in one herd, harbouring genes mainly related to the specific selection strategy of the farmer. The FROH values obtained considering ROH greater than 16 Mb, varied from 0.004 to 0.325, with the highest FROH average value of 0.136. Among mendelian heritable diseases, the Haplotype Cholesterol Deficiency was the one with the largest proportion of carrier animals, i.e. 5.6%. A herd-tailored process to assist farmers in genomic management of reproduction was released. The ROH distribution within herd, together with the genotype frequencies for disease, fertility and production mendelian traits, suggest that similar directional selection is occurring across herds. This study released to each farmer the genomic make-up of their herd used jointly with the gEBV estimated by their national breeders’ association (ANAFIBJ) for herd reproductive management. HIGHLIGHTS All females of 7 herds have been genotyped with the GGP 100K SNP chip. Genomic information on all females can be used by farmers in the process to manage reproduction, selection and genetic herd variability. The availability of genomic information on the whole herd allowed to release to each farmer the genomic make-up of their herd. The ROH distribution together with the genotype frequencies for disease, fertility and production mendelian traits, made it possible to identify genomic regions under selection according to farmer strategies.


Introduction
The Holstein breed is widely recognised as the most productive dairy cattle in the world.The Holstein is the most wide-spread dairy cattle breed in Italy with 9,552 farms, 1,130,734 lactating cows and an average production of 10,710 kg of milk.In particular, in the Lombardy region located in the northern flat, which is suited to intensive dairy farming, there are 2,759 farms with a total of 566,583 animals (ANAFIBJ -National statistics 2022).
The selection index of the Holstein breed has evolved through the years, with an initial emphasis on increasing milk yield per cow, followed by a shift towards milk components and functional and health traits (Egger-Danner et al. 2015).However, as in all selection programs, the intense selection practiced over the years may led to a loss of genetic variability and to an increase in inbreeding in the population.The need to control the increase in inbreeding even in large populations under selection has been discussed for a long time (Mcdaniel 2001;Weigel 2001).In the second half of the last century, there was a motivation to introduce new molecular tools to integrate traditional phenotypic selection programs (Henderson 1975) with the use of information of loci and QTL regions that contain genes capable of influencing economically important traits in animal production (Georges et al. 1995;Andersson 2001).
In the genomics era occurred in the last decade, the paradigm of animal breeding has changed significantly.Current genotyping techniques make it possible to determine the genotype of an animal at hundreds of thousands of markers known to be associated with phenotypic variability at very low cost and use this information to select animals even without any performance available.This is the principle of genomic selection proposed by Meuwissen et al. (2001) based on the use of Single Nucleotide Polymorphisms (SNP) as markers.
The SNP are biallelic genomic markers very frequent in the genome of any individual (approximately one per 100 base pairs) (Collins et al. 1997).The SNP genotyping technology is nowadays a routine process in cattle breeding, both for males and females, producing a large number of genotyping information, allowing as such to implement efficient selection programs also for traits with low heritability (Boichard et al. 2015), and to develop comprehensive mating plans that make use of all the genomic information available to the breeder female herd.
One of the most important elements necessary to perform selection in livestock populations is the existence of genetic variability.Indeed, the occurred selection of superior animals over time has resulted in a loss of genetic diversity.This may cause a reduced response to selection and an increase of the frequency of homozygous loci (Dickerson and Hazel 1944).
Before the advent of genomics, the study of inbreeding was based on pedigree information which, however, has limitations: 1) the value of inbreeding depends on the quality and completeness of the pedigree data (Oliehoek and Bijma 2009); 2) it does not consider the genetic variability between full siblings due to recombination during meiosis leading to an underestimation and/or overestimation of inbreeding (Hill and Weir 2011); 3) the comparison between genomic and pedigree information, showed that the frequency of misidentified bulls can be as high as 13.9% (Wiggans et al. 2012), affecting as such the inbreeding values based on this information.With the genomic selection and the development of high-density SNP arrays, it has become possible to obtain more accurate estimates of genome-wide inbreeding and relatedness (Engelsma et al. 2012).Genomic inbreeding can be calculated from a genome-wide relationship matrix (GRM) between individuals (Hayes et al. 2009), or as ratio between the length of the genome where homozygous markers form Run Of Homozygosity (ROH) (McQuillan et al. 2008) and the total genome length analysed (FROH).The length of the ROH provides also information on whether a ROH segment is the result of recent (long ROH) versus more distant (short ROH) autozygosity events (Pemberton et al. 2012).Additionally, identification of genes annotated in the ROH can provide insights on the selection occurring in the population (MacLeod et al. 2009;Purfield et al. 2012).
Due to intense selection, in dairy populations the inbreeding rates (Charlesworth and Willis 2009) and the frequency of deleterious alleles (Ouborg et al. 2010) have increased significantly over the years with a possible effect on productivity due to inbreeding depression (Falconer and Mackay 1983;Keller and Waller 2002).In Friesian cattle deleterious effects due to this phenomenon have been described on productive and functional traits (Martikainen et al. 2018;Doekes et al. 2019).
For these reasons, fostering genetic variability in the herd and controlling inbreeding is considered a priority in dairy farms and in management of cattle populations under selection.A farm-driven project funded by the Regione Lombardia in the EC EIP-AGRI Rural Development Program 2014-2020 framework is bringing genomics into the management of female replacements through the genotyping of all the animals in the herds.The project is named GENORIP: 'GENOmic' tool for the management of reproduction in dairy cattle and for the control of inbreeding.The project aimed to release a process to integrate the genomic management of herd reproduction and to manage inbreeding and genetic variability using dense SNP genotyping data.
This study is part of the activities of the GENORIP and aimed to investigate the genetic variability, the genomic inbreeding and the allele frequencies of hereditary monogenic traits in the females of seven Italian Holstein large dairy cattle herds.

Ethics statement
The study was approved by the OPBA of the University of Milan (Protocol number 160_2019), in accordance with the Directive 2010/63/EU of the European Parliament and of the Council of 22 September 2010, updating Directive 86/609/EEC of 1986 on the protection of animals used for scientific purposes.
Animal sampling, DNA extraction, quality control and genotyping A total of 3,953 Holstein cows were sampled from seven different herds located in the Lombardy region of Italy.These herds were chosen for their different sizes and management practices, in order to provide some examples of dairy cattle farms in the Lombardy region.The farms differ, in fact, in their structures and available technologies (such as milking parlours vs. milking robots) and size (120 milking heads vs 600 milking cows).
Animals were sampled using ear Tissue Sampling Units (TSU) for adult individuals and bioptic ear tags for newborn calves.The collected samples were then classified in a project structured database and stored at the University of Milan tissue repository Animal Bio-Arkive (Longeri et al. 2021).
The DNA was extracted from ear tissue using the Quick-DNA TM Miniprep Kit of Zymo Research according to the manufacturer's protocol (Zymo Research Corporation).Cows were genotyped with NEOGEN's GGP Bovine 100K, consisting of approximately 100,000 SNPs, with an average SNP spacing of about 29 kb.
All samples had a call rate value >95%.Only SNPs located on the 29 autosomes annotated according to the ARS-UCD1.2bovine genome assembly (n. 89,762) were considered in this study to perform all analyses.In order to avoid bias, we excluded SNPs detecting the same mutation: they were more than 600 SNPs on the autosomes.

Analysis of population structure
Principal component analysis (PCA) was used to determine genetic diversity within and among herds and has been performed using SNP & Variation Suite (SVS) v8.9 (Golden Helix Inc., Bozeman, MT, USA).The 2-D graphical visualisation of PCA results was obtained using the 'ggplot2' R library (Wickham 2016).

Genotype frequencies for health, phenotypic and productive traits
The GGP Bovine 100K chip releases a large number of SNPs genotypes associated to mendelian hereditary traits, such as genetic disorders and mutations related to phenotypic and productive traits or haplotypes linked to fertility traits.The allele and genotypic frequencies for these loci were estimated using an inhouse R script.

Runs of homozygosity detection
Runs of Homozygosity (ROH) were obtained for each individual using the consecutive run method of the 'detectRUNS' library of R software (Biscarini et al. 2019).The parameters used were: (i) minimum number of 30 SNPs/ROH; (ii) a minimum length of 1 Mb for the identified ROH, to avoid the detection of short and common ROH across the genome due to Linkage Disequilibrium; (iii) a maximum distance of 1 Mb between consecutive SNPs to eliminate the bias in detection due to the density occurrence of SNPs; (iv) no missing SNPs as well as no heterozygous genotypes presence in ROH definition.
The ROH distribution per herd was evaluated separately using five classes of ROH length (<2 Mb, [2][3][4][4][5][6][7][8][8][9][10][11][12][13][14][15][16].Descriptive statistics relative to the total number of ROH, the ROH average number per individuals and, the average length of ROH were calculated. The 'detectRUNS' library was also used to obtain: i) the graphical representations (Manhattan plots) for the percentage of occurrence of SNPs in ROH, estimated by counting the number of times that each SNP falls inside a ROH over the total number of individuals; ii) the ROH_islands, identified as peaks in Manhattan plot where SNPs are inside a ROH in more than 50% (chosen threshold) of the cows as discussed and suggested by Schiavo et al. (2021).

Gene annotation of ROH_islands and functional analyses
All ROH_islands were annotated with the genes downloaded from the NCBI online Database (NCBI Annotation Release: 106).Only genes with an official gene name were considered.Database for Annotation, Visualisation, and Integrated Discovery (DAVID) v6.8 (DAVID online Database) was used to perform a gene ontology (GO) functional annotation and KEGG pathway analyses.
Additionally, the CattleQTLdb database (AnimalQTLdb) was used to identifyusing the 'Search by associated gene' optionthe QTL overlapping the ROH_islands.

Inbreeding coefficient
The genomic inbreeding coefficients (F ROH ) based on ROH were calculated for each cow as: where L ROH is the total length of ROH proper of each individual genome, and L AUTO is the total genome length covered by the used SNP dataset (2,487,916,500 bp in this study).F ROH were calculated for each of the five classes of ROH length (L ROH ) previously defined.

Results and discussion
Table 1 reports the number of individuals sampled and genotyped in each farm, and the average gEBV for several traits released by ANAFIBJ for all genotyped females.All paternity and maternity consistency was verified based on the genome data by ANAFIBJ to solve inconsistencies due to incorrect genealogy, i.e. errors in sire or in maternal grandsire registration.The proportion of these inconsistences varied from 8% to 45%.

Population structure
A first sight of the genetic variability of the 7 herds is provided by the graphical representation of PCA shown in Figure 1.As visualised in Figure 1, within herd PCAs, cows cluster clearly in separate groups in Herd_3 (PCA1 ¼ 15.88%, PCA2 ¼ 13.94%), Herd_4 (PCA1 ¼ 1 5.02%, PCA2 ¼ 12.63%), Herd_6 (PCA1 ¼ 19.67%, PCA2 ¼ 14.44%), and to some extent also in Herd_2 (PCA1 ¼ 16.81%, PCA2 ¼ 16.16%) and Herd_5 (PCA1 ¼ 13.11%, PCA2 ¼ 12.40%), while it appears to exist more homogeneity among cows for Herd_7 (PCA1 ¼ 14.79%, PCA2 ¼ 13.08%) and Herd_1 (PCA1 ¼ 15.24%, PCA2 ¼ 12.22%).The cow clustering is expected to reflect the choices made by farmers in terms of use of sires for reproduction.More specifically it is likely that the variability shown by PCAs depends on the sire origin as system/country of selection scheme (e.g.USA, CAN, NLD, DEU, FRA, ITA) is mediated by the AI centres selling the semen to farmers.Discussion with farmers (partners of the project) on this topic disclosed a different approach in sire choice: some farmers rely on the technical advice (and semen) from a unique AI centre, some others select sires personally across all available on the market, also taking advantage of the information deriving from the mating plans offered by the farmer association.Only one, Herd_4 is selecting sires based on a herd genomic selection on females already applied for some years.All gEBVs here reported (Table 1) are based on the breeders' reproductive choices that were made without taking into account the genomic information of the females in the herd, with the exception of Herd_4.
To provide a rational for the cows' clustering in PCAs, we investigated the variability of sires used in farms as number of daughters from same sire/maternal grandsire (i.e. a bull being both sire and maternal grandsire in the same herd), within herd and among herds (Table 2).
In relation to the herd size, Herd_3 and Herd_6 use the lowest number of bulls: each sire has, in fact, an average 10.5 and 8.1 daughters, respectively.Herd_7 is the one with the largest number of sires (3.2 daughters per sire, on average) accounting for the herd size.The large size of Herd_3 somehow affects the possibility to use a large number of sires in the breeding plans, maintaining a high genetic level of the group of males: using a large number of sires to decrease the number of daughters per male, would in fact diminish the average genetic value of the reproducers.To avoid decreasing too much the genetic level of the bulls Herd_3 accept to have larger groups of daughters per sire, if compared to other herds.Nevertheless, the 5 clusters visible for Herd_3 in Figure 1 indicate that the farmer in selecting sires, in addition to the selection goal, was also paying attention to the genetic variability: the EBVs of the cows for PFT (selection index for Productivity, Functionality and Type) are in fact comparable to those of other herds, with a higher average value for milk gEBV (1,052 kg).On the contrary, the large spectrum of sires used (173 sires and 186 maternal grand sires, the largest in all herds), was widening the genetic variability, but only at a very low extent also reducing the average genetic values of used sires.It is noteworthy that Herd_7 is the one with the lowest proportion of bulls being at the same time sires and maternal grand sires of females in the herd.This indicates a fast change in sires used in reproduction and a great attention to reduce genetic inbreeding.Herd_4 has even fewer sires being also maternal grand sires, only 38%, an indication consistent with the applied selection plan of the farmer, based on genomic selection of females for several years.Herd_4 farmer is using young sires of the most recent generations as much as possible; the selection is addressed to prioritise improvement of functional traits, with particular emphasis on fertility and longevity.The mating plan based on genomic information used by Herd_4 appears very successful, both for the very good gEBVs of females, greater than all other herds for functional and production traits (Somatic Cell Count, Udder Health, Longevity and Fertility), the selection indexes (PFT, IES$, ICS-PR), and for the maintenance of genetic variability.Herd_2 is under a very successful introduction of precision farming (Automatic Milking Systems), never used genomic selection at farm female level and is systematically applying the same breeding goal in the herd for the past several years.The impact of GENORIP on this farm was positive as it allowed the farmer a fast step forward in matching the technology available in the herd with the genomic information to manage female reproduction in the herd.
When we compared the number of common bulls across herds, this can contribute to the explanation why Herd_6 is clustering separately from others: the number of sires in Herd_6, when compared to others, is a maximum 15 in common with maximum 2 other farms, while others have in common 21 to 59 sires up to 4 herds.The mating plan in Herd_6 is in fact fully relying on the technical advice of a unique AI centre with all bulls deriving from its selection program, while other herds acquire bulls from various AI centres.

Genotype frequencies for productive traits
Regarding genes linked to milk production and quality (Table 3), the AA SNP mutation of a-S1-Casein is close to 100% in all herds, while the GG SNP mutation wasn't identified: as a key to compare to other studies the SNP allele A mutation correspond to the B variant and the G SNP mutation correspond to the C variant as usually reported (Sanchez et al. 2020).In Holstein the effect of B-and C-variants for a-S1-Casein were identified (Poulsen et al. 2013), with the B variant linked to increase in milk yield and the positive effect of C variant on curd coagulation time and curd firmness rate (Bovenhuis et al. 1992).The frequency of the B variant indicates that the genetic potential of our herds is for milk production as it is closely to be fixed in the population.
Also for the b-Casein locus the most frequent genotype is AA in all herds (ranging from about 84% to 95%), whilst the BB genotype was found with a low genotype frequency only in Herd_1 (1.07%) and in Herd_7 (1.10%).b-Caseins show numerous genetic variants that result in different quality characteristics in milk.The most common variants are A1 and A2.A1A1 is the less frequent genotype variant of b-Casein in all herds; while the A1A2 and A2A2 vary according with herds, ranging from 37.4% to 51.7% for A1A2 and from 28.8% to 58.0% for A2A2.The molecular difference between the two proteins is related to a mutation resulting in an amino acid change (proline vs histidine) at position 67 of b-Casein (Ginger and Grigor 1999).The amino acid change was associated with a different gastric digestion of caseins.Indeed, during the enzymatic digestion of A1 casein, an opioid peptide (BCM-7) is released, which is not released in the digestion of A2 variants (Brooke-Taylor et al. 2017).In recent years, there has been an increased focus on the b-Casein A2 allele as some studies have suggested that the b-Casein A2 allele is better tolerated by the human population (He et al. 2017).To date, however, no relationship has been found between the consumption of cow's milk with the A1 allele for b-Casein and disease incidence.In addition, the A1 variant improves rennet coagulation properties compared to the A2 variant (Dinc et al. 2013;Ketto et al. 2017).Interest in marketing dairy products, with improved health impact, has opened the market to  milk selected for its b-Casein A2 content only (Mendes et al. 2019).As Table 2 shows, it is evident that two herds, in particular Herd_5 and Herd_6, have a proportion of A2A2 genotypes >50%.These breeders, in fact, select for this genotype while the others generally have higher values for the heterozygous A1A2 genotype.The higher proportion of heterozygous genotype variants is in agreement with those reported by other authors for Holstein breed (Massella et al. 2017).
For b-Lactoglobulin, depending on the herd, the AA or AB genotypes are the most frequent genotypes: i.e. the AA genotype frequencies is higher in Herd_5 and Herd_7 (50% and 43.86%, respectively); in all other herds the higher frequency has been found for the AB genotype (ranging from 46.81% in Herd_4 to 50.95% in Herd_2).The b-Lactoglobulin is the major serum protein in cow's milk, accounting for about 50% the total amount of milk proteins and the B variant has been associated with a higher casein content, resulting in a higher cheese yield (van den Berg et al. 1992;Stasio and Mariani 2000).
The frequencies of the six K-Casein genotypes (AA, AB, AE, BB, BE, and EE) had the same pattern in all herds: the most represented genotype was AB, ranging from about 36.97% in Herd_3 to 50.45% in Herd_ 5. Instead, the EE genotype has been registered in a very low number of cows (not one EE was found in Herd_2 and Herd_5) (Table 3).Studies in the literature show that in Holstein both A and B alleles are the most frequent (Prinzenberg et al. 1999;Farrell et al. 2004) and the E allele the least frequent (Caroli et al. 2000).In fact, negative effects on coagulum formation during cheesemaking have been observed in milk produced by individuals carrying the E allele variant (Caroli et al. 2000;Comin et al. 2008).As Table 3 shows, heterozygous are generally the most widespread, while the BB variant, which has intermediate values, positively influences the production of cheeses, such as Parmigiano Reggiano and Grana Padano, by increasing cheese yields, as shown in the study of Mariani et al. (Mariani et al. 1976).
The heterozygous cows (AG) at marker linked to milk YellowFat feature, are still present with very low frequencies in Herd_3 (0.47%), Herd_4 (0.91%) and Herd_7 (0.36%).The AA genotypes causes a characteristic yellow colour of fat in tissues and milk, due to carotenoids depositions in adipose tissue.(Yang et al. 1992).
Finally, for the Holstein cows here analysed, at Lactoferrin locus, we found a similar distribution of the three genotype frequencies across herds, with a higher frequency of AA, mainly in Herd_4 (about 67%).AA genotype was associated with a low milk SCC values (Wojdak-Maksymiec et al. 2006).
Regarding meat traits, we observed high variability at all analysed loci (Table 3), as the Holstein breed is not selected for these traits.

Genotype frequencies for reproductive traits and disease
Table 4 shows the genotype frequencies of genes and haplotypes influencing bovine fertility.Animals bearing mutations affecting reproduction efficiency were found for all haplotypes, with different frequencies in the seven herds.The haplotype HH5 was the one with the largest number of carriers in all herds, with about 10% of females being carriers in Herd_7.Instead, high carrier haplotype frequencies were found in some herds only for specific HH.HH4 carriers are counted only for Herd_3 and Herd_7.The genotype frequencies for COQ9-rs109301586, STAT3 and 5, Leptin_2F, and PKP2_988 markers loci (all mainly involved in embryo development) are similar across all herds, counting a higher proportion of one or both homozygous genotypes, except for STAT3_25402 and STAT5_13319, for which the most frequent genotype was BB (Table 4).
In this study, the highest carrier frequency was found for HCD (haplotype cholesterol deficiency), that represents an economic loss for the farmer.Animals homozygous for this disease show the first clinical signs between 1 and 5 months of age with decreased appetite, weight loss, diarrhoea and subsequent death with a frequency in the German Friesian population of 4.2% (Kipp et al. 2016).In our study only one animal from Herd_3 was affected.In previous study, it was shown that animals carrying HCD had significantly higher protein yields than non-carriers, but it is still unclear how HCD affects cheese yields (Cole et al. 2016).
The presence of a high proportion of carrier haplotypes/SNPs may be linked to the breeders' choice to concentrate on the selection of productive traits, which are negatively correlated with certain reproductive traits, or the use of bulls carrying some haplotype.Some of these haplotypes influence heifer conception rate, cow conception rate, milk, and protein (Cole et al. 2016).The increased consideration of these haplotypes in cattle selection criteria could lead breeders to a gradual improvement in herd fertility, reducing the losses associated with it.
Regarding genetic diseases, this study has shown the presence of carriers with low frequencies (ranging from 0.07 to 4.28, Table 3) for BLAD, GSDV and RP1 in all herds.Higher carrier frequencies were observed for Brachyspina (up to 7.22% in Herd_2).
The proportion of BLAD-carrier animals ranged from 0.07 to 0.66 which, compared with other studies such as in Brazilian Holsteins that found a carrier frequency of 5.7% (Ribeiro et al. 2000), is very low.
Avoiding carrier-to-carrier mating, and thus identifying heterozygous cows, would be a way of managing the reproduction and presence of female carriers in the herd.

Runs of homozygosity detection
A total of 458,267 ROH was identified in all cows of the seven herds.The count of ROH (per herd) reflects the size of herd sampling (correlation ¼ 0.997).At the individual level, the average number of ROH ranged from 7 (Herd_1) to 251 (Herd_3), with a similar total mean ROH length close to 2.6 Mb, except for Herd_5 that had on average longer ROH, close to 3 Mb (Table 5, Figure 2A).Herd_2 showed both the lowest mean number of ROH per individual (106.8) and lowest total genome length (average value) covered by ROH (11.3%) (Table 5).Differences among cows were identified also considering the total length of the genome covered by the ROH (sum of all ROH per animal, Figure S1).
Same selection occurring across herds may have affected same regions of the genome; at genomic level some evidence may be related to the fact that across herds ROH are found in largely overlapping genomic regions among females of different herds.Over the years selection may have affected same regions where genes involved in expression of traits under selection are annotated (Zhang et al. 2015).
The ROH were found for all classes of length (Figure 2B), with shorter regions (<2 Mb), being the most frequent classes of length (about 50%), even if this proportion may be slightly overestimated according to results from (Feren cakovi c et al. 2013), study however based on a 50K SNP chip and not on a 100K SNP chip array.Contrariwise, a small number of ROH    longer than 16 Mb were mapped in all herds (observed frequencies ranging from 0.05% -Herd_5 to 0.30% -Herd_5) with a maximum of two ROH longer than 16 Mb per individual.
Finally, ROHs were found over all chromosomes: there was no evidence of a relationship between the chromosome's length and mean ROH length, as shown in Figure 2C.
The graphs in Figure 3 shows the proportion of SNP in ROH segments across all the autosomes (Manhattan Plots) for all Herds.
ROH_islands (i.e.SNP within ROH with frequency value greater than 50% as herein before defined) were detected for the Herds 2, 4, 5, 6, and 7 and are listed in Table 6.On chr 7 two very close regions were identified for Herd_6.ROH_islands identified on chr 10 and 20 are identified in three and two herds, respectively.These two genomic regions overlap to that identified in other Holstein cows bred in Italy (Mastrangelo et al. 2018).All ROH_islands except one mapping on chr 7, harboured annotated genes (n. 68).
According to the Animal QTL database, the genes lying within the ROH_Island located on chr 10 (region shared by cows of three different herds) are mainly associated with reproduction traits (e.g.fertility traits) and morphology traits (e.g.Udder and Conformation traits).Among the genes annotated within the ROH_ Island reported in Table 6, the ERBB4 and MKRN3 genes affects udder and fertility traits.In detail, ERBB4 was identified as the hub gene of the network that regulates udder growth and development and seems to affect the genes' expression that are involved in the udder involution and that promote mammary gland remodelling (Xuan et al. 2022), whereas MKRN3 controls the initiation of puberty (Abreu et al. 2015) and inhibits the reproductive axis (Abreu et al. 2020).The ROH region on chr 20 is under selection in Herd_ 4 and Herd_5 and includes the PELO gene.This region

Inbreeding coefficient
As reported in Table 7, the F ROH values varied from 0.004 (Herd_1) to 0.325 (Herd_3), with an overall average value of F ROH ranging from 0.113 to 0.136.These values are comparable with the genomic inbreeding calculated in the US Holstein by Lozada-Soto et al. (2022) and in Italian Holstein as reported by Dadousis et al. (2022).As shown in Figure 4, the distribution of the F ROH values calculated per each class of ROH length differed among the herds, ranging from 0.113 (Herd_1) to 0.136 (Herd_4).For the two greater classes of ROH, representing the most recent genomic inbreeding, the average values (per farm) were between 0.012 and 0.023 for ROH of 8-16 Mb and 0.008 and 0.010 for ROH > 16 Mb (Table S2).We want to highlight that the maximum proportion of cows with ROH > 16 Mb was identified in Herd_5 (25.5%).In other herds this proportion is lower than 17% with a minimum value in Herd_ 3 of 5.7%.The inbreeding coefficients here obtained for each class of length were lower (except for the class of length <2 Mb) than calculated in Italian Holstein breed (Ablondi et al. 2022).Considering this overview, we can easily deduct that the inbreeding is taken under control in all farms in the last decades by the farmers, applying breeding strategies aimed to maintain genetic diversity among cows.

Conclusions
In the last decade, genomic selection has been very successful and rapidly adopted in the genetic improvement plans of large dairy cattle populations, such as the Italian Friesian breed.The introduction of genomic selection in selection schemes made the improvement for low heritable traits, as functional or health traits, much more efficient.In addition, the availability of SNP genotypes also on females is making it possible for farmers to customise herd breeding goals by implementing efficient selection, especially for functional and health traits, and develop comprehensive mating plans that exploit all the information available to the breeder.
The use of SNP genotypes on females can be extended to optimise the herd mating plans also to manage herd genetic variability, control inbreeding at a genomic level and for specific selection for mendelian monogenic traits.
In this study, the analysis of genotypes produced by the GENORIP project provided a snapshot of the genetic variability, of the genomic inbreeding, as well as the presence of mendelian genetic variants linked to traits of interest in seven Holstein dairy herds for a total of 3,953 animals.A particularity of the project was to genotype all females present in the herds and the newborn female calves along the three years of the project duration, allowing as such to validate all the genealogical information.In addition, the availability of gEBV for all females makes it possible to evaluate the selection program adopted by farmers based on the sire side.
The knowledge of both milk properties and carriers of unfavourable traits is useful to farmers who process their milk into cheese or regarding payment according to the protein composition of the milk, and to implement mating plan, respectively.

Figure 1 .
Figure 1.Graphical representation of PCA results both for each herd and for all individual together.

Figure 2 .
Figure 2. Graphical representation of ROH statistics per herd: (A) relationship between number and averaged total length (Mb) of ROH proper of each cow; (B) frequencies of ROH for each class of length together with details on the > 16 Mb ROH class of length.

Figure 3 .
Figure 3. Manhattan plot of the proportion of SNPs in identified Runs of Homozygosity (ROH), along all the autosomes, for all analysed Herds.
has been recently found under selection also in US Holstein and Jersey by (Lozada-Soto et al. 2022) who characterised the ROH in several dairy cattle populations.In their study this ROH region resulted wider respect to the one here identified.

Figure 4 .
Figure 4. Descriptive statistics of total F ROH (Table7) and graphical representations (Boxplots) of F ROH calculated in concordance with the five Runs of Homozygosity (ROH) classes of length.

Table 1 .
Sampling and milk destination of each considered herd together with descriptive statistics (minimum, maximum and mean) of gEBV (genomic EBV) values calculated by ANAFBJ.

Table 2 .
Number of sire and maternal grandsires used as reproducers in each herd and shared among them.

Table 3 .
Genotype frequencies for productive traits (i.e.milk protein and meat variant of monogenic related loci) in each herd.

Table 6 .
Runs of Homozygosity (ROH)_Islands detected and shared in at least 50% of cows together with the annotated genes annotated and associated traits.