Skip to main content

Genomic diversity and population structure of the indigenous Greek and Cypriot cattle populations



The indigenous cattle populations from Greece and Cyprus have decreased to small numbers and are currently at risk of extinction due to socio-economic reasons, geographic isolation and crossbreeding with commercial breeds. This study represents the first comprehensive genome-wide analysis of 10 indigenous cattle populations from continental Greece and the Greek islands, and one from Cyprus, and compares them with 104 international breeds using more than 46,000 single nucleotide polymorphisms (SNPs).


We estimated several parameters of genetic diversity (e.g. heterozygosity and allelic diversity) that indicated a severe loss of genetic diversity for the island populations compared to the mainland populations, which is mainly due to the declining size of their population in recent years and subsequent inbreeding. This high inbreeding status also resulted in higher genetic differentiation within the Greek and Cyprus cattle group compared to the remaining geographical breed groups. Supervised and unsupervised cluster analyses revealed that the phylogenetic patterns in the indigenous Greek breeds were consistent with their geographical origin and historical information regarding crosses with breeds of Anatolian or Balkan origin. Cyprus cattle showed a relatively high indicine ancestry. Greek island populations are placed close to the root of the tree as defined by Gir and the outgroup Yak, whereas the mainland breeds share a common historical origin with Buša. Unsupervised clustering and D-statistics analyses provided strong support for Bos indicus introgression in almost all the investigated local cattle breeds along the route from Anatolia up to the southern foothills of the Alps, as well as in most cattle breeds along the Apennine peninsula to the southern foothills of the Alps.


All investigated Cyprus and Greek breeds present complex mosaic genomes as a result of historical and recent admixture events between neighbor and well-separated breeds. While the contribution of some mainland breeds to the genetic diversity pool seems important, some island and fragmented mainland breeds suffer from a severe decline of population size and loss of alleles due to genetic drift. Conservation programs that are a compromise between what is feasible and what is desirable should focus not only on the still highly diverse mainland breeds but also promote and explore the conservation possibilities for island breeds.


Cattle is a species that is present worldwide and well adapted to diverse environments, and has evolved into the most important species for production purposes since its domestication. The presence of cattle was invaluable in the evolution of human societies with a great impact on agricultural development, diet alterations, cultural heritage and social structure [1, 2]. Molecular evidence [3] has indicated two domestication events in cattle from which the taurine (Bos taurus) and indicine (Bos indicus) species arose, sharing the aurochs (Bos primigenius) as common ancestor ~ 250,000 years ago [4]. During the introduction of domesticated cattle into Europe, the Mediterranean coast played a crucial role [2, 5]. Therefore, Greece and Cyprus, which are closely located to the core region of the domestication of cattle, i.e. in the near East, represented an important crossroad for the dispersion of human groups and their herds from Anatolia towards Europe [6,7,8].

Historically, the Southern Balkan peninsula has been characterized by the free movement of humans and animals, especially in the areas near the current borders, from almost the Neolithic throughout Antiquity, the Roman, Byzantine and Ottoman empire times to nearly 40 years before the present [6, 9]. Pastoralism, the seasonal movement of herds and people to exploit the grazing lands across the Balkan region, is a common practice since centuries [10, 11]. Migration events that enhanced the gene flow among the domesticated cattle populations, genetic drift, physical isolation due to geographic barriers [12] during the above historic periods, and human weak selection pressures led to the formation of well-adapted local cattle breeds in rather marginal and harsh environments [6, 13].

The farming and breeding practices associated with these local breeds differ substantially from those in high-performing breeds [14]: (i) they are fairly undifferentiated and unselected; (ii) pedigree records are incomplete or do not exist; (iii) breeding associations either do not exist or are recently established for conservation purposes only; (iv) accordingly, there are no classical breed standards but a breed is defined by common origin; (v) systematic and standardized records of traits do not exist; and (vi) the infrastructure necessary for such recording is rudimentary. Therefore, these local indigenous populations of South East European cattle do not meet all the conditions to be referred as breeds. However, to avoid confusion between the terms breed, strain and population, we will use the term ‘breed’ only with the adjectives ‘indigenous’, ‘local’ or ‘rare’ or without an adjective. These local breeds have been considered to be under constant risk mainly since the 1970s and onwards. In fact, in some extreme cases, the current population size of some local breeds consists of only a few animals. This global trend [15,16,17] reflects also, in generic terms, the situation of Greek local breeds.

In Greece and Cyprus, economic and social conditions as well as geomorphological and climatic reasons have not allowed the development of high-yield local cattle populations since 1950. After the 1960s, with the implementation of artificial insemination, many indigenous populations were crossbred with highly selected commercial breeds [18, 19]. For example, Brachyceros, which is a Greek autochthonous breed with short horns resembling the Albanian and Buša cattle, was crossed with Swiss Brown in some regions [20]. In the middle of the 20th century, eight indigenous cattle breeds were reported in Greece. Nowadays, four of these are considered as officially extinct (Tinos, Andros, Chios, Kerkyra), three as threatened (Brachyceros, Katerini and Sykia) and one (Kea) as a rare breed [21, 22]. The aforementioned breeds were mainly bred in mountain regions and/or on islands with poor infrastructure. The latter combined with geographical and natural barriers as well as the absence of artificial insemination led to reproductive isolation, fragmentation, and gradual depletion of the genetic diversity in these breeds [23].

Genetic diversity studies during the past years have been conducted in a large number of domestic cattle breeds [13, 24,25,26,27,28,29,30]. Improving our knowledge on the genetic diversity within and among local breeds is considered an issue of crucial importance for enhancing their efficient use with regard to sustainable animal farming in a harsh and less intensified environment, and for implementing further conservation programs [31]. Since the indigenous breeds exhibit adaptive ability to their local environment and remarkable longevity, the genetic pool of unselected local breeds can represent a valuable source of genes [32]. However, with a few exceptions [33, 34], little research has been carried out on the genetic diversity and genetic relationships of indigenous cattle from South East Europe with regard to Greece and Cyprus.

Thus, using the single nucleotide polymorphism (SNP) array technology, our objectives were: (i) to obtain unbiased estimations of the neutral genetic diversity of the Greek cattle populations, which represents the first comprehensive genome-wide analysis for these populations; (ii) to evaluate within-breeds’/populations’ different sources of genetic variance as well as their distinctiveness level; (iii) to predict recent admixture patterns of the highly selected and competitive breeds with the non-selected and heterogeneous indigenous breeds from Greece and Cyprus; (iv) to predict historical admixture patterns in cattle breeds in Greece and their expansion route towards the southern foothills of the Alps; and (v) to build an objective basis for the implementation of conservation programs for breeder’s associations, and national and international bodies, as a solution towards the uncontrolled mating and outcrossing of certain rare breeds under high risk of extinction.



Hair roots or blood were sampled from 285 individuals originating from different indigenous cattle populations. Sampling areas are in Fig. 1, which as all subsequent figures, uses the same color code specified for each geographical group in Additional file 1 Table S1. All the samples were collected by trained personnel according to the best veterinary practices. An ethical permission was given by the ethics committee of the Agricultural University of Athens. Briefly, the following local breeds from Greece and Cyprus were sampled in our analysis: (i) from mainland Greece: Greek Brachyceros breed (GRB; n = 97), Katerini breed (KTR; n = 20), Prespa cattle (PRG; n = 10), Rodope cattle (ROG; n = 12), Sykia breed (SYK; n = 16), (ii) from the islands: Kea breed (KEA; n = 97), Agathonisi cattle (AGT; n = 6), Crete cattle (CRT; n = 11), Kastelorizo cattle (KAS; n = 4), Nisyros cattle (NSY; n = 7) and Cyprus cattle (CYP; n = 5). A detailed description of the local cattle breeds studied here is provided in Additional file 2 and see Additional file 1: Table S2. The samples collected in this study were complemented by whole-genome genotypes for GRB (n = 19) and CYP (n = 9) reported by Flori et al. [34] and SYK cattle (n = 5) reported by Verdugo et al. [35] (see Additional file 1: Table S1). In addition, for comparison purposes, we included genetic information of 104 international breeds, based on genetic, historical and geographical criteria. More precisely the large dataset of the above-selected breeds fell into eight main geographic groups (Minor Asia, South East Europe, East Podolian, Tyrrhenian (Apennin-Sicily-Sardinia-Corse), Alpine, France, Iberian and North West European breeds) plus one outgroup that included Gir (GIR), Yak (YAK), and N’Dama (NDA). We included all the above breeds in our study because the geographical origin of some of them have an obvious proximity to Greece and Cyprus or because they have more or less affected the genetic pool of local breeds from Greece and Cyprus through long-term crossbreeding events during the past [36] and possibly till present. For example, the long-term crossing of the indigenous short horned (GRB) cattle and some breeds from the Alpine group (e.g. OBV, BBV and TGV) led to the formation of the KEA breed. However, the GRB and the short horned Buša cattle sampled in the neighboring Balkan countries (South East European group) probably share the same origin, which is why we sampled 18 additional Buša cattle individuals from North Macedonia (see Additional file 1 Table S1). The KTR and SYK breeds are hypothesized to share ancestry with the East Podolian steppe geographic cattle group. The AGT and NSY breeds are hypothesized to share ancestry with some strains of Podolian steppe and Anatolian origin, whereas CYP and KAS might share Anatolian and zebu origins. In addition to the Bos taurus breeds that come from Europe and Minor Asia, in our phylogenetic analyses, we used as outgroups YAK from Mongolia, NDA from West Africa representing African Bos taurus and GIR cattle originating from India but bred in Brazil, representing Bos indicus cattle. These 115 breeds and 10 geographical breed groups used in our analyses are described in Additional file 1: Table S1. We believe that the best way to display the grouping of breeds is by geographical origin (color code) without taking previous genetic knowledge into account, and to improve visualization, we added a symbol for some breeds to indicate possible or known genetic similarities with other breeds or groups. More specifically, these symbols combine the color of the group from which they originate geographically with the color of the group with which they have some genetic similarity. Thus, the breeds BURL (from the Tyrrhenian group), and PUST and PIN (from the Alpine group) display an additional symbol with a color that indicates a known influence from the North West group, the breeds KEA (from Greece), and CABA, AGER, and SBRU (from the Tyrrhenian group) display an additional color that indicates influence from the Alpine group, and some local breeds from Italy (RMG, MCH, CALV, CHI, MARE, and PODO) display an additional color that indicates possible podolian influence [37].

Fig. 1

Origin of the breeds used in the analyzed dataset. Special square marks represent the influence of East-Podolian (grey), Alpine (green) and North-West (olive green) groups

DNA isolation and SNPs

DNA isolation was performed using a commercial kit (QIAamp DNA MiniKit, QIAGEN) according to the manufacturer’s instructions. For SNP genotyping, the Illumina BovineSNP50 BeadChip array was used following standard procedures ( Furthermore, quality controls were applied to obtain a high-quality dataset: (i) only the animals with call rate higher than 0.95 were included; (ii) SNPs that mapped to unknown or sex chromosomes were removed; (iii) SNPs that were genotyped for less than 90% of the samples were removed; (iv) SNPs with a MAF lower than 0.025 and showing departures from Hardy–Weinberg equilibrium within breed (P ≤0.01) were removed. The dataset, at this stage, consisted of 46,678 autosomal SNPs genotyped in 3457 individuals.

Haplotyping and unified additive relationships (UAR)

The Beagle software package (v 5.0) was used for imputation of missing genotypes [38] and haplotype phasing [39]. To improve the efficiency of phasing and imputation, we considered genotyping data of all available bovine animals in the in-house database, which includes a large number of pairs and trios from other projects. Genome-wide relationships between individuals were estimated using a unified additive relationship (UAR) matrix [40] implemented in the R package snpReady [41] and applied to 46,678 SNP genotypes of 3457 animals.

Diversity, phylogeny and population structure analyses require samples of representative and least related animals in each breed. To retain the most representative animals, we excluded outlier individuals (i.e. erroneously sampled and/or admixed animals) by multivariate outlier analysis [42], and then decreased the level of relationship by successive exclusion of highly related animals. Both the multivariate outlier analysis and the reduction of familial structures within breeds rely on the genome-wide additive genetic relationships stored in the UAR matrix. Finally, the dataset used in subsequent diversity and phylogenetic analyses included the 2858 most representative and unrelated animals. The starting and optimized sample sizes for each breed are listed in columns N and Nd in Additional file 1: Table S1.

Haplotype diversity

To design the Illumina Bovine SNP50 BeadChip (Illumina), only five taurine breeds and one indicine breed were considered [43]. Therefore, since this array could contain biased sets of pre-ascertained SNPs, we adopted a 4-SNP-block approach as described previously to reduce this ascertainment bias [13, 14]. Specifically, 4-SNP blocks (haplotypes) that spanned less than 150 kb and had an inter-marker distance shorter than than 50 kb were defined, leading to a compromise between the maximum number of SNPs and the minimum recombination probability within the block [14]. In total, 5756 SNP blocks were taken into account as multi-allelic markers and their haplotypes as alleles in the subsequent unbiased allelic diversity and heterozygosity analyses. Hereafter, SNP blocks are also referred to as multi-allelic markers.

Genetic diversity

The following population genetic parameters were estimated: total number of alleles (nA), mean number of alleles per block (mA), allelic richness (AR), observed (HO) and expected heterozygosity (HE), number of private alleles (npA; alleles observed only in one population), frequency of private alleles (fpA) and number of common alleles (ncA, observed in all subpopulations) as described previously [13, 44, 45]. In addition, we used the number of semi-private alleles (nspA), defined as the alleles observed in two populations only, to further assess the allelic diversity and to consider any genetic and geographical closeness of the included breed pairs. Inbreeding coefficients of each animal i (Fi) were estimated from the diagonal elements of the UAR matrix as Fi= UAR(i,i) − 1. To improve the presentation and discussion of the summary statistics related to diversity, we standardized and then plotted these statistics onto a map with a tessellated projection using the R-script available with the package Tess (

To assess the subpopulation-differentiation, we used the DEST estimator, which is analogous to the classical GST for multi-allelic loci but is unbiased and more suitable when the level of gene diversity is high [46]. In this approach, we used the dataset of the genotypes for 5756 multi-allelic SNP blocks in 115 breeds.

Past effective population size based on linkage disequilibrium

We applied a linkage disequilibrium (LD) based method as implemented in the SneP tool [47] to estimate the historical and recent effective population size (Ne) for all the cattle breeds with sample size larger than 8. The estimation was performed on the SNPs with minimum and maximum distances equal to 20,000 and 10,000,000 bp, respectively, and by applying a recombination rate correction [48] and a sample size correction. The most recent effective population is represented by Ne5, (i.e. five generations ago), the effective population size in pre-industrial times (i.e. 50 generations or 250 years ago) is represented by Ne50, and in times close to domestication (10,000 years ago) by Ne2000. To improve the presentation and discussion of the effective population size across time and space, we standardized and then plotted Ne5, Ne50 and Ne2000 onto a map with a tessellated projection using the R-package Tess as described above.

Phylogeny and population structure

We applied four phylogenetic and population structure analyses to infer relationships between animals and breeds. Two of these analyses rely on bi-allelic SNP genotypes and two on multi-allelic SNP-block genotypes. In addition, two of these analyses represent supervised clustering and two represent unsupervised clustering.

Supervised phylogeny of 115 cattle breeds

To elucidate the phylogenetic relationships among the studied populations, we implemented maximum likelihood (ML) methods using the TreeMix program [49]. In this approach, we used the dataset of SNPs genotyped in 115 breeds (Table 1), and set YAK as outgroup to root the tree. The second supervised approach used the allele frequencies of 5756 haplotype blocks to estimate the Nei’s unbiased DA-distances [50]. Then, the DA-distance matrix of 115 breeds was used to reconstruct the neighbor-net with the SplitsTree4 software and the neighbor-joining tree with the FigTree 1.4 software ( [51]. We set YAK as outgroup to root the neighbor-joining tree.

Table 1 Parameters of genetic diversity in the 115 examined cattle breeds with 46,678 SNPs

Unsupervised population structure analyses

Similar to the above analysis, first we examined the population structure based on SNP genotypes. For this purpose, we conducted clustering with the program Admixture 1.23 [52] under the assumption that the number of clusters is equal to K, with K ranging from 1 to 115, i.e. the number of breeds plus 1. Since the admixture analysis does not need an outgroup, we used the dataset without YAK. To assess the quality of clustering and thus infer the most likely K, we performed 10 cross-validations [53] and estimated the cross-validation error for each K. To illustrate the results of the admixture analyses, we used the Python package Pong [54]. The second unsupervised approach used the proportion of genome-wide shared SNP-block alleles (PS) among all pairs of 2858 animals. The PS matrix was then converted into an allele sharing distance matrix (DPS) as described by Bowcock et al. [55]. Multidimensional scaling (MDS) [56] was used to project the multidimensional DPS distance matrix onto a two-dimensional (2D) plane for the 115 breeds. For a better illustration and interpretation of the MDS results, for each breed, we calculated the mean MDS-coordinates of all the individuals, which correspond to the center of each breed symbol (circle), and then we estimated the standard deviation (SD) around that center. Specifically, for each breed, we calculated the spatial distance of each individual to the group center by applying the Pythagorean theorem, assuming that the hypotenuse is the distance between the center of the breed symbol and the position per individual. Then, we estimated the SD of these spatial distances. Finally, the SD was used as the radius around the breed center symbol, as a proxy for spatial dispersion of animals of each breed. For visualization purposes, plot dimensions were proportionally adjusted in R, considering a 1-inch (= 0.254 cm) length as the longest radius.

We calculated the mean DPS within breed to show the level of breed differentiation. Mean DPS values were standardized and then plotted as a tessellated projection onto a map, using the R-package Tess as described above.

D-statistics analysis

To investigate the historical admixture between taurine and indicine cattle, the D-statistics [57] was computed by using the qpDstats tool of the AdmixTool software package [58]. In brief, for a set of three populations P1, P2 and P3, and an outgroup O that fits to this phylogeny: (((P1,P2),P3).O), the numbers of shared alleles between P1 and P3 (BABA) and, P2 and P3 (ABBA) are calculated by assuming that allele “A” represents the ancestral allele and allele “B” the derived allele. The significant excess of either “ABBA” or “BABA” indicates admixture between populations P2 and P3, or P1 and P3, respectively. We used YAK as outgroup O and Bos indicus (GIR) as P3. We selected Highland cattle that originate from the most northwestern part of Europe (Scotland) as P2 and, thus as a Bos taurus breed with the lowest, if any, admixture level with Bos indicus. All other cattle breeds were tested as P1. To present the gradient of Bos indicus genes in European taurine cattle, we standardized and then plotted the D-values onto a map with a tessellated projection using the R-package Tess as described above. The D-values with Z > |3| were considered as significant and indicated on the map.


Genetic diversity

In total, 590 common alleles were detected among the 115 breeds studied, which represents only 0.7% of the total number (80,720) of alleles. All estimators of genetic diversity for the breeds studied here were highly differentiated among the predefined geographical groups (Table 1). The geographic group of Minor Asia displays the highest average value for almost all the estimators of allelic diversity used, i.e. for: total number of alleles (nA), number of alleles per haplotype block (mA), number of private (npA), number of semi-private alleles (nspA) and allelic richness (AR) (Table 1). The heterozygosity estimates based on multi-allelic SNP blocks and on bi-allelic SNPs are highest for the South East European Buša breeds. The only diversity parameter, which indicates a higher diversity in the central Europe group than in the Minor Asian group, is the observed heterozygosity estimator based on bi-allelic SNPs (HO[SNP]) (Table 1 and Fig. 2). Intermediate values were found for the breeds in the Greek and Cyprus group for mean number of observed alleles (\( \overline{nA} \) = 37,879), mean number of private (\( \overline{npA} \) = 58.6) and semi-private alleles (\( \overline{nspA} \) = 82.6). The Greek and Cyprus group showed the highest mean inbreeding coefficient (\( \overline{F} \) = 0.178) and relatively low values for: mean allelic richness (\( \overline{AR} \) = 3.54), mean observed (\( \overline{{H_{O} }} \) = 0.656) and mean expected heterozygosity (\( \overline{{H_{E} }} \) = 0.641). Within the Greek and Cyprus cattle breeds, CRT exhibited the lowest values for all genetic diversity parameters and, at the same time, the highest frequency of private alleles (fpA = 0.376) as well as the highest average inbreeding coefficient (F = 0.457). For GRB, most of the diversity estimates had the highest values but the frequency of private alleles (fpA = 0.021) and the inbreeding coefficient (F = 0.108) were low. Generally, all the analyzed island populations except the KEA breed, which was sampled on the island of Kea and on mainland, had very high levels of inbreeding and very low diversity parameters (Table 1). GRB and the South East Europe Buša cattle shared similar values for almost all diversity parameters. The tessellated projection of diversity statistics provided strong support for a high allelic diversity in breeds from Anatolia and part of South East Europe. Based on Additional file 3: Figure S1, a southeast to northwest gradient of genetic diversity could be inferred, but it is interrupted by the genetic diversity parameters of the Greek island breeds. However, if the Greek island breeds are excluded, this possible southeast to northwest gradient of genetic diversity remains consistent (see Additional file 3 Figure S1). To illustrate a possible ascertainment bias of the SNP chip data, we present the tessellated projection of the observed heterozygosity estimations based on multi-allelic blocks (HO) and bi-allelic SNPs (HO[SNP]) side-by-side in Fig. 2. This shows that HO[SNP] suggests a high level of genetic diversity in some Alpine and North West European breeds, whereas HO highlights breeds from South East Europe and Anatolia as having the highest level of diversity. The ascertainment bias of SNP chip data was further highlighted by the scatterplot of HO[SNP] versus HO in Fig. 2. Both HO[SNP] and HO are estimators of the true diversity. Therefore, the diversity of the breeds placed above the overall trend line (e.g. North West Europe) is overestimated by HO[SNP] and the diversity of the breeds placed below this line (e.g. Minor Asia) is underestimated by HO[SNP].

Fig. 2

Tessellated projection and value distribution plot of observed heterozygosity estimated based on multi-allelic SNP-blocks (HO) and bi-allelic SNPs (HO[SNP])

Effective population size Ne

The weighted mean of the current effective population size (Ne5) is relatively small and ranges from 28 in the Iberian to 67 in the French geographic group. The Greek and Cyprus group showed a small average value (Ne5 = 40). Only GRB showed a larger Ne5(Ne5 = 85) than the weighted average of the group (Table 1). Going back in the past, the effective population size increased faster in the Buša group than in other European cattle groups. Consequently, during the pre-industrial time (Ne50 250 years ago), the effective population size was clearly larger for the South East Europe Buša (Ne50 = 399) group than for other European breed groups. Only the Minor Asian group showed a larger value (Ne50 = 549) during pre-industrial times. A comparable trend was also observed for effective population size 10,000 years ago (Ne2000). The values for the GRB and Buša groups were comparable and differed from those for the other Greek cattle breeds. The 11 indigenous CRT animals sampled on the island of Crete showed the highest inbreeding level and the smallest effective population size (Ne5, Ne50 and Ne2000) in the entire dataset (Table 1) and see Additional file 4: Figure S2. The tessellated projections provide strong support to the observation that as we go back in time the effective population size increases in regions closer to the domestication center and decreases in North West Europe (see Additional file 4: Figure S2). As Ne is inversely correlated to the extent of LD, our results suggest that the level of LD is high in the fragmented breeds under extinction pressure, which is most probably caused by uncontrolled inbreeding.

Genetic differentiation

Pairwise population differentiation values as estimated by the unbiased estimator developed for multi-allelic markers, DEST, are in Additional file 5: Table S4). Τhe DEST values between 112 cattle breeds (YAK, GIR and NDA excluded) ranged from 0.001 to 0.462. The breeds and breed groups with high allelic diversity and high heterozygosity showed a very low level of differentiation. For example, five breeds of the Minor Asia group showed an average differentiation to each other of only 0.006 and the Buša breeds a value of 0.050. On the other hand, the Greek and Cyprus breeds were highly differentiated (\( \overline{{D_{EST} }} \) = 0.250), followed by the Iberian group (\( \overline{{D_{EST} }} \) = 0.189) and the East Podolian group (\( \overline{{D_{EST} }} \) = 0.180) see Additional file 5 Table S4). Similar trends were also observed when breeds from one pre-defined group were compared to all the other breeds, i.e. the level of differentiation was lowest for the Buša breeds (0.113) and highest for the Greek and Cyprus breeds (0.224).

Within the Greek and Cyprus group, low pairwise DEST values were obtained between two mainland breeds i.e. GRB-SYK (= 0.073) and high values were obtained between two island breeds i.e. CRT-AGT (DEST = 0.413). It is remarkable that the CRT breed showed the highest differentiation level among all the investigated breeds (\( \overline{{D_{EST} }} \) = 0.391). Again, the GRB and Buša breeds shared comparable DEST values.

In addition, we used average allele sharing distance (DPS) between animals within a breed as a measure of the breed-level differentiation. As shown in Additional file 6: Figure S3), the highest DPS values were observed for most of the Buša and Anatolian breeds together with the Greek GRB and the Italian PODO breeds, whereas the lowest values were observed for the Greek island breeds.

Genetic distances and phylogenetic neighbor net

We estimated Nei’s genetic distance DA based on multi-allelic SNP blocks and present the values as a neighbor-net (Fig. 3) and as a neighbor-joining tree routed by YAK (see Additional file 7 Figure S4). In Fig. 3 to improve the visibility of the main part of the tree, we shortened the branch length for YAK. Among the Greek and Cyprus group, the island breeds CYP, AGT and KAS were placed close to Bos indicus. This is also the case for the Anatolian breeds, which form a cluster with the aforementioned breeds. Interestingly, the GRB breed from the mainland Greek group, is positioned within the cluster of the Buša breeds with short branches. On the opposite, for the island breeds CRT, AGT, KAS and NSY, long branches result from a high inbreeding level. KTR and SYK, which represent Greek podolian cattle, are placed between the TRG and Italian podolic breeds. KEA is the only Greek breed that clusters together with some breeds from the Alpine region, as a result of crossbreeding that occurred between indigenous Greek cattle and some Alpine cattle breeds. This cluster includes also some Italian breeds with Brown-Swiss influence (CABA, AGER and SBRU). In addition, we used the allele counts of bi-allelic SNPs and reconstructed the phylogeny with the TreeMix program (Fig. 4) by using YAK to root the tree. Compared to other European cattle breeds, CYP, KAS and AGT were closer to the root of the tree and thus, closer to GIR, the Bos indicus representative. The aforementioned Greek breeds are gradually followed by the cattle breeds from Minor Asia, Greece, Bulgaria, North Macedonia, Kosovo and Serbia. In agreement with the neighbor-net, KEA is the only Greek breed that clusters in the Alpine cluster.

Fig. 3

Neighbor-network based on pairwise Nei’s DA genetic distances among 115 breeds. Special square marks represent the influence of East-Podolian (grey), Alpine (green) and North-West (olive green) groups. Dotted lines indicate the shortened branch length of Yak to improve visibility

Fig. 4

Maximum likelihood (ML) tree inferred from genome-wide allele frequency data by methods implemented in the TREEMIX program. The ML dendrogram of the relationships between the examined cattle populations was rooted with the Mongolian yak as an outgroup breed

The remaining Buša breeds sampled along the Ionian-Adriatic route, i.e. Albania, Montenegro and Dalmatia, were placed after two clusters of Italian breeds. Our phylogenetic analyses, both with multi-allelic (Fig. 3) and bi-allelic markers (Fig. 4), do not aggregate the so-called Podolian or gray steppe cattle breeds (TRG, KTR, SYK, HRI, HRP, and UKP) in a single separate cluster, instead they are scattered along the phylogenetic tree and are positioned closer to their geographic neighbor than to the hypothetical steppe cattle or Podolian group. On the one hand, TRG, KTR, and SYK are positioned between the Anatolian and some of the Greek and North Macedonian breeds. On the other hand, HRI, HRP and UKP do not form an own cluster but are placed among some of the Buša neighbor breeds (Fig. 4). The Italian-podolian breeds form a separate cluster, which is placed between some of the Buša and other Italian breeds and the East podolian breeds. Finally, in the phylogenetic tree reconstructed based on SNP allele counts (TreeMix), for some samples with a high inbreeding level (e.g. CRT and BHB) long branches are observed whereas for less differentiated and highly diverse breeds (e.g. TRG and RMB) short or no branches are found.

Assessment of population structure using unsupervised heuristic and unsupervised model-based methods

We used multi-allelic SNP blocks to estimate the allele sharing distance matrix (DPS) among 2858 animals and projected these by multidimensional scaling (MDS) on a two-dimensional (2D) plane. The MDS projection of 115 breeds is shown in Fig. 5. Αlong the first dimension of MDS (MDS1), we observe that the Anatolian, Greek and Cyprus breeds have an intermediate position between Bos indicus and the remaining European cattle breeds. CYP, AGT and KAS cluster together with the Anatolian breeds, except TRG. Subsequently, the mainland Greek breeds (SYK, KTR, PRG, ROG, and GRB) and the Nisyros island breed (NSY) cluster in the geographic region corresponding to some of the Buša breeds (RHS, MKB, PRE, RMD, SHD, and BHB; South East geographic group), some of the breeds from the East Podolian geographic group (UKP, and HRP) and some Italian breeds (Tyrrhenian group) of South and Central Italy (SINI, MOSA, MOSI, and RSIC) including all Italian podolic breeds (PODO, MARE, RMG, CALV, CHI, and MCH). CRT is isolated from the other two main Greek breed clusters. KEA is positioned in the geographic region of some Alpine and Italian breeds, showing a closer relationship to these breeds than to its own (Greek) geographic cluster. The second dimension of MDS (MDS2) separates the breeds of North Europe and Alpine geographic area from the breeds of central, west and southern Europe.

Fig. 5

MDS projection of the estimated allele sharing distance matrix (DPS) among 115 breeds using multi-allelic SNP blocks. The position of each breed is represented with a circle in which the centre is the average position of all animals of the breed with a radius equal to the SD (standard deviation)

The Admixture analysis presents the second unsupervised clustering method applied to 2832 animals (26 Mongolian yak excluded) using the genotypes of bi-allelic SNPs. The lowest cross-validation error (cv error = 0.5358) was determined at K = 78. Figure 6 presents the clustering at K = 4, 10, 24 and 78. The reasons for choosing K = 78 are that: (i) with a K larger than 4, animals from highly selected breeds start to separate into distinct clusters; (ii) the clusters at K = 10 represent the number of pre-defined geographical breed groups plus 1, i.e. nine groups of European Bos taurus plus NDA and GIR; (iii) the cross-validation error decreases almost linearly from K = 1 to K = 24; and (iv) the clustering with the lowest cross-validation error is at K = 78. Additional file 8: Figure S5 presents the Admixture results with the default color assignment performed by the Pong package for K = 2 till 26. All breeds from Minor Asia and most of the South East geographic groups, as well as GRB, showed a high level of complex admixture even at K = 78. All highly selected breeds and isolated breeds formed their own cluster. GRB and KTR, and SYK to a lesser degree, shared the same pattern with most of the breeds belonging to the South East group with evidence of shared ancestry between these and the Anatolian breeds. The highly inbred and differentiated island breeds assigned to separate clusters. The same was also observed for small-sized mainland populations, such as PRG and ROG. At K = 24, a significant proportion of the Alpine ancestry (OBV, and BBV) was estimated for KEA, which is also retained at the most probable true K = 78.

Fig. 6

Model-based clustering of the estimated membership fractions of individuals from all examined breeds. Each cluster is indicated by a different color and each individual is shown as a vertical bar, at K = 4, K = 10, K = 24 and K = 78. The K-value of 78 represents the lowest cross-validation error (cv = 0.5358). The heights of the colored segments are proportional to genotype memberships. The names of the breed groups are given according to their geographical distributions as shown in Fig. 1. For a full definition of the breed groups (see Additional file 1: Table S1). From left to right within each group are presented the following breeds: Minor Asia: ATER, ATBC, ATSR, ATSY, TRG; Greece and Cyprus: CYP, AGT, CRT, NSY, GRB, KAS, KEA, PRG, ROG, SYK, KTR; South East Europe: RHS, MKB, SRB, PRE, RMB, SHB, DGB, DBB, MAB, LKB, SKB, MNB, BHB, HRB; East Podolian: HRI, HRP, UKP; Tyrrhenian: PODO, CINI, MOSI, RSIC, MOSA, SARD, SBRU, CORS, AGER, MARE, CHI, MPIS, CALV, MCH, RMG, GARF, PONT, MODE, CABA, REGG, PMT, BURL; Alpine: PRDO, OVAR, REND, BPUS, PUST, SIC, PIN, TGV, MWF, OBV, BBV, DFV, FGV, VOG, ABO, MON, TAR; France: RDBI, SAL, AUB, LIM, CHR, PAR, BAQ, GAS; Iberian: MNRQ, MALL, NGAN, CANA, MARI, ALEN, BAR, MARO, SYG; North West Europe: BPN, NOR, MAN, BBB, LKF, HF, GNS, JSY, HER, SHR, KRY, DXT, GLW, AAN, HGL, NRC, SERC, FJL, FIAY, FINE, FINW, FINN, YARO

D-statistics analysis

The historical admixture between taurine and indicine cattle was confirmed for most of the cattle breeds that originate from Anatolia, Cyprus, Greek and South East Europe and most of the Central and South Tyrrhenian breeds. The KEA, AGER and SBRU breeds, which are influenced by the Alpine group, as well as the MPIS and SARD breeds deviate from the above described general geographic trend. In Fig. 7, the values of all these significant D-statistics (Z < |3|) are shown with a red spot on a tessellated map to highlight the gradient of the Bos indicus introgression from the Southeast to the Northwest direction. The results provide strong evidence that support the influence of Bos indicus along the Balkan continental route and along the Mediterranean route up to North Italy.

Fig. 7

Tessellated projection of D-statistics values. Red dots represent a significant influence of Bos indicus (Z > |3|)


In this study, we focused on the genome-wide structure of 11 indigenous cattle populations, 10 from Greece and one from Cyprus. The dataset included all the indigenous cattle populations that are reared in Greece and Cyprus and considered to be under constant risk mainly since the 1970s and onwards. An essential reason for this situation is the increased use of sires from the so-called cosmopolitan breeds mostly from the mainland. Sampling included as many individuals as possible from the whole Greek territory (Fig. 1). Some island indigenous breeds (CYP, AGT, KAS, and NSY) were represented by very small sample numbers (Table 1). However, due to the aforementioned limited population sizes and limitations regarding the herds (i.e. highly related, almost feral and difficult to handle), it was not realistic to consider that the number of samples can be increased. For the remaining Greek breeds, a sufficient number of 10 or more of the most representative individuals of the local breeds were sampled. Therefore, this study represents the most complete and up-to-date dataset available for the indigenous Greek cattle and integrated into the most complete Minor Asia and European sample collection.

Improving our knowledge of the genetic diversity within and among local breeds is an issue of crucial importance for implementing further conservation programs that are necessary for sustainable development and future animal farming in changing environments [34]. This could be particularly important for traditional unselected breeds that cover a geographical area close to the domestication center [32].

In spite of the ascertainment bias of the BovineSNP50 chip data [13, 43], it was highly informative for the analyzed Greek and Cyprus cattle populations as reflected by the consistent values obtained for various parameters of genetic diversity levels. These parameters indicated a more profound loss of genetic diversity as well as a higher risk of inbreeding in the Greek island populations than the Greek mainland breeds. Experimental studies [59] have demonstrated that a higher level of nucleotide diversity is associated with a stronger selection response under stressful conditions. Furthermore, Vilas et al. [60] reported that the high adaptive potential of a population is better indicated by the allelic diversity of neutral markers than by expected heterozygosity. Here, we investigated cattle breeds, whose geographical origin ranged from the cattle domestication center to the most western and the most northern parts of Europe. The overall diversity estimators suggest higher allelic diversity in breeds from the region of South East Europe up to Anatolia. However, this general trend is interrupted by some highly fragmented and inbred Greek cattle breeds (see Additional file 3: Figure S1), especially by island breeds such as CRT, AGT, NSY and KAS which are represented by samples of the last remaining indigenous animals. In complete accordance with the strong effect of genetic drift, the corresponding average frequency of private alleles (fpA; Table 1) reached the highest value in the Greek island breeds. Comparable high fpA-values were observed only in genetically isolated island breeds such as MNRQ and MALL (Table 1). We observed a general trend that characterizes the fragmented breeds at risk of extinction, i.e. that they are highly inbred with a small number of high-frequency private alleles. This is reflected in the high positive correlation between F and fpA (r(F, fpA) = 0.72) and low negative correlation between F and npA (r(F, npA) = − 0.21).

As outlined in the Background section, all indigenous breeds of South East Europe, Greece, Cyprus, and Anatolia are either under weak artificial selection or not under any kind of coordinated artificial selection. However, there are substantial differences with respect to isolation. The estimated DEST levels (see Additional file 5: Table S4) suggest a substantial genetic differentiation of the Greek and Cyprus group compared to the remaining groups. This high differentiation is attributed to the genetic drift and inbreeding in the highly fragmented and naturally isolated island breeds (CRT, AGT, NSY, and KAS) and is confirmed by the low differentiation between animals within these breeds (see Additional file 6: Figure S3).

However, empirical studies [59] provide evidence that, rather than inbreeding, genetic drift led to a reduction in diversity in comparable scenarios. Moreover, purging in stressful environments (i.e. natural selection) can maintain a higher diversity level than expected by inbreeding.

In general, each island of the Aegean represents a different entity with a variety of environmental and socioeconomic factors that affect the livestock populations [61]. However, the cattle populations from the Greek islands have constantly suffered from declining population sizes during recent decades. In the case of the CRT breed, even 60 years ago, the statistical data for livestock numbers of indigenous cattle in Crete identified only 149 of the 1954 (7.6%) recorded animals as indigenous or local, which is a rather small percentage [62]. This statistic also documents a very small number of animals for such a large island (one cow per 4.23 km2). For the whole territory of the Dodecanese islands, on which the AGT, NSY, and KAS breeds are raised, the same source identifed 2434 of the 6210 recorded animals (39%) as local. Today, the population size of all the above-mentioned cattle breeds is very small (see Additional file 1: Table S3 and Additional file 2). In addition to small population sizes, genetic drift is strengthened by geographical isolation that further hampers animal exchange. Thus, in such cases, inbreeding and genetic drift are unavoidable factors that shape the diversity of these populations. More generally, there is an opposite trend between the economic significance of cattle and of small ruminants when the regions along the line from the Aegean islands toward the Greek mainland, Balkan and Middle Europe are compared.

In contrast to the island breeds, the lowest levels of genetic differentiation (DEST) were observed for the breed groups from Minor Asia and South East Europe and for some mainland Greek breeds [e.g. GRB, PRG, SYK, and KTR (see Additional file 5: Table S4)]. This low differentiation is accompanied by a high level of allelic diversity including a large number of low-frequency private and semiprivate alleles. Low differentiation and high diversity are probably the consequences of a low artificial selection pressure combined with low genetic drift in effectively large (Table 1) and less isolated populations.

All supervised clustering methods applied in this study clearly reveal the level of genetic differentiation. The highly differentiated breeds show long branches in the neighbor network (Fig. 3), neighbor joining tree (see Additional file 7 Figure S4) and maximum likelihood tree (Fig. 4). This is independent of the evolutionary sources of the differentiation, i.e. the artificially selected and artificially isolated breeds (e.g. beef breeds from the UK) and the naturally isolated and naturally selected breeds (e.g. Greek island breeds) show very long branches. However, breeds that are under strong random sampling (island breeds) show even longer branches especially in the ML analysis. On the contrary, for local breeds with a low differentiation and high diversity level, the three methods reveal short branches or even no branches (Figs. 3 and 4) and see Additional file 7: Figure S4. The unsupervised clustering methods provide a partly different tale. Most of the animals from the artificially unselected breeds of South East Europe, Greece, Cyprus and Anatolia remain unclustered with the admixture approach, even at very high K-values (Fig. 6) and see Additional file 8 Figure S5. The CRT breed which has the highest inbreeding level (0.457) among the 114 breeds investigated here was the first unselected breed that was clearly clustered with the admixture approach at K = 7. The K-value of 78 is associated with the smallest cross-validation error and should represent the true number of clusters in our design. Among the aforementioned unselected breeds, we observed clear clustering for the island breeds and the most inbred Buša breed BHB only, which are characterized by high inbreeding. We observed a similar trend regarding clustering for the Greek KTR, SYK, ROG and PRG populations but in different proportions, which reflects their isolation status (see Additional file 2 and (Fig. 1). The remaining unselected breeds and/or animals remained as an unresolved mixture even at the K-value with the lowest cross-validation error.

The MDS projection places almost all the unselected breeds from Minor Asia, Greece and Cyprus, and geographic breed groups from South East Europe and East Podolian, in the overlapping space with suggestive trends (Fig. 5), i.e. Anatolian and Cyprus breeds and the Greek island populations KAS and AGT close to the root of the tree and to GIR. The remaining breeds that were sampled in Greece were placed among the above-mentioned island breeds, the South East European Buša and some Italian local breeds. Thereby the breeds that are geographic neighbor of each other overlap in the MDS presentation. Only the KEA breed, which is the historical product of a cross between indigenous animals similar to the modern-day breeds of Brachyceros and some Alpine breeds (OBV, BBV and TGV), is placed between the Alpine and Tyrrhenian breeds, which is also confirmed by the admixture results (see K = 8 to 26 in Figure S5 [see Additional file 8 Figure S5].

Estimating the LD-based effective population size (see Additional file 4: Figure S2) for different time intervals during the evolution of cattle provides an interesting insight into their demographic history. For example, the pattern of Ne at the time of domestication (~ 2000 generations ago) suggests that compared to western European breeds, modern South East European cattle breeds had larger founder population sizes as a result of their proximity to the center of domestication. In fact, up to ~ 50 generations ago, cattle breeds from southeastern Europe had larger Ne than those from northwestern Europe. However, the current Ne indicates that the population size of the South East European cattle breeds such as Buša and GRB, has decreased substantially in the last 50 generations or so. Furthermore, this observation implies that such populations have accumulated a small number of common but long haplotypes, which contribute to the high LD. The industrial revolution, modern breeding practices and changes in customers’ preference have led to the replacement of local cattle breeds with commercial cattle breeds, which mostly originate from northwestern Europe. This factor, coupled with uncontrolled breeding in the remaining fragmented populations due to the absence of modern reproduction and breeding management practices contributed to population shrinkage of most of the previously highly diverse cattle breeds from southeastern Europe.

Based on Fig. 5 and Additional file 8: Figure S5, and based on the results of the previous studies (e.g. [34, 63]), we assume that introgression of indicine ancestry into Anatolian and some Mediterranean cattle occurred. To test this putative indicine introgression, which could have occurred along the migration route from Anatolia, Greece and South East Europe to the southern foothills of the Alps and North West Europe, we calculated the D-statistics. For all the breeds from these areas (Minor Asia, Greece and Cyprus, South East Europe and East Podolian) except KEA and also all the Italian Podolian breeds and some other breeds from South and central Italy, we obtained significant D-statistic values (Fig. 7). The first breed from the southern Alps for which no significant Bos indicus influence was found is SIC (Z < |3|), with a D-value close to zero. The D-values clearly decrease as the spatial distance to the origin of Bos indicus increases. Interestingly, three breeds from the East Podolian roup also show a significant indicine introgression, among which UKP that was sampled east of the Carpate Mountains. As discussed recently by Verdugo et al. [35], indicine introgression started ~ 4000 years ago and may have been stimulated by the onset of a period of increased aridity known as the 4.2-thousand-year abrupt climate change event. The increased level of diversity in the breeds from South East Europe, and in the Anatolian and some Greek cattle breeds could be, also, the result of various demographic events, including the introgression of Bos indicus alleles. It should be noted that we cannot distinguish the proportion of diversity caused by introgression from that caused by other evolutionary forces, e.g. low pressure of artificial selection. We also note that some estimators of allelic diversity, such as the number of private and semi-private alleles measures the proportion of unique alleles not present even in GIR or other neighbor breeds. The fact that we sampled and genotyped some of the close neighboring breeds intrinsically reduced the number of private and semi-private alleles. However, although GIR and many neighbor breeds from South East Europe and Greece were sampled (Fig. 1), the numbers of private and semi-private alleles remained larger than for breeds from other parts of Europe.

Conservation of genetic diversity for sustainable development must be understood as a global long-term task. The repeatedly confirmed evolutionary trade-off hypothesis [59, 64] suggests that increased fitness in the environment of selection is accompanied by a decrease in fitness in other environments. This trade-off is applied for high-performance breeds that are adapted to benign (temperate) environments and for breeds that are adapted for production in stressful environments. In spite of the above-mentioned trade-off, replacement of well-adapted local breeds by high-performance breeds and/or animals from their crosses is forced by the free-market mechanism that intrinsically causes the extinction of local breeds [65]. This replacement process has either already reached the terminal stage or is close to reaching it in many regions [14]. If we consider animal diversity as a globally shared-resource and as a prerequisite for sustainable development in the changing environments, then gradual depletion of neutral diversity present in local breeds under soft selection pressure is a kind of “tragedy of the commons” [65, 66]. Therefore, conservation breeding programs must be seen as a regulated long-term exploitation of common resources.

As already discussed above, most probably, because of an abrupt climate change event, farmers started to cross taurine cattle with Bos indicus [4, 9] during the early Bronze Age. Such crosses increased the already high level of diversity of local soft-selected breeds. The subsequent long-term adaptation to local environments shaped the bovine mosaic genomes with an unknown but low proportion of indicine alleles in modern-day cattle breeds from Anatolia up to the southern foothills of the Alps. Therefore, conservation of a high diversity level in these partly fragmented breeds could offer valuable genetic resources for future human needs and future abrupt or gradual climate-change events.


The phylogenetic patterns derived from genome-wide information were quite consistent with the geographic positions and the historical information regarding crossbreeding between Greek and Cyprus cattle. All investigated Cyprus and Greek breeds present a complex mosaic genome of historical and recent admixture events between neighbor and well-separated breeds. While the contribution of some mainland breeds to the pool of genetic alleles seems important, some island and fragmented mainland cattle strains suffer from a severe decline in population size and a loss of alleles due to strong bottlenecks and genetic drift. However, in spite of a markedly reduced genetic diversity level in most island breeds, these show high fertility and longevity in stressful environments. Conservation programs that are a compromise between what is feasible and what is desirable should focus not only on the still highly diverse mainland breeds but also promote and explore the conservation possibilities for island breeds.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.


  1. 1.

    Reed CA. Animal domestication in the prehistoric near east. Science. 1959;130:1629–39.

    CAS  PubMed  Google Scholar 

  2. 2.

    Howe T. Domestication and breeding of livestock. Horses, mules, asses, cattle, sheep, goats and swine. In: Campbell GL, editor. The Oxford handbook of animals in classical thought and life. Oxford: Oxford University Press; 2014.

    Google Scholar 

  3. 3.

    Loftus RT, MacHugh DE, Bradley DG, Sharp PM, Cunningham P. Evidence for two independent domestications of cattle. Proc Natl Acad Sci USA. 1994;91:2757–61.

    CAS  PubMed  Google Scholar 

  4. 4.

    Pitt O, Sevane N, Nicolazzi EL, MacHugh DE, Park SDE, Colli L, et al. Domestication of cattle: two or three events? Evol Appl. 2018;12:123–36.

    PubMed  PubMed Central  Google Scholar 

  5. 5.

    Sampson A. Τhe mesolithic hunter-gatherers in the Southeastern Mediterranean and their contribution in the Neolithisation of the Aegean. Archaeol Cult. 2018;1:11–36.

    Google Scholar 

  6. 6.

    Ripoll MP. The Knossos fauna and the beginning of the neolithic in the Mediterranean islands. In: Efstratiou N, Karetsou A, Ntinou M, editors. The neolithic settlement of Knossos in Crete: New evidence for the early occupation of Crete and the Aegean islands. Philadelphia: INSTAP Academic Press; 2013. p. 133–69.

    Google Scholar 

  7. 7.

    Peters J, Arbuckle BS, Pöllath N. Subsistence and beyond: Animals in neolithic Anatolia. In: Özdoğan M, Başgelen N, Kuniholm P, editors. The Neolithic in Turkey. 10500-5200 BC: environment, settlement, flora, fauna, dating, symbols of belief, with views from North, South, East, and West. Istanbul: Archaeology and Art Publications; 2014. p. 135–203.

  8. 8.

    Arbuckle BS, Kansa SW, Kansa E, Orton D, Cakirlar C, Gourichon L, et al. Data sharing reveals complexity in the westward spread of domestic animals across neolithic Turkey. PLoS ONE. 2014;9:e99845.

    PubMed  PubMed Central  Google Scholar 

  9. 9.

    Lenstra JA, Felius M. Genetic aspects of domestication. In: Garrick DJ, Ruvinsky A, editors. In The genetics of cattle. 2nd ed. Wallingford: CABI; 2014. p. 19–32.

    Google Scholar 

  10. 10.

    Chang C. Pastoral transhumance in the Southern Balkans as a social ideology: Ethnoarcheological research in Northern Greece. Am Anthropol. 1993;95:687–703.

    Google Scholar 

  11. 11.

    Hadjigeorgiou I. Past, present and future of pastoralism in Greece. Pastor Res Pol Pract. 2011;1:24.

    Google Scholar 

  12. 12.

    Hedrick PW. Genetics of populations. 3rd ed. Sudbury: Jones and Bartlett Publishers; 2005.

    Google Scholar 

  13. 13.

    Simčič M, Smetko A, Sölkner J, Seichter D, Gorjanc G, Kompan D, et al. Recovery of native genetic background in admixed populations using haplotypes, phenotypes, and pedigree information—using Cika cattle as a case breed. PLoS ONE. 2015;10:e0123253.

    PubMed  PubMed Central  Google Scholar 

  14. 14.

    Bett RC, Okeyo MA, Malmfors B, Johansson K, Agaba M, Kugonza DR, et al. Cattle breeds: extinction or quasi-extant? Resources. 2013;2:335–57.

    Google Scholar 

  15. 15.

    Taberlet P, Valentini A, Rezaei HR, Naderi S, Pompanon F, Negrini R, Ajmone-Marsan P. Are cattle, sheep, and goats endangered species? Mol Ecol. 2008;17:275–84.

    CAS  PubMed  Google Scholar 

  16. 16.

    Medugorac I, Veit-Kensch CE, Ramljak J, Brka M, Marković B, Stojanović S, et al. Conservation priorities of genetic diversity in domesticated metapopulations: a study in taurine cattle breeds. Ecol Evol. 2011;1:408–20.

    PubMed  PubMed Central  Google Scholar 

  17. 17.

    Ramljak J, Bunevski G, Bytyqi H, Marković B, Brka M, Ivanković A, et al. Conservation of a domestic metapopulation structured into related and partly admixed strains. Mol Ecol. 2018;27:1633–50.

    PubMed  Google Scholar 

  18. 18.

    Zervas N, Boyazoglou JG. La préservation et mise en valeur des populations locales en Méditerranée: état actuel du cheptel et des méthodes d’élevage en Grèce. Paris: Société d’EthnoZoologie; 1977.

    Google Scholar 

  19. 19.

    Constantinou A. Ruminant livestock genetic resources in Cyprus. Anim Genet Resour Inf. 1985;4:1–8.

    Google Scholar 

  20. 20.

    Mason IL. A world dictionary of livestock breeds, types and varieties. Wallingford: CABI; 1988.

    Google Scholar 

  21. 21.

    Bizelis I. Current situation of autochthonous breeds in Greece. 2015. Accessed 15 Mar 2019.

  22. 22.

    Domestic Animal Diversity Information System (DAD-IS). FAO. 2019. Accessed 13 Mar 2019.

  23. 23.

    FAO. The second report on the state of the world’s animal genetic resources for food and agriculture. 2015. Accessed 10 Feb 2018.

  24. 24.

    Rothammer S, Seichter D, Förster M, Medugorac I. A genome-wide scan for signatures of differential artificial selection in ten cattle breeds. BMC Genomics. 2013;14:908.

    PubMed  PubMed Central  Google Scholar 

  25. 25.

    Decker JE, McKay SD, Rolf MM, Kim J, Molina Alcalá A, Sonstegard TS, et al. Worldwide patterns of ancestry, divergence and admixture in domesticated cattle. PLoS Genet. 2014;10:e1004254.

    PubMed  PubMed Central  Google Scholar 

  26. 26.

    Ramljak J, Ivankovic A, Veit-Kensch CE, Förster M, Medugorac I. Analysis of genetic and cultural conservation value of three indigenous Croatian cattle breeds in a local and global context. J Anim Breed Genet. 2011;128:73–84.

    CAS  PubMed  Google Scholar 

  27. 27.

    Simčič M, Lenstra JA, Baumung R, Dovč P, Cepon M, Kompan D. On the origin of the Slovenian Cika cattle. J Anim Breed Genet. 2013;130:487–95.

    PubMed  Google Scholar 

  28. 28.

    Gautier M, Laloë D, Moazami-Goudarzi K. Insights into the genetic history of French cattle from dense SNP data on 47 worldwide breeds. PLoS ONE. 2010;6:e13038.

    Google Scholar 

  29. 29.

    Kukučková V, Moravčíková N, Ferenčaković M, Simčič M, Mészáros G, Sölkner J, et al. Genomic characterization of Pinzgau cattle: genetic conservation and breeding perspectives. Conserv Genet. 2017;18:839–910.

    Google Scholar 

  30. 30.

    Mastrangelo S, Ciani E, Ajmone Marsan P, Bagnato A, Battaglini L, Bozzi R, et al. Conservation status and historical relatedness of Italian cattle breeds. Genet Sel Evol. 2018;50:35.

    PubMed  PubMed Central  Google Scholar 

  31. 31.

    Groeneveld LF, Lenstra JA, Eding H, Toro MA, Scherf B, Pilling D, et al. Genetic diversity in farm animals—a review. Anim Genet. 2010;41:6–31.

    PubMed  Google Scholar 

  32. 32.

    Medugorac I, Medugorac A, Russ I, Veit-Kensch CE, Taberlet P, Luntz B, et al. Genetic diversity of European cattle breeds highlights the conservation value of traditional unselected breeds with high effective population size. Mol Ecol. 2009;18:3394–410.

    PubMed  Google Scholar 

  33. 33.

    Cymbron T, Freeman AR, Malheiro MI, Vigne JD, Bradley DG. Microsatellite diversity suggests different histories for Mediterranean and Northern European cattle populations. Proc Biol Sci. 2005;272:1837–43.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. 34.

    Flori L, Moazami-Goudarzi K, Alary V, Araba A, Boujenane I, Boushaba N, et al. A genomic map of climate adaptation in Mediterranean cattle breeds. Mol Ecol. 2019;28:1009–29.

    CAS  PubMed  Google Scholar 

  35. 35.

    Verdugo MP, Mullin VE, Scheu A, Mattiangeli V, Daly KG, Maisano Delser P, et al. Ancient cattle genomics, origins, and rapid turnover in the Fertile Crescent. Science. 2019;365:173–6.

    CAS  PubMed  Google Scholar 

  36. 36.

    Papadopoulos D. The cattle of the Kea Island. Thessaloniki: Editions of the Agricultural Bank of Greece (in Greek); 1946.

    Google Scholar 

  37. 37.

    Appuhn Κ. Ecologies of beef: eighteenth-century epizootics and the environmental history of early modern Europe. Environ Hist. 2010;15:268–87.

    Google Scholar 

  38. 38.

    Browning B, Browning S. Genotype imputation 631 with millions of reference samples. Am J Hum Genet. 2016;98:116–26.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Browning S, Browning B. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. Am J Hum Genet. 2007;81:1084–97.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Yang J, Benyamin B, McEvoy BP, Gordon S, Henders AK, Nyholt DR, et al. Common SNPs explain a large proportion of the heritability for human height. Nat Genet. 2010;42:565–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Granato ISC, Galli G, de Oliveira Couto EG, Bandeira e Souza M, Freitas Mendonça L, Fritsche-Neto R. SnpReady: a tool to assist breeders in genomic analysis. Mol Breed. 2018;38:102.

    Google Scholar 

  42. 42.

    Filzmoser P, Garrett RG, Reimann C. Multivariate outlier detection in exploration geochemistry. Comput Geosci. 2005;31:579–87.

    CAS  Google Scholar 

  43. 43.

    Matukumalli LK, Lawley CT, Schnabel RD, Taylor JF, Allan MF, Heaton MP, et al. Development and characterization of a high density SNP genotyping assay for cattle. PLoS ONE. 2009;4:e5350.

    PubMed  PubMed Central  Google Scholar 

  44. 44.

    El Mousadik A, Petit RJ. High level of genetic differentiation for allelic richness among populations of the argan tree [Argania spinosa (L) Skeels] endemic to Morocco. Theor Appl Genet. 1996;92:832–9.

    PubMed  Google Scholar 

  45. 45.

    Nei M. Molecular evolutionary genetics. New York: Columbia University Press; 1987.

    Google Scholar 

  46. 46.

    Jost L. GST and its relatives do not measure differentiation. Mol Ecol. 2008;17:4015–26.

    Google Scholar 

  47. 47.

    Barbato M, Orozco-terWengel P, Tapio M, Bruford MW. SNeP: a tool to estimate trends in recent effective population size trajectories using genome-wide SNP data. Front Genet. 2015;6:109.

    PubMed  PubMed Central  Google Scholar 

  48. 48.

    Sved JA, Feldman MW. Correlation and probability methods for one and two loci. Theor Popul Biol. 1973;4:129–32.

    CAS  PubMed  Google Scholar 

  49. 49.

    Pickrell JK, Pritchard JK. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 2012;8:e1002967.

    CAS  PubMed  PubMed Central  Google Scholar 

  50. 50.

    Nei M, Tajima F, Tateno Y. Accuracy of estimated phylogenetic trees from molecular data. J Mol Evol. 1983;19:153–70.

    CAS  PubMed  Google Scholar 

  51. 51.

    Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23:254–67.

    CAS  PubMed  PubMed Central  Google Scholar 

  52. 52.

    Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated animals. Genet Res. 2009;19:1655–64.

    CAS  Google Scholar 

  53. 53.

    Alexander DH, Lange K. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinf. 2011;12:246.

    Google Scholar 

  54. 54.

    Behr AA, Liu KZ, Liu-Fang G, Nakka P, Ramachandran S. Pong: fast analysis and visualization of latent clusters in population genetic data. Bioinformatics. 2016;32:2817–23.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Bowcock AM, Ruiz-Linares A, Tomfohrde J, Minch E, Kidd JR, Cavalli-Sforza LL. High resolution of human evolutionary trees with polymorphic microsatellites. Nature. 1994;368:455–7.

    CAS  PubMed  Google Scholar 

  56. 56.

    Kruskal JB, Wish M. Quantitative applications in the social sciences: multidimensional scaling. Thousand Oaks: SAGE Publications, Inc.; 1978.

    Google Scholar 

  57. 57.

    Green RE, Krause J, Briggs AW, Maricic T, Stenzel U, Kircher M, et al. A draft Ssequence of the Neandertal genome. Science. 2010;328:710–22.

    CAS  PubMed  PubMed Central  Google Scholar 

  58. 58.

    Patterson N, Moorjani P, Luo Y, Mallick M, Rohland N, Zhan Y, et al. Ancient admixture in human history. Genetics. 2012;192:1065–93.

    PubMed  PubMed Central  Google Scholar 

  59. 59.

    Ørsted M, Hoffmann AA, Sverrisdóttir E, Nielsen KL, Kristensen TN. Genomic variation predicts adaptive evolutionary responses better than population bottleneck history. PLoS Genet. 2019;15:e1008205.

    PubMed  PubMed Central  Google Scholar 

  60. 60.

    Vilas A, Pérez-Figueroa A, Queasada H, Caballero A. Allelic diversity for neutral markers retains a higher adaptive potential for quantitative traits than expected heterozygosity. Mol Ecol. 2015;24:4419–32.

    PubMed  Google Scholar 

  61. 61.

    Spilanis G, Kizos A. The Atlas of the islands (in Greek). Mytilene: Editions of the University of the Aegean; 2015.

    Google Scholar 

  62. 62.

    Statistique nationale de La Grèce. Résultats du recensement de l’Agriculture de la Grèce: Année 1951, 1961, 1991, 2001. Athènes: Imprimerie Nationale.

  63. 63.

    Upadhyay M, Bortoluzzi C, Barbato M, Ajmone-Marsan P, Colli L, Ginja C, et al. Deciphering the patterns of genetic admixture and diversity in southern European cattle using genome-wide SNPs. Evol Appl. 2019;12:951–63.

    CAS  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Richard E, Lenski RE. An experimental test of evolutionary trade-offs during temperature adaptation. Proc Natl Acad Sci USA. 2007;104:8649–54.

    Google Scholar 

  65. 65.

    Cebeci AF, Ince H, Mercan MA. The intrinsic fallacy of market mechanism and private property rights in alleviating the tragedy of the commons. Appl Econ. 2020.

    Article  Google Scholar 

  66. 66.

    Hardin G. The tragedy of the commons. Science. 1968;162:1243–8.

    CAS  PubMed  Google Scholar 

  67. 67.

    Medugorac I, Graf A, Grohs C, Rothammer S, Zagdsuren Y, Gladyr E, et al. Whole-genome analysis of introgressive hybridization and characterization of the bovine legacy of Mongolian yaks. Nat Genet. 2017;49:470–5.

    CAS  PubMed  Google Scholar 

  68. 68.

    Decker JE, Pires JC, Conant GC, et al. Resolving the evolution of extant and extinct ruminants with high-throughput phylogenomics. Proc Natl Acad Sci USA. 2009;106:18644–9.

    CAS  PubMed  Google Scholar 

  69. 69.

    Iso-Touru T, Sahana G, Guldbrandtsen B, Lund MS, Vilkki J. Genome-wide association analysis of milk yield traits in Nordic Red Cattle using imputed whole genome sequence variants. BMC Genet. 2016;17:55.

    CAS  PubMed  PubMed Central  Google Scholar 

  70. 70.

    Upadhyay MR, Chen W, Lenstra JA, Goderie CRJ, MacHugh DE, Park SDE, et al. Genetic origin, admixture and population history of aurochs (Bos primigenius) and primitive European cattle. Heredity. 2017;118:169–76.

    CAS  PubMed  Google Scholar 

  71. 71.

    Pitcairn A. The agricultural resources of Cyprus: the effect of natural and other factors on development. Cyprus Agric J. 1935;30:6–18.

    Google Scholar 

  72. 72.

    Bevan W. Notes on agriculture in Cyprus and its products. 1919. Accessed 19 Jul 2019.

  73. 73.

    Cypriot ministry of agriculture, rural development and environment. Census of the local cattle breed (in Greek); 2016.$file/%CE%95%CE%A0%CE%99%CE%A3%CE%9A%CE%9F%CE%A0%CE%97%CE%A3%CE%97%20%CE%9D%CE%A4%CE%9F%CE%A0%CE%99%CE%91%CE%A3%20%CE%A6%CE%A5%CE%9B%CE%97%CE%A3%20%CE%92%CE%9F%CE%9F%CE%95%CE%99%CE%94%CE%A9%CE%9D%202016.pdf?OpenElement. Accessed 19 Jul 2019.

  74. 74.

    Hatziolos B. The issue of animal production in Greece (in Greek). Athens: Greek Publishing Company; 1941.

    Google Scholar 

  75. 75.

    Keller C. Studien über die Haustiere der Mittelmeer-Inseln. Zürich: Kommissions-Verlag Von Georg & Cie; 1911.

    Google Scholar 

  76. 76.

    Papadopoulos D. The cattle of the Kea Island (in Greek). Thessaloniki: Editions of the Agricultural Bank of Greece; 1946.

    Google Scholar 

  77. 77.

    Manetti C. L’Anatolia. Firenze: R Bemporad & Figlio; 1922.

    Google Scholar 

  78. 78.

    Dimitriadis IN. The cattle of West Macedonia (in Greek). Thessaloniki: Editor Triantafyllou M; 1933.

  79. 79.

    Karantounias A. Specific zootechnia, cattle and buffalo breeding (in Greek). Athens: Editions of the Highest Agricultural School of Athens; 1967.

    Google Scholar 

  80. 80.

    Keller O. Die antike Tierwelt. Leipzig: Wilhem Engelmann; 1909.

    Google Scholar 

  81. 81.

    Kazoglou Y. Development of pastures’ management and livestock supporting tools for the ecosystems of the Prespa Municipality (Prespa National Park, Ladopotamos Valley, Krystallopigi Municipality). Deliverables of the second reporting period January-March 2015 (in Greek), under the frame of “Development of research and technological development innovation projects (AgroETAK)”. KYPE 3324/102 (MIS: 453350). Hellenic Agricultural Organization “DEMETER”; 2015.

  82. 82.

    Tsaprailis A, Kazoglou Y. Extensive breeding of the Greek shorthorn cattle breed: high quality products from the natural pastures to the consumer’s plate (in Greek with summary in English). In Proceedings of the Food Expo 2017 Meat and Products: 18–20 March 2017. Athens; 2017.

  83. 83.

    Vezzani V. Rodi e il suo Problema Zootecnico. Milano: Natura; 1929.

    Google Scholar 

  84. 84.

    Kazoglou IE, Xega N, Logotheti Α, Doleson F. Rare bovine breeds in the transboundary Prespa Park (in Greek with summary in English). In: Sidiropoulou A, Mantzanas K, Ispikoudis I, editors. In: Proceedings of the 7th Panhellenic Rangeland Congress Conference Grassland Science: 14–16 October 2010; Xanthi; 2010.

  85. 85.

    Grünenfelder HP, Trivizaki S. Description of Prespa cattle (Mistrece). St. Gallen: SAVE Foundation, Monitoring Institute for Rare Breeds and Seeds in Europe Report; 2014.

  86. 86.

    Psaltis G. Our cattle strains (in Greek). Athens: Editor Apatsidis N; 1931.

  87. 87.

    La Hatziolos B. race bovine de Sykia (Chalcidique) (in Greek and French). Athens: Ministère de l’Agriculture Report; 1935.

    Google Scholar 

Download references


Part of the samples was collected and financed by the BUSHALIVE project (FAO, PO 299432). We acknowledge all breeders and breeding organizations that provided samples. Exchange of researchers from Greece and Germany for the current study was also funded by the German Academic Exchange Service (DAAD) under the running project entitled “Different forces behind the development of genomic diversity in indigenous cattle strains and breeds in Greece and Germany”. Hellenic Foundation for Research and Innovation (HFRI) supported this work with a fellowship to DP under the HFRI Ph.D. Fellowship frame (Fellowship Number: 1301). The authors would also like to acknowledge members of AMALTHIA (Stamatina Trivizaki, Kostas Papaioannou, Yiannis Kazoglou, Vasilios Lekkas, Achilleas Tsaprailis, Stefano Dellepiane) for their useful information and their support during the implementation of the project.


Part of the samples was collected and financed by the BUSHALIVE project (FAO, PO 299432). Exchange of researchers from Greece and Germany for the current study was funded by the German Academic Exchange Service (DAAD) under the running project entitled “Different forces behind the development of genomic diversity in indigenous cattle strains and breeds in Greece and Germany” (project code 57418683, IKYDA program between the DAAD and the Greek State Scholarship Foundation (I.K.Y.)).

Author information




DP carried out data analysis, interpreted data and drafted the manuscript. PK and GPL assisted in data interpretation and drafting the manuscript. EK assisted with field work, data analysis and critically revised the manuscript. MU carried out the D-statistics analyses, interpreted data and revised the manuscript. DS and IR coordinated the genotyping of samples and provided data. GB contributed samples and genotypes. NK contributed samples and information on characteristics of the studied Greek breeds. IB designed the study, provided samples and critically revised the manuscript. IM conceived and guided this study, assisted with field work, contributed analysis tools and genotype data, analyzed data and critically revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Ivica Medugorac.

Ethics declarations

Ethics approval and consent to participate

Ethical permission for hair root and blood sampling was given by the ethics committee of the Agricultural University of Athens. All samples were collected by trained personnel according to best veterinary practice.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Sample description, group allocation, RGB (color) code assigned to the pre-defined groups within this paper, breed names, breed code, number of sampled and genotyped individuals (N), number of genotyped and unrelated individuals used to estimate the diversity parameters (Nd), current breeding purposes as well as sporadic or recent past breeding purposes in parenthesis, breed origin, and source of the samples or genotypes used in this study or from previous studies [17, 25, 28, 30, 34, 35, 67,68,69,70]. Table S2. Phenotypic, productive and reproductive traits. Description of phenotype, productive and reproductive traits of the 11 Greek and Cypriot analyzed breeds. Table S3. Evolution of bovine population sizes on the islands during the last 60 years. Statistical information regarding the evolution of cattle population on the Greek island based on [62].

Additional file 2.

Detailed description of the indigenous Greek and Cypriot cattle. Compilation of information on Greek indigenous cattle breeds from diverse sources [16, 33, 71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87]. This includes also information from records in Greek, which were formerly not accessible to the international scientific community due to language barriers.

Additional file 3: Figure S1.

Tessellated projection. Spatial geographic presentation of the herein estimated diversity parameters (mA, AR, HE, HE(SNP), npA, nspA).

Additional file 4: Figure S2.

Tessellated projection. Spatial geographic presentation of the estimated effective population number (Ne5, Ne50, Ne2000).

Additional file 5: Table S4.

Genetic differentiation index. Pairwise values for DEST, per analyzed breed. Maximum and minimum values within each of the predefined geographical groups are depicted in bold letters in gray rectangle. Maximum and minimum values considering all 115 breeds are in bold letters in green and yellow, respectively.

Additional file 6: Figure S3.

Tessellated projection. Spatial geographic presentation of the estimated allele sharing distance matrix (DPS) among breeds using multi-allelic SNP-blocks.

Additional file 7: Figure S4.

Phylogenetic tree. Neighbor-joining tree based on Nei’s genetic distance DA using multi-allelic SNP blocks. Mongolian yak (YAK) was used as a root. Dotted lines indicate the reduced length of YAK and GIR to improve visibility. Special square marks represent the influence of East-Podolian (grey), Alpine (green) and North-West (olive green) group.

Additional file 8: Figure S5.

Admixture analysis. Population structure of the 114 studied breeds/populations presented for each of K inferred ancestral population value, using a genome‐wide set of bi-allelic SNPs. All the analyzed individuals are presented as colored vertical bars, divided at most K segments and with proportional heights, according to their genotype membership. In the depicted plots (K = 2 to K = 26), a different color is presented for each cluster. From left to right within each group the presented breeds are as follows: Minor Asia: ATER, ATBC, ATSR, ATSY, TRG; Greece and Cyprus: CYP, AGT, CRT, NSY, GRB, KAS, KEA, PRG, ROG, SYK, KTR; South East Europe: RHS, MKB, SRB, PRE, RMB, SHB, DGB, DBB, MAB, LKB, SKB, MNB, BHB, HRB; East Podolian: HRI, HRP, UKP; Tyrrhenian: PODO, CINI, MOSI, RSIC, MOSA, SARD, SBRU, CORS, AGER, MARE, CHI, MPIS, CALV, MCH, RMG, GARF, PONT, MODE, CABA, REGG, PMT, BURL; Alpine: PRDO, OVAR, REND, BPUS, PUST, SIC, PIN, TGV, MWF, OBV, BBV, DFV, FGV, VOG, ABO, MON, TAR; France: RDBI, SAL, AUB, LIM, CHR, PAR, BAQ, GAS; Iberian: MNRQ, MALL, NGAN, CANA, MARI, ALEN, BAR, MARO, SYG; North West Europe: BPN, NOR, MAN, BBB, LKF, HF, GNS, JSY, HER, SHR, KRY, DXT, GLW, AAN, HGL, NRC, SERC, FJL, FIAY, FINE, FINW, FINN, YARO.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Papachristou, D., Koutsouli, P., Laliotis, G.P. et al. Genomic diversity and population structure of the indigenous Greek and Cypriot cattle populations. Genet Sel Evol 52, 43 (2020).

Download citation