Evolution of inbreeding: a gaze into five Italian beef cattle breeds history

Introduction

Inbreeding is the main consequence of mating individuals that are related, through common ancestry, to a degree that exceeds that of two individuals from the same population extracted at random (Kardos et al., 2016). Minimizing inbreeding is an aspect of paramount importance in cattle breeding to avoid the phenotypic expression of detrimental alleles in the offspring as well as to ensure the maintenance of genetic diversity (Howard et al., 2017). Inbred animals display chromosome segments that are identical-by-descent (IBD) and generate long runs of homozygosity (Kardos et al., 2016). Although inbreeding coefficients can be calculated from pedigrees comprising several generations, in recent years the advent of high throughput genotyping techniques has make it possible to estimate molecular inbreeding coefficients based on the characterization of the genome-wide patterns of homozygosity. One key advantage of this latter approach is that it captures ancient inbreeding accumulated in the base population and it is less affected by parentage errors (Howard et al., 2017). However, several studies (Bjelland et al., 2013; Pryce et al., 2014) have revealed the distortion caused by ancient contributions to inbreeding dilutes up to a degree in which pedigree-based analyses and genomic analyses may not differ that much.

Inbreeding is often reported in domestic animal populations as a measurement of one or several coefficients of inbreeding at a particular time point. However, assessing the magnitude of inbreeding on a continuous temporal scale is much more informative because it captures its tendency and predicted behaviour. Between 1960 and 2000, inbreeding coefficients of US dairy breeds, such as Ayrshire, Brown Swiss, Guernsey, Holstein, and Jersey, went from 0% (base population) to 4.5–6% in just four decades (Weigel, 2001). According to Weigel (2001), this increase in inbreeding was not associated with effective population size and, more likely, it was the result of the intensity of genetic selection as well as of the extensive use of a reduced number of elite sires. Similarly, Mc Parland et al. (2007) investigated the evolution of inbreeding in the Charolais, Limousine, Hereford, Angus, and Simmental beef cattle breeds as well as in the Holstein-Friesian dairy cows raised in Ireland. They found that, between 1960 and 2004, overall inbreeding increased from 0.10–0.25% to 0.5–2%, and over the last decade (1994–2004) the annual rate of increase in inbreeding was 0.06–0.13%. These and other studies suggest that inbreeding is accumulating rapidly in cattle breeds due to efficient genetic selection programs and reproductive management (Weigel, 2001).

Genetic selection in Italian beef cattle is implemented by the National Association of Italian Beef Cattle Breeders (ANABIC) and aimed to improve meat production, precocity, growth ability, and muscle development (Sbarra, 2011). Three of the five main Italian beef cattle breeds, Marchigiana (MAR), Chianina (CHI), and Romagnola (ROM) are highly specialised in beef production, and the other two, Maremmana (MRM) and Podolica (POD), are considered as rustic (Sarti et al., 2019). The current selection program, based on the traditional quantitative approach, has achieved a remarkable improvement of growth, daily weight and muscularity gain (Sbarra, 2011). Moreover, cattle are somatically well-developed with a correct morphology and light skeletal system (Rovelli et al., 2020a). However, the intensity of the selection in the five breeds is lower in the rustic MRM and POD than in the three specialised ones; moreover, the two rustic breeds register a low amount of young bulls/year in performance test (Fioretti et al., 2020). In the current work, we aimed to characterize the historical evolution of inbreeding in five Italian beef cattle breeds (CHI, MAR, ROM, MRM, and POD) in the period comprised between 2002 and 2019 by using both molecular and genealogical estimates of inbreeding coefficients. Our goal was to test whether intensive selection performed in the last twenty years has resulted in a significant increase of inbreeding levels in these five populations.

Materials & Methods

Cattle sampling and genotyping

This work comprised 3,581 young bulls belonging to five of the main Italian beef cattle breeds: CHI (909), MAR (879), ROM (904), MRM (334), and POD (555). The number of young bulls and their year of birth are reported in Table S1. Samples were collected by staff from the ANABIC at the genetic station of San Martino in Colle (Perugia, Italy) during 1985–2019. Geographical distribution of these breeds in Italy and the pictures of representative individuals from each breed can be found in Fig. 1.

During the performance test, blood samples were taken from the jugular veins of each bull. These samples were collected in EDTA K₃ coated vacuum tubes and stored at −20 °C (Rovelli et al., 2020b). Genomic DNA was purified with the GenElute Blood Genomic DNA kit (Sigma Aldrich, St. Louis, MO, USA). The purification method was previously described in Sarti et al. (2019). Genomic DNA samples from the 3,581 bulls were genotyped with the GeneSeek Genomic Profiler Bovine LDv4 33K chip (Illumina Inc., San Diego, CA, USA), which contains 30,111 SNPs, and processed at the Agrotis laboratory (LGS, Cremona, Italy). Standard multi-sample protocols and reagents were used according to the instructions of the manufacturer (Khatkar et al., 2012).

The positions of single nucleotide polymorphisms (SNPs) were mapped using the ARS-UCD_1.2 bovine genome assembly (Zorc, Ogorevc & Dovc, 2019). The PLINK software v1.9 (Chang et al., 2015) was used to update the names and positions of SNP markers. The SNPs that did not match the following criteria were removed before performing population diversity analyses: (i) SNPs with call rates less than 90%, (ii) SNPs with minor allele frequencies less than 5%, (iii) SNPs with more than 1% missing genotypes, and (iv) SNPs displaying highly significant deviations (P-value < 10⁻³) from the Hardy-Weinberg equilibrium (Amaral, Pavão & Gama, 2020). Another pruning step was performed to remove SNPs with high linkage disequilibrium (LD) using the command –indep 50 5 2 (Manunza et al., 2016) of the PLINK v1.9 software (Chang et al., 2015) as recommended in a previous publication (Howrigan, Simonson & Keller, 2011). This second step is necessary because stretches of SNPs with low MAF and genomic regions with many SNPs and strong LD often yield erroneous estimates of the effective population size (Manunza et al., 2016). The results of these pruning steps are shown in Table S2.

Results

The average PCI from 2002 to 2019 ranged from 99.29 to 99.91% (MAR), from 99.17 to 99.89% (CHI), from 98.80 to 99.92% (ROM), from 99.01 to 99.87% (MRM), and from 99.02 to 99.86% (POD). The removal of genotyped animals with PCI less than 90% resulted in the exclusion of less than 1.5% of the sample, as most of the genotyped animals had PCI greater than 90%. The estimates of Ho, He, F_hat2, and F_ped averaged across years (2002-2019) in the studied breeds fluctuated between 0.340–0.401, 0.348–0.392, −0.121–0.072, and 0.000–0.068, respectively (Tables 1 and 2). The F_ped coefficients were higher in the rustic POD and MRM breeds, probably because the depth of the pedigree (in our dataset) is higher than in the MAR, CHI, and ROM breeds. Moreover, POD was the breed that displayed the highest F_hat2 coefficient, followed by ROM. The analysis of the evolution of F_hat2, F_ped, and Ho in the five breeds (Figs. 2 and 3) evidenced that the observed heterozygosity, in the five studied populations, remained constant throughout the years, with slightly higher values in the rustic breeds. With regard to F_hat2 and F_ped, we observed some fluctuations across years that were particularly accentuated for F_hat2 which showed a stable or increasing trend depending on the breed under consideration. In any case, these yearly oscillations in the magnitude of inbreeding were not very important. The MRM and POD breeds lacked data in one and three years, respectively. These missing values are due to the fact that in Southern Italy, between 2004–2006, there was a “Bluetongue” epidemia which caused a temporary cessation of the activities of the POD selection center. For the same reason, MRM selection was temporarily suspended in 2013. The annual rates of increase in inbreeding (F_hat2 and F_ped) are displayed in Figs. 4 and 5. It can be seen that in general F_hat2 increases slightly but steadily in all five breeds, with averaged overall increasing rates of 0.17–0.34% between 2002–2019. In contrast, F_ped remained quite stable across time with averaged overall increasing rates of −0.04–0.08%. The effective population size (N_e) estimated t generations ago (from 13 to 235) is shown in Fig. 6. It is apparent that N_e decreases markedly across generations. Thirteen generations ago, N_e was lower than 300 for most breeds with the only exception of POD cattle (N_e = 498). In contrast, 235 generations ago N_e oscillated between 1887–3257, which is 7.48 times larger than current values.

Discussion

In this study, we have evaluated the variation of inbreeding and diversity parameters in five Italian beef cattle breeds across a window of approximately 20 years. We did not measure F_roh because the number of SNPs was too small (21,000–23,000 valid SNPs) to map runs of homozygosity in a reliable manner. Importantly, in our study the base generation of each breed was composed of animals born before 2002. We did so because there was a very small number of bulls representing each one of the years comprised between 1985 and 2002, so they were merged together in a single group. We have observed that the averaged (across years) coefficients of inbreeding F_hat2 and F_ped in the five Italian breeds are generally lower than 0.07 (Table 2, Figs. 2–5), while homozygosity was approximately 0.64–0.70 (data not shown). The F_hat2 coefficient is closely related to F_is and can be interpreted (also in the case of F_ped computed F_is) in breeding policy terms: in fact, positive and high values mean that mating between close relative are not –or cannot be –avoided. In smaller populations, like those at issue, breeders are implementing a more rigorous breeding policy to limit inbreeding within the herd. The MRM and POD breeds displayed the largest F_ped coefficients but this was expected because the depth of the pedigree, in our study, is larger than for CHI, MAR, or ROM. Another reason explaining the larger F_ped (and lower AR) observed in the rustic breeds (MRM and POD) relies on the fact that in these breeds artificial insemination (AI) is little spread. Bulls, especially in the past, were the offspring of animals born in the same farm. Furthermore, farms always remained quite isolated and poorly genetically connected to each other, at least until the foundation of test stations (Moioli, Napolitano & Catillo, 2004). On the other hand, in the specialized breeds (MAR, CHI, and ROM), the implementation of AI involved the use of unrelated lines thus avoiding inbreeding to a great extent (higher AR). We have observed slight discrepancies between F_hat2 and F_ped values displayed in Table 2, but this outcome is probably explained by the fact that these two coefficients have different properties (Alemu et al., 2021). Indeed, F_ped indicates the probability that two homologous alleles in an individual are identical by descent, as defined by Malécot (1948), and it ranges from 0 to 1. In contrast, F_hat2 is very similar to the method-of-moments F coefficient measured with the –het command of PLINK v1.9 (Purcell et al., 2007) which estimates the reduction of heterozygosity (or the excess of homozygosity) associated with inbreeding (Alemu et al., 2021). In consequence, F_hat2 can take negative values (Alemu et al., 2021). Despite these conceptual differences, comparison of F_ped and F_hom (similar to F_hat2) coefficients calculated in a pedigree of 245 Holstein cattle with whole-genome sequence data showed a good agreement between both parameters (Alemu et al., 2021). In our study, both coefficients indicate that none of the Italian breeds under study is significantly inbred.

We have also observed a lack of relationship between population size and the magnitude of inbreeding coefficients (Table 2, Figs. 2–5). For instance, MRM (F_hat2 = 0.014, F_ped = 0.061) and ROM (F_hat2 = 0.023, F_ped = 0.023) have population sizes of approximately 11,000–12,000 individuals, while MAR has a census almost five times larger but fairly comparable levels of inbreeding (F_hat2 = 0.012, F_ped = 0.018). The amount of inbreeding is mostly determined by the demographic history of populations rather than by their current size. In local breeds undergoing strong demographic reductions, genetic drift can be quite intense thus increasing homozygosity and the occurrence of matings between related individuals. For instance, Chillingham cattle (a breed that lives in Northumberland earldom, England) are currently represented by a herd of 50 males and 50 females, which has remained reproductively closed in the last 300–350 years (Visscher et al., 2001). This herd was formed by five bulls and eight cows in 1947, and the average number of males and females per generation has been three and 15, respectively (Visscher et al., 2001). Calculation of F_is inbreeding coefficient in Chillingham cattle yielded a value of 0.92, which is extremely high (Williams et al., 2016). In strong contrast, the population sizes of the five Italian breeds analyzed in the current work have remained relatively stable in the last 10 years and, to the best of our knowledge no genetic bottlenecks have been recorded. Interestingly, Weigel (2001) also observed a lack of correlation between the magnitude of inbreeding and the population size of five major US dairy breeds (Ayrshire, Brown Swiss, Guernsey, Holstein, and Jersey). Noteworthy, breeds with small (Jersey) and large (Holstein) population sizes had comparable levels of inbreeding, probably because in both populations the intensity of selection was considerably high (Weigel, 2001). We have compared the F_hat2 values measured in MAR, ROM, CHI, MRM, and POD breeds with F_hom values reported by Mastrangelo et al. (2018) in the same breeds. As previously explained, these two inbreeding coefficients, generated with the PLINK v1.9 software (Purcell et al., 2007) are quite comparable and both measure the reduction in heterozygosity associated with inbreeding. We observed that in general the F_hom coefficients measured by Mastrangelo et al. (2018) in the same five breeds analyzed by us are higher (in the 0.066–0.118 range) than the F_hat2 coefficients estimated in the current work. Biscarini et al. (2020) reported an F_ped coefficient for the MRM breed of 0.049, which is similar to the result obtained by us (Table 2). In another study, Moioli, Napolitano & Catillo (2004) described Wright F_is coefficients of 0.106 and 0.138 in the POD and MRM breeds, and Santana Jr et al. (2012) measured an F_ped = 0.013 in Brazilian MAR cattle. One potential reason for the discrepancy between our results and those of Mastrangelo et al. (2018) is that we have estimated allele and genotype frequencies with much more accuracy because our sample size is much larger i.e., Mastrangelo et al. (2018) used samples sizes of 21–24 animals while we have used sample sizes of 334-909 animals. This can be particularly important when measuring the excess of homozygosity in SNP markers with low or very low minimum allele frequencies.

Despite the fact that we have detected fluctuations in the magnitudes of F_hat2 and F_ped during the last 20 years, these oscillations were not very important (Figs. 2 and 3). Changes in homozygosity and inbreeding coefficients across time might be partly explained by the fact that the number of breeders, and particularly sires, is not constant across years. We have also observed a slight increase of F_hat2 in the five breeds during 2002–2020, particularly in MRM, POD, and CHI, increasing 0.01–0.02% per year, while F_ped remained quite stable, with an annual rate of increase of 0.003–0.004% per year (Figs. 4 and 5). Mc Parland et al. (2007) measured the evolution of the inbreeding in five Irish cattle and observed that in the 1994–2004 period F_ped remained constant or decreased in Angus, Charolais, and Limousine populations, while in Hereford, Holstein-Friesian, and Simmental a yearly increase of F_ped (between 0.06–0.13%) was detected. In another study, Weigel (2001) also reported increases of F_ped of 0.10–0.15% per year. In the case of the five Italian breeds under study, such increases are much more modest, because the ANABIC genetic program is designed to minimize inbreeding. A high number of bulls is used in natural reproduction, rather than a few elite sires, and each one of them is mated with few cows, because this approach ensures an offspring with low inbreeding coefficients. The matings are programmed and designed to minimize inbreeding, using the less related pairs of breeders and allowing a maximum inbreeding increment of 5% in each mating (Sbarra, 2011).

By using the SNeP (Barbato et al., 2015) and Endog v4.8 (Gutiérrez & Goyache, 2005; Gutiérrez, Goyache & Cervantes, 2009) tools, we have detected a sustained and marked decline in the effective population size of the five breeds under investigation. Mastrangelo et al. (2018) reported a very similar tendency in the same five populations. Effective population size is a complex parameter defining the size of a Wright–Fisher ideal population generating the same rate of inbreeding and variance of gene frequencies detected in the real population under investigation (Crow & Kimura, 1970). In principle, selection and reproductive management, particularly in AI schemes in which a reduced number of bulls mate with a large number of cows, are expected to reduce N_e. Although, N_e cannot be equated to a coefficient of inbreeding or diversity because it depends on many variables (Wang, Santiago & Caballero, 2016), it is remarkable that this strong N_e decline (Fig. 6) was not accompanied by a substantial reduction in Ho or He in the last 20 years (Figs. 2 and 3). Increased variance in family size associated with the upward trend to use a reduced number of elite sires as breeders could be one of the major reasons for this progressive decline of N_e. However, it should be also noticed that estimates of N_e historical trajectories with the SNeP tool (Barbato et al., 2015) are sometimes unreliable, particularly in the most recent and oldest generations (Corbin et al., 2012). Besides, N_e estimates are strongly affected by data manipulation factors (e.g., choice of the minimum allele frequency threshold) employed in the analysis (Corbin et al., 2012; Barbato et al., 2015) and one of the main tenets of the coalescent is that no SNP can be reliably sampled after 4N_e generations in the past. So, our estimates of N_e historical trajectories should be interpreted with caution.

Conclusions

The low level of inbreeding found in this study confirms the success of the Italian beef cattle selection program carried out by ANABIC, which aimed to minimize inbreeding. We have observed that the annual rate of increase of inbreeding in the five Italian cattle under study are lower than what has been reported in several dairy breeds from the United States of America and Ireland, probably because of factors related with reproductive management (high number of breeding bulls, matings programmed to minimize inbreeding, etc.). We have also detected a strong decrease of the effective size that is not accompanied by marked reductions of diversity or substantially increased inbreeding. These N_e estimates should be interpreted with caution due to the inherent difficulty of measuring this complex parameter in a reliable manner.

Supplemental Information