ORIGINAL RESEARCH published: 22 January 2021 doi: 10.3389/fgene.2020.610353 Genetic Diversity of Bubalus bubalis in Germany and Global Relations of Its Genetic Background Antonia Noce 1 , Saber Qanbari 1 , Rayner González-Prendes 2 , Julia Brenmoehl 1 , María Gracia Luigi-Sierra 3 , Michael Theerkorn 4 , Marc-André Fiege 5 , Heike Pilz 6 , Adrian Bota 7 , Livia Vidu 8 , Csaba Horwath 9 , László Haraszthy 10 , Pencho Penchev 11 , Yordanka Ilieva 11 , Tzonka Peeva 11 , Wolfgang Lüpcke 12 , René Krawczynski 13 , Klaus Wimmers 1 , Manfred Thiele 14 and Andreas Hoeflich 1* 1 Leibniz-Institute for Farm Animal Biology, Dummerstorf, Germany, 2 Animal Breeding and Genomics Group, Wageningen University & Research, Wageningen, Netherland, 3 Department of Animal Genetics, Center for Research in Agricultural Genomics, Campus Universitat Autonoma de Barcelona, Bellaterra, Spain, 4 Hof am Meer GmBH, Stadland, Germany, 5 Gut Darss, GmbH & Co. KG, Born, Germany, 6 Wiesenburger Land eG, Wiesenburg, Germany, 7 Research and Development Station for Buffalos Şercaia, Şercaia, Romania, 8 University of Agronomic Sciences and Veterinary Medicine of Bucharest, Bucharest, Romania, 9 Comuna Baciu sat Mera, Cluj, Romania, 10 Pro Vértes, Csákvár, Hungary, 11 Bulgarian National Association for Development of Buffalo Breeding, Shumen, Bulgaria, 12 Higher School in Agribusiness and Development of Regions, Agricultural University Plovdiv, Tarnovo, Bulgaria, 13 Energiequelle GmbH, Zossen, Germany, 14 German Buffalo Edited by: Association, Penig, Germany Muniyandi Nagarajan, Central University of Kerala, India Reviewed by: This is the first study to explore the genetic diversity and population structure of domestic Tingxian Deng, water buffalo (Bubalus bubalis) in Germany and their potential relations to herds in other Institute of Buffalo (CAAS), China Liguo Yang, parts of Europe or worldwide. To this end, animals from different herds in Germany, Huazhong Agricultural University, Bulgaria, Romania, and Hungary were genotyped and compared to genotypes from China other populations with worldwide distribution and open to the public. The pilot study *Correspondence: analyzed population structure, phylogenetic tree, and inbreeding events in our samples. Andreas Hoeflich hoeflich@fbn-dummerstorf.de In buffalos from Germany, a mixed genetic make-up with contributions from Bulgaria (Murrah breed), Romania, and Italy was found. All in all, a high degree of genetic diversity Specialty section: was identified in European buffalos, and a novel genotype was described in Hungarian This article was submitted to Evolutionary and Population Genetics, buffalos by this study. We demonstrate that European buffalos stand out from other a section of the journal buffalo populations worldwide, supporting the idea that buffalos have not completely Frontiers in Genetics disappeared from the European continent during the late Pleistocene. The high genetic Received: 25 September 2020 Accepted: 30 November 2020 diversity in European buffalos seems to be an excellent prerequisite for the establishment Published: 22 January 2021 of local breeds characterized by unique traits and features. This study may also be Citation: considered as an initial step on the way to genome characterization for the sustainable Noce A, Qanbari S, development of the buffalo economy in Germany and other parts of Europe in the future. González-Prendes R, Brenmoehl J, Luigi-Sierra MG, Theerkorn M, Keywords: Bubalus bubalis, European populations, genotyping, genetic diversity, population structure, run of Fiege M-A, Pilz H, Bota A, Vidu L, homozygosity Horwath C, Haraszthy L, Penchev P, Ilieva Y, Peeva T, Lüpcke W, Krawczynski R, Wimmers K, Thiele M and Hoeflich A (2021) Genetic INTRODUCTION Diversity of Bubalus bubalis in Germany and Global Relations Water buffalo (Bubalus bubalis) is a multipurpose animal, producing meat, milk, leather products, of Its Genetic Background. and dung. In developing countries, buffalos are used as draught animals, providing more than 40% Front. Genet. 11:610353. of farm labor (Borghese, 2005). Due to their wallowing behavior, buffalos are less susceptible to doi: 10.3389/fgene.2020.610353 ectoparasites and related diseases and suitable for grazing in the swamps (Michelizzi et al., 2010). Frontiers in Genetics | www.frontiersin.org 1 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany Buffalos are also more efficient in the digestion process due to in refugial areas. Because rock reliefs, dated back to ≈16,000 their longer rumen and produce a lower average of methane BC in France, clearly display buffalos (Krawczynski, 2013). This, emission (Moss et al., 2000; Kawashima et al., 2006; Calabró in fact, would be an indication that water buffalos in Europe et al., 2013). Compared to cattle, the buffalo is less selective for survived the ice age. feed quality (Michelizzi et al., 2010) and longer living, granting Domesticated water buffalos in Europe are referenced during husbandry at lower costs and with less frequent animal turnover the centuries VI–XII C.E. (Maymone, 1942; Alexiev, 1998) and (Borghese, 2005; Hegde, 2019; Hoffman and Valencak, 2020). had been settled in Italy, Bulgaria, Romania, and all the other Buffalo meat is a rich source of proteins, fatty acids as omega 3- Balkan countries, where they have been preserved and bred 6 fatty acids, iron, and characterized by lower concentrations of for centuries until today (Ferrara, 1964). European buffalos are cholesterol if compared to cattle (Roth and Myers, 2004). Milk river type (Borghese, 2011). The Balkans (East Europe) are an from buffalo contains biliverdin, bioactive pentasaccharides, and essential point of the river buffalo’s historical migration route. gangliosides, which are not present in bovine milk (Abd El-Salam Riverine buffalo moved to Southwestern Asia from the Indian and El-Shibiny, 2011) and more fat in total but lower cholesterol domestication area, reached Egypt and Turkey, and arrived compared to the cattle. Buffalo milk is further characterized by in East Europe and Italy during the seventh century (Zeuner, higher protein content and big casein micelles; both features 1963; Clutton-Brock, 1999). Some animals likely returned then are related to the higher yields of cheese produced from the to Egypt, Turkey, and Bulgaria with Crusaders returning, same amounts of buffalo versus dairy cow milk (Zicarelli, 2004; and spread into the other Balkan regions during the 12th Martini et al., 2018). century (Ferrara, 1964; Borghese, 2005). Since 1980 buffalos The domestic water buffalo is classified into two major have also been introduced to Germany, France, Spain, Portugal, categories: swamp buffalo (Bubalus bubalis carabanensis, Luxembourg, The Netherlands, Switzerland, and other parts 2n = 48) and river buffalo (Bubalus bubalis bubali, 2n = 50). of Europe. Nowadays, the only officially recognized breeds in While their taxonomical status is still being debated, Europe are the Mediterranean Italian buffalo and the Bulgarian B. bubalis is supposed to be descended from Indian Murrah. Selecting animals for milk traits was initiated in Italy. wild buffalo (Bubalus arnee) domesticated ≈5,000 years Crossing with other populations was avoided; a herd book with ago (Cockrill, 1984). Nevertheless, the origin of the two animals, performance, and morphological traits was established subspecies into which water buffalo are divided is still the to maintain its unique genetic identity (Iannuzzi and Di Meo, object of study (Colli et al., 2018b), especially to solve the 2009). The primary purpose was the production and marketing debate about the occurrence of a single versus two separate of milk and milk derivatives such as mozzarella cheese, one of the domestication events (Lau et al., 1998; Kierstein et al., 2004; “pasta filata” cheeses known worldwide (Zicarelli, 2001). Kumar et al., 2007). While local Bulgarian buffalos were a significant European Little is known about the history, genetic diversity, and Mediterranean population raised for draft power, meat, and milk performance of European buffalo populations. The presence of purposes, they started to be crossed with Murrah breed since water buffalo in Europe is dated from the Pleistocene until the 1972, mainly to improve their performance for milk production. warm period before the last Eem-Interglacial (≈125,000 years A selection program was initiated to develop a typical Bulgarian ago), when buffalos were also present in Central Europe Murrah breed for milk with high butterfat content (Borghese, (Zahariev et al., 1986; Alexiev, 1998). Although hunting 2005; Borghese and Moioli, 2016). (Martin, 1984) and climatic changes, particularly during the late In Germany, 7,614 buffalos, spread over 16 different Pleistocene (Lima-Ribeiro and Diniz-Filho, 2017), temporarily or regions, are recorded for the year 2020 by an internal locally displaced and shrank the European buffalo populations, communication from the German Buffalo Breeder’s Association. the presence of refugial regions in Southern and Eastern The development of the single nucleotide polymorphism (SNP) Europe may have prevented the extinction of the Bubalus genotyping assay specific for river buffalo (Iamartino et al., 2017), species from the European continent (Krawczynski, 2010). The Axiom Buffalo Genotyping Array 90K from Affymetrix, now R presence of buffalos in South-Eastern Europe during the early enables genome analysis of this comparably novel farm animal Neolithikum (9,000–7,000 BC) was suggested (Bökönyi, 1957) species and could be useful for the initiation of genetic selection but generally dismissed (Bartosiewicz, 1999). Based on bone in German buffalo. A subsequent release of the first assembly finds in Austria dated to the Atlantikum (7,000–4,000 BC), of the water buffalo genome (Low et al., 2019) further enabled other authors have also discussed the presence of buffalos researchers to expand their knowledge about this on a molecular during the Holocene (beginning 12,000 BC) (Pucher, 1988). scale. In this study, we sought to explore the genetic diversity However, defined criteria for the unambiguous distinction of and structure of domestic water buffalo populations in Germany members from the Bovini tribus (in particular, Bos versus and Eastern Europe. Bubalus) have been lacking until now (Pucher, 1988). According to Haarmann (2011), the ancient Greeks have used different terms for aurochs (Bos primigenius), wisent (Bos bonasus), MATERIALS AND METHODS and buffalo (Bubalus sp.). This linguistic approach supports the hitherto controversial bone finds from the same area at Ethics Statement about the same time. It is unclear whether water buffalos All sampling occurred as part of commercial buffalo breeding repopulated Europe after the ice age or whether they survived programs and strictly adhered to national and international laws. Frontiers in Genetics | www.frontiersin.org 2 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany FIGURE 1 | Geographic map showing the locations with specific coordinates the buffalo samples were collected from, genotyped in the current study. The number of samples for each location are reported in Table 1. Samples Collection, DNA Purification, when merging genotypes, we retained only those SNPs that were and Genotyping common in the three data sets. We further removed SNPs which mapped on sexual and mitochondrial bases, as well as those with Ear tag test samples from 285 female and male buffalos were unknown chromosome coordinates to perform an autosomal collected randomly from ten farms during the years 2018 analysis. The remaining panel of 36,759 SNPs was then updated and 2019. Importantly, all breeders were asked to provide to the last version of the water buffalo genome (UOA_WB_1) samples from unrelated animals. The farms are distributed in through the commands –update-map, –update-chr, –update- Central and Eastern Europe (Figure 1), as listed in Table 1. alleles of PLINK V1.7 software (Purcell et al., 2007). The three DNeasy Blood & Tissue Kit (Qiagen, Germany) was used to datasets were then merged using the merge module in PLINK. extract genomic DNA according to the instructions of the Individually, the quality control filters were applied to remove manufacturer. Purification and concentration of DNA were variants with a minor allele frequency lower than 0.05 or a rate of measured with a NanoDrop2000 before the normalization to missing genotyping >10%. The final dataset consisted of 36,014 the required concentration (30 ng/µL). The quality control and genotyping using the Axiom Buffalo Genotyping Array 90K R from Affymetrix1 were performed by ATLAS Biolabs GmbH (Berlin, Germany). Allele calling was carried out using Axiom TABLE 1 | River buffalo populations analyzed in the present study (n: samples Analysis Suite software V4.0.1 (Applied Biosystems by Thermo per population). Fisher Scientific) following the pipeline for Affymetrix Axiom Populations label City Country n genotyping workflow (Nicolazzi et al., 2014) and using the last version of buffalo genome assembly (UOA_WB_1) as the 1 Bul_Noce Varna Bulgaria 58 reference (Low et al., 2019). 2 Ger_Born Born Germany 28 3 Ger_Jüt Jüterbog 27 Datasets Updating and Merging 4 Ger_Stad Stadland 26 5 Ger_Wies Wiesenburg 28 We retrieved previously reported datasets from Colli et al. 6 Hun_HT Földes Hungary 19 (2018a) and Deng et al. (2019a) available in the Dryad Repository 7 Hun_Csak Csákvar 17 (Table 2).2 To solve the incongruity between the marker positions 8 Hun_Tisz Tiszataj 19 1 http://www.affymetrix.com 9 Rom_Mera Mera Romania 16 2 https://datadryad.org/resource/doi:10.5061/dryad.h0cc7; https://doi.org/10.5061/ 10 Rom_Serc Sercaia 47 dryad.310pf05 Total number 285 Frontiers in Genetics | www.frontiersin.org 3 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany TABLE 2 | Populations of buffalo with publicly available genotype data included in 2009) with default parameters. This software modulates the the present genome analysis (n: samples per population). probability of the observed genotypes considering the ancestry Populations label Geographical origin n References proportion and population allele frequencies. The best number of clusters (K-value) is estimated by the model as that reporting the 1 Bul_Colli Bulgaria 11 Colli et al., 2018a lower cross-validation error. As an additional approach for the 2 Ita_Colli Italy 15 analysis of the population structure, a phylogenetic tree was built 3 Mozamb Mozambique 7 using the R package APE (Paradis and Schliep, 2019) based on 4 Rivbr_Ana Turkey 15 computed Wright’s Fixation Index (FST ) (Wright, 1965) matrix 5 Rivbr_Mur Brazil 15 obtained with Arlequin 3.5 (Excoffier and Lischer, 2010). 6 RivCo Colombia 12 7 RivEg Egypt 15 8 Rivir_Aza, Iran 27 Run of Homozygosity Rivir_Khu, Analysis of run of homozygosity (ROH) was exclusively Rivir_Maz conducted on Central and East Europe buffalo populations 9 RivPk_in_Mur India 5 genotyped in this study (Table 1) and based on a panel of 61,813 10 RivPk_Azk, Pakistan 28 autosomal SNPs available after the application of quality control RivPk_Kun, filters, described above, to eliminate sexual, mitochondrial, RivPk_Nil and unknown chromosome coordinate SNPs, as well as those 11 Rom_Colli Romania 9 with a minor allele frequency lower than 0.01 and rate of 12 Ita_Deng Italy 35 Deng et al., 2019a missing genotyping >10%. The detection of ROH was performed using PLINK V1.7 (Purcell et al., 2007), a sliding window of 1,000 kb was designed to detect regions matching the following SNPs and 477 individuals to perform the subsequent genetic parameters: minimum size of 50 SNPs (–homozyg-snp 50), population and structure analyses (Supplementary Table 1). minimum density of SNPs of 1 SNP every 100 kb (–homozyg- density 100). We allowed one heterozygous SNP (–homozyg- Genetic Relationship and Population window-het 1), a distance of homozygosity SNPs within the Structure window of 250 kb (–homozyg-gap 250), and two missing SNPs A multidimensional scaling (MDS) analysis was performed to (–homozyg-window-missing 2) per ROH. The identified ROH explore the genetic relationship between four different buffalo were cataloged in four categories according to its length: 1–5 Mb, populations. To this end, a symmetric matrix of the identity-by- 5–15 Mb, 15–30 Mb, and >30 Mb. Afterward, from ROH, we state (IBS) distances, for all pairs of individuals, was generated also calculated the inbreeding coefficient (FROH ) using R package in PLINK V1.7 software (Purcell et al., 2007) based on the “DetectRUN” (Biscarini et al., 2018) and applying the formula: proportion of alleles shared and visualized using ggplot R P package (Wickham, 2016). Observed (Ho) and Expected (He) LROH FROH = heterozygosity for ten buffalo populations genotyped in this study Lgenome (Bulgaria, Germany, Hungary, and Romania) were also estimated using PLINK V1.7 (Purcell et al., 2007). The maximum- 6LROH is the sum of the length of ROH for each individual, and likelihood approach was used to estimate the population Lgenome is the length of the genome analyzed, which was around structure through ADMIXTURE V1.3 Software (Alexander et al., 2.65 Gb in buffalo. FIGURE 2 | Methodical workflow chart summarizing the analysis pipeline followed in the current study. QCF = Quality control filters (MAF < 0.05, -geno 0.01, -mind 0.01). Frontiers in Genetics | www.frontiersin.org 4 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany FIGURE 3 | (A) Multi-Dimensional Scaling (MDS) plot based on genome-wide identity-by-state pairwise inferred with PLINK V1.7, including 36.014 SNPs from 477 individuals. This plot shows the genetic distance between all the 22 populations analyzed worldwide and published in part (Colli et al., 2018a; Deng et al., 2019a). The percentage of variances captured for each dimension is reported in brackets. The dashed circles mark the populations that cluster together. (B) Admixture analysis results at K = 4 and 6, including 22 buffalo populations worldwide. K values indicate the number of ancestries estimation (clusters), K = 6 was the best fitting solution for the lowest cross-validation error reported (CV = 0.61856). Populations are divided by a vertical black line, and it is partitioned into K colored segments that represent the population’s estimated membership fractions in K clusters. The 22 populations are localized in Germany (Ger_Wies, Wiesenburg, n = 28: Ger_Born, Born, n = 28; Ger_Stad, Stadland, n = 26; Ger_Jüt, Jüterbog, n = 27), Italy (Ita_Deng according to Deng et al., 2019a, n = 35; Ita_Colli according to Colli et al., 2018a, n = 15), Mozambique (Mozamb according to Colli et al., 2018a, n = 7), Bulgaria (Bul_Colli according to Colli et al., 2018a, n = 11), Bulgaria (Bul_Noce, n = 58), Hungary (Hun_HT, Földes, n = 19; Hun_Csak, Csákvar, n = 17; Hun_Tisz, Tiszataj, n = 19), Romania (Rom_Colli according to Colli et al., 2018a, n = 9; Rom_Mera, Mera, n = 16; Rom_Serc, Sercaia, n = 47), and also according to Colli et al., 2018a: Brazil (n = 15), Colombia (n = 12), Egypt (n = 15), Turkey (n = 15), Iran (n = 27), India (n = 5), and Pakistan (n = 28). Population labels and the number of individuals of both analyses are reported in Tables 1, 2. Frontiers in Genetics | www.frontiersin.org 5 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany FIGURE 4 | (A) Multidimensional scaling (MDS) plot based on European buffalo populations available in this work. The dashed circles mark the populations that cluster together. (B) Admixture analysis results at K = 4 and 6 exclusively in European buffalo populations. K = 6 was the best fitting solution for the lowest cross-validation error reported (CV = 0.59712). Abbreviations and sample numbers are specified in Figure 3. The frequency of occurrence of SNPs into an ROH was Graphical visualization was performed using R package qqman analyzed, identifying the genomic regions most commonly (Turner, 2014) by constructing Manhattan plots. For the associated with ROH. Regions were selected containing the top identification of genomic regions associated with ROH hotspots, 1% of most common ROH-associated SNPs (ROH hotspots or we used NCBI map viewer of the water buffalo UOA_WB_1. ROH islands) (Purfield et al., 2017; Mastrangelo et al., 2018). Functional annotation was performed by the use of DAVID Frontiers in Genetics | www.frontiersin.org 6 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany TABLE 3 | Observed (Ho) and expected (He) heterozygosity of buffalo populations Murrah buffalos as the only exception of a European population. de novo genotyped distributed in Bulgaria, Germany, Hungary, and Romania Four different groups could be distinguished in the rest of the (abbreviations are defined in Table 1). European population: the Hungarian cluster, which included Populations Ho He all three farms, the Romanian cluster, which showed one farm separated from the others (Rom_Serc), the Italian cluster, which Bul_Noce 0.41 0.40 included samples from Mozambique, and a German cluster from Ger_ Born 0.39 0.39 Ger _Jüt 0.39 0.35 Brandenburg (Ger_Jüt), which was separated from the other Ger _Stad 0.38 0.38 populations. The remaining animals from Germany were mixed Ger _Wies 0.40 0.38 and scattered between Romanian and Italian populations. The Hun_HT 0.37 0.34 admixture analysis carried out on genetic information from 22 Hun_Csak 0.35 0.35 buffalo populations available with global distribution (Figure 3B) Hun_Tisz 0.35 0.32 revealed a higher diversity among European populations than Rom_Mera 0.38 0.37 groups from India, Pakistan, the Middle East, and Latin America. Rom_Serc 0.38 0.36 K = 4 evidenced that the Bulgarian buffalos form the only European population that shared a high average of alleles with the non-European groups. The Romanian group showed mixed V6.8 (Huang et al., 2009) and PANTHER (Thomas et al., 2003) ancestry, except for Rom_Serc. Low admixture levels were software. The analysis pipeline applied in this study is provided found in the Hungarian group. A high admixture of buffalo was as a flow chart (Figure 2). found in all German farms, noticeably with a different structure in the buffalo from Brandenburg (Ger_Jüt). In Figure 4, we aimed to focus on the relative genetic distance and structure RESULTS exclusively among European populations (N = 14). The MDS plot (Figure 4A) showed that after excluding Mozambique buffalos, Analysis of Genetic Diversity and the Italian samples were presented as a more homogenous Population Structure and distinct cluster. The same appeared for the two Bulgarian The MDS plot (Figure 3A) of 22 buffalo populations (Tables 1, 2) populations. Furthermore, animals from two German farms with different geographic distributions revealed a clear separation (Lower Saxony – Ger_Stad, Mecklenburg-Western Pomerania – of European buffalo from other populations distributed Ger_Born) had relations to the Italian cluster, while animals worldwide. Indian-Pakistan and Latino-American buffalos from Saxony (Ger_Wies) appeared to be related more closely to were presented in a single cluster together with Bulgarian the Romanian buffalos. Buffalos from Brandenburg (Ger_Jüt) FIGURE 5 | Phylogenetic tree indicating the genetic distance of 22 buffalo populations, analyzed and published in part (Colli et al., 2018a; Deng et al., 2019a), with worldwide distributions based on the fixation index (FST) matrix obtained from 36,014 SNPs in 477 individuals using Arlequine software. Abbreviations and sample numbers are specified in Figure 3. Frontiers in Genetics | www.frontiersin.org 7 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany FIGURE 6 | The mean sum of runs of homozygosity (ROH) per population within each ROH length category in samples collected from buffalos in Bulgaria (Bul_Noce), Romania (Rom_Mera, Rom_Serc), Hungary (Hun_Csak, Hun_Tisz, Hun_HT), and Germany (Ger_Born, Ger_Jüt, Ger_Stad, Ger_Wies). Data include 61,813 SNPs genotyped in 285 samples (Mb: megabase); all other abbreviations and sample numbers are specified in Figure 3. could be positioned closer to the Bulgarian cluster. From K = 6 of worldwide. The phylogenetic tree based on the FST index matrix admixture analysis including exclusively European populations (Supplementary Table 2) was consistent with the Admixture (Figure 4B), three principal components were distinguished analysis (Figure 5). Two principal groups defined the difference in Central and Southeast Europe: the Mediterranean pattern, between the Italian Mediterranean and the Murrah breeds. including Italian and German populations, the Hungarian- Along with the Italian group, Mozambique and three German Romanian, and the Bulgarian ones. Buffalos from farms herds, including Mecklenburg-Western Pomerania (Ger_Born), in Germany displayed several mixed components, Lower Lower Saxony (Ger_Stad), and Brandenburg (Ger_Jüt) clustered, Saxony buffalos (G_Stad) had a higher proportion of Italian while Brandenburg formed a separate branch. In the second Mediterranean background compared to herds in Mecklenburg- group, Murrah breeds from Pakistan and India clustered with Western Pomerania (G_Born) and Saxony (G_Wies). Saxonian Latino America, Middle East, and Bulgarian samples. Hungarian buffalos (G_Wies) showed a consistent portion of alleles from and Romanian buffalos stood out, while animals from Saxony Romania, while the comparison at K = 4 and 6 revealed that (Ger_Wies) grouped with animals from Romanian clusters. buffalos from Brandenburg (G_Jüt) mainly contained Bulgarian genetic background. Run of Homozygosity The results of the heterozygosity metric ranged from In order to make better use of all information and to learn Bulgaria to Hungary (Hung_Csak, Hung_Tisz) with 0.41–0.35 about the demographic past and/or recent history of the buffalo for observed values (Ho) and 0.40–0.35/0.32 for expected populations of Central and Southeast Europe, we performed values (He), respectively. Similar heterozygosity values were an ROH analysis based on 285 samples genotyped in this found in the Hungarian and Romanian samples, lowest in study (Table 1), for which a panel of 61,813 autosomal SNPs relation to the other population groups (from 0.35 to 0.38 was available (Figure 6). We found 10,797 ROH regions Ho and from 0.32 to 0.36 He), while German buffalos had distributed in four categories (0–5, 5–15, 15–30, >30 Mb) with the highest values of diversity, after the Bulgarian buffalo patterns different for each population. Among German buffalos, population, between groups, and within its group (Table 3). Ger_Born had the highest average of ROH in 0–5 Mb size To support the interpretation of the Admixture and MDS category, indicative of ancient inbreeding, which decreased in all analyses, we also measured the population differentiation due the other categories. Ger_Jüt showed the opposite trend, with the to genetic structure by an FST pairwise distance analysis of highest average in the long ROH category (>30 Mb), suggestive 36,014 SNPs in all individuals of 22 populations distributed of the small effective population size and more recent inbreeding. Frontiers in Genetics | www.frontiersin.org 8 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany FIGURE 7 | Violin Plot reporting the proportion of autosome covered in run of homozygosity (ROH) aka Genomic inbreeding coefficients (FROH ) in Buffalo populations from Bulgaria (Bul_Noce), Romania (Rom_Mera, Rom_Serc), Hungary (Hun_Csak, Hun_Tisz, Hun_HT), and Germany (Ger_Born, Ger_Jüt, Ger_Stad, Ger_Wies). The black line indicates the median of FROH within the breed. Abbreviations and sample numbers are specified in Figure 3. The highest average of both short and long ROH categories was Figure 6, where both Romanian buffalo populations presented found in Ger_Stad. Buffalos from Saxony (Ger_Wies) showed a different trend of ROH categories. Hungarian buffalos had an average of ROH quite similar in all categories except a high ROH in their genome and high inbreeding due to a decrease in the 5–15 Mb category. Figure 7 displays an average small population size reflected in the different values within the inbreeding coefficient calculated from the individual sum of farms Csákvar (Hun_Csak) and Tiszataj (Hun_Tisz) with the ROH tracts per population, and in Table 4, the corresponding highest FROH (0.144 and 0.173), Földes (Hun_HT) with lower values are reported. The Bulgarian buffalo population had the FROH (0.121). The different values of FROH obtained within the lowest average of homozygosity in the genome (116.93 Mb) and, German group were an indication of the high variability between consequently, the lowest level of inbreeding (FROH = 0.047). farms (Table 4). Romanian samples from Mera village had 185.31 Mb of the genome in ROH with an FROH of 0.074, lower than the second farm (Rom_Serc) with 252.31 Mb and FROH 0.102. The buffalo Localizing ROH Hotspots and Gene Set population in Mera (Rom_Mera) is distributed in smaller herds Enrichment Analyses from multiple private farms. The second Romanian population In order to detect the ROH hotspots, Manhattan plots were (Rom_Serc) studied here is located in Sercaia and is managed by built to identify genomic regions frequently associated with ROH the Institute of Research and Development for Buffalo Breeding. (Figure 8) and the SNP locus with the highest frequency (%) Buffalos from Sercaia were characterized by ancient inbreeding, in those ROH (Table 5) for each German buffalo population. showing the highest proportion of ROH in the category 5–15 Mb. The thresholds to define the regions containing the top 1% These values are in agreement with the results presented in of most common ROH-associated SNPs (ROH hotspots or Frontiers in Genetics | www.frontiersin.org 9 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany TABLE 4 | Averages of individuals ROH in the genome and coefficient of FROH per population. Populations ID Geographic average sum of average FROH location ROH (Mb) Bul_Noce Bulgary 116.93 0.047 Ger_Born Born 192.15 0.078 Ger_Jüt Jüterbog 270.62 0.109 Ger_Stad Stadland 272.96 0.110 Ger_Wies Wiesenburg 146.96 0.059 Hun_Csak Csákvar 356.91 0.144 Hun_Tisz Tiszataj 427.16 0.173 Hun_HT Földes 299.01 0.121 Rom_Mera Mera 182.45 0.074 Rom_Serc Sercaia 252.31 0.102 ROH, run of homozygosity; FROH , inbreeding coefficient; Mb, Megabases; all other abbreviations are explained in Table 1. ROH islands) were 0.38, 0.37, 0.42, and 0.29 for Mecklenburg- Western Pomerania, Brandenburg, Lower Saxony, and Saxony, respectively. With their genomic coordinates and after functional annotation (Supplementary Table 3), we were able to reveal potential candidate genes under directional selection among the different populations (Table 6 and Figure 9). We then examined genes co-localized with the ROH hotspots to identify candidate genes present in potential genomic regions that have undergone selection in each population, and where no significant (P ≤ 0.05) enrichment was found. All the hotspot SNPs in each population were in the short ROH length, likely suggestive of ancient inbreeding or selection events (Table 5). The populations from Brandenburg (Ger_jüt) and Lower Saxony (Ger_Stad) had similar levels of inbreeding, as also found for the populations from Mecklenburg-Western Pomerania (Ger_Born) and Saxony (Ger_Wies; Table 4). The distribution of hotspot SNPs were different in every population (Figures 8A–D), indicative of no directional selection in the groups. The highest number of genes co-localized with ROH regions were found in Saxony (83) and Lower Saxony (78), followed by Brandenburg (56) and Mecklenburg-Western Pomerania (53) (Tables 5, 6 and Figure 9) involving 11–16 biological processes (Figure 10). In particular, only in buffalos from Saxony we identified genes essential for the immune system (Il18bp, Rhog), and mammalian autophagy FIGURE 8 | Genomic distribution of ROH islands in German buffalo samples. (Atg16l2) in chromosome 16 with 26, 21, and 15 hotspots The x-axis represents the SNP genomic coordinate chromosome-wise, and SNPs, respectively. the y-axis visualizes the frequency (%) of overlapping ROH shared among individuals. The red line indicates the threshold to define the regions containing the top 1% of most common ROH-associated SNPs (ROH hotspots or ROH islands). Populations are presented as DISCUSSION (A) Mecklenburg-Western Pomerania (Ger_Born), (B) Brandenburg (Ger_Jüt), (C) Lower Saxony (Ger_Stad), and (D) Saxony (Ger_Wies). Abbreviations and European Buffalos Stand Out From sample numbers are specified in Figure 3. Others This is the first study aimed to explore the genetic background of German buffalos in relation to other European populations of Indo-Pakistani background in the genome of Mediterranean (Colli et al., 2018a; Deng et al., 2019a). Our knowledge about the buffalo (Colli et al., 2018b). Here, we compared 22 populations genetic relationships of European breeds, especially in Central worldwide. The separation between Murrah (India-Pakistan, Europe and Germany are still poor. A recent study from Latino America, and the Middle East) and Mediterranean Colli et al. (2018b) reported a high genetic variability within (European populations) breeds was clearly evident. European Eastern European populations and identified recognizable traces buffalo populations are characterized by higher genetic diversity Frontiers in Genetics | www.frontiersin.org 10 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany TABLE 5 | Number of significant SNPs in ROH regions, genomic distribution, and of the Italian Mediterranean, Hungarian, and Bulgarian genetic genes mapped in the same regions for each German buffalo population analyzed background, while buffalos from the National breeding facility in this work (Chr: chomosome). (Rom_Serc) have a defined genotype. In this national breed, a Populations Chr with SNPs (total Genes (total ROH size higher level of FROH was found than in Mera buffalos. Moreover, hotspot SNPs number) number) range the Mera herds (0.074) showed the highest proportion of long size ROH that confirms the small population size but is also Ger_Born 2, 3, 5, 10, 17 679 53 0–1.4 Mb indicative of recent inbreeding events. It may be due to the Ger_Jüt 1, 3, 4, 5, 7, 949 56 0–1.9 Mb 10, 13, 14, 16, application of selection programs and the sign of environmental 20, 23 adaptation of the Mediterranean Carpathian buffalo in Rom_Serc Ger_Stad 1-5, 12, 15 899 78 0–2.5 Mb population, while in Mera the animals are less subjected to Ger_Wies 1,3, 5,6, 12, 16 635 83 0–1.4 Mb predefined or joint selection programs. According to Karpati (1997), Mediterranean buffalos were introduced to Hungary, by the Turks, during the 16th century. The animals are distributed than the Murrah breeds, which is related to their shorter distance in small private farms used for family husbandry but mainly kept from the place of domestication. This could be a consequence in national parks as a protected reserve (Borghese, 2005). In the of the natural buffalo repopulation in Europe after the Ice Age MDS plot, the Hungarian cluster, including three populations or by the distribution of domesticated buffalo by different routes (Földes, Csákvar, and Tiszataj), showed a lower admixture (Bartosiewicz, 1999). Admixture analysis revealed that European with Balkan populations (Figure 4B). The comparably high and non-European populations do not share the same gene pool. level of inbreeding coefficient and a high proportion of ROH, Only Bulgarian Murrah buffalos were profoundly different from especially in short length, is indicative of ancient inbreeding the other European animals. In K = 4, a large amount of alleles events in Hungarian buffalos. Notably, the genetic background were common with Murrah breeds from India-Pakistan but not in Hungarian buffalo is comparably homogenous compared in K = 6 (Figure 3B), indicating the development of a different to all genotypes published before or presented here and may genetic background. In fact, the values of observed and expected represent a valuable resource for local economy or breeding heterozygosity (He = 0.40) were comparable to those previously programs worldwide in the future. Clearly, and maybe as reported by Colli et al. (2018b) in Murrah breeds from India- a surprise, the genetic background in Hungarian buffalo is Pakistan, underlining their similarity and indicating an isolate- different from the population in Mera, which represents a former breaking effect, the mixing of previously isolated populations. Hungarian village. In all three Hungarian farms included in From the data provided here, we have strong evidence that this study, the buffalos are used for meat production only European buffalos stand out from other buffalo populations and serve particularly for the establishment of local marketing worldwide. We may thus assume ancient origins of European together with other regional products. Under current natural buffalos (excluding Bulgaria), supporting the idea of Krawczynski breeding, there is no selection for milk production at present (2010) postulating that buffalos did not wholly disappear from in the herds included, which may explain the lack of recent the European continent during the late Pleistocene. inbreeding. Bulgarian Murrah buffalos, instead, had the lowest proportion of ROH in each length category and the lowest FROH values between Central and East European populations. Different Genotypes in European Indeed, the Bulgarian is the biggest population analyzed in Buffalos and a Novel One in Hungary this study, and the samples were collected from two different Four clear clusters could be distinguished in the MDS locations, which also could be related to comparably low levels plot when including samples from five European countries of inbreeding. In the German populations, genomic relationships only, representing the geographical distribution of Bulgarian, with Bulgarian Murrah, Italian Mediterranean breeds, and Hungarian, Romanian, and the Italian population. German Romanian genetics were identified. Buffalos in Germany had populations, however, were mixed and located in between the highest observed and expected heterozygosity values within the Italian and Romanian groups. In the Romanian cluster, all Central Eastern European populations studied here. The samples collected from the National breeding facility (Rom_Serc) observed genetic diversity is a significant factor for their potential were separated from the other farm located in a former adaptation to the local environment and provides a valuable Hungarian village (Mera), now part of Romania. The origin substrate for an upcoming selection program. Drift rather than of Romanian buffalo is still unclear. According to Borghese, it mutation has likely created this variability due to the short period is a Mediterranean type adapted to the cold climate and local in generations and the relatively small effective population sizes environment and is classified as a Mediterranean Carpathian (Wang, 2005). Buffalos from Brandenburg (Ger_Jüt) stood out breed (Borghese, 2005). Because of the massive decrease of their from all other German populations studied here. As the only population size during the last 20 years, they are nowadays in an German population, clear genetic relations to Bulgarian buffalos endangered status, distributed in small farms in some villages, can be postulated for the Jüterbog herd. In fact, part of the and almost 98% in Transylvania used for local food production Bulgarian Murrah’s gene pool that was visible in the admixture and draft strength (Matiuti et al., 2020). Our results revealed at K = 4, and disappeared at K = 6, is likely an effect of genetic that Mera buffalos (Rom_Mera), which were distributed in drift after isolation over time. The Brandenburg buffalo herd multiple smaller familiar farms, were characterized by a mixture showed the highest deviation of heterozygosity values, with He Frontiers in Genetics | www.frontiersin.org 11 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany TABLE 6 | Shared and unique genes overlapping ROH regions with the highest frequency SNPs in the German populations. Populations N. of Genes Name of Genes Ger_Born,Ger_Jüt,Ger_Stad 4 RAB38, GRM5, TMEM135, CTSC Ger_Born, Ger_Jüt, Ger_Wies 3 CA10, MFSD14B, ANXA10 Ger_Born, Ger_Stad 6 TOM1L1, MMD, LOC102416224, LOC102390200, HLF, STXBP4 Ger_Born, Ger_Wies 14 TAF1D, LOC112585128, PANX1, DEUP1, ELAVL2, SLC36A4, HEPHL1, MED17, LOC112585257, FAT3, LOC102389890, C5H11orf54, VSTM5, CEP295 Ger_Jüt, Ger_Stad 2 ME3, PRSS23 Ger_Jüt, Ger_Wies 3 NLGN1, KIF2B, NAALADL2 Ger_Stad, Ger_Wies 14 GPR17, LOC102408414, HS6ST1, AMMECR1L, PLEKHB2, UGGT1, IMP4, PTPN18, ERCC3, AMER3, LRRTM4, LOC112581251, TUBGCP5, BIN1 Ger_Born 26 FUT9, GALNTL6, MYO1B, FSTL5, LOC102389548, MFAP3L, FHL5, CAVIN2, NDUFAF4, MMS22L, MTNR1B, TMEFF2, LOC102389468, UFL1, PALLD, AADAT, CBR4, NEK1, MANEA, LOC102409600, SH3RF1, GPR63, KLHL32, CLCN3, NABP1, TMEM100 Ger_Jüt 44 NSMCE3, FAM189A1, ZNF25, ZPLD1, DIAPH3, ZNF37A, BMS1, LOC102392518, LOC112581575, TJP1, BTBD3, SIM1, RASGEF1A, LOC112581582, APBA2, FXYD4, FZD4, TASP1, ZNF248, LOC102406297, C23H10orf143, SPTLC3, GRIK2, MCHR2, GLRX3, MGMT, LOC102407296, ASCC3, LOC102400568, TARSL2, BICC1, SNRPA1, OTUD7A, CSGALNACT2, RET, LOC112578510, EBF3, PCSK6, HNRNPF, CHRM3, LOC102395743, ISM1, TM2D3, LOC102397911 Ger_Stad 52 SLC6A15, CPA6, MYBL1, AHSA2P, TYR, CCT4, VPS54, LOC102392047, PPP1R42, C15H8orf34, MGAT4C, LRRIQ1, CSMD1, SULF1, DNAJC5B, DLG2, VCPIP1, NTS, PREX2, TCF24, FAM161A, CSPP1, ARFGEF1, PDE7A, MYOM2, ARHGEF10, COPS5, WDPCP, CLN8, COMMD1, VXN, CRH, UGP2, OCA2, EHBP1, MCMDC2, ADHFE1, TSPAN19, TRIM55, KBTBD11, LOC112577753, DLGAP2, OTX1, MDH1, ALX1, B3GNT2, RRS1, NIPA1/2, ERICH1, TMEM17, RASSF9, LOC102407794 Ger_Wies 49 NIPA1, CLPB, NAB1, ART1, PGAP2, LOC102413639, CADM2, CHRNA10, NEMP2, ARHGEF4, LAMTOR1, FCHSD2, NUP98, RNF121, STIM1, PDE2A, SAP130, MFSD6, LOC112577983, ARHGEF17, PHOX2A, LOC102407599, ATG16L2, LIMS2, IWS1, WDR33, MYO7B, RELT, CCDC115, STARD10, ARAP1, SFT2D3, P2RY6, TRIM21, LOC102390798, FAM168B, P2RY2, POLR2D, RHOG, IL18BP, FOLR2/INPPL1, RRM1, CYFIP1, MAP3K2, TRIM68, PROC, ANAPC15, NUMA1, IZUMO3 being lower than Ho, suggesting Hardy-Weinberg disequilibrium in 11 biological processes. While Saxony buffalos (Ger_Wies), and inbreeding (Mayo, 2008). The higher values of inbreeding the population with the lowest inbreeding level and the lowest with an increasing trend of ROH average from short to long size number of hotspot SNPs (635), showed their distribution may further suggest recent inbreeding and isolation without any limited to 6 chromosomes but to 83 genes involved in 15 admixture since the introduction of the buffalos in Germany. biological processes. Overlapping ROH signals affecting genes involved in the immune system (Il18bp, Rhog) and autophagy (Atg16l2) were only found in Saxony. It is reported that Overlapping ROH Signals in German the gene Atg16l2 (autophagy related 16 like 2) influences Buffalo Genome the adaptation of the immune system to the recovery from The control of inbreeding in a population is necessary to mastitis in Danish Holstein cattle (Welderufael et al., 2018). avoid depression of phenotype traits (Ouborg et al., 2010; From these very initial interpretations of the potential effects Fontanesi et al., 2015). Considering that ROH regions are not of local selection for functional traits, we may get an idea randomly present in the genome (McQuillan et al., 2008), the of how genomic selection could be used to improve animal detection of the genomic regions frequently covered with ROH health and performance in the future. Certainly, additional is useful. This helps to understand whether they influence studies of gene expression, metabolism, or animal health are any quantitative trait, especially in a population of small needed to link SNP markers with particular phenotypes. We population size, such as local breeds (Biscarini et al., 2015). are aware of the limitations of our study. Only four different The negative effect of inbreeding depression in milk production farms in Germany are included that represent less than 15% has been reported in Egyptian (Khattab and Kawthar, 2007), of the German buffalo population. Moreover, the selection South Iranian (Mirhabibi et al., 2007), and Brazilian buffalo of farms in Eastern Europe was based on potential familiar populations (Santana et al., 2011). We found several differences relationships and can be seen as a start only. For German among German buffalo populations in ROH analysis in terms and other European countries, higher coverage, including of the number of hotspot SNPs, candidate genes in the same more herds and additional countries to further exploit the genomic coordinates, and biological processes in which they full setting of genetic diversity in Europe, and an in-depth could be involved. In the Brandenburg (Ger_Jüt) population, analysis of signatures of past selection combining several metrics the most distinguished, with a high level of inbreeding, we together with ROH is required. Finally, in the genotyping found the highest number of hotspot SNPs (949), distributed array designed for buffalo, only 30% of SNPs from European over 11 chromosomes and concerned only 57 genes involved breeds (Italian Mediterranean) have been included, which could Frontiers in Genetics | www.frontiersin.org 12 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany FIGURE 9 | Venn diagram reporting the number of common and unique genes present in overlapping ROH regions with the highest frequency SNPs in different buffalo populations from Germany (gene names are referred by Table 6). FIGURE 10 | Classification of genes with hotspot SNPs in ROH genome regions. Populations are presented as (A) Mecklenburg-Western Pomerania (Ger_Born), (B) Brandenburg (Ger_Jüt), (C) Lower Saxony (Ger_Stad), (D) Saxony (Ger_Wies). Abbreviations and sample numbers are specified in Figure 3. Frontiers in Genetics | www.frontiersin.org 13 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany be a further limitation to obtain genomic information about are not reviewed by ethical committees. Written informed European buffalo populations. Accordingly, we may initiate consent was obtained from the owners for the participation of new genome sequencing initiatives to develop more informative their animals in this study. SNP arrays to improve further functional genome analysis and genome-based selection in European buffalo. AUTHOR CONTRIBUTIONS CONCLUSION AN, KW, and AH conceived the project and designed the experiments. AN performed the laboratory part. AN, SQ, RG-P, This work aimed to contribute to an improved understanding and ML-S performed the data analysis. AN, SQ, RG-B, JB, RK, of Central and Southeast European buffalo populations from and AH contributed to the interpretation of results and wrote a genomic perspective. We defined the genetic characteristics the manuscript. MiT, M-AF, HP, AB, LV, CH, LH, PP, YI, TP, of European buffalo populations in Romania, Bulgaria, and WL, MaT, and AH collected or provided samples and contributed Hungary. Together with the Italian Mediterranean breed, these to data production. All authors read, made corrections, and populations contribute to the high degree of genetic diversity of approved the final version of the manuscript. European buffalos. Incidentally, we may have genetic evidence for a novel or original Hungarian breed, characterized by outmost uniformity in three different farms with no clear relations to ACKNOWLEDGMENTS any other genotype described so far. Since all farms are in the close neighborhood, follow-on studies may help to unravel the This study would not have been possible without the help of many origin of the Hungarian population. This novel genetic identity buffalo breeders. We express our gratitude to Family Henrion found in Hungarian buffalos thus adds one additional specific (Bobalis Agrargesellschaft mbH in Jüterbog), Dr. Möhring breed to the existing genetic diversity in Europe, which so far (Gut Darss GmbH & Co. KG in Born), the staff in Stadland consisted only of the Italian Mediterranean, Bulgarian Murrah, (Hof am Meer, Stadland/Norderschwei), and all helpers from and Romanian genetics. Wiesenburger Land eG in Wiesenburg. We further want to This is the first genetic characterization of buffalos in thank the managers and helpers in Hungary: Mr. Mihály Bodnár, Germany, where farming started only about 40 years ago and Tiszatáj Fundation (H- 4450 Tiszalök Rákóczi u. 14), Mr. Levente therefore is in its infancy. We have identified familiar relations Viszló Pro Vértes Fundation (H-8083 Csákvár, Kenderesi u. to all genotypes identified in European buffalos so far with Geszner ház), and Mrs. Gabriella Bona Bihar Fundation (H- the exception of Hungary. Importantly and based only on a 4177, Földes Kálóháti Tanya), the staff in Sercaia (Romania) and comparably small number of herds, a diverse genetic make-up is all the families in Mera (Romania) that supported our project. available in European and German buffalos, which can be seen Furthermore, the local support in Bulgaria is thanked here. as an excellent prerequisite for the development of Buffalo-based Finally, we would like to thank Dr. Taina Figueiredo Cardoso for bioeconomy in Europe based on local breeds characterized by helpful advice when performing genetic diversity analyses. unique features. SUPPLEMENTARY MATERIAL DATA AVAILABILITY STATEMENT The Supplementary Material for this article can be found Our data set has been uploaded in the Dryad Repository, online at: https://www.frontiersin.org/articles/10.3389/fgene. https://datadryad.org/stash/share/eT_rs3441UHCMlN4qpbRi8u 2020.610353/full#supplementary-material 6sFPHdX-ZtDPuVS5TpSI and https://datadryad.org/resource/ doi:10.5061/dryad.h0cc7. Supplementary Table 1 | A descriptive summary of three SNP panels used in this study. Supplementary Table 2 | Wright’s Fixation Index (FST) matrix computed with ETHICS STATEMENT 36,014 SNPs in 22 buffalo populations. Supplementary Table 3 | List of potential candidate genes and relative functional Ethical review and approval was not required for the animal annotation found in ROH hotspots region in German buffalo populations: (a) study because sampling occurred as part of livestock breeding Mecklenburg-Western Pomerania, (b) Brandenburg, (c) Lower programs which are granted by international law and therefore Saxony, (d) Saxony. REFERENCES Alexander, D. H., Novembre, J., and Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664. doi: 10.1101/ Abd El-Salam, M. H, and El-Shibiny, S. (2011). A comprehensive review on the gr.094052.109 composition and propertiesof buffalo milk. Dairy Sci. Technol. 91, 663–699. Alexiev, A. (1998). The Water Buffalo. Bulgaria: St. Kliment Ohridski University doi: 10.1007/s13594-011-0029-2 Press. Frontiers in Genetics | www.frontiersin.org 14 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany Bartosiewicz, L. (1999). “The emergence of holocene faunas in the carpathian basin: Iannuzzi, L., and Di Meo, G. (2009). “Water buffalo,” in Genome Mapping and a review,” in The Holocene History of the European Vertebrate Fauna, Ed. N. Genomics in Domestic Animals, eds N. E. Cockett, and C. Kole, (Germany: Benecke, (Rahden: Modern Aspects of Research) 6, 73–90. Springer-Verlag). Biscarini, F., Cozzi, P., Gaspa, G., and Marras, G. (2018). detectRUNS: Detect Karpati, L. (1997). Bufaloes in hungary. Bufalo Newlett. 7:9. Runs of Homozygosity and Runs of Heterozygosity in Diploid Genomes. Kawashima, T., Sumamal, W., Pholsen, P., Chaithiang, R., and Kurihara, M. (2006). Available at: http://orca.cf.ac.uk/id/eprint/108906 (accessed March 16, Comparative study on energy and nitrogen metabolisms between Brahman 2020). cattle and swamp buffalo fed with low quality diet. JARQ 40, 183–188. doi: Biscarini, F., Nicolazzi, E. L., Stella, A., Boettcher, P. J., and Gandini, G. (2015). 10.6090/jarq.40.183 Challenges and opportunities in genetic improvement of local livestock breeds. Khattab, A. S., and Kawthar, A. M. (2007). Inbreeding and it is effects on some Front. Genet. 6:33. doi: 10.3389/fgene.2015.00033 productive and reproductive traits in a herd of Egyptian buffalos. Ital. J. Anim. Bökönyi, S. (1957). Mitteleuropäische vorgeschichtliche Büffelfunde. Folia Sci. 6, 275–278. doi: 10.4081/ijas.2007.s2.275 Archaeol. 9, 39–35. Kierstein, G., Vallinoto, M., Silva, A., Schneider, M. P., Iannuzzi, L., and Brenig, B. Borghese, A. (2005). Buffalo production and research. FAO Ed. REU Tech. Ser. 67, (2004). Analysis of mitochondrial D-loop region casts new light on domestic 1–315. water buffalo (Bubalus bubalis) phylogeny. Mol. Phylogenet. Evol. 30, 308–324. Borghese, A. (2011). Situation and perspectives of buffalo in the World. Europe and doi: 10.1016/s1055-7903(03)00221-5 Macedonia. Mac. J. Anim. Sci. 1, 281–296. Krawczynski, R. (2010). “Zur historischen verbreitung des wasserbüffel (Bubalus Borghese, A. and Moioli, B. (2016). Buffalo Breeds and Management Systems. Bubalis L., 1758) in Eropa,” in Wasserbüffel in Der Landschaftspflege ed. J. Northern Greece: Animal Production Research Institute. Hoffman, R. Krawczynski, and H. G. Wagner, (Lexxion: Tagungsband), 17–26. Calabró, S., Infascelli, F., Tudisco, R., Musco, N., Grossi, M., Monastra, G., et al. doi: 10.1111/j.1439-0507.1988.tb03852.x (2013). Estimation of in vitro methane production in buffalo and cow. Buffalo Krawczynski, R. (2013). “Die potentiell natürliche megafauna Europas,” in Bull. 32, 924–927. Nationalpark-Jahrbuch Unteres Odertal 9 Nationalparkstiftung Unteres Odertal, Clutton-Brock, J. (1999). A Natural History of Domesticated Mammals. Cambridge: ed. A. Vössing, (Schwedt: Schloss Criewen), 29–40. Cambridge University Press, 238. Kumar, S., Nagarajan, M., Sandhu, J. S., Kumar, N., Behl, V., and Nishanth, G. Cockrill, W. R. (1984). “Water buffalo,” in Evolution of Domesticated Animals, ed. (2007). Mitochondrial DNA analyses of Indian water buffalo support a distinct I. L. Mason, (London: Longman), 52–63. genetic origin of river and swamp buffalo. Anim. Genet. 38, 227–232. doi: Colli, L., Milanesi, M., Vajana, E., Iamartino, D., Bomba, L., Puglisi, F., 10.1111/j.1365-2052.2007.01602.x et al. (2018b). New insights on water buffalo genomic diversity and post- Lau, C. H., Drinkwater, R. D., Yusoff, K., Tan, S. G., Hetzel, D. J., and domestication migration routes from medium density SNP chip data. Front. Barker, J. S. (1998). Genetic diversity of Asian water buffalo (Bubalus Genet. 9:53. doi: 10.3389/fgene.2018.00053 bubalis): mitochondrial DNA D-loop and cytochrome b sequence Colli, L., Milanesi, M., Vajana, E., Iamartino, D., Bomba, L., Puglisi, F., et al. variation. Anim. Genet. 29, 253–264. doi: 10.1046/j.1365-2052.1998.00 (2018a). Data From: New insights on Water Buffalo Genomic Diversity and 309.x Post-Domestication Migration Routes from Medium Density SNP Chip Data, Lima-Ribeiro, M. S., and Diniz-Filho, J. A. F. (2017). Climate change, Philadelphia, PA: Dryad. doi: 10.5061/dryad.h0cc7 human overkill, and the extinction of megafauna: a macroecological Deng, T., Aixin, L., Jiajia, L., Guohua, H., Tingzh, Y., Shenhe, L., et al. (2019a). Data approach based on pattern-oriented modelling. Evol. Ecol. Res. 18, From: Genome-wide SNP Data Revealed the Extent of Linkage Disequilibrium, 97–121. Persistence of Phase and Effective Population Size in Purebred and Crossbred Low, W. Y., Tearle, R., Bickhart, D. M., Rosen, B. D., Kingan, S. B., Swale, T., et al. Buffalo Populations, Philadelphia, PA: Dryad. doi: 10.5061/dryad.310pf0 (2019). Chromosome-level assembly of the water buffalo genome surpasses 5 human and goat genomes in sequence contiguity. Nat. Commun. 10:260. Deng, T., Liang, A., Liu, J., Hua, G., Ye, T., Liu, S., et al. (2019b). Genome-Wide Martin, P. S. (1984). “Prehistoric overkill: the global model,” in Quaternary SNP data revealed the extent of linkage disequilibrium, persistence of phase and Extinctions: A Prehistoric Revolution, eds P. S. Martin, and R. G. Klein, (Tucson: effective population size in purebred and crossbred buffalo populations. Front. University of Arizona Press), 354–403. Genet. 9:688. doi: 10.3389/fgene.2018.00688 Martini, M., Altomonte, I., and Salari, F. (2018). Native milk fat globule size and its Excoffier, L., and Lischer, H. E. L. (2010). Arlequin suite ver 3.5: a new influence on the natural creaming properties of buffalo milk. Buffalo Bull. 37, series of programs to perform population genetics analyses under Linux and 481–488. Windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010. Mastrangelo, S., Sardina, M. T., Tolone, M., Di Gerlando, R., Sutera, A. M., 02847.x Fontanesi, L., et al. (2018). Genome-wide identification of runs of homozygosity Ferrara, B. (1964). The present situation of buffalo breeding in Italy. Acta Med. Vet. islands and associated genes in local dairy cattle breeds. Animal 12, 2480–2488. Napoli. 10, 325–355. doi: 10.1017/s1751731118000629 Fontanesi, L., Scotti, E., Samore, A. B., Bagnato, A., and Russo, V. (2015). Matiuti, M., Matiuti, C. L., Garleac, C., and Hutu, I. (2020). Particularities of buffalo Association of 20 candidate gene markers with milk production breeding in Romania. One Health Int. J. 6, 1–5. and composition traits in sires of Reggiana breed, a local dairy Maymone, B. (1942). The Buffalo in Italy. Ann. Ist. Sper. Zoot. 3, 5–66. cattle population. Livest. Sci. 176, 14–21. doi: 10.1016/j.livsci.2015. Mayo, O. (2008). A century of hardy–weinberg equilibrium. Twin Res. Hum. Genet. 03.022 11, 249–256. doi: 10.1375/twin.11.3.249 Haarmann, H. (2011). Das Rätsel der Donauzivilisation: Die Entdeckung der ältesten McQuillan, R., Leutenegger, A. L., Abdel-Rahman, R., Franklin, C. S., Pericic, M., Hochkultur Europas. Munich: C. H. Beck Verlag, 1–286. Barac-Lauc, L., et al. (2008). Runs of homozygosity in European populations. Hegde, N. G. (2019). Buffalo husbandry for sustainable development of small Am. J. Hum. Genet. 83, 359–372. farmers in india and other developing countries. Asian J. Res. Anim. Vet. Sci. Michelizzi, V. N., Dodson, M. V., Pan, Z., Amaral, M. E., Michal, J. J., McLean, 3, 1–20. D. J., et al. (2010). Water buffalo genome science comes of age. Int. J. Biol. Sci. Hoffman, J. M., and Valencak, T. G. (2020). A short life on the farm: aging and 6, 333–349. doi: 10.7150/ijbs.6.333 longevity in agricultural, large-bodied mammals. Geroscience 42, 909–922. doi: Mirhabibi, S., Manafiazar, G., Qaravisi, S., and Mahmoodi, B. (2007). Inbreeding 10.1007/s11357-020-00190-4 and its effect on some productive traits in buffalos of South Iran. Ital. J. Anim. Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and Sci. 6, 372–374. doi: 10.4081/ijas.2007.s2.372 integrative analysis of large gene lists using DAVID bioinformatics resources. Moss, A. R., Jouany, J. P., and Newbold, J. (2000). Methane production by Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211 ruminants: its contribution to global warming. Ann. De Zootech. 49, 231–253. Iamartino, D., Nicolazzi, E. L., Van Tassell, C. P., Reecy, J. M., Fritz-Waters, E. R., doi: 10.1051/animres:2000119 Koltes, J. E., et al. (2017). Design and validation of a 90K SNP genotyping assay Nicolazzi, E. L., Iamartino, D., and Williams, J. L. (2014). AffyPipe: an open- for the water buffalo (Bubalus bubalis). PLoS One 12:e0185220. doi: 10.1371/ source pipeline for Affymetrix Axiom genotyping workflow. Bioinformatics 30, journal.pone.0185220 3118–3119. doi: 10.1093/bioinformatics/btu486 Frontiers in Genetics | www.frontiersin.org 15 January 2021 | Volume 11 | Article 610353 Noce et al. Genetic Diversity of Bubalus b. in Germany Ouborg, N. J., Pertoldi, C., Loeschcke, V., Bijlsma, R., and Hedrick, P. W. (2010). from mastitis in danish holstein cows. Front. Genet. 9:141. doi: 10.3389/fgene. Conservation genetics in transition to conservation genomics. Trends Genet. 26, 2018.00141 177–187. doi: 10.1016/j.tig.2010.01.001 Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. New York, NY: Paradis, E., and Schliep, K. (2019). ape 5.0: an environment for modern Springer-Verlag. phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528. doi: Wright, S. (1965). The interpretation of population-structure by F-statistics with 10.1093/bioinformatics/bty633 special regard to systems of mating. Evolution 19, 395–420. doi: 10.1111/j.1558- Pucher, E. V. (1988). Erstnachweis des Europäischen Wildesels (Equus hydruntinus 5646.1965.tb01731.x Regalia, 1907) im Holozän Österreichs. Ann. Naturhist. Mus. Wien 92/B, 31–48. Zahariev, Z. I., Aleksiyev, A., and Nikolova, S. D. (1986). Morphological and Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A., Genetics Survey of Water Buffalos. Bivali: Zemizdat. Bender, D., et al. (2007). PLINK: a tool set for whole-genome association Zeuner, F. E. (1963). A History of Domesticated Animals. New York, NY: and population-based linkage analyses. Am. J. Human Genet. 81, Hutchinson; Harper & Row. 559–575. Zicarelli, L. (2001). La bufala Mediterranea Italiana: esempio di una razza Purfield, D. C., McParland, S., Wall, E., and Berry, D. P. (2017). The distribution autoctona in espansione. Sci. Tecn. Latt. Cas. 52, 279–284. of runs of homozygosity and selection signatures in six commercial meat sheep Zicarelli, L. (2004). Buffalo milk: its properties, dairy yield and mozzarella breeds. PLoS One 12:e0176780. doi: 10.1371/journal.pone.0176780 production. Vet. Res. Commun. 28(Suppl. 1), 127–135. doi: 10.1023/b:verc. Roth, J., and Myers, P. (2004). Bubalus Bubalis. Avaliable online at: 0000045390.81982.4d http://animaldiversity.ummz.umich.edu/site/accounts/information/Bubalus_ bubalis.html (accessed April 20, 2012). Conflict of Interest: M-AF, MiT, HP, and RK are employed by the companies Santana, M. L. Jr., Aspilcueta-Borquis, R. R., Bignardi, A. B., Albuquerque, L. G., Gut Darss GmbH & Co., Hof am Meer GmBH, Wiesenburger Land eG, and and Tonhati, H. (2011). Population structure and effects of inbreeding on Energiequelle GmbH, respectively. milk yield and quality of Murrah buffalos. J. Dairy Sci. 94, 5204–5211. doi: 10.3168/jds.2011-4377 The remaining authors declare that the research was conducted in the absence of Thomas, P. D., Campbell, M. J., Kejariwal, A., Mi, H. Y., Karlak, B., Daverman, any commercial or financial relationships that could be construed as a potential R., et al. (2003). PANTHER: a library of protein families and subfamilies conflict of interest. indexed by function. Genome Res. 13, 2129–2141. doi: 10.1101/gr.77 2403 Copyright © 2021 Noce, Qanbari, González-Prendes, Brenmoehl, Luigi-Sierra, Turner, S. D. (2014). qqman: an R package for visualizing GWAS results using Q-Q Theerkorn, Fiege, Pilz, Bota, Vidu, Horwath, Haraszthy, Penchev, Ilieva, Peeva, and manhattan plots. bioRxiv [Preprint] doi: 10.1101/005165 [Preprint], Lüpcke, Krawczynski, Wimmers, Thiele and Hoeflich. This is an open-access article Wang, J. (2005). Estimation of effective population sizes from data on genetic distributed under the terms of the Creative Commons Attribution License (CC BY). markers. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 360, 1395–1409. doi: The use, distribution or reproduction in other forums is permitted, provided the 10.1098/rstb.2005.1682 original author(s) and the copyright owner(s) are credited and that the original Welderufael, B. G., Løvendahl, P., de Koning, D.-J., Janss, L. L. G., and Fikse, W. F. publication in this journal is cited, in accordance with accepted academic practice. No (2018). Genome-wide association study for susceptibility to and recoverability use, distribution or reproduction is permitted which does not comply with these terms. Frontiers in Genetics | www.frontiersin.org 16 January 2021 | Volume 11 | Article 610353