Genotyping

Genotypes were provided by the National Friesian Italian Breeds Association (ANAFI) for a total of 44 Italian Holstein-Friesian cattle aged 12 to 15 mo all coming from a unique herd, the University experimental station. All individuals had a marker call rate greater than 98 %. Animals used in this study were part of an experimental population raised at University of Milan’s experimental farm, Azienda Agraria Didattico Sperimentale Angelo Menozzi. Animals were genotyped for 30,125 SNP using the GeneSeek GGP Bovine LD v4 array. After excluding SNPs that were not assigned to a bovine chromosome (BTA) or that were assigned to BTA X or mitochondrial DNA, 27,365 SNP remained.

Principal component analysis (PCA)

Additional information available for each individual included its pedigree, Genomic Productivity, Functionality, Type (GPFT) genomic selection index, and batch of genotyping. The PCA was conducted within SNP with the Variation Suite Golden Helix v8.4 software (SVS) (Golden Helix Inc., Bozeman, MT). The procedure used in SVS considered 20 principal components calculating for each of them a relative eigenvalue; the first principal component (eigenvalue=1.188) was plotted against the second (eigenvalue=0.983).

Inbreeding coefficients

The inbreeding estimate of F _IS (or ƒ) was calculated using SVS, defined as 1-(HET_OBS/HET_EXP). This coefficient is equivalent to Wright’s within-subpopulation fixation index with values in the range of -1 to +1. The whole genome inbreeding coefficients (F _ROH) were calculated according to⁹:

Where: L _ROH: is the total length of all ROH in the genome of an individual while.

L_AUTO: refers to the specified length of the autosomal genome covered by SNP (2,505,649,802 bp in the current study).

Pedigree-based coefficients of inbreeding (F _PED) were derived from Pedigree Viewer¹⁰. The pedigree included records for 2,760 individuals up to 19 generations as calculated by Pedigree Viewer that orders genealogical information in discrete generations. F _PED and F _ROH were compared using linear regression and Pearson’s correlation coefficient to test their similarity and to highlight differences in the two methods of estimation in a group of animals of the same herd.

Runs of homozygosity

Runs of homozygosity were estimated for each individual using SVS v8.4. This program does not rely on sliding windows to identify ROH, but instead the algorithm works continuously across an entire chromosome by examining every possible run for a match with the specified criteria¹¹. The following criteria to define ROH were used: 1) Two missing SNP were allowed; 2) One heterozygous genotype was allowed; 3) The minimum number of SNP that constituted the ROH was set to 30; 4) A minimum length of 1 Mb; 5) A maximum gap between consecutive SNP of 1 Mb. Runs of homozygosity were classified into five classes (<2, 2 to 4, 4 to 8, 8 to 16, and >16 Mb) using the same nomenclature as reported by other authors^4,5. The incidence of common runs per SNP was plotted using the Genome Browse tool of Golden Helix, and subsequently aligned to the BTAU 5.0.1 bovine assembly to identify genes consistent with SNP in the top 1% and 5% of runs. The online STRING database classified the network amongst these genes¹², while gene ontology (GO) was performed through the Protein Analysis Through Evolutionary Relationships (PANTHER) classification system (Release 13.1), (http://www.pantherdb.org/pathway/).

Principal component analysis (PCA)

Individuals were grouped based on their batch of genotyping (n= 6), sire (n= 18), sire country of origin (n= 6), and class of GPFT (n= 5). The different batches of genotyping did not cluster whatsoever, indicating that the time of genotyping did not affect individual genotypes (data not shown). On the contrary, there was a strong clustering when individuals are identified by their sire and their sire’s country of origin (Figure 1). Sires 1 and 5 are especially distanced from the others, while those originating from Canada, Italy, and the United States tend to group as well, albeit separately by origin. Individuals were subsequently classified by their GPFT value in increments of 500, and no clustering was observed based upon those assigned categories.

Figure 1 Scatter plot of principal component analysis (PCA) values based on individual SNP genotypes. A) PCA values grouped by sire. B) PCA values grouped by sire origin. C) PCA values grouped by Productivity, Functionality, Type (GPFT) genomic selection index classes 

Inbreeding coefficients

Wright’s inbreeding coefficient estimate of F _IS was calculated for each individual. The majority of animals (n= 31) displayed negative F _IS, while a smaller proportion (n= 13) possessed positive coefficients (Figure 2). The highest and lowest F _IS values were 0.057 and -0.077, respectively. The average F _IS existed at -0.018 shoving an increase in heterozygosity due to an outbreeding mating scheme used by the farmer. The majority of individuals possessed F _PED coefficients between 0.050 and 0.070 (n= 31); values ranged from 0.044 to 0.116. The whole genome inbreeding coefficient estimate of F _ROH was higher overall than the inbreeding coefficients calculated from pedigree. The average F _PED was 0.064, while the mean F _ROH >1 Mb equaled 0.083.

Figure 2 Regression of individual inbreeding coefficient. A) Regression of inbreeding calculated by pedigree (F _PED) on inbreeding calculated on run of homozigosity (F _ROH). R² = 0.205, red line indicates regression line (F _ROH = 0.042 + 0.621 * F _PED ). B) Regression of F _PED on F _IS. R² = 0.214, red line indicates regression line (F _is= -0.077 + 0.913*F _PED ) 

Runs of homozygosity

With the ROH criteria set at 1 Mb and minimum of 30 SNP, the SVS software identified a total of 1,950 runs within the population, ranging from 1 to 22 Mb in length (Figure 3). All 44 animals displayed ROH, with the average number of ROH per animal being 44. The largest number of ROH for an individual was 64, while the least number of ROH for an individual existed at 25. The average ROH length was 4.66 Mb, and the greatest amount of ROH existed in the 4-8 Mb range (Table 1). The number of ROH per chromosome was greatest for BTA 10 (122 runs), while lowest for BTA 27 (24 runs). The longest ROH also existed on BTA 10 at over 22 Mb in length, contrary to other results identifying the longest ROH on BTA 8^13-15. The average number of SNPs falling into a ROH was consistent among ROH length category, ranging from 42 (<2 Mb) to 153 (>16 Mb) SNP. ROH length appeared to be relatively proportional to chromosome size, with longer runs appearing on longer chromosomes.

Figure 3 Number of Run of Homozigosity (ROH) according to their class of length 

The genomic regions most commonly associated with ROH were identified by selecting the top 1% and 5% of SNPs most frequently observed (Figure 4, Table 2). The incidence of common runs per SNP indicated that the genomic distribution of ROH was non-uniform across chromosomes. Quantitative trait loci (QTL) corresponding with regions of homozygosity housing the top 5% of SNPs were identified using the Animal Quantitative Trait Loci (QTL) Database (AnimalQTLdb) (Table 2)¹⁶.

Figure 4: Incidence of SNPs in ROH identified by SVS. The red and blue lines indicate the adopted threshold for the top 1% and 5% of observations, respectively 

Principal component analyses

The PCA is a useful tool that allows farmers to identify genetically different individuals within their herd based on genotypes. The graphical representation of individuals in fact represent an immediate easy to interpret tool that does not require any specific skill by farmers except the concept that closed points are more similar and distant ones are different. Well-differentiated females can be identified in this herd, mostly due to sire, although a group appearing genetically similar exists belonging to a variety of sires. Sires 1 and 5 are responsible for the out-groups we see, with their derivations being Canada and Italy, respectively. Therefore, as proposed by the farmer involved in this study, this visual is beneficial to determine the a-posteriori analysis of his mating strategies and whether the sires he had chosen are helping to his attempts to limit the reduction in genetic variability of the herd. The PCA is also beneficial in sire selection, as the farmer can ascertain if sires deriving from specific countries have an impact on the genetic distribution in the population. Likewise, it can be visualized how genomic values, in this case GPFT, coincide with different individuals. GPFT of similar value did not cluster, indicating that both sire and sire origin did not impact GPFT in the females of this herd. Sire 1 contributed to individuals with GPFT values ranging from 1,500 to 3,000 and sire 5 to animals with values from 1,000 to 2,500, for example.

Inbreeding coefficients

This result here obtained are concordant with another study performed in Italian Holstein, where the average F _PED was 0.044 and the F _ROH >1 Mb was 0.116⁵. Across all 44 animals, the correlation between F _PED and F _ROH was significantly different from 0 with r equal to 0.453 (P<0.01). Intense and accurate artificial selection practiced over many years has resulted in high rates of genetic gain; however, the high rates of gain have been accompanied by large increases in inbreeding¹⁷. Additionally, as population sizes decrease, the probability of mating among relatives increases, especially in small herds or local breeds¹³. Producers have largely depended on pedigrees to estimate the coefficient of inbreeding within the herd, but this number can only increase with each generation. Pedigree relatedness gives an expected, not actual, proportion of genomic identity by descent among individuals and it is anticipated that genotype-based estimates provide greater accuracy on relatedness¹⁸. In practice, pedigree information is difficult to obtain, potentially unreliable, and rarely assessed for inbreeding arising from common ancestors¹⁹. When the pedigree-based coefficients of this herd was examined, it was possible that individuals to be up to a tenth inbred. But, in the absence of pedigree data, the extent of a genome under ROH may be used to infer aspects of recent population history, even from relatively few samples^3,9. McQuillan et al.⁹ revealed that F _ROH correlates strongly with the inbreeding coefficient estimated from pedigree (r= 0.86), and it has been found that F _ROH values are preferable to F _PED, as they are thought to be a more realistic reflection of inbreeding level^15,19. However, it is important to note that the setting of the parameters used to derive ROH is crucial to account for the effects of SNP density correctly²⁰. In this dataset, it is visible a linear relationship between F _ROH and F _PED. The Pearson’s correlation coefficient provides a positive correlation at r= 0.453. The correlation in this work is in concordance with other studies using cattle that were genotyped at similar density^15,20, yet lower than others that were genotyped at medium to high density^3,5. It is possible then that the correlation coefficient may be affected partially by the small number of animals in this study as such as for the use of a low-density SNP chip. Additionally, the F _ROH was larger than the mean F _PED in 35 animals, implying that pedigree based inbreeding coefficients could be underestimated.

The overall low F _IS values are indicative of a low level of inbreeding in the population and of a relatively high number of individuals in a heterozygous state. This may be related to the result of the specific mating strategy implemented by the University herd manager co-author of this work, as he stated his goal was to increase genetic variability and preserve low inbreeding. Therefore, this genomic information is providing unique feedback as to the successful efforts in maintaining inbreeding within this population. Dairy breeds especially have been under more intensive selection and may be related to the recent increase of consanguineous mating resulting from small number of high genetic merit sires used for artificial insemination⁵. Because of this, it is important to understand the inbreeding occurring at a genomic level, as well as to determine whether or not the genetic variability within one’s herd is being well maintained.

Runs of homozygosity and associated genes

The most commonly observed QTL within these regions reported in Table 2 included milk protein percentage, milk fat, calving ease, PTA type, and stature. After aligning the SNP to the BTAU 5.0.1 reference, 260 annotated genes were identified to be consistent with the top 5% of SNP, while 37 genes remained in the top 1%. Genes in the top 1% resided on chromosomes 5 and 10; no annotated genes were identified with the SNP on BTA 16. STRING subsequently identified networks amongst these genes, with 142 genes identified to be involved in some kind of network (Figure 5). These genes also underwent gene ontology in PANTHER. Genes within the top 5% corresponded with 12 biological processes (Figure 6) and 70 pathways (data not shown). The pathway that included the greatest number of genes (n= 12) was the gonadotropin-releasing hormone receptor pathway, while other pathways sharing genes were the 5HT2 type receptor mediated signaling pathway, endothelin signaling pathway, inflammation mediated by chemokine and cytokine signaling pathway, integrin signaling pathway, and Wingless-related integration site (Wnt) signaling pathway. The respective biological process most commonly shared among genes (n= 116) was “cellular processes,” specifically cell communication, while another large portion (n= 79) has a role in metabolic processes (Figure 6).

Figure 5 Network of genes corresponding with SNP in the top 5% of incidence of common runs obtained by STRING online database 

Figure 6 Biological processes of genes corresponding with SNP in the top 5% of incidence of common runs 

The practices of intense selection of sires, artificial insemination, and embryo transfer have been featured heavily in some breeds, reducing effective population sizes, genetic diversity and affecting levels of homozygosity³. The deleterious effects associated with boosted homozygosity arising from inbreeding are predisposed to reduce the genetic gains, implicating in a clear loss of genetic variability¹⁵. ROH give insight to this variability, as they are continuous homozygous segments that are common in individuals and populations. The ability of these segments to give a glimpse into a population’s genetic events makes them a useful tool in providing information about the evolution of the population over time, therefore enabling producers to maintain diversity and fitness within their livestock breed²¹. Furthermore, ROH provide useful information about the genetic relatedness among individuals, helping to minimize the inbreeding rate and also helping to expose deleterious variants in the genome²¹. Long ROH arise as a result of recent inbreeding, while shorter ROH can indicate more distant ancestral effects such as breed founder effects³. In point of fact, the presence of segments longer than 10 Mb is traceable to inbreeding from recent common ancestors that occurred only five generations ago²², and 66 % of the animals comprised in this study presented at least one ROH extending over 10 Mb. This signifies that the majority of individuals in this population derived from recent common ancestors and is the product of recent inbreeding, consistent with what we see from the inbreeding coefficients.

ROH were also identified as a means to provide insight to the conserved genomic regions within the herd. ROH less than 5 Mb were recently characterized as being short²³. These shorter ROH can be related to a more distant ancestral positive selection effect due to recombination events from repeated meiosis breaking long chromosomes into segments and therefore reducing their size^13,24. Lower SNP density, such as that used in this study, tends to inflate ROH³. Given that we see primarily short ROH in this herd (average ROH=4.66 Mb), even with the potential inflation, it is possibly to say with confidently that the level of inbreeding within this population has been well maintained, is low, and is recent based upon the presence of long ROH in the majority of individuals yet the overall short ROH throughout the herd. This result identifies with a separate study that established that the Italian Holstein show higher number of ROH segments related to ancient consanguinity¹³.

Genes within the extended homozygous regions corresponding to SNPs in the top 5% of common runs were found to be largely involved in survival processes. This indicates that basic biological processes, such as cell communication and metabolism, have been maintained in selection efforts. QTL within these regions, however, show that these genomic regions coincide with beneficial production facets such as milk protein percentage, milk fat, and calving ease, although these may be specific to the Holstein-Friesian breed due to stringent selection efforts over time. Additionally, genes in the top 1% (BTA 5, 10, and 16) also correlated with survival and developmental processes. Such processes included cell proliferation (GNG2, ADAM10, ALDH1A2), metabolic function (SCG5, PTGDR, RORA, MYO1E, LIPC, ALDH1A2, MYO5A, GALK2), organ development (ANKS1B, MYO1E, CCNB2, ALDH1A2, PRTG, MAPK6), and immunity (NEDD4). No genes were identified on BTA 16 due to poor annotation, but several genes on BTA 5 and 10 have been well characterized in different phenotypes in cattle. PARP1 binding protein (PARPBP) on BTA 5 has been associated with fat percentage in the third lactation stage in Jersey cattle, while ankyrin repeat and sterile alpha motif domain containing 1B (ANKS1B), also on BTA 5, has been associated with bovine spongiform encephalopathy (BSE) incidence in cattle families^25,26. Meanwhile, genes on BTA 10 were largely found to be involved in fat metabolism and feed intake. Secretogranin V (SCG5) and G protein subunit gamma 2 (GNG2), for example, have been linked to the regulation of hormone metabolic processes involved in feed intake in Holstein^27,28, while transcription factor 12 (TCF12) was found to be associated with lipid and organoleptic traits in European Bos taurus breeds²⁹. GNG2 was also found to be in a region associated with susceptibility to Mycobacterium avium ssp. Paratuberculosis (Map) infection³⁰. Lastly, FERM domain containing 6 (FRMD6) was identified to be in a region of positive selection in the N’Dama breed³¹. The regions on BTA 10 and 16 harboring the top 1% of SNP were also concordant with other studies. Studies involving Holstein were especially similar, with positive selection occurring on BTA 10 at 50-60 Mb and a putative QTL identified at 34.7-56.9 on BTA 10 for somatic cell count and non-return rate at 90 days (paternal effect)^20,32,33. Another study identified a high proportion of ROH on BTA 16, while a separate found selection at regions 24.7-26.7 and 26.5-28.5 on BTA 16^3,20.