Discussion

The cattle commercial panel used to genotype buffaloes allowed us to infer the homology and the divergence between the species. However, a clear bias in the comparative analyses between the genetic groups was observed because the commercial panel is composed of polymorphic markers of B. taurus, especially with taurine breeds (Illumina BovineHD Genotyping BeadChip Data Sheet). This fact was most evident in the analyses of Hierarchical Clusters, because as the MAF limit values decreased, the bias effect of the selection of markers also decreased, resulting in genetic clustering that coincided with the expected zoological classification between the Bubalus and the species of Bos genus (Vaughan 1986). The genetic similarity/divergence analyses, in fact, inferred about the distance of B. taurus taurus to other bovids and the inference across other groups was prejudiced. For example, PCA distances indicated that B. taurus indicus was closer to buffalo than to taurine cattle breeds due to marker selection bias. A similar effect was found in cervidae genotyped with cattle SNPs, indicating that Bos taurus indicus were closer to cervids than to Bos taurus taurus (Kasarda et al. 2015). Taurine was used as reference to avoid erroneous interpretations of the results, so the distances were taken from them, not between other groups. However, PC1 fully discriminated buffalo group from the taurine group, so buffaloes can also be considered a reference for comparison (similarity of the genetic groups from buffaloes). It implies that the interpretation of results apart from buffaloes contains little bias and it can be considered valid. It is an important result for the forthcoming chromosomal rearrangements proposed.

The PCA allowed us to observe the discrimination within the genus Bos and between the genera Bos and Bubalus, but not the discrimination between B. depressicornis and Bubalus. Moreover, when Bos taurus taurus was assumed as a reference for comparison, the furthest group were the buffalo, followed by Bos javanicus and the Asian B. taurus indicus and, finally, the African breed. Hierarchical Clusters (Figure 2 - HC0.5 and HC0.3), and admixture analysis (Figure 3 - k = 2) revealed greater genomic similarities between buffaloes and Asian cattle than to European cattle. These results indicate that buffaloes may actually be closer to B. javanicus followed by B. taurus indicus and Bos taurus taurus. This relation can be attributed to the analogous evolution process of these groups, such as the adaptation of the cattle and buffalo to the same Indian environment (Hoffpauir 1982).

The geographical origins of buffaloes and B. javanicus could contribute to explaining the genetic similarity because both of them originated in South-East Asia (Lau et al. 1998; Ropiquet et al. 2008). After its origins in Southeast Asia, buffaloes spread to northern China and western India, the natural habitat of zebu (Lau et al. 1998). Colli et al. (2018) suggest that two independent domestication events occurred in Indo-Pakistani region for water buffalo and close to the China/Indochina border for swamp buffalo. However, there are two current explanations of zebu origins in India (MacHugh et al. 1997). The first has Auroque (Bos primigenius primigenius) as a common ancestor, which despite being Asian, it first originated Bos taurus in Europe and after, it diversified into the zebu in India. The second and the most accepted is that zebu developed independently in South Asia from another subspecies of the Auroque, the Bos primigenius namadicus.

The high r² and the low haplotypic diversity observed in the centromeric regions of buffalo submetacentric chromosomes are strong evidence of a fusion between the chromosomes of a common ancestor that had 29 acrocentric autosomes (2n = 60), and not a rupture of the submetacentric chromosomes in this supposed ancestral to originate the 29 acrocentric chromosomes (Iannuzzi et al. 2009). The hypothesis of a presumed ancestral with 29 autosomes and the greater genomic similarity between Asian Bovidae corroborate the theory that the common ancestor to these bovids was possibly Bos primigenius namadicus.

Chromosome structural evidence also indicate similarity across species. For example, B. javanicus (2n = 56) have two Robertsonian translocations, while buffaloes have 5 to 6 (2n = 50 and 2n = 48) (Iannuzzi et al. 2003; Ropiquet et al. 2008). Although translocations are different among these species, they always have involved orthologs to BTA 1, 2, 28 and 29 (Ropiquet et al. 2008). Since the translocation process always affects the same groups of chromosomes and occurs between centromeres, the predisposition for translocation may not be only random and structural but it may also play a role in their composition, such as the number of tandem repeats of sequences that may cause instability and consequent centromere displacement during meiosis (Berend et al. 2003; Purgato et al. 2015). Other evidence is that the most common chromosomal aberration in cattle (intraspecific) is the translocation of 1:29 (De Lorenzi et al. 2012; Yimer and Rosnina 2014).

In this way, the ancestor 2n = 60 (common to both groups) possibly had variations/mutations that predisposed translocations with a frequency that allowed the appearance of different species by reproductive isolation. Another evidence about the similarity between zebu and buffalos is the presence of acrocentric Y sex chromosome, unlike the European cattle that have metacentric Y sex chromosome (Iannuzzi et al. 2001; Yindee et al. 2010). Although these results suggest that Asian cattle may have a common origin independent from the European, the possibilities of crossbreeding between the ancestors of these groups at different periods are plausible in both theories, making it more difficult to analysis the evidence.

The similarity between Asian Bovidae observed can also be attributed to natural selection processes that converged to generate similar functions among the species because these are subject to the same tropical environment and this may have favored the predominance of certain alleles (Vickrey et al. 2015). In addition, the similarities obtained with higher frequency alleles are also possibly related to the convergence of artificial selection for domestication and production.

The introgression analysis allowed to discriminate between Bubalus species and B. depressicornis, when the number of clusters was equal to 8, indicating a slight B. depressicornis introgression in the water buffalo. This introgression possibly originated from the centromeric regions of the four submetacentric chromosomes common to these species, considering that these regions are more conserved/preserved than other chromosomal regions (Gallagher et al. 1999; Schneider et al. 2016). The buffaloes and Bos javanicus were clustered together when cluster number was equal to 2. Also, the Asian zebu had a greater buffalo and B. javanicus introgression when compared to the European animals, unlike the African zebu that presented greater introgression of the European clustering. This introgression may possibly have resulted from the fact that Boran is a recent breed, originated from mixed taurine breeds (Rege 1999). The possible reasons for the Indian zebu to present greater Asian Bovidae introgression are, in addition to crossbreeding, the homology of the species and the possible convergent evolution of these groups as discussed above. B. javanicus was the last to differentiate from buffaloes when the number of clusters increased, this similarity was observed in the other analyses. B. javanicus also had buffalo introgressions and a slightly higher zebu proportion when the number of clusters was lower than 5.

Decker et al. (2014) reported zebu introgression into B. javanicus and considered these animals as belonging to the same cluster, in a study considering B. javanicus, Eurasian and African B. taurus, for cluster number equal to 3. In our study, when the cluster number was greater than or equal to 5, no relevant introgressions were observed between the Bubalis and Bos genera and the distinction between the taurines breeds was observed when the number of cluster was lower to 4, while more than 5 clusters were necessary to separate zebus. This result strongly indicates that crosses between the ancestors of these Asian Bovidae of both genera occurred, as it was previously considered. In addition, although no algorithm was used to select the best number of clusters, eight clusters (k = 8) showed to be the best option for these analyses, since these correctly discriminated all genetic groups (species and breeds) in consonance to a biological perspective.

The r² levels were similar to reported to buffalo species by Deng et al. (2018) and Mokhber et al. (2019). The buffalo specific SNP panel enabled a better inference of LD among markers, especially between chromosomes than the one resulted from the cattle SNP panel. The LD between the markers of the homologous structures allowed not only to infer which structures were related such as the FISH (Fluorescent in Situ Hybridization) cytogenetic probes (Iannuzzi et al., 2003), but also to orient their positioning such as in the original buffalo assembly (Low et al. 2019), and possible evolution.

The use of LD to verify the relationship of the homologous structures still had the advantage of being easily obtained by panels of markers that are used in large scale for other genomic studies, without the necessity to design specific probes with single functionality. Comparing to sequencing data, the LD approach allows to use many individuals in the analysis at a low-cost; whereas in sequencing, only one or few individuals are used (despite this method allows greater number of variants). The disadvantage is that LD estimates can be affected by several factors such as history and population structure, effective population size, sample size, markers density and distribution, and strict filtering of SNPs (Hayes et al., 2003; Yan et al., 2009; Bohmanova et al., 2010).

LD methodology also detected a correlation of a 0.7 Mb region of BTA21 with BTA 1 and 27. Although it appears to be a translocation in buffaloes, Utsunomiya et al. (2016) reported that this region is a possible assembling error of the UMD3.0 bovine reference.

Our rearrangement showed high agreement with UOA_WB_1 for SNP order and chromosome length according to the correlations and linear regression analyses. These results were expected since, according to Low et al. (2019), the UOA_WB_1 was scaffolded in an order conserve synteny with the homologous Bos taurus genome (UMD3.1). These authors also observed good agreement for the chromosome sizes and proportion of sequences aligned to corresponding homologous B. taurus chromosomes, comparing UOA_WB_1 and UMD3.1.

Some misassembled SNPs comparing UOA_WB_1 and UMD3.1 were previously reported by Utsunomiya et al. (2016). Additional possible misassembled SNPs were also detected in our study. Possible misassembled SNPs can be detected in new regions due to the assessment bias of the SNP panels as well as the different polymorphisms in the buffalo genome. In general, all the results confirm that the buffalo genome assembly was more accurate than the rearrangement of the UMD3.1 to predict the right marker position. However, the rearrangement allowed to include many SNPs given as unknown position (1,036) in UOA_WB_1, despite a few of them were considered as potential misplaced SNPs. The main difference between the genome assemblies is not related to segment translocations, but is due to the quality of sequence methods, since the UMD3.1 used Sanger (https://www.ncbi.nlm.nih.gov/assembly/GCF_000003055.6/) and while UOA_WB_1 used PacBio technology (https://www.ncbi.nlm.nih.gov/assembly/GCF_003121395.1/). These results indicate that the improvement in sequence technologies used for cattle as well as new cattle assemblies, would improve the quality of the rearrangement. Although we have not tested, the results also suggest that the improvement in cattle assembly can also help any specific buffalo assemblies, mainly in misassembled regions.

The r² estimates were similar at the telomeres of the five studied chromosomes and regions with intense LD depression were not observed. This result indicated that there were no large deletions (> 4 Mb) in these buffalo chromosomes compared to the cattle homologs, as was observed by Low et al. (2019) with sequence data. However, significant differences were observed when comparing the centromeric and telomeric regions. It indicates a fusion process in common ancestors to originate buffalo submetacentric chromosomes, instead of the rupture of the ancestral chromosomes to generate the acrocentric chromosomes in cattle. The BBU5 assembled with UAO_WB_1 had no high LD in the centromere region, as we observed in our rearrangement. This high LD in the expected region indicates that the rearrangement is the better way to study LD levels in the submetacentric chromosomes.

High LD levels in the centromeric and pericentromeric regions have been observed in plant and animal species (Smith et al. 2005; Shi et al. 2010; Talbert and Henikoff 2010). The possible explanation is that the probability of crossing-over is very low in these regions or that there is no crossing-over in the centromere (Shi et al. 2010; Talbert and Henikoff 2010). Higher LD levels and less diversity were observed in the aforementioned studies, when compared to the present study. It possibly happens because the translocations in buffaloes are more recent. The species has only five submetacentric chromosomes, while the other species have all the autosomes as acrocentrics (Smith et al. 2005; Shi et al. 2010). Therefore, the low recombination rate and the reduced diversity in the centromere may serve as indicators of the previous translocation processes.

The different LD patterns observed in the centromeres of the buffalo submetracentric chromosomes may be due to the order of occurrence of the translocations during evolution time. The ascertainment bias of the markers in this region and the current location of the centromere may also interfere in the LD behavior and haplotype diversity. The ascertainment bias effect was minimized by adjusting the physical distance, using a LD decay function. Five LD peak patterns were observed on the submetacentric chromosomes of water buffaloes:

1. High disequilibrium (r² ≈1) in the centromeric region and small LD extension in the pericentromeric region (BBU 2);

2. High disequilibrium in the centromeric region (r² ≈1) and slightly higher LD extension in the pericentromeric region (BBU 3);

3. Medium disequilibrium in the centromeric region (r² ≈0.4) and greater LD extension region (BBU 4);

4. Low disequilibrium (r² ≈0.2) and high extension of the LD region (BBU 1);

5. Extension of LD very high in the pericentromeric region and, medium, in the centromeric region (BBU 5).

Considering that, possibly, the last chromosome to be formed was BBU5, because it is the only one of the 5 submetacentric of B. bubalis that is missing in B. depressicornis (Gallagher et al. 1999), the LD extension in the centromeric region may be related to recent events in these regions. Likewise, Hayes et al. (Hayes et al. 2003) also correlated the high extension LD between markers to recent events in the history of a population. Therefore, the chronological order of the translocation events was possibly BBU2, BBU3, BBU4, BBU1, and BBU5. Thus, a hypothetical model on the evolution stages was assumed on the LD behavior and diversity in the chromosome arms, pericentromeric and centromeric regions (Figure 9).

• Download figure
• Open in new tab
• Download powerpoint

Figure 9

Hypothetical model on linkage disequilibrium (r²) level behavior in chromosome arms (ARM), pericentromeric (PE) and centromeric (CE) regions throughout generations after the translocation process forms a new chromosome. The linkage disequilibrium standards were not observed in stages 1 and 5.

Initially, the r² peak tends to be low in the centromeric regions due to the low divergence in this region, increasing over the generations with the increase in diversity (stages 1 and 2). After increasing the diversity in the region, LD (stage 3) also increases. Subsequently, the loci tend to enter in equilibrium over generations, thus decreasing the LD, but much slower than the other chromosomal regions. This equilibrium starts from the outermost centromeric regions toward the center (stage 4). At the end, the low recombination rate of the centromere is balanced with the level of diversity and the region has higher averages than the other chromosome regions, however at a below level of the previous stage (stage 5).

This sequence was conceived without considering the events of genetic bottlenecks that tend to increase LD in both regions, which tend to be dissipated more quickly in the regions away from the centromeres.

Diversity is fundamental to obtain LD estimates while diversity is propagated in the population through recombination. After the occurrence of a translocation, the telomeres of the first new submetacentric chromosome (in the first individual) can still suffer crossover with the non-translocated homologous chromosomes, although at a much lower rate, while the centromere remains with diversity and r² equal to 0. The possibility of reproduction among buffalo subspecies is indicative of this effect. After mutations in the centromere regions (diversity), the estimated LD is already higher due to the lower recombination rate.

The diversity analysis showed that BBU5 had less diversified pericentromeric and centromeric regions besides a greater pericentromeric extension of this conservation. BBU1 also presented a similar profile for these regions, however with larger magnitudes for haplotype diversity. These characteristics also indicate that these chromosomes were the last to be fused, agreeing with the results observed in the LD analyses. With UAO_WB_1 we observed the lowest haplotype diversity levels in extra-centromeric regions for BBU4 and BBU5. Thus, the rearranged genome cattle assembly was considered more accurate and used as model to elucidate evolution process in the centromeric regions.

The r²a and PERC values obtained for the buffalo submetacentric chromosomes were close to the variation obtained for Holstein and Angus (Mulder et al. 2012; Berry et al. 2014) and for Nellore and Gir (Carvalheiro et al. 2014; Boison et al. 2015), for imputations using 10% of the density of the original marker. This result indicated that the cattle reference for arranging the markers for each orthologous structure could also be used for the buffalo genome. In fact, imputation accuracy is not expected to increase significantly by the rearrangement of the adjacent regions, but to decrease abruptly when the rearrangement is inadequate, reflecting the effect of the misplaced position on the formation of the haplotype blocks (Utsunomiya et al. 2016).

Thus, BBU 1 was observed to have the greatest increase in imputation accuracy, especially by correcting the bovine reference between the BTA 21 and BTA27 regions. Regarding the other chromosomes studied, a marked reduction of the imputation accuracy in the centromeric and pericentromeric regions would be observed if the structures had been rearranged inadequately. However, a slight increase in accuracy was observed, indicating a small improvement in the construction of haplotype blocks. This result indicated that the proposed rearrangement was adequate and could be used as the approximate genomic reference for buffalo in studies that are dependent on this coordinates. The better performance of the rearrangements when compared to UAO_WB_1 in many scenarios as well as the small trend of this assembly to decrease the imputation accuracy, according to the number of markers, indicated that rearrangement is the best option for studies in the centromeres, in agreement to the LD and haplotype diversity scan analysis.

The water buffalo has genomic homology to cattle, and its five submetacentric chromosomes come from chromosomal fusions of a common ancestor closer to Asian Bovidae than to taurine. This genomic homology between the species enabled us to use the reference cattle assembly to recreate a buffalo genomic reference by rearranging the submetacentric chromosomes. When using the bovine genomic reference, the rearrangement of the buffalo submetacentric chromosomes could be done by SNP BTA calculations. The centromere of these rearranged chromosomes had the expected profile for the centromere of submetacentric chromosomes exhibiting high linkage disequilibrium and low haplotype diversity, enabling us to hypothesize about evolutionary-genetic events. The imputation accuracy revealed that the proposed rearrangement for these chromosomes was adequate and that the approximate genomic reference for water buffalo can be used in scientific studies. The SNP-based approach is a cheap and easy solution to transfer any cattle assembly technology to buffalo, beyond to help place and order many unplaced SNP/segments in a buffalo assembly. Moreover, the proposed approach can be applied to elucidate other chromosomal rearrangement events in other species and possibly better understand the evolution relationship among them.