Results and discussion

The main aim of this study was to analyze estimates of inbreeding derived from ROH in three important Italian local cattle breeds. Moreover, genotypes from Italian Holstein were also included in these analyses in order to compare results among breeds.

We used a definition of ROH as tracts of homozygous genotypes that were >4 Mb in length identified with a minimum number of 40 SNPs. In fact, the density of SNP panel used to generate the data for ROH identification is an important factor that strongly affects autozygosity estimates. Ferenčaković et al. (Reference Ferenčaković, Solkner and Curik2013b) showed that the 50K panel revealed an abundance of small segments and overestimated the numbers of segments 1 to 4 Mb long, suggesting that it is not sensitive enough for the precise determination of small segments. We estimated mean F _{ROH>4 Mb} values allowing one, two and three heterozygous SNPs and paired t-tests were conducted within each cattle breed. In fact, considering that genotyping errors in SNP chip data do occur, it seems reasonable to allow some heterozygous calls, especially for long segments that are more frequent in cattle populations (Ferenčaković et al., Reference Ferenčaković, Solkner and Curik2013b) than in human species (Kirin et al., Reference Kirin, McQuillan, Franklin, Campbell, McKeigue and Wilson2010). The results showed different values depending on whether one, two and three heterozygous genotypes were allowed (Table 1). The differences between F _ROH estimated using one and two heterozygous SNPs were very small in all breeds and did not have important effects on estimates of inbreeding levels, with the highest value of 0.003 units in Italian Holstein and Modicana (Table 1). The highest different values of F _ROH were observed when one and three heterozygous SNPs were compared, with the highest value of 0.007 units for the same above mentioned breeds. Ferenčaković et al. (Reference Ferenčaković, Solkner and Curik2013b) suggested that for long ROH (which can have >5000 to 6000 SNPs) some heterozygous calls must be allowed, especially with high-density chip, but at the same time, the number of allowable heterozygous calls should be limited. In fact, the same authors showed that allowing certain minimum numbers of heterozygous SNPs leads to inaccurate ROH calls, in particular at the ends of ROH. Marras et al. (Reference Marras, Gaspa, Sorbolini, Dimauro, Ajmone-Marsam, Valentini, Williams and Macciotta2014), in a study of ROH using medium-density chip, reported that when heterozygous SNPs were allowed, the number of longer ROH increased dramatically, and preferred not to use them in the ROH. Therefore, considering that in our study medium-density SNP data were used, and that the longest segment was below 2000 SNPs, only one heterozygous SNP was allowed in the ROH in order to avoid underestimation of long ROH.

Table 1 Comparison of inbreeding derived from runs of homozygosity (F _ROH) values obtained by allowing different numbers of heterozygous (het) single nucleotide polymorphisms (SNPs)

^a,b,cDifferent letters indicate statistical significance within the same breed (P<0.001, paired t-test).

We analyzed animals from four Italian cattle breeds with different inbreeding background and selection histories. In particular, Cinisara and Modicana are two ancient Sicilian breeds that are not subject to breeding programs (Mastrangelo et al., Reference Mastrangelo, Saura, Tolone, Salces-Ortiz, Di Gerlando, Bertolini, Fontanesi, Sardina, Serrano and Portolano2014), whereas Reggiana is characterized by limited selection program. For this breed, only few studies have been carried out so far with the aim to identify associations with production traits that might be useful to refine selection and conservation programs (Fontanesi et al., Reference Fontanesi, Scotti, Samorè, Bagnato and Russo2015). Holstein dairy cattle has dominated the milk production industry over decades. 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 increase of inbreeding (Rodríguez-Ramilo et al., Reference Rodríguez-Ramilo, Fernández, Toro, Hernández and Villanueva2015).

A total of 3661 ROH were identified among the four breeds. All individuals of Italian Holstein displayed at least two ROH, whereas in the local breeds there were individuals that did not show ROH >4 Mb. In all breeds, except for Reggiana, the number of ROH per chromosome was greater in BTA1 and BTA2, and tended to decrease with chromosome length. The maximum size of ROH was 112.65 Mb and was found on BTA8 in Cinisara breed. Kim et al. (Reference Kim, Cole, Huson, Wiggans, Van Tassell, Crooker, Liu, Da and Sonstegard2013) showed similar results in Holstein cow with the maximum size of ROH of 87.13 Mb on BTA8. Modicana and Italian Holstein breeds showed the longest ROH on BTA9 (89.61 and 70.11 Mb, respectively), whereas the Reggiana breed on BTA4 (102.18 Mb). Modicana breed showed the highest MN_ROH per individual and the highest value of F _{ROH>4 Mb} (11.03 and 0.055, respectively), whereas Reggiana breed showed the lowest ones (7.15 and 0.035, respectively) (Table 2). L_ROH values indicated low variation among the four breeds showing that this value is not a good descriptor of ROH as reported by other authors (Marras et al., Reference Marras, Gaspa, Sorbolini, Dimauro, Ajmone-Marsam, Valentini, Williams and Macciotta2014). The comparison of ROH is not straightforward since different studies used different criteria in particular for the minimum length of ROH and the minimum number of SNPs involved in ROH. Furthermore, the number of SNPs, density of the SNP chip and selection criteria for SNPs used to determine the genomic inbreeding can have a huge effect on these values (Bjelland et al., Reference Bjelland, Weigel, Vukasinovic and Nkrumah2013). Ferenčaković et al. (Reference Ferenčaković, Hamzić, Gredler, Solberg, Klemetsdal, Curik and Sölkner2013a) found higher number of ROH in four analyzed cattle breeds probably because of the shorter length considered to define the ROH (>1 Mb). Similar results of F _ROH>4 Mb were reported by Ferenčaković et al. (Reference Ferenčaković, Solkner and Curik2013b) using a 50K panel for Pinzgauer (0.037) and Tyrol Grey (0.042) local cattle breeds, and by Marras et al. (Reference Marras, Gaspa, Sorbolini, Dimauro, Ajmone-Marsam, Valentini, Williams and Macciotta2014) in Marchigiana (0.046) beef cattle breed. Differences among breeds existed also for the ROH length. Figure 1 showed the total number of ROH and the total lengths of genome in ROH for each individual of the four breeds. Considerable differences among animals and breeds have been found. The individuals of Italian Holstein breed showed high number of short ROH segments. Similar results were showed for Reggiana breed with some extreme animals with segments covering 400 Mb and more of genome, and with a number of ROH per individual >25. The Sicilian breeds showed analogous results between them with the total length of ROH characterized by the presence of large segments. S_ROH varied among breeds (Figure 2). The highest average S_ROH was 132 Mb in Cinisara, whereas the lowest one was 90 Mb in Reggiana. Considering the median values, the highest one was found in Italian Holstein, whereas the lowest one was found in Reggiana. The average reported S_ROH values were lower than the ones reported in other studies (Purfield et al., Reference Purfield, Berry, McParland and Bradley2012; Ferenčaković et al., Reference Ferenčaković, Hamzić, Gredler, Solberg, Klemetsdal, Curik and Sölkner2013a). The three most homozygous animals present in our dataset were from Cinisara (676.9 Mb), Modicana (681.2 Mb) and Reggiana (725.2 Mb) with almost a quarter of their genome classified as ROH. In all breeds, most ROH segment coverage was in the shorter length categories (4 to 8 Mb), in particular Modicana (51%) and Italian Holstein (50%) (Table 3). In fact, as reported in studies of ROH in human (Kirin et al., Reference Kirin, McQuillan, Franklin, Campbell, McKeigue and Wilson2010) and cattle populations (Ferenčaković et al., Reference Ferenčaković, Hamzić, Gredler, Solberg, Klemetsdal, Curik and Sölkner2013a; Marras et al., Reference Marras, Gaspa, Sorbolini, Dimauro, Ajmone-Marsam, Valentini, Williams and Macciotta2014) longer ROH were found less frequently than shorter ones. The expected length of autozygous segments that are IBD follows an exponential distribution with mean equal to 1/2g Morgans, where g is the number of generations since the common ancestor (Howrigan et al., Reference Howrigan, Simonson and Keller2011). Therefore, considering that 16 Mb segments are expected to present inbreeding up to three generations ago, recent inbreeding was observed in the studied local breeds due to the higher frequencies of ROH in this length category (Table 3), whereas the short ROH segments observed in Italian Holstein (4 Mb) was related to more ancient inbreeding, occurring 12.5 generation ago (about 75 years ago). However, the findings suggest that the local breeds experienced both recent and ancient inbreeding events, since that some animals lacked such long ROH, whereas other showed long segments. The results also indicated that these breeds have not recently been extensively crossed with other ones otherwise the long ROH would have broken down.

Figure 1 Relationship between the total number of runs of homozygosity (ROH) >4 Mb and the total length (kb) of genome in such ROH for individuals from each breed. Each dot represents an individual.

Figure 2 Box plots of within-breed average and median sum of all ROH segments per individual. ROH=runs of homozygosity; CIN=Cinisara; MOD=Modicana; HOL=Holstein; REG=Reggiana.

Table 2 Descriptive statistics for runs of homozygosity (ROH) for each cattle breed

MN_ROH=mean number of ROH per individual with minimum and maximum value in brackets; F _{ROH>4 Mb}=mean ROH-based inbreeding coefficient with standard deviation and minimum and maximum value in brackets; L_ROH=average length of ROH in Mb with minimum and maximum value in brackets; SNPs=minimum and maximum number of single nucleotide polymorphisms (SNPs) involved in ROH.

Table 3 Descriptive statistics of the number and the frequency distribution of runs of homozygosity (ROH) in different ROH length categories (Mb) for each cattle breed

n ROH=number of ROH; Freq=relative frequency of ROH on different ROH length categories.

One of the main advantages of genomic coefficients is the availability of chromosomal inbreeding coefficients. F _ROHBTA estimates were reported in Figure 3. In general, for each breed, the F _ROHBTA values followed the same pattern as those computed for the whole genome. Higher F _ROHBTA values were found on BTA28 (for Cinisara), BTA16 (Modicana), BTA26 (Italian Holstein) and BTA23 (Reggiana), whereas for all breeds the lowest one was found in BTA5. In a previous study on Italian Holstein, Gaspa et al. (Reference Gaspa, Marras, Sorbolini, Ajmone Marsan, Williams, Valentini, Dimauro and Macciotta2014) identified an interesting region of ~2 Mb on BTA26 that harbors some genes involved in the metabolism of mammary gland. Similar values were reported by Marras et al. (Reference Marras, Gaspa, Sorbolini, Dimauro, Ajmone-Marsam, Valentini, Williams and Macciotta2014) in Italian Simmental and local Marchigiana cattle breeds.

Figure 3 Distribution of inbreeding coefficient estimates for each chromosome (F _ROHBTA) calculated as the proportion of BTA in ROH over the length of the BTA covered by the involved SNPs. ROH=runs of homozygosity; CIN=Cinisara; MOD=Modicana; HOL=Holstein; REG=Reggiana.

In the absence of pedigree information, the origin of ROH could also be explained using other indicators, as LD and N _e . In fact, another explanation for ROH is the lack of recombination in a specific region. Pairwise r ² values were averaged over all autosomes and plotted as a function of genomic distance between markers (Figure 4). The highest level of r ² was found in Italian Holstein, whereas the lowest one in Cinisara. The extent of LD was used to estimate current and past N _e that is an important parameter for the assessment of genetic diversity and helps to explain how population evolved (Tenesa et al., Reference Tenesa, Navarro, Hayes, Duffy, Clarke, Goddard and Visscher2007). In the four breeds, the highest N _e (estimated five generation ago) was observed in Cinisara (94.58), whereas the lowest one was observed in Modicana (59.84) (Table 4). For Sicilian breeds, the N _e estimates based on LD were substantially higher than those reported in a previous study (Mastrangelo et al., Reference Mastrangelo, Saura, Tolone, Salces-Ortiz, Di Gerlando, Bertolini, Fontanesi, Sardina, Serrano and Portolano2014) calculated from the rates of F and f. Different estimates for N _e were also reported in Iberian pigs with complete and accurate pedigree records, where N _e calculated from the rates of molecular F and f were 17 and 10, respectively (Saura et al., Reference Saura, Fernández, Rodríguez, Toro, Barragán, Fernández and Villanueva2013), whereas N _e estimate using information from LD and recombination rate was 36 (Saura et al., Reference Saura, Tenesa, Woolliams, Fernández and Villanueva2014). Therefore, the discrepancies were due to the different used methods. In fact, as for the pedigree-based methods, the different molecular methods may give divergent results depending on the sampling strategy or the parameters used to compute N _e (Leroy et al., Reference Leroy, Mary-Huard, Verrier, Danvy, Charvolin and Danchin-Burge2013). These methods differ also in terms of time scale investigated and the amount of available information. The rates of F and f only give estimates of N _e based on limited time period, and taking into account the year of birth of individuals (that in local breeds as Cinisara and Modicana may be incorrect) may result in biased estimates. LD-based method uses more information, leads to an accurate estimate (Waples and Do, Reference Waples and England2010; Waples and England, Reference Waples and Do2011; Saura et al., Reference Saura, Woolliams, Tenesa, Fernández and Villanueva2015), with the possibility of investigating the change of N _e over time, as LD between loci at a specific recombination distance reflects the ancestral N _e 1/2c generations ago (Hayes et al., Reference Hayes, Visscher, McPartlan and Goddard2003), if the population grows linearly over time. However, it should be underlined that some parameters, as density and frequency of SNP pairs and distribution of MAF, affect the estimations of LD (Ober et al., Reference Ober, Malinowski, Schlather and Simianer2013) and then of N _e . Moreover, the methods used to convert physical distances between SNP pairs to genetic distance may result in different estimated N _e values (García-Gámez et al., Reference García-Gámez, Sahana, Gutiérrez-Gil and Arranz2012). Estimate of N _e obtained in this study for Italian Holstein was closed to those previously published for other Holstein population (Rodríguez-Ramilo et al., Reference Rodríguez-Ramilo, Fernández, Toro, Hernández and Villanueva2015). In general, the breed with the highest average inbreeding coefficient had the lowest N _e , as in Modicana breed. Moreover, LD and N _e were influenced by the recent history of selection. In fact, the strong selection for milk production and artificial insemination in Holstein and the highest inbreeding in Modicana have led to a reduction in the N _e .

Figure 4 Linkage disequilibrium across the genome as a function of genomic distance (Mb). CIN=Cinisara; MOD=Modicana; HOL=Holstein; REG=Reggiana.

Table 4 Effective population size (Ne) estimated from linkage disequilibrium values for each cattle breed

In Table 5 the average inbreeding and coancestry molecular coefficients estimated using different approaches were reported. Cinisara presented the highest values for all F coefficients (F _GRM, F _HOM and F _{MOL i} ); Modicana showed the lowest values for F _GRM and F _HOM and the highest value for f _{MOL ij} (Table 5). Italian Holstein breed showed the lowest values of f _{MOL ij} and F _{MOL i} . Estimates of inbreeding coefficients depend on the used methods. In fact, F coefficients estimated using allele frequencies (F _HOM and F _GRM) showed considerable variation among breeds respect to F _ROH and F _{MOL i} . In all breeds, f _{MOL ij} and F _{MOL i} values were much higher than the other coefficients because these two methods (that are obtained on a SNP-by-SNP basis) do not discriminate alleles that are IBD or identical by status (IBS) (Rodríguez-Ramilo et al., Reference Rodríguez-Ramilo, Fernández, Toro, Hernández and Villanueva2015). However, these estimates computed from SNP array data were strongly correlated with genealogical estimates, represent a useful alternative to genealogical information for measuring and maintaining genetic diversity and are very accurate in predicting genealogical coancestry (Gómez-Romano et al., Reference Gómez-Romano, Villanueva, de Cara and Fernández2013; Saura et al., Reference Saura, Fernández, Rodríguez, Toro, Barragán, Fernández and Villanueva2013). Spearman’s rank correlation between F _ROH and the other genomic inbreeding estimated measures was calculated (Table 6). High correlation was found between F _HOM and F _ROH ranged from 0.83 in Reggiana to 0.95 in Cinisara and Modicana. The correlations among F _ROH and other inbreeding estimates (F _GRM, F _HOM and F _{MOL i} ) were generally lower ranged from 0.45 (F _{MOL i,}−F _ROH) in Cinisara to 0.17 (F _GRM−F _ROH) in Modicana (Table 6). High correlation between F _HOM and F _ROH (0.84) was also reported by Zhang et al. (Reference Zhang, Young, Wang, Sun, Wolc and Dekkers2014) in a study on pig in which ROH >5 Mb after LD pruning were detected, whereas really different values (0.06, 0.35 and 0.61) were reported by Zhang et al. (Reference Zhang, Calus, Guldbrandtsen, Lund and Sahana2015) in three cattle breeds. Ferenčaković et al. (Reference Ferenčaković, Hamzić, Gredler, Solberg, Klemetsdal, Curik and Sölkner2013a) reported high correlation between F _HOM and F _ROH based on short segments (ROH >1 and >2 Mb). The poor correlation reported in our study between F _GRM and F _ROH was according to other studies (Marras et al., Reference Marras, Gaspa, Sorbolini, Dimauro, Ajmone-Marsam, Valentini, Williams and Macciotta2014; Zhang et al., Reference Zhang, Calus, Guldbrandtsen, Lund and Sahana2015). Zavarez et al. (Reference Zavarez, Utsunomiya, Carmo, Neves, Carvalheiro, Ferenčaković, O’Brien, Curik, Cole, Van Tassell, da Silva, Sonstegard, Sölkner and Garcia2015) in a study on autozygosity using high-density SNPs, showed that the correlation between F _GRM and F _ROH decreased from 0.74 per ROH >0.5 Mb to 0.41 per ROH >16 Mb, probably due to the properties of the G matrix which is based on individual loci, whereas F _ROH is based on chromosomal segments. A higher correlation between F _{MOL i} and F _ROH were reported by Gómez-Romano et al. (2014) in Austrian Brown Swiss cattle (0.76) and Rodríguez-Ramilo et al. (Reference Rodríguez-Ramilo, Fernández, Toro, Hernández and Villanueva2015) in Spanish Holstein breed (0.88). However, while the alternative used estimates of inbreeding and coancestry coefficients could not distinguish between recent and ancient inbreeding, F _ROH provided the direct estimated level of autozygosity in the current populations and allowed us to detect recent inbreeding (up to three generations ago) in the local cattle breeds, in particular for Cinisara and Modicana ones. In fact, in the Sicilian farming system, natural mating is the common practice for local breeds, and the exchange of animal among flocks is quite unusual, with an increase of inbreeding within the population due to uncontrolled mating of related individuals (Mastrangelo et al., Reference Mastrangelo, Sardina, Riggio and Portolano2012). As pedigree data were unavailable for animals in this study, comparison of genomic and pedigree inbreeding coefficients was not possible. However, the strong correlation between the pedigree inbreeding coefficient and the sum of ROH reported by several authors (Purfield et al., Reference Purfield, Berry, McParland and Bradley2012; Ferenčaković et al., Reference Ferenčaković, Solkner and Curik2013b) suggests that in absence of animal’s pedigree data, the extent of a genome under ROH may be used to infer aspects of recent population history even from relatively few samples. It should be underlined that the occurrence of ROH in an individual may be the result of inbreeding events but they may also be present in outbreed populations as result of other phenomena. In fact, an increased frequency of common extended haplotypes can also be a consequence of selection pressure on genomic regions involved in functional roles (Gaspa et al., Reference Gaspa, Marras, Sorbolini, Ajmone Marsan, Williams, Valentini, Dimauro and Macciotta2014), but as reported above, Sicilian cattle breeds are not subject to selection programs, therefore the presence of ROH in these two breeds was only due to inbreeding effect. Moreover, recent studies showed that the genomic estimates of inbreeding can be used to calculate the effects of inbreeding on performance and fitness traits. Pryce et al. (Reference Pryce, Haile-Mariam, Goddard and Hayes2014), in a study on the identification of genomic regions associated with inbreeding depression in Holstein cattle breed, showed that long ROH (>60 SNPs or 3.5 Mb), as those identified in our breeds, were associated with a reduction in milk yield, independently of the proportion of the genome that was homozygous. Therefore, our results showed the necessity of implementing conservation programs to preserve the local breeds in order to avoid further loss of genetic distinctiveness.

Table 5 Estimated mean of genomic inbreeding and coancestry coefficients for each cattle breed

F _GRM=inbreeding coefficient based on genomic relationship matrix; F _HOM=inbreeding coefficient based on the difference between observed v. expected number of homozygous genotypes; F _{MOL i} =molecular inbreeding coefficient of individual i; f _{MOL ij} =molecular coancestry coefficient between individuals i and j.

Table 6 Correlation between F _ROH and other genomic inbreeding coefficients for each cattle breed

F _HOM=inbreeding coefficient based on the difference between observed v. expected number of homozygous genotypes; F _ROH=inbreeding coefficient based on the runs of homozygosity; F _GRM=inbreeding coefficient based on genomic relationship matrix; F _{MOL i} =molecular inbreeding coefficient of individual i.

*P<0.05, **P<0.01, ***P<0.001.

Selection and mating strategies have been proposed in the past for controlling inbreeding and coancestry. The best know strategy to achieve these goals is optimizing the contributions of the parents to minimize global coancestry in their offspring (Fernández et al., Reference Fernández, Meuwissen, Toro and Mäki-Tanila2003). Recently, measures of coancestry based on IBD segments (de Cara et al., Reference de Cara, Villanueva, Toro and Fernández2013) and on shared segments of the genome (Bosse et al., Reference Bosse, Megens, Madsen, Crooijmans, Ryder, Austerlitz, Groenen and de Cara2015) have been proposed as good balance between maintaining diversity and fitness, with a higher fitness than managing with molecular coancestry and higher diversity than managing with genealogical coancestry. Therefore, determining the occurrence of IBD segments in potential parents, thereby measuring their relatedness and coancestry, can be used to minimize the occurrence of long ROH in the offspring. The availability of genome-wide genotyping platforms allows us now to study populations from a more detailed perspective, providing information on the genetic status and on their evolution across time.