Fine-scale analysis of six beef cattle breeds revealed patterns of their genomic diversity

Abstract

There has been a significant decrease in genetic diversity on a global scale for several livestock species, including beef cattle. Therefore, the aim of this study was to deliver fine-scale analysis of the genomic structure and to describe the diversity status of six beef cattle breeds, with a specific focus on genomic inbreeding and the trend of effective population size. The genome-wide data of 639 animals belonging to cosmopolitan (Aberdeen Angus, Charolais, Hereford and Limousin) and local (Piedmontese and Romagnola) cattle breeds were included due to their importance for the local and global beef industry in Europe. In addition, the local dual-purpose Slovak Pinzgau was used to extend the study due to increasing interest from beef cattle breeders in Slovakia. The results showed that more than 3% of the Aberdeen Angus, Hereford and Romagnola genomes, less than 1% of the Limousin, Charolais and Slovak Pinzgau genomes and 0.38% of the Piedmontese genome could have been significantly affected by the mating of relatives over the last three generations. The estimates of effective population size indicated a linear decline within each of the analysed populations. The decrease was mainly evident in populations of the Slovak Pinzgau (7.39 animals per generation), Limousin, Piedmontese and Charolais cattle (average 5 animals per generation). Due to the fact that the decrease in effective population size may be accompanied by a decrease in breed competitiveness with regard to a changing production environment, genetic diversity should be continuously monitored to preserve the long-term potential of the analysed breeds.

Highlights

• Decreased diversity across beef breeds is evident

• Monitoring genetic diversity is needed to prevent inbreeding depression

• Use of genomic data can increase the accuracy of breeding schemes

Keywords:

• Autozygosity
• cattle
• genomic data
• local populations
• sustainable management

Introduction

In recent decades, there has been a significant decline in livestock genetic diversity for several species, including cattle (Kasarda et al. 2017, 2019; Kukučková et al. 2017). This genetic erosion is most likely a consequence of the specialisation of production favouring the world’s most popular and high-producing breeds. It should be stressed that an intense selection of specific sire lines or maternal families continually leads to a reduction in genetic diversity even today. Consequently, monitoring genetic diversity has a crucial role in the sustainable management of breeding programmes, especially local breeds of dairy and beef cattle.

Nowadays, single nucleotide polymorphism (SNP) microarrays are the most popular and cost-effective tool to identify the genetic predisposition of individuals as well as populations with regard to preserving the gene pool of valuable local animal genetic resources throughout the world. Based on the fact that SNP microarrays capture information about thousands of genetic markers distributed throughout the genome, they can provide comprehensive insight into the genomic structure and diversity status of a particular population. Moreover, the rapid advances in genome sequencing and high-throughput DNA techniques have led to new and more accurate measures of genetic diversity, coancestry and genomic inbreeding (Neuditschko et al. 2012; Ferenčaković, Hamzić, et al. 2013; Ferenčaković, Sölkner, et al. 2013; Coates et al. 2018). Most of the methodological approaches for identifying and estimating genetic relationships within and across populations are based on genome-wide average of the estimated number of identity-by-descent (IBD) or alternatively identity-by-state (IBS) fragments among individuals. Particular fragments located in the genome of two or more individuals are IBS if they have identical nucleotide sequence. If such genomic fragments come from a common ancestor without recombination, then they are also IBD because of the common ancestral origin in individuals. Due to factors such as genetic drift, mating of close relatives or selection, IBD alleles in a homozygous state that have accumulated in a particular area of a genome form continuous homozygous stretches with various length described as runs of homozygosity (ROHs) (Mastrangelo et al. 2017). In cattle, a ROH with a specific length can reflect the proportion of shared autozygosity inherited from different generations of common ancestors in the past (Howrigan et al. 2011; Curik et al. 2014). Given this feature, the identification of ROH distributions in the genome is recognised as a key reference method, mainly for calculating genomic inbreeding (F_ROH) (McQuillan et al. 2008).

Genetic distances as the molecular equivalent of the relatedness coefficient are usually expressed as differences between populations in term of the number of mutations or differences in allele frequencies. The genetic distance matrix concerning the identification of paternal lines or maternal families can be depicted simply by coordinate analysis (discriminant analysis of principal components [DAPC] method) or construction of phylogenetic trees resulting from an IBD matrix, Nei’s genetic distances or Wright’s F_ST index (Neuditschko et al. 2012; Chang et al. 2015). Cluster analysis represents another way to determine genetic differentiation within and across populations. It allows determining the membership probability of each individual to particular populations or breeds and identifying the degree of admixture among them at the genomic level. Segments of individual genomes that belong to each cluster are usually estimated by a Bayesian algorithm, which allows applying different model-based clustering approaches (with or without admixture, spatial correlation, a priori information of population origin, etc.) (Toro et al. 2014; Jombart and Collins 2015). Concerning inter-population genetic diversity, genomic-based signatures of recent or ongoing selection and overall breeding history of particular populations, it is also necessary to understand the level of linkage disequilibrium (LD). The extent of LD patterns in the genome may also be used to infer the ancestral and recent effective population size (N_eLD) (Flury et al. 2010).

Current genetic studies are predominantly focussed on cosmopolitan dairy cattle breeds. Nevertheless, the availability of shared genomic data from various resources allows researchers to apply sophisticated methodologies on other breeds that are the main interest for cattle breeders’ society. Such studies significantly augment knowledge about the regional and global state of cattle biodiversity. This information is crucial to conserve specific genomic segments of breeds while simultaneously improving their performance. Aberdeen Angus, Charolais, Hereford, Limousin, Romagnola and Piedmontese breeds were selected due to their importance for the local and global beef industry and comparison of the gene pool among them. In addition, the local dual-purpose Slovak Pinzgau breed was used to extend the study based on the increasing interest of the Slovakian cattle industry in beef production with higher added value in the form of labelled production. Because of the lack of young fattening cattle in Slovakia, extensive import from countries with a long tradition in the production of young cattle for contract fattening is realised. Due to this, ‘new’ breeds are imported into Slovakia (e.g. Piedmontese and Romagnola). The novelty of the study lies in the fact that this is one of the first studies concerning the gene pool of various beef breeds at the genome-wide level in Slovakia as well as in the Central and Eastern Europe.

This study aimed to deliver fine-scale analysis of the genomic structure and to describe the diversity status of six beef cattle breeds, with a specific focus on genomic inbreeding and the trend of effective population size.

Materials and methods

Analysed populations and quality control of SNP data

For this study, the genome-wide data of 639 animals belonging to seven breeds (Aberdeen Angus, Charolais, Hereford, Limousin, Piedmontese, Romagnola, Slovak Pinzgau) were selected. The dataset comprised new and previously published data. The sample size of each breed and genotyping data sources is provided in Table 1. The obtained genome-wide data consisted of representative animals; closely related animals were avoided if possible.

Table 1. Source of genome-wide data and sample size of analysed breeds.

Breed	Region of origin	Sample size	Genotyping Microarray	Data source
Aberdeen Angus	UK	90	BovineSNP50v1 BeadChip	Decker et al. (2009)
Charolais	France	53	BovineSNP50v1 BeadChip	Decker et al. (2009)
		71	International Dairy & Beef Chip	This study
Hereford	UK	98	BovineSNP50v1 BeadChip	Decker et al. (2009)
Limousin	France	100	BovineSNP50v1 BeadChip	Decker et al. (2009)
		17	International Dairy & Beef Chip	This study
Piedmontese	Italy	30	BovineSNP50v1 BeadChip	Decker et al. (2009)
Romagnola	Italy	29	BovineSNP50v1 BeadChip	Decker et al. (2009)
Slovak Pinzgau	Slovakia	151	BovineSNP50v2 BeadChip	This study

The quality control of genotyping data was performed with PLINK 1.9 software (Chang et al. 2015). In first step, the consensus map file was prepared to consistently merge datasets from different sources (Table 1). Overall, 36909 SNP common loci were available before SNP pruning. Only autosomal SNPs with known chromosomal positions according to the latest bovine genome assembly (ARS-UCD1.2) were left for analyses. The minimum call rate for animals and SNP loci was set to 90%. According to previous studies, the level of minor allele frequency (MAF) can bias subsequent analyses (e.g. Qanbari et al. 2010); hence, the MAF threshold was set to 1%. After SNP pruning, the final dataset consisted of 36842 SNP high-informative markers and 634 animals.

Genomic autozygosity and inbreeding coefficient

The observed (H_o) and expected (H_e) heterozygosity and average frequency of minor allele across SNPs was calculated with PLINK 1.9 software (Chang et al. 2015). This software was also used to scan the ROH distribution across the genome of analysed breeds. The ROH stretches were defined as regions with a minimum of 30 SNPs in a run, a minimum length of run 1000 kilobase pairs (kbp), the maximum distance between consecutive SNPs in a window 1000 kbp and a lower density limit of 1 SNP per 100 kbp. The ROH distribution was evaluated separately for five classes, depending on the ROH length (1–2 megabase pair [Mbp], 2–4 Mbp, 4–8 Mbp, 8–16 Mbp and > 16 Mbp). For the class ROH > 16 Mbp, one heterozygous call was allowed. Missing genotypes were enabled for ROH 4–8 Mbp (1 call), ROH 8–16 Mbp (2 calls) and ROH > 16 Mbp (4 calls). The total length of detected ROHs, average ROH length (Mbp) and the sum of all ROHs were calculated separately for each length category and population. Subsequently, the autosome-specific proportion of ROHs longer than 1 Mbp for each breed was determined by the R package detectRUNS (Biscarini et al. 2018).

The genomic inbreeding coefficient (F_ROH) – defined as the proportion of the genome expressed by ROHs with a specific minimum length relative to the total autosomal genome covered by SNP markers (McQuillan et al. 2008) – was calculated according to the formula: $$F_{\mathit{ROH}}=\mathrm{ }\frac{\sum L_{\mathit{ROH}}}{L_{\mathit{autosome}}},$$ where $\sum L_{\mathit{ROH}}$ is the total length of ROHs in the genome of an individual that contains a specific number of continuous homozygous SNPs and $L_{\mathit{autosome}}$ is the length of the autosomal genome expressed by all SNPs. Ferenčaković, Hamzić, et al. (2013) demonstrated that the bovine 50k genotyping microarray is not adequately sensitive to accurately estimate the inbreeding level based on the ROHs 1–4 Mbp. Thus, genomic inbreeding was only calculated for three categories defined by the ROHs with a minimum length of 4, 8 and 16 Mbp (F_ROH4-8Mb, F_ROH8-16Mb and F_ROH>16Mb, respectively).

LD and effective population size

The level of LD across syntenic SNP pairs was expressed by the parameter $r_{LD}^2$ (Hill and Robertson 1968), based on the Pearson’s product-moment correlation coefficient calculated between frequencies of genotypes of two loci using an algorithm implemented in the SNeP programme (Barbato et al. 2015). This approach does not require phasing of haplotypes prior to LD analysis and is much more efficient in case of high-density SNP data or whole-genome sequences. LD values obtained by this approach are identical to the results produced by the analysis of phased haplotype frequencies.

The effective size of a real population X (N_eLD) with known LD values for a certain interval length was defined as the size of a hypothetical ideal population that in equilibrium would have the same level of LD for the same interval length as found in the real finite population (Makina et al. 2015). N_e estimation was based on its functional relationship [E($r_{LD}^2$)≈f(c, N_eLD)], with $r_{LD}^2$ and the recombination rate (c) expressed as the physical position (distance) of SNP markers in the genome (1 centimorgan [cM] ∼ 1 Mbp). This relationship can be presented as follows: $${Ne}_{LD}=\frac{1}{kc}\left[\frac{1}{r_{LD}^2-\frac{1}{n_g}}-\alpha \right],$$ where α is a factor that considers the occurrence of mutations, k is a factor for the inheritance model (e.g. k = 4 in case of biallelic loci), n_g is the sample size of gametes (twice the number of individuals selected from the population) and c is the genetic distance between two loci in Morgans. Assuming constant linear growth of N_eLD with time (t) represented by past generations, the historical effective population size can be expressed as a function of time and physical genetic distance between two loci (Hayes et al. 2003). Longer genetic distances are usually used to estimate the N_eLD in recent generations. The N_eLD estimates for first nine generations are not considered reliable (Corbin et al. 2012); hence, the recent effective population size (N_eLD0) and its confidence interval were calculated by linear regression of N_eLD estimates for 10–60 generations ago (N_eLD10 to N_eLD60).

Genomic structure and differentiation within and across populations

The genetic structure of populations and genetic relationships on an intra- as well as inter-population level were estimated based on calculation of Wright’s F_ST index (Weir and Cockerham 1984), Nei’s genetic distances D_a (Nei 1972), DAPC (Jombart and Ahmed 2011) and a model-based Bayesian clustering approach (Pritchard et al. 2000). Before analysis of genetic differentiation within and across populations, LD-based pruning was done to remove SNP loci that were in LD with each other using PLINK’s pairwise threshold model (indep-pairwise 50 5 0.05). This procedure removed 29580 SNP loci that exceeded an r² threshold of 0.05, after which 7108 remained for subsequent analyses.

Wright’s F_ST index was used mainly to estimate the genetic differentiation degree among analysed populations. Generally, the F_ST index is a measure of population fragmentation expressing the reduction of heterozygosity in subpopulations due to genetic drift. The F_ST index can be defined as the correlation between gametes within subpopulations proportional to randomly selected gametes in the total population: $$F_{ST}=\frac{H_T-H_S}{H_T},$$ where H_T is expected heterozygosity of the total population and H_S is the average expected heterozygosity in subpopulations. The F_ST index calculation was supported by bootstrapping (1000 replications) implemented in the R package StAMPP (Pembleton et al. 2013).

Nei’s genetic distance was defined by Nei (1972) as the difference in genetic constitution between the two populations. This theory assumes that if two populations with a small genetic distance between them are similar, they originated with high probability from common ancestors. Nei’s genetic distances among analysed populations were calculated with the StAMPP package (Pembleton et al. 2013).

DAPC implemented in the Adegenet package (Jombart and Ahmed 2011) was used to determine the degree of genetic differentiation among populations. This multivariate method is usually used to identify and describe clusters that produce genetically related individuals. If the information of the individual’s membership to a certain group or population is unknown, DAPC allows estimating the individual’s origin using a Bayesian-based algorithm. Moreover, this approach allows estimating the probability of an individual’s membership to a particular population and determining the proportion of individual alleles in the population structure. In this study, DAPC was based on the predefined populations according to the animal’s origin. The optimal number of principal components (PCs) reflecting the highest proportion of variance in the database was tested based on the α-score (Jombart and Collins 2015).

The model-based Bayesian clustering approach implemented in Structure 2.3.6 (Pritchard et al. 2000) was used to analyse the systematic difference in allele frequencies among breeds due to different ancestry and to infer the optimal number of clusters represented by all animals in relation to the individual admixture proportion. The default parameters of an admixture model and correlated allele frequencies based on burn-in period of 10000 followed by 100000 Markov chain Monte Carlo (MCMC) iterations were run to perform analysis. The maximum number of clusters was set to 10, where each run was repeated 10 times. The optimal value of K was selected based on the log probability of delta K using the web-based platform Structure Harvester, according to Evanno et al. (2005).

Results and discussion

Analysed populations and quality control of SNP data

After quality control of SNP data, the database included genotypes for 36842 SNP markers covered 2.50 gigabase pairs (Gbp) of the autosomal genome, with an average spacing of 67.95 ± 72.03 kbp. The expected high genotyping call rate at level 99.54% is consistent with previous studies in cattle (e.g. McClure et al. 2018; Sasaki et al. 2018) and reflects the fact that the applied microarray is primarily intended to genotype highly informative SNP loci that are uniformly distributed across the genome of major cattle breed types. It has been shown that a bovine microarray can be partially suitable for genotyping genetically related non-model species, but with a significant decrease in the number of successfully genotyped SNPs (Decker et al. 2009; Kharzinova et al. 2015; Kasarda et al. 2017; More et al. 2019).

The distribution of SNP loci, as well as MAF frequency, was non-uniform across the autosomal genome (Figure 1). The highest number of markers was found on BTA 1 (2452), while the lowest was indicated for autosome 28 (669). Even if the MAF of SNP loci varied across autosomes, the average MAF (0.34 ± 0.13) was relatively homogenous among the analysed breeds and ranged from 0.24 ± 0.15 (Romagnola) to 0.34 ± 0.14 (Charolais). These values are consistent with previous studies in beef cattle (Gautier et al. 2010, Mastrangelo et al. 2018).

Figure 1. Distribution of single nucleotide polymorphisms (SNPs) (left y axis) and average minor allele frequency (MAF) frequency (right y axis) across autosomes.

The observed (H_o) and expected (H_e) heterozygosity were calculated to describe the genetic diversity level and quantify the degree of increase in homozygosity. The observed average H_o (0.33 ± 0.16) and H_e (0.38 ± 0.18) indicated a slight decrease in genetic diversity across the analysed populations. These data were also confirmed with the average Wright’s F_IS index (0.08 ± 0.05). A closer look at the populations showed that the decrease of diversity was more evident in Charolais cattle (F_IS = 0.16 ± 0.02), while the Piedmontese revealed the highest level of genetic variability (F_IS = −0.02 ± 0.02) (Table 2). The lower genome-wide genetic diversity observed for cosmopolitan cattle compared to local ones may be a result of more intense selection pressure on particular genotypic combinations controlling desired phenotypic traits and generally higher genomic uniformity of cosmopolitan breeds (Kukučková et al. 2017). The state of genome-wide genetic diversity in the analysed breeds should be continuously monitored to avoid a reduction in reproductive ability, response to selection and adaptation potential to changing environmental conditions in the future.

Table 2. Genetic diversity indices for each breed.

Breed	Ho ± SD	He ± SD	MAF ± SD	FIS ± SD
Aberdeen Angus	0.316 ± 0.163	0.368 ± 0.200	0.245 ± 0.152	0.022 ± 0.073
Charolais	0.343 ± 0.150	0.412 ± 0.138	0.336 ± 0.143	0.155 ± 0.021
Hereford	0.327 ± 0.157	0.366 ± 0.181	0.257 ± 0.149	0.030 ± 0.094
Limousin	0.340 ± 0.158	0.377 ± 0.141	0.282 ± 0.135	0.077 ± 0.023
Piedmontese	0.348 ± 0.166	0.373 ± 0.177	0.258 ± 0.145	−0.023 ± 0.018
Romagnola	0.317 ± 0.181	0.389 ± 0.224	0.235 ± 0.153	−0.009 ± 0.038
Slovak Pinzgau	0.349 ± 0.151	0.353 ± 0.153	0.261 ± 0.143	−0.011 ± 0.034
Total	0.334 ± 0.162	0.377 ± 0.177	0.268 ± 0.149	0.047 ± 0.079

F_IS: fixation index; H_o: observed heterozygosity; H_e: expected heterozygosity; MAF: minor allele frequency; SD: standard deviation.

Genomic autozygosity and inbreeding coefficient

Distribution of ROH stretches in the genome

ROHs were identified in the genome of all analysed animals in the database. A detailed description of the ROH number and length for each category and breed is given in Table 3. The highest level of homozygosity was found in the Hereford population, where the ROH stretches covered on average 15.84% (average 396.29 Mbp) of the autosomal genome. On the other hand, the lowest proportion of ROHs was identified in the Piedmontese genome, where the homozygous stretches accounted for only 1.48% (average 36.91 Mbp). The proportion of ROHs detected corresponds to the average number of ROHs within breeds. The highest ROH number was found in the Hereford cattle (70.06), whereas the lowest was identified in Piedmontese cattle (7.83). The observed proportion of genomic autozygosity, mainly in the case of Aberdeen Angus and Hereford, corresponds to previous studies in beef cattle (Purfield et al. 2012; Zavarez et al. 2015; Lashmar et al. 2018). However, such a relatively high proportion of ROHs has been found in several dairy breeds (Purfield et al. 2012; Kim et al. 2015; Mastrangelo et al. 2016). In both dairy and beef cattle, the high level of genomic autozygosity is most likely a result of the intensive selection and repeated use of sires or animals from the same paternal line or maternal families in breeding practice (Peripolli et al. 2018). A higher level of genomic autozygosity in Hereford and Angus among cosmopolitan breeds may be the result of historical constitution of particular populations, island origin and thus low genetic fragmentation of founder population (Wright 1940). The lower level of observed genomic autozygosity in Limousin and Charolais cattle is mainly comparable to dual-purpose breeds (Marras et al. 2015; Szmatoła et al. 2016; Kukučková et al. 2017).

Table 3. The length (Mbp) and number (in brackets) of homozygous stretches for each run of homozygosity (ROH) category and breed.

Breed	ROH category	Average ± SD	Range	Proportion in genome (%)
Aberdeen Angus	1–2 Mbp	361.98 ± 148.12 (63.56 ± 13.83)	98.74–845.62 (23.00–98.00)	14.47
	2–4 Mbp	351.60 ± 148.04 (57.50 ± 13.51)	93.40–835.23 (20.00–90.00)	14.06
	4–8 Mbp	278.89 ± 145.38 (31.07 ± 11.62)	76.79–770.69 (12.00–58.00)	11.15
	8–16 Mbp	168.42 ± 121.08 (11.34 ± 6.77)	8.39–624.64 (1.00–31.00)	6.73
	> 16 Mbp	86.58 ± 93.27 (3.29 ± 3.28)	0.00–456.86 (0.00–15.00)	3.46
Charolais	1–2 Mbp	93.04 ± 48.48 (21.08 ± 7.06)	14.07–233.55 (7.00–48.00)	3.72
	2–4 Mbp	87.29 ± 48.58 (17.74 ± 6.51)	7.32–231.72 (3.00–44.00)	3.49
	4–8 Mbp	70.57 ± 44.41 (8.41 ± 3.67)	4.52–221.00 (1.00–20.00)	2.82
	8–16 Mbp	40.47 ± 38.91 (2.79 ± 2.22)	0.00–177.04 (0.00–10.00)	1.62
	> 16 Mbp	21.84 ± 29.45 0.84 ± 1.05)	0.00–141.53 (0.00–5.00)	0.87
Hereford	1–2 Mbp	396.29 ± 201.89 (70.06 ± 22.30)	117.72–1410.24 (24.00–128.00)	15.84
	2–4 Mbp	385.46 ± 201.35 (63.79 ± 2 1.63)	112.82–1404.27 (22.00–124.00)	15.41
	4–8 Mbp	307.34 ± 193.96 34.60 ± 1 5.484)	81.81–328.77 (12.00–90.00)	12.29
	8–16 Mbp	187.20 ± 165.21 (13.12 ± 9.25)	11.45–1133.59 (1.00–57.00)	7.48
	> 16 Mbp	93.29 ± 121.65 (3.67 ± 4.03)	0.00–779.11 (0.00–22.00)	3.73
Limousin	1–2 Mbp	89.67 ± 52.69 (17.47 ± 5.41)	19.88–333.95 (8.00–34.00)	3.59
	2–4 Mbp	84.66 ± 52.73 (14.50 ± 5 .17)	14.00–330.77 (5.00–31.00)	3.38
	4–8 Mbp	66.96 ± 50.11 7.09 ± 3.60)	0.00–308.61 (0.00–20.00)	2.68
	8–16 Mbp	43.20 ± 46.20 (2.87 ± 2.51)	0.00–261.34 (0.00–14.00)	1.73
	> 16 Mbp	22.66 ± 34.84 (0.89 ± 1.23)	0.00–216.55 (0.00–7.00)	0.91
Piedmontese	1–2 Mbp	36.91 ± 30.01 (7.83 ± 4.03)	8.25–158.87 (3.00–20.00)	1.48
	2–4 Mbp	34.66 ± 29.35 (6.55 ± 3.50)	8.25–155.37 (3.00–18.00)	1.39
	4–8 Mbp	25.94 ± 28.43 (2.76 ± 2.53)	0.00–143.30 (0.00–12.00)	1.04
	8–16 Mbp	16.24 ± 22.74 1.07 ± 1.31)	0.00–110.09 (0.00–6.00)	0.65
	> 16 Mbp	9.57 ± 21.01 (0.41 ± 0.87)	0.00–96.33 (0.00–4.00)	0.38
Romagnola	1–2 Mbp	230.73 ± 74.98 (36.17 ± 8.43)	100.05–396.04 (22.00–56.00)	9.22
	2–4 Mbp	225.01 ± 74.29 (32.89 ± 8.09)	96.20–387.61 (21.00–53.00)	8.99
	4–8 Mbp	195.52 ± 69.77 (19.24 ± 5.89)	68.01–362.70 (8.00–34.00)	7.82
	8–16 Mbp	137.84 ± 61.30 (8.90 ± 3.78)	22.46–286.28 (2.00–19.00)	5.51
	> 16 Mbp	79.00 ± 50.67 (3.14 ± 2.05)	0.00–178.42 (0.00–7.00)	3.16
Slovak Pinzgau	1–2 Mbp	73.79 ± 53.67 13.61 ± 7.13)	0.00–347.57 (0.00–37.00)	2.95
	2–4 Mbp	70.22 ± 52.93 11.49 ± 6.43)	0.00–344.42 (0.00–36.00)	2.81
	4–8 Mbp	56.51 ± 49.85 (5.93 ± 3.98)	0.00–337.58 (0.00–19.00)	2.26
	8–16 Mbp	37.32 ± 42.85 (2.46 ± 2.33)	0.00–331.60 (0.00–16.00)	1.49
	>16 Mbp	21.08 ± 34.66 (0.84 ± 1.20)	0.00–263.62 (0.00–7.00)	0.84

Mbp: megabase pair; SD: standard deviation.

Figure 2 illustrates the distributions of ROH stretches (%) with a minimum length of 1 Mbp for each autosome. The results showed that the ROHs distribution, like MAF frequency, is not uniform and depends on the analysed breed. In Aberdeen Angus, Hereford and Romagnola cattle, the largest proportion of ROHs from all stretches detected was found on BTA1 (6.62%, 6.28% and 7.44%, respectively), in Limousin and Piedmontese cattle on BTA12 (9.86% and 9.69%, respectively), in Charolais cattle on BTA6 (6.19%) and in Slovak Pinzgau cattle on BTA26 (9.40%). Given that ROH stretches arise as a result of selection pressure, when only individuals with the desired phenotype are selected for breeding – and thus certain genotype combinations are indirectly selected – the different proportion of ROHs in the genome is most likely a result of different breeding strategies used in grading-up for each breed. In addition to these breed-specific ROH stretches, there were common signals of selection in beef cattle on BTA2, BTA6, BTA7, BTA13 and BTA20 close to the genes responsible for coat colour, muscle development, double muscling, meat tenderness and intramuscular fat content (Moravčíková et al. 2019).

Figure 2. Autosome-specific proportion of run of homozygosity (ROH) stretches with a minimum length 1 megabase pair (Mbp).

Figure 3 shows the differences in the number of ROH stretches (4–8 Mbp, 8–16 Mbp and > 16 Mbp), depending on the length of the genome covered. The results indicate that the number and length of the genome vary considerably, depending on the breed as well as the evaluated animal. Mészáros et al. (2015) reported that such a difference between individuals is most likely a consequence of their distances from common ancestors.

Figure 3. Relationships between the number and length of run of homozygosity (ROH) stretches (4–8 Mbp, 8–16 Mbp and > 16 Mbp). Mbp: megabase pair.

Overall, across all ROH length categories, there were 15008 ROH stretches in the Aberdeen Angus, 18154 stretches in the Hereford, 5010 in the Limousin, 6103 in the Charolais, 540 in the Piedmontese, 5184 in the Slovak Pinzgau and 2910 in the Romagnola genome. One of the most commonly used applications of this analysis is the determination of ancestral or recent inbreeding in a population (Curik et al. 2014). Very long autozygous stretches originate from common ancestors born a few generations ago, while short ROHs are likely derived from more distant ancestors or founders of a given breed. However, short ROH stretches can sometimes persist in the population for a very long time mainly due to a low recombination rate (Curik et al. 2017). Howrigan et al. (2011) also found that the expected length of autozygous stretches follows an exponential distribution, with a mean equal to 1/2 g Morgans (g is the number of generations from a common ancestor). In cattle, Ferenčaković, Sölkner, al. (2013) revealed that ROH > 1 Mbp date back ∼ 50 generations, > 2 Mbp ∼ 25 generations, > 4 Mbp ∼ 12.5 generations, > 8 Mbp ∼ 6 generations and > 16 Mbp reflects ∼ 3 generations of ancestors.

Across the analysed breeds, the total length of ROH mainly comprised short fragments (1–2 Mbp): Aberdeen Angus – 31.76%, Hereford – 37.82%, Limousin – 40.80%, Charolais – 41.44%, Piedmontese – 42.04%, Slovak Pinzgau – 39.64% and Romagnola – 37.59%. The ROH stretches with length 4–8 Mbp accounted for on average 16.94 ± 1.58% from all of the identified stretches for a particular length category. The lowest proportion across all stretches detected occurred for the longest ROH stretches (> 16 Mbp): Aberdeen Angus – 1.64%, Hereford – 1.98%, Limousin – 2.08%, Charolais – 1.65%, Piedmontese – 2.22%, Slovak Pinzgau – 2.45 and Romagnola – 3.13%. The obtained results showed that on average more than 3% of the Aberdeen Angus, Hereford and Romagnola genomes, up to 1% of the Limousin, Charolais and Slovak Pinzgau genomes and 0.38% of the Piedmontese genome could have been significantly affected by the mating of relatives over the last three generations.

Genomic inbreeding

Assuming that the ROH stretches with different lengths represent genomic autozygosity derived from ancestors coming from different generations (Curik et al. 2014), the inbreeding coefficient was calculated for each individual separately based on three ROH length categories (F_ROH4-8Mb, F_ROH8-16Mb, F_ROH>16Mb). Table 4 summarises F_ROH statistics for individual breeds and length categories. On average, ROH stretches with length 4–8 Mbp accounted for 5.72%, while ROHs > 16 Mbp covered 1.91% of the scanned autosomal genomes. The obtained average value of genomic inbreeding showed that the selection was most intense in the population of Aberdeen Angus, Hereford and Romagnola cattle (recent genomic inbreeding at level 3.5%). These results are comparable with previous studies of Purfield et al. (2012) and Peripolli et al. (2017), who found a similarly high intensity of inbreeding in these breeds. This outcome is the result of non-random preferential mating of sires representing competitive lineages and thus higher probability of co-selection of close relatives in subsequent generations.

Table 4. Genomic inbreeding (%) derived from run of homozygosity (ROH) stretches with different length (4–8 Mbp, 8–16 Mbp and > 16 Mbp) for each breed.

Breed	Inbreeding coefficient	Average ± SD	Lower CI	Upper CI	Range
Aberdeen Angus	FROH4-8Mb	11.15 ± 5.81	9.93	12.37	3.07–30.81
	FROH8-16Mb	6.73 ± 4.84	5.72	7.75	0.34–24.97
	FROH>16Mb	3.46 ± 3.73	2.68	4.24	0.00–18.26
Charolais	FROH4-8Mb	2.82 ± 1.78	2.50	3.14	0.18–8.83
	FROH8-16Mb	1.62 ± 1.56	1.34	1.90	0.00–7.08
	FROH>16Mb	0.87 ± 1.18	0.66	1.09	0.00–5.66
Hereford	FROH4-8Mb	12.29 ± 7.75	10.73	13.84	3.27–53.12
	FROH8-16Mb	7.48 ± 6.60	6.16	8.81	0.46–45.31
	FROH>16Mb	3.73 ± 4.86	2.75	4.70	0.00–31.14
Limousin	FROH4-8Mb	2.68 ± 2.00	2.31	3.04	0.00–12.34
	FROH8-16Mb	1.73 ± 1.85	1.39	2.07	0.00–10.45
	FROH>16Mb	0.91 ± 1.39	0.65	1.16	0.00–8.66
Piedmontese	FROH4-8Mb	1.04 ± 1.13	0.61	1.47	0.00–5.73
	FROH8-16Mb	0.65 ± 0.91	0.30	1.00	0.00–4.40
	FROH>16Mb	0.38 ± 0.84	0.06	0.70	0.00–3.85
Romagnola	FROH4-8Mb	7.82 ± 2.79	6.76	8.88	2.72–14.50
	FROH8-16Mb	5.51 ± 2.45	4.58	6.44	0.90–11.44
	FROH>16Mb	3.16 ± 2.03	2.39	3.93	0.00–7.13
Slovak Pinzgau	FROH4-8Mb	2.26 ± 1.99	1.94	2.58	0.00–13.49
	FROH8-16Mb	1.49 ± 1.71	1.21	1.77	0.00–13.26
	FROH>16Mb	0.84 ± 1.39	0.62	1.07	0.00–10.54

SD: standard deviation; CI: confidence interval.

On the other hand, there was a relatively good trend of recent inbreeding (up to 1%) for the Limousin, Charolais and Slovak Pinzgau cattle (Table 4). The lowest level of genomic inbreeding was found in Piedmontese cattle (F_ROH>16Mb = 0.38%). This finding confirms the previous study of Marras et al. (2015), who reported that the majority of the autozygosity of the Piedmontese genome originated from remote ancestors, so the current population is not significantly affected by the increase in inbreeding.

LD and effective population size

Previous studies have indicated that markers with a very low frequency of minor allele may bias LD estimates (Qanbari et al. 2010). Therefore, the LD extent was expressed by the parameter $r_{LD}^2,$ which is less sensitive to differences in allele frequencies (Bohmanova et al. 2010). Considering that $r_{LD}^2$ has limited parametric space and point estimates of effective population size are infinite if $r_{LD}^2$ is equal to 0 or 1 (Uimari and Tapio 2011), only $r_{LD}^2$ values in the range 0.01–0.99 were used to estimate LD and N_eLD.

The density of SNP markers in the database was consistent with previous studies that have assessed LD in cattle based on the 50k SNP microarray (Sargolzaei et al. 2008). However, the average distance between adjacent SNPs was higher compared to other studies (Flury et al. 2010; Beghain et al. 2013) mainly due to the reduced total number of SNP markers after SNP pruning. Descriptive statistics of the observed LD level for each population are shown in Table 5.

Table 5. Average level of linkage disequilibrium in the genome of the studied breeds.

Breed	No. of SNP pairs	${\mathit{r}}_{\mathit{LD}}^2$ ± SD	SNP pairs with ${\mathit{r}}_{\mathit{LD}}^2$>0.3 (%)
Aberdeen Angus	2450579	0.149 ± 0.077	4.23
Charolais	1320646	0.119 ± 0.031	8.71
Hereford	2419523	0.162 ± 0.088	5.06
Limousin	1995006	0.112 ± 0.043	5.52
Piedmontese	2636210	0.103 ± 0.057	1.44
Romagnola	2535148	0.143 ± 0.053	5.68
Slovak Pinzgau	1330045	0.183 ± 0.056	0.63

SD: standard deviation; SNP: single nucleotide polymorphism; r²_LD; linkage disequilibrium.

Aberdeen Angus had the highest average $r_{LD}^2$ value, whereas the Slovak Pinzgau cattle had the lowest observed value (Table 5). The obtained LD values between adjacent SNP markers across evaluated breeds were lower compared to other breeds: for example, Swiss Eringer (0.24; Flury et al. 2010), Holstein (0.2; Bohmanova et al. 2010) or Blonde d´Aquitaine (0.21; Beghain et al. 2013). Although the highest proportion of markers with $r_{LD}^2$ > 0.3 was found in the Charolais genome (8.71%), this value was significantly lower compared to other breeds: for example, Brown Swiss (31.49%) or Nelore (21.51%) (Kasarda et al. 2016). The low level of LD in the genome of the studied breeds may be a consequence of more frequent use of natural services compared to cosmopolitan dairy or dual-purpose breeds. This phenomenon was previously demonstrated in the artificially unselected Italian breed Cinisara (Mastrangelo et al. 2014).

To determine the effect of inter-marker physical distance on LD, all $r_{LD}^2$ values were analysed according to their physical position in the genome depending on the autosome and the breed. All pairs of SNP markers were ordered into 100 groups according to the distance between them, from 50 to 5000 kbp (bin size 50 kbp). Figure 4 shows a rapid decrease in $r_{LD}^2$ values depending on the increasing inter-marker distance. The variability of the LD values within the chromosomes was relatively high, especially in the case of shorter distances (about 250 kb). Arias et al. (2009) found that the recombination distance in the bovine genome decreases directly in proportion to the length of the chromosomes. This phenomenon suggests that the recombination rate increases in proportion to the length of the chromosomes. In other words, longer chromosomes will show lower LD compared with shorter ones (Lee et al. 2011). Many studies have also confirmed a decrease in LD depending on genetic distance (Flury et al. 2010; Qanbari et al. 2010; Kasarda et al. 2016).

Figure 4. The relationship between linkage disequilibrium (LD) and inter-marker distance.

The level of LD between markers located close to each other in the genome is usually used to estimate the trend of ancient N_e (50 kb ∼ 150 generations ago), while LD between markers separated by a much higher distance reflects the current effective population size (4000 kb ∼ 12.5 generations ago). However, estimation of N_eLD more than 60 generations ago based on the inter-marker physical distances (1 Mbp = 0.01 Morgans) may be significantly biased (Flury et al. 2010).

For all evaluated breeds, there was a linear decrease in effective population size over the last 60 generations (average R² of 0.99) (Figure 5). The significant decrease was mainly evident in populations of the Slovak Pinzgau (7.39 animals per generation), Limousin, Piedmontese and Charolais cattle (average 5 animals per generation). In Romagnola, there was a decrease of 2.5 animals per generation. The lowest decrease was detected in the Hereford (1.5 animals per generation) and Aberdeen Angus (1.9 animals per generation) populations.

Figure 5. The trend of effective population size for the analysed breeds over the last 60 generations.

The estimates of current effective population size (Aberdeen Angus, N_e = 52.04; Hereford, N_e = 55.89; Limousin, N_e = 11.51; Charolais, N_e = 5.98; Piedmontese, N_e = 14.06; Romagnola, N_e = 22.97; and Slovak Pinzgau, N_e = 45.18) indicate that the analysed populations are at risk of significant loss of genetic diversity. According to Franklin (1980), populations with less than 50 animals may be considered to be at immediate risk of extinction. Other authors also refer to this value as the threshold for long-term survival of animal populations. For example, Meuwissen (2009) and FAO (2015) stated that a minimum effective population size of 100 individuals is needed to ensure its long-term viability. Specifically, small populations are subject to a much faster loss of genetic variability due to inbreeding or genetic drift compared to populations with high numbers of individuals. In addition, the significant decrease in effective population size as one of the most important genetic diversity parameters may be accompanied by a decrease in breed competitiveness with regard to a changing production environment. Therefore, in order to preserve the long-term genetic potential of the analysed breeds, it is necessary to continuously monitor the trend of effective population size and level of genomic inbreeding.

Genomic structure and differentiation within and across populations

Wright’s F_ST index and Nei’s genetic distances, which are the most commonly used parameters in population genetics to evaluate the relationships on inter-population levels, were calculated to determine the degree of genetic similarity among the analysed breeds. The low average value of genetic distances observed among populations (D_a = 0.07 ± 0.02) indicates a decrease in the genetic differentiation degree among the studied breeds. While particular breeds formed genetically differentiated groups (Figure 6), analysis of genetic differentiation based on the Wright’s statistics confirmed that all evaluated populations of cattle may be considered as separate, i.e. genetically differentiated groups (p < .0001). Simply put, genetically similar groups show less variable F_ST values, while in groups with different genotype frequencies, the F_ST value increase (Barendse et al. 2009). The observed F_ST values ranged from 0.04 to 0.13 (Table 6).

Figure 6. Genetic differentiation within and between evaluated populations based on discriminant analysis of principal components (DAPC) (A), first discriminant function of DAPC (B) and individual membership probabilities suggested by a model-based Bayesian clustering approach (C).

Table 6. Genetic relationships among populations based on the Wright’s statistics (under the diagonal) and Nei’s genetic distances (above the diagonal).

	Aberdeen Angus	Hereford	Limousin	Charolais	Piedmontese	Romagnola	Slovak Pinzgau
Aberdeen Angus		0.080	0.059	0.058	0.066	0.097	0.075
Hereford	0.114		0.063	0.055	0.068	0.098	0.072
Limousin	0.082	0.087		0.034	0.038	0.070	0.069
Charolais	0.079	0.076	0.044		0.045	0.074	0.076
Piedmontese	0.088	0.090	0.045	0.052		0.061	0.053
Romagnola	0.131	0.132	0.090	0.092	0.076		0.096
Slovak Pinzgau	0.102	0.099	0.086	0.073	0.070	0.122

The highest genetic similarity was found between Limousin and Charolais cattle (Table 6), two breeds with a French origin. The origin of the Charolais breed is derived from the autochthonous regional breed of yellow cattle distributed mainly in the Charolles area, surrounded by the Seine, Loire, Rhône and Allier Rivers. Although Charolais cattle were originally horned, in recent years selection has focussed on breeding genetically hornless animals. The Limousin breed originated in the Limousin and Marche regions of southwestern France. Since 1900, Limousin have been bred in the direction of beef cattle, with a high slaughter yield, feed conversion efficiency, growth and muscle development. Limousin muscling is intermediate to that of British cattle breeds such as Aberdeen Angus and Hereford and the extreme double muscling found in the European Belgian Blue and Piedmontese breeds. After Charolais, the Limousin cattle is the second most common breed used for commercial crossbreeding in Slovakia.

The high genetic similarity between the Piedmontese and Charolais populations is most likely explained by the fact that in the 19^th century Charolais bulls were used in the grading up of Piedmontese cattle (Cattle Network – EAAP Working Group 2005). The Piedmontese cattle come from north-western Italy, from the Savoy Alps. The character that has determined the recent genetic evolution of the Piedmontese breed is the development of the so-called ‘double-muscle’ factor, i.e. extremely developed muscular masses not only in the thigh or the buttocks but also in the neck. Due to this factor, Piedmontese cattle have gained the attention of beef cattle breeders around the world, including in North and South America. The Piedmontese breed showed partially higher level of genetic similarity with Slovak Pinzgau cattle compared to the other analysed breeds. The Slovak Pinzgau cattle is one of two national Slovakian cattle breeds and therefore is considered to be an important animal genetic resource. The Slovak Pinzgau belongs to the group of Alpine cattle breeds and is bred mainly for its good longevity, fertility and adaptability on poor production conditions.

The analyses showed that Romagnola and Hereford breeds may be considered the most genetically differentiated populations. In addition, as can be seen in Figure 6(A), Romagnola was the most genetically distant group in this study. The Romagnola breed, which is from the Emilia-Romagna region of Italy, belongs to the Podolic group of grey cattle. The Romagnola have mainly been used as a draft breed, but since the mechanisation of agriculture in the mid-20^th century, it has been bred mainly for beef production.

There was relatively high genetic distance between the Aberdeen Angus and Hereford cattle, both of which originated from the UK. The Hereford is the most widespread beef cattle breed throughout the world. The Hereford breed originated in the county of Herefordshire, in the West Midlands of England. Many Hereford lines have been used for cross-breeding with other breeds to increase economically desired traits. These endeavours have led to changes in the breed’s genome. However, some lines have been kept separate, and these have retained traits of the earlier breed, such as hardiness and thriftiness. The original Hereford is now treated as an animal genetic resource of value for genetic conservation. The Aberdeen Angus breed, a Scottish breed of small beef cattle, was developed in the early part of the 19^th century from the polled and predominantly black cattle of Aberdeenshire and Angus in north-eastern Scotland. Due to easy calving, excellent maternal characteristics, good fertility, longevity and resistance to harsh climatic conditions, it has spread throughout the world.

Relationships among individuals at the intra- and inter-population levels were also evaluated by DAPC. This method is considered to be suitable, especially in the case of a high-density genotyping microarrays, because it allows the analysis of multidimensional data through several dimensions. This method has been successfully applied in various studies focussed on the assessment of population structure in model and non-model animal species (Kasarda et al. 2017; Kukučková et al. 2017; Moravčíková et al. 2017). In the first step of DAPC, 1000 PCs and 10 discriminant functions were used. Subsequently, the test of optimal PC number using α-score (Jombart and Collins 2015) was performed to successfully classify individuals into clusters. The results of this analysis showed that the first 60 PCs, which corresponded to 50.13% of the total variability in the dataset, were sufficient to explain genetic relationships among the analysed breeds. The six discriminant functions described 35.51% of the total variability conserved across animals. The first two principal components (PC1 = 10.86% and PC2 = 7.27%) and the first discriminant function are presented in Figure 6 for visualisation. DAPC demonstrated the distribution of individuals into clusters based on their geographical origin, except Piedmontese cattle (Figure 6(A)). Consistent with expectations, the high degree of genetic similarity was found between the Charolais and Limousin populations with a French origin and Aberdeen Angus and Hereford cattle from the UK, respectively. The Slovak Pinzgau and Romagnola populations formed separate genetic groups. The Piedmontese population was clustered together with Charolais and Limousin cattle (Figure 6(B)) probably as a consequence of cross-breeding between cosmopolitan and local breeds in Italy (Mastrangelo et al. 2018).

The Bayesian clustering and the method proposed by Evanno et al. (2005) – which evaluates the second order of change of the likelihood with respect to ΔK – confirmed the local maxima of K = 7, consistent with previous approaches (Figure 6(C)). The model-based Bayesian clustering approach confirmed that the frequency of alleles varied consistently within each population and individual membership probabilities reflected the origin and unique gene pool of each breed. The membership probability was less than 95% in Hereford, Limousin, Charolais and Piedmontese genomes (Table 7). The analysis revealed the Aberdeen Angus as the breed with the highest contribution to the genetic makeup of Hereford cattle (2.5%). The Limousin cattle gene pool was most affected by Charolais cattle (1.5%) most likely as a result of similar geographical origin. The highest proportion of admixture was found in the Piedmontese genome, where Romagnola cattle was recognised as the breed with the highest contribution to its gene pool (11.6%).

Table 7. The membership proportion of individual breeds.

	Slovak Pinzgau	Aberdeen Angus	Piedmontese	Romagnola	Hereford	Limousin	Charolais
Aberdeen Angus	0.007	0.951	0.005	0.004	0.017	0.006	0.009
Hereford	0.014	0.025	0.013	0.005	0.929	0.006	0.007
Limousin	0.010	0.014	0.012	0.007	0.009	0.933	0.015
Charolais	0.004	0.020	0.008	0.004	0.018	0.035	0.911
Piedmontese	0.019	0.010	0.782	0.116	0.009	0.046	0.018
Romagnola	0.005	0.003	0.031	0.951	0.001	0.005	0.004
Slovak Pinzgau	0.966	0.005	0.012	0.006	0.006	0.003	0.001

Bold values mean the highest membership probabilities identified within particular population/group.

Conclusions

Current genetic studies are predominantly oriented towards cosmopolitan dairy cattle breeds. Nevertheless, the availability of shared genomic data sources allows researcher to apply sophisticated analyses on other breeds that are the main focus of interest for cattle breeders’ society to enlarge knowledge of the regional and global state of cattle biodiversity and to conserve specific genomic segments in breeds while simultaneously improving their performance.

The observed proportion of genome autozygosity corresponded to recent genomic inbreeding, which are expected to be about 3.5% in populations of Aberdeen Angus, Hereford and Romagnola as a result of non-random preferential mating of sires representing competitive lineages and thus higher probability of co-selection of close relatives in subsequent generations. In the Limousin, Charolais and Slovak Pinzgau breeds, there was a good trend of inbreeding (up to 1%) , with the Piedmontese cattle presenting the lower value (0.38%). However, the trend of effective population size across all breeds has declined over time. Over the past 60 generations, there has been a significant decrease in the number of animals. This decline has been most significant in the Limousin, Piedmontese, Charolais and Romagnola cattle, all of which have a current effective population size of less than 25 individuals. For the Hereford, Aberdeen Angus and Slovak Pinzgau breeds, the analysis revealed an effective population size of approximately 50 individuals. One of the reasons for this significant decrease in the effective size of the populations in the current generation is most likely the repeated use of animals from the preferred paternal lines in breeding practice to keep only desired phenotypic traits. It also follows that in order to preserve the gene pool of these breeds for future generations, it is necessary to continually monitor them. Genetic differentiation analysis confirmed that all evaluated populations may be considered genetically unique, i.e. genetically differentiated groups. In the context of gene flow among populations, an observed degree of genetic admixture may be attributed to common geographical origin of particular breeds.

Ethical approval statement

All institutional and international guidelines for the animal hair samples collection were followed.

Disclosure statement

The authors state that they have no conflict of interest.

Additional information

Funding

This work was supported by the Slovak Research and Development Agency [Projects No. APVV-17-0060 and APVV-17-0508], VEGA [1/0742/17] and by the Ministry of Agriculture of the Czech Republic [Project No. QK1910059 and MZE-RO0718-V003].