Evaluation of Genetic Diversity and Structure of Turkish Water Buffalo Population by Using 20 Microsatellite Markers

Simple Summary In the present study, twenty microsatellite loci were tested to assess and analyze the genetic diversity between and within 17 different populations of Turkish water buffalo. The total number of animals sampled was 837, collected from six geographical regions: Marmara Region (MRM), Black Sea Region (BSR), Aegean Region (AER), Central Anatolia Region (CAR), Eastern Anatolia Region (EAR) and Southeastern Anatolia Region (SAR). All studied microsatellites markers showed allelic polymorphism. In this study, the results indicated a definite genetic diversity among the Turkish water buffalo populations which indicates the existence of at least two major clusters. Abstract The present study was aimed to investigate the genetic diversity among 17 Turkish water buffalo populations. A total of 837 individuals from 17 provincial populations were genotyped, using 20 microsatellites markers. The microsatellite markers analyzed were highly polymorphic with a mean number of alleles of (7.28) ranging from 6 (ILSTS005) to 17 (ETH003). The mean observed and expected heterozygosity values across all polymorphic loci in all studied buffalo populations were 0.61 and 0.70, respectively. Observed heterozygosity varied from 0.55 (Bursa (BUR)) to 0.70 (Muş (MUS)). It was lower than expected heterozygosity in most of the populations indicating a deviation from Hardy–Weinberg equilibrium. The overall value for the polymorphic information content of noted microsatellite loci was 0.655, indicating their suitability for genetic diversity analysis in buffalo. The mean FIS value was 0.091 and all loci were observed significantly deviated from Hardy–Weinberg Equilibrium (HWE), most likely based on non-random breeding. The 17 buffalo populations were genetically less diverse as indicated by a small mean FST value (0.032 ± 0.018). The analysis of molecular variance (AMOVA) analysis indicated that about 2% of the total genetic diversity was clarified by population distinctions and 88 percent corresponded to differences among individuals. The information produced by this study can be used to establish a base of national conservation and breeding strategy of water buffalo population in Turkey.


Introduction
Water buffaloes have been reported to be of great importance to the lives of farmers and thus to the economies of many countries worldwide [1]. The number of water buffaloes in the world has increased rapidly over the past few decades, and according to Food and Agriculture Organization (FAO) statistics, there are about 208 million buffaloes in the world. Most of the world's buffaloes live in Asia (96.79%), Africa (1.68%), the Americas (1.23%) and Europe (0.22%) [1,2].
According to 1974 FAO statistics, at that time, there was one million buffalo heads in Turkey. From 1984 to 1998, there has been a decrease in the buffalo breeding population of water buffalo populations [12,17,18] The genetic identification of Turkish water buffalos takes importance for the conservation of genetic diversity in indigenous breeds. The aims of this study were (I) to assess the genetic diversity within and between the Turkish water buffalo population; (II) to estimate the level of inbreeding, using 20 microsatellite loci; and (III) to identify the genetic relationship and describe geographical and genetic distinction between different water buffalo populations at different sites in Turkey.

Materials and Methods
This experiment was conducted at the National Anatolian Water Buffalo Breeding Project breeders buffalo in 6 different geographical regions in Turkey. The experimental protocol was approved by the Republic of Turkey Ministry of Agriculture and Forestry Pendik Veterinary Control Institute Local Ethics Committee (AEC approval number: 12/2013). This study was conducted between June 2013 and June 2015.

Sampling and DNA Extraction
In this study, a total of 837 blood samples from unrelated Turkish water buffalo (Bubalus bubalis) individuals were collected at random from 17 different populations located in 6 different geographical regions (17 cities, 37 districts and 119 villages) (Supplementary Figure S1). These regions have been chosen to represent the expected Turkish water buffalo populations. Most of these samples are local populations belonging to exceedingly small farms.
Blood samples were collected in vacutainer tubes that include EDTA (0.5 mM, pH 8.0) and stored at +4 • C until DNA extraction. Around 8 mL of blood per buffalo was taken, and the purification of genomic DNA was performed from 2 mL blood samples, using High Pure PCR Template Preparation Kit (Roche Molecular Systems Inc., Pleasanton, CA, USA), following the manufacturer's protocol. All DNA extraction and polymerase chain reaction (PCR) amplification were performed by Geometry Biotechnology (http://www.genometri.com.tr/, accessed date 10 August 2020) in Istanbul.

Microsatellites Amplifications and Analysis
A total of 20 heterologous buffalo microsatellite loci were studied from a panel of 30 markers recommended by the International Society for Animal Genetics (ISAG) and Food and Agriculture Organization of the United Nations (FAO) working group for biodiversity study. The criterion for selection of the microsatellite loci was based on the high polymorphism information content value (PIC) and the number of exhibited alleles of the loci [12]. Information on the twenty-microsatellite investigated is presented in Supplementary Table S1. These microsatellites were analyzed to estimate various genetic diversity parameters. The microsatellite markers analysis was performed by using an Applied Biosystems 3130 Genetic Analyzer. The microsatellite markers were grouped in five sets of fluorescent-labeled primers. Five primer pairs were performed in each set for multiplex amplification. The forward primer for each locus was labeled with one of the four fluorescent dyes FAM, HEX and TAMRA (Applied Biosystems, Foster City, CA, USA) (Supplementary Table S1). The polymerase chain reaction analyses were carried out in a T100TM Thermal Cycler (Bio-Rad), using the primers listed in Supplementary Table S1. The reaction mixture was composed of genomic DNA (50 ng), 200 mM dNTPs, 2.0 mM MgCl 2 , 1X PCR buffer, 10 pmol forward and reverse primers and Taq DNA polymerase (0.5 u/sample). The PCR was performed with a total reaction volume of 25 µL, using the following thermal conditions, 94 • C for 10 min, followed by 32 cycles of 94 • C for 1 min, 55 • C for 30 s, 56 • C for 3 min and a final extension at 60 • C for 1 h. Amplified DNA was controlled in 1% agarose gel, using SYBR™ Safe DNA Gel Stain.
After agarose gel electrophoresis for 20 microsatellites, allele size was identified on all samples with an ABI Prism ® 3130 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA), using the GeneScan ® Analysis Software (Applied Biosystems, Foster City, CA USA), which check different alleles via size comparison with standard DNA size

Data Analysis
Allele frequencies, the total number of alleles, the mean number of alleles (N a ), the number of effective alleles (N e ), allelic richness (R S ), polymorphic information content for each locus (PIC), observed heterozygosity (H O ), expected heterozygosity (H E ), Shannon's information index (I), Hardy-Weinberg equilibrium and null allele frequencies, using Genetix v4.05 [21], FSTAT v2.9.3.2 [22], POPGENE Version 1.31 [23] and GenAlEx Version 6.5 [24]. Wright's F statistics (F ST , F IS and FIT), as proposed by Weir and Cockerham [25], were analyzed by using Genetix ® software [21]. Nei's gene diversity (H T ), diversity within populations (H S ), diversity between populations (D ST ) and coefficient of gene differentiation (G ST ) values were analyzed with FSTAT v2.9.4 [22]. Exact tests for deviation from the Hardy-Weinberg Equilibrium (HWE) and partitioning of genetic diversity using analysis of molecular variance (AMOVA) were analyzed, using the AR-LEQUIN v. 3.5.2.2 [26]. Pairwise genetic distances (Reynold's genetic distance) and Nei's unbiased D AS genetic distances (Nei,[27]) were calculated by using the Populations v 1.2.30 software. Neighbor-net dendrogram constructed from Reynold's genetic distances by using SplitsTree v4. 16.0 [28]. Principal components analysis (PCA) was calculated for the 20 microsatellites, using NTSYSpc V2.10q software [29]. Genetic diversity and the degree of admixture of Turkish water buffalo populations were analyzed by using the Bayesian clustering procedure of STRUCTURE ver. 2.3 [30]. Twenty replicate runs were calculated for each K between 1 and 20, with a burn-in period of 1,000,000 iterations, followed by 500,000 iterations of the Markov chain Monte Carlo algorithm. To determine the most possible groups (K) that best fit the data; we used the STRUCTURE HARVESTER [31], which implements the Evanno method [32]. Evanno's method was carried out to determine the suitable number of clusters using ∆K, due to the rate of change in the log probability of the data. The program CLUMPP ver. 1.1 [33] was used to align the 20 repetitions of each K. CLUMPP software, an online web-based program, was performed for collating the outcomes produced by the program STRUCTURE. The clustering pattern was applied in the CLUMPP program and visualized by way of the software DISTRUCT software version 1.1 [34].

Results
In the present study, 20 microsatellite markers were used to analyze the relationships within and among 17 Turkish water buffalo populations. A total of 837 individuals were sampled from 17 provinces, 37 districts and 119 villages belonging to six geographical regions that represent the most important sites of water buffalo breeding.

Genetic Diversity
A summary of statistic results for genetic diversity is shown in Table 1 and Supplementary Table S2. The mean number of alleles (N a ), the number of effective alleles (N e ), observed (H O ) and unbiased expected (H E ) heterozygosity, Shannon's information index (I) and the deficit of heterozygotes (F IS ) values for individual subpopulations and overall population are presented in Table 1. The properties of the analyzed microsatellite, along with the genetic variation statistics, were listed in Supplementary Table S2. For the entire 17 Turkish water buffalo populations, a total of 190 alleles were found in 837 animal genotypes for the 20 microsatellite markers.
The mean number of alleles (N a ) was the highest in Central Anatolia (Kayseri (KAY) and Sivas (SVS)) Region buffalo individuals (8.80) and the lowest in Black Sea Region (Giresun (GIR)) buffalo individuals (5.70). The number of effective alleles per locus (N e ) and allelic richness (R s ) values are a measure of the genetic diversity. The effective number of alleles per population (N e ) showed lower values, which varied from 2.87 (Bursa (BUR)) to 4.14 (Muş (MUS)). Additionally, the R S means ranged from 5.03 to 6.83 in the Marmara, Black Sea and Aegean Regions buffaloes, being lower than that in the Central Anatolia, Eastern Anatolia and Southeastern Anatolia Regions buffaloes. The observed (H O ) and unbiased expected (H E ) heterozygosity per population varied from 0.55 (BUR) to 0.70 (MUS) and 0.58 (BUR) to 0.73 (MUS), respectively ( Table 1). The H O and H E values in the Marmara, Black Sea and Aegean Regions buffaloes were lower than the Central Anatolia, Eastern Anatolia and Southeastern Anatolia Regions' buffaloes. Genetic diversity of the Turkish water buffalo population of Central Anatolia, Eastern Anatolia and Southeastern Anatolia Regions was observed to be higher than in other regions (Marmara, Aegean and the Black Sea), but the differences were not statistically significant. It is thought that the reason for this difference may be due to the geographical proximity of the buffalo to the domestication center, the fact that the buffalo breeding in these regions consists of small family businesses and there is no selection in this population. The inbreeding coefficient (F IS ) within the populations varied between 0.017 and 0.183. Additionally, the SVS and Amasya (AMS) populations have shown that the lowest and highest genetic diversity, respectively. All of the F IS values were determined to be positive and statistically significant (p < 0.05, p < 0.001) ( Table 1). The average F IS value, which describes the excess or deficit of heterozygotes within subpopulations was 0.117 (p < 0.05) and therefore different from zero. Shannon's information index (I) across populations ranged from 1.15 (BUR) to 1.58 (MUS).   When the results are evaluated based on microsatellite loci, a total number of alleles per locus varied from 6 (ILSTS005) to 17 (ETH003), while the mean number of alleles per locus (N a ) ranged between 4.29 and 12.59 for the same loci (Supplementary Table S2). N e ranged from 1.80 (CSSM033) to 7.74 (CSSM047) for the Turkish water buffalo population based on 20 polymorphic microsatellite loci. The value of allelic richness (R s ) ranged from 3.99 (ILSTS005) to 12.08 (ETH003), with a mean of 6.87. PIC and the Shannon information index are another measure of genetic variability indicating the informativeness of the assessed loci. The polymorphic information content (PIC) was analyzed for each locus and varied from 0.412 (CSSM033) to 0.859 (CSSM047), which has the highest number of alleles per locus in the current study. A total of 16 of these 20 microsatellite loci had polymorphic information content (PIC) values greater than (0.5), which make them useful in genetic diversity studies. Additionally, four loci (CSSM033, ILSTS005, CSSM032 and CSSM029) showed moderate polymorphism (PIC > 0.40) (Supplementary Table S2). Shannon's information index (I) across populations ranged from 0.91 (CSSM033) to 2.32 (CSSM045). The value of gene flow (Nm*) between the 17 subpopulations was positive and varied between 2.05 (CSSME070) and 26.87 (CSSM029) for different microsatellite loci. This confirms that Turkish buffalo samples were exchanged between the 17 subpopulations.
The values of fixation indexes (F IS , F ST and F IT ) for the overall populations are given in Supplementary Table S2. Most of the markers had positive values for F IS (Supplementary  Table S2), showing a deficiency in heterozygosity. The F IS index was negative for ILSTS005, CSSM036 and CSSM029 markers indicating a high frequency of heterozygotes in these loci.  Supplementary Table S2.
The average coefficient of gene distinction (G ST ) over the 20 loci was 0.030 ± 0.021 (p < 0.01). The G ST values for single loci varied from −0.001 for CSSM029 to 0.098 for CSSME070. The gene differential coefficient G ST (3%) indicated that most of the total genetic variation was due to intra-population difference and only a few existed among populations, which implied these Turkish water buffalo populations had relatively less genetic diversity and distinctiveness. Nei's gene diversity index (H T ) for loci ranged from 0.44 (CSSM033) to 0.88 (CSSM045), with an average of 0.69.
Eighteen of the 20 microsatellite loci found a highly important departure from Hardy-Weinberg equilibrium (HWE) (p < 0.001) in the whole population, whereas the other two loci showed different significant differences (p < 0.01, p < 0.05) (Supplementary Table S2). But, when considering populations separately, many markers per population were in Hardy-Weinberg disequilibrium (p > 0.05). The number of these markers ranged between 5 loci inİstanbul-Çatalca (IST) to 14 loci in Tekirdag (TEK) and GIR populations. All the studied populations performed highly significant departure (p < 0.001, p < 0.01, p < 0.05) from the Hardy-Weinberg equilibrium when considering all microsatellite loci. This departure and the high positive mean values of F IS (Supplementary Table S2) may remark the presence of heterozygote deficiencies, which could be the consequences of an uncontrolled mating between populations. The presence of null alleles, defined as non-amplifying alleles, because of mutations at PCR priming sites, causes overestimation of both F ST and genetic distance values. The null allele frequencies varied from 0.017 (CSSM036) to 0.241 (BMC1013). We identified only one locus (BMC1013) that had potential null alleles at high frequency (r ≥ 0.2) in at least one breed. It has been observed that the frequencies of other microsatellite loci are lower than (0.20) (Supplementary Table S2 Table S5) of Turkish water buffalo, showed an overall genetic differentiation F ST of (0.032 ± 0.018) and pairwise F ST values ranging from 0.0000 (SVS vs. KAY; KAY vs. Bitlis (BIT)) to 0.0866 (BUR vs. Çorum (COR)) (varying from Indian red, white and blue colors in Figure 1). Significant genetic variation was observed after sequential Bonferroni correction in 118 out of 136 population pairs (Supplementary  Table S5). (SIN)), (between I and II) from only the Aegean Region's buffaloes (AFY) and (II) from the Marmara, Black Sea, Central Anatolia, Eastern Anatolia and Southeastern Anatolia Regions' buffaloes in Turkey (Figure 2). The phylogeny of Reynold's distances (Supplementary Table S4) was similar to that performed using Nei's DAS distances (Supplementary  Table S3).   Table S4) was similar to that performed using Nei's D AS distances (Supplementary Table S3).
The analysis of molecular variance (AMOVA) test was performed to evaluate genetic variability is distributed within and among populations. We calculated possible structures by composing and contrasting different population groups. We performed the analysis under two hypotheses: For Hypothesis (I), the AMOVA analyses results showed that most of the molecular variation occurred within individuals (88.33%), while it represented 0.55% among geographic groups, 2.69% among populations within groups and 8.42% among individuals within populations (Table 2). Variance components among groups, among populations within groups and within individuals were significant (p < 0.001) for all the studied loci (Table 2), implementing significant geographical distribution in Animals 2021, 11, 1067 8 of 17 studied buffalo populations. Furthermore, the variance component among individuals within populations was significant (p < 0.05). For Hypothesis (II), AMOVA analyses outcomes indicated that the variation among groups, among populations within groups, among individuals within populations and within individuals was 2.15%, 2.02%, 8.34% and 87.49%, respectively. Variance components among groups, among populations within groups, among individuals within populations and within individuals were significant (p < 0.001) for all the studied loci implementing significant Reynold's genetic distances distribution in studied Turkish water buffalo populations (Table 2).  The analysis of molecular variance (AMOVA) test was performed to evaluate genetic variability is distributed within and among populations. We calculated possible structures by composing and contrasting different population groups. We performed the analysis under two hypotheses: For Hypothesis (I), the AMOVA analyses results showed that most of the molecular variation occurred within individuals (88.33%), while it represented 0.55% among geographic groups, 2.69% among populations within groups and 8.42% among individuals within populations (Table 2). Variance components among groups, among populations within groups and within individuals were significant (p < 0.001) for all the studied loci (Table 2), implementing significant geographical distribution in studied buffalo populations. Furthermore, the variance component among individuals within populations was significant (p < 0.05). For Hypothesis (II), AMOVA analyses outcomes indicated that the variation among groups, among populations within groups, among individuals within populations and within individuals was 2.15%, 2.02%, 8.34% and 87.49%, respectively. Variance components among groups, among populations within groups, among individuals within populations and within individuals were significant (p < 0.001) for all the studied loci implementing significant Reynold's genetic distances distribution in studied Turkish water buffalo populations ( Table 2).
The principal component analysis (PCA) of 17 Turkish water buffalo populations due to 20 microsatellite markers is shown in Figure 3. The 44.46% of the total genetic variation present between the populations was explained by the first three axes of the PCA test. This is an acceptable fit, given the small amount of variability from a large number of samples and microsatellite alleles used in the analysis. The PCA analysis classified Turkish water buffalo populations into two basic cluster involving different region populations to some extent, i.e., Cluster I is BSR buffaloes (GIR, SIN, AMS, TOK and COR); Cluster II is the first group, namely BSR, CAR, EAR and SAR buffaloes (  The principal component analysis (PCA) of 17 Turkish water buffalo populations due to 20 microsatellite markers is shown in Figure 3. The 44.46% of the total genetic variation present between the populations was explained by the first three axes of the PCA test. This is an acceptable fit, given the small amount of variability from a large number of samples and microsatellite alleles used in the analysis. The PCA analysis classified Turkish water buffalo populations into two basic cluster involving different region populations to some extent, i.e., Cluster I is BSR buffaloes (GIR, SIN, AMS, TOK and COR); Cluster II is the first group, namely BSR, CAR, EAR and SAR buffaloes (Samsun (SAM), DYB, MUS, BIT, SIV and KAY), and the second group, namely MRM and BSR buffaloes (BAL, TEK, IST and Düzce (DUZ)) ( Figure 3). The Aegean Region's buffaloes (AFY) are between the II and III groups. The first component (PC1), which was responsible for 22.42% of the genetic variation, separated the BUR population from all the other studied populations. The second component, which represented 35.87% of the genetic variations, separated the DUZ population (Figure 3).  The genetic structure of each population was identified regarding admixture level for each water buffalo individual using a correlated allele frequencies model implemented within the STRUCTURE software. The results of Delta K (∆K = 92.42) represented that the optimal number of genetic clusters indicating most like ancestral breeds was at K = 2 ( Figure 4A). The value indicates that the studied water buffalo populations were better identified by two genetic clusters instead of 17 populations ( Figure 4B). The two clusters' backgrounds were made up of IST, DUZ, TEK, Balıkesir (BAL), BUR, SAM, DYB, MUS, BIT, KAY and SVS in the first (red color), and AFY, GIR, AMS, TOK, COR and SIN in the second (dark blue color) cluster ( Figure 4B). Approximately 51% of the individuals were classified within their source cluster assuming a threshold of q ≥ 0.700, whereas, for more stringent threshold q values, only~38% (q ≥ 0.900),~24% (q ≥ 0.950) and 0.23% (q ≥ 0.999) of the individuals were correctly assigned. There were several heterogeneous populations (AFY, SAM and BAL) with less than 48-30% of the individuals assigned to their source cluster (for q ≥ 0.700). Consequently, the SAM population evidenced two clusters more than an intermediate position between the two reference populations.

Discussion
The studies of genetic diversity play a significant role in developing breeding strategies and programs for livestock. In order to retain genetic variation, breeding strategies that increase effective population size minimizing the genetic drift effect should be performed. Microsatellite markers in combination with recent scientific applications of statistic showed a powerful tool for the conservation of native breeds [35].
The advantage of using microsatellite DNA polymorphism to estimate genetic diversity among breed and among associated populations has been researched in livestock such as water buffalo [8,[36][37][38][39][40][41][42], cattle [43][44][45][46], sheep [47][48][49] and goat [35,50,51]. To determine genetic diversity among water buffalo populations, several microsatellite studies have been published [12,17,18] in Turkey, but a study of genetic diversity of Turkish water buffalo that has collected 837 individuals from six geographical regions has not been performed until now. In the present study, we wanted to provide basic data on the genetic variation and population structure of Turkish water buffalo populations, using 20 microsatellite markers to provide a foundation for a more comprehensive genetic resource protection and genetic management.

Discussion
The studies of genetic diversity play a significant role in developing breeding strategies and programs for livestock. In order to retain genetic variation, breeding strategies that increase effective population size minimizing the genetic drift effect should be performed. Microsatellite markers in combination with recent scientific applications of statistic showed a powerful tool for the conservation of native breeds [35].
The advantage of using microsatellite DNA polymorphism to estimate genetic diversity among breed and among associated populations has been researched in livestock such as water buffalo [8,[36][37][38][39][40][41][42], cattle [43][44][45][46], sheep [47][48][49] and goat [35,50,51]. To determine genetic diversity among water buffalo populations, several microsatellite studies have been published [12,17,18] in Turkey, but a study of genetic diversity of Turkish water buffalo that has collected 837 individuals from six geographical regions has not been performed until now. In the present study, we wanted to provide basic data on the genetic variation and population structure of Turkish water buffalo populations, using 20 microsatellite markers to provide a foundation for a more comprehensive genetic resource protection and genetic management.

Genetic Diversity of Turkish Water Buffalo
All measures of genetic diversity: the mean number of alleles (N a ), the number of effective alleles (N e ), allelic richness (R s ), Shannon's Information Index (I) and polymorphic information content values (PIC) showed that most of the studied loci were highly informative, indicating high polymorphism across the loci, thus suggesting the suitability of these microsatellite markers for genetic diversity studies in Turkish water buffalo populations. The values of genetic diversity parameters were higher compared with a similar study of Asian water buffalo population (Pakistan, Thai, Indonesia, Egypt and China) [42,[52][53][54][55], South American buffalo populations (Brazil and Colombia) [10,38,56], North American buffalo populations (Cuba) [36,56,57], southeast of Central Europe (Romania) [58] and Southern Europe water buffalo populations (Italia and Greece) [41]. Molecular genetic parameters (N a , N e , PIC and HO) obtained from this study were higher than previous research in Turkish water buffalo [12,17,18]. These different studies, both the loci and the number of loci involved were different from each other. Thus, differences in the results may partly be attributable to the differences in the loci employed. These results indicated that microsatellites used in the present study have a high confidence to reveal genetic diversity for Turkish water buffalo populations.
The allele diversity for Turkish water buffalo populations was lower than that of the Colombian buffalo population studied with 10 microsatellite loci [10] but much higher than that reported for Cuban, Pakistan, Romanian, Egyptian, Iranian and Brazilian water buffalo populations [52,56,[58][59][60]. Turkish water buffalo populations showed a drastic low number of the effective number of alleles (even lower than half) than the observed mean number of alleles. This is due to the very low frequency of most of the alleles at each locus and a very few alleles might have contributed a major part of the allelic frequency at each locus. Allelic richness was considerably high in the CAR's, EAR's and SAR's buffaloes, indicating high genetic polymorphism as expected heterozygosity (>0.6912).
Another measure of genetic variability is expected heterozygosity where the maximum expected heterozygosity (0.7280) was showed in the MUS population and the minimum (0.5833) was shown in the BUR population. This study results indicate the 17 tested Turkish water buffalo population have substantial amount of genetic diversity, when compared to some other water buffalo breeds around the world where, Pakistan water buffalo breeds such as, Nili (H E = 0.  [14]), some of Turkish water buffalo populations (H E = 0.5359; Özkan Ünal et al. [12]), Asian swamp buffalo (H E = 0.50; Barker et al. [42]) and Chinese buffalo (H E = 0.53, Zhang et al. [55]). On the other hand, the buffalo populations tested in our study showed less genetic diversity when compared to Iranian indigenous buffalo populations (H E = 0.75; Darestani et al. [59]) and Iraqi buffalo populations (H E = 0.86; Jaayid and Dragh, [61]), Indian river buffalo breeds (H E = 0.71-0.78; Kumar et al. [39]), and African buffalo (H E = 0.759; Van Hooft et al. [62]).
The estimate of inbreeding value (F IS ) shows the excess of homozygosity in a subpopulation and, with reference to molecular markers, informs if a pattern of reduction in diversity based on several causes exists. The positive F IS values showed heterozygotes deficiency within populations. This deficit might be because of inbreeding and the Wahlund effect. The mean F IS value of the Turkish water buffalo populations was 0.1170 (p < 0.05) and is similar to that obtained by Uffo et al. [36] in Cuban water buffalos, but lower than what was described by Ángel-Marín et al. [10], using 10 microsatellite markers in Colombian buffalo herds; other authors describe a lower value of F IS [56][57][58][59]. The F IS in Turkish water buffalo can be considered higher compared with other populations (F IS = 0.1170) like the results reported by Shokrollahi et al. [9], who found a value of 0.047 in the Iranian river buffalo. The highly significant (p < 0.001) F IS value (0.1820) showed a rather high inbreeding degree within the population. The heterozygote deficiency performed in the AMS population could be because of the higher rate of inbreeding, the population subdivision (Wahlund effect) and the presence of "null alleles" (non-amplifying alleles). Another possible reason for the high values of inbreeding could be due to their small population size, a small number of breeding males and their limited geographical area of dispersion [63].
The F IT values, which measure the heterozygosity loss of the individual concerning the overall population, were 0.119 (p < 0.05), indicating that there is a general lack of heterozygous individuals in the Turkish water buffalo populations of 12%. In the study that the F IT value was considerably lower than the Pakistan buffalo populations [52], Cuban water buffalo breed [36] and Turkish water buffalo populations [17,18], and higher than the values reported by Popa et al. [58], Darestani et al. [59], Marrero et al. [56] and Acosta et al. [57].
The presence of null alleles, defined as non-amplifying alleles, due to mutations at PCR priming sites, causes overestimation of both F ST and genetic distance values. The null allele frequencies in the studied microsatellite loci were below 20% except for BMC1013 loci (24.09%). The lowest and highest null allele frequencies were 0.168 (CSSM036) and 0.2409 (BMC1013), respectively. Taking this value into consideration, it has been implemented that the studied 19 loci can be safely used in paternity tests (except for BMC1013 loci).
Mean value of DST indicating genetic diversity between populations, GST, which is an important indicator of the relative magnitude of genetic differentiation, and HT, described as total genetic diversity, values were found as 0.021, 0.030 and 0.694, respectively. The D ST value obtained from this study can be considered as an indicator of low genetic diversity between 17 Turkish water buffalo populations. The average G ST value obtained from the overall loci pointed out that 3% of total genetic diversity resulted from the differences between the populations. In all other respects, it can be said that 97% of genetic variation is caused by the difference between individuals. The G ST values were considerably lower than Ángel-Marín et al. [10] in Colombian buffalo herds. All studied loci showed a significant deviation from the Hardy-Weinberg Equation (p > 0.05, p > 0.01, p > 0.001).

Genetic Structure of Turkish Water Buffalo
Population differentiation was analyzed by estimation of the F ST index. In regard to all pairwise differences (Slatkins linearized F ST ) in this study, the distribution of F ST represented low genetic divergence (0.000 < F ST < 0.0866) among populations in general. The F ST comparison values obtained were significant in 118 pairwise calculations (p < 0.05; p < 0.01; p < 0.001). The highest level of differentiation was obtained between COR-BUR and TOK-BUR populations (F ST > 0.0866), and the lowest between KAY-BIT and KAY-SVS populations (F ST = 0.000), respectively. Thus, the average proportion of genetic differentiation among breeds was 3.2%. This value is lower than the 16.8% found in another genetic study on Asian buffalo [42], Indian buffalo (3.4%) [39], Turkish water buffalo populations (6.2%) [17], Cuban and Brazilian buffaloes (7.5%) [56] and Pakistan buffalo (7%) [52], but higher than Chinese buffalo (2.8%) [16] and Iranian buffalo (1%) [59]. In the study, it was shown that the F ST value differences between the buffalo populations belonging to the provinces that are geographically close to each other are quite low but important. The reason for this is thought to be due to the male bull changes between nearby provinces.
Neighbor An AMOVA was implemented to analyze the relative contribution of different factors to the observed genetic variation, with each factor considered in a separate analysis, i.e., six groups according to the geographical prevalence (MRM, BSR, AER, CAR, EAR and SAR), three groups according to Reynold's genetic distances distribution: The first group is the BSR populations (GIR, AMS, TOK, COR and SIN); the second group is the AER population (AFY); third group comprises the MRM populations (IST, TEK, BAL and BUR), BSR buffaloes (SAM, DUZ), CAR buffaloes (KAY, SVS), EAR buffaloes (MUS, BIT) and SAR buffaloes (DYB). The AMOVA analysis data implemented that the majority of the obtained variance is because of differences among individuals within populations. Most of the variation is obtained within the individuals (88.33% Hypothesis I and 87.49% Hypothesis II), yet the differences among groups show only 0.55% and 2.15% of the variation, respectively. Among groups, among populations within groups and among individuals within populations were also an important origin of variation (p < 0.001, p < 0.05), though fundamentally smaller than the within individual's component.
In this study, the analysis with the STRUCTURE program showed that Turkish water buffalo populations were grouped into two major lineages when K = 2 ( Figure 4B). Cluster I included IST, BAL, BUR, TEK, DUZ, SVS, KAY, BIT, MUS and DYB populations, and Cluster II gathered GIR, AMS, TOK, COR and SIN populations, while other water buffalo populations (AFY and SAM) appeared to be the contact zone between both clusters, as individuals had mixed lineages. The STRUCTURE analysis results support the neighbornet dendrogram results, as well as F ST , PCA and genetic distance results. Our results provide a broad perspective on the extant genetic variation and population structure of Turkish water buffalo populations.
Identification of within and between populations genetic diversity is a prerequisite for well-structured and sustainable animal breeding and conservation programs. The study, which was carried out by using 20 microsatellite markers, showed that within population genetic diversity was higher than between population diversity. This situation can be seen as an opportunity in terms of breeding programs and genetic conservation programs for these populations. Hence, it can be concluded that the Turkish water buffalo population possessed a considerable amount of genetic diversity due to the low pressure of artificial selection and the possibility of random mating; however, Turkish buffalo populations require a scientific production system in order to improve the production, without losing the significant genetic structure of these economically important animals.

Conclusions
In this paper, we have provided noticeably powerful data on genetic diversity and population structure of Turkish water buffalo populations, which might be helpful for similar studies. Our results suggested the relatively low but statistically significant genetic diversity of 17 Turkish water buffalo populations and brought insight into the structure of the analyzed populations. This study is the first assessment of the molecular genetic diversity of six geographical regions' water buffalo populations in Turkey. The findings on the genetic structure of Turkish water buffalo populations in the present study will have significant implications in formulating the future strategies for conservation and breeding programs.