Evolution of inbreeding: a gaze into five Italian beef cattle breeds history Giacomo Rovelli1 ,2 , Maria Gracia Luigi-Sierra2 , Dailu Guan2 ,3 , Fiorella Sbarra4 , Andrea Quaglia4 , Francesca Maria Sarti1 , Marcel Amills2 ,5 and Emiliano Lasagna1 1 Department of Agricultural, Food and Environmental Sciences (DSA3), University of Perugia, Perugia, Italy 2 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Universitat Autónoma de Barcelona, Bellaterra, Barcelona, Spain 3 Department of Animal Science, University of California, Davis, Davis, CA, United States of America 4 National Association of Italian Beef-Cattle Breeders (ANABIC), San Martino in Colle, Perugia, Italy 5 Departament de Ciència Animal i dels Aliments, Universitat Autónoma de Barcelona, Bellaterra, Barcelona, Spain ABSTRACT In the last decades, intensive selection programs have led to sustained increases of inbreeding in dairy cattle, a feature that might have adverse consequences on the viability and phenotypic performance of their offspring. This study aimed to determine the evolution of inbreeding of five Italian beef cattle breeds (Marchigiana, Chianina, Romagnola, Maremmana, and Podolica) during a period of almost 20 years (2002– 2019). The estimates of Ho, He, Fhat2 , and Fped 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. Moreover, annual rates of increase of the estimated inbreeding coefficients have been very low (Fhat2 = 0.01–0.02%; Fped = 0.003–0.004%). The use of a high number of bulls combined with strategies implemented by the Association of Italian Beef Cattle Breeders ANABIC to minimize inbreeding might explain these results. Despite the fact that diversity and inbreeding have remained quite stable during the last two decades, we have detected a sustained decrease of the population effective size of these five breeds. Such results should be interpreted with caution due to the Submitted 20 May 2021 inherent difficulty of estimating Ne from SNPs data in a reliable manner. Accepted 3 August 2021 Published 27 September 2021 Subjects Agricultural Science, Genetics, Veterinary Medicine, Zoology Corresponding authors Marcel Amills, marcel.amills@uab.cat Keywords Bos taurus, Genetic diversity, Coancestry, Effective population size, Single nucleotide Emiliano Lasagna, polymorphism emiliano.lasagna@unipg.it Academic editor Mudasir Ahmad Syed INTRODUCTION Additional Information and Inbreeding is the main consequence of mating individuals that are related, through Declarations can be found on common ancestry, to a degree that exceeds that of two individuals from the same population page 13 extracted at random (Kardos et al., 2016). Minimizing inbreeding is an aspect of paramount DOI 10.7717/peerj.12049 importance in cattle breeding to avoid the phenotypic expression of detrimental alleles Copyright in the offspring as well as to ensure the maintenance of genetic diversity (Howard et 2021 Rovelli et al. al., 2017). Inbred animals display chromosome segments that are identical-by-descent Distributed under (IBD) and generate long runs of homozygosity (Kardos et al., 2016). Although inbreeding Creative Commons CC-BY 4.0 coefficients can be calculated from pedigrees comprising several generations, in recent OPEN ACCESS years the advent of high throughput genotyping techniques has make it possible to estimate How to cite this article Rovelli G, Luigi-Sierra MG, Guan D, Sbarra F, Quaglia A, Sarti FM, Amills M, Lasagna E. 2021. Evolution of in- breeding: a gaze into five Italian beef cattle breeds history. PeerJ 9:e12049 http://doi.org/10.7717/peerj.12049 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. Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 2/17 Figure 1 Geographical spread of the studied breeds in the different Italian regions. Photo credit: AN- ABIC. Map of Italy adapted from https://www.d-maps.com/carte.php?num_car=14531&lang=it. Full-size DOI: 10.7717/peerj.12049/fig-1 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 K3 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 Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 3/17 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−3 ) 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. Data analysis Calculation of molecular inbreeding coefficient In the five studied breeds, inbreeding coefficient Fhat 2 was calculated for each bull with the PLINK v1.9 software (Chang et al., 2015). The –ibc command of PLINK v1.9 (Chang et al., 2015) was used to compute Fhat 2 . The formula used to calculate Fhat 2 is as follows: Ohom − Ehom Fhat 2 = 1 − Ehom where Ohom is the observed number of homozygotes and Ehom is the expected number of homozygotes. The mean of the inbreeding coefficient was calculated per year (18 levels) for each one of the five studied breeds. Animals born between 1985 and 2002 were merged into a single group because the number of genotyped individuals born before 2002 is very scarce. The PROC REG v14.1 tool (SAS Inst. Inc., Cary, NC) was used to estimate the annual rate of increase in Fhat 2 by fitting a linear regression and considering the 2002–2019 period (Sall, 1981). The R software v4.0.3 (R Core Team, 2018) was used to perform a box plot to represent graphically the inbreeding coefficient variation per year for each of the studied breeds. Calculation of a pedigree –based inbreeding coefficient The Endog software v4.8 (Gutiérrez & Goyache, 2005; Gutiérrez, Goyache & Cervantes, 2009) was used to calculate the pedigree inbreeding coefficient (Fped ), which is defined as the probability that an individual has two IBD alleles (Ferenčaković et al., 2017). In addition to the default variables proposed by the software, we also considered for Fped estimation the average relatedness (AR) coefficient. This parameter is defined as the probability of an allele, chosen randomly from the entire population, to belong to a given animal, so AR can be understood as the representation of the animal in the entire pedigree regardless of the knowledge of such pedigree (Gutiérrez, Goyache & Cervantes, 2009). The depth and completeness of the pedigree are key when estimating inbreeding coefficients, because an incomplete pedigree will lead to an underestimation of the mean inbreeding. We calculated a pedigree completeness index (PCI) for each animal included in the pedigree using the method developed by MacCluer et al. (1983) and implemented in Sargolzaei (2014). The depth of the pedigree varied across breeds, since we have considered Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 4/17 only the generations with all known ancestors. For the rustic breeds (MRM and POD), the inbreeding coefficient (Fped ) was calculated considering four ancestral generations. In contrast, the MAR and CHI breeds were represented by individuals from three ancestral generations while in ROM only genealogical data from two generations were available. As previously said, the mean Fped per year was computed for each of the five studied breeds, merging in one single group the animals born between 1985 and 2002. The annual rate of increase in Fped coefficient was estimated by fitting a linear regression using PROC REG v14.1 (SAS Inst. Inc., Cary, NC) through the time period from 2002 to 2019 (Sall, 1981). This linear regression was plotted with the R software v4.0.3 (R Core Team, 2018). Estimation of genetic diversity and historic effective population size trends The Arlequin software v3.5.2.2 (Excoffier & Lischer, 2010) was used to estimate within- population diversity, by calculating observed (Ho) and expected (He) heterozygosities subsequently corrected over the number of usable SNPs. Historical trends in effective population size (Ne ) were estimated with the SNeP software (Barbato et al., 2015) using default settings and a correction to adjust linkage disequilibrium (LD) r 2 values for small sample sizes. The same index was also calculated through the individual increase in inbreeding, using the software Endog v4.8 (Gutiérrez & Goyache, 2005; Gutiérrez, Goyache & Cervantes, 2009). The formula used to estimate Ne from LD (Corbin et al., 2012), with SNeP software, was:   1  h 1 i − α NT (t ) = (4f (ct )) E r 2 |C adj t Where NT (t ) is the effective population size estimated t generations ago in the past, ct is the recombination rate t generations ago in the past, r 2 adj is the linkage disequilibrium estimation adjusted for sampling bias, and α is a constant. 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, Fhat 2 , and Fped 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 Fped 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 Fhat 2 coefficient, followed by ROM. The analysis of the evolution of Fhat 2 , Fped , 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 Fhat 2 and Fped , we observed Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 5/17 Table 1 Mean, standard deviation, and range (minimum and maximum values) for expected and ob- served heterozygosity. Breed Hoa He b x̄ ± sdc ranged x̄ ± sdc ranged Marchigiana 0.350 ± 0.022 0.349–0.351 0.350 ± 0.019 0.349–0.352 Chianina 0.361 ± 0.024 0.359–0.362 0.356 ± 0.020 0.354–0.359 Romagnola 0.343 ± 0.021 0.340–0.345 0.356 ± 0.020 0.354–0.359 Maremmana 0.383 ± 0.025 0.380–0.384 0.391 ± 0.021 0.388–0.392 Podolica 0.399 ± 0.025 0.397–0.401 0.391 ± 0.021 0.388–0.392 Notes. a Ho, observed heterozygosity b He, expected heterozygosity c x̄ ± sd, mean and standard deviation of Ho and He d range, minimum and maximum values of Ho and He. Table 2 Mean, standard deviation, and range (minimum and maximum values) for different inbreed- ing coefficients. Breed F hat 2 a F ped b x̄ ± sdc ranged x̄ ± sdc ranged Marchigiana 0.012 ± 0.005 −0.051–0.071 0.018 ± 0.002 0.002–0.048 Chianina 0.018 ± 0.005 −0.051–0.023 0.024 ± 0.003 0.002–0.053 Romagnola 0.023 ± 0.010 −0.121–0.072 0.007 ± 0.001 0.000–0.016 Maremmana 0.014 ± 0.003 −0.061–0.044 0.062 ± 0.006 0.039–0.067 Podolica 0.025 ± 0.006 −0.110–0.047 0.061 ± 0.006 0.032–0.068 Notes. a Fhat 2 : molecular inbreeding coefficient. b Fped : pedigree-based inbreeding. c x̄ ± sd, mean and standard deviation of Fhat 2 and Fped . d range: minimum and maximum values of Fhat 2 and Fped . some fluctuations across years that were particularly accentuated for Fhat 2 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 (Fhat 2 and Fped ) are displayed in Figs. 4 and 5. It can be seen that in general Fhat 2 increases slightly but steadily in all five breeds, with averaged overall increasing rates of 0.17–0.34% between 2002–2019. In contrast, Fped remained quite stable across time with averaged overall increasing rates of −0.04–0.08%. The effective population size (Ne ) estimated t generations ago (from 13 to 235) is shown in Fig. 6. It is apparent that Ne decreases markedly across generations. Thirteen generations ago, Ne was lower than 300 for most breeds with the only exception of POD cattle (Ne = 498). In contrast, 235 generations ago Ne oscillated between 1887–3257, which is 7.48 times larger than current values. Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 6/17 Figure 2 Molecular (F hat 2 ) and pedigree (F ped ) inbreeding coefficients and observed heterozygosity of the specialised breeds (Marchigiana, Chianina, and Romagnola) measured in a period of almost 20 years (2002–2019). Ho: observed heterozygosity; F hat 2 : molecular inbreeding coef- ficient; F ped : pedigree-based inbreeding. The dots represent the mean value and the whiskers are the standard deviation. The base generation, rep- resented as 1985–2002, is composed by animals born between 1985 and 2002 (they have been merged in a single group because data from few bulls were available before 2002). Full-size DOI: 10.7717/peerj.12049/fig-2 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 Froh 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 Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 7/17 Figure 3 Molecular (F hat 2 ) and pedigree (F ped ) inbreeding coefficients and observed heterozygosity of the rustic breeds (Maremmana and Podolica) measured in a period of almost 20 years (2002–2019). Ho: observed heterozygosity; F hat 2 : molecular inbreeding coefficient; F ped : pedigree-based inbreeding. The dots represent the mean value and the whiskers are the standard deviation. The base generation, represented as 1985–2002, is composed by animals born between 1985 and 2002 (they have been merged in a single group because data from few bulls were available before 2002). Full-size DOI: 10.7717/peerj.12049/fig-3 averaged (across years) coefficients of inbreeding Fhat 2 and Fped 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 Fhat 2 coefficient is closely related to Fis and can be interpreted (also in the case of Fped computed Fis ) 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 Fped 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 Fped (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 Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 8/17 Figure 4 Overall percentual increase of F hat 2 and F ped for the specialised breeds (Marchigiana, Chian- ina, and Romagnola). F hat 2 : molecular inbreeding coefficient; F ped : pedigree-based inbreeding. The red line represents the regression of the coefficients over time expressed in years. The numbers next to the fig- ure correspond to the overall increase of inbreeding across all analyzed years. Full-size DOI: 10.7717/peerj.12049/fig-4 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 Fhat 2 and Fped 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, Fped 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, Fhat 2 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, Fhat 2 can take negative values (Alemu et al., 2021). Despite these conceptual differences, comparison of Fped and Fhom (similar to Fhat 2 ) coefficients calculated in a pedigree of 245 Holstein cattle with whole-genome sequence data showed a good agreement between both Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 9/17 Figure 5 Overall percentual increase of F hat 2 and F ped for the rustic breeds (Maremmana and Podolica). F hat 2 : molecular inbreeding coefficient; F ped : pedigree-based inbreeding. The red line represents the regression of the coefficients over time expressed in years. The numbers next to the figure correspond to the overall increase of inbreeding across all analyzed years. Full-size DOI: 10.7717/peerj.12049/fig-5 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 (Fhat 2 = 0.014, Fped = 0.061) and ROM (Fhat 2 = 0.023, Fped = 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 (Fhat 2 = 0.012, Fped = 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 Fis 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 Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 10/17 Figure 6 Effective population size (Ne ) across generation (GenAgo) for each Italian beef cattle breed. Full-size DOI: 10.7717/peerj.12049/fig-6 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 Fhat 2 values measured in MAR, ROM, CHI, MRM, and POD breeds with Fhom 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 Fhom 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 Fhat 2 coefficients estimated in the current work. Biscarini et al. (2020) reported an Fped 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 Fis coefficients of 0.106 and 0.138 in the POD and MRM breeds, and Santana Jr et al. (2012) measured an Fped = 0.013 in Brazilian MAR cattle. One potential reason for the discrepancy Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 11/17 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 Fhat 2 and Fped 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 Fhat 2 in the five breeds during 2002–2020, particularly in MRM, POD, and CHI, increasing 0.01–0.02% per year, while Fped 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 Fped remained constant or decreased in Angus, Charolais, and Limousine populations, while in Hereford, Holstein-Friesian, and Simmental a yearly increase of Fped (between 0.06–0.13%) was detected. In another study, Weigel (2001) also reported increases of Fped 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 Ne . Although, Ne 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 Ne 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 Ne . However, it should be also noticed that estimates of Ne 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, Ne 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 Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 12/17 coalescent is that no SNP can be reliably sampled after 4Ne generations in the past. So, our estimates of Ne 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 Ne estimates should be interpreted with caution due to the inherent difficulty of measuring this complex parameter in a reliable manner. ACKNOWLEDGEMENTS The authors wish to thank the ANABIC (S. Martino in Colle, Perugia, Italy) for their cooperation in this study, selecting and providing SNPs data and phenotypic information of Italian beef cattle breeds. The authors also wish to thank the editor and the three referees for their valuable comments on the manuscript and their constructive suggestions. ADDITIONAL INFORMATION AND DECLARATIONS Funding This work was supported by the projects: ‘‘Italian Biodiversity Environment Efficiency Fitness’’—I-BEEF 1 and 2 —2014–2020 and 2020–2023. PSRN: Support for the conservation, use and sustainable development of genetic resources in agriculture (Sub- measure 10.2). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Grant Disclosures The following grant information was disclosed by the authors: Italian Biodiversity Environment Efficiency Fitness’’—I-BEEF 1 and 2 —2014-2020 and 2020-2023. PSRN. Competing Interests The authors declare there are no competing interests. Author Contributions • Giacomo Rovelli performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 13/17 • Maria Gracia Luigi-Sierra and Dailu Guan performed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. • Fiorella Sbarra performed the experiments, authored or reviewed drafts of the paper, and approved the final draft. • Andrea Quaglia and Francesca Maria Sarti analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. • Marcel Amills and Emiliano Lasagna conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft. Animal Ethics The following information was supplied relating to ethical approvals (i.e., approving body and any reference numbers): Blood extraction from young bulls is a routine procedure performed by trained veterinarians from the National Association of Italian Beef Cattle Breeds (ANABIC) for purposes that are completely unrelated to this research project. These blood samples were retrieved well before (2002–2019) our research project began and they were kept by ANABIC. The data were collected in accordance with the FAO guidelines for the characterization of animal genetic resources. Animal management, animal husbandry and trait recording followed the criteria for the assessment of animal welfare as identified and defined by the Welfare Quality Project (Vapnek & Chapman, 2011; Food & Agriculture Organization of the United Nations FAO, 2012). In summary, ethical permission from the University of Perugia was not required. Data Availability The following information was supplied regarding data availability: All the data supporting the results are available in the article or in the Supplemental Files. The raw genotypic and pedigree data available at figshare: Rovelli, Giacomo; Luigi Sierra, Maria Gracia; GUAN, DAILU; Sbarra, Fiorella; Quaglia, Andrea; Sarti, Francesca Maria; et al. (2021): Inbreeding trends in five Italian beef cattle breeds: a perspective of almost two decades (2002-2020). figshare. Dataset. https://doi.org/10.6084/m9.figshare.14625864.v1. Supplemental Information Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.12049#supplemental-information. REFERENCES Alemu SW, Kadri NK, Harland C, Faux P, Charlier C, Caballero A, Druet T. 2021. An evaluation of inbreeding measures using a whole-genome sequenced cattle pedigree. Heredity 126(3):410–423 DOI 10.1038/s41437-020-00383-9. Amaral AJ, Pavão AL, Gama LT. 2020. Genomic tools for the conservation and genetic improvement of a highly fragmented breed—The Ramo Grande cattle from the Azores. Animals 10(6):1089 DOI 10.3390/ani10061089. Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 14/17 Barbato M, Orozco-ter Wengel P, Tapio M, Bruford MW. 2015. SNeP: a tool to estimate trends in recent effective population size trajectories using genome-wide SNP data. Frontiers in Genetics 6:109 DOI 10.3389/fgene.2015.00109. Biscarini F, Mastrangelo S, Catillo G, Senczuk G, Ciampolini R. 2020. Insights into genetic diversity, runs of homozygosity and heterozygosity-rich regions in marem- mana semi-feral cattle using pedigree and genomic data. Animals 10(12):2285 DOI 10.3390/ani10122285. Bjelland DW, Weigel KA, Vukasinovic N, Nkrumah JD. 2013. Evaluation of inbreeding depression in Holstein cattle using whole-genome SNP markers and alternative measures of genomic inbreeding. Journal of Dairy Science 96(7):4697–4706 DOI 10.3168/jds.2012-6435. Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. 2015. Second- generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 4(1):s13742–015 DOI 10.1186/s13742-015-0047-8. Corbin LJ, Liu AYH, Bishop SC, Woolliams JA. 2012. Estimation of historical effective population size using linkage disequilibria with marker data. Journal of Animal Breeding and Genetics 129(4):257–270 DOI 10.1111/j.1439-0388.2012.01003.x. Crow JF, Kimura M. 1970. An introduction to population genetics theory. An introduction to population genetics theory. CABI: Cab direct press. Excoffier L, Lischer HEL. 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular Ecology Resources 10:564–567 DOI 10.1111/j.1755-0998.2010.02847.x. Ferenčaković M, Sölkner J, Kapš M, Curik I. 2017. Genome-wide mapping and estima- tion of inbreeding depression of semen quality traits in a cattle population. Journal of Dairy Science 100(6):4721–4730 DOI 10.3168/jds.2016-12164. Fioretti M, Negrini R, Biffani S, Quaglia A, Valentini A, Nardone A. 2020. Demographic structure and population dynamics of Maremmana cattle local breed after 35 years of traditional selection. Livestock Science 232:103903 DOI 10.1016/j.livsci.2019.103903. Food and Agriculture Organization of the United Nations (FAO). 2012. Phenotypic characterization of animal genetic resources. Rome. Gutiérrez JP, Goyache F. 2005. A note on ENDOG: a computer program for Analysis pedigree information. Journal of Animal Breeding and Genetics 122:172–176 DOI 10.1111/j.1439-0388.2005.00512.x. Gutiérrez JP, Goyache F, Cervantes I. 2009. User’s Guide-ENDOG v 4.6 a computer pro- gram for monitoring genetic variability of populations using pedigree information. Juan Pablo Gutiérrez García press. Available at http:// www.ucm.es/ info/ prodanim/ html/ JP_Web_archivos/ EN_Us_G_.pdf (accessed on 17 April 2021). Howard JT, Pryce JE, Baes C, Maltecca C. 2017. Invited review: Inbreeding in the genomics era: inbreeding, inbreeding depression, and management of genomic variability. Journal of Dairy Science 100(8):6009–6024 DOI 10.3168/jds.2017-12787. Howrigan D, Simonson M, Keller M. 2011. Detecting autozygosity through runs of homozygosity: a comparison of three autozygosity detection algorithms. BMC Genomics 12:460 DOI 10.1186/1471-2164-12-460. Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 15/17 Kardos M, Taylor HR, Ellegren H, Luikart G, Allendorf FW. 2016. Genomics ad- vances the study of inbreeding depression in the wild. Evolutionary Applications 9(10):1205–1218 DOI 10.1111/eva.12414. Khatkar MS, Moser G, Hayes BJ, Raadsma HW. 2012. Strategies and utility of imputed SNP genotypes for genomic analysis in dairy cattle. BMC Genomics 13(1):538 DOI 10.1186/1471-2164-13-538. MacCluer JW, Boyce AJ, Dyke B, Weitkamp LR, Pfenning DW, Parsons CJ. 1983. Inbreeding and pedigree structure in Standardbred horses. Journal of Heredity 74(6):394–399 DOI 10.1093/oxfordjournals.jhered.a109824. Malécot G. 1948. Mathematics of heredity. Mathematics of heredity press. Manunza A, Noce A, Serradilla JM, Goyache F, Martínez A, Capote J, Delgado JV, Jordana J, Muñoz E, Molina A, Landi V, Pons A, Balteanu V, Traoré A, Vidilla M, Sanchez-Rodruguez M, Sànchez A, Cardoso TF, Amills M. 2016. A genome- wide perspective about the diversity and demographic history of seven Spanish goat breeds. Genetic Selection Evolution 48(1):1–9 DOI 10.1186/s12711-016-0229-6. Mastrangelo S, Ciani E, Marsan PA, Bagnato A, Battaglini L, Bozzi R, Carta A, Catillo G, Cassandro M, Casu S, Ciampolini R, Crepaldi P, D’Andrea M, DiGerlando R, Fontanesi L, Longeri M, Macciotta NP, Mantovani R, Marletta D, Matassino D, Mele M, Pagnacco G, Pieramati C, Portolano B, Sarti FM, Tolone M, Pilla F. 2018. Conservation status and historical relatedness of Italian cattle breeds. Genetics Selection Evolution 50(1):1–16 DOI 10.1186/s12711-018-0406-x. Mc Parland S, Kearney JF, Rath M, Berry DP. 2007. Inbreeding trends and pedigree analysis of Irish dairy and beef cattle populations. Journal of Animal Science 85(2):322–331 DOI 10.2527/jas.2006-367. Moioli B, Napolitano F, Catillo G. 2004. Genetic diversity between Piedmontese, Maremmana, and Podolica cattle breeds. Journal of Heredity 95(3):250–256 DOI 10.1093/jhered/esh032. Pryce JE, Haile-Mariam M, Goddard ME, Hayes BJ. 2014. Identification of genomic regions associated with inbreeding depression in Holstein and Jersey dairy cattle. Genetics Selection Evolution 46(1):1–14 DOI 10.1186/s12711-014-0071-7. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, DeBakker PIW, Daly MJ, Sham PC. 2007. PLINK: a tool set for whole-genome association and population-based linkage analyses. The American Journal of Human Genetics 81(3):559–575 DOI 10.1086/519795. R Core Team. 2018. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Available at https:// www.r-project. org (accessed on 17 April 2021). Rovelli G, Ceccobelli S, Perini F, Demir E, Mastrangelo S, Conte G, Abeni F, Marletta D, Ciampolini R, Cassandro M, Bernabucci U, Lasagna E. 2020a. The genetics of phenotypic plasticity in livestock in the era of climate change: a review. Italian Journal of Animal Science 19(1):997–1014 DOI 10.1080/1828051X.2020.1809540. Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 16/17 Rovelli G, Ceccobelli S, Perini F, Sbarra F, Quaglia A, Sarti FM, Lasagna E. 2020b. Genome analysis in five Italian beef cattle breeds. Acta Fytotechnica et Zootechnica 23(4):112–115 DOI 10.15414/afz.2020.23.mi-fpap.112-115. Sall J. 1981. SAS regression applications. SAS Institute, Cary, NC (EUA). Santana Jr ML, Oliveira PS, Eler JP, Gutiérrez JP, Ferraz JBS. 2012. Pedigree analysis and inbreeding depression on growth traits in Brazilian Marchigiana and Bonsmara breeds. Journal of Animal Science 90(1):99–108 DOI 10.2527/jas.2011-4079. Sargolzaei M. 2014. SNP1101 User’s Guide. Version 1. Guelph, Ontario, Canada: HiggsGene Solutions Inc. Sarti FM, Ceccobelli S, Lasagna E, DiLorenzo P, Sbarra F, Pieramati C, Giontella A, Panella F. 2019. Influence a single nucleotide polymorphisms in some candidate genes related to the performance traits in Italian beef cattle breeds. Livestock Science 230:103834 DOI 10.1016/j.livsci.2019.103834. Sbarra F. 2011. Genetics of autochthonous Italian beef cattle breeds. PhD Thesis, University of Padua (Italy) Available at http:// paduaresearch.cab.unipd.it (accessed on 26 April 2021). Vapnek J, Chapman M. 2011. Legislative and regulatory options for animal welfare. Rome: Food and Agriculture Organization of the United Nations (FAO). Visscher PM, Smith D, Hall SJ, Williams JL. 2001. A viable herd of genetically uniform cattle. Nature 409(6818):303–303 DOI 10.1038/35053160. Wang J, Santiago E, Caballero A. 2016. Prediction and estimation of effective population size. Heredity 117(4):193–206 DOI 10.1038/hdy.2016.43. Weigel KA. 2001. Controlling inbreeding in modern breeding programs. Journal of Dairy Science 84:E177–E184 DOI 10.3168/jds.S0022-0302(01)70213-5. Williams JL, Hall SJ, Del Corvo M, Ballingall KT, Colli L, Marsan PA, Biscarini F. 2016. Inbreeding and purging at the genomic level: the Chillingham cattle reveal extensive, non-random SNP heterozygosity. Animal Genetics 47(1):19–27 DOI 10.1111/age.12376. Zorc M, Ogorevc J, Dovc P. 2019. The new bovine reference genome assembly pro- vides new insight into genomic organization of the bovine major histocom- patibility complex. Journal of Central European Agriculture 20(4):1111–1115 DOI 10.5513/JCEA01/20.4.2679. Rovelli et al. (2021), PeerJ, DOI 10.7717/peerj.12049 17/17