Genetic Diversity of Bubalus bubalis in Germany and Global Relations of Its Genetic Background

This is the first study to explore the genetic diversity and population structure of domestic water buffalo (Bubalus bubalis) in Germany and their potential relations to herds in other parts of Europe or worldwide. To this end, animals from different herds in Germany, Bulgaria, Romania, and Hungary were genotyped and compared to genotypes from other populations with worldwide distribution and open to the public. The pilot study analyzed population structure, phylogenetic tree, and inbreeding events in our samples. 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 was identified in European buffalos, and a novel genotype was described in Hungarian buffalos by this study. We demonstrate that European buffalos stand out from other buffalo populations worldwide, supporting the idea that buffalos have not completely disappeared from the European continent during the late Pleistocene. The high genetic diversity in European buffalos seems to be an excellent prerequisite for the establishment of local breeds characterized by unique traits and features. This study may also be considered as an initial step on the way to genome characterization for the sustainable development of the buffalo economy in Germany and other parts of Europe in the future.

This is the first study to explore the genetic diversity and population structure of domestic water buffalo (Bubalus bubalis) in Germany and their potential relations to herds in other parts of Europe or worldwide. To this end, animals from different herds in Germany, Bulgaria, Romania, and Hungary were genotyped and compared to genotypes from other populations with worldwide distribution and open to the public. The pilot study analyzed population structure, phylogenetic tree, and inbreeding events in our samples. 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 was identified in European buffalos, and a novel genotype was described in Hungarian buffalos by this study. We demonstrate that European buffalos stand out from other buffalo populations worldwide, supporting the idea that buffalos have not completely disappeared from the European continent during the late Pleistocene. The high genetic diversity in European buffalos seems to be an excellent prerequisite for the establishment of local breeds characterized by unique traits and features. This study may also be considered as an initial step on the way to genome characterization for the sustainable development of the buffalo economy in Germany and other parts of Europe in the future.

INTRODUCTION
Water buffalo (Bubalus bubalis) is a multipurpose animal, producing meat, milk, leather products, and dung. In developing countries, buffalos are used as draught animals, providing more than 40% of farm labor (Borghese, 2005). Due to their wallowing behavior, buffalos are less susceptible to ectoparasites and related diseases and suitable for grazing in the swamps (Michelizzi et al., 2010).
Buffalos are also more efficient in the digestion process due to their longer rumen and produce a lower average of methane emission (Moss et al., 2000;Kawashima et al., 2006;Calabró et al., 2013). Compared to cattle, the buffalo is less selective for feed quality (Michelizzi et al., 2010) and longer living, granting husbandry at lower costs and with less frequent animal turnover (Borghese, 2005;Hegde, 2019;Hoffman and Valencak, 2020). Buffalo meat is a rich source of proteins, fatty acids as omega 3-6 fatty acids, iron, and characterized by lower concentrations of cholesterol if compared to cattle (Roth and Myers, 2004). Milk from buffalo contains biliverdin, bioactive pentasaccharides, and gangliosides, which are not present in bovine milk (Abd El-Salam and El-Shibiny, 2011) and more fat in total but lower cholesterol compared to the cattle. Buffalo milk is further characterized by higher protein content and big casein micelles; both features are related to the higher yields of cheese produced from the same amounts of buffalo versus dairy cow milk (Zicarelli, 2004;Martini et al., 2018).
The domestic water buffalo is classified into two major categories: swamp buffalo (Bubalus bubalis carabanensis, 2n = 48) and river buffalo (Bubalus bubalis bubali, 2n = 50). While their taxonomical status is still being debated, B. bubalis is supposed to be descended from Indian wild buffalo (Bubalus arnee) domesticated ≈5,000 years ago (Cockrill, 1984). Nevertheless, the origin of the two subspecies into which water buffalo are divided is still the object of study (Colli et al., 2018b), especially to solve the debate about the occurrence of a single versus two separate domestication events (Lau et al., 1998;Kierstein et al., 2004;Kumar et al., 2007).
Little is known about the history, genetic diversity, and performance of European buffalo populations. The presence of water buffalo in Europe is dated from the Pleistocene until the warm period before the last Eem-Interglacial (≈125,000 years ago), when buffalos were also present in Central Europe (Zahariev et al., 1986;Alexiev, 1998). Although hunting (Martin, 1984) and climatic changes, particularly during the late Pleistocene (Lima-Ribeiro and Diniz-Filho, 2017), temporarily or locally displaced and shrank the European buffalo populations, the presence of refugial regions in Southern and Eastern Europe may have prevented the extinction of the Bubalus species from the European continent (Krawczynski, 2010). The presence of buffalos in South-Eastern Europe during the early Neolithikum (9,000-7,000 BC) was suggested (Bökönyi, 1957) but generally dismissed (Bartosiewicz, 1999). Based on bone finds in Austria dated to the Atlantikum (7,000-4,000 BC), other authors have also discussed the presence of buffalos during the Holocene (beginning 12,000 BC) (Pucher, 1988). However, defined criteria for the unambiguous distinction of members from the Bovini tribus (in particular, Bos versus 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), and buffalo (Bubalus sp.). This linguistic approach supports the hitherto controversial bone finds from the same area at about the same time. It is unclear whether water buffalos repopulated Europe after the ice age or whether they survived in refugial areas. Because rock reliefs, dated back to ≈16,000 BC in France, clearly display buffalos (Krawczynski, 2013). This, in fact, would be an indication that water buffalos in Europe survived the ice age.
Domesticated water buffalos in Europe are referenced during the centuries VI-XII C.E. (Maymone, 1942;Alexiev, 1998) and had been settled in Italy, Bulgaria, Romania, and all the other Balkan countries, where they have been preserved and bred for centuries until today (Ferrara, 1964). European buffalos are river type (Borghese, 2011). The Balkans (East Europe) are an essential point of the river buffalo's historical migration route. Riverine buffalo moved to Southwestern Asia from the Indian domestication area, reached Egypt and Turkey, and arrived in East Europe and Italy during the seventh century (Zeuner, 1963;Clutton-Brock, 1999). Some animals likely returned then to Egypt, Turkey, and Bulgaria with Crusaders returning, and spread into the other Balkan regions during the 12th century (Ferrara, 1964;Borghese, 2005). Since 1980 buffalos have also been introduced to Germany, France, Spain, Portugal, Luxembourg, The Netherlands, Switzerland, and other parts of Europe. Nowadays, the only officially recognized breeds in Europe are the Mediterranean Italian buffalo and the Bulgarian Murrah. Selecting animals for milk traits was initiated in Italy. Crossing with other populations was avoided; a herd book with animals, performance, and morphological traits was established to maintain its unique genetic identity (Iannuzzi and Di Meo, 2009). The primary purpose was the production and marketing of milk and milk derivatives such as mozzarella cheese, one of the "pasta filata" cheeses known worldwide (Zicarelli, 2001).
While local Bulgarian buffalos were a significant European Mediterranean population raised for draft power, meat, and milk purposes, they started to be crossed with Murrah breed since 1972, mainly to improve their performance for milk production. A selection program was initiated to develop a typical Bulgarian Murrah breed for milk with high butterfat content (Borghese, 2005;Borghese and Moioli, 2016).
In Germany, 7,614 buffalos, spread over 16 different regions, are recorded for the year 2020 by an internal communication from the German Buffalo Breeder's Association. The development of the single nucleotide polymorphism (SNP) genotyping assay specific for river buffalo (Iamartino et al., 2017), Axiom R Buffalo Genotyping Array 90K from Affymetrix, now enables genome analysis of this comparably novel farm animal species and could be useful for the initiation of genetic selection in German buffalo. A subsequent release of the first assembly of the water buffalo genome (Low et al., 2019) further enabled researchers to expand their knowledge about this on a molecular scale. In this study, we sought to explore the genetic diversity and structure of domestic water buffalo populations in Germany and Eastern Europe.

Ethics Statement
All sampling occurred as part of commercial buffalo breeding programs and strictly adhered to national and international laws.  Table 1.

Samples Collection, DNA Purification, and Genotyping
Ear tag test samples from 285 female and male buffalos were collected randomly from ten farms during the years 2018 and 2019. Importantly, all breeders were asked to provide samples from unrelated animals. The farms are distributed in Central and Eastern Europe (Figure 1), as listed in Table 1. DNeasy Blood & Tissue Kit (Qiagen, Germany) was used to extract genomic DNA according to the instructions of the manufacturer. Purification and concentration of DNA were measured with a NanoDrop2000 before the normalization to the required concentration (30 ng/µL). The quality control and genotyping using the Axiom R Buffalo Genotyping Array 90K from Affymetrix 1 were performed by ATLAS Biolabs GmbH (Berlin, Germany). Allele calling was carried out using Axiom Analysis Suite software V4.0.1 (Applied Biosystems by Thermo Fisher Scientific) following the pipeline for Affymetrix Axiom genotyping workflow (Nicolazzi et al., 2014) and using the last version of buffalo genome assembly (UOA_WB_1) as the reference (Low et al., 2019).

Datasets Updating and Merging
We retrieved previously reported datasets from Colli et al. (2018a) and Deng et al. (2019a) available in the Dryad Repository ( Table 2). 2 To solve the incongruity between the marker positions when merging genotypes, we retained only those SNPs that were common in the three data sets. We further removed SNPs which mapped on sexual and mitochondrial bases, as well as those with unknown chromosome coordinates to perform an autosomal analysis. The remaining panel of 36,759 SNPs was then updated to the last version of the water buffalo genome (UOA_WB_1) through the commands -update-map, -update-chr, -updatealleles of PLINK V1.7 software (Purcell et al., 2007). The three datasets were then merged using the merge module in PLINK. Individually, the quality control filters were applied to remove variants with a minor allele frequency lower than 0.05 or a rate of missing genotyping >10%. The final dataset consisted of 36,014  SNPs and 477 individuals to perform the subsequent genetic population and structure analyses (Supplementary Table 1).

Genetic Relationship and Population Structure
A multidimensional scaling (MDS) analysis was performed to explore the genetic relationship between four different buffalo populations. To this end, a symmetric matrix of the identity-bystate (IBS) distances, for all pairs of individuals, was generated in PLINK V1.7 software (Purcell et al., 2007) based on the proportion of alleles shared and visualized using ggplot R package (Wickham, 2016). Observed (Ho) and Expected (He) heterozygosity for ten buffalo populations genotyped in this study (Bulgaria, Germany, Hungary, and Romania) were also estimated using PLINK V1.7 (Purcell et al., 2007). The maximumlikelihood approach was used to estimate the population structure through ADMIXTURE V1.3 Software (Alexander et al., 2009) with default parameters. This software modulates the probability of the observed genotypes considering the ancestry proportion and population allele frequencies. The best number of clusters (K-value) is estimated by the model as that reporting the lower cross-validation error. As an additional approach for the analysis of the population structure, a phylogenetic tree was built using the R package APE (Paradis and Schliep, 2019) based on computed Wright's Fixation Index (F ST ) (Wright, 1965) matrix obtained with Arlequin 3.5 (Excoffier and Lischer, 2010).

Run of Homozygosity
Analysis of run of homozygosity (ROH) was exclusively conducted on Central and East Europe buffalo populations genotyped in this study ( Table 1) and based on a panel of 61,813 autosomal SNPs available after the application of quality control filters, described above, to eliminate sexual, mitochondrial, and unknown chromosome coordinate SNPs, as well as those with a minor allele frequency lower than 0.01 and rate of 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 parameters: minimum size of 50 SNPs (-homozyg-snp 50), minimum density of SNPs of 1 SNP every 100 kb (-homozygdensity 100). We allowed one heterozygous SNP (-homozygwindow-het 1), a distance of homozygosity SNPs within the window of 250 kb (-homozyg-gap 250), and two missing SNPs (-homozyg-window-missing 2) per ROH. The identified ROH were cataloged in four categories according to its length: 1-5 Mb, 5-15 Mb, 15-30 Mb, and >30 Mb. Afterward, from ROH, we also calculated the inbreeding coefficient (F ROH ) using R package "DetectRUN" (Biscarini et al., 2018) and applying the formula: L ROH is the sum of the length of ROH for each individual, and L genome is the length of the genome analyzed, which was around 2.65 Gb in buffalo.   The frequency of occurrence of SNPs into an ROH was analyzed, identifying the genomic regions most commonly associated with ROH. Regions were selected containing the top 1% of most common ROH-associated SNPs (ROH hotspots or ROH islands) (Purfield et al., 2017;Mastrangelo et al., 2018).
Graphical visualization was performed using R package qqman (Turner, 2014) by constructing Manhattan plots. For the identification of genomic regions associated with ROH hotspots, we used NCBI map viewer of the water buffalo UOA_WB_1. Functional annotation was performed by the use of DAVID  V6.8 (Huang et al., 2009) and PANTHER (Thomas et al., 2003) software. The analysis pipeline applied in this study is provided as a flow chart (Figure 2).

Analysis of Genetic Diversity and Population Structure
The MDS plot ( Figure 3A) of 22 buffalo populations (Tables 1, 2) with different geographic distributions revealed a clear separation of European buffalo from other populations distributed worldwide. Indian-Pakistan and Latino-American buffalos were presented in a single cluster together with Bulgarian Murrah buffalos as the only exception of a European population. Four different groups could be distinguished in the rest of the European population: the Hungarian cluster, which included all three farms, the Romanian cluster, which showed one farm separated from the others (Rom_Serc), the Italian cluster, which included samples from Mozambique, and a German cluster from Brandenburg (Ger_Jüt), which was separated from the other populations. The remaining animals from Germany were mixed and scattered between Romanian and Italian populations. The admixture analysis carried out on genetic information from 22 buffalo populations available with global distribution ( Figure 3B) revealed a higher diversity among European populations than groups from India, Pakistan, the Middle East, and Latin America. 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 ancestry, except for Rom_Serc. Low admixture levels were found in the Hungarian group. A high admixture of buffalo was 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 exclusively among European populations (N = 14). The MDS plot ( Figure 4A) showed that after excluding Mozambique buffalos, the Italian samples were presented as a more homogenous and distinct cluster. The same appeared for the two Bulgarian populations. Furthermore, animals from two German farms (Lower Saxony -Ger_Stad, Mecklenburg-Western Pomerania -Ger_Born) had relations to the Italian cluster, while animals from Saxony (Ger_Wies) appeared to be related more closely to the Romanian buffalos. Buffalos from Brandenburg (Ger_Jüt)  (Table 3).
To support the interpretation of the Admixture and MDS analyses, we also measured the population differentiation due to genetic structure by an F ST pairwise distance analysis of 36,014 SNPs in all individuals of 22 populations distributed worldwide. The phylogenetic tree based on the F ST index matrix (Supplementary Table 2) was consistent with the Admixture analysis ( Figure 5). Two principal groups defined the difference between the Italian Mediterranean and the Murrah breeds. Along with the Italian group, Mozambique and three German herds, including Mecklenburg-Western Pomerania (Ger_Born), Lower Saxony (Ger_Stad), and Brandenburg (Ger_Jüt) clustered, while Brandenburg formed a separate branch. In the second group, Murrah breeds from Pakistan and India clustered with Latino America, Middle East, and Bulgarian samples. Hungarian and Romanian buffalos stood out, while animals from Saxony (Ger_Wies) grouped with animals from Romanian clusters.

Run of Homozygosity
In order to make better use of all information and to learn about the demographic past and/or recent history of the buffalo populations of Central and Southeast Europe, we performed an ROH analysis based on 285 samples genotyped in this study (Table 1), for which a panel of 61,813 autosomal SNPs was available (Figure 6). We found 10,797 ROH regions distributed in four categories (0-5, 5-15, 15-30, >30 Mb) with patterns different for each population. Among German buffalos, Ger_Born had the highest average of ROH in 0-5 Mb size category, indicative of ancient inbreeding, which decreased in all the other categories. Ger_Jüt showed the opposite trend, with the highest average in the long ROH category (>30 Mb), suggestive of the small effective population size and more recent inbreeding.

Localizing ROH Hotspots and Gene Set Enrichment Analyses
In order to detect the ROH hotspots, Manhattan plots were built to identify genomic regions frequently associated with ROH (Figure 8) and the SNP locus with the highest frequency (%) in those ROH ( Table 5) for each German buffalo population. The thresholds to define the regions containing the top 1% of most common ROH-associated SNPs (ROH hotspots or  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 (Atg16l2) in chromosome 16 with 26, 21, and 15 hotspots SNPs, respectively.

European Buffalos Stand Out From Others
This is the first study aimed to explore the genetic background of German buffalos in relation to other European populations (Colli et al., 2018a;Deng et al., 2019a). Our knowledge about the genetic relationships of European breeds, especially in Central Europe and Germany are still poor. A recent study from Colli et al. (2018b) reported a high genetic variability within Eastern European populations and identified recognizable traces of Indo-Pakistani background in the genome of Mediterranean buffalo (Colli et al., 2018b). Here, we compared 22 populations worldwide. The separation between Murrah (India-Pakistan, Latino America, and the Middle East) and Mediterranean (European populations) breeds was clearly evident. European buffalo populations are characterized by higher genetic diversity  (Bartosiewicz, 1999). Admixture analysis revealed that European and non-European populations do not share the same gene pool. Only Bulgarian Murrah buffalos were profoundly different from the other European animals. In K = 4, a large amount of alleles were common with Murrah breeds from India-Pakistan but not in K = 6 ( Figure 3B), indicating the development of a different genetic background. In fact, the values of observed and expected heterozygosity (He = 0.40) were comparable to those previously reported by Colli et al. (2018b) in Murrah breeds from India-Pakistan, underlining their similarity and indicating an isolatebreaking effect, the mixing of previously isolated populations.
From the data provided here, we have strong evidence that European buffalos stand out from other buffalo populations worldwide. We may thus assume ancient origins of European buffalos (excluding Bulgaria), supporting the idea of Krawczynski (2010) postulating that buffalos did not wholly disappear from the European continent during the late Pleistocene.

Different Genotypes in European Buffalos and a Novel One in Hungary
Four clear clusters could be distinguished in the MDS plot when including samples from five European countries only, representing the geographical distribution of Bulgarian, Hungarian, Romanian, and the Italian population. German populations, however, were mixed and located in between the Italian and Romanian groups. In the Romanian cluster, samples collected from the National breeding facility (Rom_Serc) were separated from the other farm located in a former Hungarian village (Mera), now part of Romania. The origin of Romanian buffalo is still unclear. According to Borghese, it is a Mediterranean type adapted to the cold climate and local environment and is classified as a Mediterranean Carpathian breed (Borghese, 2005). Because of the massive decrease of their population size during the last 20 years, they are nowadays in an endangered status, distributed in small farms in some villages, and almost 98% in Transylvania used for local food production and draft strength (Matiuti et al., 2020). Our results revealed that Mera buffalos (Rom_Mera), which were distributed in multiple smaller familiar farms, were characterized by a mixture of the Italian Mediterranean, Hungarian, and Bulgarian genetic background, while buffalos from the National breeding facility (Rom_Serc) have a defined genotype. In this national breed, a higher level of F ROH was found than in Mera buffalos. Moreover, the Mera herds (0.074) showed the highest proportion of long size ROH that confirms the small population size but is also indicative of recent inbreeding events. It may be due to the application of selection programs and the sign of environmental adaptation of the Mediterranean Carpathian buffalo in Rom_Serc population, while in Mera the animals are less subjected to 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 in small private farms used for family husbandry but mainly kept in national parks as a protected reserve (Borghese, 2005). In the MDS plot, the Hungarian cluster, including three populations (Földes, Csákvar, and Tiszataj), showed a lower admixture with Balkan populations (Figure 4B). The comparably high level of inbreeding coefficient and a high proportion of ROH, especially in short length, is indicative of ancient inbreeding events in Hungarian buffalos. Notably, the genetic background in Hungarian buffalo is comparably homogenous compared to all genotypes published before or presented here and may represent a valuable resource for local economy or breeding programs worldwide in the future. Clearly, and maybe as a surprise, the genetic background in Hungarian buffalo is different from the population in Mera, which represents a former Hungarian village. In all three Hungarian farms included in this study, the buffalos are used for meat production only and serve particularly for the establishment of local marketing together with other regional products. Under current natural breeding, there is no selection for milk production at present in the herds included, which may explain the lack of recent inbreeding. Bulgarian Murrah buffalos, instead, had the lowest proportion of ROH in each length category and the lowest F ROH values between Central and East European populations. Indeed, the Bulgarian is the biggest population analyzed in this study, and the samples were collected from two different locations, which also could be related to comparably low levels of inbreeding. In the German populations, genomic relationships with Bulgarian Murrah, Italian Mediterranean breeds, and Romanian genetics were identified. Buffalos in Germany had the highest observed and expected heterozygosity values within all Central Eastern European populations studied here. The observed genetic diversity is a significant factor for their potential adaptation to the local environment and provides a valuable substrate for an upcoming selection program. Drift rather than mutation has likely created this variability due to the short period in generations and the relatively small effective population sizes (Wang, 2005). Buffalos from Brandenburg (Ger_Jüt) stood out from all other German populations studied here. As the only German population, clear genetic relations to Bulgarian buffalos can be postulated for the Jüterbog herd. In fact, part of the Bulgarian Murrah's gene pool that was visible in the admixture at K = 4, and disappeared at K = 6, is likely an effect of genetic drift after isolation over time. The Brandenburg buffalo herd showed the highest deviation of heterozygosity values, with He being lower than Ho, suggesting Hardy-Weinberg disequilibrium and inbreeding (Mayo, 2008). The higher values of inbreeding with an increasing trend of ROH average from short to long size may further suggest recent inbreeding and isolation without any admixture since the introduction of the buffalos in Germany.

Overlapping ROH Signals in German Buffalo Genome
The control of inbreeding in a population is necessary to avoid depression of phenotype traits (Ouborg et al., 2010;Fontanesi et al., 2015). Considering that ROH regions are not randomly present in the genome (McQuillan et al., 2008), the detection of the genomic regions frequently covered with ROH is useful. This helps to understand whether they influence any quantitative trait, especially in a population of small population size, such as local breeds (Biscarini et al., 2015). The negative effect of inbreeding depression in milk production has been reported in Egyptian (Khattab and Kawthar, 2007), South Iranian (Mirhabibi et al., 2007), and Brazilian buffalo populations (Santana et al., 2011). We found several differences among German buffalo populations in ROH analysis in terms of the number of hotspot SNPs, candidate genes in the same genomic coordinates, and biological processes in which they could be involved. In the Brandenburg (Ger_Jüt) population, the most distinguished, with a high level of inbreeding, we found the highest number of hotspot SNPs (949), distributed over 11 chromosomes and concerned only 57 genes involved in 11 biological processes. While Saxony buffalos (Ger_Wies), the population with the lowest inbreeding level and the lowest number of hotspot SNPs (635), showed their distribution limited to 6 chromosomes but to 83 genes involved in 15 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 the gene Atg16l2 (autophagy related 16 like 2) influences the adaptation of the immune system to the recovery from mastitis in Danish Holstein cattle (Welderufael et al., 2018). From these very initial interpretations of the potential effects of local selection for functional traits, we may get an idea of how genomic selection could be used to improve animal health and performance in the future. Certainly, additional studies of gene expression, metabolism, or animal health are needed to link SNP markers with particular phenotypes. We are aware of the limitations of our study. Only four different farms in Germany are included that represent less than 15% of the German buffalo population. Moreover, the selection of farms in Eastern Europe was based on potential familiar relationships and can be seen as a start only. For German and other European countries, higher coverage, including more herds and additional countries to further exploit the full setting of genetic diversity in Europe, and an in-depth analysis of signatures of past selection combining several metrics together with ROH is required. Finally, in the genotyping array designed for buffalo, only 30% of SNPs from European breeds (Italian Mediterranean) have been included, which could 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).

CONCLUSION
This work aimed to contribute to an improved understanding of Central and Southeast European buffalo populations from a genomic perspective. We defined the genetic characteristics of European buffalo populations in Romania, Bulgaria, and Hungary. Together with the Italian Mediterranean breed, these populations contribute to the high degree of genetic diversity of 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 any other genotype described so far. Since all farms are in the close neighborhood, follow-on studies may help to unravel the origin of the Hungarian population. This novel genetic identity found in Hungarian buffalos thus adds one additional specific breed to the existing genetic diversity in Europe, which so far consisted only of the Italian Mediterranean, Bulgarian Murrah, and Romanian genetics. This is the first genetic characterization of buffalos in Germany, where farming started only about 40 years ago and therefore is in its infancy. We have identified familiar relations to all genotypes identified in European buffalos so far with the exception of Hungary. Importantly and based only on a comparably small number of herds, a diverse genetic make-up is available in European and German buffalos, which can be seen as an excellent prerequisite for the development of Buffalo-based bioeconomy in Europe based on local breeds characterized by unique features.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because sampling occurred as part of livestock breeding programs which are granted by international law and therefore are not reviewed by ethical committees. Written informed consent was obtained from the owners for the participation of their animals in this study.

AUTHOR CONTRIBUTIONS
AN, KW, and AH conceived the project and designed the experiments. AN performed the laboratory part. AN, SQ, RG-P, and ML-S performed the data analysis. AN, SQ, RG-B, JB, RK, and AH contributed to the interpretation of results and wrote the manuscript. MiT, M-AF, HP, AB, LV, CH, LH, PP, YI, TP, WL, MaT, and AH collected or provided samples and contributed to data production. All authors read, made corrections, and approved the final version of the manuscript.

ACKNOWLEDGMENTS
This study would not have been possible without the help of many buffalo breeders. We express our gratitude to Family Henrion (Bobalis Agrargesellschaft mbH in Jüterbog), Dr. Möhring (Gut Darss GmbH & Co. KG in Born), the staff in Stadland (Hof am Meer, Stadland/Norderschwei), and all helpers from Wiesenburger Land eG in Wiesenburg. We further want to thank the managers and helpers in Hungary: Mr. Mihály Bodnár, Tiszatáj Fundation (H-4450 Tiszalök Rákóczi u. 14), Mr. Levente Viszló Pro Vértes Fundation (H-8083 Csákvár, Kenderesi u. Geszner ház), and Mrs. Gabriella Bona Bihar Fundation (H-4177, Földes Kálóháti Tanya), the staff in Sercaia (Romania) and all the families in Mera (Romania) that supported our project. Furthermore, the local support in Bulgaria is thanked here. Finally, we would like to thank Dr. Taina Figueiredo Cardoso for helpful advice when performing genetic diversity analyses.