Genomic study of population structure and phylogenetic relationship of Khuzestani buffaloes

The science of phylogenetics, and specially the subfield of molecular systematics, has grown exponentially by amount of publications, general interest and amount of available genetic data. Modern phylogenomic studies use large genomic and transcriptomic resources. however a comprehensive molecular phylogeny of animals, including the newest types of data for all phyla, remains elusive.
Genome sequences and SNP chips are now accessible for many of species across the animal phylogeny, conveying key features of animal genome evolution into sharper emphasis. The field of animal evolutionary genomics has focused on recognizing the diversity genomic features, rebuilding the history of evolutionary variations in animal genomes and testing hypotheses about the evolutionary relationships of animals. One of the most important buffaloes in Iran are Khuzestani buffaloes, which have a wide distribution. Given the fact that they may be divided into several ecotypes and sub-populations, these may cause problems in genomic wide association study. One of the challenges that always threatens the accuracy of genomic association and genetic variation studies is the existence of sub-populations within the populations. The probable relationship between genetic groups among Khuzestani buffaloes by whole genome markers has not been studied. This study is an effort to understand and provide a detailed insight into the population structure and genomic phylogeny of Khuzestani buffaloes.
Material and Blood or hair samples from 121 Khuzestani buffaloes from different regions of Khuzestan province including Ahvaz (n=28), Dezfoul (n=19), Shadegan (n=6), Shosh (n=4), Shushtar (n=9), Dashtazadegan (n=24) and Kermanshah (n=8) samples were collected. The samples were chosen from herds that were registered by Iranian Animal Breeding Center (IABC) and were under milk recording system. DNA of the samples were isolated and sent to Parco Tecnologico Padano in Itally. All samples were genotyped by Axiom@Buffalo 90K bead array of Affymetrix Company. PLINK software was used for pre-processing of raw genotype data based on minor allele frequency (MAF) ≤0.01, Hardy Weinberg disequilibrium (p-value ≤ 10-6) and call rate ≤ 0.05. Thereafter, a total of 64,709 SNPs remaind for the analysis. The genetic distance for each pair of individuals was calculated using R software (ape package). The Neighbor-joining algorithm was used to plot trees based on those distances. Also, the phylogenetic diagram is derived based on the distance matrix and using R software (Phyclust package). To identify the optimal number of clusters, K-means were sequentially increased with K values, and different clustering responses were compared using Bayesian Information Criteria (BIC). FST values (theta) for autosomal chromosomes and counties were calculated by unbiased method. Then, the average genetic distance of each area was calculated with other sampling areas. Principal Component Analysis (PCA) and Discriminant Analysis of Principal Components (DAPC) were performed for genetic clustering of individuals. Adegenet R package was used for PCA and DAPC analysis.
In our study Unrooted Neighbor-Joining Tree Plot grouped the individuals mainly in one cluster depending on the origin of group or subpopulation. Dezfoul, Shushtar and Dashtazadegan individuals formed their individual clusters and due to genetic relatedness, Ahvaz and Kermanshah clusters were formed almost overlapping with other. Along with the closeness of Dezfoul among themselves, they also showed closeness to Shosh. Kermanshah group lied between Ahvaz, Dashtazadegan and Dezfoul buffaloes.
Results obtained from DAPC analysis confirmed with those obtained from PCA analysis. Individuals were correctly assigned to their respective clusters. Based on phylogenic result, it is evident that buffaloes of Ahwaz are affected by buffaloes of other province’s regions, because they are in the center. Buffaloes of Shushtar, Dezful and Dasht-E-azadegan have lower emigration rate in comparison with Ahwaz buffaloes. This result also confirmed by principal component analysis. To survey the genomic structure and patterns of Khuzestani buffaloe populations, the Fst values (unbiased Theta) for every SNPs were calculated. This criterion reflects the distribution of diversity between groups and within the group. Its high level indicates that major part of heterozygosity is due to differences in the groups. In other hand, the low level of it, suggests that the diversity is due to the individuals within the group, and the populations or groups are associated. The highest value of Fst was obtained for buffaloes of Dashte-E-azadegan (0.0159) and Shadegan (0.0147) because of its special geographical condition, and the lowest value of Fst was observed for the buffaloes of Ahwaz (0.0087) and Shosh (0.0086). The results of statistical tests of population genetic showed a high genetic diversity in buffaloes of Khuzestan. Graphical view of clustering and assignment test suggested at least three subpopulations for Khuzestani buffaloes. According to the overall results and existence of great SNP diversity, it is possible to move toward construction of base population based on this results.
In this study a genome-wide insight into the genomic phylogeny of Khuzestani buffaloes was provided. The results of nucleotide genetic diversity and genetic distance between buffaloes are worthy of attention in subsequent studies such as genomic wide association study, genomic selection and formation of Base populations of Khuzestani buffalo.