Genetic diversity of BoLA-DRB3 in South American Zebu cattle populations

Background Bovine leukocyte antigens (BoLAs) are used extensively as markers of disease and immunological traits in cattle. However, until now, characterization of BoLA gene polymorphisms in Zebu breeds using high resolution typing methods has been poor. Here, we used a polymerase chain reaction sequence-based typing (PCR-SBT) method to sequence exon 2 of the BoLA class II DRB3 gene from 421 cattle (116 Bolivian Nellore, 110 Bolivian Gir, and 195 Peruvian Nellore-Brahman). Data from 1416 Taurine and Zebu samples were also included in the analysis. Results We identified 46 previously reported alleles and no novel variants. Of note, 1/3 of the alleles were detected only in Zebu cattle. Comparison of the degree of genetic variability at the population and sequence levels with genetic distance in the three above mentioned breeds and nine previously reported breeds revealed that Zebu breeds had a gene diversity score higher than 0.86, a nucleotide diversity score higher than 0.06, and a mean number of pairwise differences greater than 16, being similar to those estimated for other cattle breeds. A neutrality test revealed that only Nellore-Brahman cattle showed the even gene frequency distribution expected under a balanced selection scenario. The FST index and the exact G test showed significant differences across all cattle populations (FST = 0.057; p <  0.001). Neighbor-joining trees and principal component analysis identified two major clusters: one comprising mainly European Taurine breeds and a second comprising Zebu breeds. This is consistent with the historical and geographical origin of these breeds. Some of these differences may be explained by variation of amino acid motifs at antigen-binding sites. Conclusions The results presented herein show that the historical divergence between Taurine and Zebu cattle breeds is a result of origin, selection, and adaptation events, which would explain the observed differences in BoLA-DRB3 gene diversity between the two major bovine types. This allelic information will be important for investigating the relationship between the major histocompatibility complex and disease, and contribute to an ongoing effort to catalog bovine MHC allele frequencies according to breed and location. Electronic supplementary material The online version of this article (10.1186/s12863-018-0618-7) contains supplementary material, which is available to authorized users.


Background
Domestic bovines comprise two major breed types: humpless (Bos taurus) and humped (Bos indicus or Bos taurus indicus). Based on archeological and molecular evidence, it is estimated that the Taurine and Zebu lineages divergence occurred from a common ancestor between 330,000 and 2 million years ago, respectively [1][2][3][4]. Since then, the two types of cattle have accumulated many genetic variations resulting in a high genetic distance; this has contributed to highly differentiated phenotypes and adaptation. Currently, Zebu breeds are found worldwide in tropical and subtropical regions; indeed, beef production in South America's warm regions is based mainly on Zebu cattle (Nellore, Brahman, and Gir breeds).
The bovine MHC, named bovine leukocyte antigen (BoLA), comprises three DRB class II copies, with the BoLA-DRB3 gene being the most highly expressed and polymorphic [5]. Polymorphic sites within the BoLA-DRB3 gene are located mainly at peptide-binding sites; these polymorphisms are mostly maintained by balancing selection [6]. Several lines of evidence show that polymorphisms influence both the magnitude and epitope specificity of antigen-specific T cell responses to infectious diseases [7]. In this sense, the BoLA-DRB3 study has used SBT or NGS assays to examine genetic diversity of BoLA-DRB3 in large samples. With this in mind, the aim of the present study was to examine the genetic diversity of the BoLA-DRB3 gene in three Zebu breeds from Bolivia and Peru, and to compare the results with those reported for Taurine breeds. The results may contribute to an ongoing effort to catalog bovine MHC allele frequency according to breed and location.
Cattle were handled by veterinarians from local farms, RIKEN, and Universidad Autónoma Gabriel René Moreno, in strict accordance with good animal practice following the Universidad Austral de Chile Institutional guidelines. Veterinarian Collaborator explained about the sample collection to farmer in verbal way, and their future use in the research project. This study was approved by Committees on the Ethics of Animals for Research at the National from University of LA PLATA (Argentine: Certificate date May 26th, 2014), Universidad Autónoma Gabriel Moreno (Bolivia, and Universidad Austral de Chile (Certificate No. 153-2014). BoLA-DRB3 typing using PCR-sequence based typing (SBT) BoLA-DRB3 alleles were genotyped using PCR-SBT. Briefly, exon 2 of BoLA-DRB3 was amplified by PCR as described by Takeshima et al. [19], using primers DRB3FRW and DRB3REV [28]. The PCR fragments were purified using an ExoSAP-IT PCR Product Purification Kit (USB Corp., Cleveland, OH) and sequenced using the ABI PRISM BigDye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems, Foster City, CA). Raw sequence data were analyzed using Assign 400ATF ver. 1. 0.2.41 software (Conexio Genomics, Fremantle, Australia).

Measures of genetic variability
Allele frequencies and the number of alleles (N a ) were obtained by direct counting, and 95% confidence intervals (CI) for allele frequencies were calculated using the binomial distribution implemented in R with the binom. confint function (http://cran.r-project.org/web/packages/ binom/). The distribution of alleles across major groups of breeds was represented by a Venn plot using the R package 'VennDiagram'. The observed (H o ) and expected unbiased (H e ) heterozygosity for the BoLA-DRB3 locus were estimated as described previously [29] using ARLEQUIN 3.5 software for population genetic analyses [30]. Potential deviations from Hardy-Weinberg equation (HWE) were estimated for each breed by F IS statistics [31] using the exact test included in GENEPOP 4.0 [32]. The Ewens-Watterson-Slatkin exact test for neutrality was estimated as described by Slatkin [33] and implemented in ARLEQUIN 3.5 software.

Genetic structure and differentiation
Genetic structure and genetic differentiation among breeds were assessed using Wright's F statistics (F ST ) and the previously described variance-based method [31]. These parameters were estimated using ARLE-QUIN 3.5 and GENEPOP 4.0. F ST values were represented graphically using the pairFstMatrix.r function in the R statistical environment.

Principal component analysis (PCA)
To condense the genetic variation at the BoLA-DRB3 locus, allele frequencies were used to perform principal component analysis (PCA) as described previously [34] using PAST software [35]. Nei's standard genetic distance (Ds) [36] and Nei's unbiased genetic distance (D A ) [37] were calculated from allele frequencies, and cluster analysis was performed using the UPGMA [38] and the neighbor-joining (NJ) [39] algorithms. For these analyses, only pure breeds and a Nellore-Brahman mixed population were included due to limitations inherent in the dendrograms. CIs for the groupings were estimated by bootstrap re-sampling of data using 1000 replications.
Genetic distances and trees were calculated using POP-ULATIONS 1.2.28 software [40], and trees were constructed using TREEVIEW [41].
Genetic diversity at the nucleotide and amino acid sequence levels Nucleotide diversity (π) and pairwise comparison of nucleotide substitutions between alleles were calculated using ARLEQUIN 3.5. The mean number of nonsynonymous (dN) and synonymous (dS) nucleotide substitutions per site was estimated for each pair using the Jukes-Cantor's formula as described by Nei and Gojobori [42]. These parameters were estimated using MEGA 5 [43].

Results
Distribution of BoLA-DRB3 alleles among analyzed zebu cattle breeds PCR-SBT genotyping allowed to identify 46 previously reported BoLA-DRB3 alleles, while no new variants were observed ( Table 2) and all genotyping result were shown in additional table (Additional file 1). The N a ranged from 19 BoLA-DRB3 alleles in Gir to 33 variants in Nellore-Brahman (Table 2). Of these, six alleles reported in the IPD-MHC public database (https://www.ebi.ac.uk/ ipd/mhc/) were from Bos indicus and 26 alleles were from B. taurus; 14 alleles were shared among these two major bovine types. Furthermore, a Venn diagram was constructed using data obtained in this study and from previous reports [12][13][14][15][16][17][18][19][20][21][22]; these data, that included 14 breeds and cross-breeds, were then used to evaluate the distribution of the BoLA-DRB3 allele ( Fig. 1.). To achieve this, data were grouped in terms of geographical origin, i.e., British, European Continental, American Creole, and Zebu cattle breeds. The plot clearly shows that the Zebu group, that comprised Bolivian Gir, Bolivian Nellore, Peruvian Nellore-Brahman and Philippine Brahman breeds, harbors 21 out of the 90 alleles identified in the four cattle groups (Fig. 1); these 21 alleles represent about 30% of the 66 alleles detected in the Zebu group. The N a within each Zebu population with frequencies > 5% varied; from five in Bolivian Nellore and Peruvian Nellore-Brahman to ten in Bolivian Gir. Some of the high-frequency (> 5%) alleles were present in at least two out of three Zebu populations; these alleles were BoLA-DRB3*2801, *2201, *3101, and *3601 ( Table 2). These common alleles account for 65.95% of all alleles sampled in Bolivian Nellore, 67.18% in Peruvian Nellore-Brahman, and 85.00% in Bolivian Gir (Fig. 2).
Nucleotide and amino acid diversity of BoLA-DRB3 alleles in analyzed zebu breeds Genetic diversity was examined at the DNA level using two parameters: π and the number of pairwise differences (NPD). The results obtained were then compared Maillard et al., [26] Takeshima et al., [18] with those previously reported for ten other breeds ( Table 3). The π ranged from 0.068 in the Nellore × Brahman population to 0.078 in the Bolivian Gir population, while the mean NPD varied from 16.91 in the Nellore × Brahman to 19.43 in Bolivian Gir. The 13 cattle breeds analyzed by PCR-SBT varied in terms of the number and composition of BoLA-DRB3 alleles. However, with the exception of Shorthorn and Native Philippine cattle, the values observed for the π and NPD indices were similar to those estimated for other cattle breeds (Table 3).
We also calculated the average number for dn and ds in Bolivian Nellore, Peruvian Nellore × Brahman, and Bolivian Gir. As shown in Table 4, the values obtained for these populations (0.037 to 0.042 for ds and 0.103 to 0.108 for dn) when analyzing the whole of exon 2 suggested a dn/ds difference range for this genome region of 0.058 to 0.066. The dn/ds differences were more prominent when only the antigen-binding sites (ABS) were analyzed (dn/ds > 0.132; Table 3).
Gene diversity, HWE, and neutrality testing of BoLA-DRB3 in the analyzed zebu populations Genetic diversity within the three analyzed Zebu populations was estimated according to the N a and the observed and expected heterozygosities. As mentioned above, the N a ranged from 19 in Bolivian Gir to 33 in Peruvian Nellore-Brahman, while the values of expected heterozygosity ranged from 0.86 in Nellore-Brahman to 0.92 in Bolivian Gir ( Table 4). The observed heterozygosity in Peruvian Nellore-Brahman and Bolivian Gir varied from 0.76 to 0.88, respectively (Table 4). Regarding HWE, two out of three populations were in equilibrium, while Nellore × Brahman showed a significant deviation from the theoretical values (Table 4). This deviation can be explained by a significant excess of homozygous animals (F IS = 0.113, p value = 0.033). To assess balancing selection, Slatkin's exact neutrality test was performed ( Table 4). The BoLA-DRB3 gene frequency profile in Gir cattle showed an even distribution (p = 0.008) despite the finding that its genotype proportions did not deviate from the theoretical proportions. This even gene frequency distribution is consistent with the theoretical proportion expected under recent balancing selection as opposed to positive or neutral selection (p > 0.025). Conversely, the Nellore and Nellore × Brahman gene frequency distributions were compatible with neutral selection (Table 4).

Genetic structure and differentiation of BoLA-DRB3 in zebu breeds
To examine population genetic structure and differentiation levels among cattle breeds, we calculated two parameters: the F ST index and the exact G test. There were significant differences in the F ST parameter across all cattle populations (F ST = 0.057; p value ≤0.0001); pairwise comparison ranged from 0.022 (between Philippine Brahman and Philippine Native) to 0.133 (between Peruvian Nellore × Brahman and Chilean Hereford) (Fig. 3). However, there was around 5% differentiation among Zebu populations. The exact G test for differentiation indicated that, on average, differences in allele frequency distribution among populations were highly significant (exact p value ≤0.0001). Furthermore, when all pairs of breeds were tested, significant results (p value ≤0. 0001) were also obtained for all comparisons, indicating a significant level of genetic structure for BoLA-DRB3 among breeds. In conclusion, when Zebu breeds were included in the analysis, the average F ST value was 1 to 5% higher than that obtained when analyzing taurine breeds alone [12] or a population of the same breeds [22], respectively.

Genetic differentiation of BoLA-DRB3 alleles in zebu and taurine breeds
To assess the genetic differentiation of BoLA-DRB3 alleles between Zebu cattle and nine breeds previously studied by SBT, we performed two types of analysis: dendrograms and PCA. First, BoLA-DRB3 allele frequencies were used to generate Nei's D A and D S genetic distances matrices; data from Bolivian Nellore, Bolivian Gir, Peruvian Brahman × Nellore, and nine previously reported breeds (Japanese Jersey, Japanese Shorthorn, Japanese Holstein, Japanese Black, Chilean Hereford, Philippine Native, Philippine Brahman, Yacumeño, and Hartón del Valle; [12][13][14][15][16][17][18][19][20][21][22]) were used. Dendrograms based on distance were constructed using UPGMA and NJ algorithms. NJ cluster analysis based on Nei' D A or D S genetic distances showed two major clusters: one comprising mainly European taurine breeds and Yacumeño Creole cattle and one comprising Zebu breeds, Hartón del Valle, and Philippine Native breeds (Fig. 4a). Japanese Jersey or Chilean Hereford was located outside of these clusters. UPGMA cluster analysis showed a similar trend, with Jersey and Hereford breeds located in separated clades (data not shown).
Second, we used allele frequencies of BoLA-DRB3 from the 12 breeds mentioned above to perform PCA. The results are reported in Fig. 4b, which shows the first and second principal components (PCs) for BoLA-DRB3 gene frequency. The first two components accounted cumulatively for 44.68% of the variability in the data. The first PC accounted for 28.34% of the total variance, and showed a clear differentiation pattern between Zebu (positive values) and Taurine breeds (negative values), whereas native breeds were located near the origin of the plot. The first PC was determined mainly by 26   BoLA class II molecule binds various peptides derived from antigens via five antigen binding pockets named pocket 1, pocket 4, pocket 6, pocket 7 and pocket 9 [44]. To assess whether observed differences in allelic frequency between Taurine and Zebu populations is reflected within amino acid motifs in each pocket, we analyzed protein pockets implicated in antigen-binding function by PCA. With respect to taurine breeds, Zebu breeds were clearly differentiated in component 2 of pocket 1 due to a higher frequency of GVFT and VGFT and a lower frequency of VVFT and GMFT amino acid motifs (Fig. 5). Component 1 of pockets 4 and 7 allowed differentiation of Zebu (positive values) from Taurine breeds (negative values). GLDEREY and GLDRREV amino acid motifs in pocket 4, and DYWIR, DFWFR, and EYWIR amino acid motifs in pocket 7, provided the major contribution and explain the distribution observed in PCA plots. With the exception of Philippine Brahman, Zebu breeds showed a similar trend in pocket 9, which was accounted for by QYS, QFA, and END amino acid motifs. Finally, component 1 of pocket 6 only discriminated Nellore from Nellore crossbred by means of the HH motif.

Discussion
Domestic cattle were derived from the now-extinct wild aurochs (Bos primigenius) [45]. Zooarcheological evidence revealed that this wild species was widespread throughout Eurasia and North Africa, even though, events of bovine domestication occurred only few times in the human history and in precise geographic regions (called domestication center). In this sense, different local wild auroch subspecies populations would have been domesticated in five domestication centers: three Fig. 3 a Genetic distance between pairs of populations estimated by Wright's F statistics (F ST ) (below) and Nei's D A distance (above). b Graphical representation of calculated F ST values between population pairs using an R-function: pairFstMatrix.r. HeCh = Chilean Hereford; HoJa = Japanese Holstein; ShoJa = Japanese Shorthorn; JeJa = Japanese Jersey; WaJa = Japanese Black; Ya = Yacumeño; HV = Hartón del Valle; NaPh = Philippine Native; BrPh = Philippine Brahman; GirBo = Bolivian Gir; NeBo = Bolivian Nellore; and NexBrPe = Peruvian Nellore × Brahman for B. taurus (southwest Asia, north Africa, and northeastern China) at about 10,000 Years B.P., and two for B. indicus (Indus Valley and southern India) at around 8000 Years B.P. Mitochondrial DNA (mtDNA) has been extensively used to reconstruct the ancestry and domestic history of cattle. Analysis of modern and ancient mtDNAs showed that auroch diversity had an uneven distribution with multiple mtDNA clades and a complex branching structure [46,47]. Since domestication process occurred in geographically restricted areas, only a fraction of the whole wild genetic variation was retained in each domestication center (sampling effect).
Regarding MHC Class II genes, trans-species theory supports the hypothesis that MHC alleles preceded the origin of related species and survived through the speciation process [48]. Phylogenetic analysis of Bovidae DRB sequences also supports this theory [49,50]. Bearing this in mind, it is expected that domestic cattle BoLA-DRB3 diversity was present in the auroch population before the divergence of major bovine types. Just as mtDNA haplogroups, BoLA-DRB3 alleles could have shown an On the other hand, genetic structure could have been diluted due to migration and gene introgression events occurred since domestication. Most of the detected alleles in Zebu were reported in the IPD-MHC database as Taurine, while a slightly smaller number of variants was shared among these two major bovine types, and only six were described as B. indicus. In contrast, Venn diagram, used to evaluate the distribution of the BoLA-DRB3 allele among group of breeds, clearly shows that the Zebu group, comprised of Bolivian Gir, Bolivian Nellore, Peruvian Nellore-Brahman and Philippine Brahman breeds, harbors 30% of Zebu privative alleles. This data further supports the hypothesis that BoLA-DRB3 diversity has an uneven distribution.
MHC class II gene diversity can be maintained by some form of balancing or overdominant selection [51,52]. To assess this hypothesis, allele diversity, gene frequency profile, HWE and Slatkin's exact neutrality tests were performed. In this sense, between 19 and 33 alleles were detected within the analyzed Zebu breeds, resulting in values of expected heterozygosity higher than 0.85. Furthermore, common alleles accounted for 65.95% of all alleles sampled in Bolivian Nellore, 67.18% in Peruvian Nellore-Brahman, and 85.00% in Bolivian Gir. Despite the fact Zebu breeds exhibit higher resistance to tropical pathogens; these results demonstrate that analyzed populations showed allele diversity and gene frequency distribution within the range reported for Taurine breeds [12,13,22].
Regarding HWE, none of the analyzed populations showed a significant excess of heterozygotes animals that could reflect overdominant selection, while only Gir gene frequency profile showed an even distribution consistent with the theoretical proportion expected under recent balancing selection. Even though this evolutionary force could not be detected in most cattle populations studied until now [12,22], this even gene distribution was observed in Yacumeño and Japanese Black breeds. Furthermore, several of the analyzed breeds showed neutrality test p values ranging from 0.025 to 0.1 which, though not significant, could reflect the same bias towards balancing selection. Conversely, Nellore and Nellore × Brahman gene frequency distributions were compatible with neutral selection.
Observed diversity at allele level was also present when analyzing the nucleotide and amino acid variability. The 13 cattle breeds analyzed by PCR-SBT varied in terms of the number and composition of BoLA-DRB3 alleles. However, with the exception of Shorthorn and Native Philippine cattle, the values observed in Zebu breeds for the π and NPD indices were similar to those estimated for other cattle breeds [12,[20][21][22]. As expected, the dn/ds differences were more prominent when only the ABS were analyzed, which is in concordance with the biological function of these sites. Of note, the average number for dn and ds across the entire exon 2, as well as the ABS, in Zebu breeds was slightly smaller than the average reported for previously analyzed breeds [12,[20][21][22].
Genetic structure and differentiation of BoLA-DRB3 in Zebu breeds were assessed by calculating F ST index and the exact G test. Both estimations revealed significant differences across all cattle populations. Noteworthy, when Zebu breeds were included in the analysis, the average F ST value was 1 to 5% higher than that obtained when analyzing taurine breeds alone [12] or populations of the same Taurine breed [22], respectively. Moreover, dendrograms based on Nei' D A or D S genetic distance were constructed using UPGMA and NJ algorithms. These trees showed similar relationship between breeds: two major clusters revealing topologies consistent with the historical and geographical origin. In addition, PCA results are in agreement with the overall clustering observed after NJ or UPGMA tree construction, which showed a clear divergence between Taurine and Zebu populations. Also, native breeds occupy an intermediate position in both analyses, probably due to some level of Zebu gene introgression [53,54]. Finally, both dendrograms and PCA showed that Chilean Hereford and Japanese Jersey were the most divergent populations among B. taurus breeds. This topology could be partially explained by the ancient genetic structure present among and within B. primigenius subspecies; and the sampling effect during the domestication process. Furthermore, subsequent selection, local adaptation and breed origin could also have contributed to the observed genetic differentiation at BoLA-DRB3 gene.
PCA based on allelic frequencies of amino acid motif of pockets implicated in antigen-binding function showed a clear differentiation between Taurine and Zebu breeds, as consequence of an enrichment of particular amino acid motif in specific pocket. These differences could be consistent with the widely reported Zebu cattle high disease tolerance/susceptibility to tropical infectious disease, such as ticks and intestinal parasites [47,55]. Further studies involving structural modeling and molecular simulation of the BoLA-DRβ protein are needed to elucidate whether these differences play a role in the function of the resultant protein.

Conclusion
The present study is the first to use PCR-SBT to examine the allelic distribution of the BoLA-DRB3 gene in Zebu populations. The results show that the historical divergence between Taurine and Zebu cattle breeds is due to origin, selection, and adaptation events, which explains the observed differences in BoLA-DRB3 gene diversity between these two major bovine types.

Additional file
(supported the genotyping cost and analyzing cost, and interpretation of data), and C(25450405)) (supported the genotyping cost and analyzing cost, interpretation of data, and in writing the manuscript), from the Japan Society for the Promotion of Science (JSPS).

Availability of data and materials
The datasets used and/or analyzed in the current study are available from the corresponding author upon reasonable request. Ethics approval All animals were handled by veterinarians from RIKEN, Universidad Austral de Chile and LAVET, in strict accordance with good animal practice following the Universidad Austral de Chile Institutional guidelines. This study was approved by the Committee on the Ethics of Animals for Research at the National University of LA PLATA (Certificate date May 26th, 2014) and by the Committee on the Ethics of Animals for Research at Universidad Austral de Chile (Certificate No. 153-2014). The animal has derived from private owners and the consent was taken in verbal form, because each farmer commonly contacts with veterinarian in same way.

Competing interests
The authors declare that they have no competing interests.

Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.