Genetic position of Hungarian Grey among European cattle and identification of breed-specific markers

Abstract

Hungarian Grey is an indigenous cattle breed that is one of the national symbols of Hungary. However, genetic description of the Hungarian Grey cattle has not yet been conducted based on whole-genome screening. Using the GeneSeek high-density Bovine SNP (single nucleotide polymorphism) 150 K BeadChip, we sampled the genome of 36 Hungarian Grey, 12 Maremmana, 13 Hungarian Fleckvieh and 5 Holstein-Friesian cattle for population studies and used data of 139 other cattle from an additional dataset created on European cattle breeds (Upadhyay et al.2017. Heredity 118, 169–176). The performance of a multidimensional scaling plot showed that Hungarian Grey clustered independently from other European cattle. The number and total length of runs of homozygosity (ROH) is similar or slightly below the value of other European cattle; F_ROH coefficients (proportion of the autosomal genome covered by ROH) are similar to Maremmana and Maronesa. The frequency of ROH does not show increased values as it can be noticed in Heck and Maltese. These results indicate that the Hungarian Grey cattle have been successfully maintained avoiding negative genetic effects, and reflect the uniqueness among European cattle. The identification of breed-specific loci has been aimed at differentiating Hungarian Grey (n = 136 in this case) from other cattle breeds (n = 169). Ten loci (−log₁₀P > 5) were identified as markers capable for differentiation of Hungarian Grey. These markers are located on chromosomes 6, 14, 15, 16, 20 and 24.

Implications

Genetic characterisation via whole-genome scan gives feedback to the breeders about the status of this breed. The analysis provides markers to support competitiveness of the Hungarian Grey breed.

Introduction

There are several theories regarding the origins of the Hungarian Grey cattle. One of the theories is that Hungarian Grey arrived in the Pannonian Basin with the 9th-century Hungarian immigration from the east. Another theory leans towards the opinion that domestication of Hungarian Grey took place in the Pannonian Basin, during the reign of House of Árpád Kings, from aurochs (Bos Primigenius), the ancestors of domestic cattle (Bartosiewicz, 1996). The first written document referring to the Hungarian Grey cattle as ‘magnus cornuotes boves Hungaricos’ dates back to the 16th century (Bartosiewicz, 1997). In former times (10th to 13th centuries), most cattle in Hungary had been small, brachyceros-type animals (Bartosiewicz, 1997). Written evidence for westward export (Germany, Italy, Moravia) of Hungarian cattle is known from the 14th century (Miskulin, 1905). The long-horned larger cattle have been introduced to the Balkan countries and the Carpathian Basin by the Turkish Ottoman Empire during the 16th- to 18th-century campaigns (Bartosiewicz, 1996 and 1997).

The breed was used as a draught animal, but it has been bred also for its beef quality. Phenotypic traits of this breed, such as strong pigmentation, long dark eyelashes and well-developed dewlap, show similarities to ancient North African pictures from Egypt and the Sahara region. The Turkish army occupied the Near-East and the North African region long before the European invasion, so it seems that cattle of southern origin have been introduced to our continent later (16 to 18th centuries).

In 1884, the overwhelming majority (78%) of the 4.9 million cattle in Hungary were registered as Hungarian Grey, and half of the 6.7 million still belonged to this breed in 1900 (Mattesz, 1927; Tormay, 1901).

After World War I, when Hungary lost about 72% of its territory, in the remaining part of the former Hungary, by 1925, there were about 321 000 Hungarian Grey cattle (Bartosiewicz, 1997). In the following years, decline of the breed persisted. The World War II severely disrupted the breeding activity, and many herds were destroyed. Later, during the socialist era, as a result of the rescue programme initiated by Imre Bodó and supported by the Hungarian government (1961), ca. 200 purebred cows and 6 bulls were saved (Bodó et al., 1996). The breeding plan used a special rotational mating scheme of disposable bulls to avoid inbreeding (Bodó, 1990). Breeding activity was highly stimulated by the Association of Hungarian Grey Cattle Breeders, which was founded in 1991. Due to small but permanent subsidies by the state, the Hungarian Grey herds started to increase again. In 2011, the Hungarian Grey population amounted to 7000 cows (Bodó, 2011).

As mentioned in a within-breed investigation (Zsolnai et al., 2014), the Hungarian Grey is considered to belong to the ‘Podolic’ cattle group (Negrini et al., 2007), named also Podolian or West Ukrainian; however, the origin of its members has been questioned by several authors (Maróti-Agóts, 2011; Manzone, 2011; Filippini, 2011). The genetic origin of the Hungarian Grey cattle has not yet been definitively elucidated. Mitochondrial DNA studies in Hungarian Grey and Italian Grey cattle showed some similarities to other European breeds; in addition, the T1 haplogroup, located phylogeographically in Africa, was not observed in the Hungarian Grey cattle (Maróti-Agóts, 2011; Lancioni et al., 2016).

A few Maremmana bulls were brought to Hungary at the beginning of the 20th century (Bodó, 2011), and the hereditary 1:29 chromosome translocation was introduced by one of them. The abnormality that is linked to embryo losses is absent in purebred Hungarian Grey animals. This rare chromosome disorder was eradicated from the affected herds by investigating more than 800 cattle and culling the carriers of both sexes (Kovács, 1989; Zsolnai et al., 2014).

The objectives of the present paper were to compare the genetic status of Hungarian Grey to other European cattle in the hope to get insight into the history of the breed, like in the case of Maltese (Lancioni et al., 2016), and to initialise molecular trademark development, already performed in Mangalitza pig (Szántó-Egész et al., 2016; Zsolnai et al., 2017), via looking for differences between Hungarian Grey and European cattle breeds.

Material and methods

Samples

We used the genomic data (150 K SNP (single nucleotide polymorphism), Geneseek Genotyping BeadChip) of blood samples of 136 Hungarian Grey, 12 Maremmana, 13 Hungarian Fleckvieh and 5 Holstein-Friesian samples and 139 animals from 38 different cattle populations examined using a 770 K SNP chip dataset (Illumina BovineHD BeadChip, which contains 777 692 SNPs) (Upadhyay et al., 2017). Altogether 305 animals were analysed (number of Hungarian Greys was 36 in population studies and 136 in the search of breed-specific alleles). Deoxyribonucleic acid was isolated from samples using a simple protocol (Zsolnai et al., 2003).

Animals for this study were selected by the Association of Hungarian Grey Cattle Breeders to represent all the recorded, paternal lines of the herd. Maremmana samples were provided by the University of Veterinary Medicine Budapest; Fleckvieh and Holstein-Friesian samples were used earlier in other research projects (Anton et al., 2018).

For a detailed comparison of Hungarian Grey to other breeds, nine population were selected; Nelore was an outgroup; choices of Boskarin, Busha, Chianina, Heck, Maltese, Maremmana, Maronesa and Romanian Grey were based on their similar phenotypes and their positions relative to Hungarian Grey after principal component analysis.

Data evaluation

Data evaluation followed our previous survey (Bâlteanu et al., 2019). A series of quality control procedures was conducted on raw data using SNP & Variation Suite (SVS) software, v.8.8.1 (Golden Helix, Bozeman, MT, USA): monomorphic markers and unmapped SNPs as well as those with a call rate <95% were eliminated from the dataset. In addition, we removed SNPs with a minor allele frequency <0.05. Duplicated samples (identical by descent value >0.95) and individuals with a genotype call rate <95% were removed. After filtering, the final dataset included 305 animals and 126 150 SNPs. Pairwise genetic distances (F _ST) of the populations were calculated using SVS software v.8.8.1.

The proportions of mixed ancestry and population structure were evaluated with the ADMIXTURE software v.1.3 (Alexander et al., 2009) using default parameters. ADMIXTURE calculates maximum likelihood estimates of individual ancestries based on the data provided by multiple loci (Alexander et al., 2009). We evaluated different number of clusters (K-value from 2 to 10) by considering the mixed ancestry model. The optimal K-value was determined by taking into account the estimates of cross-validation errors (Alexander and Lange, 2011). PLINK software v.1.9 (Purcell et al., 2007) was used to calculate observed (H _o) and expected (H _e) heterozygosities as well as to build a multidimensional scaling plot using a genome-wide identity-by-state pairwise distances matrix (--mds-plot 2 and --cluster options). The --het command of PLINK (Purcell et al., 2007) was used to compute the method-of-moments relatedness F-coefficient.

The detection of ROH (runs of homozygosity) was carried out with PLINK (Purcell et al., 2007). The minimum length of a ROH was set to 1 Mb in order to minimise the detection of spurious ROH. To make sure that the length of ROH is not affected by low SNP density, the minimum number of SNPs that constituted a ROH was set to 50 considering the calculation method proposed by Lencz et al. (2007):

$$l = \frac{{{lo{g_{e\;\frac{{\alpha }\over{{{n_{s\cdot{n_i}}}}}}}}}\over{{lo{g_e}\;\left( {1 - het} \right)}}$$

where n _s is the number of SNPs per individual, n _i is the number of individuals, α is the percentage of false-positive ROH (set to 0.05) and het is the mean SNP heterozygosity across all SNPs. The density of SNPs was set to 1 SNP for each 100 Kb, and a maximum distance of 1000 Kb was allowed between two consecutive SNPs. The scanning window contained 50 SNPs, and the maximum number of missing SNPs per window was set to five with allowance for one heterozygous SNP.

Each ROH was classified based on its physical length into four size categories: 1 to ≤5 Mb, 5 to ≤15 Mb, 15 to ≤30 Mb, and >30 Mb. For each ROH category, the mean sum of ROH per breed was calculated by summing the lengths of all ROH in a given individual for each one of the categories under consideration. The inbreeding coefficient derived from ROH genomic coverage (F _ROH) was calculated by dividing total ROH length per individual by total genome length (2715 Mb) for each individual.

For the identification of loci in Hungarian Grey v. non-Hungarian Grey comparisons, the names of breeds were recoded to 1 and 0, respectively. For the correction of population structure, genomic kinship matrix was used in a multilocus mixed model (Segura et al., 2012). The used model was

$$y = X\beta + Zu + e$$

where y is the phenotypic value, X is the matrix of fixed effects composed of SNPs and covariates (age and farm), Z is the matrix of random animal effects, e means residual effects, and β and u are vectors representing coefficients of fixed and random effects, respectively.

Results

Analysis of diversity and population structure

On the multidimensional scaling plot, our Maremmana, Hungarian Fleckvieh and Holstein-Friesian samples positioned themselves accordingly to the Maremmana, Fleckvieh and Holstein-Friesian groups originating from a public database (Upadhyay et al., 2017). Hungarian Grey keeps a well-separated position from other cattle (Figure 1, Supplementary Figure S1). Descriptive statistics of genetic diversity of Hungarian Grey and the selected nine cattle breeds are shown in Table 1. Observed and expected heterozygosities ranged between 0.353 and 0.502. While Romanian Grey showed high levels of heterozygosity (H _o = 0.502, H _e = 0.373), Busha showed the lowest level of observed heterozygosities (H _o = 0.380, H _e = 0.383). Hungarian Grey, Maronesa and Romanian Grey have the lowest proportion of the genome covered by ROH (F _ROH values are 0.093, 0.092 and 0.061, respectively). The coefficients of pairwise genetic differentiations (Supplementary Figure S2) are ranging from 0.062 (Hungarian Grey–Maremmana) to 0.492 (Nelore–Maltese). If Nelore is taken out from the comparison, the maximal value is 0.295 (Heck–Maltese). Admixture analyses was performed from K = 2 to 10 (Supplementary Figure S3); the most probable K-value of the selected 10 populations was K = 5. Data of Upadhyay et al. (2017) have been completed with all Hungarian Grey animals (Figure 2).

Figure 1 Multidimensional scaling plot depicting the relationships between Hungarian Grey and other cattle populations. The following codes are used: C1 = first component (eigenvalue = 7.681); C2 = second component (eigenvalue = 2.797).

Table 1 Diversity statistics of Hungarian Grey cattle and nine other breeds

F _ROH = proportion of the autosomal genome covered by runs of homozygosity; H _o = observed heterozygosities; H _e = expected heterozygosities; F = method-of-moments relatedness coefficient.

Figure 2 Admixture analysis of cattle populations for a range of K-values (2 to 4). Each individual is represented by a single column divided into K coloured segments, where K is the number of assumed clusters. Populations are separated by white lines.

Figure 3 Number and total length of runs of homozygosity (ROH) in select European cattle populations. The number of ROH estimated in each individual genome (y-axis) is plotted against total ROH size (i.e. number of megabases covered by ROH in each genome, x-axis).

Analysis of runs of homozygosity

We characterised the length, distribution and frequency of ROH in 10 cattle populations (Figure 3). The number of ROH does not show a major departure between the Hungarian Grey cattle and selected breeds with regard to the distribution of individuals according to ROH number and genomic coverage. Hungarian Grey does not show very long ROH segments (>700 Mb) as can be seen in case of Maltese, Heck and Busha individuals. In a comparative analysis of Hungarian Grey (Figure 4), the mean of ROH at medium ROH (5 to 15 Mb) is above the value of Nelore, Maremmana, Maronese, Maltese and Romanian Grey. At 15 to 30 Mb ROH category, Hungarian Grey has the second lowest; at >30 Mb, it has the lowest value. A high level of homozygosity can be observed in Busha, Boskarin and especially Maltese and Heck, reflected also by their F _ROH values (Table 1).

Figure 4 Classification of runs of homozygosity (ROH) identified in select cattle populations based on their size (x-axis) and mean sum of ROH (y-axis, measured in megabases) within each ROH category and averaged per breed.

Identification of loci capable of breed differentiation

During the search for breed-specific loci, we compared the genetic profiles of Hungarian Grey with the available 39 other cattle breeds. We identified 10 candidate breed-specific markers (−log₁₀P > 5) on chromosomes 6, 14, 15, 16, 20 and 24 (Table 2). Using these 10 loci, a multidimensional scaling plot gives two distinct groups (Figure 5): Hungarian Grey and non-Hungarian Grey animals. On chromosome 20, the highest −log₁₀P value (18.96) is far above the value of other SNP hits.

Table 2 List of breed-specific loci and genomic location of Hungarian Grey

FDR = false discovery rate.

Figure 5 Multidimensional scaling plot of Hungarian Grey and European cattle populations. The following codes are used: C1 = first component (eigenvalue = 95.423); C2 = second component (eigenvalue = 14.286).

On chromosome 20, rs43349802, an intronic variant of ATP6V0E1 (ATPase H + Transporting V0 Subunit E1) gene responsible for adiposity, obesity and skeletal muscle development in humans and in mouse, is among the candidate genes associated with mid-test metabolic weight in Hereford and Simmental × Angus cattle (Seabury et al., 2017). ATP6V0E1 is the top positive hubbing gene for oleic acid content in Nelore cattle (Oliveira et al., 2019). Upon investigating transcription during sepsis, ATP6V0E1 was among the differentially expressed genes (Zhang et al., 2017) in three enriched KEGG pathways, such as oxidative phosphorylation, phagosome and epithelial cell signalling in Helicobacter pylori infection.

Discussion

Based on a principal component analysis of mtDNA, among the 18 podolian cattle breeds, two main groups were found (Di Lorenzo et al., 2018), where Maremmana and Chianina belonged to a six-member group, while Slavonian Syrmian, Istrian cattle, Piemontese, Katerini, Calvana, Bianca di Val Padana, Podolsko, Mucco Pisano, Ukrainian, Bulgarian, Hungarian and Romanian Grey were in the other group. However, Hungarian Grey was markedly positioned as a separated breed from the above-listed 11 breeds. In our multidimensional scaling analysis of SNP data, the Hungarian Grey animals also represent a distinct group among the selected cattle (Figure 1) and among the whole dataset (Supplementary Figure S1) of Upadhyay et al. (2017). The three Hungarian Grey animals, outstanding from the main group of that breed, are known as animals with Maremmana ancestry according to their pedigree.

The coefficients of genetic differentiation (F _ST) show that the closest, but moderately differentiated, group to Hungarian Grey are Maremmana (F _ST = 0.086) and Busha (F _ST = 0.067). The latter is known to have Hungarian Grey influence (van Vuure, 2005). The remaining seven breeds display a large genetic differentiation from Hungarian Grey (F _ST > 0.100).

Population structure analyses provide evidence of a composite origin of Busha, Chianina and Maronese cattle (Supplementary Figure S3; 5 ≤ K ≤ 7). At K values >7, Hungarian Grey samples started to split into two sub-clusters, which required further investigation in the herdbook. The three Hungarian Grey animals – displaying Maremmana influence at K-value >3 – have recorded Maremmana ancestry in the herdbook. The Maremmana ancestry of these animals was also detected by an analysis of 18 microsatellites using over 15 000 Hungarian Grey and 32 Maremmana samples where <1% of Hungarian Greys had been proved to have Maremmana ancestry (unpublished data).

On the reproduced figure of Upadhyay et al. (2017), Nelore represents a distinct group (at K = 4). Breeds from the database show similar admixture as reported by Upadhyay et al. (2017), while Hungarian Grey remains in distinct cluster (Figure 2) in agreement with the amplified fragment length polymorphism analysis of Negrini et al., (2007). A mixed origin of Heck, Busha and several other breeds can be noticed, as well as a low portion of mixing between Hungarian Grey and other breeds (excluding Nelore), especially in Maremmana, Busha, Chianina, Romanian Grey and Heck. We noticed reciprocal colours in Hungarian Grey and other European cattle; ocher and yellow portions in Hungarian Grey individuals; and dark blue portions in other European cattle. Such mixing events might have happened via official trading of individuals or during the movement of animal masses in the 16th century when thousands of Hungarian Greys were driven on foot to German fairs and Italy (Strassburg, Augsburg, Nürnberg and Venice) (Miskulin, 1905; Takáts,1927).

An analysis of homozygosity runs revealed that long ROH segments cannot be found in Hungarian Grey (Figure 3), and the mean ROH at 15 to 30 Mb and at >30 Mb was low (Figure 4), indicating a successful avoidance of inbreeding during the implementation of herdbook-based mating plans (Bodó et al., 1996).

SNPs performing well in differentiating Hungarian Greys from other breeds (Figure 5) are the base of developing molecular trademark (Szántó-Egész et al., 2016; Zsolnai et al., 2017) for the beef produced by Hungarian Grey breeders. Similar achievements have been described in cattle such as the inclusion of seven breeds into the quest for breed-specific SNPs (Czech et al., 2018); testing of different identification methods of a small, cost-effective set of SNPs for breed assignment (Kumar et al., 2019); or finding trait-associated loci (Mastrangelo et al., 2016). Our results also demonstrate that Hungarian Grey animals form a distinct group, markedly distinguishable from other European cattle. However, origins of the breed still remain unrevealed; so, further investigation is needed contrasting the samples with Turkish and other eastward-located breeds to follow up the hypothesis of the Turkish migration route of Hungarian Grey.