Genomic inbreeding estimation in small populations: evaluation of runs of homozygosity in three local dairy cattle breeds

Abstract

In the local breeds with small population size, one of the most important problems is the increase of inbreeding coefficient (F). High levels of inbreeding lead to reduced genetic diversity and inbreeding depression. The availability of high-density single nucleotide polymorphism (SNP) arrays has facilitated the quantification of F by genomic markers in farm animals. Runs of homozygosity (ROH) are contiguous lengths of homozygous genotypes and represent an estimate of the degree of autozygosity at genome-wide level. The current study aims to quantify the genomic F derived from ROH (F_ROH) in three local dairy cattle breeds. F_ROH values were compared with F estimated from the genomic relationship matrix (F_GRM), based on the difference between observed v. expected number of homozygous genotypes (F_HOM) and the genomic homozygosity of individual i (F_{MOL i}). The molecular coancestry coefficient (f_{MOL ij}) between individuals i and j was also estimated. Individuals of Cinisara (71), Modicana (72) and Reggiana (168) were genotyped with the 50K v2 Illumina BeadChip. Genotypes from 96 animals of Italian Holstein cattle breed were also included in the analysis. We used a definition of ROH as tracts of homozygous genotypes that were >4 Mb. Among breeds, 3661 ROH were identified. Modicana showed the highest mean number of ROH per individual and the highest value of F_ROH, whereas Reggiana showed the lowest ones. Differences among breeds existed for the ROH lengths. The individuals of Italian Holstein showed high number of short ROH segments, related to ancient consanguinity. Similar results showed the Reggiana with some extreme animals with segments covering 400 Mb and more of genome. Modicana and Cinisara showed similar results between them with the total length of ROH characterized by the presence of large segments. 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 estimated F coefficients were generally lower ranged from 0.45 (F_{MOL i}−F_ROH) in Cinisara to 0.17 (F_GRM−F_ROH) in Modicana. On the basis of our results, recent inbreeding was observed in local breeds, considering that 16 Mb segments are expected to present inbreeding up to three generations ago. Our results showed the necessity of implementing conservation programs to control the rise of inbreeding and coancestry in the three Italian local dairy cattle breeds.

Implications

In the local breeds with small population size, one of the most important problems is the increase of inbreeding that leads to different negative effects as a reduction in phenotypic values. The current study aims to quantify the genomic inbreeding derived from runs of homozygosity (ROH) (F _ROH) in three Italian local dairy cattle breeds. According to ROH results, recent inbreeding was well detected in the investigated local dairy cattle breeds. Our results showed the necessity of implementing conservation programs to preserve the local breeds in order to avoid further loss of genetic distinctiveness. Therefore, determining the occurrence of identical by descent segments in potential parents, thereby measuring their relatedness and coancestry, can be used to minimize the occurrence of long ROH in the offspring.

Introduction

Animal genetic resources must be preserved because of their contribution to human livelihood, now and in the future (Toro et al., 2011). Most local livestock breeds are the result of a particular adaptation to production systems environmentally conditioned, and in many cases no other breed could survive in the same habitat if the local breed goes extinct. In addition, such local populations might harbor specific genetic variants that are worth retaining and that might be used to recover the loss of genetic diversity that occurs in mainstream breeds because of very intensive selection on production traits (Fernández et al., 2011). Apart from that, these populations represent local culture, history and tradition and are often linked to traditional products of farm animals (milk, meat, eggs, etc.).

Typically, local breeds are small populations and their size put them at risk of extinction. Consequently, the genetic diversity stored in each of them should be treated with great care and management strategies that insure the viability, and maintenance of the population should be implemented. Selection programs in local breeds with small population size are limited by the low number of animals (families) and by the need to control inbreeding (Fontanesi et al., 2015), which represent one of the most important problems. The individual inbreeding coefficient (F) is defined as the proportion of an individual’s genome that is autozygous, that has homozygous identical by descent (IBD) status, or equivalently the probability of a randomly sampled locus in the genome to be autozygous (Ferenčaković et al., 2013a). The increase of F leads to different negative effects as reduction in phenotypic values for some traits, reduction of genetic variance, higher frequency of homozygous genotypes with the reduction of individual performance (inbreeding depression) and lower population viability (Ouborg et al., 2010). Therefore, to avoid inbreeding depression, an accurate and sensitive estimation of F is very important, especially in local breeds/populations. Traditionally, F is estimated on the basis of pedigree information but in most cases this is unavailable or inaccurate. Moreover, the probabilistic approach of pedigree analysis does not take into account the stochastic nature of recombination (McQuillan et al., 2008). Recently, with the availability of high-density single nucleotide polymorphism (SNP) arrays, F can be estimated accurately in absence of pedigree information (Allendorf et al., 2010). There are two categories of genomic inbreeding measures based on genome-wide SNPs. The first category is based on marker-by-marker estimates such as the diagonal elements of the genomic relationship matrix (GRM) (VanRaden et al., 2011), the canonical estimate based on excess SNP homozygosity in PLINK (Purcell et al., 2007) and molecular coancestry estimates (Caballero and Toro, 2002). The second one is based on runs of homozygosity (ROH) detection. ROH are contiguous lengths of homozygous genotypes that are present in an individual due to parents transmitting identical haplotypes to their offspring (Gibson et al., 2006). Nowadays, F estimated from ROH (F _ROH) is considered to be the most powerful method of detecting inbreeding effects among several alternative estimates of inbreeding (Keller et al., 2011). F _ROH provided a good measure of individual genome-wide autozygosity and allows to distinguish between recent and ancient inbreeding (McQuillan et al., 2008). Because recombination events interrupt long chromosome segments, long ROH (~10 Mb) arise as result of recent inbreeding (up to five generation ago), while shorter ROH (~1 Mb) can indicate more distant ancestral effect (up to 50 generation ago) such as breed founder effects (Howrigan et al., 2011). Therefore, estimate of F using ROH is particularly appealing as the number of generations of inbreeding and the history of recent selection events can be inferred from the extend and frequency of ROH regions (Purfield et al., 2012). Although ROH from high-throughput genotyping analyses have been studied extensively in humans, these estimates are rare in cattle, particular in local breeds, and in other livestock species (Purfield et al., 2012; Ferenčaković et al., 2013a; Silió et al., 2013; Pertoldi et al., 2014).

The current study aims to quantify the genomic inbreeding derived from ROH in three economically important Italian local dairy cattle breeds, Cinisara, Modicana and Reggiana, characterized by the same breeding goals but different selection histories. Moreover, genotypes from Italian Holstein, the most important dairy cattle breed reared in Italy, were also included in these analyses in order to compare results among breeds.

Material and methods

Breeds, genotypes and quality control

A total of 407 individuals were used for the analyses. DNA samples belonged to four different cattle breeds: Cinisara (71), Modicana (72), Reggiana (168) and Italian Holstein (96). For these breeds pedigree data were not available. Sampling was carried out in several farms and individuals were selected on the basis of information supplied by the farmers to avoid, as much as possible, closely related animals. The Cinisara, Modicana and Reggiana are three economically important local breeds with small population size (number of reared animals <4000). Cinisara and Modicana are two cattle breeds well adapted to the harshness of Sicilian marginal mountain areas and their economic importance lies on the traditional production systems of two typical ‘pasta filata’ cheeses: Caciocavallo Palermitano and Ragusano PDO (Protected Designation of Origin), respectively. Recently, Mastrangelo et al. (2014) reported the effective population size values estimated from rate of F per year (19 and 12) and from rate of coancestry (f) (four and eight individuals) in Cinisara and Modicana cattle breeds, respectively. Reggiana is a local cattle breed reared in the province of Reggio Emilia in Northern Italy specialized for the production of a niche brand of Parmigiano-Reggiano PDO cheese.

All animals were genotyped for 54 609 SNPs using Bovine SNP50K v2 BeadChip (Illumina Inc., San Diego, CA, USA). Data quality control was performed separately for each breed. We excluded all SNPs not assigned to a bos taurus chromosome (BTA) or assigned to chromosomes X and Y. Markers were filtered according to quality criteria that included: (i) call frequency (⩾0.95), (ii) minor allele frequency (MAF⩾0.01) and (iii) Hardy-Weinberg equilibrium (P-value=0.001). SNPs that did not satisfy these quality criteria were excluded. Moreover, considering that high linkage disequilibrium (LD) can lead to detection of ROH that are not truly IBD, LD pruning was also performed before the ROH call to increase power, as suggested by Purcell et al. (2007) and applied by several authors (Howrigan et al., 2011; Bjelland et al., 2013). Therefore, unlinked SNPs were selected using -indep option of PLINK with the following parameters: 50 SNPs/window, a shift of five SNPs between windows and r ² threshold of 0.5. A total of 38 937 SNPs in Cinisara, 32 179 SNPs in Modicana, 29 483 SNPs in Reggiana and 27 586 SNPs in Italian Holstein cattle breeds were retained after quality control and were used to estimate F _ROH. The main difference for the number of SNPs used for each breed, in particular the highest number of SNPs used for Cinisara, was due to different values of LD among breeds. In fact, Cinisara showed the lowest value of LD and, therefore, the lowest number of excluded SNPs.

Run of homozygosity calling option

F _ROH were calculated as the proportion of genome in ROH over the overall length of the genome covered by the involved SNPs (2 541 174 kb) using the PLINK whole-genome association analysis toolset (Purcell et al., 2007). The following criteria were used to define the ROH: (i) the minimum number of SNPs included in the ROH was fixed to 40; (ii) the minimum length that constituted the ROH was set to 4 Mb; (iii) two missing SNPs were allowed in the ROH; (iv) minimum density of one SNP every 100 kb; (v) maximum gap between consecutive SNPs of 1 Mb. Moreover, the number of allowed heterozygous SNPs was set to different values: from one to three. Mean F _ROH values obtained allowing different numbers of heterozygous SNPs were compared within the same breed using paired t-tests. The mean number of ROH per individual per breed (MN_ROH), the average length of ROH (L_ROH) and the sum of all ROH segments per animal (S_ROH) were estimated. The distribution of S_ROH within breed was assessed using box plots. In addition, chromosomal (BTA) F _ROH (F _ROHBTA) values were also estimated for each breed, as F _ROHBTA=L_ROHBTA/L_BTA (Silió et al., 2013), in which L_ROHBTA is the total length of an individual’s ROH in each BTA and L_BTA is the length of each chromosome covered by the involved SNPs (Supplementary Table S1). ROH were classified into three classes (4 to 8, 8 to 16 and >16 Mb) using the same nomenclature reported by other authors (Ferenčaković et al., 2013a; Marras et al., 2014) except for two classes (<2 and 2 to 4 Mb), which were not considered in our study. The number and percentage of ROH within each ROH length category for breed were also determined.

Genomic inbreeding analyses

Alternative estimates of inbreeding and coancestry coefficients were also calculated. In particular: (1) the values of the diagonal elements of the GRM (F _GRM) proposed by VanRaden et al. (2011); (2) the genomic inbreeding coefficient based on the difference between observed v. expected number of homozygous genotypes (F _HOM) using PLINK (Purcell et al., 2007); (3) the molecular coancestry coefficient (f _MOLij ) between individuals i and j (Caballero and Toro, 2002); (4) the molecular inbreeding coefficient (F _MOLi ) of individual i, calculated as F _{MOL i} =2 f _{MOL ii} −1 (f _{MOL ii} is the molecular self-coancestry). Spearman’s rank correlation among different genomic inbreeding measures was calculated.

Effective population size

The effective population sizes (N _e ) were calculated as N _e =(1/4c)×(1/r ²−1) (Sved, 1971) where r ² (the squared correlation coefficient of allele frequencies at pair of loci) is the value of LD and c the genetic distance in Morgans between SNPs. Physical distances between SNP pairs were converted to genetic distances with the assumption of 1 cM~1 Mb. Each genetic distance c corresponds to a value of t generation in the past, and this value was calculated as t=1/(2c), assuming a linear population growth (Hayes et al., 2003). All pairwise combinations of SNPs were estimated using LD plot function in Haploview v 4.2 software (Barrett et al., 2005). For this analysis, markers were filtered according to quality criteria reported above, except for LD pruning; in fact N _e estimates could be biased if calculated from LD pruned SNPs. A total of 44 875 SNPs in Cinisara, 42 687 SNPs in Modicana, 35 720 SNPs in Reggiana and 41 596 SNPs in Italian Holstein cattle breeds were used. For each chromosome, pairwise r ² was calculated for SNPs between 0 and 50 Mb apart. To visualize the LD pattern per chromosome, r ² values were stacked and plotted as a function of inter-marker distance categories.

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. (2013b) 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., 2013b) than in human species (Kirin et al., 2010). 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. (2013b) 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. (2014), 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., 2014), 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., 2015). 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., 2015).

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. (2013) 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., 2014). 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., 2013). Ferenčaković et al. (2013a) 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. (2013b) using a 50K panel for Pinzgauer (0.037) and Tyrol Grey (0.042) local cattle breeds, and by Marras et al. (2014) 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., 2012; Ferenčaković et al., 2013a). 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., 2010) and cattle populations (Ferenčaković et al., 2013a; Marras et al., 2014) 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., 2011). 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. (2014) 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. (2014) 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., 2007). 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., 2014) 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., 2013), whereas N _e estimate using information from LD and recombination rate was 36 (Saura et al., 2014). 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., 2013). 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, 2010; Waples and England, 2011; Saura et al., 2015), 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., 2003), 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., 2013) 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., 2012). 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., 2015). 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., 2015). 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., 2013; Saura et al., 2013). 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. (2014) 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. (2015) in three cattle breeds. Ferenčaković et al. (2013a) 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., 2014; Zhang et al., 2015). Zavarez et al. (2015) 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. (2015) 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., 2012). 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., 2012; Ferenčaković et al., 2013b) 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., 2014), 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. (2014), 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., 2003). Recently, measures of coancestry based on IBD segments (de Cara et al., 2013) and on shared segments of the genome (Bosse et al., 2015) 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.

Conclusion

This study has reported for the first time the genome-wide inbreeding estimate using ROH in three Italian local cattle breeds. The obtained results highlight differences in detection and in distribution of ROH among breeds. In particular, Cinisara and Modicana breeds showed long ROH segments and the presence of inbreeding due to recent consanguineous mating. Therefore, our results showed the necessity of implementing conservation programs with the aim to control the level of inbreeding. The control of coancestry would restrict inbreeding depression, the probability of losing beneficial rare alleles and therefore the risk of extinction for these local cattle breeds, and may be crucial for implementing genetic improvement programs. Breeders should be aware of this situation, and breeding systems should be designed to foster and maintain genetic variation in these populations. Avoiding mating among relatives, together with other actions (e.g. sires/dams ratio, balanced progeny sizes) are strategies to control the increase of inbreeding.